function BBOAckley

% Biogeography-based optimization (BBO) software for minimizing a general function
 
Maxgen = 30; % generation count limit
OPTIONS.popsize = 50; % total population size
OPTIONS.numVar = 20; % number of SIVs in each island
OPTIONS.Resolution = 0.1172; % resolution of each SIV (this gives 9 bits for the Ackley function)

RandSeed = round(sum(100*clock)); % initialize random number generator
[MaxParValue, MinParValue, Population] = Init(OPTIONS); % initialize the population
Population = CostFunction(OPTIONS, Population, MinParValue, MaxParValue);
% Sort the population from most fit to least fit
Population = PopSort(Population);
% Compute the average cost
AverageCost = ComputeAveCost(Population);
% Initialize arrays, display info to screen
MinCost = [Population(1).cost];
AvgCost = [AverageCost];
disp(['The best and mean of Generation # 0 are ', num2str(MinCost(end)), ' and ', num2str(AvgCost(end))]);

Imax = 1; % max immigration rate for each island
Emax = 1; % max emigration rate, for each island
PopSize = OPTIONS.popsize; % max species count, for each island

% Begin the optimization loop
for GenIndex = 1 : Maxgen
    % Map cost values to species counts.
    [Population] = GetSpeciesCounts(Population, PopSize);
    % Compute immigration rate and emigration rate for each species count.
    % lambda(i) is the immigration rate for habitat i.
    % mu(i) is the emigration rate for habitat i.
    [lambda, mu] = GetLambdaMu(Population, Imax, Emax, PopSize);
    % Now use lambda and mu to decide how much information to share between habitats.
    for k = 1 : PopSize
        % Probabilistically input new information into habitat k
        for j = 1 : OPTIONS.numVar
            if rand < lambda(k)
                % Pick a habitat from which to obtain an SIV
                RandomNum = rand * sum(mu);
                Select = mu(1);
                SelectIndex = 1;
                while (RandomNum > Select) & (SelectIndex < OPTIONS.popsize)
                    SelectIndex = SelectIndex + 1;
                    Select = Select + mu(SelectIndex);
                end
                Island(k,j) = Population(SelectIndex).chrom(j);
            else
                Island(k,j) = Population(k).chrom(j);
            end
        end
    end
    % Replace the habitats with their new versions.
    for k = 1 : PopSize
        Population(k).chrom = Island(k,:);
    end
    % Make sure each individual is legal.
    Population = Bound(OPTIONS, Population, MinParValue, MaxParValue);
    % Calculate cost
    Population = CostFunction(OPTIONS, Population, MinParValue, MaxParValue);
    % Sort from best to worst
    Population = PopSort(Population);
    % Sort from best to worst
    Population = PopSort(Population);
    % Compute the average cost
    [AverageCost, nLegal] = ComputeAveCost(Population);
    % Display info to screen
    MinCost = [MinCost Population(1).cost];
    AvgCost = [AvgCost AverageCost];
    disp(['The best and mean of Generation # ', num2str(GenIndex), ' are ',...
            num2str(MinCost(end)), ' and ', num2str(AvgCost(end))]);
end
Conclude(Population, nLegal, MinCost, AvgCost);
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Population] = GetSpeciesCounts(Population, P)

% Map cost values to species counts.
% This loop assumes the population is already sorted from most fit to least fit.
for i = 1 : length(Population)
    if Population(i).cost < inf
        Population(i).SpeciesCount = P - i;
    else
        Population(i).SpeciesCount = 0;
    end
end
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [lambda, mu] = GetLambdaMu(Population, I, E, P)
% Compute immigration rate and extinction rate for each species count.
% lambda(i) is the immigration rate for individual i.
% mu(i) is the extinction rate for individual i.
for i = 1 : length(Population)
    lambda(i) = I * (1 - Population(i).SpeciesCount / P);
    mu(i) = E * Population(i).SpeciesCount / P;
end
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [MaxParValue, MinParValue, Population] = Init(OPTIONS)
% Initialize population
MinParValue = 1;
MaxParValue = floor(1 + 2 * 30 / OPTIONS.Resolution);
for popindex = 1 : OPTIONS.popsize
    chrom = floor(MinParValue + (MaxParValue - MinParValue + 1) * rand(1,OPTIONS.numVar));
    Population(popindex).chrom = chrom;
end
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Population] = Bound(OPTIONS, Population, MinParValue, MaxParValue)
% Make sure each population member is within constraints
for i = 1 : OPTIONS.popsize
    for k = 1 : OPTIONS.numVar
        Population(i).chrom(k) = max(Population(i).chrom(k), MinParValue);
        Population(i).chrom(k) = min(Population(i).chrom(k), MaxParValue);
    end
end
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Population] = CostFunction(OPTIONS, Population, MinParValue, MaxParValue)
% Compute the cost of each member in Population
popsize = OPTIONS.popsize;
p = OPTIONS.numVar;
for popindex = 1 : popsize
    sum1 = 0;
    sum2 = 0;
    for i = 1 : p
        gene = Population(popindex).chrom(i);
        x = (gene - MinParValue) / (MaxParValue - MinParValue) * 2 * 30 - 30;
        sum1 = sum1 + x^2;
        sum2 = sum2 + cos(2*pi*x);
    end
    Population(popindex).cost = 20 + exp(1) - 20 * exp(-0.2*sqrt(sum1/p)) - exp(sum2/p);    
end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Population, indices] = PopSort(Population)
% Sort the population members from best to worst
popsize = length(Population);
Cost = zeros(1, popsize);
indices = zeros(1, popsize);
for i = 1 : popsize
    Cost(i) = Population(i).cost;
end
[Cost, indices] = sort(Cost, 2, 'ascend');
Chroms = zeros(popsize, length(Population(1).chrom));
for i = 1 : popsize
    Chroms(i, :) = Population(indices(i)).chrom;
end
for i = 1 : popsize
    Population(i).chrom = Chroms(i, :);
    Population(i).cost = Cost(i);
end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [AveCost, nLegal] = ComputeAveCost(Population)
% Compute the average cost of all legal individuals in the population.
% OUTPUTS: AveCost = average cost
%          nLegal = number of legal individuals in population
% Save valid population member fitnesses in temporary array
Cost = [];
nLegal = 0;
for i = 1 : length(Population)
    if Population(i).cost < inf
        Cost = [Cost Population(i).cost];
        nLegal = nLegal + 1;
    end
end
% Compute average cost.
AveCost = mean(Cost);
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Conclude(Population, nLegal, MinCost, AvgCost)
% Output results of population-based optimization algorithm.
% Count the number of duplicates
Skip = [];
NumDups = 0;
for i = 1 : length(Population)-1
    if ~isempty(find(Skip, i)), continue, end
    Chrom1 = sort(Population(i).chrom);
    for j = i+1 : length(Population)
        Chrom2 = sort(Population(j).chrom);
        if isequal(Chrom1, Chrom2)
            if isempty(find(Skip, j))
                Skip = [Skip j];
            end
            if (NumDups == 0)
                NumDups = 2;
            else
                NumDups = NumDups + 1;
            end
        end
    end
end  
disp([num2str(NumDups), ' duplicates in final population.']);
disp([num2str(nLegal), ' legal individuals in final population.']);
% Display the best solution
Chrom = sort(Population(1).chrom);
disp(['Best chromosome = ', num2str(Chrom)]); 
% Plot some results
close all;
plot([0:length(MinCost)-1], MinCost, 'r');
hold on;
plot([0:length(AvgCost)-1], AvgCost, 'b:');
legend('minimum cost', 'average cost');
xlabel('Generation');
ylabel('Cost');
return;