#OK to post homework
#Blair Seidler, 2021-04-25 Assignment 25

with(combinat):
Help:=proc(): print(`HookSet(L,cell)`,`OneStepGNW(L)`,`GNW(L)`): end:

# 1.

#HookSet(L,cell): inputs a shape L and a cell [i,j], with 1≤ i ≤nops(L), 1 ≤ j ≤ L[i],
#and outputs set set of all cells that are (weakly) to the right in the same row, and
#(weakly) down in the same column. For example,
#HookSet([5,5,3,2,1],[2,2]); should output the set {[2,2],[2,3],[2,4],[2,5],[3,2],[4,2]}

HookSet:=proc(L,cell) local i,j,k,H:
i:=cell[1]: j:=cell[2]:
if i<1 or i>nops(L) or j<1 or j>L[i] then RETURN(FAIL): fi:
H:={seq([i,k],k=j..L[i])}:
for k from i+1 to nops(L) while L[k]>=j do
  H:=H union{[k,j]}:
od:
H:
end:

# 2.

#OneStepGNW(L): that implements one iteration of the beautiful Greene-Nijenhuis-Wilf algorithm,
#by starting at cell [1,1] and "walking" to a corner cell, by at each step, (if you are currently
#at location cell), picking UNIFORMLY at random a member of
#HookSet(L,cell) minus {cell}
#and going there, until you hit a corner cell. Returns that position [i,j].

#NOTE: Starting at [1,1] does not produce a uniform distribution. The GNW algorithm
#says to start at a uniformly random cell of L, so my code does that.

OneStepGNW:=proc(L) local row,col,cur,H,r:
row:=LD(L):
col:=rand(1..L[row])():
cur:=[row,col]:
H:=HookSet(L,cur) minus {cur}:
while nops(H)>0 do
  cur:=H[rand(1..nops(H))()]:
  H:=HookSet(L,cur) minus {cur}:
od:
cur:
end:

# 3.

#GNW(L): does the same as randomSYT(L).
GNW:=proc(L) local T,L1,i,n,cell:
n:=add(L):
L1:=L:
T:=[seq([1$L[i]],i=1..nops(L))]:
for i from n to 2 by -1 do
  cell:=OneStepGNW(L1):
  T[cell[1]][cell[2]]:=i:
  L1[cell[1]]:=L1[cell[1]]-1:
  if nops(L1[-1])=0 then L1:=op(1..nops(L1)-1,L1): fi:
od:
T:
end:

#Convince yourself that GNW([10,10,10]) is MUCH faster than randomSYT([10,10,10]).
#How long does it take to get GNW([9$9]);

#Both were really fast on [10,10,10], and even on [9$9].
#at [15$9], there was a significant difference (yes, I was restarting between trials)
#time(randomSYT([15 $ 9]));
#                             33.390
#time(GNW([15 $ 9]));
#                             0.046


# 4.

# syt([a1,a2,..., ak])=n!/Product(hij, over all the n cells [i,j] of L)
# syt([a1,..., ak])= Delta(a1+k-1,a2+k-2, ..., ak)*(a1+...+ak)!/((a1+k-1)!(a2+k-2)!....ak!)
#Where Delta(x1,...xk) is the discriminant, i.e the product of xi-xj over all 1 ≤ i ≤ j ≤ k.

# These two formulas are really just two ways of writing the same thing.
# First observation (a1+...+ak)! and n! are clearly identical.

#Of the remaining factors, the first formula has all of the hook lengths in the denominator.

#The RDB (YF) formula has a discriminant in the numerator and a bunch of factorials in the denom.
#The factorials are the factorials of the first-column hook lengths.
#The discriminant in the numerator is actually all of the "missing" hook lengths from each row.

#Example: in [5,4,2], the first-row hook lengths are (7,6,4,3,1).
#The terms of the discriminant involving a1 are
#((a1+k-1)-(a2+k-2))((a1+k-1)-(a3))=(7-5)(7-2)=(5)(2).
# Now we note that (5)(2)/7! is exactly 1/(7*6*4*3*1)
#The second-row hook lengths are (5,4,2,1), and the other term of the disciminant is 3.
#The last row will never have any missing hook lengths, so it will just be (ak)!


#### From C25.txt ####
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:

#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:


#### From C24.txt ####
#C24.txt: Maple code for Lecture 24#OK to post homework

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:



#From George Spahn, 4/18/2021, Assignment 23

#2

# SYTpairs(n) inputs a positive integer n and outputs the set of pairs of 
# standard Young tableaux of the same shape [S,T] each with n cells

SYTpairs:=proc(n) local s,L,T,i,j:

L:={}:
for s in Pn(n) do:
 T:=SYT(s):
 L := L union {seq(seq([T[i],T[j]] ,j=1..nops(T)) ,i=1..nops(T))}:
od:
L:
end: 

###############################################################################
#3


#RS1m(Y,m): ONE STEP OF THE RS algorithm
# modified to also return the row of the insertion
RS1m:=proc(Y,m) local i,Y1:

Y1:=Ins(Y,1,m):
if Y1[2]=0 then
  RETURN(Y1[1],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],i-1:
end:

#RS1mOld(Y,m): ONE STEP OF THE RS algorithm
# modified to also return the row of the insertion
RS1mOld:=proc(Y,m) local i,Y1:

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

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

#RS(pi):  the pair [S,T] that is the output
#of the Robinson-Scheastead mapping
RS:=proc(pi) local Y1,Y2,i,row:
Y1:=[[pi[1]]]:
Y2:=[[1]]:
for i from 2 to nops(pi) do
 Y1,row:=RS1m(Y1,pi[i]):
 if row = nops(Y2)+1 then 
  Y2:=[op(Y2),[i]]:
 else 
  Y2[row]:=[op(Y2[row]),i]:
 fi:
od:
[Y1,Y2]:
end:

#############################################################################
#4

invert1 := proc(T1,row) local T,a,i,j,c:
T:=T1:
i:=row:
a := T[i][-1]:
if nops(T[i]) = 1 then 
 T:=[op(1..nops(T)-1, T)]:
else
 T[i]:=[op(1..nops(T[i])-1, T[i])]:
fi:
while i>1 do 
 i:=i-1:
 for j from 2 to nops(T[i]) while T[i][j] < a do od:
 c:=T[i][j-1]:
 T[i][j-1]:=a:
 a:=c:
od:
T,a:
end:

# InvRS(P) inputs a pair [S,T] of standard Young tableaux of the same shape 
# and outputs a permutation of 1,2,...n
InvRS:=proc(P) local S,T,perm,n,x,j,a,i:
S:=P[1]:
T:=P[2]:
perm := []:

n:= max( seq(max(T[i]),i=1..nops(T)) ):
for x from n to 1 by -1 do

 for j from 1 to nops(T) while not(x in T[j]) do od:
 S,a:=invert1(S,j):
 perm:=[a,op(perm)]:  

od:
perm:
end:

##############################################################################
#old stuff

#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:

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

#InsOld(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
InsOld:=proc(Y,i,m) local 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:

#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 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:

#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,j:
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