#C26.txt: maple code for Lecture 26, April 26, 2021
Help26:=proc() print(` GNW(L), detpart(L,x,m), Sc(L,x,m)`):
print(`ScY(L,x,m),delta(x,m), sytC(L), Is2First(Y) `):
print(`Est2First(L,K), ekn(k,n,x), hkn(k,n,x)`):
end:

#Elementary symmetic polynomial
#e_k(x[1],...,x[n]):
#coeff of t^(n-k) in(t+x[1])*(t+x[2])*...*(t+x[n])
#weight enumerator of all subests of k elements in {1, ...,n}
#weight({2,3,5})=x[2]*x[3]*x[5]:
#ekn(2,3,x)= weight({{1,2},{1,3},{2,3}))=
#x[1]*x[2]+x[1]*x[3]+x[2]*x[3]
#Take S any subest of k elements of {1,...,n}
#Case 1: n does not belong to S
#ekn(k,n-1,x)
#Case 2: n does belong
#x[n]*ekn(k-1,n-1,x):
#ekn(k,n,x)=
ekn:=proc(k,n,x) option remember:
if k=0 then
  RETURN(1):
fi:
if k>n then
 RETURN(0):
fi:
expand(ekn(k,n-1,x)+x[n]*ekn(k-1,n-1,x)):
end:


#Full homog. symmetic polynomial
#h_k(x[1],...,x[n]):
#coeff of t^k 1/(1-x[1]*t)/(1-x[2]*t)...
#weight enumerator of all subests of 
#1<= i1<=i2 <=i3<= ...<=ik <=n
#weight({2,3,5})=x[2]*x[3]*x[5]:
#hkn(2,3,x)= weight({[1,2],[1,3],[2,3],[1,1],[2,2],[3,3]))=
#x[1]*x[2]+x[1]*x[3]+x[2]*x[3]+x[1]^2+x[2]^2+x[3]^2
#Take S any subest of k elements of {1,...,n}
#Case 1: n does not belong to S
#hkn(k,n-1,x)
#Case 2: n does belong
#x[n]*hkn(k-1,n,x):
#hkn(k,n,x)=
hkn:=proc(k,n,x) option remember:
if k=0 then
  RETURN(1):
fi:
if n=0 then
  RETURN(0):
fi:
expand(hkn(k,n-1,x)+x[n]*hkn(k-1,n,x)):
end:



 


#Is2First(Y): Inputs a Standard Young tableau Y
#decides whether 2 is in the row 
Is2First:=proc(Y): 
if Y=[] or nops(Y[1])=1 then RETURN(false): fi:
evalb(Y[1][2]=2):
end:

#Est2First(L,K): estimates the probability of the
#event "2 is in the first row" by doing a public opinion
#poll among K random Standard Young tableaux (using GNW)
Est2First:=proc(L,K) local i:
evalf(coeff(add(Is2First(GNW(L)),i=1..K),true,1)/K):
end:

#sytC(L): The number of standard Young tableaux of shape L
#using Young-Frobenius
sytC:=proc(L) local i,k:
k:=nops(L):
add(L)!/mul((L[i]+k-i)!,i=1..k)
*delta([seq(L[i]+k-i,i=1..k)],k)
:
end:

#ScY(L,x,m): The Schur polynomial in x[1],..., x[m]
#of the shape L, according to the determinant formula
#rediscovered by Yuxuan
ScY:=proc(L,x,m): normal(detpart(L,x,m)/delta(x,m)):end:

#FROM YUXUAN
#Yuxuan Yang, April 24th, Assignment 25

with(combinat):

#1
HookSet:=proc(L,cell) local i,j,S:
S:={}:
for i from cell[1] to nops(L) while L[i]>=cell[2] do
 S:=S union {[i,cell[2]]}:
od:
for j from cell[2] to L[cell[1]] do
 S:=S union {[cell[1],j]}:
od:

S:
end:
 
#2
OneStepSNW:=proc(L) local m,n,i,j,H,k,allpos,pos:
allpos:=[seq(seq([n,m],m=1..L[n]),n=1..nops(L))]:
pos:=allpos[rand(1..nops(allpos))()]:
i:=pos[1]:
j:=pos[2]:
for k from 1 to add(L) do
H:=HookSet(L,[i,j]):
if nops(H)=1 then 
 RETURN([i,j]):
fi:
pos:=H[rand(2..nops(H))()]:
i:=pos[1]:
j:=pos[2]:
od:
pos:
end:

#3

