% Funny hamiltonian for practicing searches (simulated annealing, genetic
% algorithm, etc.)
clc

%rand('state',0);  % Restart the random number generator if you wish
k = 0;  % Number of permutations found so far
N = 16; % Number of sites in the model
Nocc = ceil(N/2); % Number of occupied sites

Jmatrix = initialize_hamiltonian(N); % Set up the interactions
Nperms = nchoosek(N,Nocc);  % How many possible permutations
fprintf('Search space: %d\n',Nperms)
energies = zeros(1,Nperms);  % Initialize the array for energies
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);           % Convert the number to binary
    occupation = str2num(binstr'); % Convert it to a string
    occupation = [zeros(N-length(occupation),1); occupation]; % Pad with extra zeroes
    if Nocc~=sum(occupation) % Skip over cases with the
        continue             % wrong number of occupied sites
    end
    k = k + 1;
    energies(k) = compute_energy(Jmatrix,occupation');
end

hist(energies,20)