function BinaryGA

% Binary genetic algorithm, based on Goldberg's book.

% maximize 2^10
popsize = 10; % population size
lchrom = 30; % length of each chromosome (number of bits)
numGenes = 1; % number of genes per chromosome
pcross = 0.6; % probability of crossover
pmutation = 0.033; % probability of mutation (per bit)

% Ackley test function
% popsize = 50; % population size
% numGenes = 20; % number of genes per chromosome
% lchrom = 9 * numGenes; % 9 bits per gene, 20 genes per chromosome
% pcross = 1.0; % probability of crossover
% pmutation = 0.0; % probability of mutation (per bit)

maxgen = 15; % max # of generations to run GA

chromosome = zeros(lchrom, 1);
individual = struct('chrom',chromosome, ... % chromosome = genotype = bit string
    'x',zeros(1,numGenes), ... % phenotype = parameters
    'fitness',0, ... % objective function value
    'parent1',0, ... % first parent
    'parent2',0, ... % second parent
    'xsite',0); % crossover point

% Define "oldpop" and "newpop" as an array of individuals
clear oldpop newpop
for i = 1 : popsize
    oldpop(i) = individual;
    newpop(i) = individual;
end

% initialize the population
[oldpop] = initpop(lchrom, popsize, numGenes);
[sumfitness, minfitness, maxfitness, avgfitness] = statistics(oldpop);
initreport(popsize, lchrom, maxgen, pcross, pmutation, maxfitness, avgfitness, minfitness);
maxfitnessArr = [maxfitness]; % array of maximum fitness values per generation
avgfitnessArr = [avgfitness]; % array of average fitness values per generation
% run the GA for the specified max # of generations
for gen = 1 : maxgen
    [newpop] = generation(oldpop, sumfitness, pcross, lchrom, pmutation, numGenes); % select, crossover, mutate
    [sumfitness, minfitness, maxfitness, avgfitness] = statistics(newpop); % collect statistics
    maxfitnessArr = [maxfitnessArr maxfitness];
    avgfitnessArr = [avgfitnessArr avgfitness];
    report(gen, maxfitness, avgfitness, minfitness); % output to the screen
    oldpop = newpop; % replace the old population with the new population
end

% plot results
close all;
gen = 0 : maxgen;
plot(gen,maxfitnessArr, gen,avgfitnessArr,'--');
legend('max fitness', 'avg fitness');
xlabel('generation');
ylabel('fitness');
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [obj] = objfunc(x, lchrom, numGenes)
% fitness function evaluation of x^10
coef = 2^lchrom - 1;
n = 10;
obj = (x / coef)^n;

% Compute the Ackley cost function
% sum1 = 0;
% sum2 = 0;
% for numVar = 1 : numGenes
%     sum1 = sum1 + x(numVar)^2;
%     sum2 = sum2 + cos(2*pi*x(numVar));
% end
% obj = 20 + exp(1) - 20 * exp(-0.2*sqrt(sum1/numGenes)) - exp(sum2/numGenes);
% obj = 25 - obj;

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function value = decode(chrom, numGenes)
% decode a chromosome into its parameters (phenotypes)

% x^10 test function
value = 0;
powerof2 = 1;
for j = 1 : length(chrom)
    if chrom(j) == 1
        value = value + powerof2;
    end
    powerof2 = 2 * powerof2;
end

% Ackley - map the numVar genes in the chromosome to numVar values between -30 and +30
% BitsPerGene = length(chrom) / numGenes; % number of bits per gene
% ndx = 0;
% for numVar = 1 : numGenes % each gene is mapped to a phenotype
%     value(numVar) = 0;
%     for j = 1 : BitsPerGene
%         ndx = ndx + 1;
%         value(numVar) = value(numVar) + chrom(ndx) * 2^(BitsPerGene-j);
%     end
%     value(numVar) = -30 + 60 * value(numVar) / (2^BitsPerGene - 1);
% end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [oldpop] = initpop(lchrom, popsize, numGenes)
% initialize a random population
for j = 1 : popsize
    for i = 1 : lchrom
        if rand < 0.5
            oldpop(j).chrom(i) = 0;
        else
            oldpop(j).chrom(i) = 1;
        end
    end
    oldpop(j).x = decode(oldpop(j).chrom, numGenes);
    oldpop(j).fitness = objfunc(oldpop(j).x, lchrom, numGenes);
    oldpop(j).parent1 = 0;
    oldpop(j).parent2 = 0;
    oldpop(j).xsite = 0;
end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [individual] = select(pop, sumfitness)
% Select a single individual via roulette wheel selection
partsum = 0;
j = 0;
randomnumber = rand * sumfitness;
while (partsum <= randomnumber) & (j <= length(pop))
    j = j + 1;
    partsum = partsum + pop(j).fitness;
end
individual = j;
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [child1, child2, jcross] = crossover(parent1, parent2, pcross, lchrom, pmutation)
% Cross two parent chromosomes to get two child chromosomes
if rand < pcross 
    jcross = floor(lchrom * rand + 1 - eps); % random # betwen 1 and lchrom
else
    jcross = lchrom;
end
% first exchange: 1 to 1, and 2 to 2
child1 = zeros(lchrom, 1);
child2 = zeros(lchrom, 1);
for j = 1 : jcross
    child1(j) = mutation(parent1(j), pmutation);
    child2(j) = mutation(parent2(j), pmutation);
end
% second exchange: 1 to 2, and 2 to 1
for j = jcross + 1 : lchrom
    child1(j) = mutation(parent2(j), pmutation);
    child2(j) = mutation(parent1(j), pmutation);
end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [allele] = mutation(alleleval, pmutation)
% mutate a bit with some specified probability
if rand < pmutation
    allele = ~alleleval;
else
    allele = alleleval;
end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [newpop] = generation(oldpop, sumfitness, pcross, lchrom, pmutation, numGenes)
% select, crossover, and mate until the entire population is replaced
for j = 1 : 2 : length(oldpop)
    % pick a pair of mates
    mate1 = select(oldpop, sumfitness);
    mate2 = select(oldpop, sumfitness);
    % crossover and mutate
    [newpop(j).chrom, newpop(j+1).chrom, jcross] = ...
        crossover(oldpop(mate1).chrom, oldpop(mate2).chrom, pcross, lchrom, pmutation);
    % decode, evaluate fitness, and record parentage
    for i = j : j + 1
        newpop(i).x = decode(newpop(i).chrom, numGenes);
        newpop(i).fitness = objfunc(newpop(i).x, lchrom, numGenes);
        newpop(i).parent1 = mate1;
        newpop(i).parent2 = mate2;
        newpop(i).xsite = jcross;
    end
end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [sumfitness, minfitness, maxfitness, avgfitness] = statistics(pop)
% compute the statistics of the population
sumfitness = 0;
minfitness = inf;
maxfitness = -inf;
for j = 1 : length(pop)
    sumfitness = sumfitness + pop(j).fitness;
    minfitness = min(minfitness, pop(j).fitness);
    maxfitness = max(maxfitness, pop(j).fitness);
end
avgfitness = sumfitness / length(pop);
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function initreport(popsize, lchrom, maxgen, pcross, pmutation, maxfitness, avgfitness, minfitness)
% output an initial report to the screen
disp('    Simple Genetic Algorithm Parameters');
disp('    -----------------------------------');
disp(['   Population size (popsize)         =   ', num2str(popsize)]);
disp(['   Chromosome length (lchrom)        =   ', num2str(lchrom)]);
disp(['   Maximum # of generations (maxgen) =   ', num2str(maxgen)]);
disp(['   Crossover probability (pcross)    =   ', num2str(pcross)]);
disp(['   Mutation probability (pmutation)  =   ', num2str(pmutation)]);
disp(' ');
disp('    Initial Generation Statistics');
disp('    -----------------------------');
disp(['   Initial population maximum fitness = ', num2str(maxfitness)]);
disp(['   Initial population average fitness = ', num2str(avgfitness)]);
disp(['   Initial population minimum fitness = ', num2str(minfitness)]);
disp(' ');
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function report(gen, maxfitness, avgfitness, minfitness)
% output the GA progress to the screen
disp(['Generation ', num2str(gen)]);
disp(['   Maximum fitness = ', num2str(maxfitness)]);
disp(['   Average fitness = ', num2str(avgfitness)]);
disp(['   Minimum fitness = ', num2str(minfitness)]);
disp(' ');
return