GNW:=proc(L) local i,j,pos,YT,n,L1:
n:=add(L):
if nops(L)=1 then
RETURN([[seq(i,i=1..L[1])]]):
fi:
pos:=OneStepSNW(L):
L1:=L:
if L1[pos[1]]>1 then
L1[pos[1]]:=L1[pos[1]]-1:
YT:=GNW(L1):
YT[pos[1]]:=[op(YT[pos[1]]),n]:
RETURN(YT):
else
L1:=[op(1..nops(L1)-1,L1)]:
YT:=GNW(L1):
YT:=[op(YT),[n]]:
RETURN(YT):
fi:
end:




#5
delta:=proc(x,n) local j,i:
product(product(x[i]-x[j],j=i+1..n),i=1..n-1):
end:

guesssc:=proc(L,x,m) 
expand(Sc(L,x,m)*delta(x,m)):
end:

with(linalg):
detpart:=proc(L,x,m) local i,j,M,degree,size,xnew:
size:=max(m,nops(L)):
degree:=[0$size]:
xnew:=[1$size]:
for i from 1 to m do
degree[i]:=degree[i]+m-i:
xnew[i]:=x[i]:
od:
for i from 1 to nops(L) do
degree[i]:=degree[i]+L[i]:
od:
M:=[seq([seq(xnew[i]^degree[j],j=1..size)],i=1..size)]:
expand(det(matrix(M))):
end:

checksc:=proc(L,x,m)
evalb(guesssc(L,x,m)=detpart(L,x,m)):
end:

#I checked several examples. The formula should work.
#The formula is in fact a generalized vandermonde det over the regular vandermonde det.
#6 no idea how to prove.




















#C25.txt: Maple code for Lecture 25 of Dr.Z.'s Experimental Mathematica class, Spring 2021

Help25:=proc(): print(`Wt(L,x),Sc(L,x,m), randomSYT(L)`):end:

#Wt(L,x): The weight according to x[1],..., x[m] of
#the tableau L
Wt:=proc(L,x) local i,j: mul(mul(x[L[i][j]],j=1..nops(L[i])),i=1..nops(L)):end:

ScStupid:=proc(L,x,m) local S,i: S:=TAB(L,m): add(Wt(S[i],x),i=1..nops(S)):end:



#Sc(L,x,m): the sum  of all the weights of tableaux of shape L with entries
#{1,2,...,m}  x is the symbol for
#the weight. These are called Schur polynomials
#
Sc:=proc(L,x,m)  local L1,B,i,S,b,y,S1:
option remember:

if L=[] then
  RETURN(1):
fi:

if m<nops(L) then
  RETURN(0):
fi:

if nops(L)=1 and m=1 then
 RETURN(x[1]^L[1]):
fi:


B:=NewBox([seq(L[i]-L[i+1],i=1..nops(L)-1),L[nops(L)]]):

S:=0:

for b in B do
 L1:=[seq(L[i]-b[i],i=1..nops(L))]:
 if L1[-1]=0 then
  L1:=[op(1..nops(L1)-1,L1)]:
fi:

 

S:=expand(S+ x[m]^add(b)*Sc(L1,x,m-1)):

od:
S:

end:

#THANKS TO Blair Seidler
#randomSYT(L): inputs an integer partition L (aka shape) and outputs, uniformly at random,
#one of the syt(L) Standard Young tableaux of that shape
randomSYT:=proc(L) local L1,die,S1,n,i:
n:=add(L):
if n=1 then RETURN([[1]]): fi:
die:=[]:
for i from 1 to nops(L)-1 do
  if L[i]>L[i+1] then
    L1:=[op(1..i-1,L),L[i]-1,op(i+1..nops(L),L)]:
    die:=[op(die),syt(L1)]:
  else
    die:=[op(die),0]:
  fi:
od:

if L[-1]>1 then
  L1:=[op(1..nops(L)-1,L),L[i]-1]:
else
  L1:=[op(1..nops(L)-1,L)]:
fi:
die:=[op(die),syt(L1)]:
i:=LD(die):
L1:=[op(1..i-1,L),L[i]-1,op(i+1..nops(L),L)]:
if L1[-1]=0 then L1:=[op(1..nops(L1)-1,L1)]: fi:
S1:=randomSYT(L1):
AP
(S1,i,n):
#S:

end:

#LD(L): Inputs a list of positive integers L (of n:=nops(L) members)
#outputs an integer i from 1 to n with the prob. of i being
#proportional to L[i]
#For example LD([1,2,3]) should output 1 with prob. 1/6
#output 2 with prob. 1/3
#output 3 with prob. 3/6=1/2
LD:=proc(L) local n,i,su,r:
n:=nops(L):
r:=rand(1..convert(L,`+`))():
su:=0:

for i from 1 to n do
 su:=su+L[i]:
  if r<=su then
    RETURN(i):
  fi:
od:
end:



#OLD STUFF
Help24:=proc(): print(`SYTpairs(n), RS(pi), InvRS(P) , Words(m,n)`): 
print(`APP(Y,b,m), Box(L) , NewBox(B) TAB(L,m) `):
end:

#APP(Y,b,m): inputs a tableaux Y with entires {1,...,m-1} and
#and a vector b of either the same length of Y or one longer and appends
#b[i] m's to the end of row i for i=1..nops(Y) and if nops(b)=nops(Y)+1 a new row of b[nops(b)] m's. For example
#APP([[1,1,2],[2,2]],[1,1,1],3) should be
#[[1,1,2,3],[2,2,3],[3]]
APP:=proc(Y,b,m) local i,Y1:
Y1:=[seq([op(Y[i]),m$b[i]],i=1..nops(Y))]:
if nops(b)=nops(Y)+1 then
 Y1:=[op(Y1),[m$b[-1]]]:
fi:
Y1:
end:

#TAB(L,m): the set of all tableaux of shape L with entries
#{1,2,...,m}
TAB:=proc(L,m)  local L1,B,i,S,b,y,S1:
option remember:

if L=[] then
  RETURN({[]}):
fi:

if m<nops(L) then
  RETURN({}):
fi:

if nops(L)=1 and m=1 then
 RETURN({[[1$L[1]]]}):
fi:


B:=NewBox([seq(L[i]-L[i+1],i=1..nops(L)-1),L[nops(L)]]):

S:={}:

for b in B do
 L1:=[seq(L[i]-b[i],i=1..nops(L))]:
 if L1[-1]=0 then
  L1:=[op(1..nops(L1)-1,L1)]:
fi:

 S1:=TAB(L1,m-1):

S:=S union {seq(APP(y,b,m),y in S1)}:

od:
S:

end:

#NewBox(B): inputs a list of NON-NEGATIVE integers and outputs
#the set of all lists [i1,i2,...] i1<=B[1], i2<=B[2] ,...
NewBox:=proc(B) local k,B1,S1,i,v:
option remember:
if B=[]  then
 RETURN({[]}):
fi:

k:=nops(B):
B1:=[op(1..k-1,B)]:

S1:=NewBox(B1):
{seq(seq([op(v),i], v in S1),i=0..B[k])}:

end:

#Box(L): Inputs a list L (where k=nops(L)) of non-negative integers L outputs
#the  LIST of all LISTS [a1,a2,..., ak] in LEX ORDER such that ai is in {0,1,.., L[i])
#In particular the n-dim UNIT cube would be Box([1$n]);
Box:=proc(L) local k,i,L1,B1,j:
option remember:

if not type(L,list) then
  print(L, `should be a list `):
 RETURN(FAIL):
fi:

k:=nops(L):
if k=0 then 
 RETURN([[]]):
fi:
 if not ({seq(type(L[i],integer),i=1..nops(L))}={true} and min(op(L))>=0) then
   print(`bad input`):
  RETURN(FAIL):
 fi:




L1:=[op(1..k-1,L)]:

B1:=Box(L1):
 [ seq(seq([op(B1[j]),i],i=1..L[k]),j=1.. nops(B1))]:
end:

#C23.txt:  Maple code for Lecture 23

Help23:=proc(): print(` Ins(Y,i,m), RS1(Y,m) , RSleft(pi), AllSYT(n) `): end:

#Ins(Y,i,m): Inputs a (partially filled tableux Y,
#a row i, between 1 and nops(Y)+1 and an integer m
#and tries to place it LEGALLY at row i
#(if it is larger than all the current entries of row i
#and does not stick out, i=1 OR nops(Y[i])<nops(Y[i-1])
#and outputs the new tableau and either 0 or the new bumpbee
Ins:=proc(Y,i,m) local Y1,j:
if not (1<=i and i<=nops(Y)+1) then
 RETURN(FAIL):
fi:

if i=nops(Y)+1 then
  RETURN([op(Y),[m]],0):
fi:

if i=1 then 
  if m>=Y[1][-1] then
    RETURN([[op(Y[1]),m], op(2..nops(Y),Y)],0):
   fi:
fi:

for j from 1 to nops(Y[i]) while m>Y[i][j] do od:

if j=nops(Y[i])+1  and (i=1 or nops(Y[i])<nops(Y[i-1])) then

#corrected after class , Maple does not like op(m..n,L) if m-n>=2
if j+1<=nops(Y[i]) then
  RETURN ([op(1..i-1,Y),[op(1..j-1,Y[i]),m,op(j+1..nops(Y[i]),Y[i])],
op(i+1..nops(Y),Y)],0):
else

RETURN ([op(1..i-1,Y),[op(1..j-1,Y[i]),m],
op(i+1..nops(Y),Y)],0):
fi:


else
RETURN([op(1..i-1,Y),[op(1..j-1,Y[i]),m,op(j+1..nops(Y[i]),Y[i])],
op(i+1..nops(Y),Y)],Y[i][j]):
fi:

