% SpeedyGA is a vectorized implementation of a Simple Genetic Algorithm in Matlab
% version 1.2.1
% Copyright (C) 2007,2008  Keki Burjorjee

%  Licensed under the Apache License, Version 2.0 (the "License"); you may
%  not use this file except in compliance with the License. You may obtain 
%  a copy of the License at  
%
%  http://www.apache.org/licenses/LICENSE-2.0 
%
%  Unless required by applicable law or agreed to in writing, software 
%  distributed under the License is distributed on an "AS IS" BASIS, 
%  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 
%  See the License for the specific language governing permissions and 
%  limitations under the License. 

% Created and tested under Matlab 7 (R14). 


% If a parameter is passed to the speedyGA function, it is used to set the
% seed of the random number generator. Otherwise, the seed is made to be
% dependent upon the system clock. 
% Sample usage: speedyGA(), speedyGA(27463)

function speedyGA(varargin)
len=640                 % The length of bitstrings
popSize=500             % The size of the population
maxGens=500             % The maximum number of generations allowed for a run
probCrossover=1         % The probability of crossing over. 
                        % (Note: 1-probCrossover is the probability of
                        % reproduction without crossover)
probMutation=1/popSize % The mutation probability (per bit)
sigmaScalingFlag=1      % Sigma Scaling is described on pg 168 of M. Mitchell's
                        % GA book. It often improves GA performance.
sigmaScalingCoeff=1     % Higher values => less fitness pressure 
crossoverType=1         % 0 => no crossover
                        % 1 => 1pt crossover
                        % 2 => uniform crossover
visualizationFlag=1     % 0 => don't visualize bit frequencies
                        % 1 => visualize bit frequencies
makeVideoFlag=1                        
aviFileName=['speedyGArun_crossoverType=' num2str(crossoverType)]; 

%pregenerate crossover and mutation masks
maskReposFactor=5;
crossmasks=[];
if crossoverType==0
    crossmasks=false(popSize,len);
elseif crossoverType==2
    crossmasks=rand(popSize,(len+1)*maskReposFactor)<0.5;
end
mutmasks=rand(popSize,(len+1)*maskReposFactor)<probMutation;

% preallocate vectors for recording the average and maximum fitness in each
% generation
avgFitnessHist=zeros(1,maxGens+1);
maxFitnessHist=zeros(1,maxGens+1);

% If a parameter was passed to the vectorGA function, it is used to set the
% seed of the random number generator. Otherwise, the seed is made to be
% dependent upon the system clock
if nargin==1
    rand('state',varargin{:});
else
    rand('state',sum(100*clock));
end
    
% the population is a popSize by len matrix of randomly generated boolean
% values
pop=rand(popSize,len)<.5;
if visualizationFlag==1
    figure(1);
    set(gcf, 'position', [360 173 624 527])
end
if visualizationFlag==1 & makeVideoFlag==1
    delete(aviFileName);
    aviObj=avifile(aviFileName);
    aviObj.Quality=100;
    aviObj.fps=8;
    aviObj.compression='Cinepak';
end
for gen=0:maxGens
        
     fitnessVals=royalRoads(pop);
    
    %obtain the fitness value and index of the fittest individual 
    [maxFitnessHist(gen+1), maxIndex]=max(fitnessVals);
    
    %%%%%%%%%%%%%%% display population information
    % create a string of the bestIndividual
    bestIndiv=[];
    for i=1:len
        bestIndiv=[bestIndiv num2Str(pop(maxIndex,i))];
    end
    avgFitnessHist(gen+1)=mean(fitnessVals);
    
    % display the generation number, the average Fitness of the population,
    % the maximum fitness of an individual in the population, and an 
    % individual with the maximum fitness
    display(['gen=' num2str(gen,'%.3d') '   avgFitness=' ...
        num2str(avgFitnessHist(gen+1),'%3.3f') '   maxFitness=' ...
        num2str(maxFitnessHist(gen+1),'%.3d') '  bestIndiv=' bestIndiv]);
    
    % Conditionally perform bit-frequency visualization
    if (visualizationFlag==1)
        fig=figure(1)
        bitFreqs=sum(pop)/popSize;
        plot(1:len,bitFreqs, '.');
        axis([0 len 0 1]);
        title([' \newline \newline Generation = ' num2str(gen) ', Average Fitness = ' sprintf('%0.3f', avgFitnessHist(gen+1))]);
        ylabel('Frequency of the bit 1');
        xlabel('Locus \newline \newline ');
        drawnow;
        if makeVideoFlag==1
            F=getFrame(fig);
            aviObj=addFrame(aviObj, F);
        end
    end    

    % if sigma scaling is on then perform it
    if(sigmaScalingFlag)
        sigma=std(fitnessVals);
        if sigma~=0;
            fitnessVals=1+(fitnessVals-mean(fitnessVals))/...
            (sigmaScalingCoeff*sigma);
            fitnessVals(fitnessVals<=0)=realmin;
        else
            fitnessVals=1;
        end
    end        
    
    
    % Normalize the fitness values and then create an array with the 
    % cumulative normalized fitness values (the last value in this array
    % will be 1)
    cumNormFitnessVals=cumsum(fitnessVals/sum(fitnessVals));
    
    % Preallocate an array for the post-selection population 
    newPop=zeros(popSize,len);
    % Perform fitness proportional selection using Stochastic Universal
    % Sampling
    SUSMarkers=rand(1,1)+[1:popSize]/popSize;
    gt1=SUSMarkers>1;
    SUSMarkers(gt1)=SUSMarkers(gt1)-1;
    [temp newPopIndices]=histc(SUSMarkers,[0 cumNormFitnessVals]);
    pop=pop(newPopIndices,:);

    %%%%%%%%%%%%%%% Perform uniform recombination and reproduction    
    % generate pairs of mating indices
    matingIndices=ceil(rand(popSize,2)*popSize);
    % create a new population using the first mate from each mating pair
    newPop=pop(matingIndices(:,1),:);
    % create a new population using the second mate from each mating pair
    secondMate(:,:)=pop(matingIndices(:,2),:);
    % create crossover masks for the entire population
    if crossoverType==0
        masks=crossmasks;
    elseif crossoverType==1
        masks=zeros(popSize, len);
        temp=ceil(rand(popSize,1)*(len-1));
        for i=1:popSize
            masks(i,1:temp(i))=1;
        end
    else    
        temp=floor(rand*len*(maskReposFactor-1));
        masks=crossmasks(:,temp+1:temp+len);
    end
    % generate an array of booleans indicating which first mates to leave 
    % untouched
    reprodIndices=rand(popSize,1)<1-probCrossover;
    masks(reprodIndices,:)=false;
    % implement uniform recombination with reproduction
    newPop(logical(masks))=secondMate(logical(masks));
    % implement mutation
    temp=floor(rand*len*(maskReposFactor-1));
    masks=mutmasks(:,temp+1:temp+len);
    pop=xor(newPop,masks);    
end
figure(2)
hold off
plot([0:maxGens],avgFitnessHist,'k-');
hold on
plot([0:maxGens],maxFitnessHist,'c-');
title('Average and Maximum Fitness Per Generation')
xlabel('Generations')
ylabel('Fitness')
if makeVideoFlag==1
    aviObj=close(aviObj)
end

% onemax
function fitness=oneMax(pop)
fitness=sum(pop,2)';    

% The royal roads function. The chromosome length (i.e. len) 
% should be a multiple of 8
function fitness=royalRoads(pop)
[popSize len]=size(pop);
fitness=zeros(popSize,1);
for i=1:8:len
    temp=sum(pop(:,i:i+7),2);
    temp=double(temp==8);
    fitness=fitness+temp;
end
fitness=fitness';