# -------------------------------------------------------------------------
# HeckeLIS code for Maple (Version of 01-06-08)
# Companion to the paper: "Longest increasing subsequences, 
#                         Plancherel-type measure and and the Hecke
#                         Insertion algorithm"
#                         by Hugh Thomas and Alexander Yong
# Report errors to Alexander Yong (ayong@math.umn.edu).
# Maple compatibility: Tested for Maple 7
# -------------------------------------------------------------------------

# -------------------------------------------------------------------------
# QUICK INSTRUCTIONS:
#
# Purpose: This software implements codes to study the Plancherel-Hecke
#          measure, and in particular computes Heckeshape of a word w.
#          Main procedures are Hecke (which actually computes the insertion
#          tableau associated to a word), and patience_ties_allowed. Other
#          procedures are included but not discussed here.
#
# Usage:   Start maple and load HeckeLIS.v0.1.txt
#          Words are inputed as an row vector, e.g., [3,4,1,5,4]
#
# Example: Type Hecke([3,4,1,5,4],5) to compute the Hecke insertion
#          of [3,4,1,5,4]. The "5" is the maximum number of rows of the
#          output. It can be set to some large number, but the alphabet 
#          size suffices. (This is only really relevant for very large 
#          computations.) 
#
#          The output is [[[1, 4, 5], [3, 5]], [[1, 3], [2, 2]]]. The
#          [[1, 4, 5], [3, 5]] is the insertion tableau, read by rows
#          i.e.
#
#          145
#          35 
#
#          The [[1, 3], [2, 2]] are the row lengths (3,2)
#
#          You can generate a random word from W_{n,q} by using
#          randword_repeatallowed(n,q); 
#          and use this with Hecke.
#
#          You can use the data [[1, 3], [2, 2]] to plot the shape
#          with plot([[1, 3], [2, 2]]);
#
#          To play patience with ties allowed, type 
#          e.g., patience_ties_allowed([3,5,1,4,2]);
#
#          MC_patience_ties_allowed(100); does some Monte Carlo
#          simulation on pile sizes (here with 100 samples).
# -------------------------------------------------------------------------

# Code to play patience with ties allowed, using the greedy strategy
# decklist is a list of integers, possibly with repeats
# outputs a list of column sizes 

`patience_ties_allowed`:=proc(decklist)

	local ii, numcols, gameboard, gameboardarray, currentend,
              failed, currentcol, lastentryindex:

	# gameboard is a list of lists
	gameboard:=[]:
	currentend:=0;

	gameboardarray:=array(1..20000);
	for ii from 1 to nops(decklist) do

	    failed:=1:
	    #special case, inserting the very first number
	    if currentend=0 then
	       gameboardarray[1]:=[op(ii,decklist)]:		
	       failed:=0:
	       currentend:=1:
	    fi:
	    # otherwise find the first column whose top card (the last
	    # entry in the list) is <= the current number to insert
	    currentcol:=1:
	    while failed=1 and currentcol<=currentend do
		lastentryindex:=nops(gameboardarray[currentcol]):  
		if(op(lastentryindex,gameboardarray[currentcol])>=op(ii,decklist)) then
		        failed:=0:
			gameboardarray[currentcol]:=[op(gameboardarray[currentcol]),op(ii,decklist)]:
		fi:
		currentcol:=currentcol+1:
	    od:

	    # if no such columns exist, create a new column with the
	    # current entry to be inserted 
	    if failed=1 then
	       gameboardarray[currentend+1]:=[op(ii,decklist)]:
	       currentend:=currentend+1:
	    fi:
	od:
	
	for ii from 1 to currentend do
	    gameboard:=[op(gameboard),gameboardarray[ii]]:
	od:
	return(gameboard);
end:

# tool to convert a permutation in S_52 to a word with repeats corresponding
# to a real deck

`convertperm52`:=proc(permin)
	
	local ii, permout:

	permout:=[]:
	for ii from 1 to nops(permin) do
	    permout:=[op(permout), ceil(op(ii,permin)/4)]:
	od:
	return(permout);
end:

# For Monte Carlo simulation with patience_ties_allowed

`MC_patience_ties_allowed`:=proc(numsamples)

	local jj, RR, summaryarray, is, avgnum, TT:

	avgnum:=0;
	summaryarray:=array(1..13);
	for jj from 1 to 13 do
	    summaryarray[jj]:=0;
	od:

	for jj from 1 to numsamples do
	    RR:=randperm[combinat](52);
	    RR:=convertperm52(RR);
	    TT:=patience_ties_allowed(RR);
	    is:=nops(TT):
	    avgnum:=avgnum+is;
	    summaryarray[is]:=summaryarray[is]+1:
	od:	    
	avgnum:=avgnum*1.0/numsamples;
	return([summaryarray, avgnum]);
end:

# constructs a uniformly chosen random composition
# x_1+..+x_q=n
# where each summand x_i is at least one

`randcomposition`:=proc(n,q, upperbnd)

	local ii, bigarray, compn, ku, indextochange, temp:

	#randomize():
	compn:=[]:
	bigarray:=array(1..upperbnd);
	for ku from 1 to upperbnd do
	    bigarray[ku]:=1:
	od:

	for ii from 1 to n-q do
	    indextochange:=floor(stats[random,uniform](1)*q)+1:
	    temp:=bigarray[indextochange]+1:
	    bigarray[indextochange]:=temp:
	od:
	
	for ii from 1 to q do
	    compn:=[op(compn),bigarray[ii]]:
	od:

	return(compn);
end:

`randcomposition_norepeat`:=proc(n,q, upperbnd)

	local ii, bigarray, compn, ku, indextochange, temp:

	#randomize():
	compn:=[]:
	bigarray:=array(1..upperbnd);
	for ku from 1 to upperbnd do
	    bigarray[ku]:=0:
	od:

	for ii from 1 to n do
	    indextochange:=floor(random[stats][uniform](1)*q)+1:
	    
	    temp:=bigarray[indextochange]+1:
	    bigarray[indextochange]:=temp:
	od:
	
	for ii from 1 to q do
	    compn:=[op(compn),bigarray[ii]]:
	od:

	return(compn);
end:

`randword`:=proc(n,q, upperbnd)

	local themultperm, randcomp, ii,jj, therandperm, newmultperm:

	themultperm:=[]:
	randcomp:=randcomposition_norepeat(n,q, upperbnd):
	for ii from 1 to nops(randcomp) do
	    for jj from 1 to op(ii,randcomp) do
		themultperm:=[op(themultperm),ii]:
	    od:
	od:
	therandperm:=combinat[ 'randperm'](n):
	newmultperm:=[]:
	for ii from 1 to n do
	    newmultperm:=[op(newmultperm),op(op(ii,therandperm), themultperm)]:
	od:
	return(newmultperm);
end:

`randword_repeatallowed`:=proc(n,q)

	local themultperm, randarray, ii:

	themultperm:=[]:
	randarray:=[stats[random,uniform](n)*q]:
	for ii from 1 to n do
	    themultperm:=[op(themultperm),floor(op(ii,randarray))+1]:
	od:
	return(themultperm);
end:
	
`randmultperm`:=proc(n,q, upperbnd)

	local themultperm, randcomp, ii,jj, therandperm, newmultperm:

	themultperm:=[]:
	randcomp:=randcomposition(n,q, upperbnd):
	for ii from 1 to nops(randcomp) do
	    for jj from 1 to op(ii,randcomp) do
		themultperm:=[op(themultperm),ii]:
	    od:
	od:
	therandperm:=combinat[ 'randperm' ](n):
	newmultperm:=[]:
	for ii from 1 to n do
	    newmultperm:=[op(newmultperm),op(op(ii,therandperm), themultperm)]:
	od:
	return(newmultperm);
end:

`w0w0`:=proc(wordin, q)

	local ii, ans:
	ans:=[]:
	for ii from 1 to nops(wordin) do
	    ans:=[op(ans),q-op(ii,wordin)+1]:
	od:
	return(ans):
end:

# For ties allowed, look at pile sizes

`MC_piles`:=proc(n,q,numsamples, upperbnd)

	local ii, jj, RR, summaryarray, is, avgnum, maxpilesize, TT, kkk,
              plotarray, variance_list, variance_est:

	maxpilesize:=-1:
	avgnum:=0;
	variance_list:=[]:
	summaryarray:=array(1..q);
	for jj from 1 to q do
	    summaryarray[jj]:=0;
	od:

	for jj from 1 to numsamples do
	    RR:=randword_repeatallowed(n,q):
	    TT:=patience_ties_allowed(RR);
	    avgnum:=avgnum+nops(TT);
	    variance_list:=[op(variance_list),nops(TT)]:
	    maxpilesize:=max(maxpilesize,nops(TT));
	    for kkk from 1 to nops(TT) do
		summaryarray[kkk]:=summaryarray[kkk]+nops(op(kkk,TT)):
	    od:
	od:	    
	
	for kkk from 1 to q do
	    summaryarray[kkk]:=summaryarray[kkk]*1.0/numsamples;
	od:	
	
	avgnum:=avgnum*1.0/numsamples:
	plotarray:=[];
	for kkk from 1 to maxpilesize do
	    plotarray:=[op(plotarray),[kkk,summaryarray[kkk]]]:
	od:

	# estimate variance
	variance_est:=0:
	for ii from 1 to nops(variance_list) do
	    variance_est:=variance_est+(op(ii,variance_list)-avgnum)^2:
	od:
	variance_est:=variance_est/(nops(variance_list)-1):
	return([summaryarray,plotarray,avgnum, variance_est]);
end:


`Hecke`:=proc(wordin, arraymax)

	local ii, bigarray, numberrows, insertedyet, rowtotry, 
              currentlabel, currentlist,
              possible, temp, tt, found, returnlist,
              anothertemp, arraydimfact, yy, plotshapelist;
	bigarray:=array(1..arraymax): 
	numberrows:=0:

	for ii from 1 to nops(wordin) do
	      # special case, first number inserted
	      if numberrows=0 then
		 bigarray[1]:=[op(1,wordin)]:
		 numberrows:=numberrows+1:
	      fi:
	      
	      insertedyet:=0:
	      rowtotry:=1:
	      currentlabel:=op(ii,wordin):
	      while insertedyet=0 do
		    if(numberrows=rowtotry-1) then
		       bigarray[rowtotry]:=[currentlabel]:
		       insertedyet:=1:
		       numberrows:=numberrows+1:
		    fi:
		    currentlist:=bigarray[rowtotry]:

		    if insertedyet=0 then
		     if(currentlist[nops(currentlist)]=currentlabel) then
			insertedyet:=1:
		     elif (currentlist[nops(currentlist)]<currentlabel) then
			# add new number to end of row if possible

			possible:=0:
			if(rowtotry=1) then
			      possible:=1:

			else 
			     # check if the previous row blocks putting in the new entry
			     if (nops(bigarray[rowtotry-1])>nops(bigarray[rowtotry])) then
				temp:=bigarray[rowtotry-1]:
				if(temp[nops(bigarray[rowtotry])+1]<currentlabel) then
					possible:=1:
				fi:
			     fi:
			fi: 
			
			if possible=1 then
			   bigarray[rowtotry]:=[op(currentlist),currentlabel]:
			fi:
			insertedyet:=1:			
		     else

			# find the smallest label larger than currentlabel and bump
			found:=0:
			tt:=1:
			while found=0 do
			      if(currentlist[tt]>currentlabel) then
				found:=1:
			      else
				tt:=tt+1:
			      fi:
			od:

			# now determine if we can replace 

			possible:=1:
			if rowtotry=1 then
			   if(tt=1) then 
				    possible:=1:
			   else
				if(currentlist[tt-1]>=currentlabel) then
				    possible:=0:
				fi:
			   fi:
			else
		           temp:=bigarray[rowtotry-1]:
			   if(temp[tt]>=currentlabel) then
				possible:=0:
			   fi:
			   if(tt=1) then
				    possible:=1:
			   else
				if(currentlist[tt-1]>=currentlabel) then
			   	    possible:=0:
			        fi:
			   fi:
			fi:

			temp:=currentlist[tt]:
			if(possible=1) then
			   anothertemp:=array(1..nops(currentlist)):
			   arraydimfact:=nops(currentlist):
			   for yy from 1 to nops(currentlist) do
			       anothertemp[yy]:=op(yy,currentlist):
			   od:
			   anothertemp[tt]:=currentlabel:
			   currentlist:=[]:
			   for yy from 1 to arraydimfact do
			       currentlist:=[op(currentlist),anothertemp[yy]]:
			   od:
			   bigarray[rowtotry]:=currentlist:
			fi:
			# in any case bump the label to the next row
			   currentlabel:=temp:
			   rowtotry:=rowtotry+1:
		     fi:
		    fi: 
	      od:
	od:
	returnlist:=[]:
	for ii from 1 to numberrows do
	    returnlist:=[op(returnlist),bigarray[ii]]:
	od:
	plotshapelist:=[]:
	for ii from 1 to numberrows do
	    plotshapelist:=[op(plotshapelist),[ii,nops(op(ii,returnlist))]]:
	od:
	return([returnlist,plotshapelist]);

end: