Code covered by the BSD License

# One Rank Cuckoo Search (ORCS) algorithm

### Ahmed Tawfik (view profile)

01 May 2013 (Updated )

An improved cuckoo search optimization algorithm. Well competent and easy to use.

...
function [BestNest, BestFitness, EvaluationsCount] = ...
ORCS(FitnessFunction, Dimensions, LowerBound, UpperBound, Tolerance, EvaluationsMaximum)

% One Rank Cuckoo Search (ORCS) algorithm
% Implementation of the ORCS optimization algorithm.
%
% Usage: [BestNest, BestFitness, EvaluationsCount] = ...
% 	ORCS(FitnessFunction, Dimensions, LowerBound, UpperBound, Tolerance, EvaluationsMaximum)
%
% Example: [BestNest, BestFitness, EvaluationsCount] = ORCS(@Sphere, 5, -100, 100, 0, 25000)
% Example: [BestNest, BestFitness, EvaluationsCount] = ORCS(@Sphere, 5, -100 * ones(1,5), 100 * ones(1,5), 0, 25000)
%
% Written by Ahmed Twafik.
% Last Updated, 3 May 2013.
%
% One Rank Cuckoo Search Algorithm with Application to Algorithmic Trading Systems Optimization.
% International Journal of Computer Applications 64(6):30-37, February 2013."
%

if nargin < 6
help ORCS
return;
end

% Normalize inputs
LowerBound = NormalizeBounds(LowerBound, Dimensions);
if isempty(LowerBound)
help ORCS
return;
end

UpperBound = NormalizeBounds(UpperBound, Dimensions);
if isempty(UpperBound)
help ORCS
return;
end

% ORCS defaults (change, if required)
PopulationSize = 15;
Pa = 0.25; % Probability of abandon (discovery rate)
OrT = 0.1; % One Rank ratio update trigger, as a percent of maximum no. of evaluations

% Initialize
Beta = 1.5;
sigma = ((gamma(1 + Beta) * sin(pi * Beta / 2) / (gamma((1 + Beta) / 2) * Beta * 2^((Beta - 1) / 2))))^(1.0 / Beta);

orRatio = 1; % One Rank ratio
orDelta = 1.0 - 0.5 / Dimensions; % One Rank ratio-update delta

% Initialize nests (solutions)
nests = repmat(LowerBound, PopulationSize, 1) + ...
rand(PopulationSize, Dimensions) .* repmat(UpperBound - LowerBound, PopulationSize, 1);

% Evaluate
progress.Fitness = [];
progress.Evaluations = [];
[nests, fitness, evaluationsCount, progress] = EvaluateNests(nests, LowerBound, UpperBound, FitnessFunction, 0, progress);

% Get best
[bestFitness,bestNestIndex] = min(fitness);
bestNest = nests(bestNestIndex,:);

% Iterations
while (evaluationsCount < EvaluationsMaximum && bestFitness > Tolerance && !EqualNests(nests))
newNests = GetNewNests(nests, bestNest, Beta, sigma); % Get new nests (Levy flights)
if (orRatio <= rand)
[nests, fitness, bestNest, bestFitness, evaluationsCount, progress, orRatio] = ...
GetBestNest(newNests, LowerBound, UpperBound, FitnessFunction, evaluationsCount, progress, nests, fitness ...
, bestNest, bestFitness, OrT, EvaluationsMaximum, orRatio, orDelta); % Get best
newNests = EmptyNests(nests, Pa); % Empty nests (discovery rate)
else
newNests = EmptyNests(newNests, Pa); % Empty new nests (discovery rate)
end
[nests, fitness, bestNest, bestFitness, evaluationsCount, progress, orRatio] = ...
GetBestNest(newNests, LowerBound, UpperBound, FitnessFunction, evaluationsCount, progress, nests, fitness ...
, bestNest, bestFitness, OrT, EvaluationsMaximum, orRatio, orDelta); % Get best
end

% Display
BestNest = bestNest;
BestFitness = bestFitness;
EvaluationsCount = evaluationsCount;

fprintf('Best fitness is %g\r\n', BestFitness);

% Draw
figure;
if(progress.Fitness(end) > 0)
semilogy(progress.Evaluations, progress.Fitness, '-r', 'LineWidth',2);
else
plot(progress.Evaluations, progress.Fitness, '-r', 'LineWidth',2);
end
grid on;
title('ORCS (Fitness / Evaluations)');
xlabel('Evaluations');
ylabel('Fitness');

