Metropolis Algorithm

import numpy
import pylab

size = 16
n = 2**size

T = numpy.arange(0.2,5,.1)

B = 1/(T)

isingarray = numpy.zeros((n,size))
Hvals = numpy.zeros(n)
m = n

for col in range (0,size,1):
    m = m/2
    boolean = True
    for row in range (0,n,1):
        if row%m == 0:
            boolean = not boolean
        if boolean == True:
            isingarray[row,col] = 1
        else:
            isingarray[row,col] = -1

asize = numpy.int(numpy.sqrt(size))

H = numpy.zeros(n)

Z = numpy.zeros(B.size)
ZHsum = numpy.zeros(B.size)
ZH2sum = numpy.zeros(B.size)

avgH = numpy.zeros(B.size)
avgH2 = numpy.zeros(B.size)

for row in range (0,n,1):
    testarray = numpy.reshape(isingarray[row,:], (asize,asize))
    for i in range (0,asize,1):
        H[row] += -(testarray[0,i]*testarray[asize-1,i])
        H[row] += -(testarray[i,0]*testarray[i,asize-1])


A = numpy.eye(asize**2, asize**2,1)
A[numpy.arange(asize**2) % asize == asize-1] = 0
A += numpy.eye(asize**2, asize**2, asize)

H += -(numpy.dot(isingarray,A)*isingarray).sum(axis=1)

for i in range (0,B.size):
    for k in range (0,H.size):
        Z[i] += numpy.exp(-B[i]*H[k])
        ZHsum[i] += numpy.exp(-B[i]*H[k])*H[k]
        ZH2sum[i] += numpy.exp(-B[i]*H[k])*(H[k])**2

    avgH[i] = ZHsum[i]/Z[i]
    avgH2[i] = ZH2sum[i]/Z[i]

C = (avgH2-avgH**2)/((T**2)*size)

#Begin Metropolis Algorithm

metroarray = numpy.zeros(B.size)


for j in range(0,B.size,1):

    metro = numpy.reshape(isingarray[numpy.random.randint(0,65537),:], (asize,asize))

    metrosize = numpy.zeros(10000)
    waitperiod = metrosize.size*1.5

    for count in range (numpy.int(metrosize.size-waitperiod),metrosize.size,1):

        mH = 0.
        mHp = 0.

        for i in range (0,asize,1):
            mH += (metro[0,i]*metro[asize-1,i])
            mH += (metro[i,0]*metro[i,asize-1])
            for k in range (0, asize-1,1):
                mH += (metro[i,k]*metro[i,k+1])
                mH += (metro[k,i]*metro[k+1,i])
        mH = mH * -1

        metro2 = metro

        choice = numpy.random.randint(0,4,1)
        choice2 = numpy.random.randint(0,4,1)

        metro2[choice,choice2] = metro2[choice,choice2]*-1

        for i in range (0,asize,1):
            mHp += (metro2[0,i]*metro2[asize-1,i])
            mHp += (metro2[i,0]*metro2[i,asize-1])
            for k in range (0, asize-1,1):
                mHp += (metro2[i,k]*metro2[i,k+1])
                mHp += (metro2[k,i]*metro2[k+1,i])
        mHp = mHp * -1

        if mHp < mH:
            metro = metro2
        if numpy.random.rand() < numpy.exp(-B[j]*(mH-mHp)):
            metro = metro2

        mH = 0.

        for i in range (0,asize,1):
            mH += (metro[0,i]*metro[asize-1,i])
            mH += (metro[i,0]*metro[i,asize-1])
            for k in range (0, asize-1,1):
                mH += (metro[i,k]*metro[i,k+1])
                mH += (metro[k,i]*metro[k+1,i])
        mH = -1*mH
        if count >= 0:
            metrosize[count] += mH
            print count
    metroarray[j] = numpy.sum(metrosize)/metrosize.size