#Nathan Fox
#Homework 27
#I give permission for this work to be posted online

#Read procedures from class/last homework
read(`C27.txt`):

#Import packages
with(`plots`):

#Help procedure
Help:=proc():
    print(` OneStepIsing(M, z) , Ising(m, n, z, K) `):
    print(` DrawMat(M, h) , PatternSearch(m, n, zmin, zmax, zinc, K, h) `):
end:

##PROBLEM 1##

#OneStepIsing(M, z): inputs a matrix M (of size m by n, say) whose
#entries are {1, -1}, picks, uniformaly at random one location
#[i,j] (1<=i<=m, 1<=j<=n) and flips the M[i][j] entry (i.e.
#multiplies it by -1), getting a new matrix (that only differs from
#M in the [i,j] entry) definitely if Gibbs(M',z)/Gibbs(M,z) is >=1,
#and does the flipping with probability Gibbs(M',z)/Gibbs(M,z) if
#it is less 1 (using raR()). Otherwise it does nothing, and goes
#to the next step.
OneStepIsing:=proc(M, z) local i, j, rat, M2:
    i:=rand(1..nops(M))():
    j:=rand(1..nops(M[1]))():
    M2:=M:
    M2[i][j]:=-M[i][j]:
    rat:=Gibbs(M2, z)/Gibbs(M, z):
    if raR() > rat then
        return M:
    else
        return M2:
    fi: 
end:

##PROBLEM 2##

#Ising(m, n, z, K): starts with an m by n {-1,1} matrix taken
#uniformly at random, a parameter z, and a very large positive
#integer K (say 30000), and only outputs the FINAL matrix.
Ising:=proc(m, n, z, K) local M, i:
    M:=RandU(m, n):
    for i from 1 to K do
        M:=OneStepIsing(M, z):
    od:
    return M:
end:

##PROBLEM 3##

#DrawMat(M, h): enters a matrix (given a list of lists) with
#entries in {-1,1}, and "draws" it on a grid with resolution
#"pixel-size" h, where the pixel at location [i,j] is black if
#M[i][j]=1, and white if it is M[i][j]=-1
#
#NOTE: to conform with standard matrix conventions, entry [1,1] is
#in the upper left (as opposed to the lower left, which would
#conform with standard graphics conventions)
DrawMat:=proc(M, h) local i, j, m, ret, cols:
    ret:=NULL:
    cols[1]:=`black`:
    cols[-1]:=`white`:
    m:=nops(M):
    for i from 1 to m do
        for j from 1 to nops(M[i]) do
            ret:=ret, polygonplot(
                [[h*(j-1), h*(m-i)], [h*(j-1), h*(m-i-1)], [h*j, h*(m-i-1)], [h*j, h*(m-i)]],
                color=cols[M[i][j]], filled=true):
        od:
    od:
    return display(ret, axes=none):
end:

##PROBLEM 4##

PatternSearch:=proc(m, n, zmin, zmax, zinc, K, h) local z:
    return display(seq(display(DrawMat(Ising(m, n, z*zinc, K), h), textplot([0, -2*h, sprintf("z=%f", z*zinc)])), z=zmin/zinc..zmax/zinc), insequence=true, axes=none):
end:

#PatternSearch(15, 15, 1/20, 5, 1/20, 30000, 3) took about 20
#minutes.  For small values of z, there's a checkerboard pattern
#(low energy).  As z increases, it heads toward large regions of
#a solid color (high energy). See hw27.gif for an animation.
#PatternSearch(5, 5, 1/20, 5, 1/20, 3000, 3) takes less time.
#For small values of z, it appears to result in a random scattered
#matrix (low energy). For high z values, it appears to result in
#either all white or all black (high energy).  It goes completely
#solid around z=1.6 or so.  See hw27_5.gif for an animation.