Flat histograms everywhere

So now that we have our Wang-Landau idea we can apply it to different problems. The key idea is that we have an underlying distribution $\mathbb{G}$ which is the logarithm of the number of objects (up to an additive constant) and an increment $f$ that decreases as the quality of our $\mathbb{G}$ distribution improves. We then set our transition probabilities to be $\begin{align*} P(\varphi \to \tau) &= \min\bigg\{ 1, \exp(\mathbb{G}(\varphi) - \mathbb{G}(\tau) ) \bigg\} \end{align*}$ Using these transitions we then get a distribution that is (roughly) flat in the size of the objects.

There is, however, nothing stopping us from having $\mathbb{G}(\varphi)$ being a 1-d distribution in the size of $\varphi$ . Exactly the same ideas will work if we define $\mathbb{G}(\varphi)$ to be a 2-d distribution, say in the size of $\varphi$ and another statistic. Lets see a simple example.

Flat in more than 1 variable

A (now) classic application of SAWs is to study polymers adsorbing at a surface. Without getting into details of the physics, to study the problem, we need to count / sample SAWs in a half-space according to their length, and the number of vertices lying in the boundary of that half-space.

It is easy to tweak our code to keep the SAW in the half-space (if you try to step out of it, reject the move) and also to keep track of the number of “visits” to the boundary (just update a counter). The code for keeping the two-dimensional distribution of histograms, $\mathbb{G}$ and counts is a little more cumbersome. To make the code (a bit) easer to read I’ll use \texttt{defaultdict}, but in production code it would much faster to use 2d arrays.

Lets take a look:

from random import random, randint
from math import exp, log
from collections import defaultdict
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

N=20
S=10**5

class AdsSAW:
  xMoves=[1,-1,0,0]
  yMoves=[0,0,1,-1]

  def __init__(self):
    self.length=0
    self.walk=[(0,0)]
    self.occupied=set()
    self.visits=[0]

  def lv(self): #tuple of length,visits
    return( self.length, self.visits[-1])

  def doMove(self):
      xn,yn=self.walk[-1] #last vertex
      c=randint(0,3); xn += AdsSAW.xMoves[c]; yn += AdsSAW.yMoves[c]

      if(self.length==0): #cannot delete
          self.addStep( (xn,yn) )
      else:
          if( (xn,yn)==self.walk[-2] ):
              self.deleteStep() #immediate reversal
          else:
              self.addStep( (xn,yn) )

  def deleteStep(self):
    if( log(random()) > G[self.length,self.visits[-1]]-G[self.length-1,self.visits[-2]] ):
        return()
    x,y=self.walk.pop()
    self.occupied.remove( (x,y) )
    self.visits.pop()
    self.length-=1

  def addStep(self, p):
    if(p[1]<0): #stepping outside the half-space
      return
    v=self.visits[-1] 
    if(p[1]==0):#update number of visits
      v+=1
    if(self.length==N or log(random()) > G[self.length,self.visits[-1]]-G[self.length+1,v] ):
        return
    if( p not in self.occupied ):
        self.visits.append(v)        
        self.walk.append(p)
        self.occupied.add(p)
        self.length+=1

counts=defaultdict(int) #total samples
H = defaultdict(int) #running histogram
G = defaultdict(float) #G-distribution
f = 1 #G-incrementer

def flatten():
  global f,H,G
  r=min(H.values())/max(H.values())
  if(r>0.9): #time to reduce f
    f=f/2
    H = defaultdict(float) #reset histogram
    g0 = min(G.values()) 
    G = defaultdict(float, {k: X-g0 for k,X in G.items()} ) #avoid overflows
    print("Histogram is flat = ",r, "\tReduce f to ",f)
  else:
    print("Histogram not flat enough", r, "\tf left as ", f)

phi=AdsSAW()

while(f >=0.01):
  for s in range(S):
    phi.doMove()
    if( s%17==0 ):
      lv=phi.lv()
      counts[ lv ]+=1
      G[ lv ]+=f
      H[ lv ]+=1

  flatten()

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
x,y = zip(*counts.keys() )
z=counts.values()
w = list(map( lambda X:0.5, x))
b = list(map( lambda X:0, x))
ax.bar3d(x,y,b,w,w,z)
ax.set_xlabel('length'); ax.set_ylabel('visits')
plt.show()

for k in range(0,5):
  for j in range(0,k+1):
    print( " + {:f}*a^{:d}".format( exp(G[k,j]-G[0,0]), j ), end='' )
  print()

One thing that you notice immediately when running this code is that it takes much longer to converge — this is because our random walk through the state-space needs to difuse over a 2d space of (length,visits) rather than just over a 1d space of (length). Anyway - you can see that we’ve reduced the max length and we (finally) get a nice flat histogram output.

2d flat histogram

Of course we should really do something more meaningful with this data — like compute specific heats or…

Polygons and BFACF

Rather than going further into analysing SAWs and their interactions, lets take a look at polygons and the BFACF algorithm which is used to study them while keeping their topology fixed. If we think about trying to apply our Berretti-Sokal to simulate a polygon we immedately realise that there are no ends to add to or delete from.

BFACF — named for Berg, Foster (1981) and Arago de Varvalho, Caracciolo and Frohlich (1983) — acts on an edge and neighbouring face of a polygon:

BFACF moves

So — the basic move consists of

pick and edge and a neighbouring face
look at the edges around that face — if occupied, vacate them, and if vacant, occupy them.
if the result is a legal SAP, then keep it, else reject it and keep current configuration.

Each move changes the total length of the polygon by either $+2, 0$ or $-2$ . Traditionally we would attach a Metropolis-Hastings step to each move with a parameter $\beta$ , so that if a move changes the length by $0,-2$ then we accept it, while if it increases the length by $+2$ we accept it with probability $\beta^2$ and otherwise reject it. Detailed balance then gives us $\begin{align*} P(\varphi_n \to \tau_{n-2}) \cdot \pi(\varphi_n) &= P(\tau_{n-2} \to \varphi_n) \cdot \pi(\tau_{n-2}) \\ \pi(\varphi_n) \cdot \frac{1}{4n} &= \frac{\beta^2}{4(n-2)} \cdot \pi(\tau_{n-2}) \end{align*}$ where the factor of $4n$ comes from chooing one of the $n$ edges and one of the 4 faces adjacent to it. Hence we have $\begin{align*} \pi(\varphi_n) &\propto n \cdot \beta^{2n} \end{align*}$ The main reason that we use this algorithm so much in topological applications rather than the much faster pivot algorithm, is that it has been shown by Buks van Rensburg (1991) that BFACF conserves topology (Buks and I (2010) generalised this to a couple of other lattices). More precisely the ergodicity classes of the algorithm are exactly the knot types. ie - if you start with a figure-8 knot, then BFACF will allow you to sample uniformly from the space of all figure-8 knots. Further — this also applies to links — start with a hopf-link and you’ll sample uniformly from the space of all hopf-links.

The coding of this algorithm has a couple of annoying technicalities that we need to sort out. Lets start at the beginning:

We need to be able to navigate around our polygon — if I give you a vertex, you need to be able to find the next vertex and previous vertex around the polygon.
We need to be able to pick an edge uniformly at random and investigate the adjacent face(s).

We can sort these problems by storing

the set of vertices
a map that tells us the edges coming into and out from each vertex
we encode the direction of edges just as $0,1,2$ as $+x,+y,+z$ and $3,4,5$ as $-x,-y,-z$ (respectively).

This means we can very easily determine

is a given vertex occupied (just set the set)
select a vertex at random (python will happily sample from a set — a bit harder in c++)
given a direction, what is the reverse direction (just $\pm 3$ )
given a vertex, what is the direction of the outgoing edge, and so what are the previous and next vertices (look up the incoming and outgoing edges, and step along them)
given a vertex what are the 4 adjacent faces (given by the 4 vectors normal to the outgoing edge) and pick one of them at random.

This last one is a bit more tedious — but if our edge is in direction $d$ , then a vector in direction $v$ is normal to it provided $d \mod 3 \neq v \mod 3$ .

Lets look at what happens when we try a move.

Pick a vertex uniformly at random — call it $x$ .
Find the edge outgoing from $x$ and call it $eXtoY$ — we need its direction.
Find the next vertex along the polygon — call it $y$ .
We’ll also need the vertex previous to $x$ and the vertex after $y$ — call them $px$ and $ny$ .
Pick a direction uniformly at random normal to $eXtoY$ — this determines the face adjacent to the edge from $x$ to $y$ and so the face on which we’ll do out move. Call this direction $direction$ , and its reverse $rdirection$ .
Let $X$ be the vertex one unit in direction $direction$ from $x$ , and
let $Y$ be the vertex one unit in direction $direction$ from $y$ .

Here is a picture of the story so far:

BFACF pre-move

At present we’ve made no assumptions about $X,Y$ and $px,ny$ . But lets walk through what can happen:

If $X=px$ and $Y=py$ then our move will reduce length by $2$
If $X \neq px$ and $Y\neq py$ then our move will increase length by $2$
If $X=px$ but $Y\neq py$ then the move will not change length, and
If $Y=ny$ but $X\neq px$ then the move will not change length.

These 4 possibilities are shown (in order) below

Four possible BFACF states

So in order to do the above moves we must (in order)

we don’t need to check self-avoidance, connect $px \to ny$ and delete vertices $x,y$ .
we need to check that both $X,Y$ are occupied — if they are then reject since self-avoidance violated. Else add $X$ and $Y$ and connect $x \to X \to Y \to y$ .
check if $Y$ is occupied — if so then reject since self-avoidance violated. Else, delete $x$ and connect $px \to Y \to y$ .
check if $X$ is occupied — if so then reject since self-avoidance violated. Else, delete $y$ and connect $x \to X \to ny$ .

Oof! The resulting code is a bit messy. The same ideas work on the square-lattice, triangular lattice, FCC-lattice and (with some more work) BCC-lattice.

Anyway — we still need to make sure we put in our Metropolis-Hastings accept-reject moves. And (finally) we get the following code:

from collections import defaultdict
from random import randint, random, sample
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

BETA=0.05
N=100
S=10**6

#ENUWSD=012345 
#so reverse = dir+/-3
xMove = [1,0,0,-1,0,0]
yMove = [0,1,0,0,-1,0]
zMove = [0,0,1,0,0,-1]

# edge i -> {a,i,b}, so direction of a,b != direction of i or reverse of i.
def UConf(direction):
  # d -> not d or (d+/-3), ie c%3!=d%3
  c = randint(0,5)
  while(c%3==direction%3):
    c = randint(0,5)
  return(c)

def step(p, d):
  return( ( p[0]+xMove[d],p[1]+yMove[d],p[2]+zMove[d]) )

def rstep(p, d):
  return( ( p[0]-xMove[d],p[1]-yMove[d],p[2]-zMove[d]) )

def reverseDirection(d):
  if(d>=3):
    return(d-3)
  else:
    return(d+3)

def diff(p,p0):
  return( (p[0]-p0[0], p[1]-p0[1], p[2]-p0[2]) )

class SAP:
  def __init__(self):
    self.length=4
    self.vertices=set() 
    self.vertices.add( (0,0,0) )
    self.vertices.add( (1,0,0) )
    self.vertices.add( (1,1,0) )
    self.vertices.add( (0,1,0) )

    self.edges=defaultdict(list)
    self.edges[ (0,0,0) ]=(4 ,0)
    self.edges[ (1,0,0) ]=(0 ,1)
    self.edges[ (1,1,0) ]=(1 ,3)
    self.edges[ (0,1,0) ]=(3 ,4)

  def preparePlot(self):
    p0 = sample(self.vertices,1)[0] # a vertex at random
    p=p0
    xs=[0];ys=[0];zs=[0]
    for k in range(0,self.length):
      p=step(p, self.edges[p][1] )
      xs.append(p[0]-p0[0])
      ys.append(p[1]-p0[1])
      zs.append(p[2]-p0[2])
    return(xs,ys,zs)

  def printV(self):
    p0 = sample(self.vertices,1)[0] # a vertex at random
    p=p0
    print( "[", diff(p,p0), end='')
    for k in range(0,self.length):
      p=step(p, self.edges[p][1] )
      print( ",", diff(p,p0), end='' )
    print("]")

  def printVE(self):
    p0 = sample(self.vertices,1)[0] # a vertex at random
    p=p0
    print("**********************************")
    print( diff(p,p0), self.edges[p] )
    for k in range(0,self.length):
      p=step(p, self.edges[p][1] )
      print( diff(p,p0), self.edges[p] )

  def BFACF(self):
    x = sample(self.vertices,1)[0] # a vertex at random
    eXtoY = self.edges[x] # edges in and out from it
    y = step(x,eXtoY[1]) # next vertex
    px = rstep(x,eXtoY[0]) # previous vertex
    ny = step( y, self.edges[y][1] ) # next^2 vertex

    #pick direction to push edge x->y
    direction = UConf(eXtoY[1])
    rdirection = reverseDirection(direction)    
    XX = step(x,direction)
    YY = step(y,direction)

    # check if new XX,YY are px or ny
    # use that to compute delta-length
    if(XX==px):
      if(YY==ny):
        delta = -2;
      else:
        delta=0;
    else: # X!=px
      if(YY==ny):
        delta=0; 
      else:
        delta=2;

    if(delta==0):
      if(self.length==4): 
        return
      if(XX==px): #connect px->YY->y, delete x if legal
        if(YY in self.vertices):
          return #not self-avoiding
        else:
          self.vertices.remove(x); del self.edges[x]
          self.vertices.add(YY)
          self.edges[px]=(self.edges[px][0],eXtoY[1])
          self.edges[YY]=(eXtoY[1], rdirection)
          self.edges[y]=(rdirection,self.edges[y][1])
      else: #must have YY=ny, connect x->XX->ny, delete y if legal
        if(XX in self.vertices):
          return #not self-avoiding
        else:
          self.vertices.remove(y); del self.edges[y]
          self.vertices.add(XX)
          self.edges[x]=(self.edges[x][0],direction)
          self.edges[XX]=(direction, eXtoY[1])
          self.edges[ny]=(eXtoY[1],self.edges[ny][1])
    elif(delta==2):
      if(XX in self.vertices or YY in self.vertices or self.length>=N):
        return #not self-avoiding
      if(random()<BETA): #connect x->XX->YY->y
        self.vertices.add(XX); self.vertices.add(YY)
        self.edges[x]=(self.edges[x][0],direction)
        self.edges[XX]=(direction,eXtoY[1])
        self.edges[YY]=(eXtoY[1],rdirection)
        self.edges[y]=(rdirection, self.edges[y][1])
        self.length+=2
    else: #delta=-2
      if(self.length==4): 
        return
      #connect px->ny, delete x,y
      self.vertices.remove(x); del self.edges[x]
      self.vertices.remove(y); del self.edges[y]
      self.edges[px]=(self.edges[px][0],eXtoY[1])
      self.edges[ny]=(eXtoY[1],self.edges[ny][1])
      self.length-=2

phi=SAP()
counts=[0 for k in range(0,N+1)]

for s in range(0,S):
  phi.BFACF()
  if(s%17==0):
    counts[phi.length]+=1

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
xs,ys,zs = phi.preparePlot()
ax.plot(xs,ys,zs)
ax.scatter(xs,ys,zs)
ax.set_xlabel('x'); ax.set_ylabel('y');  ax.set_zlabel('z')
plt.show()

plt.bar(range(0,N+1),counts)
plt.show()

We can take a look at one of the 3d polygons it produces:

Sample 3d SAP

and the sample histogram:

Sample 3d SAP

This isn’t so great. Lets flatten it.

Flattened BFACF

We don’t need to make too many changes — just import the wang-landau stuff we did earlier and make sure we put the Metropolis-Hastings steps in there when the length changes. One annoying thing to take care of is that we only have polygons of even lengths, and the smallest one is the unit square which has length $4$ . So to avoid the odd-even problem, we’ll use half-length and to take care of $\mathbb{H}(0)=\mathbb{H}(1)=0$ we need to compute the max and min of our histogram more carefully.

Doing all this gives

from collections import defaultdict
from random import randint, random, sample
from math import log,exp
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

N=50
S=10**5

#ENUWSD=012345 
#so reverse = dir+/-3
xMove = [1,0,0,-1,0,0]
yMove = [0,1,0,0,-1,0]
zMove = [0,0,1,0,0,-1]

# edge i -> {a,i,b}, so direction of a,b != direction of i or reverse of i.
def UConf(direction):
  # d -> not d or (d+/-3), ie c%3!=d%3
  c = randint(0,5)
  while(c%3==direction%3):
    c = randint(0,5)
  return(c)

def step(p, d):
  return( ( p[0]+xMove[d],p[1]+yMove[d],p[2]+zMove[d]) )

def rstep(p, d):
  return( ( p[0]-xMove[d],p[1]-yMove[d],p[2]-zMove[d]) )

def reverseDirection(d):
  if(d>=3):
    return(d-3)
  else:
    return(d+3)

def diff(p,p0):
  return( (p[0]-p0[0], p[1]-p0[1], p[2]-p0[2]) )

class SAP:
  def __init__(self):
    self.halfLength=2
    self.vertices=set() 
    self.vertices.add( (0,0,0) )
    self.vertices.add( (1,0,0) )
    self.vertices.add( (1,1,0) )
    self.vertices.add( (0,1,0) )

    self.edges=defaultdict(list)
    self.edges[ (0,0,0) ]=(4 ,0)
    self.edges[ (1,0,0) ]=(0 ,1)
    self.edges[ (1,1,0) ]=(1 ,3)
    self.edges[ (0,1,0) ]=(3 ,4)

  def preparePlot(self):
    p0 = sample(self.vertices,1)[0] # a vertex at random
    p=p0
    xs=[0];ys=[0];zs=[0]
    for k in range(0,self.halfLength*2):
      p=step(p, self.edges[p][1] )
      xs.append(p[0]-p0[0])
      ys.append(p[1]-p0[1])
      zs.append(p[2]-p0[2])
    return(xs,ys,zs)

  def printV(self):
    p0 = sample(self.vertices,1)[0] # a vertex at random
    p=p0
    print( "[", diff(p,p0), end='')
    for k in range(0,self.halfLength*2):
      p=step(p, self.edges[p][1] )
      print( ",", diff(p,p0), end='' )
    print("]")

  def printVE(self):
    p0 = sample(self.vertices,1)[0] # a vertex at random
    p=p0
    print("**********************************")
    print( diff(p,p0), self.edges[p] )
    for k in range(0,self.halfLength*2):
      p=step(p, self.edges[p][1] )
      print( diff(p,p0), self.edges[p] )

  def BFACF(self):
    x = sample(self.vertices,1)[0] # a vertex at random
    eXtoY = self.edges[x] # edges in and out from it
    y = step(x,eXtoY[1]) # next vertex
    px = rstep(x,eXtoY[0]) # previous vertex
    ny = step( y, self.edges[y][1] ) # next^2 vertex

    #pick direction to push edge x->y
    direction = UConf(eXtoY[1])
    rdirection = reverseDirection(direction)    
    XX = step(x,direction)
    YY = step(y,direction)

    # check if new XX,YY are px or ny
    # use that to compute delta-length
    if(XX==px):
      if(YY==ny):
        delta = -2;
      else:
        delta=0;
    else: # X!=px
      if(YY==ny):
        delta=0; 
      else:
        delta=2;

    if(delta==0):
      if(self.halfLength==2): 
        return
      if(XX==px): #connect px->YY->y, delete x if legal
        if(YY in self.vertices):
          return #not self-avoiding
        else:
          self.vertices.remove(x); del self.edges[x]
          self.vertices.add(YY)
          self.edges[px]=(self.edges[px][0],eXtoY[1])
          self.edges[YY]=(eXtoY[1], rdirection)
          self.edges[y]=(rdirection,self.edges[y][1])
      else: #must have YY=ny, connect x->XX->ny, delete y if legal
        if(XX in self.vertices):
          return #not self-avoiding
        else:
          self.vertices.remove(y); del self.edges[y]
          self.vertices.add(XX)
          self.edges[x]=(self.edges[x][0],direction)
          self.edges[XX]=(direction, eXtoY[1])
          self.edges[ny]=(eXtoY[1],self.edges[ny][1])
    elif(delta==2):
      if(XX in self.vertices or YY in self.vertices or self.halfLength>=N):
        return #not self-avoiding
      if( log(random()) > G[self.halfLength]-G[self.halfLength+1] ):
        return
      else: #connect x->XX->YY->y
        self.vertices.add(XX); self.vertices.add(YY)
        self.edges[x]=(self.edges[x][0],direction)
        self.edges[XX]=(direction,eXtoY[1])
        self.edges[YY]=(eXtoY[1],rdirection)
        self.edges[y]=(rdirection, self.edges[y][1])
        self.halfLength+=1
    else: #delta=-2
      if(self.halfLength==2): 
        return
      if( log(random()) > G[self.halfLength]-G[self.halfLength-1] ):
        return
      #connect px->ny, delete x,y
      self.vertices.remove(x); del self.edges[x]
      self.vertices.remove(y); del self.edges[y]
      self.edges[px]=(self.edges[px][0],eXtoY[1])
      self.edges[ny]=(eXtoY[1],self.edges[ny][1])
      self.halfLength-=1

counts=[0 for i in range(N+1)] #total samples
H = [0 for k in range(N+1)] #running histogram
G = [0 for k in range(N+1)] #G-distribution
f = 1 #G-incrementer

def flatten():
  global f,H,G
  ## because H[0]=H[1]=0 more care required
  mnH=10**10
  mxH=-1
  for k in range(2,N):
    if(H[k]<mnH):
      mnH=H[k]
    if(H[k]>mxH):
      mxH=H[k]

  r=mnH/mxH
  if(r>0.9): #time to reduce f
    f=f/2
    H = [0 for k in range(N+1)] #reset histogram
    g0 = G[2] #shortest half-length
    G = list(map(lambda X: X-g0, G) ) #avoid overflows
    print("Histogram is flat = ",r, "\tReduce f to ",f)
  else:
    print("Histogram not flat enough", r, "\tf left as ", f)


phi=SAP()
counts=[0 for k in range(0,N+1)]

while(f >=0.01):
  for s in range(S):
    phi.BFACF()
    if( s%17==0 ):
      counts[ phi.halfLength ]+=1
      G[ phi.halfLength ]+=f
      H[ phi.halfLength ]+=1

  flatten()

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
xs,ys,zs = phi.preparePlot()
ax.plot(xs,ys,zs)
ax.scatter(xs,ys,zs)
ax.set_xlabel('x'); ax.set_ylabel('y');  ax.set_zlabel('z')
plt.show()

plt.bar(range(0,N+1),counts)
plt.show()

and a sample unknotted SAP:

A sample 3d unknot

and the sample histogram:

Histogram for 3d SAPs

This is much better.

Now — if we want to sample a different knot type, we just need to change the initialisation stuff to encode a trefoil or a figure-8 or a whatever. This can be a bit tedious since one needs to encode the vertices and edges of a small knot.

Lets move on to a Hopf-link.

Hopf hopf hopf

The biggest change that happens here is that we need to make sure we keep the two components of the hopf-link separated. Because of this we’ll keep two buckets of vertices — one for each component. We will, however, need to be a bit careful when we (now) pick a vertex uniformly at random to do our BFACF move. If we have $n$ vertices in the first component and $k$ vertices in the second, then to pick one uniformly at random, we first choose a random integer in $[1,n+k]$ . If it is in $[1,n]$ then we pick randomly from the first component, and otherwise we pick randomly from the second. Otherwise our code is very similar.

One other pain is coding in the link itself and its edges.

from collections import defaultdict
from random import randint, random, sample
from math import log,exp
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

N=200
S=10**5

#ENUWSD=012345 
#so reverse = dir+/-3
xMove = [1,0,0,-1,0,0]
yMove = [0,1,0,0,-1,0]
zMove = [0,0,1,0,0,-1]

# edge i -> {a,i,b}, so direction of a,b != direction of i or reverse of i.
def UConf(direction):
  # d -> not d or (d+/-3), ie c%3!=d%3
  c = randint(0,5)
  while(c%3==direction%3):
    c = randint(0,5)
  return(c)

def step(p, d):
  return( ( p[0]+xMove[d],p[1]+yMove[d],p[2]+zMove[d]) )

def rstep(p, d):
  return( ( p[0]-xMove[d],p[1]-yMove[d],p[2]-zMove[d]) )

def reverseDirection(d):
  if(d>=3):
    return(d-3)
  else:
    return(d+3)

def diff(p,p0):
  return( (p[0]-p0[0], p[1]-p0[1], p[2]-p0[2]) )



class Hopf:
  def __init__(self):
    self.halfLength=[4,4]
    self.vertices=[set(),set()]
    self.edges=defaultdict(list)

    self.insertVertexEdge(0, (0,0,0), (4, 0));
    self.insertVertexEdge(0, (1,0,0), (0, 0));
    self.insertVertexEdge(0, (2,0,0), (0, 1));
    self.insertVertexEdge(0, (2,1,0), (1, 1));
    self.insertVertexEdge(0, (2,2,0), (1, 3));
    self.insertVertexEdge(0, (1,2,0), (3, 3));
    self.insertVertexEdge(0, (0,2,0), (3, 4));
    self.insertVertexEdge(0, (0,1,0), (4, 4));

    self.insertVertexEdge(1, (1,1,0), (2, 2));
    self.insertVertexEdge(1, (1,1,1), (2, 0));
    self.insertVertexEdge(1, (2,1,1), (0, 0));
    self.insertVertexEdge(1, (3,1,1), (0, 5));
    self.insertVertexEdge(1, (3,1,0), (5, 5));
    self.insertVertexEdge(1, (3,1,-1), (5, 3));
    self.insertVertexEdge(1, (2,1,-1), (3, 3));
    self.insertVertexEdge(1, (1,1,-1), (3, 2));

  def insertVertexEdge(self, component, v, e):
      self.vertices[component].add(v)
      self.edges[v]=e

  def totalHL(self):
      return( self.halfLength[0]+self.halfLength[1] )

  def preparePlot(self):
    p0 = sample(self.vertices[0],1)[0] # a vertex at random
    p=p0
    xs=[0];ys=[0];zs=[0]
    for k in range(0,self.halfLength[0]*2):
      p=step(p, self.edges[p][1] )
      xs.append(p[0]-p0[0])
      ys.append(p[1]-p0[1])
      zs.append(p[2]-p0[2])

    p=sample(self.vertices[1],1)[0]
    xt=[p[0]-p0[0]];yt=[p[1]-p0[1]];zt=[p[2]-p0[2]]
    for k in range(0,self.halfLength[1]*2):
      p=step(p, self.edges[p][1] )
      xt.append(p[0]-p0[0])
      yt.append(p[1]-p0[1])
      zt.append(p[2]-p0[2])
    return([xs,ys,zs,xt,yt,zt])

  def printV(self, component):
    p0 = sample(self.vertices[component],1)[0] # a vertex at random
    p=p0
    print( "[", diff(p,p0), end='')
    for k in range(0,self.halfLength[component]*2):
      p=step(p, self.edges[p][1] )
      print( ",", diff(p,p0), end='' )
    print("]")

  def printVE(self, component):
    p0 = sample(self.vertices[component],1)[0] # a vertex at random
    p=p0
    print("**********************************")
    print( diff(p,p0), self.edges[p] )
    for k in range(0,self.halfLength[component]*2):
      p=step(p, self.edges[p][1] )
      print( diff(p,p0), self.edges[p] )

  def BFACF(self):
    #pick a vertex from a component 
    if(randint(1,self.halfLength[0]+self.halfLength[1]) <= self.halfLength[0]):
        comp=0
    else:
        comp=1

    x = sample(self.vertices[comp],1)[0] # a vertex at random
    eXtoY = self.edges[x] # edges in and out from it
    y = step(x,eXtoY[1]) # next vertex
    px = rstep(x,eXtoY[0]) # previous vertex
    ny = step( y, self.edges[y][1] ) # next^2 vertex

    #pick direction to push edge x->y
    direction = UConf(eXtoY[1])
    rdirection = reverseDirection(direction)    
    XX = step(x,direction)
    YY = step(y,direction)

    # check if new XX,YY are px or ny
    # use that to compute delta-length
    if(XX==px):
      if(YY==ny):
        delta = -2;
      else:
        delta=0;
    else: # X!=px
      if(YY==ny):
        delta=0; 
      else:
        delta=2;

    if(delta==0):
      if(self.halfLength[comp]==2): 
        return
      if(XX==px): #connect px->YY->y, delete x if legal
        if(YY in self.vertices[0] or YY in self.vertices[1]):
          return #not self-avoiding
        else:
          self.vertices[comp].remove(x); del self.edges[x]
          self.vertices[comp].add(YY)
          self.edges[px]=(self.edges[px][0],eXtoY[1])
          self.edges[YY]=(eXtoY[1], rdirection)
          self.edges[y]=(rdirection,self.edges[y][1])
      else: #must have YY=ny, connect x->XX->ny, delete y if legal
        if(XX in self.vertices[0] or XX in self.vertices[1]):
          return #not self-avoiding
        else:
          self.vertices[comp].remove(y); del self.edges[y]
          self.vertices[comp].add(XX)
          self.edges[x]=(self.edges[x][0],direction)
          self.edges[XX]=(direction, eXtoY[1])
          self.edges[ny]=(eXtoY[1],self.edges[ny][1])
    elif(delta==2):
      if(XX in self.vertices[0] or YY in self.vertices[0] or XX in self.vertices[1] or YY in self.vertices[1] or self.totalHL()>=N):
        return #not self-avoiding
      if( log(random()) > G[self.totalHL()]-G[self.totalHL()+1] ):
        return
      else: #connect x->XX->YY->y
        self.vertices[comp].add(XX); self.vertices[comp].add(YY)
        self.edges[x]=(self.edges[x][0],direction)
        self.edges[XX]=(direction,eXtoY[1])
        self.edges[YY]=(eXtoY[1],rdirection)
        self.edges[y]=(rdirection, self.edges[y][1])
        self.halfLength[comp]+=1
    else: #delta=-2
      if(self.halfLength==2): 
        return
      if( log(random()) > G[self.totalHL()]-G[self.totalHL()-1] ):
        return
      #connect px->ny, delete x,y
      self.vertices[comp].remove(x); del self.edges[x]
      self.vertices[comp].remove(y); del self.edges[y]
      self.edges[px]=(self.edges[px][0],eXtoY[1])
      self.edges[ny]=(eXtoY[1],self.edges[ny][1])
      self.halfLength[comp]-=1

counts=[0 for i in range(N+1)] #total samples
H = [0 for k in range(N+1)] #running histogram
G = [0 for k in range(N+1)] #G-distribution
f = 2 #G-incrementer

def flatten():
  global f,H,G
  ## because H[0]=H[1]=0 more care required
  mnH=10**10
  mxH=-1
  for k in range(8,N):
    if(H[k]<mnH):
      mnH=H[k]
    if(H[k]>mxH):
      mxH=H[k]
  print("Min hist = ", mnH, "\tMax hist = ",mxH)
  r=mnH/mxH
  if(r>0.5): #time to reduce f
    f=f/2
    H = list(map(lambda X: 0, H) ) #reset histogram
    g0 = G[8] #shortest halfLength.
    G = list(map(lambda X: X-g0, G) ) #avoid overflows
    print("Histogram is flat = ",r, "\tReduce f to ",f)
  else:
    print("Histogram not flat enough", r, "\tf left as ", f)


phi=Hopf()

while(f >0.1):
  for s in range(S):
    phi.BFACF()
    if( s%17==0 ):
      counts[ phi.totalHL() ]+=1
      G[ phi.totalHL() ]+=f
      H[ phi.totalHL() ]+=1

  flatten()

#make sure we get a reasonably long sample
while(phi.totalHL()<N*0.5): 
    phi.BFACF()

print( phi.halfLength )
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
comps = phi.preparePlot()
ax.plot(comps[0],comps[1],comps[2],'bo-')
ax.plot(comps[3],comps[4],comps[5],'ro-')
ax.set_xlabel('x'); ax.set_ylabel('y');  ax.set_zlabel('z')
plt.show()


plt.bar(range(0,N+1),counts)
plt.show()

and a sample hopf-link:

Sample Hopf link

So we can see here that a “typical” hopf link wants to maximise the size of one link while the other shrinks. Indeed, if we do our simulations carefully and let stuff run for a while then we’ll find that a typical hopf link of total length $n$ , the component sizes grow as $\begin{align*} \langle \max\{l_1,l_2\} \rangle_n &= n - o(n)\\ \langle \min\{l_1,l_2\} \rangle_n &= o(n). \end{align*}$

The sample histogram isn’t too bad — especially given this was run only for a few minutes (on my desktop computer):

Sample Hopf histogram

Flatter Hopf

One way to overcome the issue of “big-link-small-link” is to flatten our simulation in the length of each component. Of course, it will take quite a bit longer to flatten as our random-walk over $(l_1,l_2)$ space is now diffusing in a 2d space instead of a 1d-line. But it is not too hard to combine our adsorbing-SAW code with our BFACF-hopf code to get something that works.

from collections import defaultdict
from random import randint, random, sample
from math import log,exp,floor
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

N=50
S=10**5

#ENUWSD=012345 
#so reverse = dir+/-3
xMove = [1,0,0,-1,0,0]
yMove = [0,1,0,0,-1,0]
zMove = [0,0,1,0,0,-1]

# edge i -> {a,i,b}, so direction of a,b != direction of i or reverse of i.
def UConf(direction):
  # d -> not d or (d+/-3), ie c%3!=d%3
  c = randint(0,5)
  while(c%3==direction%3):
    c = randint(0,5)
  return(c)

def step(p, d):
  return( ( p[0]+xMove[d],p[1]+yMove[d],p[2]+zMove[d]) )

def rstep(p, d):
  return( ( p[0]-xMove[d],p[1]-yMove[d],p[2]-zMove[d]) )

def reverseDirection(d):
  if(d>=3):
    return(d-3)
  else:
    return(d+3)

def diff(p,p0):
  return( (p[0]-p0[0], p[1]-p0[1], p[2]-p0[2]) )



class Hopf:
  def __init__(self):
    self.halfLength=[4,4]
    self.vertices=[set(),set()]
    self.edges=defaultdict(list)

    self.insertVertexEdge(0, (0,0,0), (4, 0));
    self.insertVertexEdge(0, (1,0,0), (0, 0));
    self.insertVertexEdge(0, (2,0,0), (0, 1));
    self.insertVertexEdge(0, (2,1,0), (1, 1));
    self.insertVertexEdge(0, (2,2,0), (1, 3));
    self.insertVertexEdge(0, (1,2,0), (3, 3));
    self.insertVertexEdge(0, (0,2,0), (3, 4));
    self.insertVertexEdge(0, (0,1,0), (4, 4));

    self.insertVertexEdge(1, (1,1,0), (2, 2));
    self.insertVertexEdge(1, (1,1,1), (2, 0));
    self.insertVertexEdge(1, (2,1,1), (0, 0));
    self.insertVertexEdge(1, (3,1,1), (0, 5));
    self.insertVertexEdge(1, (3,1,0), (5, 5));
    self.insertVertexEdge(1, (3,1,-1), (5, 3));
    self.insertVertexEdge(1, (2,1,-1), (3, 3));
    self.insertVertexEdge(1, (1,1,-1), (3, 2));

  def insertVertexEdge(self, component, v, e):
      self.vertices[component].add(v)
      self.edges[v]=e

  def totalHL(self):
      return( self.halfLength[0]+self.halfLength[1] )

  def deltaHL(self):
      return( abs(self.halfLength[0]-self.halfLength[1] ))

  def currentHL(self):
      return( (self.halfLength[0], self.halfLength[1]) )

  def futureHL(self, c, v):
      if(c==0):
          return( self.halfLength[0]+v, self.halfLength[1] )
      else:
          return( self.halfLength[0], self.halfLength[1]+v )

  def preparePlot(self):
    p0 = sample(self.vertices[0],1)[0] # a vertex at random
    p=p0
    xs=[0];ys=[0];zs=[0]
    for k in range(0,self.halfLength[0]*2):
      p=step(p, self.edges[p][1] )
      xs.append(p[0]-p0[0])
      ys.append(p[1]-p0[1])
      zs.append(p[2]-p0[2])

    p=sample(self.vertices[1],1)[0]
    xt=[p[0]-p0[0]];yt=[p[1]-p0[1]];zt=[p[2]-p0[2]]
    for k in range(0,self.halfLength[1]*2):
      p=step(p, self.edges[p][1] )
      xt.append(p[0]-p0[0])
      yt.append(p[1]-p0[1])
      zt.append(p[2]-p0[2])
    return([xs,ys,zs,xt,yt,zt])

  def printV(self, component):
    p0 = sample(self.vertices[component],1)[0] # a vertex at random
    p=p0
    print( "[", diff(p,p0), end='')
    for k in range(0,self.halfLength[component]*2):
      p=step(p, self.edges[p][1] )
      print( ",", diff(p,p0), end='' )
    print("]")

  def printVE(self, component):
    p0 = sample(self.vertices[component],1)[0] # a vertex at random
    p=p0
    print("**********************************")
    print( diff(p,p0), self.edges[p] )
    for k in range(0,self.halfLength[component]*2):
      p=step(p, self.edges[p][1] )
      print( diff(p,p0), self.edges[p] )

  def BFACF(self):
    #pick a vertex from a component 
    if(randint(1,self.halfLength[0]+self.halfLength[1]) <=self.halfLength[0]):
        comp=0
    else:
        comp=1

    x = sample(self.vertices[comp],1)[0] # a vertex at random
    eXtoY = self.edges[x] # edges in and out from it
    y = step(x,eXtoY[1]) # next vertex
    px = rstep(x,eXtoY[0]) # previous vertex
    ny = step( y, self.edges[y][1] ) # next^2 vertex

    #pick direction to push edge x->y
    direction = UConf(eXtoY[1])
    rdirection = reverseDirection(direction)    
    XX = step(x,direction)
    YY = step(y,direction)

    # check if new XX,YY are px or ny
    # use that to compute delta-length
    if(XX==px):
      if(YY==ny):
        delta = -2;
      else:
        delta=0;
    else: # X!=px
      if(YY==ny):
        delta=0; 
      else:
        delta=2;

    if(delta==0):
      if(self.halfLength[comp]==2): 
        return
      if(XX==px): #connect px->YY->y, delete x if legal
        if(YY in self.vertices[0] or YY in self.vertices[1]):
          return #not self-avoiding
        else:
          self.vertices[comp].remove(x); del self.edges[x]
          self.vertices[comp].add(YY)
          self.edges[px]=(self.edges[px][0],eXtoY[1])
          self.edges[YY]=(eXtoY[1], rdirection)
          self.edges[y]=(rdirection,self.edges[y][1])
      else: #must have YY=ny, connect x->XX->ny, delete y if legal
        if(XX in self.vertices[0] or XX in self.vertices[1]):
          return #not self-avoiding
        else:
          self.vertices[comp].remove(y); del self.edges[y]
          self.vertices[comp].add(XX)
          self.edges[x]=(self.edges[x][0],direction)
          self.edges[XX]=(direction, eXtoY[1])
          self.edges[ny]=(eXtoY[1],self.edges[ny][1])
    elif(delta==2):
      if(XX in self.vertices[0] or YY in self.vertices[0] or XX in self.vertices[1] or YY in self.vertices[1] or self.totalHL()>=N):
        return #not self-avoiding
      if( log(random()) > G[self.currentHL()]-G[self.futureHL(comp,1)] ):
        return
      else: #connect x->XX->YY->y
        self.vertices[comp].add(XX); self.vertices[comp].add(YY)
        self.edges[x]=(self.edges[x][0],direction)
        self.edges[XX]=(direction,eXtoY[1])
        self.edges[YY]=(eXtoY[1],rdirection)
        self.edges[y]=(rdirection, self.edges[y][1])
        self.halfLength[comp]+=1
    else: #delta=-2
      if(self.halfLength==2): 
        return
      if( log(random()) > G[self.currentHL()]-G[self.futureHL(comp,-1)] ):
        return
      #connect px->ny, delete x,y
      self.vertices[comp].remove(x); del self.edges[x]
      self.vertices[comp].remove(y); del self.edges[y]
      self.edges[px]=(self.edges[px][0],eXtoY[1])
      self.edges[ny]=(eXtoY[1],self.edges[ny][1])
      self.halfLength[comp]-=1

counts=defaultdict(int) #total samples
H = defaultdict(int) #running histogram
G = defaultdict(float) #G-distribution
f = 2 #G-incrementer

def flatten():
  global f,H,G
  r=min(H.values())/max(H.values())
  if(r>0.5): #time to reduce f
    f=f/2
    H = defaultdict(float) #reset histogram
    g0 = G[0,0] 
    G = defaultdict(float, {k: X-g0 for k,X in G.items()} ) #avoid overflows
    print("Histogram is flat = ",r, "\tReduce f to ",f)
  else:
    print("Histogram not flat enough", r, "\tf left as ", f)

phi=Hopf()

while(f >0.1):
  for s in range(S):
    phi.BFACF()
    if( s%17==0 ):
      counts[ phi.currentHL()]+=1
      G[ phi.currentHL() ]+=f
      H[ phi.currentHL() ]+=1

  flatten()

s=0
while(phi.totalHL()<N*0.75 or phi.deltaHL()>4):
    phi.BFACF(); s+=1
    if( s%17==0 ):
      counts[ phi.currentHL()]+=1
      G[ phi.currentHL() ]+=f
      H[ phi.currentHL() ]+=1

print( phi.halfLength )
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
comps = phi.preparePlot()
ax.plot(comps[0],comps[1],comps[2],'bo-')
ax.plot(comps[3],comps[4],comps[5],'ro-')
ax.set_xlabel('x'); ax.set_ylabel('y');  ax.set_zlabel('z')
plt.show()


fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
x,y = zip(*counts.keys() )
z=counts.values()
w = list(map( lambda X:0.5, x))
b = list(map( lambda X:0, x))
ax.bar3d(x,y,b,w,w,z)
ax.set_xlabel('l1'); ax.set_ylabel('l2')
plt.show()

Here is our sample hopf-link (we tweaked the code to make sure it produced one where the component lengths were similar:

Nicely equal Hopf

The sample histogram isn’t too bad — again this was run only for a few minutes:

Nice Hopf histogram

Flatter-delta to get the Hopf we want

To improve the performance of the above — especially the time taken to difuse over the $(l_1,l_2)$ parameter space — I want to change how we are “bucketing” the data.

For many purposes we are actually just interested in Hopf links in which both components have the same (or simular) size. We don’t really care about links in which one component is big and the other is small. So we don’t really need to collect all possible $(l_1,l_2)$ pairs separately and have the algorithm try to get equal numbers of samples in all those buckets. Instead, lets partition our state-space into total length (or total half-length) and then some simple measure of the difference in lengths between the components. So one could do something like $\begin{align*} (l_1+l_2, |l_1-l_2|) \end{align*}$ but then this still gives $O(n^2)$ buckets. We could also do something like $\begin{align*} (l_1+l_2, \lfloor \sqrt{ |l_1-l_2| } \rfloor) \end{align*}$ so reduce the number of buckets, while still giving us plenty of samples in which the components are similar lengths. But (arguably) even better is $\begin{align*} (l_1+l_2, \lfloor \log(1 + |l_1-l_2| ) \rfloor) \end{align*}$ This gives a $O(n\log n)$ number of buckets for our algorithm to difuse over. We just need to tweak how we store $\mathbb{G}$ and $\mathbb{H}$ .

from collections import defaultdict
from random import randint, random, sample
from math import log,exp,sqrt,floor
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

N=200
S=10**5

#ENUWSD=012345 
#so reverse = dir+/-3
xMove = [1,0,0,-1,0,0]
yMove = [0,1,0,0,-1,0]
zMove = [0,0,1,0,0,-1]

# edge i -> {a,i,b}, so direction of a,b != direction of i or reverse of i.
def UConf(direction):
  # d -> not d or (d+/-3), ie c%3!=d%3
  c = randint(0,5)
  while(c%3==direction%3):
    c = randint(0,5)
  return(c)

def step(p, d):
  return( ( p[0]+xMove[d],p[1]+yMove[d],p[2]+zMove[d]) )

def rstep(p, d):
  return( ( p[0]-xMove[d],p[1]-yMove[d],p[2]-zMove[d]) )

def reverseDirection(d):
  if(d>=3):
    return(d-3)
  else:
    return(d+3)

def diff(p,p0):
  return( (p[0]-p0[0], p[1]-p0[1], p[2]-p0[2]) )



class Hopf:
  def __init__(self):
    self.halfLength=[4,4]
    self.vertices=[set(),set()]
    self.edges=defaultdict(list)

    self.insertVertexEdge(0, (0,0,0), (4, 0));
    self.insertVertexEdge(0, (1,0,0), (0, 0));
    self.insertVertexEdge(0, (2,0,0), (0, 1));
    self.insertVertexEdge(0, (2,1,0), (1, 1));
    self.insertVertexEdge(0, (2,2,0), (1, 3));
    self.insertVertexEdge(0, (1,2,0), (3, 3));
    self.insertVertexEdge(0, (0,2,0), (3, 4));
    self.insertVertexEdge(0, (0,1,0), (4, 4));

    self.insertVertexEdge(1, (1,1,0), (2, 2));
    self.insertVertexEdge(1, (1,1,1), (2, 0));
    self.insertVertexEdge(1, (2,1,1), (0, 0));
    self.insertVertexEdge(1, (3,1,1), (0, 5));
    self.insertVertexEdge(1, (3,1,0), (5, 5));
    self.insertVertexEdge(1, (3,1,-1), (5, 3));
    self.insertVertexEdge(1, (2,1,-1), (3, 3));
    self.insertVertexEdge(1, (1,1,-1), (3, 2));

  def insertVertexEdge(self, component, v, e):
      self.vertices[component].add(v)
      self.edges[v]=e

  def totalHL(self):
      return( self.halfLength[0]+self.halfLength[1] )

  def deltaHL(self):
      return( abs(self.halfLength[0]-self.halfLength[1] ))

  def currentHL(self):
      return( (self.halfLength[0]+self.halfLength[1], floor(log(1+abs(self.halfLength[0]-self.halfLength[1]) ))) )

  def futureHL(self, c, v):
      if(c==0):
        return( (self.halfLength[0]+self.halfLength[1]+v, floor(log(1+abs(self.halfLength[0]+v-self.halfLength[1]) ))) )
      else:
        return( (self.halfLength[0]+self.halfLength[1]+v, floor(log(1+abs(self.halfLength[0]-self.halfLength[1]-v) ))) )

  def preparePlot(self):
    p0 = sample(self.vertices[0],1)[0] # a vertex at random
    p=p0
    xs=[0];ys=[0];zs=[0]
    for k in range(0,self.halfLength[0]*2):
      p=step(p, self.edges[p][1] )
      xs.append(p[0]-p0[0])
      ys.append(p[1]-p0[1])
      zs.append(p[2]-p0[2])

    p=sample(self.vertices[1],1)[0]
    xt=[p[0]-p0[0]];yt=[p[1]-p0[1]];zt=[p[2]-p0[2]]
    for k in range(0,self.halfLength[1]*2):
      p=step(p, self.edges[p][1] )
      xt.append(p[0]-p0[0])
      yt.append(p[1]-p0[1])
      zt.append(p[2]-p0[2])
    return([xs,ys,zs,xt,yt,zt])

  def printV(self, component):
    p0 = sample(self.vertices[component],1)[0] # a vertex at random
    p=p0
    print( "[", diff(p,p0), end='')
    for k in range(0,self.halfLength[component]*2):
      p=step(p, self.edges[p][1] )
      print( ",", diff(p,p0), end='' )
    print("]")

  def printVE(self, component):
    p0 = sample(self.vertices[component],1)[0] # a vertex at random
    p=p0
    print("**********************************")
    print( diff(p,p0), self.edges[p] )
    for k in range(0,self.halfLength[component]*2):
      p=step(p, self.edges[p][1] )
      print( diff(p,p0), self.edges[p] )

  def BFACF(self):
    #pick a vertex from a component 
    if(randint(1,self.halfLength[0]+self.halfLength[1]) <=self.halfLength[0]):
        comp=0
    else:
        comp=1

    x = sample(self.vertices[comp],1)[0] # a vertex at random
    eXtoY = self.edges[x] # edges in and out from it
    y = step(x,eXtoY[1]) # next vertex
    px = rstep(x,eXtoY[0]) # previous vertex
    ny = step( y, self.edges[y][1] ) # next^2 vertex

    #pick direction to push edge x->y
    direction = UConf(eXtoY[1])
    rdirection = reverseDirection(direction)    
    XX = step(x,direction)
    YY = step(y,direction)

    # check if new XX,YY are px or ny
    # use that to compute delta-length
    if(XX==px):
      if(YY==ny):
        delta = -2;
      else:
        delta=0;
    else: # X!=px
      if(YY==ny):
        delta=0; 
      else:
        delta=2;

    if(delta==0):
      if(self.halfLength[comp]==2): 
        return
      if(XX==px): #connect px->YY->y, delete x if legal
        if(YY in self.vertices[0] or YY in self.vertices[1]):
          return #not self-avoiding
        else:
          self.vertices[comp].remove(x); del self.edges[x]
          self.vertices[comp].add(YY)
          self.edges[px]=(self.edges[px][0],eXtoY[1])
          self.edges[YY]=(eXtoY[1], rdirection)
          self.edges[y]=(rdirection,self.edges[y][1])
      else: #must have YY=ny, connect x->XX->ny, delete y if legal
        if(XX in self.vertices[0] or XX in self.vertices[1]):
          return #not self-avoiding
        else:
          self.vertices[comp].remove(y); del self.edges[y]
          self.vertices[comp].add(XX)
          self.edges[x]=(self.edges[x][0],direction)
          self.edges[XX]=(direction, eXtoY[1])
          self.edges[ny]=(eXtoY[1],self.edges[ny][1])
    elif(delta==2):
      if(XX in self.vertices[0] or YY in self.vertices[0] or XX in self.vertices[1] or YY in self.vertices[1] or self.totalHL()>=N):
        return #not self-avoiding
      if( log(random()) > G[self.currentHL()]-G[self.futureHL(comp,1)] ):
        return
      else: #connect x->XX->YY->y
        self.vertices[comp].add(XX); self.vertices[comp].add(YY)
        self.edges[x]=(self.edges[x][0],direction)
        self.edges[XX]=(direction,eXtoY[1])
        self.edges[YY]=(eXtoY[1],rdirection)
        self.edges[y]=(rdirection, self.edges[y][1])
        self.halfLength[comp]+=1
    else: #delta=-2
      if(self.halfLength==2): 
        return
      if( log(random()) > G[self.currentHL()]-G[self.futureHL(comp,-1)] ):
        return
      #connect px->ny, delete x,y
      self.vertices[comp].remove(x); del self.edges[x]
      self.vertices[comp].remove(y); del self.edges[y]
      self.edges[px]=(self.edges[px][0],eXtoY[1])
      self.edges[ny]=(eXtoY[1],self.edges[ny][1])
      self.halfLength[comp]-=1

counts=defaultdict(int) #total samples
H = defaultdict(int) #running histogram
G = defaultdict(float) #G-distribution
f = 2 #G-incrementer

def flatten():
  global f,H,G
  r=min(H.values())/max(H.values())
  if(r>0.5): #time to reduce f
    f=f/2
    H = defaultdict(float) #reset histogram
    g0 = min(G.values())
    G = defaultdict(float, {k: X-g0 for k,X in G.items()} ) #avoid overflows
    print("Histogram is flat = ",r, "\tReduce f to ",f)
  else:
    print("Histogram not flat enough", r, "\tf left as ", f)

phi=Hopf()


while(f > 0.3):
  for s in range(S):
    phi.BFACF()
    if( s%17==0 ):
      counts[ phi.currentHL()]+=1
      G[ phi.currentHL() ]+=f
      H[ phi.currentHL() ]+=1

  flatten()

print( phi.totalHL() )
print( len(phi.vertices[0]), len(phi.vertices[1]) )

s=0
while(phi.totalHL()<N*0.75 or phi.deltaHL()>4 ):
    phi.BFACF(); s+=1
    if( s%17==0 ):
      counts[ phi.currentHL()]+=1
      G[ phi.currentHL() ]+=f
      H[ phi.currentHL() ]+=1

print( phi.totalHL() )
print( len(phi.vertices[0]), len(phi.vertices[1]) )


fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
comps = phi.preparePlot()
ax.plot(comps[0],comps[1],comps[2],'bo-')
ax.plot(comps[3],comps[4],comps[5],'ro-')
ax.set_xlabel('x'); ax.set_ylabel('y');  ax.set_zlabel('z')
plt.show()


fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
x,y = zip(*counts.keys() )
z=counts.values()
w = list(map( lambda X:0.5, x))
b = list(map( lambda X:0, x))
ax.bar3d(x,y,b,w,w,z)
ax.set_xlabel('length'); ax.set_ylabel('log(delta)')
plt.show()

Another sampled Hopf

And a sample-histogram:

Another Hopf histogram