end:

#RS1(Y,m): ONE STEP OF THE RS algorithm
RS1:=proc(Y,m) local i,Y1:

Y1:=Ins(Y,1,m):
if Y1[2]=0 then
  RETURN(Y1[1]):
fi:

#continued after class
for i from 2 to nops(Y)+1 while Y1[2]<>0 do
Y1:=Ins(Y1[1],i,Y1[2]):
od:
Y1[1]:
end:

with(combinat):

#RSleft(pi): The left tableau,S of the pair (S,T) that is the output
#of the Robinson-Scheastead mapping
RSleft:=proc(pi) local Y,i:
Y:=[[pi[1]]]:
for i from 2 to nops(pi) do
 Y:=RS1(Y,pi[i]):
od:
Y:
end:

#AllSYT(n): The set of all Standard Young Tableau with n cells
AllSYT:=proc(n) local S:
S:=Pn(n):
{seq(op(SYT(S[i])),i=1..nops(S))}:
end:

#OLD STUFF
#C22.txt
Help22:=proc(): print(` Pn(n), AP(Y,i,m), SYT(L) , AllSYT(n) `): end:

#Ap(Y,i,m):write a procedure that inputs a Standard Young Tableau
#and an integer i between 1 and nops(Y)+1 and inserts the
#entry m at  the end of the i-th row, or if i=nops(Y)+1 make
#a new row with 1 cell occupied with m
AP:=proc(Y,i,m):

if i=nops(Y)+1 then
 RETURN([op(Y),[m]]):
fi:


if i>1 and nops(Y[i])=nops(Y[i-1]) then
  RETURN(FAIL):
fi:

[op(1..i-1,Y),[op(Y[i]),m],op(i+1..nops(Y),Y)]:


end:

#SYT(L): inputs an integer partion (given as a weakly decreasing
#list of POSITIVE integers OUTPUTS the LIST of
#Standard Young Tableax of shape L (n=Number of boxes of L)
#(A way of filling 1, ..., n in the boxes such that horiz. and
#vertically they are increasing
SYT:=proc(L) local n,L1,S,S1,i,j:
option remember:
n:=add(L):

if n=1 then 
  RETURN([[[1]]]):
fi:

S:=[]:

for i from 1 to nops(L)-1 do
 if L[i]>L[i+1] then
  L1:=[op(1..i-1,L),L[i]-1,op(i+1..nops(L),L)]:
  S1:=SYT(L1):
 S:=[op(S), seq(AP(S1[j],i,n),j=1..nops(S1))]:
 fi:
od:


 if L[-1]>1 then
  L1:=[op(1..nops(L)-1,L), L[nops(L)]-1 ]:
  else
  L1:=[op(1..nops(L)-1,L)]:
 fi:

 S1:=SYT(L1):
 S:=[op(S), seq(AP(S1[j],nops(L),n),j=1..nops(S1))]:

S:

end:



#syt(L): inputs an integer partion (given as a weakly decreasing
#list of POSITIVE integers OUTPUTS the NUMBER of
#Standard Young Tableax of shape L (n=Number of boxes of L)
#(A way of filling 1, ..., n in the boxes such that horiz. and
#vertically they are increasing
syt:=proc(L) local n,L1,S,S1,i,j:
option remember:
n:=add(L):

if n=1 then 
  RETURN(1):
fi:

S:=0:

for i from 1 to nops(L)-1 do
 if L[i]>L[i+1] then
  L1:=[op(1..i-1,L),L[i]-1,op(i+1..nops(L),L)]:
  S:=S+syt(L1):
 
fi:
od:


 if L[-1]>1 then
  L1:=[op(1..nops(L)-1,L), L[nops(L)]-1 ]:
  else
  L1:=[op(1..nops(L)-1,L)]:
 fi:

 S:=S+syt(L1):
S:
end:

##FROM C11.txt
#Pnk(n,k): The LIST of integer partitions of n with largest part k
Pnk:=proc(n,k) local k1,L,L1:
option remember:

if not (type(n,integer) and type(k,integer) and n>=1 and k>=1 )then
 RETURN([]):
fi:
 
if k>n then
 RETURN([]):
fi:

if k=n then
 RETURN([[n] ]):
fi:

L:=[]:

for k1 from min(n-k,k) by -1 to 1 do
 L1:=Pnk(n-k,k1):
 L:=[op(L), seq([k, op(L1[j])],j=1..nops(L1))]:
od:

L:

end:

#Pn(n): The list of integer partitions of n in rev. lex. order
Pn:=proc(n) local k:option remember:[seq(op(Pnk(n,n-k+1)),k=1..n)]:end:
#END FROM C11.txt


#END OLD STUFF