function picoGA(varargin)
LEN=60  %the length of bitstrings
POPSIZE=200 % the size of the population
GENERATIONS=100 % the maximum number of generations allowed for a run
CROSSOVER_PROB=.8% the probability of crossing over
PMUT=.5/LEN %the mutation probability per bit

% If a parameter was passed to the TinyGA 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 size(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=generateRandomBooleanMatrix(POPSIZE,LEN,.5);

for gen=1:GENERATIONS
    
    fitnessVals=zeros(1,POPSIZE);
    %calculate fitness values of chromosomes in the population
    for i=1:POPSIZE
        for j=1:LEN
            fitnessVals(i)=fitnessVals(i)+double(pop(i,j));
        end        
    end
    
    %obtain the fitness value and index of the fittest individual 
    [maxFitness, maxIndex]=max(fitnessVals);
    
    %%%%%%%%%%%%%%% display population information
    % create a string of the bestIndividual
    bestIndiv=[];
    for i=1:LEN
        bestIndiv=[bestIndiv num2Str(pop(maxIndex,i))];
    end
    
    % 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(mean(fitnessVals),'%3.3f') '   maxFitness=' ...
        num2str(maxFitness,'%.3d') '  bestIndiv=' bestIndiv]);
    if sum(fitnessVals==LEN)>0
        display('Maximum Found!');
        return;
    end
    
    %%%%%%%%%%%%%%%% Perform Fitness Proportional Selection
    %fitnessVals=fitnessVals-min(fitnessVals)+realmin;
    
    % 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=false(POPSIZE,LEN);
    % perform fitness proportional selection over the current population
    for i=1:POPSIZE
        k=rand(1,1);
        
        % calculate the index of the the chromosome (in the old population) 
        % which will be added to the new population.
        % index is an array of POPSIZE+1 booleans with one and only one 
        % true value (the rest are all false). This true value is at the 
        % index of the individual 
        % which needs to be added to the newPop
        index=[cumNormFitnessVals>=k false ]&[ true  cumNormFitnessVals<k];
        
        %Uses logical indexing to add the individual to the newPop
        newPop(i,:)=pop(index(1:POPSIZE),:);
    end
    
    pop=newPop;
    %%%%%%%%%%%%%%% 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 uniform crossover masks for the entire population
    masks=generateRandomBooleanMatrix(POPSIZE,LEN,.5);
    % generate an array of booleans indicating which first mates to leave 
    % untouched
    reprodIndices=generateRandomBooleanMatrix(POPSIZE,1,1-CROSSOVER_PROB);
    masks(reprodIndices,:)=false;
    % implement uniform recombination with reproduction
    newPop(logical(masks))=secondMate(logical(masks));
    % implement mutation
    masks=generateRandomBooleanMatrix(POPSIZE,LEN,PMUT);
    pop=xor(newPop,masks);    
end
display(['Maximum not found after ' num2str(GENERATIONS) ' generations']);
    
% returns rows x cols array of booleans in which the probability of a 1 is
% pTrue
function randomBooleanMatrix=generateRandomBooleanMatrix(rows,cols,pTrue)
for i=1:rows
    for j=1:cols
        randomBooleanMatrix(i,j)=rand(1,1)<pTrue;
    end
end