%
% opt.m
% Neil Gershenfeld  (c) 11/18/95
%
clear all
nbits = 100;
bonds = randn(1,nbits);
min_energy = -sum(abs(bonds));

%
% genetic algorithms solution
%
nsteps = 300;
population = 100;
beta = [10 1 0.1 0.01];
for temp = 1:length(beta)
   spins = [];
   energy = [];
   newspins = [];
   for n = 1:population
      spins = [spins; (2 * (rand(1,nbits) > 0.5) - 1)];
      end
   for i = 1:nsteps
      %
      % evaluate fitness
      %
      for n = 1:population
         energy(i,n) = - sum(bonds .* (spins(n,:) .* ...
            [spins(n,(2:nbits)) spins(n,1)]));
         end
      prob = exp(-beta(temp)*(energy(i,:)-min_energy));
      prob = prob/sum(prob);
      prob = cumsum(prob);
      %
      % reproduce
      %
      for n = 1:population
         %
         % draw from distribution
         %
         index1 = 1+sum((rand) > prob);
         index2 = 1+sum((rand) > prob);
         %
         % crossover
         %
         crossover = ceil(rand*nbits);
         newspins(n,(1:crossover)) = spins(index1,(1:crossover));
         newspins(n,((crossover+1):nbits)) = ...
            spins(index2,((crossover+1):nbits));
         %
         % mutate
         %
         bit = ceil(nbits * rand);
         newspins(n,bit) = -newspins(n,bit);
         end
      spins = newspins;
      end
   plot(min(energy'),'b')
   hold on
   end
xlabel('step')
ylabel('energy')
hold on
plot([0 nsteps],[min_energy min_energy],'r')
hold off
print -deps optganew.eps

return

%
% simulated annealing solution
%
nsteps = 5000;
nrates = 5;
rates = [0.0001 0.001 0.01 0.1];
starting_spins = 2 * (rand(1,nbits) > 0.5) - 1;
for rate = 1:length(rates)
   spins = starting_spins;
   for step = 1:nsteps
      old_energy = - sum(bonds .* (spins .* [spins(2:nbits) spins(1)]));
      energy(step,rate) = old_energy;
      %
      % flip a random bit
      %
      bit = ceil(nbits * rand);
      spins(bit) = -spins(bit);
      new_energy = - sum(bonds .* (spins .* [spins(2:nbits) spins(1)]));
      delta_energy = new_energy - old_energy;
      %
      % check to see if we need to reject it
      %
      if ((delta_energy > 0) & (rand(1) > exp(-delta_energy*rates(rate)*step)))
         spins(bit) = -spins(bit);
         end;
      end
   end
plot(energy,'b')
xlabel('step')
ylabel('energy')
hold on
plot([0 nsteps],[min_energy min_energy],'r')
hold off
print -deps optsanew.eps
figure