# Percolation Model
# Kellen Myers & Bobby DeMarco, 2009
# Released under Kellen's License http://math.rutgers.edu/~kellenm/license.html

Help:=proc(): print(`PerP(M,N,p), States(M) `): 
print(`IsGoodState(M) `):
end:


HelpOld:=proc(): print(` LD(p,x) , GenM(A,B,p), Fatten1(V,S) , Fatten(v,S)`): 
print(`RO(M), MC(A,B,p,K) , Wt(M,p)`):
end:


LD:=proc(p,x) local P,N,L,i,L1,j,r,ra:


if {seq(evalb( coeff(p,x,i)>=0)  ,i=ldegree(p,x)..degree(p,x))}<>
{true}  then
 RETURN(FAIL):
fi:


if {seq(type(coeff(p,x,i),rational)  ,i=ldegree(p,x)..degree(p,x))}<>
{true}  then
 RETURN(FAIL):
fi:

P:=expand(p/subs(x=1,p));

N:=lcm( seq( denom( coeff(P,x,i) ) ,i=ldegree(P,x)..degree(P,x))):


L:=[]:

for i from ldegree(P,x) to degree(P,x) do
 if coeff(P,x,i)<>0 then
   L:=[op(L), [i, N*coeff(P,x,i)]]:
  fi:
od:

ra:=rand(1..N):

L1:=[seq( L[i][2], i=1..nops(L))]:

L1:=[  seq(  add(L1[j],j=1..i)   , i=1..nops(L1))]:

r:=ra():

for j from 1 while L1[j]<r do   od:

L[j][1]:

end:



#GenM(A,B,p): a random A by B 0-1 matrix  with Pr(1)=p
#(indep.)
GenM:=proc(A,B,p) local a,b:
[seq( [ seq( LD(1-p+p*x,x), b=1..B)], a=1..A)]:
end:








#expands S by one 1
Fatten1:=proc(v,S) local i:

for i from 1 to nops(v) do
 if v[i]=1 and (member(i-1,S)  or member( i+1,S)) and not member(i,S)  then
     RETURN(S union {i}):
 fi:
od:

FAIL:

end:

#Fatten(v,S): inputs a vector of 0-1 v, and a set S of indices 
#with S[i]=1, expands it to the set of all 1's reachable from S
#horiz.
Fatten:=proc(v,S) local S1,S2,S3:

S1:=S:

S2:=Fatten1(v,S1):

while S2<>FAIL and S2<>S1 do
S3:=Fatten1(v,S2):

 S1:=S2:
 S2:=S3:

od:

S1:

end:

  #RO(M): the set of places where M[nops(M)] has
#a 1 reachable from M[1]
RO:=proc(M) local m,S,a,M1,S1:
m:=nops(M):

if m=1 then
S:={}:

for a from 1 to nops(M[1]) do
 if M[1][a]=1 then
    S:=S union {a}:
 fi:
od:
RETURN(S):
fi:

M1:=[op(1..m-1,M)]:

S1:=RO(M1):


S:={}:

for a from 1 to nops(M[m]) do
 if M[m][a]=1  and member(a,S1) then
    S:=S union {a}:
 fi:
od:

Fatten(M[m],S):

end:



#IsGood(M): Given a 0-1 matrix decides whether there is a
#"continuous" path of 1 from the first floor to the nops(M)'s floor
IsGood:=proc(M)
evalb(RO(M)<>{}):
end:

#MC(A,B,p,K):runs IsGood(M) K times applied to random A by B 0-1 matrices
MC:=proc(A,B,p,K) local c,M,i:
c:=0:

for i from  1 to K do
 M:=GenM(A,B,p):
 if IsGood(M) then
    c:=c+1:
 fi:
od:

c/K:

end:

#Wt(M,p): the prob. of the 0-1 matrix M
Wt:=proc(M,p) local i,j,c:

c:=1:

for i from 1 to nops(M) do
 for j from 1 to nops(M[i]) do
   if M[i][j]=1 then
      c:=c*p:
   else
       c:=c*(1-p):
   fi:
 od:
od:

c:

end:


###end old stuff
##############start new stuff for April 2, 2009

#PerP(M,N,p): the exact expression (as a poly. in p)
#of the prob. of (simplified) percolation from
#the bottm row to the top row
PerP:=proc(M,N,p) local c,i,v:
c:=0:
for i from 0 to 2^(M*N)-1 do
 v:=convert(i,base,2):
  v:=[op(v),0$(M*N-nops(v))]:
  v:=[seq( [op(N*j+1..N*(j+1),v)],j=0..M-1)]:
 
if IsGood(v) then
  c:=expand(c+Wt(v,p)):
fi:

od:
c:
end:


#States(M): all the vectors of size M
#with {0,1,2} with at least one 2
States:=proc(M) local S,i,v:

S:={}:

for i from 0 to 3^M-1 do
v:=convert(i,base,3):
v:=[op(v),0$(M-nops(v))]:

if IsGoodState(v) then

S:=S union {v}:

fi:

od:

S:

end:

IsGoodState:=proc(v) local i: 

if not member(2,convert(v,set)) then
 RETURN(false):
fi:

for i from 1 to nops(v)-1 do
  if {v[i],v[i+1]}={1,2} then
      RETURN(false):
  fi:
od:

true:
end: