double & GlobalMin,

NumericVector & GlobalParams,

NumericVector & f,

NumericMatrix & Foods,

int & unchanged

) {

double GlobalMinOld=GlobalMin;

for (int i=0;i<f.size(); i++)

if (f[i] < GlobalMin) {

GlobalMin = f[i];

GlobalParams = Foods(i,_);

}

if (GlobalMin == GlobalMinOld) ++unchanged;

else unchanged = 0;

return;

int index,

NumericVector & fitness,

NumericVector & f,

IntegerVector & trial,

Function & fun,

NumericMatrix & Foods,

const NumericVector & lb,

const NumericVector & ub

) {

for (int j=0; j<lb.size(); j++)

Foods.at(index,j) = runif(1, lb.at(j), ub.at(j))[0];

f[index] = as<double>(fun(Foods(index,_)));

fitness[index] = CalculateFitness(f[index]);

trial[index] = 0;

return;

double & GlobalMin,

NumericVector & GlobalParams,

NumericVector & fitness,

NumericVector & f,

IntegerVector & trial,

Function & fun,

NumericMatrix & Foods,

const NumericVector & lb,

const NumericVector & ub

) {

int param2change, neighbour;

double ObjValSol, FitnessSol;

NumericVector solution(Foods.ncol());

for (int i=0;i<Foods.nrow();i++) {

// Random parameter to change

param2change = (int)(unif_rand()*Foods.ncol());

// Random neighbour to select

neighbour = i;

while (neighbour == i)

neighbour = (int)(unif_rand()*Foods.nrow());

// Suggesting new solution

solution = Foods(i,_);

solution[param2change] = Foods(i,param2change) +

(Foods(i, param2change) - Foods(neighbour, param2change))*(unif_rand()-.5)*2;

// Truncating

if (solution[param2change] < lb.at(param2change))

solution[param2change] = lb.at(param2change);

if (solution[param2change] > ub.at(param2change))

solution[param2change] = ub.at(param2change);

// Comparing current solution with new one

ObjValSol = as<double>(fun(solution));

FitnessSol = CalculateFitness(ObjValSol);

if (FitnessSol > fitness[i]) {

Foods(i,_) = solution;

fitness[i] = FitnessSol;

f[i] = ObjValSol;

trial[i] = 0;

} else {

trial[i]+=1;

}

}

return;

NumericMatrix & Foods,

NumericVector & fitness,

NumericVector & prob

) {

double maxfit = fitness[0];

for (int i=0;i<Foods.nrow();i++)

if (fitness[i] > maxfit) maxfit = fitness[i];

for (int i=0;i<Foods.nrow();i++)

prob[i] = (0.9*((fitness[i] + 1e-40)/(maxfit + 1e-40))) + 0.1;

return;

double & GlobalMin,

NumericVector & GlobalParams,

NumericVector & fitness,

NumericVector & f,

IntegerVector & trial,

NumericVector & prob,

Function & fun,

NumericMatrix & Foods,

const NumericVector & lb,

const NumericVector & ub) {

int param2change, neighbour;

double ObjValSol, FitnessSol;

NumericVector solution(Foods.ncol());

int t = 0, i=0;

while (t < Foods.nrow()) {

// Randomly choose a food source

if (unif_rand() < prob[i]) {

t++;

// Random parameter to change

param2change = (int)(unif_rand()*Foods.ncol());

// Random neighbour to select

neighbour = i;

while (neighbour == i)

neighbour = (int)(unif_rand()*Foods.nrow());

// Suggesting new solution

solution = Foods(i,_);

solution[param2change] = Foods(i,param2change) +

(Foods(i, param2change) - Foods(neighbour, param2change))*(unif_rand()-.5)*2;

// Truncating

if (solution[param2change] < lb.at(param2change))

solution[param2change] = lb.at(param2change);

if (solution[param2change] > ub.at(param2change))

solution[param2change] = ub.at(param2change);

// Comparing current solution with new one

ObjValSol = as<double>(fun(solution));

FitnessSol = CalculateFitness(ObjValSol);

if (FitnessSol > fitness[i]) {

Foods(i,_) = solution;

fitness[i] = FitnessSol;

f[i] = ObjValSol;

trial[i] = 0;

} else {

trial[i]+=1;

}

} else { /* if */

i++;

if (i==Foods.nrow()) i=0;

}

} /* while */

return;

NumericVector & fitness,

NumericVector & f,

IntegerVector & trial,

NumericVector & prob,

Function & fun,

NumericMatrix & Foods,

const NumericVector & lb,

const NumericVector & ub,

int limit

) {

int maxtrialindex = 0;

for(int i=0;i<Foods.nrow();i++) {

if (trial[i] > trial[maxtrialindex])

maxtrialindex = i;

}

// If it has reach the max, then init again

if (trial[maxtrialindex]>=limit)

init(maxtrialindex, fitness, f, trial, fun, Foods, lb, ub);

return;

NumericVector & par,

Function & fn,

const NumericVector & lb,

const NumericVector & ub,

int FoodNumber = 20,

int limit = 100,

int maxCycle = 1000,

int criter = 50 // , double tol=1e-10

) {

// Initialize:

NumericMatrix Foods(FoodNumber, par.size());

NumericVector f(FoodNumber), fitness(FoodNumber), prob(FoodNumber);

IntegerVector trial(FoodNumber);

NumericMatrix ans(maxCycle, par.size());

// Initializing

NumericVector GlobalParams = clone(par);

double GlobalMin = as<double>(fn(GlobalParams));

// Should be distributed equally

for (int i=0;i<FoodNumber;i++) {

for (int j=0;j<par.size();j++)

Foods.at(i,j) = lb.at(j) + (ub.at(j) - lb.at(j))/(FoodNumber-1.0)*i;

// Checking if it is defined

f[i] = as<double>(fn(Foods(i,_)));

if (NumericVector::is_na(f[i]))

stop("Undefined value for -fn-. Check the function's domain.");

fitness[i] = CalculateFitness(f[i]);

trial[i] = 0;

}

int unchanged = 0;

MemorizeBetsSource(GlobalMin, GlobalParams, f, Foods, unchanged);

ans(0,_) = GlobalParams;

unchanged = 0;

// double OldGlobalMin = GlobalMin;

int i=0;

while (++i<maxCycle) {

SendEmployedBees(GlobalMin, GlobalParams, fitness, f, trial, fn, Foods, lb, ub);

CalculateProbabilities(Foods, fitness, prob);

SendOnlookerBees(GlobalMin, GlobalParams, fitness, f, trial, prob, fn, Foods, lb, ub);

MemorizeBetsSource(GlobalMin, GlobalParams, f, Foods, unchanged);

// Storing and breaking

ans(i,_) = GlobalParams;

if (unchanged>=criter) break;

// If not, then send them out again

SendScoutBees(fitness, f, trial, prob, fn, Foods, lb, ub, limit);

}

return List::create(

_["Foods"] = Foods,

_["f"] = f,

_["fitness"]= fitness,

_["trial"] = trial,

_["value"] = GlobalMin,

_["par"] = GlobalParams,

_["counts"] = i,

_["hist"] = ans(Range(0,i),_)

);