print(`This is a Maple package companion to the paper titled`):
print(`"An Inhomogeneous Multispecies TASEP on a Ring",`):
print(`by Arvind Ayyer & Svante Linusson.`):
print(`The package implements the explicit formulas for the Markov chains`):
print(`for the multispecies TASEP and mulitiline queues described in the paper.`):
print(``):
print(`Type Help(); to view available procedures,`):
print(`or Help(procedure_name); to view details of the procedure.`):

#Version of June 1, 2012.

Help := proc()

if nops([args])=0 then
   print(`Procedures available are:`):
   print(`MultiSpeciesMatrix, MultiSpeciesSteadyState,`):
   print(`MultiPermtoMultiline, MultilinetoMultiPerm,`):
   print(`MultilineMatrix, MultilineSteadyState`):
   print(`CoupeMatrix, CoupeSteadyState`):
fi:

if nops([args])=1 and op(1,[args])=`MultiSpeciesMatrix` then
   print(`MultiSpeciesMatrix(m,x): returns the transition matrix of the`):
   print(`multispecies TASEP with m_1 1's, m_2 2's, and so on.`):
   print(`For example, try MultiSpeciesMatrix([1,1,1],x);`):
fi:

if nops([args])=1 and op(1,[args])=`MultiSpeciesSteadyState` then
   print(`MultiSpeciesSteadyState(m,x,p) returns the stationary distribution p(C)`):
   print(`of the multispecies TASEP with m_1 1's, m_2 2's, and so on.`):
   print(`For example, try MultiSpeciesSteadyState([1,1,1],x,p);`):
fi:

if nops([args])=1 and op(1,[args])=`MultiPermtoMultiline` then
   print(`MultiPermtoMultiline(alpha):  given a multiline queue returns the bully path`):
   print(`projection to a multispecies configuration and positions of the covered zeros.`):
   print(`For example, try MultilinetoMultiPerm([[1, 0, 0], [0, 1, 1], [1, 1, 1]]);`):
fi:

if nops([args])=1 and op(1,[args])=`MultilinetoMultiPerm` then
   print(`MultilinetoMultiPerm(pi): given a multispecies configuration`):
   print(`returns all multiline queues which project via bully paths to pi`):
   print(`For example, try MultilinetoMultiPerm([1,1,2,2,3]);`)
fi:

if nops([args])=1 and op(1,[args])=`MultilineMatrix` then
   print(`MultilineMatrix(m,x): returns the transition matrix of the`):
   print(`multiline queue which projects via bully paths to m_1 1's, m_2 2's, and so on.`):
   print(`For example, try MultilineMatrix([1,1,1],x);`):
fi:

if nops([args])=1 and op(1,[args])=`MultilineSteadyState` then
   print(`MultilineSteadyState(m,x,p) returns the stationary distribution p(C)`):
   print(`of the multiline queue which projects via bully paths to m_1 1's, m_2 2's, and so on.`):
   print(`For example, try MultilineSteadyState([1,1,1],x,p);`):
fi:

if nops([args])=1 and op(1,[args])=`CoupeMatrix` then
   print(`CoupeMatrix(m,x): returns the transition matrix of the`):
   print(`three species coupe process which projects via bully paths to m_1 1's, m_2 2's, and m_3 3's.`):
   print(`For example, try CoupeMatrix([1,1,1],x);`):
fi:

if nops([args])=1 and op(1,[args])=`CoupeSteadyState` then
   print(`CoupeSteadyState(m,x,p) returns the stationary distribution p(C)`):
   print(`of the multispecies TASEP which projects via bully paths to m_1 1's, m_2 2's,  and m_3 3's.`):
   print(`For example, try CoupeSteadyState([1,1,1],x,p);`):
fi:

end:

with(LinearAlgebra):

#################################################GRAPH CREATION############################################################
#MultiSpeciesGraph(config) creates the relevant graph for the closed system with edges(0,...,n-1)^d and vertices for transition elements assuming the higher colour can be exchanged with the lower colour in the bulk (there is no boundary)
MultiSpeciesGraph := proc(config,x) local i,j,k,Ve2,Ed,Ve,tem,sites:

sites := add(i,i in config):
Ve := combinat[permute]([seq(i$(config[i]),i=1..nops(config))]):
Ed := {}:

