function minE = annealing(N, numiter, rates, myrnd)

% This function simulates a cooling 1D spin glass using the annealing optimisation technique
% 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
% 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 the spins using our 'Random Generator' from a previous class
% We choose 'startwith' randomly
Spins = randomGenerator(N,myrnd);
Spins = ((Spins > 0.5) - 0.5) * 2;

% Set the vector of Energies
E = zeros(1,numiter);

for i = 1:length(rates),

	rate = rates(i);
	S = Spins;
	
	% Calculate the Energy for one step 
	E(1) = - (sum(J .* S .* [S(2:N) S(1)]));

	for t = 2:numiter,

		% Choose a spin randomly
		flipIndex = ceil(rand(1) * N);
	
		% Flip the spin from -1 to 1 or 1 to -1
		S(flipIndex) = -S(flipIndex);
	
		% Calculate the new Energy for one step 
		E(t) = - (sum(J .* S .* [S(2:N) S(1)]));

		% Compare newE to previous E
		DeltaE = E(t)-E(t-1);
	
		% Move with probability if Energy increased
		% Move always if Celta Energy decreased
		if (DeltaE > 0)
			probability = exp(-rate * t * DeltaE);
		else
			probability = 1;
		end
	
		% Reverse previous flip and keep the same E if rejected
		if (rand(1) > probability)
			S(flipIndex) = -S(flipIndex);
			E(t) = E(t-1);
		end
	
	end

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

% calculate the minimum possible Energy
minE = -sum(abs(J));
plot([1 numiter],[minE minE],'color',[1 0.4 0.4]);

hold off;

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

