function [en_hist, it_count] = ga_minimize(H, Npop, mut_r, max_its, BFmin)
% GA_MINIMIZE attempts to find a minimum using a genetic algorithm
% H is the Hamiltonian, i.e, the coupling matrix (diagonal entries are
% on-site energies
% indv is a 2D array, each row is an individual, a string of 1's and 0's

en_hist = ones(max_its,Npop)+BFmin; % History of energies over the iteraction of the GA
                               % initialize to something larger than BFmin
it_count = 1;                  % Counter for number of generations (or iterations)

N = length(H);                 % Number of sites in the lattice
indv = ceil(rand(Npop,N)-.5);  % Initialize the array, random 0's & 1's
ga_ens = zeros(1,Npop);        % Energies of the GA's current population
fitness = zeros(Npop);         % Fitness scores of the GA's current population

while (1) % Break condition for the while loop is at the end of the loop

    % Compute energy for each individual in the population
    for i = 1:Npop
        ga_ens(i) = compute_energy(H,indv(i,:));
    end
	minE = min(ga_ens);
    en_hist(it_count,:) = ga_ens; 
    
    % Determine fitness scores for each individual
    for i = 1:Npop
        fitness(i) = 1/(1+ga_ens(i)-minE);
	end
	
	% Order the individuals and their fitness scores
	[dummy, order] = sort(ga_ens);  % Somewhat obtuse but this code section
	indv = indv(order,:);           % just orders the individuals and their
	fitness = fitness(order);       % scores in ascending order of energy
	
	% This block of code peforms the mating. Individuals with better
    % fitness scores are selected prefentially to be parents. Mating
    % is done using a simple one-point crossover
	new_indv = indv; % Make a copy before mating
	fitcum = cumsum(fitness)/sum(fitness); % This array is used in a "weighted" selection
                                           % of parents that favors those with better
                                           % fitness scores
	for i = 1:Npop  % Loop over each individual, make one new child for each individual
         rnd = rand(1);
         [a,p1] = max(rnd<fitcum); % Select the first parent, p1
         rnd = rand(1);
         [a,p2] = max(rnd<fitcum); % Select the first parent, p2
         split = ceil((N-2)*rand(1)+1); % select the place to "split," i.e., the crossover point
         new_indv(i,:)=[indv(p1,1:split-1) indv(p2,split:end)]; % Splice together the new individual
	end
	
	% Make random mutations 
	for i = 1:Npop
        muts = [ceil(rand(1,N)-1+mut_r)];  % Make an array of sites to mutate
        new_indv(i,:) = mod(muts+new_indv(i,:),2);  % At chosen mutation sites, switch 
	end                                             % 1's to 0's and 0's to 1's
	
	indv=new_indv;  % Copy children to current population 

    % Exit the function if we have found the minimum or reached the 
    % maximum number of allowed iterations (generations)
    if  it_count == max_its || any(en_hist(it_count,:) == BFmin) 
        break
    end
    it_count = it_count + 1;
end