for i from 1 to nops(Ve) do
	if Ve[i][1] < Ve[i][sites] then
		Ed := Ed union {[x[Ve[i][1]], Ve[i],[Ve[i][sites],op(2..sites-1,Ve[i]),Ve[i][1]]] }:
	fi:

	for j from 1 to sites-1 do
		if Ve[i][j] > Ve[i][j+1] then
			Ed := Ed union { [x[Ve[i][j+1]],Ve[i], [op(1..j-1,Ve[i]), Ve[i][j+1], Ve[i][j], op(j+2..sites,Ve[i])]] }:
		fi:
	od:
od:

return Ve,Ed:

end:
	

#rotate1(config) returns the config translated by 1 step assuming it is on the ring
rotate1 := proc(config) local i:
return [config[nops(config)],op(1..nops(config)-1,config)]:
end:

#MultiSpeciesMatrix2(config,x) returns the transitions of the Markov chain
MultiSpeciesMatrix2 := proc(config,x) local i,j,k,M,Ve,Ed:

(Ve,Ed) := MultiSpeciesGraph(config,x):
M := Array(1..nops(Ve),1..nops(Ve)):
for i from 1 to nops(Ed) do
    for j from 1 to nops(Ve) do
        for k from 1 to nops(Ve) do
    	    if Ed[i][2] = Ve[j] and Ed[i][3] = Ve[k] then
	       M[k,j] := M[k,j] + Ed[i][1]:
	    fi:
	od:
    od:
od:

return Matrix(M):

end:

#MultiSpeciesMatrix(config,x) returns the transition matrix
MultiSpeciesMatrix := proc(config,x) local i,j,M:

M:=MultiSpeciesMatrix2(config,x):
for i from 1 to op(1,M)[1] do
    M[i,i] := -add(M[j,i],j=1..i-1)-add(M[j,i],j=i+1..op(1,M)[1]):
od:

return M:

end:

#MultiSpeciesSteadyState(config,x,p) returns the steady state
MultiSpeciesSteadyState := proc(config,x,p) local i,j,M,V,V2,eqs,vars,soln,n:

n := nops(config):
M := MultiSpeciesMatrix(config,x):
V2 := MultiSpeciesGraph(config,x)[1]:
V := Vector([seq(p[op(V2[i])],i=1..nops(V2))]):
vars := convert(V,set):
eqs := convert(M.V,set) union {vars[nops(vars)]=mul(x[i]^binomial(n-i,2),i=1..n-1)}:
eqs := eqs union {seq(p[op(V2[i])]=p[op(rotate1(V2[i]))],i=1..nops(V2))}:

soln := solve(eqs,vars):
return {seq(op(1,soln[i])=factor(op(2,soln[i])),i=1..nops(soln))}:

end:

#################################################MULTI LINE PROCESS FOR PERMUTATIONS#################################################
#MultilinetoPerm(C) given a multiline process returns the Ferrari-Martin permutation
MultilinetoPerm := proc(C) local i,j,k,pos,pos2,n,p,C2:

n := nops(C)+1:
C2 := C:
p := [0$n]:
for i from 1 to n-1 do
    pos := seq(`if`(C2[i,k]=1,k,op({})),k=1..n):
    for j from i to n-2 do
    	pos2 := Findnext1(C2,j,pos):
	C2[j,pos]:=0:
	pos := pos2:
    od:
    C2[n-1,pos]:=0:
    p[pos]:=i:
od:

p[seq(`if`(p[i]=0,i,op({})),i=1..n)]:=n:
return p:

end:

#PermtoMultiline(P) given a permutation returns all the multiline configs
PermtoMultiline := proc(P) local i,j,n,M,x:

n := nops(P):
M := MultilineVert(n):
return {seq(`if`(MultilinetoPerm(M[i])=P,M[i],op({})),i=1..nops(M))}:

end:

#Findnext1(C,m,k) finds the 1 in row m+1 corresponding to the 1 in C[m,k]
Findnext1 := proc(C,m,k) local i,j,n:

if C[m,k]=0 then
   ERROR(`C[m,k] is not equal to 1`):
fi:

n := nops(C)+1:
for i from k to n while C[m+1,i]=0 do od:
if i=n+1 then
   for i from 1 to k-1 while C[m+1,i]=0 do od:
fi:
return i:

end:

#Clockseq(C,k) returns the list of positions for the clock to strike on the multiline
#starting with position k on the lowest
Clockseq := proc(C,k) local i,j,n,S:

n := nops(C)+1:
S := [0$(n-2),k]:
j := k:
for i from n-1 to 2 by -1 do
    if C[i,j]=1 then
       S[i-1]:=j:
    else
       S[i-1]:=`if`(j<n,j+1,1):
    fi:
    j := S[i-1]:
od:

return S:

end:


#MultilineVertPerm(n) returns the n-1 multiline process configs and transition rates
MultilineVertPerm := proc(n) local i,j,k,c,C,configs:

C:=[0$(n-1)]:
for i from 1 to n-1 do
    C[i] := combinat[permute]([1$i,0$(n-i)]):
od:

configs := combinat[permute]([seq(i$(n-1),i=1..binomial(n,floor(n/2)))],n-1):
configs := {op(configs)} minus {seq(seq(`if`(c[i]>binomial(n,i),c,op({})),i=1..n-1),c in configs)}:
configs := {seq([seq(C[j][configs[i][j]],j=1..n-1)],i=1..nops(configs))}:

return configs:

end:

#MultilineGraphPerm(n,x) returns the n-1 multiline process configs and transition rates
MultilineGraphPerm := proc(n,x) local i,j,k,c,C,configs,Ed,cseq,cseqm1:

configs := MultilineVertPerm(n):
Ed:={}:
for c in configs do
    for i from 1 to n do
    	cseq := Clockseq(c,i):
	cseqm1 := [seq(`if`(cseq[j]>1,cseq[j]-1,n),j=1..n-1)]:
	Ed := {op(Ed), [x[MultilinetoPerm(c)[i]],
			c,
			[seq(`if`(c[j][cseqm1[j]]=0 and c[j][cseq[j]]=1, [seq(`if`(k=cseqm1[j],1,`if`(k=cseq[j],0,c[j][k])),k=1..n)],c[j]),j=1..n-1)]]}:
    od:
od:

return configs,Ed:

end:

#MultilineMatrixPerm2(n,x) returns the transitions of the Markov chain
MultilineMatrixPerm2 := proc(n,x) local i,j,k,M,Ve,Ed:

(Ve,Ed) := MultilineGraphPerm(n,x):
M := Array(1..nops(Ve),1..nops(Ve)):
for i from 1 to nops(Ed) do
    for j from 1 to nops(Ve) do
        for k from 1 to nops(Ve) do
    	    if Ed[i][2] = Ve[j] and Ed[i][3] = Ve[k] then
	       M[k,j] := M[k,j] + Ed[i][1]:
	    fi:
	od:
    od:
od:

return Matrix(M):

end:

#MultilineMatrix(n,x) returns the transition matrix
MultilineMatrix := proc(n,x) local i,j,M:

M:=MultilineMatrixPerm2(n,x):
for i from 1 to op(1,M)[1] do
    M[i,i] := -add(M[j,i],j=1..i-1)-add(M[j,i],j=i+1..op(1,M)[1]):
od:

return M:

end:

#rotate2(config) returns the config translated by 1 step assuming it is on the ring
rotate2 := proc(config) local i,c,n:
n:=nops(config)+1:
return [seq([c[n],seq(c[i],i=1..n-1)],c in config)]:
end:


#MultilineSteadyStatePerm(n,x,p) returns the steady state
MultilineSteadyStatePerm := proc(n,x,p) local i,j,M,V,V2,eqs,vars,soln:

M := MultilineMatrixPerm(n,x):
V2 := MultilineVertPerm(n,x):
V := Vector([seq(p[op(V2[i])],i=1..nops(V2))]):
vars := convert(V,set):
eqs := convert(M.V,set) union {vars[1]=mul(x[i]^binomial(n-i,2),i=1..n-1)}:
eqs := eqs union {seq(p[op(V2[i])]=p[op(rotate2(V2[i]))],i=1..nops(V2))}:

soln := solve(eqs,vars):
return {seq(op(1,soln[i])=factor(op(2,soln[i])),i=1..nops(soln))}:

end:


#conjss(n,x,p) returns the conjectured steady state for rates [x,1,...,1]
conjss := proc(n,x,p) local i,j,V,M,S:

V := MultilineVertPerm(n,x):
#S := {seq(p[op(V[i])]=x^(binomial(n-1,2)-nops(MultilinetoMultiPerm(V[i])[2][1])),i=1..nops(V))}:
S := Vector([seq(x^(binomial(n-1,2)-nops(MultilinetoMultiPerm(V[i])[2][1])),i=1..nops(V))]):
M := MultilineMatrixPerm(n,[x,1$(n-1)]):
return convert(simplify(M.S),set):

end:

#In the case when x[1]=q and all other parameters are set to 1 and when there is one species of each type, 
#then we have the following formula for the partition function
partfn:=(n,q)->n*mul(subs({x[1]=1,seq(x[i]=q,i=2..n)},homsym(j+2,n-1-j,x)),j=1..n-2);

#################################################MULTI LINE PROCESS FOR MULTIPERMUTATIONS#################################################
#MultilinetoMultiPerm(C) given a multiline configuration returns the Ferrari-Martin multipermutation
MultilinetoMultiPerm := proc(C) local i,j,k,pos,pos2,n,r,p,C2,ptype,vert2,vert,prev:

n := nops(C[1]):
r := nops(C):
ptype := [add(`if`(C[1][i]=1,1,0),i=1..n),seq(add(`if`(C[j][i]=1,1,0),i=1..n)-add(`if`(C[j-1][i]=1,1,0),i=1..n),j=2..r)]:
C2 := C:
p := [0$n]:
vert := [{}$(r-1)]:
prev := {}:
for i from 1 to r do
    for j from 1 to ptype[i] do
        for pos from 1 while C2[i][pos]=0 do od:
	prev := {op(prev),[i,pos]}:
    	for k from i to r-1 do
    	    (pos2,vert2) := Findnext2(C2,k,pos,prev):
	    prev := {op(prev),[k+1,pos2]}:
	    C2[k,pos]:=0:
	    pos := pos2:
	    vert[i] := {op(vert[i]),op(vert2)}:
	od:
	C2[r,pos]:=0:
    	p[pos]:=i:
    od:
od:

#print(prev):

#p[seq(`if`(p[i]=0,i,op({})),i=1..n)]:=n:
return p,vert:

end:

#Findnext2(C,m,k,prev) finds the 1 in row m+1 corresponding to the 1 in C[m,k]
Findnext2 := proc(C,m,k,prev) local i,j,n,vert:

if C[m,k]=0 then
   ERROR(`C[m,k] is not equal to 1`):
fi:

n := nops(C[1]):
vert := {}:
for i from k to n while C[m+1,i]=0 or member([m+1,i],prev) do 
    if not member([m+1,i],prev) then
        vert:={op(vert),[m+1,i]}:
    fi:
od:
if i=n+1 then
   for i from 1 to k-1 while C[m+1,i]=0 or member([m+1,i],prev) do 
       if not member([m+1,i],prev) then
       	  vert:={op(vert),[m+1,i]}:
       fi:
   od:
fi:
return i,vert:

end:

#MultiPermtoMultiline(P) given a permutation returns all the multiline configs
MultiPermtoMultiline := proc(P) local i,j,n,M,x:

n:=convert(P,multiset):
n:=[seq(n[i][2],i=1..nops(n))]:
M := MultilineVert(n):
return {seq(`if`(MultilinetoMultiPerm(M[i])[1]=P,M[i],op({})),i=1..nops(M))}:

end:

#Clockseq2(C,k) returns the list of positions for the clock to strike on the multiline
#starting with position k on the lowest
Clockseq2 := proc(C,k) local i,j,n,r,S:

r := nops(C):
n := nops(C[1]):
S := [0$(r-1),k]:
j := k:
for i from r to 2 by -1 do
    if C[i,j]=1 then
       S[i-1]:=j:
    else
       S[i-1]:=`if`(j<n,j+1,1):
    fi:
    j := S[i-1]:
od:

return S:

end:

#MultilineVert(T) returns the multiline process configs with T[1] 1's, etc, 
MultilineVert := proc(T) local i,j,k,n,c,C,configs:

n := add(T[i],i=1..nops(T)):
C:=[0$(nops(T))]:
for i from 1 to nops(T) do
    j := add(T[k],k=1..i):
    C[i] := combinat[permute]([1$j,0$(n-j)]):
od:

j := max(seq(nops(C[i]),i=1..nops(C))):
configs := combinat[permute]([seq(i$(nops(T)),i=1..j)],nops(C)):
configs := {op(configs)} minus {seq(seq(`if`(c[i]>nops(C[i]),c,op({})),i=1..nops(C)),c in configs)}:
configs := {seq([seq(C[k][configs[i][k]],k=1..nops(C))],i=1..nops(configs))}:
return configs:

end:

#MultilineGraph(T,x) returns the multiline process configs and transition rates
MultilineGraph := proc(T,x) local i,j,k,c,C,configs,Ed,cseq,cseqm1,n,r,val,val2:

n := add(T[i],i=1..nops(T)):
r := nops(T):
configs := MultilineVert(T):
Ed:=[]:
for c in configs do
    for i from 1 to n do
    	cseq := Clockseq2(c,i):
	cseqm1 := [seq(`if`(cseq[j]>1,cseq[j]-1,n),j=1..r)]:
	val := MultilinetoMultiPerm(c):
	if val[1][i]=r then
#	   val2 := {seq(op(val[2][j]),j=1..r-1)} intersect {seq([j,i],j=1..r-1)}:
#	   val2 := max(seq(
	   if {seq(op(val[2][j]),j=1..r-1)} intersect {seq([j,i],j=1..r-1)} ={} then
	      val := r-1:
	   else
	      val := r-2:
	      
	   fi:
	else
	   val := val[1][i]:
	fi:
	Ed := [op(Ed), [x[val],
			c,[seq(`if`(c[j][cseqm1[j]]=0 and c[j][cseq[j]]=1, 
				[seq(`if`(k=cseqm1[j],1,`if`(k=cseq[j],0,c[j][k])),k=1..n)],c[j]),j=1..r)]]]:
    od:
od:

return configs,Ed:

end:

#MultilineMatrix2(T,x) returns the transitions of the Markov chain
MultilineMatrix2 := proc(T,x) local i,j,k,M,Ve,Ed:

(Ve,Ed) := MultilineGraph(T,x):
M := Array(1..nops(Ve),1..nops(Ve)):
for i from 1 to nops(Ed) do
    for j from 1 to nops(Ve) do
        for k from 1 to nops(Ve) do
    	    if Ed[i][2] = Ve[j] and Ed[i][3] = Ve[k] then
	       M[k,j] := M[k,j] + Ed[i][1]:
	    fi:
	od:
    od:
od:

return Matrix(M):

end:

#MultilineMatrix(T,x) returns the transition matrix
MultilineMatrix := proc(T,x) local i,j,M:

M:=MultilineMatrix2(T,x):
for i from 1 to op(1,M)[1] do
    M[i,i] := -add(M[j,i],j=1..i-1)-add(M[j,i],j=i+1..op(1,M)[1]):
od:

return M:

end:

#rotate2(config) returns the config translated by 1 step assuming it is on the ring
rotate2 := proc(config) local i,c,n:
n:=nops(config[1]):
return [seq([c[n],seq(c[i],i=1..n-1)],c in config)]:
end:


#MultilineSteadyState(T,x,p) returns the steady state
MultilineSteadyState := proc(T,x,p) local i,j,M,V,V2,eqs,vars,soln,n,N:

n := nops(T):
N := add(T[i],i=1..n):
M := MultilineMatrix(T,x):
V2 := MultilineVert(T,x):
V := Vector([seq(p[op(V2[i])],i=1..nops(V2))]):
vars := convert(V,set):
eqs := convert(M.V,set) union {vars[1]=mul(x[i]^add(add(`if`(V2[1][j,k]=0,1,0),k=1..N),j=i+1..n),i=1..n)}:
eqs := eqs union {seq(p[op(V2[i])]=p[op(rotate2(V2[i]))],i=1..nops(V2))}:

soln := solve(eqs,vars):
return {seq(op(1,soln[i])=factor(op(2,soln[i])),i=1..nops(soln))}:

end:


#conj2SteadyState(T,x,p) returns the conjectured steady state for rates [x,1,...,1]
conj2SteadyState := proc(T,x,p) local i,j,r,V,M,S:

r := nops(T):
V := MultilineVertR(T,x):
S := {seq(p[op(V[i])]=x^(binomial(r-1,2)-nops(MultilinetoMultiPerm(V[i])[2][1])),i=1..nops(V))}:
#return S:

S := Vector([seq(x^(binomial(r-1,2)-nops(MultilinetoMultiPerm(V[i])[2][1])),i=1..nops(V))]):
M := MultilineMatrix(T,[x,1$(r-1)]):
#return convert(simplify(M.S),set):
return simplify(M.S):

end:

#wtmulline(L) guesses the stationary weight of the multiline config
wtmulline := proc(L,x) local i,j,r,P,vert,wt:

r := nops(L):
wt := mul(x[i]^add(add(`if`(L[j,k]=0,1,0),k=1..nops(L[1])),j=i+1..r),i=1..r):
(P,vert) := MultilinetoMultiPerm(L):
for i from 1 to r do
    wt := wt*mul((x[j]/x[i])^(add(`if`(vert[i][k][1]=j,1,0),k=1..nops(vert[i]))),j=i+1..r):
od:

return wt:

end:

#wtmulperm(L) guesses the stationary weight of the multiperm using above
wtmulperm := proc(P,x) local i,j,L,vert,wt:

L := MultiPermtoMultiline(P):
return simplify(add(wtmulline(i,x),i in L)):

end:

#################################################Multiline Coupe Process###############################################
#CoupeGraph(T,x) returns the multiline coupe process configs and transition rates
#Works only for 3 species
CoupeGraph := proc(T,x) local i,j,k,c,c2,C,configs,Ed,cseq,cseqm1,n,r,val,val2:

n := add(T[i],i=1..nops(T)):
r := nops(T):
if r <> 3 then
   ERROR(`Only three species`):
fi:

configs := MultilineVert(T):
Ed:=[]:
for c in configs do
    for i from 1 to n do
    	c2 := CoupeTrans(c,i):
	val := MultilinetoMultiPerm(c)[1][i]:
	if c2<> c then
		Ed := [op(Ed), [x[val],c,c2]]:
	fi:
    od:
od:

return configs,Ed:

end:

#CoupeTrans(C,k) returns the new multiline coupe by making a coupe transition
#starting with position k on the lowest
CoupeTrans := proc(C,k) local i,j,n,r,S,pi,cmid,cmidp1,cend,kp1:

pi := MultilinetoMultiPerm(C)[1]:
n := nops(pi):
r := nops(C):
if r <> 3 then
   ERROR(`Only three species`):
fi:

S := C:
if pi[`if`(k>1,k-1,n)]=1 or pi[k]=3 or (pi[k]=2 and pi[`if`(k>1,k-1,n)]=2) then
   return C:
fi:

if S[2,`if`(k>1,k-1,n)] = 0 then
   S[2,k]:=0:
   S[2,`if`(k>1,k-1,n)] := 1:
fi:


#Regular jump
if S[1,k] = 1 then 
   if S[1,`if`(k>1,k-1,n)] = 0 then
      S[1,k]:=0:
      S[1,`if`(k>1,k-1,n)] := 1:
   fi:
#fi:
#Pulling jump
elif S[1,k] = 0 then 
   for cend from k to n while pi[cend]=pi[k] do od:
   if cend=n+1 and pi[1]=pi[k] then
      for cend from 1 to n while pi[cend]=pi[k] do od:
   fi:
   cend := cend-1:

   #Front seat is not the back seat
   if cend <> k then
      if cend>k and k<n then
      	 for i from k+1 to cend do
      	     if S[1,i]=1 then
	     	S[1,i]:=0:
		S[1,i-1]:=1:
	     fi:		
	 od:
      else
      	 for i from k+1 to n do
      	     if S[1,i]=1 then
	     	S[1,i]:=0:
		S[1,i-1]:=1:
	     fi:		
	 od:
      	     if S[1,1]=1 then
	     	S[1,1]:=0:
		S[1,n]:=1:
	     fi:			 
      	 for i from 2 to cend do
      	     if S[1,i]=1 then
	     	S[1,i]:=0:
		S[1,i-1]:=1:
	     fi:		
	 od:
      fi:

   #Front seat is same as the back seat
   else
	#Next coupe
	kp1 := `if`(k<n,k+1,1): 
	if pi[kp1]=1 or pi[kp1]=2 then
	   for cend from kp1 to n while pi[cend]=pi[kp1] do od:
   	   if cend=n+1 and pi[1]=pi[kp1] then
      	   	for cend from 1 to n while pi[cend]=pi[kp1] do od:
   	   fi:
   	   cend := `if`(cend=1,n,cend-1):

	else
#         print(`here`):
	   for cend from kp1 to n while pi[cend]=pi[kp1] do od:
      	   if cend=n+1 and pi[1]=pi[kp1] then
      	      for cend from 1 to n while pi[cend]=pi[k+1] do od:
      	   fi:
	   cmid := cend-1:
	   cmidp1 := `if`(cmid<n,cmid+1,1):
	   for cend from cmidp1 to n while pi[cend]=pi[cmidp1] do od:
      	   if cend=n+1 and pi[1]=pi[cmidp1] then
      	      for cend from 1 to n while pi[cend]=pi[cmid+1] do od:
      	   fi:
	   cend := cend-1:
	fi:

      if cend=kp1  then
      	 if S[1,kp1]=1 then
	    S[1,kp1]:=0:
	    S[1,k]:=1:
	 fi:
      elif cend>kp1 and k<n-1 then
      	 for i from kp1 to cend do
      	     if S[1,i]=1 then
	     	S[1,i]:=0:
		S[1,i-1]:=1:
	     fi:		
	 od:
      else
      	 for i from k+1 to n do
      	     if S[1,i]=1 then
	     	S[1,i]:=0:
		S[1,i-1]:=1:
	     fi:		
	 od:
      	     if S[1,1]=1 then
	     	S[1,1]:=0:
		S[1,n]:=1:
	     fi:			 
      	 for i from 2 to cend do
      	     if S[1,i]=1 then
	     	S[1,i]:=0:
		S[1,i-1]:=1:
	     fi:		
	 od:
      fi:

  fi:
fi:

return S:

end:


#CoupeMatrix2(T,x) returns the transitions of the Markov chain
CoupeMatrix2 := proc(T,x) local i,j,k,M,Ve,Ed:

(Ve,Ed) := CoupeGraph(T,x):
M := Array(1..nops(Ve),1..nops(Ve)):
for i from 1 to nops(Ed) do
    for j from 1 to nops(Ve) do
        for k from 1 to nops(Ve) do
    	    if Ed[i][2] = Ve[j] and Ed[i][3] = Ve[k] then
	       M[k,j] := M[k,j] + Ed[i][1]:
	    fi:
	od:
    od:
od:

return Matrix(M):

end:

#CoupeMatrix(T,x) returns the transition matrix
CoupeMatrix := proc(T,x) local i,j,M:

M:=CoupeMatrix2(T,x):
for i from 1 to op(1,M)[1] do
    M[i,i] := -add(M[j,i],j=1..i-1)-add(M[j,i],j=i+1..op(1,M)[1]):
od:

return M:

end:


#CoupeSteadyState(T,x,p) returns the steady state
CoupeSteadyState := proc(T,x,p) local i,j,M,V,V2,eqs,vars,soln:

M := CoupeMatrix(T,x):
V2 := MultilineVert(T,x):
V := Vector([seq(p[op(V2[i])],i=1..nops(V2))]):
vars := convert(V,set):
eqs := convert(M.V,set) union {vars[1]=x[1]^T[3]}:
#mul(x[i]^binomial(n-i,2),i=1..n-1)}:
eqs := eqs union {seq(p[op(V2[i])]=p[op(rotate2(V2[i]))],i=1..nops(V2))}:

soln := solve(eqs,vars):
return {seq(op(1,soln[i])=factor(op(2,soln[i])),i=1..nops(soln))}:

end:


#################################################Schubert Polynomials###############################################
#divdiff(f,x,i) returns the divided difference of the function f(x_1,...) with respect to x[i]
divdiff := proc(f,x,i)

return factor((f-subs({x[i]=x[i+1],x[i+1]=x[i]},f))/(x[i]-x[i+1])):

end:

#homsym(n,d,x) returns the complete homogenous symmetric polynomial in variables x[1]..x[n] of degree d
homsym := proc(n,d,x) local i,j,S:

if d<0 then
	return 0:
fi:

S:= combinat[choose]([seq(x[i]$d,i=1..n)],d):
return add(mul(S[i][j],j=1..d),i=1..nops(S)):

end: