% Application of simulated annealing to our "funny hamiltonian"
clc

%rand('state',0);  % Restart the random generator if you want
k = 0; % Number of permutations found
N = 12; %Number of sites in the model
Jmatrix = initialize_hamiltonian(N); % Set up the interactions
Nocc = ceil(N/2); % Number of occupied sites
Nperms = nchoosek(N,Nocc);
fprintf('Search space: %d\n',Nperms)
energies = zeros(1,Nperms);
for i =1:2^N % In this loop is a simple-minded, brute-force method
             % for generating all possible permutations of 0's and 1's
    binstr = dec2bin(i); 
    occupation = str2num(binstr');
    occupation = [zeros(N-length(occupation),1); occupation];
    if Nocc~=sum(occupation)
        continue
    end
    occupation;
    k = k + 1;
    energies(k) = compute_energy(Jmatrix,occupation');
end

oldguess = [zeros(N-Nocc,1); ones(Nocc,1)];
k = 0; %iteration counter
kT = 6; % Temperature
dt = .99;

energy = compute_energy(Jmatrix,oldguess');
%while kT>1e-5
while k < 500
    site1=ceil(rand(1)*N); % Select two particle
    site2=ceil(rand(1)*N); % locations to switch
    while oldguess(site1)==oldguess(site2)  % Keep selecting sites till we've selected both an empty and occupied site
        site2=ceil(rand(1)*N);
    end
    newguess = oldguess;
    temp = oldguess(site1);
    oldguess(site1) = oldguess(site2);
    oldguess(site2) = temp;
    deltaE = compute_energy(Jmatrix,oldguess') - energy;
    if deltaE < 0
        newguess = oldguess;
        %disp('jump down')
        
    elseif rand(1) < exp(-deltaE/kT)
        newguess = oldguess;
        %disp('jump up')
    else
        oldguess = newguess;
        %disp('no jump')
    end
    k = k + 1;
    energy = compute_energy(Jmatrix,newguess');
    MCe(k) = energy;
    kT = kT*dt;
    newguess';
end
plot(MCe,'.-')
%min(energies)
%energy
fprintf('Global min: %15.7f\n Found min: %15.7f\n',min(energies),energy)
if (min(energies)<energy)
    disp('Not converged')
else
    disp('Converged')
end