Help:=proc(): print(`LC(p), Rect(M,N),RP(M,N,p), Neig(G,v), CC(G,v) `): end:

#LC(p): Given a prob. p (a rational number) outputs 1 with prob. p
#and 0 with prob. 1-p
LC:=proc(p) local a,b,ra:
a:=numer(p): b:=denom(p):

ra:=rand(1..b)():

if ra<=a then
 RETURN(1):
else
 RETURN(0):
fi:

end:

#Rect(M,N): the set of edges of the rectangle [0,M]x[0,N]
#For example, Rect(1,1) = { {[0,0],[1,0]},{[0,0],[0,1]}, {[1,0],[1,1]}, 
#{[0,1],[1,1]}}

Rect:=proc(M,N) local i,j:
option remember:

{seq(seq([[i,j],[i+1,j]],i=0..M-1),j=0..N),

seq(seq([[i,j],[i,j+1]],i=0..M),j=0..N-1)}:

end:

with(plots):
P:=proc(G): plot(G,axes=none):end:

#beautiful code by Dennis Hou
RP:=proc(M,N,p) local i,j,S:
S:={seq(seq(LC(p)*[[i,j],[i+1,j]],i=0..M-1),j=0..N),
seq(seq(LC(p)*[[i,j],[i,j+1]],i=0..M),j=0..N-1)}:
S minus {[[0,0],[0,0]]}:
end:


#Neig(G,v): the set of (up to four neighbors) of vertex v in G
Neig:=proc(G,v) local N,n:

N:={v+[1,0],v+[-1,0],v+[0,1],v+[0,-1]}:

for n in N do
 if  not( member([v,n],G) or member([n,v],G)) then
   N:=N minus {n}:
 fi:
od:
      
N union {v}:

end:


CC:=proc(G,v) local N1,N2,N3,n:
N1:={v}:

N2:=Neig(G,v):

while N1<>N2 do
 
N3:={ seq( op(Neig(G,n)) , n in N2) }:

N1:=N2:
N2:=N3:

od:

N2:

end: