######################################################################
##IsingPolygons.txt: Save this file as IsingPolygons.txt             #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
#and then type:  read IsingPolygons.txt<Enter>                       #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created: Feb. 2018
 
print(`Created:  Feb. 2018`):
print(` This is IsingPolygons.txt `):
print(`To count (and weight)-count lattice graphs (collection of edges) such that every vertex a degree in a prescribed set`):
print(`It is one of the packages that accompany the article `):
print(`Y`):
print(`by Manuel Kauers and Doron Zeilberger`):

print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to DoronZeil at gmail dot com `):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://sites.math.rutgers.edu/~zeilberg/  .`):

print(`---------------------------------------`):
 print(`For a list of the Checking procedures type ezraC();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the CYLINDER procedures type ezraCy();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):

ezraCy:=proc()

if args=NULL then
 print(` The CYLINDER procedures are: FollowersCy, FollowersEcy, GFnstCy, IsFoFCy,  IsFoCy, LeftLettersCy, RanWordCy, VocabCy `):
 print(``):

else
ezra(args):
fi:

end:

ezraC:=proc()

if args=NULL then
 print(` The Checking procedures are: CheckGFnst, CheckGFnstCy `):
 print(``):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: Alphabet, Conts, ContsCy, ContsE, ContsEcy, CZmn, DegWt, Followers, FollowersE, IsFo, IsFoF, LeftLetters, NuToLet `):
 print(`TM, TMz`):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: DrawWord, GFnst, GFnstTo, Vocab, RanWordCy`): 

elif nops([args])=1 and op(1,[args])=Alphabet then
print(`Alphabet(n): the alphabet of graphs in an torus of width n i.e. [0,n] times [0,k] for arbitray k. Try:`):
print(`Alphabet(3);`):

elif nops([args])=1 and op(1,[args])=CheckGFnst then
print(`CheckGFnst(n,k0, DegS): Checks GFnst(n,s,t,DegS) up to the k0-coefficient of t`):
print(`Try:`):
print(`CheckGFnst(2,4,{0,2,4});`):

elif nops([args])=1 and op(1,[args])=CheckGFnstCy then
print(`CheckGFnstCy(n,k0, DegS): Checks GFnstCy(n,s,t,DegS) up to the k0-coefficient of t`):
print(`Try:`):
print(`CheckGFnstCy(2,4,{0,2,4});`):

elif nops([args])=1 and op(1,[args])=Conts then
print(`Conts(w,n,DegS): Given a word w in the alphabet of  graphs in a strip of width n, where every vertex has a degree in`):
print(`the set DegS. Outputs all its legal continuations. Try:`):
print(`Conts([[{1,2},{0,2}]],2,{0,2,4});`):

elif nops([args])=1 and op(1,[args])=ContsCy then
print(`ContsCy(w,n,DegS): Given a word w in the alphabet of  degS-degree graphs in a CYLINDER of width n, outputs all its legal continuations. Try:`):
print(`ContsCy([[{1,2},{0,2}]],2,{0,2,4});`):

elif nops([args])=1 and op(1,[args])=ContsE then
print(`ContsE(w,n,DegS): Given a word w in the alphabet of  degS-degree polygons in a strip of width n, outputs all its legal continations that end the word. Try:`):
print(`ContsE([[{1,2},{0,2}]],2,{0,2,4});`):

elif nops([args])=1 and op(1,[args])=ContsEcy then
print(`ContsEcy(w,n): Given a word w in the alphabet of  even-degree polygons in a CYLINDER of width n, outputs all its legal continations that end the word. Try:`):
print(`ContsEcy([[{1,2},{0,2}]],2);`):


elif nops([args])=1 and op(1,[args])=CZmn then
print(`CZmn(m,n,v): the polynomial Z_{m,n}(v) done via the combinatorial method. Try: `):
print(`CZmn(2,2,v);`):

elif nops([args])=1 and op(1,[args])=DegWt then
print(`DegWt(S1,S2,S3,n,z): the weight of the transiton S1,S2,S3 where the weight is mul(z[i]#verticesOfDegree i,i=0..4). Try: `):
print(` DegWt({0,1},{0,1},{0,1},2,z); `):

elif nops([args])=1 and op(1,[args])=DrawWord then
print(`DrawWord(w): draws the word w. Try:`):
print(`DrawWord([[{2},{1,2}], [{1,2},{0,2}],{1}]); `):

elif nops([args])=1 and op(1,[args])=Followers then
print(`Followers(P,n,DegS): inputs a letter P=[S1,S2] in the DegS-degree graphs  alphabet for strips of width n,`):
print(` and outputs the set of its legal followers. Try:`):
print(`Followers([{1,2},{0,2}],2,{0,2,4});`):

elif nops([args])=1 and op(1,[args])=FollowersCy then
print(`FollowersCy(P,n,DegS): inputs a letter P=[S1,S2] in the DegS-graphs  alphabet for a CYLINDER of width n,`):
print(` and outputs the set of its legal followers. Try:`):
print(`FollowersCy([{1,2},{0,2}],2,{0,2,4});`):

elif nops([args])=1 and op(1,[args])=FollowersE then
print(`FollowersE(P,n,DegS): inputs a letter P=[S1,S2] in the DegS-graph alphabet for strips of width n, and outputs the set of "half letters"`):
print(`i.e. subsets of {1, ...,n} that end the word.`):
print(`Try: FollowersE([{1,2},{0,2}],2,{0,2,4});`):

elif nops([args])=1 and op(1,[args])=FollowersEcy then
print(`FollowersEcy(P,n,DegS): inputs a letter P=[S1,S2] in the DegS-graphs  alphabet for CYLINDERS of width n, and outputs the set of "half letters"`):
print(`i.e. subsets of {0, ...,n} that end the word.`):
print(`Try: FollowersEcy([{1,2},{0,2}],2,{0,2,4});`):


elif nops([args])=1 and op(1,[args])=GFnst then
print(`GFnst(n,s,t,DegS): The generating function enumerating DegS-graphs contained in a strip of length n with weight t^(length)*s^(number of edges)`):
print(`Try:`):
print(`GFnst(2,s,t,{0,2,4});`):

elif nops([args])=1 and op(1,[args])=GFnstCy then
print(`GFnstCy(n,s,t,DegS): The generating function enumerating even-degree polygons contained in a CYLINDER of length n`):
print(` with weight t^(length)*s^(number of edges)`):
print(`Try:`):
print(`GFnstCy(2,s,t,{0,2,4});`):

elif nops([args])=1 and op(1,[args])=GFnstTo then
print(`GFnstTo(n,s,t,DegS): the generating function in the variables s,t, for polygons where each vertex has degree`):
print(` belonging to the set DegS, whose coefficient of t^m*s^k is the number of such polygons bounded in an m by n torodial strip with`):
print(` exactly k edges. Warning: Very slow. On my computer I can only do n=0,1,2. Try:`):
print(`GFnstTo(1,s,t,{0,2,4});`):

elif nops([args])=1 and op(1,[args])=GFnstToZ then
print(`GFnstToZ(n,s,t,z): the generating function in the variables s,t, and z[0],..., z[4], for all graphs`):
print(` bounded in [0,n] by [0,m] whose coefficient of t^m*s^k*z[0]^d0*...z[4]^d4 is the number of graphs bounded in an m by n torodial strip with`):
print(` exactly k edges,  d0 vertices of degree 0, d1, vertices of degree 1, ..., d4 vertices of degree 4. `):
print(` Extremely slow. On my computer I can only do n=0,1,2. Try:`):
print(`GFnstToZ(1,s,t,z);`):

elif nops([args])=1 and op(1,[args])=IsFo then
print(`IsFo(S1,S2,S3,n,DegS): inputs subsets S1, S2, S3 of {0,1,...,n}, {1,..,n}, and {0,1,....n}, pos. integer n, and a `):
print(` a set of integers, DegS, (a subset of {0,1,2,3,4}, where `):
print(`S1 is the set of horizontal edges where i represents a horizontal edge from [a-1,i] to [a,i], `):
print(`S2 is the set of vertical edges along a generic vertical line x=a, where i represents a vertical edge from [a,i-1] to [a,i], `):
print(`S3 is the sets of horizontal edges where i represents a horizontal edge from [a,i] to [a+1,i], `):
print(`decides whether S2,S3 can follow after S1 by checking that every vertex on the line x=a has a degree that belongs to the set DegS`):
print(`(i.e. 0 or 2 or 4, and never 1 or 3). Try: `):
print(` IsFo({0,3},{3},{0,2},3,{0,2,4}); `):

elif nops([args])=1 and op(1,[args])=IsFoCy then
print(`IsFoCy(S1,S2,S3,n,DegS):  inputs subsets S1, S2, S3 of {0,1,...,n}, {1,..,n}, and {0,1,....n}, pos. integer n, and a `):
print(` a set of integers, DegS, (a subset of {0,1,2,3,4}, where `):
print(`S1 is the set of horizontal edges where i represents a horizontal edge from [a-1,i] to [a,i], `):
print(`S2 is the set of vertical edges along a generic vertical line x=a, where i represents a vertical edge from [a,i-1] to [a,i], `):
print(`(where the edge connecting [a,0] and [a,n] is denoted by 0)`):
print(`S3 is the set of horizontal edges where i represents a horizontal edge from [a,i] to [a+1,i], `):
print(`decides whether S2,S3 can follow after S1 by checking that every vertex on the line x=a has a degree in the set DegS`):
print(`(i.e. 0 or 2 or 4, and never 1 or 3). Try: `):
print(` IsFoCy({0,3},{3},{0,2},3,{0,2,4}); `):


elif nops([args])=1 and op(1,[args])=IsFoF then
print(`IsFoF(S1,S2,n,DegS): inputs subsets S1 and S2 of {1,...,n} and {0,1,..,n} and a subset of {0,1,2,3,4}, DegS.`):
print(`S1 is the set of vertical edges along the y-axis, in left-most border of an n by infinity strip and S2 is the set of horizontal edges`):
print(`between x=0 and x=1, decides whether S1 can follow S2 by checking that every vertex on the y-axis has even degree`):
print(` (i.e. 0 or 2 or 4, and never 1 or 3). Try: `):
print(` IsFoF({1,2,3,4},{0,4},4,{0,2,4}); `):
print(`IsFoF({1,2,3,4},{0},5,{0,2,4}); `):

elif nops([args])=1 and op(1,[args])=IsFoFCy then
print(`IsFoFCy(S1,S2,n,DegS):  For the CYLINDER. `):
print(`inputs subsets S1 and S2 of {1,...,n} and {0,1,..,n} and a subset of {0,1,2,3,4}, DegS.`):
print(`S1 is the set of vertical edges along the y-axis, in left-most border of an n by infinity strip and S2 is the set of horizontal edges`):
print(`between x=0 and x=1, decides whether S1 can follow S2 by checking that every vertex on the y-axis has a degree that belongs to DegS`):
print(` (i.e. 0 or 2 or 4, and never 1 or 3). Try: `):
print(` IsFoF({1,2,3,4},{0,4},4,{0,2,4}); `):
print(`IsFoFCy({1,2,3,4},{0},5,{0,2,4}); `):


elif nops([args])=1 and op(1,[args])=LeftLetters then
print(`LeftLetters(n,DegS): the set of pairs [S1,S2] where S1 is a subset of {1,....,n} and S2 is a subset of {0,...,n} denting`):
print(`sets of vertical and horizontal edges respectively such that IsFoF(S1,S2,DegS) (q.v.) is true. Try: `):
print(`LeftLetters(3,{0,2,4});`):

elif nops([args])=1 and op(1,[args])=LeftLettersCy then
print(`LeftLettersCy(n,DegS): For the CYLINDER. the set of pairs [S1,S2] where S1 and S2 are subsets of {0,...,n} denting`):
print(`sets of vertical and horizontal edges respectively (the edge connection y=0 and y=n is denoted by 0) such that IsFoFCy(S1,S2,DegS) (q.v.) is true. Try: `):
print(`LeftLettersCy(3,{0,2,4});`):


elif nops([args])=1 and op(1,[args])=NuToLet then
print(`NuToLet(k,n): inputs pos. integers k and n such that 1<=k<=2^n and outputs the letter [V,H] corresponding to it`):
print(`by converting k to binary and letting  the first n bits correpond to V and the last n bits to H. Try:`):
print(`NuToLet(3,2);`):

elif nops([args])=1 and op(1,[args])=RanWord then
print(`RanWord(k,n,DegS): a random word in the STRIP alphabet of length k in the [0,n] strip. Try:`):
print(`RanWord(10,4,{0,2,4});`):

elif nops([args])=1 and op(1,[args])=RanWordCy then
print(`RanWordCy(k,n,DegS): a random word of length k in the [0,n] in the CYLINDER alphabet. Try:`):
print(`RanWordCy(10,4,{0,2,4});`):

elif nops([args])=1 and op(1,[args])=TM then
print(`TM(n,DegS,s): the transfer matrix for graphs in the TORUS with width n [0,n], followed by the. Try:`):
print(`TM(2,{0,2,4},s);`):

elif nops([args])=1 and op(1,[args])=TMz then
print(`TMz(n,DegS,s,z): the transfer matrix for graphs in the TORUS according to s^#NumberOfEdges*degreeCount`):
print(`TMz(2,s,z);`):

elif nops([args])=1 and op(1,[args])=Vocab then
print(`Vocab(k,n,DegS): the set of words for DegS-degree graphs, bounded in an n by k rectangle. Try: `):
print(`Vocab(4,2,{0,2,4});`):

elif nops([args])=1 and op(1,[args])=VocabCy then
print(`VocabCy(k,n, DegS): the set of words for DegS-graphs bounded in an n by k CYLINDER. Try: `):
print(`VocabCy(4,2,{0,2,4});`):

else
print(`There is no ezra for`,args):
fi:
 
end:



emet:=proc(B)
if B then
 1:
else
0:
fi:
end:


IsFoF:=proc(S1,S2,n,DegS) local i,d:

if not (type(S1,set) and  type(S2,set) and  S1 minus {seq(i,i=1..n)}={} and S2 minus {seq(i,i=0..n)}={} ) then
  RETURN(FAIL):
fi:

d:=emet(member(1,S1))+emet(member(0,S2)):

if not member(d,DegS) then
  RETURN(false):
fi:

for i from 1 to n-1 do
 d:=emet(member(i,S1))+emet(member(i+1,S1))+emet(member(i,S2)):

 if   not member(d,DegS) then
    RETURN(false):
 fi:

od:

d:=emet(member(n,S1))+emet(member(n,S2)):

 if not member(d,DegS) then
    RETURN(false):
 fi:

true:

end:




IsFo:=proc(S1,S2,S3,n,DegS) local i,d:

if not (type(S1,set) and  type(S2,set) and type(S3,set) and  S1 minus {seq(i,i=0..n)}={} and S2 minus {seq(i,i=1..n)}={} 
and  S3 minus {seq(i,i=0..n)}={}) then
  RETURN(FAIL):
fi:

d:=emet(member(0,S1))+emet(member(1,S2))+ emet(member(0,S3)):

if not member(d,DegS)  then
  RETURN(false):
fi:

for i from 1 to n-1 do
 d:=emet(member(i,S1))+emet(member(i,S2))+emet(member(i+1,S2)) +emet(member(i,S3)):

 if not member(d,DegS)  then
    RETURN(false):
 fi:

od:

 d:=emet(member(n,S1))+emet(member(n,S2)) +emet(member(n,S3)):

 if not member(d,DegS)  then
    RETURN(false):
 fi:

true:

end:


LeftLetters:=proc(n,DegS) local gu1,gu2,gu,S1,S2,i:

gu1:=powerset({seq(i,i=1..n)}): gu2:=powerset({seq(i,i=0..n)}):
gu:={}:

for S1 in gu1 do
 for S2 in gu2 do
  if IsFoF(S1,S2,n,DegS) then
   gu:=gu union {[S1,S2]}:
  fi:
 od:
od:

gu:

end:

Followers:=proc(P,n,DegS)  local gu1,gu2,gu,S1,S2,i:

gu1:=powerset({seq(i,i=1..n)}): gu2:=powerset({seq(i,i=0..n)}):
gu:={}:

for S1 in gu1 do
 for S2 in gu2 do
  if IsFo(P[2],S1,S2,n,DegS) then
   gu:=gu union {[S1,S2]}:
  fi:
 od:
od:

gu:

end:



FollowersE:=proc(P,n,DegS)  local gu1,gu,S1,i:

gu1:=powerset({seq(i,i=1..n)}): 
gu:={}:

for S1 in gu1 do
  if IsFoF(S1,P[2],n,DegS) then
   gu:=gu union {S1}:
  fi:
od:

gu:

end:



Conts:=proc(w,n,DegS) local gu,gu1:

if w=[] then
  gu:=LeftLetters(n,DegS):
  RETURN({seq([gu1],gu1 in gu)}):
fi:

gu:={}:


gu:=Followers(w[-1],n,DegS):

{seq([op(w),gu1],gu1 in gu)}:

end:


ContsE:=proc(w,n,DegS) local gu,gu1:

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

gu:={}:


gu:=FollowersE(w[-1],n,DegS):

{seq([op(w),gu1],gu1 in gu)}:

end:



Vocab:=proc(k,n,DegS) local gu,gu1,i:
if k=0 then
 RETURN({[]}):
fi:

gu:=LeftLetters(n,DegS):

gu:={seq([gu1],gu1 in gu)}:

for i from 2 to k do
gu:={seq(op(Conts(gu1,n,DegS)),gu1 in gu)}:
od:

{seq(op(ContsE(gu1,n,DegS)),gu1 in gu)}:
end:

with(plots):



DrawWord:=proc(w) local pic,i,mu,mu1:

if w=[] then
 RETURN(FAIL):
fi:

if not type(w[-1],set) then
 RETURN(FAIL):
fi:


pic:=plot({[0,0]},style=point,symbol=CIRCLE, scaling=constrained,axes=none):

for i from 1 to nops(w)-1 do
 mu:=w[i][1] minus {0}:

 for mu1 in mu do
  pic:=pic, plot([[i-1,mu1-1],[i-1,mu1]]):
 od:

  mu:=w[i][2]:

  for mu1 in mu do
  pic:=pic, plot([[i-1,mu1],[i,mu1]]):
 od:

od:

mu:=w[nops(w)] minus {0}:

 for mu1 in mu do
  pic:=pic, plot([[nops(w)-1,mu1-1],[nops(w)-1,mu1]]):
 od:

display(pic);

end:




RanWord:=proc(k,n,DegS) local gu,gu1,w,i:
gu:=LeftLetters(n,DegS):

gu:={seq([gu1],gu1 in gu)}:

if gu={} then
  RETURN(FAIL):
fi:

w:=gu[rand(1..nops(gu))()]:

for i from 2 to k do
gu:=Conts(w,n,DegS):
if gu={} then
  RETURN(FAIL):
fi:

w:=gu[rand(1..nops(gu))()]:

od:

gu:=ContsE(w,n,DegS):
if gu={} then
  RETURN(FAIL):
fi:

gu[rand(1..nops(gu))()]:
end:


GFnst:=proc(n,s,t,DegS) local guL,guM,L,A,fu,eq1,L1,eq,var:
option remember:

guL:=LeftLetters(n,DegS):

guM:={seq(op(Followers(L,n,DegS)),L in guL)}:

eq:={}:
var:={}:

for L in guM do
 var:=var union {A[L]}:
 fu:=FollowersE(L,n,DegS):
 eq1:=A[L]-add(t*s^(nops(L[1])+nops(L[2])+nops(fu1)),fu1 in fu):
 fu:=Followers(L,n,DegS):
  eq1:=eq1-t*s^(nops(L[1])+nops(L[2]))*add(A[L1],L1 in fu):
 eq:=eq union {eq1}:
od:

var:=solve(eq,var):

normal(1+add(subs(var,A[L]),L in guL)):

end:




CheckGFnst1:=proc(n,k0,DegS) local gu,s,t,gu1,i,mu,lu:
gu:=Vocab(k0,n,DegS):

lu:=add(s^(add(nops(gu1[i][1])+nops(gu1[i][2]),i=1..nops(gu1)-1)+nops(gu1[nops(gu1)])),gu1 in gu):

mu:=GFnst(n,s,t,DegS):

mu:=taylor(mu,t=0,k0+1):
expand(coeff(mu,t,k0)-lu):
end:

CheckGFnst:=proc(n,k0,DegS) local k:
evalb({seq(CheckGFnst1(n,k,DegS),k=1..k0)}={0}):
end:




###start Cylinder





IsFoFCy:=proc(S1,S2,n,DegS) local i,d:

if not (type(S1,set) and  type(S2,set) and  S1 minus {seq(i,i=0..n)}={} and S2 minus {seq(i,i=0..n)}={} ) then
  RETURN(FAIL):
fi:



for i from 0 to n-1 do
 d:=emet(member(i,S1))+emet(member(i+1,S1))+emet(member(i,S2)):
 if not member(d,DegS) then
    RETURN(false):
 fi:

od:

d:=emet(member(n,S1))+ emet(member(0,S1))+emet(member(n,S2)):

 if not member(d,DegS) then
    RETURN(false):
 fi:

true:

end:


LeftLettersCy:=proc(n,DegS) local gu1,gu2,gu,S1,S2,i:

gu1:=powerset({seq(i,i=0..n)}): gu2:=powerset({seq(i,i=0..n)}):
gu:={}:

for S1 in gu1 do
 for S2 in gu2 do
  if IsFoFCy(S1,S2,n,DegS) then
   gu:=gu union {[S1,S2]}:
  fi:
 od:
od:

gu:

end:


IsFoCy:=proc(S1,S2,S3,n,DegS) local i,d:

if not (type(S1,set) and  type(S2,set) and type(S3,set) and  S1 minus {seq(i,i=0..n)}={} and S2 minus {seq(i,i=0..n)}={} 
and  S3 minus {seq(i,i=0..n)}={}) then
  RETURN(FAIL):
fi:


for i  from 0 to n-1 do
 d:=emet(member(i,S1))+emet(member(i,S2))+emet(member(i+1,S2)) +emet(member(i,S3)) :
 if not member(d,DegS) then
    RETURN(false):
 fi:

od:

 d:=emet(member(n,S1))+emet(member(n,S2)) + emet(member(0,S2)) +emet(member(n,S3)):

 if not member(d,DegS) then
    RETURN(false):
 fi:

true:

end:


FollowersCy:=proc(P,n,DegS)  local gu1,gu2,gu,S1,S2,i:

gu1:=powerset({seq(i,i=0..n)}): gu2:=powerset({seq(i,i=0..n)}):
gu:={}:

for S1 in gu1 do
 for S2 in gu2 do
  if IsFoCy(P[2],S1,S2,n,DegS) then
   gu:=gu union {[S1,S2]}:
  fi:
 od:
od:

gu:

end:



FollowersEcy:=proc(P,n,DegS)  local gu1,gu,S1,i:

gu1:=powerset({seq(i,i=0..n)}): 
gu:={}:

for S1 in gu1 do
  if IsFoFCy(S1,P[2],n,DegS) then
   gu:=gu union {S1}:
  fi:
od:

gu:

end:


ContsCy:=proc(w,n,DegS) local gu,gu1:

if w=[] then
  gu:=LeftLettersCy(n,DegS):
  RETURN({seq([gu1],gu1 in gu)}):
fi:

gu:={}:


gu:=FollowersCy(w[-1],n,DegS):

{seq([op(w),gu1],gu1 in gu)}:

end:


ContsEcy:=proc(w,n,DegS) local gu,gu1:

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

gu:={}:


gu:=FollowersEcy(w[-1],n,DegS):

{seq([op(w),gu1],gu1 in gu)}:

end:


VocabCy:=proc(k,n,DegS) local gu,gu1,i:
if k=0 then
 RETURN({[]}):
fi:

gu:=LeftLettersCy(n,DegS):

gu:={seq([gu1],gu1 in gu)}:

for i from 2 to k do
gu:={seq(op(ContsCy(gu1,n,DegS)),gu1 in gu)}:
od:

{seq(op(ContsEcy(gu1,n,DegS)),gu1 in gu)}:
end:




RanWordCy:=proc(k,n,DegS) local gu,gu1,w,i:
gu:=LeftLettersCy(n,DegS):

gu:={seq([gu1],gu1 in gu)}:


if gu={} then
  RETURN(FAIL):
fi:

w:=gu[rand(1..nops(gu))()]:



for i from 2 to k do
gu:=ContsCy(w,n,DegS):

if gu={} then
  RETURN(FAIL):
fi:

w:=gu[rand(1..nops(gu))()]:



od:

gu:=ContsEcy(w,n,DegS):

if gu={} then
  RETURN(FAIL):
fi:


gu[rand(1..nops(gu))()]:
end:


GFnstCy:=proc(n,s,t,DegS) local guL,guM,L,A,fu,eq1,L1,eq,var:
option remember:

guL:=LeftLettersCy(n,DegS):

guM:={seq(op(FollowersCy(L,n,DegS)),L in guL)}:

eq:={}:
var:={}:

for L in guM do
 var:=var union {A[L]}:
 fu:=FollowersEcy(L,n,DegS):

 eq1:=A[L]-add(t*s^(nops(L[1])+nops(L[2])+nops(fu1)),fu1 in fu):
 fu:=FollowersCy(L,n,DegS):
  eq1:=eq1-t*s^(nops(L[1])+nops(L[2]))*add(A[L1],L1 in fu):
 eq:=eq union {eq1}:
od:


var:=solve(eq,var):


normal(1+add(subs(var,A[L]),L in guL)):

end:


CheckGFnstCy1:=proc(n,k0,DegS) local gu,s,t,gu1,i,mu,lu:
gu:=VocabCy(k0,n,DegS):

lu:=add(s^(add(nops(gu1[i][1])+nops(gu1[i][2]),i=1..nops(gu1)-1)+nops(gu1[nops(gu1)])),gu1 in gu):

mu:=GFnstCy(n,s,t,DegS):

mu:=taylor(mu,t=0,k0+1):
expand(coeff(mu,t,k0)-lu):
end:

CheckGFnstCy:=proc(n,k0,DegS) local k:
evalb({seq(CheckGFnstCy1(n,k,DegS),k=1..k0)}={0}):
end:




#Alphabet(n): the alphabet of graphs in an torus of width n i.e. [0,n] times [0,k] for arbitray k. Try:
#Alphabet(3);
Alphabet:=proc(n) local gu,gu1,gu2,i:

gu:=powerset({seq(i,i=0..n)}): 

{seq(seq([gu1,gu2],gu2 in gu),gu1 in gu)}:

end:


#VecToSet(v): inputs a binary vetor, outputs the corresponding suset of {1, ..., nops(v))}
VecToSet:=proc(v) local gu,i:
gu:={}:
for i from 1 to nops(v) do
 if v[i]=1 then
  gu:=gu union {i}:
 fi:
od:
gu:
end:

#NuToLet(k,n): inputs pos. integers k and n such that 1<=k<=2^(2*(n+1)) and outputs the letter [V,H] corresponding to it
#by converting k to binary and letting  the first n+1 bits correpond to V and the last n+1 bits to H. Try:
#NuToLet(3,2);
NuToLet:=proc(k,n) local gu,i,gu1,gu2:

if not (1<=k and k<=2^(2*n+2)) then
 RETURN(FAIL):
fi:

gu:=convert(k-1,base,2):

gu:=[seq(gu[nops(gu)+1-i],i=1..nops(gu))]:

gu:=[0$(2*n+2-nops(gu)),op(gu)]:
gu1:=[op(1..n+1,gu)]: gu2:=[op(n+2..2*n+2,gu)]:
gu1:=VecToSet(gu1):
gu2:=VecToSet(gu2):

gu1:=subs({seq(i=i-1,i=1..n+1)},gu1): gu2:=subs({seq(i=i-1,i=1..n+1)},gu2):
[gu1,gu2]:

end:




#TM(n,DegS,s): the transfer matrix for graphs in the TORUS with width n [0,n], followed by the. Try:
#TM(2,{0,2,4},s);
TM:=proc(n,DegS,s) local T,i,j,gu1,gu2:
option remember:
for i from 1 to 2^(2*n+2) do
  gu1:=NuToLet(i,n):
 for j from 1 to 2^(2*n+2) do
   gu2:=NuToLet(j,n):
    if IsFoCy(gu1[2],gu2[1],gu2[2],n,DegS) then
     T[i,j]:=s^(nops(gu1[1])+nops(gu1[2])):
    else
       T[i,j]:=0:
    fi:
 od:
od:

[seq([seq(T[i,j],j=1..2^(2*n+2))],i=1..2^(2*n+2))]:
   
end:

#GFnstTo(n,s,t,DegS): the generating function in the variables s,t, for polygons where each vertex has degree belonging to the set DegS, 
#whose coefficient of t^m*s^k is the number of such polygons bounded in an m by n torodial strip with exactly k edges. Try:
#GFnstTo(1,s,t,{0,2,4});
GFnstTo:=proc(n,s,t,DegS) local gu:
gu:=Matrix(TM(n,DegS,s)):
normal(LinearAlgebra[Trace](LinearAlgebra[MatrixInverse](1-t*gu))):
end:





#DegWt(S1,S2,S3,n,z): the weight of the transiton S1,S2,S3 where the weight is mul(z[i]#verticesOfDegree i,i=0..4). Try
#DegWt({0,1},{0,1},{0,1},2,z);
DegWt:=proc(S1,S2,S3,n,z) local i,d,gu:

if not (type(S1,set) and  type(S2,set) and type(S3,set) and  S1 minus {seq(i,i=0..n)}={} and S2 minus {seq(i,i=0..n)}={} 
and  S3 minus {seq(i,i=0..n)}={}) then
  RETURN(FAIL):
fi:


gu:=1:

for i  from 0 to n-1 do
 d:=emet(member(i,S1))+emet(member(i,S2))+emet(member(i+1,S2)) +emet(member(i,S3)) :
 gu:=gu*z[d]:

od:

 d:=emet(member(n,S1))+emet(member(n,S2)) + emet(member(0,S2)) +emet(member(n,S3)):

 gu:=gu*z[d]:

gu:

end:



#TMz(n,s,z): the transfer matrix for graphs in the TORUS according to s^#NumberOfEdges*degreeCount
#TMz(2,s,z);
TMz:=proc(n,s,z) local T,i,j,gu1,gu2:
option remember:
for i from 1 to 2^(2*n+2) do
  gu1:=NuToLet(i,n):
 for j from 1 to 2^(2*n+2) do
   gu2:=NuToLet(j,n):

     T[i,j]:=s^(nops(gu1[1])+nops(gu1[2]))*DegWt(gu1[2],gu2[1],gu2[2],n,z) :


 od:
od:

[seq([seq(T[i,j],j=1..2^(2*n+2))],i=1..2^(2*n+2))]:
   
end:


#GFnstToZ(n,s,t,z): the generating function in the variables s,t, z[0], ..., z[4], for all subgraphs of [0,n] times [0,m]
#whose coefficient of t^m*s^k*prod(z[degree(v)], v all vertices). Try:
#GFnstToZ(1,s,t,z);
GFnstToZ:=proc(n,s,t,z) local gu:
gu:=Matrix(TMz(n,s,z)):
normal(LinearAlgebra[Trace](LinearAlgebra[MatrixInverse](1-t*gu))):
end:

#CZmn(m,n,v): the polynomial Z_{m,n}(v) done via the combinatorial method. Try
#CZmn(2,2,v):
CZmn:=proc(m,n,v) local gu:
gu:=Matrix(TM(m-1,{0,2,4},v)):
expand(LinearAlgebra[Trace](gu^n)):
end: