# Newman-Ziff percolation code

from random import *
from math import *
import Numeric 
from visual import *
from visual.graph import *

L=8
N=(L*L)
EMPTY=(-N-1)
p=0.5
seed(22222)

scene=display(x=0,y=400,width=400,height=400,fov=1e-8,background=color.white)

pic1=gdisplay(x=0,y=0,width=400,height=400,title="Fraction of sites in largest cluster vs p",
             xtitle="p",ytitle="Fraction of sites in largest cluster",xmax=1,xmin=0,ymin=0,ymax=1,foreground=color.black,background=color.white)
            
plot1=gcurve(gdisplay=pic1,color=color.blue)


def findroot(i):
#    print "%4d" % (ptr[i])
    if ptr[i]<0:
        return i
    else:
        ptr[i] = findroot(ptr[i])
        return ptr[i]


nn = Numeric.zeros((N,4))   
for i in range (N):
    xpos=i%L
    ypos=int(i/L)
    nn[i,0]=i-1+L*(xpos==0)
    nn[i,1]=i+1-L*(xpos==(L-1))
    nn[i,2]=i-L+N*(ypos==0)
    nn[i,3]=i+L-N*(ypos==(L-1))
   
order=Numeric.zeros((N))
for i in range (N):
    order[i] = i
for i in range (N):
    j=int(i+(N-i)*random())
    temp = order[i]
    order[i] = order[j]
    order[j] = temp
#    print "%4d\t%4d" % (i,order[i])
big =0
ptr=Numeric.zeros((N))
for i in range (N):
    ptr[i] = EMPTY
for i in range (N):
    r1 = order[i]
    s1 = order[i]
    ptr[s1] = -1
    for j in range (4): 
        s2 = nn[s1][j]
        if ptr[s2]!=EMPTY:
           r2 = findroot(s2)
           if r2!=r1: 
                if ptr[r1]>ptr[r2]:
                    ptr[r2] += ptr[r1]
                    ptr[r1] = r2
                    r1 = r2
                else:
                    ptr[r1] += ptr[r2]
                    ptr[r2] = r1
	#   print "%4d\t%4d" % (ptr[r1],order[i]))
                if -ptr[r1]>big:
                    big = -ptr[r1]
    print "%4d\t%4d" % (i+1,big)
    plot1.plot(pos=( float(i+1)/float(N) , float(big)/float(N) ))

for i in range (N):
    if i<int(p*N):
        ptr[order[i]] = -1
    else:
        ptr[order[i]]=EMPTY
for k in range (2):    
    for i in range (N):
        if ptr[i]!=EMPTY: 
            s1=i
            r1=findroot(s1);
            for j in range (4):
                s2=nn[s1,j]
                if ptr[s2]!=EMPTY:
                    r2 = findroot(s2)
                    if (r2!=r1):
                        if (ptr[r1]>ptr[r2]):
                            ptr[r2] += ptr[r1]
                            ptr[r1] = r2
                            r1 = r2
                        else:                        
                            ptr[r1] += ptr[r2]
                            ptr[r2] = r1

scene.background=color.white                            

max=0
for i in range(N):
    if EMPTY<ptr[i]<0:
        if i>max:
            max=i
box_size=1.1
for i in range (N):
#    print "%4d\t%4d" % (i,ptr[i])
    xpos=i%L
    ypos=(i/L)
#    print "%4d\t%4d" % (xpos,ypos)
    if ptr[i]>=0:
        box(pos=(box_size*xpos-L/2,box_size*ypos-L/2,0),color=((1.0*ptr[i])/(max+1),(1.0*ptr[i])/(max+1),(1.0*ptr[i])/(max+1)))
    if EMPTY<ptr[i]<0:
        box(pos=(box_size*xpos-L/2,box_size*ypos-L/2,0),color=((1.0*i)/(max+1),(1.0*i)/(max+1),(1.0*i)/(max+1)))