from numarray import *
from MLab import *
from pylab import *

def energy(spins, bonds):
    n = size(spins)
    e = bonds[:n-1] * spins[:n-1] * spins[1:]
    return - (sum(e) + bonds[n-1] * spins[n-1] * spins[0])

def anneal(spins, bonds, alpha, steps = 100):
    spins = spins.copy()
    n = size(spins)
    e = zeros((steps,1), Float64);
    for t in range(steps):
        # Energy
        E = energy(spins,bonds)

        # Flip a random spin
        r = int(n*rand())
        spins[r] = - spins[r]

        # New Energy and Delta
        En = energy(spins,bonds)
        dE = En - E
        e[t] = E

        # Reject Move?
        p = exp(-alpha * t * dE)
        if ((dE > 0) & (p < rand())):
            spins[r] = - spins[r] # Reject
    return e

def main() :
    steps = 6000
    alphas = array([0.1, 0.01, 0.001])
    N = 100
    spins = 2 * (rand(N,1) > 0.5) - 1
    bonds = randn(N,1)
    for i in range(size(alphas)):
        e = anneal(spins, bonds, alphas[i], steps)
        plot(e)
    title("Simulated Annealing");
    show()
    

if __name__ == "__main__" : main()