#OK to post homework
#Yukun Yao, Mar. 20, Assignment 15



###Problem 1### Read and understand the three new procedures in C15.txt , and experiment with them. Try to do: NickAv(4,4, 1.3 , 1.0, 1000,10000, E, E,M); ten times, and see whether they are close (when I tried it, it was a little disappointing)

#It returns [10.66840000, 9.506400000, 11.22280000, 11.39240000, 10.48120000, 10.26200000, 10.61120000, 10.66120000, 10.79320000, 11.77320000], which looks good. 



###Problem 2### Now try: NickAv(40,40, 1.3 , 1.0, 1000,10000, E, E,M); ten times. Are the outputs closer to each other?

#[876.3900000, 970.6956000, 917.2284000, 882.6540000, 886.8880000, 946.5564000, 910.4008000, 839.3156000, 882.1044000, 936.0376000], it seems not too bad. 



###Problem 3### According to Julia Yeomans's book starting with either the all 1-matrix or the alternationg 1,-1, rather than a random matrix gives better results. Write procedures NickAvONE(M,N,x,y,K1,K2,G,En,Ma), NickAvALT(M,N,x,y,K1,K2,G,En,Ma): that instead of starting with a random matrix, starts with the all 1s and the 1,-1,... matrices respectively, and see if you get better agreements between the various runs.

NickAvONE:=proc(M,N,x,y,K1,K2,G,En,Ma) local A,i,j,su:
A:=[seq([seq(1, j=1..N)], i=1..M)]:

for i from 1 to K1 do
 A:=OSM(A,M,N,x,y):
od:

su:=0:

for i from 1 to K2 do
 A:=OSM(A,M,N,x,y):
 su:=su+subs({En=Ene(A),Ma=Mag(A)},G):
od:

evalf(su/K2):
end:


NickAvALT:=proc(M,N,x,y,K1,K2,G,En,Ma) local A,i,j,su:
A:=[seq([seq((-1)^(i+j), j=1..N)], i=1..M)]:

for i from 1 to K1 do
 A:=OSM(A,M,N,x,y):
od:

su:=0:

for i from 1 to K2 do
 A:=OSM(A,M,N,x,y):
 su:=su+subs({En=Ene(A),Ma=Mag(A)},G):
od:

evalf(su/K2):
end:


#[seq(NickAvONE(4, 4, 1.3, 1.0, 1000, 10000, E, E, M), i = 1 .. 10)] returns [11.60360000, 10.66880000, 9.498000000, 11.21480000, 11.37640000, 10.54240000, 10.26000000, 10.63280000, 10.59840000, 10.81560000]

#[seq(NickAvALT(4, 4, 1.3, 1.0, 1000, 10000, E, E, M), i = 1 .. 10)] returns [11.77320000, 11.29480000, 11.14360000, 11.17080000, 10.74880000, 10.16000000, 10.75280000, 10.62120000, 10.60400000, 11.40440000]

#[seq(NickAvONE(40, 40, 1.3, 1.0, 1000, 10000, E, E, M), i = 1 .. 10)] returns [1340.871600, 1352.812000, 1375.066400, 1381.382800, 1321.116800, 1264.835600, 1384.162400, 1399.124400, 1390.297600, 1464.889200]

#[seq(NickAvALT40,40, 1.3 , 1.0, 1000,10000, E, E,M), i=1..10 )] returns [853.3140000, 830.5640000, 906.7880000, 811.3176000, 870.3404000, 834.2068000, 879.1888000, 808.0080000, 865.3692000, 848.1572000]



###Problem 4### Write a procdeure Shrink(A) that starts with a 3k by 3k matrix of {1,-1} and outputs a 3k-1 by 3k-1 matrix, that returns the "majority rule" in each 3 by 3 consecutive square. Run Nick(81,81,x,y,10000): for various numerical values of x (not too far from 1), then shrink them successively, and use PlotMat(A) to plot them. Try to reproduce the kind of pictures in Kenneth G. Wilson's seminal article

#The following one shrinks 3k by 3k to k by k. Actually I am not sure whether you meant shrink to 3k-2 by 3k-2 (but in this case we don’t need the dimensional to be a multiply of 3) or k by k. 

Shrink:=proc(A) local m, n, B, i, j, i1, j1, temp:
m:=nops(A):
n:=nops(A[1]):
if m mod 3 <> 0 or n mod 3 <> 0 then
  return fail:
