#######################################################################
#Using the metropolis algorthim to search for Folner sets
#######################################################################
#
#
#
#Folner sets are finite subsets of an infinite graph with a very low
#ratio of outgoing edges to vertices.  The inf of this ratio over all
#finite vertex subsets is called the isoperimetry constant of the graph.
#If this constant is zero, we call the graph amenable.  A 30-year old
#open question is whether the Cayley graph of a certain group, Thompson's Group
#F, is amenable.
#
#The best known upper bound for the isoperimetry constant of F is 1/2.
#This code is designed to search a graph automatically for "good" Folner sets,
#perhaps even finding one in F with ratio less than 1/2.
#
#The basic plan is to have a subset of the Cayley graph which we keep track of,
#as well as the set of its Neighbors (vertices it touches but that aren't in the set).
#We then choose a random vertex either on the boundary of the set or one of its Neighbors,
#and add/subtract that vertex according to the Metropolis algorithm.
#
#S is usually the current set of vertices, N the set of neighbors, and
#B is a subset of S consisting of its "boundary", i.e., vertices of S which are adjacent to
#vertices of N
#
#This code tests our algorithm on the Cayley graph of ZxZ (which is just the
#2-d integer lattice).  We've had some troubles getting things together for F,
#but hope to run it soon.


print(`For help type Help()`):

Help := proc() print(`For background information, call HelpB().`):
print(`For procedure names and functions, call HelpP().`): end:

HelpB := proc() print(`Using the metropolis algorthim to search for Folner sets`):
print(`Folner sets are finite subsets of an infinite graph with a very low ratio of outgoing edges to vertices.  The inf of this ratio over all finite vertex subsets is called the isoperimetry constant of the graph. If this constant is zero, we call the graph amenable.  A 30-year old open question is whether the Cayley graph of a certain group, Thompson's Group F, is amenable.`):

print(`The best known upper bound for the isoperimetry constant of F is 1/2. This code is designed to search a graph automatically for "good" Folner sets, perhaps even finding one in F with ratio less than 1/2.`):

print(`The basic plan is to have a subset of the Cayley graph which we keep track of, as well as the set of its Neighbors (vertices it touches but that aren't in the set). We then choose a random vertex either on the boundary of the set or one of its Neighbors, and use the Metropolis algorithm to decide whether to add/subtract the vertex.`):

print(`S is usually the current set of vertices, N the set of neighbors, and B is a subset of S consisting of its "boundary", i.e., vertices of S which are adjacent to vertices of N`):

print(`This code tests our algorithm on the Cayley graph of ZxZ (which is just the 2-d integer lattice).  We've had some troubles getting things together for F, but hope to run it soon.`):

end:

HelpP := proc()

print(`To run the Metropolis algorithm on the 2-d integer lattice, start with a set S of pairs of integers.  Then call Metropolis(S,K,ZZNeighbors), filling in the number of iterations you want to run for K.  The ZZNeighbors function is a function that gives the neighbors of a vertex.  One could replace it with a different neighbors function to test other graph structures on the integer lattice, e.g., connecting diagonals.`):

print(`AdjustNeighbors(S,N,v,Neighbors): Adds/removes vertex v from set S (with neighbors N) and adjusts N accordingly (Neighbors is the Neighbors function).`):

print(`ZZNeighbors(v): Gives the set of neighbors of the vertex v in the square integer lattice`):

print(`SetNeighbors(S,Neighbors): Gives the set of neighbors of the set S in the square integer lattice`):

print(`SNplot(S,N): Plots two sets of points (intended to be a set and its neighbors) in red and green.`):

print(`Weight(S,N,Neighbrs): Calculates the weight (based on the ratio of neighbors to set size) of S with neighbors N.`):

print(`OneStep(S,N,oldweight,Neighbors): One step of the metropolis algorithm using the set S with neighbors N and neighbor function Neighbors.  oldweight is the current weight of the set (so we don't have to recalculate it each time).`):

print(`Metropolis(S,K,Neighbors): Run the Metropolis algorithm K times, starting with set S, using the neighbors function Neighbors.`):


end:


#####
#AdjustNeighbors(S,N,v) 
#Adds/removes vertex v from set S (with neighbors N)
#  and adjusts neighbors set N accordingly
#  where Neighbors(v) returns all neighbors of v.
#  Assumes that v is either in B or N
AdjustNeighbors:=proc(S,N,B,v, Neighbors) local S1, N1,B1, Nv, u:
 Nv:=Neighbors(v):
 if v in S then
  S1:=S minus {v}:
  B1:= B minus {v}:
  B1:= B1 union (Neighbors(v) intersect S):
  if Nv intersect S1 = {} then
   N1:=N:
   for u in Nv do
     if Neighbors(u) intersect S1={} then N1:=N1 minus {u}: fi:
   od:
  else N1:=N union {v}:
   for u in Nv do
     if Neighbors(u) intersect S1={} then N1:=N1 minus {u} fi:
   od:
  fi:
 else
  S1:=S union {v}:
  N1:=(N union Neighbors(v)) minus S1:

  if not Nv intersect S1 = Nv then B1:= B union {v}: else B1:=B: fi:

  for u in Nv do
    if u in B and Neighbors(u) intersect N1 = {} then B1:= B1 minus {u}: fi:
  od:
  
 fi:
 [S1,N1,B1]:
end:


#Returns neighbors of v in ZxZ
ZZNeighbors:=v->{v+[1,0], v+[0,1], v-[1,0], v-[0,1]}:

#Returns the neighbors of set S
#according to the function Neighbors
SetNeighbors:=(S,Neighbors)->`union`(op(map(Neighbors,S))) minus S:


#Plots the set of points S and N in red and green repectively.
SNplot:=(S,N)->plots[display](plot(S, style = point, symbol = solidcircle, symbolsize = 10, scaling = constrained, color = red), plot(N, style = point, symbol = solidbox, symbolsize = 10, scaling = constrained, color = green)):

#Computes the Weight of the set-neighbor pair (S,N)
Weight:=proc(S,N,Neighbors):
add( nops(Neighbors(v) intersect S), v in N)/nops(S):
end:


#One step of the metropolis algorithm
OneStep:=proc(S,N,B,oldweight,Neighbors,p:=1) local a,b,v,S1,N1,B1,newweight,ratio:
 a:=nops(B):  b:=nops(N):
 v:=rand(1..a+b)():

 if v<=a then v:=B[v]:  else v:=N[v-a]: fi:

# if Neighbors=ZZNeighbors then
#  while abs(v[1])+abs(v[2]) > 20 and not v in S do
#   v:=rand(1..a+b)():
#   if v<=a then v:=B[v]:  else v:=N[v-a]: fi:   
#  od:
# fi:

 S1,N1,B1:=op(AdjustNeighbors(S,N,B,v,Neighbors)):
 
 newweight:=Weight(S1,N1,Neighbors):

 ratio:=(oldweight/newweight)^p:
 if ratio>=1 then [S1,N1,B1,newweight]:
 elif stats[random,uniform[0,1]](1) <= ratio then [S1,N1,B1,newweight]:
 else [S,N,B,oldweight]: fi:
end:


#K steps of the Metropolis algorithm, starting with S
Metropolis:=proc(S,K,Neighbors,p:=1) local N,i,L,B:
 N:=SetNeighbors(S,Neighbors):
 B := SetNeighbors(N,Neighbors) intersect S:
 L:=[S,N,B,Weight(S,N,Neighbors)]:
 for i from 1 to K do
  L:=OneStep(L[1],L[2],L[3],L[4],Neighbors,p*nops(L[1])):
 od:
L:
end: