% % 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