% NormalizeBounds -------------------------------------------------------------
function bounds = NormalizeBounds(bounds, Dimensions)

if isequal(size(bounds), [1,1])
bounds = bounds * ones(1, Dimensions);
end

if !isequal(size(bounds), [1,Dimensions])
bounds = [];
end

% EqualNests ------------------------------------------------------------------
function result = EqualNests(nests)

firstNest = nests(1,:);
for n=2:size(nests,1)
if(!isequal(nests(n), firstNest))
result = false;
return;
end
end

result = true;

% GetNewNests -----------------------------------------------------------------
function newNests = GetNewNests(nests, bestNest, Beta, sigma)

for n=1:size(nests,1)
for d=1:size(nests,2)
scale = 0.01 * randn();
lfs = randn() * sigma / abs(randn())^(1.0 / Beta); % Levy-flight sample
step = scale * lfs * (nests(n,d) - bestNest(d));
newNests(n,d) = nests(n,d) + step;
end
end

% EmptyNests ------------------------------------------------------------------
function newNests = EmptyNests(nests, Pa)

newNests = nests;
PopulationSize = size(nests, 1);

% Get 2 random-permutations of population size
rpx = randperm(PopulationSize);
rpy = randperm(PopulationSize);

for n=1:PopulationSize
for d=1:size(nests, 2)
if (rand <= Pa) % Discovered
step = rand * (nests(rpx(n),d) - nests(rpy(n),d));
newNests(n,d) = nests(n,d) + step;
end
end
end

% GetBestNest -----------------------------------------------------------------
function [nests, fitness, bestNest, bestFitness, evaluationsCount, progress, orRatio] = ...
GetBestNest(newNests, LowerBound, UpperBound, FitnessFunction, evaluationsCount, progress, nests, fitness ...
, bestNest, bestFitness, OrT, EvaluationsMaximum, orRatio, orDelta)

% Evaluate nests
[newNests, newFitness, evaluationsCount, progress] = EvaluateNests(newNests, LowerBound, UpperBound, FitnessFunction ...
, evaluationsCount, progress);

% Get best nests
for n=1:size(nests, 1)
if (newFitness(n) <= fitness(n))
nests(n,:) = newNests(n,:);
fitness(n) = newFitness(n);
end
end

% Update bestNest/orRatio
[newBestFitness,bestNestIndex] = min(fitness);
if newBestFitness < bestFitness
bestNest = nests(bestNestIndex,:);
bestFitness = newBestFitness;
elseif evaluationsCount - progress.Evaluations(end) > OrT * EvaluationsMaximum
orRatio = orRatio * orDelta;
end

% EvaluateNests ---------------------------------------------------------------
function [nests, fitness, evaluationsCount, progress] = EvaluateNests(nests, LowerBound, UpperBound, FitnessFunction ...
, evaluationsCount, progress)

PopulationSize = size(nests, 1);

% Bound
nests = Bound(nests, LowerBound, UpperBound);

% Evalute fitness function
fitness = inf * ones(PopulationSize, 1);
for n=1:PopulationSize
fitness(n) = FitnessFunction(nests(n,:));
end

% Update evaluations count
evaluationsCount = evaluationsCount + PopulationSize;

% Update progress (if fitness improved)
if size(progress.Fitness, 1) > 0
lastFitness = progress.Fitness(end);
else
lastFitness = inf;
end

bestFitness = min(fitness);
if (bestFitness < lastFitness)
progress.Fitness = [progress.Fitness bestFitness];
progress.Evaluations = [progress.Evaluations evaluationsCount];
end

% Bound -----------------------------------------------------------------------
function boundNests = Bound(nests, LowerBound, UpperBound)

boundNests = nests;

PopulationSize = size(nests, 1);
Dimensions = size(nests, 2);

bbbFactor = 1.0 - 1.0 / sqrt(Dimensions); % Bound by Best factor

for n=1:PopulationSize
for d=1:Dimensions
if (nests(n,d) < LowerBound(d) || nests(n,d) > UpperBound(d))
if (rand <= bbbFactor)
boundNests(n,d) = nests(randi(n),d); % Bound by Best nests
else
boundNests(n,d) = LowerBound(d) + rand *(UpperBound(d) - LowerBound(d)); % Bound by a random value
end
end
end
end