Help:=proc(): print(` Wt(A,z,w) , GenM(m,n,z,w,K), OneStep(m,n,A,z,w)`): end:



HelpOld:=proc(): print(`En(A), PF(m,n,z) , RealAve1(m,n,z,w,p,q)`): 
print(`RealAve2(m,n,z,w,p1,q1,p2,q2), Cor(m,n,z,w,p1,q1,p2,q2) `):
end:

#En(A): the sum of the products of all nearest-neighbors of 
#the matrix A given as a list of lists
En:=proc(A) local i,j:

 add(add(A[i][j]*A[i][j+1],j=1..nops(A[i])-1),i=1..nops(A))+

add(add(A[i-1][j]*A[i][j],i=2..nops(A)),j=1..nops(A[1])):

end:

#PF(m,n,z)
PF:=proc(m,n,z) local s,c,ab:

s:=0:

for c from 0 to 2^(m*n)-1 do
 
 ab:=convert(c,base,2):
 ab:=[op(ab),0$(m*n-nops(ab))]:
 ab:=subs(0=-1,ab):
 ab:=[seq([op((i-1)*n+1..i*n,ab)],i=1..m)]:
  
 s:=s+z^En(ab):
od:

s:

end:


#Ex: evaluate f(z)=lim (m,n->infinity) PF(m,n,z)^(1/(m*n))
#(equivalenty: limit of log(PF(m,n,z))/(m*n) )


#gPF(m,n,z,w)
gPF:=proc(m,n,z,w) local s,c,ab,ab1:

s:=0:

for c from 0 to 2^(m*n)-1 do
 
 ab:=convert(c,base,2):
 ab:=[op(ab),0$(m*n-nops(ab))]:
 ab:=subs(0=-1,ab):
 ab1:=[seq([op((i-1)*n+1..i*n,ab)],i=1..m)]:
  
 s:=s+z^En(ab1)*w^convert(ab,`+`):
od:

s:

end:

#Ex: evaluate g(z,w)=lim (m,n->infinity) PF(m,n,z,w)^(1/(m*n))
#(equivalenty: limit of log(PF(m,n,z,w))/(m*n) )


#RealAve1(m,n,z,w,p,q): the (weighted) average of A[i][j]
#amongsts all m by n matrices with {-1,1} with Boltzmann weight
RealAve1:=proc(m,n,z,w,p,q) local t,s,ab,ab1,i,j,c:

s:=0:

for c from 0 to 2^(m*n)-1 do
 
 ab:=convert(c,base,2):
 ab:=[op(ab),0$(m*n-nops(ab))]:
 ab:=subs(0=-1,ab):
 ab1:=[seq([op((i-1)*n+1..i*n,ab)],i=1..m)]:
  
 s:=s+z^En(ab1)*w^convert(ab,`+`)*ab1[p][q]:
od:

s/gPF(m,n,z,w):

end:




#RealAve2(m,n,z,w,p1,q1,p2,q2): the (weighted) average of A[p1][q1]*A[p2][q2]
#amongsts all m by n matrices with {-1,1} with Boltzmann weight
RealAve2:=proc(m,n,z,w,p1,q1,p2,q2) local t,s,ab,ab1,i,j,c:

s:=0:

for c from 0 to 2^(m*n)-1 do
 
 ab:=convert(c,base,2):
 ab:=[op(ab),0$(m*n-nops(ab))]:
 ab:=subs(0=-1,ab):
 ab1:=[seq([op((i-1)*n+1..i*n,ab)],i=1..m)]:
  
 s:=s+z^En(ab1)*w^convert(ab,`+`)*ab1[p1][q1]*ab1[p2][q2]:
od:

s/gPF(m,n,z,w):

end:

#E(XY)-E(X)E(Y)=E((X-E(X))(Y-E(Y))
Cor:=proc(m,n,z,w,p1,q1,p2,q2):

RealAve2(m,n,z,w,p1,q1,p2,q2)-
RealAve1(m,n,z,w,p1,q1)*RealAve1(m,n,z,w,p2,q2):

end:







#Wt(A,z,w)
Wt:=proc(A,z,w) local s,c,ab,ab1: 
  
 z^En(A)*w^convert(convert(A,`+`),`+`)

end:


#OneStep(m,n,A,z,w): performs one step in the Metropolis arlgorithm
#applied to the Ising model with parameters z and w (importatnce given
#by z^En(A)*w^convert(convert(A,`+`),`+`)
OneStep:=proc(m,n,A,z,w) local i,j,A1,rat,dennis:

i:=rand(1..m)():

j:=rand(1..n)():

A1:=A:

A1[i][j]:=-A[i][j]:


rat:=Wt(A1,z,w)/Wt(A,z,w):

if rat>=1 then
  RETURN(A1):
fi:

dennis:=stats[random,uniform[0,1]](1):

if dennis<rat then
 RETURN(A1):
else
 RETURN(A):
fi:

end:


#GenM(m,n,z,w,K): running Metropolis to generate "typical" spin-systems
#(matrices in {-1,1}) with the parameters z exp(Constant/T) and w (something
#to do with magenetization)
GenM:=proc(m,n,z,w,K) local A,A0,ra,i,j,bn:
ra:=rand(0..1):

A0:=[seq([seq( 2*ra()-1,i=1..n)],j=1..m)]:

for bn from 1 to K do
A0:=OneStep(m,n,A0,z,w):
od:

A0:


end: