reachability on random diagonal graphs

diagonals

#!/usr/bin/env python

from __future__ import division,print_function
from numpy import *
import scipy.sparse
import scipy.sparse.csgraph

"""
Consider approximate disjoint diagonal cycles concentric around the origin.
We can arrange that cycle C_n has radius O(n) and length O(n).
The probability of all cycles being unbroken is thus roughly

  prod_n (1-2^O(n))

which could easily be positive.  This doesn't answer the question.
"""

random.seed(37832)

def reached(n):
  """How far can we reach on an (2n+1)^2 lattice of random diagonals?"""

  # Build our random graph
  up = random.randint(2,size=4*n*n).reshape(2*n,2*n).astype(bool)
  def V(x,y):
    return (2*n+1)*x+y
  x = arange(2*n)
  x,y = x[:,None],x[None,:]
  i = V(x,y+1-up)
  j = V(x+1,y+up)
  v = ones((2*n)**2,dtype=bool)
  A = scipy.sparse.coo_matrix((v,(i.ravel(),j.ravel())),shape=((2*n+1)**2,)*2).tocsr()

  # Perform depth first search
  reached = scipy.sparse.csgraph.depth_first_order(A,i_start=V(n,n),directed=False,return_predecessors=False)
  good = zeros((2*n+1,2*n+1),dtype=bool)
  good.ravel()[reached] = 1
  return good

def plot():
  import matplotlib.pyplot as plt
  from mpl_toolkits.mplot3d import Axes3D
  from matplotlib import cm
  fig = plt.figure()

  steps = 1000
  ns = array([25,50,100,200,400,800])
  for sub,n in enumerate(ns):
    ax = fig.add_subplot(1,len(ns),sub+1,projection='3d')
    prob = zeros((2*n+1,2*n+1))
    for s in xrange(steps):
      print('step %d / %d'%(s,steps))
      prob += reached(n)
    prob /= steps
    print('center prob = %g'%(prob[n,n]))
    print('prob = \n%s'%prob)
    i = arange(2*n+1)-n
    i,j = meshgrid(i,i)
    def half(x):
      assert x.shape==(2*n+1,2*n+1)
      y = zeros_like(x)
      for i in xrange(2*n+1):
        l = min(i,2*n-i)
        j = arange(2*l+1)-l
        y[i,n+j] = x[i+j,i-j]
      return y
    print('half(prob) = \n%s'%half(prob))
    if n<=50:
      ax.plot_surface(i,j,half(prob),rstride=1,cstride=1,linewidth=0,cmap=cm.coolwarm)
    else:
      ax.plot_surface(i,j,half(prob),linewidth=0,cmap=cm.coolwarm)
    ax.set_title('n = %d'%n)
    ax.view_init(elev=5.)
    ax.w_xaxis.set_ticklabels([])
    ax.w_yaxis.set_ticklabels([])
  plt.show()

if __name__=='__main__':
  plot()