ORCS(FitnessFunction, Dimensions, LowerBound, UpperBound, Tolerance, EvaluationsMaximum)

function [BestNest, BestFitness, EvaluationsCount] = ORCS(FitnessFunction, Dimensions, LowerBound, UpperBound, Tolerance, EvaluationsMaximum)
% One Rank Cuckoo Search (ORCS) algorithm.
% A simplified implemention of the ORCS optimization (minimization) algorithm.
% Example: [BestNest, BestFitness, EvaluationsCount] = ORCS(@Sphere, 50, -100, 100, 0, 25000)
% Written by Ahmed Twafik.
% Last Updated, 30 April 2013.
% For more information, 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

% 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 evolutions

% 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
nests = LowerBound + rand(PopulationSize,Dimensions) * (UpperBound - LowerBound); % Solutions

% 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');

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

firstNest = nests(1,:);
for n=2:size(nests,1)
	for d=1:size(nests,2)
		if (nests(n,d) != firstNest(d))
			result = false;
			return;
		end
	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 || nests(n,d) > UpperBound)
			if (rand <= bbbFactor)
				boundNests(n,d) = nests(randi(n),d); % Bound by Best nests
			else
				boundNests(n,d) = LowerBound + rand *(UpperBound - LowerBound); % Bound by random
			end
		end
	end
end