function minE = genetic(N, numiter, betas, myrnd, population)

% This function simulates a cooling 1D spin glass using a genetic algorithm
% Each spin S1...SN can be either +1 or -1
% J is a random variable drawn from the Gaussian distribution with zero mean and unit variance
% We compare the searching algorithm for different cooling rates
% Input: 
% 	N = size of the spin glass
%	numiter = the number of iterations
% 	rates = a vector of cooling rates to compare
%	myrnd = a first random value for the random generator
%	population = number of sets of spin glasses
% Output: 
%	E = Energy curve of the spin glass over time

% Create the random variable J of size N
J = normrnd(0,1,1,N);

% Create sets spins using our 'Random Generator' from a previous class
for p = 1:population,
	Spins(p,:) = randomGenerator(N,myrnd+p);
	Spins(p,:) = ((Spins(p,:) > 0.5) - 0.5) * 2;
end

% calculate the minimum possible Energy
Emin = -sum(abs(J));

for i = 1:length(betas),

	S = Spins;
	
	for t = 1:numiter,

		% Calculate the Energy for each set 
		for p = 1:population,
			E(p) = - (sum(J .* S(p,:) .* [S(p,2:N) S(p,1)]));
		end
	
		% Calculate the Probability for each set 
		probability = exp(-betas(i) * (E-Emin));

		% Normalize the sum of probabilities to 1
		probability = probability/sum(probability);
		
		% Reorganize the data from 0 to 1
		probability = cumsum(probability);
		
		for p = 1:population,

			% Pick up 2 sets out ot a hundred
			set1 = sum(rand > probability) + 1;
			set2 = sum(rand > probability) + 1;
			
			% find a random crossover
			crossover = ceil(rand * (population-1));

			% Couple the 2 set an create a new one
			newS(p,:) = [S(set1,(1:crossover)) S(set2,(crossover+1:N))];
		
			% Flip one bit randomly (Mutation)
			flipIndex = ceil(rand(1) * N);
			newS(p,flipIndex) = -newS(p,flipIndex);
		end
		
		S = newS;
		Energy(t) = min(E);
	end

	% Plot
	plot(Energy,'color',[0 (1-(1/(log(i)+1))) (1/(log(i)+1))]);
	hold on;
	
end

plot([1 numiter],[Emin Emin],'color',[1 0.4 0.4]);
hold off;

xlabel('Steps');
ylabel('Energy');

