function [Bestmem,bestval,nfe,DATA] = DE33F(fcn,Xmin,Xmax,f,c,np,tol,itMax)
% Function implementing multi-dimensional
% minimization using differential evolution of R. Storn
% R. Storn, "Differential Evolution - A Simple and Efficient Heuristic
% Strategy for Global Optimization over Continuous Spaces",
% J. Global Optimization, vol. 11, p.341-359 (1997).
% The present Matlab program is based on the one by Rainer Storn
% (http://www.icsi.berkeley.edu/~storn/code.htm), which was
% improved by Ken Price, Arnold Neumaier and Jim Van Zandt
% The present program:
% (1) constructs the trial vectors more efficiently
% (2) makes the program loopless except for the main iteration loop:
% key point: operate on the entire population in "parallel"
% (3) requires only one copy of the population
% (4) uses identifier names that reflect the matrix nature of the variables
% (5) implements only scheme: DE/rand/1/bin for constructing trial vectors
% INPUT:
% fcn = string name of the cost function
% Xmin = lower limits of initial parameters
% Xmax = upper limits of initial parameters
% f = amplification factor
% c = crossover probability
% np = number of individuals in the population
% Tol = tolerance in solution norm or function value
% itMax = maximum number of iteration allowed
% OUTPUT:
% Bestmem = best ever individual
% bestval = best ever cost function value
% iter = actual number of iterations
% DATA = history of the best function value in each generation

% K. Ming Leung, 10/06/03, version 1
% 10/23/03, version 2

n = length(Xmin);                          % number of parameters
Ones = ones(1,np);
% Create initial population as an n by np matrix
POP = Xmin*Ones + rand(n,np).*((Xmax-Xmin)*Ones);
Val = feval(fcn,POP),                   % function values of population (row vector)
[bestval,ibest] = min(Val,[],2);      % find the best value
Bestmem = POP(:,ibest);              % genes of the best ever individual

% Initialization of vector and matrices:
% Matrices of 3 randomly chosen vectors to construct trial vectors
P1 = zeros(n,np);
P2 = P1;
P3 = P2;
TRIAL = zeros(n,np);                % Matrix of trial vectors for population
MASK = zeros(n,np);                % mask for selecting trial vectors

% Indices used to construct trial vectors
I1 = zeros(np);     I2 = zeros(np);     I3 = zeros(np);
Shift = [2:np 1];
% shifting indices 1 to the right
del = 1;                                     % absolute changes in function or vector norm
iter = 1;                                     % iteration counter
DATA = [];                               % history of the best function values

% Loop through the generations
while ((iter<itMax) & (del>tol))
I1 = randperm(np);                    % index to randomize vectors in population
I2 = I1(Shift);                            % index for vectors in population 2
I3 = I2(Shift);                            % index for vectors in population 3
P1 = POP(:,I1);                         % shuffled population 1
P2 = POP(:,I2);                         % shuffled population 2
P3 = POP(:,I3);                         % shuffled population 3

% randomly choose elements of trial vectors to be constructed
MASK = (rand(n,np) < c);
TRIAL = POP;                         % copy population as trial vectors
% construct trial vectors according to scheme DE/rand/1/bin
TRIAL(MASK) = P1(MASK) + f*(P2(MASK)-P3(MASK));
Tempval = feval(fcn,TRIAL);     % function values of trial vectors
Idx = (Tempval<Val);                % indices of the better ones
POP(:,Idx) = TRIAL(:,Idx);       % replace each by its better off-spring
Val(Idx) = Tempval(Idx);           % update function values of population
[bestval1,ibest] = min(Val,[],2); % find best individual of the generation
Bestmem1 = POP(:,ibest);         % find its function value
% See if new generation has a better individual than the previous generation
if bestval1 < bestval
%del = norm(Bestmem - Bestmem1); % control error of parameter vector
del = bestval - bestval1;             % control error of function value
bestval = bestval1;                     % update the best function value
Bestmem = Bestmem1;              % update parameter of best individual
end;
DATA = [DATA; iter bestval];  % collect data for analysis
iter = iter+1;
end                                            % end of loop over generations

if iter >= itMax
error('Maximum number of iteration reached without convergence!');
end

nfe = (iter-1)*np;                        % number of function evaluations