fi:
B:=[seq([seq(0, j=1..n/3)], i=1..m/3)]:
for i from 1 to m/3 do
  for j from 1 to n/3 do
    temp:=add(add(A[i1][j1], j1=3*j-2..3*j), i1=3*i-2..3*i):
    if temp > 0 then
      B[i][j]:=1:
    else
      B[i][j]:=-1:
    fi:
  od:
od:
B:
end:


#I attach the picture of shrinked Nick(300, 300, 1.002, 1.0011, 10000) with this homework. 










###### Following is from C15.txt ######



#March 12, 2020 "Virtual Class". Corrected version of March 13, 2020 
#C15.txt: Using  Metropolis to estimate macroscopic quantities (i.e. expectations) of physical quantities
#and comparing them for the actual values for small lattices
#The new procedures for C15.txt are listed in Help(). To see the procedures done in Classes 13 and 14
#type: Help13(); and  Help14();
Help:=proc(): print(` NickAv(M,N,x,y,K1,K2,G,Eng,Mag), ExactAv(M,N,x,y,G,En,Ma)  , ExactAvBrute(M,N,x,y,G,En,Ma) `): end:



#ExactAvBrute(M,N,x,y,G,En,Ma):   Like ExactAv(M,N,x,y,G,En,Ma) but by brute force
ExactAvBrute:=proc(M,N,x,y,G,En,Ma) local S,tot,su, A:
S:=Mats(M,N):

tot:=0:
su:=0:
for A in S do
tot:=evalf(tot+x^Ene(A)*y^Mag(A)):
su:=evalf(su+x^Ene(A)*y^Mag(A)*subs({En=Ene(A),Ma=Mag(A)},G)):
od:

su/tot:
end:


#ExactAv(M,N,x,y,G,En,Ma):  Corrected version of March 13, 2020
#Inputs positive integers M and N and real numbers x and y, positive integers K1 and K2
#and an expression G that is a MONOMIAL in the symbols E and M (denoting energy and magnetization, respectively), outputs the
#the EXACT average of G(En,Ma) over all 2^(M*N) possible M by N {1,-1} matrices with the Ising prob. distribution with
#temperature parameter x and external field parameter y. It uses "generating-functionlogy". If the
#Prob. gen. fun of a r.v. is f(z) then its expectation is (z*d/dz f(z))|z=1 and its r-th moment is
#(zd/z)^r) f(z) |z=1 , similarly for more than variable
#For example, to the exact average of the energy with x=1.1 and y=1 on the 4 by 4 rectangular torus do:
#ExactAv(4,4, 1.1, 1.0, En,En,Ma). 
#To get the exact average of En^2*Ma^2 do
#ExactAv(4,4, 1.1, 1.0, En^2*Ma^2,En,Ma). 
ExactAv:=proc(M,N,x,y,G,En,Ma) local d1,d2,C,i,f,g,X,Y:

#We check that G is indeed a monomial in En and Ma, i.e. is of the form C*En^d1*Ma^d2 for some number C 

d1:=degree(G,En):
d2:=degree(G,Ma):
C:=normal(G/(En^d1*Ma^d2)):

if not type(C,numeric) then
 print(G, `is not a monomial in`, En, Ma ):
 RETURN(FAIL):
fi:

#Introduce FORMAL variables X and Y to represent the energy and magnetization
#this is needed for the "generating-functionlogy" part, we need x and y to be numeric
#to reflect the probabilty distribution
#f is the exact expression in X and Y for the so-called (unnormalized) partition function for the M by N torodial rectangle
f:=ParF(M,N,X,Y):

#Let's normalize f
f:=f/subs({X=1,Y=1},f):

#Now let's use the numbers x and y to make it a numeric prob. generating function where the coeff. of the monomial X^a*Y^b
#is the exact probability that the energy of the spin configuration is a and its magnetization is b
f:=subs({X=x*X,y=y*Y},f):

#Now let's have a copy of f, let's call it g
g:=f:
#Now let's apply (X*d/dX) d1 times followed by (Y*d/dY) d2 times

for i from 1 to d1 do
 g:=expand(X*diff(g,X)):
od:


for i from 1 to d2 do
 g:=expand(Y*diff(g,Y)):
od:


#now let's plug-in X=1, Y=1 to get a number that only depends on the physical parameters x and y.
#Note that x and y can be left symbolic, if desired (but for comparing it to the output of NickAv, they have to be numeric)
#we also have to multiply by the coefficient of the monomial
#and divide by the numerator
C*subs({X=1,Y=1},g)/subs({X=1,Y=1},f):
end:

#####Old erroneous code
#ExactAvNoGood(M,N,x,y,G,En,Ma):  Inputs positive integers M and N and real numbers x and y, positive integers K1 and K2
#and an expression G that is a MONOMIAL in the symbols E and M (denoting energy and magnetization, respectively), outputs the
#the EXACT average of G(En,Ma) over all 2^(M*N) possible M by N {1,-1} matrices with the Ising prob. distribution with
#temperature parameter x and external field parameter y. It uses "generating-functionlogy". If the
#Prob. gen. fun of a r.v. is f(z) then its expectation is (z*d/dz f(z))|z=1 and its r-th moment is
#(zd/z)^r) f(z) |z=1 , similarly for more than variable
#For example, to the exact average of the energy with x=1.1 and y=1 on the 4 by 4 rectangular torus do:
#ExactAv(4,4, 1.1, 1.0, En,En,Ma). 
#To get the exact average of En^2*Ma^2 do
#ExactAvNoGood(4,4, 1.1, 1.0, En^2*Ma^2,En,Ma). 
ExactAvNoGood:=proc(M,N,x,y,G,En,Ma) local d1,d2,C,i,f,X,Y:

#We check that G is indeed a monomial in En and Ma, i.e. is of the form C*En^d1*Ma^d2 for some number C 

d1:=degree(G,En):
d2:=degree(G,Ma):
C:=normal(G/(En^d1*Ma^d2)):

if not type(C,numeric) then
 print(G, `is not a monomial in`, En, Ma ):
 RETURN(FAIL):
fi:

#Introduce FORMAL variables X and Y to represent the energy and magnetization
#this is needed for the "generating-functionlogy" part, we need x and y to be numeric
#to reflect the probabilty distribution
#f is the exact expression in X and Y for the so-called (unnormalized) partition function for the M by N torodial rectangle
f:=ParF(M,N,X,Y):

#Let's normalize f
f:=f/subs({X=1,Y=1},f):

#Now let's use the numbers x and y to make it a numeric prob. generating function where the coeff. of the monomial X^a*Y^b
#is the exact probability that the energy of the spin configuration is a and its magnetization is b
f:=subs({X=x*X,y=y*Y},f):

#Now let's apply (X*d/dX) d1 times followed by (Y*d/dY) d2 times

for i from 1 to d1 do
 f:=expand(X*diff(f,X)):
od:


for i from 1 to d2 do
 f:=expand(Y*diff(f,Y)):
od:


#now let's plug-in X=1, Y=1 to get a number that only depends on the physical parameters x and y.
#Note that x and y can be left symbolic, if desired (but for comparing it to the output of NickAv, they have to be numeric)
#we also have to multiply by the coefficient of the monomial
C*expand(subs({X=1,Y=1},f)):
end:




#NickAv(M,N,x,y,K1,K2,G,En,Ma):  Inputs positive integers M and N and real numbers x and y, positive integers K1 and K2
#and an expression G in the symbols E and M (denoting energy and magnetization, respectively), outputs the
#estimate produced by running the Metropolis algorithm  K1 in the "warm-up" stage (w/o taking records)
#followed by K2 iterations where a record of G(E,M) is taken, and added to the average. 
#For example, to get an estimate for the average of the energy-squared, try:
#NickAv(4,4,2. ,1.1, 100,1000, E^2, E,M);

NickAv:=proc(M,N,x,y,K1,K2,G,En,Ma) local A,i,su:
A:=RM(M,N):

for i from 1 to K1 do
 A:=OSM(A,M,N,x,y):
od:

su:=0:

for i from 1 to K2 do
 A:=OSM(A,M,N,x,y):
 su:=su+subs({En=Ene(A),Ma=Mag(A)},G):
od:

evalf(su/K2):
end:


with(Statistics):

###stuff from C13.txt
Help13:=proc(): print(`RM(M,N), PlotMat(A), Ene(A) , Mag(A)`):
print(`Prob(A,x,y) , Vecs(N) , Mats(M,N) ,  ParF(M,N,x,y)`):
end:

RM:=proc(M,N) local i,j,ra:
ra:=rand(0..1):
subs(0=-1,[seq([seq(ra(),j=1..N)],i=1..M)]):
end:

with(plots):

#PlotMat(A): plots the matrix A where 1 is black and -1 is white
PlotMat:=proc(A) local i,j,d:

if A[1][1]=1 then
d:=pointplot([[1,1]],axes=none,color=black,symbol=circle):
else
d:=pointplot([[1,1]],axes=none,color=red,symbol=circle):
fi:

for j from 2 to nops(A[1]) do
if A[1][j]=1 then
d:=d,pointplot([[1,j]],color=black,symbol=circle):
else
d:=d,pointplot([[1,j]],color=red,symbol=circle):
fi:
od:

for i from 2 to nops(A) do
for j from 1 to nops(A[i]) do
if A[i][j]=1 then
d:=d,pointplot([[i,j]],color=black,symbol=circle):
else
d:=d,pointplot([[i,j]],color=red,symbol=circle):
fi:
od:
od:



display(d);

end:


#Ene(A): The energy of the matrix A of {1,-1}
Ene:=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][j]*A[i+1][j]   , j=1..nops(A[1]) )  ,i=1..nops(A)-1)+
add(A[1][i]*A[-1][i],i=1..nops(A[1]))+
add(A[i][1]*A[i][-1],i=1..nops(A)):
end:

#Mag(A): the magnetization of the spin-config. A
Mag:=proc(A): add(add(A)): end:

#Vecs(N): the set of all 2^N {1,-1} lists of length N
Vecs:=proc(N) local S,v:
option remember:
if N=0 then
 RETURN({[]}):
fi:
S:=Vecs(N-1):

{seq([op(v),1], v in S),seq([op(v),-1], v in S)}:

end:




#Mats(M,N): the set of all M by N {1,-1} matrices
Mats:=proc(M,N) local S,S1,A,v:
if M=0 then
 RETURN({[]}):
fi:

S:=Mats(M-1,N):
S1:=Vecs(N):

{seq(seq([op(A),v],v in S1),A in S)}:

end:

#x=exp(k/T) (T=abs. temp. (in Kelvin)) y=exp(External Magnetication)
#x^Ene*y^Mag
#ProbA(A,x,y):Unnormalied prob.
ProbA:=proc(A,x,y): x^Ene(A)*y^Mag(A):end:

#ParF(M,N,x,y): the Gibbs partition function of M by N Ising model
ParF:=proc(M,N,x,y) local S,A:
S:=Mats(M,N):
add(ProbA(A,x,y),A in S):
end:

###end from C13.txt



###start stuff from C14.txt

Help14:=proc(): print(` Neis(c,M,N) , NewA(A,c), OSM(A,M,N,x,y), Nick(M,N,x,y,K1)`): end:

#Neis(c,M,N): Inputs a location c=[i,j] and pos. integers M and N outputs the SET
#of four neighbors of the cell [i,j], {[i \pm 1,j ], [i, j \pm 1] but with toral convention
#CODE by Alison Bu and Yukun Yao, with improvements by Charles Kenney
Neis:=proc(c,M,N) local i,j:

i:=c[1]: j:=c[2]:

{[(i-2 mod M)+1,j ],[(i mod M)+1, j ],[i, (j-2 mod N)+1 ] ,[i,(j mod N)+1]}:


end:

#NewA(A,c): The result of flipping the spin at location c in the matrix A
#followed by the change in Energy, and change in Magnetization
NewA:=proc(A,c) local i,j,s, A1,Ne,DiffE, DiffM,a:
i:=c[1]: j:=c[2]:
s:=A[i][j]:

A1:=[op(1..i-1,A),
 [op(1..j-1,A[i]), -s, op(j+1..nops(A[i]), A[i])],  
 op(i+1..nops(A),A)]:

Ne:=Neis(c,nops(A),nops(A[1])):

DiffE:=-2*s*add(A[a[1]][a[2]], a in Ne):
DiffM:=-2*s:
[A1,DiffE,DiffM]:
end:

#OSM(A,M,N,x,y): one step in the Metropolis algorithm inputs an M by N matrix of {1,-1}
#and pos. numbers x,y outputs one step in the Metropolis algorithm according to
#the importance x^Ene(A)*y^Mag(A)
OSM:=proc(A,M,N,x,y) local c,Zidong, A1,Rat,r,U:
c:=[rand(1..M)(), rand(1..N)()]:

Zidong:=NewA(A,c):
A1:=Zidong[1]:

Rat:=x^Zidong[2]*y^Zidong[3]:

r:=Sample(RandomVariable(Uniform(0,1)),1)[1]:

if r<Rat then
  A1:
else
 A:
fi:

end:

#Nick(M,N,x,y,K1): Starts with a random M by N {1,-1} matrix with physical parameters x and y
#and  a positive integer K1 runs Metropolis K1 times
Nick:=proc(M,N,x,y,K1) local A,i:
A:=RM(M,N):

for i from 1 to K1 do
 A:=OSM(A,M,N,x,y):
od:
A:
end:

#####end stuff from C14.txt