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
% POP = an n by np matrix representing the entire population
% 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 sum(Val)/np]; % collect data for analysis
  iter = iter+1;
end                                 % end of loop over generations

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

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

