# OK to post homework
# Jingyang Deng, Mar.12th, 2020, Assignment 15

# There's something wrong with my Maple and I could not use "subs", 
# so I used "algsubs" instead.

Help := proc(): print(`ExactAvNew(M, N, x, y, G, En, Ma), NickAvONE(M, N, x, y, K1, K2, G, En, Ma), NickAvALT(M, N, x, y, K1, K2, G, En, Ma), Shrink(A)`): end:
with(Statistics):
with(plots):
with(MTM):

####  prob.1  ####

# I run [seq(NickAv(4, 4, 1.3, 1.0, 1000, 10000, E, E, M), i = 1 .. 10)]; and got
# [10.65760000, 10.07760000, 11.07760000, 11.86680000, 9.775200000, 
#  10.80520000, 11.07200000, 10.44240000, 10.43080000, 10.82320000]
# In fact, it's not that close.


####  prob.2  ####

# Run ExactAv(4, 4, 1.3, 1.0, E, E, M); and the output is 37.27622212. It's far from answers above.
# The code is erroneous because we can not use ParF(M, N, X, Y) to add all the probability
# together directly. The calculation of pgf in "ExactAv" is wrong.
# A correct version of the procedure ExactAvNew(M, N, x, y, G, En, Ma) is as below:

ExactAvNew := proc(M, N, x, y, G, En, Ma) local d1, d2, C, S, A, pgf, X, Y, i, j:
d1, d2 := degree(G, En), 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:
S := Mats(M, N):
pgf := Array([seq(ProbA(A, x, y) * X^Ene(A) * Y^Mag(A), A in S)]):
pgf := add(pgf / add(algsubs(X = 1, algsubs(Y = 1, pgf)))):
for i from 1 to d1+d2 do
    if i<=d1 then
        pgf := X * diff(pgf, X):
    else
        pgf := Y * diff(pgf, Y):
    fi:
od:
C * algsubs(X = 1, algsubs(Y = 1, pgf)):
end:

# Here's a test:
# ExactAvNew(4, 4, 1.3, 1.0, E, E, M);
#                           10.79301894
# It seems that this is the right answer.

# A := [seq(NickAv(4, 4, 1.3, 1.0, 1000, 10000, E, E, M), i = 1 .. 100)];
# Mean(A);
#                         10.7736000000000
# StandardDeviation(A);
#                        0.579936991352848
# Answers derived by Metropolis algorithm are around exact answer but not close to it.
# However, as long as the sample population is large enough, the approximation is also good.


####  prob.3  ####

# Run [seq(NickAv(40, 40, 1.3, 1.0, 1000, 10000, E, E, M), i = 1 .. 10)]; and I got
# [864.6096000, 866.5960000, 874.2252000, 898.3248000, 942.3916000, 
#  902.4960000, 904.1988000, 816.5968000, 875.4144000, 895.9156000]
# It's still not that good, mainly because the population is not large enough.


####  prob.4  ####

#the initial state is ones(M, N)
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 + algsubs(En = Ene(A), algsubs(Ma = Mag(A), G)):
od:
evalf(su / K2):
end:

#the initial state is a alternate matrix
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 + algsubs(En = Ene(A), algsubs(Ma = Mag(A), G)):
od:
evalf(su / K2):
end:

# Here's a test, which shows that the initial state doesn't matter a lot.
# A := [seq(NickAvONE(4, 4, 1.3, 1.0, 1000, 10000, E, E, M), i = 1 .. 100)];
# Mean(A);
#                         10.7867120000000
# StandardDeviation(A);
#                        0.565804919316925
# A := [seq(NickAvALT(4, 4, 1.3, 1.0, 1000, 10000, E, E, M), i = 1 .. 100)];
# Mean(A);
#                         10.7236040000000
# StandardDeviation(A);
#                        0.589484679579195


####  prob.5  ####

Shrink := proc(A) local M, N, count, tmp, S, B, i, j:
B := Array(A):
M, N := seq(size(B)):
if (not M=N) or (not type(log[3](M), integer)) then RETURN(FAIL): fi:
#count(C): inputs a matrix C, outputs the dominant element (1 or -1) in C.
count := proc(C):
    if add(add(C)) > 0 then RETURN(1): else RETURN(-1): fi:
end proc:
S := [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
        tmp := [ [B[3*i-2, 3*j-2], B[3*i-2, 3*j-1], B[3*i-2, 3*j]],
                      [B[3*i-1, 3*j-2], B[3*i-1, 3*j-1], B[3*i-1, 3*j]],
                      [B[3*i, 3*j-2], B[3*i, 3*j-1], B[3*i, 3*j]] ]:
    S[i][j] := count(tmp):
    od:
od:
S:
end:

#generate a matrix whose size is 3^lg3size, shrink it susseccively and plot it.
plt := proc(lg3size, x, y) local A:
A := Nick(3^lg3size, 3^lg3size, x, y, 10000):
print(listdensityplot(A, axes = none)):
do
    A := Shrink(A):
    print(listdensityplot(A, axes = none)):
until nops(A)=1:
end:

####  code from class (revised)  ####

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

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:=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:=proc(A): add(add(A)): end:

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:=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:

ProbA:=proc(A,x,y): x^Ene(A)*y^Mag(A):end:

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:=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:=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:=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:

ExactAv:=proc(M,N,x,y,G,En,Ma) local d1,d2,C,i,f,X,Y:
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:
f:=ParF(M,N,X,Y):
f:=f/subs({X=1,Y=1},f):
f:=subs({X=x*X,y=y*Y},f):
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:
C*expand(algsubs(X=1,algsubs(Y=1,f))):
end:

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+algsubs(En=Ene(A), algsubs(Ma=Mag(A),G)):
od:
evalf(su/K2):
end: