Skip to content

Instantly share code, notes, and snippets.

@girving
Created Aug 20, 2015
Embed
What would you like to do?
reachability on random diagonal graphs
#!/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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment