function E = anneal(J)
nsteps = 5000;
rate = 0.001;

nbits = size(J,1);
spins = 2 * (rand(1,nbits) > 0.5) - 1;

for t = 1:nsteps
    E = -(spins*J*spins');
    
    r = ceil(nbits * rand);
    spins(r) = -spins(r); 
    
    En = -(spins*J*spins');
    dE = En - E;
    
    p = exp(-rate * t * dE);
    if ((dE > 0) & (p < rand()))
        pins(r) = -spins(r);
    end;
end
     
E = -(spins*J*spins');