# ------------------------------------------------------------------------- # 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)]nops(bigarray[rowtotry])) then temp:=bigarray[rowtotry-1]: if(temp[nops(bigarray[rowtotry])+1]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: