...

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.
%
% For more information, please check "Ahmed S. Tawfik, Amr A. Badr and Ibrahim F. Abdel-Rahman. 
% 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