Code covered by the BSD License

### Highlights fromDifferential Evolution

from Differential Evolution by Markus Buehren
Optimization using the evolutionary algorithm of Differential Evolution.

differentialevolution(...
```function varargout = differentialevolution(...
DEParams, paramDefCell, objFctHandle, objFctSettings, ...
objFctParams, emailParams, optimInfo)
%DIFFERENTIALEVOLUTION  Start Differential Evolution optimization.
%    BESTMEM = DIFFERENTIALEVOLUTION(DEPARAMS, ...) starts a Differential
%    Evolution (DE) optimization to minimize the cost returned by a given
%    function. For a quick start, check out and modify the functions DEMO1
%    and DEMO2.
%
%    Output arguments:
%
%    bestmem          - Best population member.
%    bestval          - Lowest evaluated cost.
%    bestFctParams    - Structure like input objFctParams containing the
%                       best parameter set.
%    nrOfIterations   - Number of iterations done.
%
%
%    Input arguments:
%
%    DEParams         - Struct with parameters for DE.
%    paramDefCell     - Cell specifying the parameters to optimize.
%    objFctHandle     - Handle to the objective function.
%    objFctSettings   - Additional settings to be passed (a cell array will
%                       be expanded using {:}). If no additional settings
%                       are needed, set objFctSettings to an empty cell: {}
%    objFctParams     - Struct with initial parameters.
%                       Parameters are used for sending E-mail
%                       and 'subtitle' are displayed and included in saved
%                       files if existing. No influence on optimization.
%
%    The structure DEParams needs to contain the following fields:
%
%    VTR            "Value To Reach" (set to empty matrix for no VTR).
%    NP             Number of population members (e.g. 10*D).
%    F              DE-stepsize F from interval [0, 2]. A good initial guess
%                   is the interval [0.5, 1], e.g. 0.8.
%    CR             Crossover probability constant from interval [0, 1]. If
%                   the parameters are correlated, high values of CR work
%                   better. The reverse is true for no correlation.
%    strategy       1 --> DE/best/1/exp (def.)   6 --> DE/best/1/bin
%                   2 --> DE/rand/1/exp          7 --> DE/rand/1/bin
%                   3 --> DE/rand-to-best/1/exp  8 --> DE/rand-to-best/1/bin
%                   4 --> DE/best/2/exp          9 --> DE/best/2/bin
%                   5 --> DE/rand/2/exp          else  DE/rand/2/bin
%                   Experiments suggest that /bin likes to have a slightly
%                   larger CR than /exp
%    maxiter        Maximum number of iterations.
%    maxtime        Maximum time (in seconds) before finishing optimization.
%                   Set to empty or Inf for no time limit.
%    maxclock       Time (as returned by function clock.m) when to
%                   finish optimization. Set to empty for no end time.
%    minvalstddev   Population is reinitialized if the standard deviation of
%                   the cost values in the population is lower than
%                   minvalstddev.
%    minparamstddev Population is reinitialized if the maximum parameter
%                   standard deviation (normalized to the parameter range)
%                   is lower than minparamstddev.
%    nofevaliter    Population is reinitialized if there was no function
%                   evaluation during the last nofevaliter iterations.
%    nochangeiter   Population is reinitialized if there was no change in
%                   the population during the last nochangeiter iterations.
%    refreshiter    Info is displayed and current state is saved every
%                   refreshiter iterations.
%    refreshtime    State is saved after refreshtime seconds.
%    refreshtime2   Additional progress information is displayed every
%                   refreshtime2 seconds (usually refreshtime2 >>
%                   refreshtime).
%    refreshtime3   Progress information is sent by E-mail every
%                   refreshtime3 seconds (usually refreshtime3 >>
%                   refreshtime2).
%    useInitParams  If one, the given parameters in struct objFctParams
%                   OR those in the fourth column of paramDefCell are used
%                   as one of the initial population members. If two,
%                   additionally the first twenty percent of the population
%                   members are computed as the given initial parameter
%                   vector plus small random noise.
%    saveHistory    Save intermediate results.
%    displayResults Draw graphs for visualization of the optimization
%                   result.
%    feedSlaveProc  Use slave process for parallel computation.
%    maxMasterEvals Maximum number of function evaluations done by the
%                   master process. Warning: Use this option with caution!
%                   If maxMasterEvals is set to a number less than the
%                   number of population members and one of the slave
%                   processes is interrupted, the optimization will be
%                   stuck!
%    slaveFileDir   Base directory for saving slave files.
%    minimizeValue  If true, the evaluation value is minimized, otherwise
%                   maximized.
%    validChkHandle Handle to a function which takes the same arguments as
%                   the objective function and decides if a given parameter
%                   set is valid (subject to a constraint) or not. The
%                   function should return 1 if the constraint is fulfilled
%                   and 0 if not. Set to empty string if no constraint is
%                   needed.
%    playSound      Play a short sound when a new best member was found.
%
%    If DEParams is empty or fields are missing, default parameters are used
%    (see function GETDEFAULTPARAMS).
%
%    The cell array paramDefCell has to contain the names of the parameters
%    to optimize, its ranges, their quantizations and the initial values.
%    Each parameter may be a real-valued scalar or column vector.
%
%    Example 1 (only scalar parameters):
%
%    paramDefCell = {
%      'useSmoothing',     [0    1],     1,   0
%      'nrOfCoefficients', [5   20],     1,  10
%      'threshold',        [0.01 1], 0.001, 0.5 }
%
%    The first cell in each row contains the name of the parameter, the
%    second a two-element row vector specifying the allowed range, the third
%    the quantization and the fourth the initial values (the fourth column
%    of the cell array may be omitted). Provide a non-empty value either in
%    objFctParams or in the fourth column of paramDefCell as initial value.
%    If both are present, a warning message is issued and the value in
%    paramDefCell is used. If objFctParams is empty and no initial
%    parameters are given in paramDefCell, the centers of the parameter
%    ranges are used as initial parameters.
%
%    Using parameter quantization allows for the use of binary variables
%    like 'useSmoothing' above as well for parameters that are of integer
%    nature, like a number of coefficients. If the quantization of a
%    parameter is set to zero, the parameter is not quantized. Using a
%    quantization grid for continuous parameters can save computational
%    effort. If saveHistory is true, all evaluated parameter vectors are
%    saved with the corresponding cost value and the same parameter value
%    will never be evaluated twice. With quantization, it is more likely
%    that a generated parameter vector was already evaluated and saved
%    before.
%
%    Example 2 (vector parameter):
%
%    paramDefCell = {'weightings', [0 1; 0 2], [0.01; 0.02], [0.5; 0.5]};
%
%    Here, the parameter weightings is defined as a two-element column
%    vector. The ranges are set to [0, 1] for the first element and [0, 2]
%    for the second. The quantizations are 0.01 and 0.02 and the initial
%    values are both 0.5.
%
%    The objective function (given as function handle objFctHandle) is
%    started as
%
%    value = objFctHandle(objFctSettings,    objFctParams) or
%    value = objFctHandle(objFctSettings{:}, objFctParams).
%
%    The second case is used if objFctSettings is a cell array, thus
%    allowing for an arbitrary number of additional input arguments. The
%    provided structure objFctParams may contain further fixed parameters
%    and/or the current parameter values. The fields with the names of the
%    parameters given in paramDefCell are overwritten by the values of the
%    current parameters. If the objective function handle is empty, the
%    distance to a randomly chosen optimal parameter vector is used as cost
%    value (for testing purposes).
%
%    Example 3 (vector parameter):
%
%    paramDefCell = {'', [0 1; 0 2], [0.01; 0.02], [0.5; 0.5]};
%
%    In this special case (one single parameter with empty name), the
%    objective function is called as
%
%    value = objFctHandle(objFctSettings,    paramVec) or
%    value = objFctHandle(objFctSettings{:}, paramVec)
%
%    with the current parameters in column vector paramVec.
%
%    When displaying an info string, the current optimization state
%    including all tested members etc. is saved in the file
%    XXX_result_YYY_ZZ.mat, where XXX is the name of the objective function,
%    YYY is the name of the current host and ZZ is a number between 1 and 50
%    (to avoid overwriting old results).
%
%    A 'time over'-file is saved at the start of the optimization. The
%    optimization is stopped if this file is deleted. Using this mechanism
%    to stop the simulation avoids to break Matlab during saving a file,
%    which can make a file unaccessible for the rest of the session and
%    leads to repeating warning messages. The name of the file to delete is
%    XXX_timeover_YYY.mat, where XXX is the name of the objective function
%    and YYY is the hostname. Result- and 'time over'-files are saved in
%    directory 'data' if existing, otherwise in the current directory.
%
%    The optimization can be performed in parallel by more than one
%    processor/computer. Function DIFFERENTIALEVOLUTION has to be started in
%    one Matlab session, function DIFFERENTIALEVOLUTIONSLAVE on one or
%    more other Matlab sessions. Function DIFFERENTIALEVOLUTION acts as
%    master and saves information about which function to evaluate and which
%    parameters to use into files in a commonly accessible directory. The
%    Distributed Computing toolbox is not used. If input parameter
%    slaveFileDir is empty, the directory differentialevolution is used (or
%    created) below the temporary directory returned by function TEMPDIR2
%    (something like C:\Documents and Settings\<UserName>\Local Settings\
%
%    Function DIFFERENTIALEVOLUTION was developed for objective functions
%    that need relatively long for one function evaluation (several seconds
%    or more). When used with objective functions that evaluate very fast,
%    memory problems could occur. When saveHistory is true, every evaluated
%    parameter vector is kept in memory. Further, the overhead for checking
%    if a parameter vector was already evaluated might be larger than a
%    function evaluation itself.
%
%    Start this function without input arguments or with only the first
%    input argument DEParams to run a demo optimization of Rosenbrock's
%    saddle. No files are saved during the demo.
%
%    This function is based on the differential evolution (DE) algorithm of
%    Rainer Storn (http://www.icsi.berkeley.edu/~storn/code.html). The core
%    evolutionary algorithm was taken from function devec3.m.
%
%    <a href="differentialevolution.html">differentialevolution.html</a>  <a href="http://www.mathworks.com/matlabcentral/fileexchange/18593">File Exchange</a>  <a href="https://www.paypal.com/cgi-bin/webscr?cmd=_s-xclick&hosted_button_id=KAECWD2H7EJFN">Donate via PayPal</a>
%
%    Markus Buehren
%
%    GETDEFAULTPARAMS, DEMO1, DEMO2, TEMPDIR2.

% get default DE parameters
DEParamsDefault = getdefaultparams;

% set text width for wrapping displayed information
textWidth = 75;

% get DE parameters from input structure
if nargin == 0 || isempty(DEParams)
DEParams = DEParamsDefault;
else
fieldNames = fieldnames(DEParamsDefault);
for k=1:length(fieldNames)
if ~isfield(DEParams, fieldNames{k})
DEParams.(fieldNames{k}) = DEParamsDefault.(fieldNames{k});
disp(textwrap2(sprintf(['Warning: Field ''%s'' not included in DEParams. ', ...
'Using default value.'], fieldNames{k}), textWidth));
end
end
end

switch nargin
case {0, 1}
% generate default parameter set for demonstration
objFctParams.parameter1 = -1;
objFctParams.parameter2 = -1;
objFctSettings          = 100;
paramDefCell = {
'parameter1', [-2 2], 0.05
'parameter2', [-2 2], 0.05};
optimInfo.title         = 'Optimization of Rosenbrock''s saddle';
emailParams             = [];
DEParams.feedSlaveProc  = 1;
DEParams.refreshiter    = 1;
DEParams.refreshtime    = 10;  % in seconds
DEParams.refreshtime2   = 20;  % in seconds
DEParams.refreshtime3   = 40;  % in seconds
DEParams.maxiter        = 100;
DEParams.maxtime        = 60;  % in seconds
rand('state', 1); % always use the same population members

case 2
error(textwrap2('Wrong number of input arguments.'));

otherwise
if ~exist('objFctSettings', 'var'), objFctSettings = {}; end
if ~exist('objFctParams',   'var'), objFctParams   = []; end
if ~exist('emailParams',    'var'), emailParams    = []; end
if ~exist('optimInfo', 'var') || isempty(optimInfo) || ~isstruct(optimInfo)
optimInfo       = [];
optimInfo.title = 'DE optimization';
end
end

% check paramDefCell
checkinputs(paramDefCell, objFctParams);

% modify paramDefCell if there are vector-valued parameters
k = 1;
parameterDimVector = [];
while k <= size(paramDefCell, 1)
parameterDim = size(paramDefCell{k,2}, 1);
if parameterDim == 1
% scalar parameter, save dimension
parameterDimVector(k,1) = 1; %#ok
k = k + 1;
if isempty(paramDefCell{1,1})
% scalar parameter with empty name
paramDefCell{1,1} = '_1';
end
else
% vector parameter, introduce new rows in paramDefCell
parameterDimVector(k:k+parameterDim-1,1) = parameterDim; %#ok
parameterName = paramDefCell{k,1};
paramDefCell = [paramDefCell(1:k,:); ...
cell(parameterDim-1, size(paramDefCell,2)); paramDefCell(k+1:end,:)];
for d = parameterDim:-1:1
paramDefCell{k+d-1, 1} = sprintf('%s_%d', parameterName, d);
for col = 2:size(paramDefCell,2)
paramDefCell{k+d-1,col} = paramDefCell{k, col}(d,:);
if col == 4 && isnan(paramDefCell{k+d-1,col}) % initial value = NaN
paramDefCell{k+d-1,col} = [];
end
end
end
k = k + parameterDimVector(k,1);
end
end

% initialize functions
getparametername (paramDefCell, parameterDimVector);
displaybestmember(paramDefCell);

% get parameter bounds
parameterBounds = cell2mat(paramDefCell(:,2));
parGridVector   = cell2mat(paramDefCell(:,3));
D     = size(parameterBounds, 1);
XVmin = parameterBounds(:, 1)';
XVmax = parameterBounds(:, 2)';

% check values
errorFound = 0;
for parNr = 1:D
if XVmin(parNr) > XVmax(parNr)
fprintf('Error: Lower bound (%g) of parameter %s is larger than upper bound (%g).\n', ...
XVmin(parNr), getparametername(parNr, 1), XVmax(parNr));
errorFound = true;
end
if parGridVector(parNr) < 0
fprintf('Error: Negative quantization values are not allowed (parameter %s).\n', ...
getparametername(parNr, 1));
errorFound = true;
end
minQuantizationUser = 1e-12;
if parGridVector(parNr) > 0 && parGridVector(parNr) < minQuantizationUser
fprintf('Error: Minimum quantization step size is %g (parameter %s).\n', ...
minQuantizationUser, getparametername(parNr, 1));
errorFound = true;
end
end
if errorFound
error('Erroneous parameters found.');
end

% compute number of possible parameter vectors
if all(parGridVector > 0)
nrOfPossibleMembers = prod(floor((diff(parameterBounds, 1, 2) + 0.5*parGridVector) ./ parGridVector) + 1);
else
nrOfPossibleMembers = inf;
end

% check parameters
DEParams.NP = min(DEParams.NP, nrOfPossibleMembers);
if (DEParams.maxiter <= 0)
error('maxiter must be greater than zero.');
end
if DEParams.displayResults && ~DEParams.saveHistory
fprintf('Warning: Optimization history can not be displayed if not saved.\n');
DEParams.displayResults = 0;
end
DEParams.refreshiter = floor(DEParams.refreshiter);

% get parameters
% TODO: pass struct DEParams to subfunctions instead of unwrapping
% everything.
NP                = DEParams.NP;
VTR               = DEParams.VTR;
refreshtime       = DEParams.refreshtime;
refreshtime2      = DEParams.refreshtime2;
refreshtime3      = DEParams.refreshtime3;
maxtime           = DEParams.maxtime;
maxclock          = DEParams.maxclock;
saveHistory       = DEParams.saveHistory;
displayResults    = DEParams.displayResults;
feedSlaveProc     = DEParams.feedSlaveProc;
slaveFileDir      = DEParams.slaveFileDir;
maxMasterEvals    = DEParams.maxMasterEvals;
playSound         = DEParams.playSound;
validChkHandle    = DEParams.validChkHandle;

% check if validChkHandle is valid
if isempty(validChkHandle)
validChkHandle = @alwaysvalid; % handle to function that always returns true
elseif ~isa(validChkHandle, 'function_handle')
if ischar(validChkHandle)
validChkHandle = str2func(validChkHandle);
else
error('validChkHandle is neither empty nor a string nor a function handle.');
end
end

% get initial parameter vector
if DEParams.useInitParams
initialMem = zeros(1,D);
if size(paramDefCell, 2) == 4
% use initial parameters from fourth column of paramDefCell
parNr = 1;
while parNr <= D
if ~isempty(paramDefCell{parNr, 4})
% check if objFctParams also contains initial value
parameterName = getparametername(parNr, 2);
index = parNr:parNr+parameterDimVector(parNr)-1;
initialMem(index) = cell2mat(paramDefCell(index,4))';

% warn if overwriting intial values in objFctParams
if isstruct(objFctParams) && isfield(objFctParams, parameterName) && ...
~isempty(objFctParams.(parameterName)) && ~all(isnan(objFctParams.(parameterName)))
if ~isequal(objFctParams.(parameterName), initialMem(index))
disp(textwrap2(sprintf(['Warning: Initial value of parameter ''%s'' in ', ...
'struct objFctParams was overwritten by settings in paramDefCell.'], ...
parameterName), textWidth));
end
elseif isnumeric(objFctParams) && ~isempty(objFctParams) && ...
~all(isnan(objFctParams))
if isscalar(objFctParams)
disp(textwrap2(sprintf(['Warning: Initial parameter value ', ...
'objFctParams was overwritten by settings in paramDefCell.'], ...
parameterName), textWidth));
else
disp(textwrap2(sprintf(['Warning: Initial parameter vector ', ...
'objFctParams was overwritten by settings in paramDefCell.'], ...
parameterName), textWidth));
end
end
end

for d=1:parameterDimVector(parNr)
index = parNr+d-1;
if isnan(initialMem(index))
parameterName = getparametername(index, 1);
disp(textwrap2(sprintf(['Warning: No initial value for parameter %s given. ', ...
'Using center of range interval as initial value.'], parameterName), textWidth));
initialMem(index) = 0.5*(XVmin(index) + XVmax(index));
end
end
parNr = parNr + parameterDimVector(parNr);
end
checkInitialMem = true;

elseif isempty(objFctParams)
% no initial parameter set given
disp('Warning: Option DEParams.useInitParams is true but no initial parameter set is given.');
DEParams.useInitParams = 0;
checkInitialMem = false;
else
% check if initial values are given by objFctParams
% and collect them
if isstruct(objFctParams)
parNr = 1;
while parNr <= D
parameterDim  = parameterDimVector(parNr);
index = parNr:parNr+parameterDim-1;
parameterName = getparametername(parNr, 2);
if isfield(objFctParams, parameterName)
initialMem(index) = objFctParams.(parameterName)';
else
disp(textwrap2(sprintf(['Warning: No initial value for parameter %s given. ', ...
'Using center of range interval as initial value.'], parameterName), textWidth));
initialMem(index) = 0.5*(XVmin(index) + XVmax(index));
end
parNr = parNr + parameterDim;
end
elseif strcmp(paramDefCell{1,1}, '_1')
initialMem = objFctParams;
for parNr = 1:D
parameterName = getparametername(parNr, 1);
if isnan(initialMem(parNr))
disp(textwrap2(sprintf(['Warning: No initial value for parameter %s given. ', ...
'Using center of range interval as initial value.'], parameterName), textWidth));
initialMem(parNr) = 0.5*(XVmin(parNr) + XVmax(parNr));
end
end
end
checkInitialMem = true;
end

% check if initial member is on quantization grid
if checkInitialMem
[ignore, quantInitialMem] = considerparametercontraints(...
objFctParams, paramDefCell, parameterDimVector, initialMem); %#ok
for parNr = 1:D
if ~isnan(initialMem(parNr)) && paramDefCell{parNr,3} > 0 && abs(initialMem(parNr) - quantInitialMem(parNr)) > eps
disp(textwrap2(sprintf(['Warning: Initial value of parameter %s (%g) ', ...
'not on quantization grid (next grid value: %g).'], getparametername(parNr, 1), ...
initialMem(parNr), quantInitialMem(parNr)), textWidth));
end
end
end
end

% display title
disp(repmat('*', 1, textWidth)); % ***************
if isstruct(optimInfo) && isfield(optimInfo, 'title')
disp(textwrap2(sprintf('Starting ''%s'' at %s.', optimInfo.title, ...
datestr(clock, 'mmm dd, HH:MM')), textWidth));
if isfield(optimInfo, 'subtitle')
disp(textwrap2(optimInfo.subtitle));
end
else
fprintf('Starting optimization at %s.\n', datestr(clock, 'mmm dd, HH:MM'));
end

% display number of possible parameter vectors
if isfinite(nrOfPossibleMembers)
fprintf('Number of possible parameter vectors: %g\n', nrOfPossibleMembers);
else
fprintf('Infinite number of possible parameter vectors.\n');
end

% quit if maxtime is less or equal zero (can be used to check if all
% initial parameter values are on the quantization grid)
if (DEParams.maxtime <= 0) || (~isempty(maxclock) && etime(clock, maxclock) > 0)
disp(repmat('*', 1, textWidth)); % ***************
disp(textwrap2(sprintf(['Function %s was left because parameter maxtime ', ...
'is less or equal zero.\n'], mfilename), textWidth));
if nargout > 0
varargout = {[], [], [], 0};
end
return
end

% get slave file directory
if isempty(objFctHandle)
% do not generate slave files in testing mode
feedSlaveProc = false;
end
if feedSlaveProc
% build name of slave file directory
if isempty(slaveFileDir)
slaveFileDir = concatpath(tempdir2, 'differentialevolution');
end
% create slave file directory if not existing
if ~exist(slaveFileDir, 'dir')
mkdir(slaveFileDir);
end
else
slaveFileDir = '';
end

% save 'time over'-file
if ~isempty(objFctHandle)
objFctName = strrep(func2str(objFctHandle), '/', '_');
else
objFctName = 'test';
end
timeOverFileName = sprintf('%s_timeover_%s.mat', objFctName, gethostname);
if exist('./data', 'dir')
timeOverFileName = ['data/' timeOverFileName];
end
info = 'Delete this file to stop parameter optimization.'; %#ok
sem = setfilesemaphore(timeOverFileName);
save(timeOverFileName, 'info');
removefilesemaphore(sem);

% initialize variables
pop                = zeros(NP,D); % population array
val                = nan  (1,NP); % cost vector
valstddev          = inf;         % cost vector standard deviation
paramstddev        = inf;         % mean parameter standard deviation
nfeval.overall     = 0;           % number of function evaluations
nfeval.last        = 0;           % number of function evaluations before current iteration
nfeval.saved       = 0;           % number of evaluations that were saved
nfeval.slave       = 0;           % number of evaluations performed by slave process
nofevalCounter     = 0;           % number of iterations without any function evaluation
noChangeCounter    = 0;           % number of iterations without any population change
bestvalhist        = [];          % history of best values
bestmemhist        = [];          % history of best member
valstddevhist      = [];          % history of cost vector standard deviation
paramstddevhist    = [];          % history of mean parameter standard deviation
allval             = [];          % vector with all computed cost values
allmem             = [];          % matrix with all tested vectors as columns
% note: parameter vectors are saved in pop as rows, in allmem as columns!

if DEParams.minimizeValue
minMaxSign =  1; % minimization is default behavior
else
minMaxSign = -1;
end

% clear persistent variables in subfunctions
clear sendmailblat getparametername timeovercheck displayprogressinfo
clear displaybestmember
saveoptimizationresult();
computeevaluationvalue();

% initialize function sendmailblat
sendEmail = sendmailblat([], [], emailParams);

% save current time
startTime           = mbtime;
nextRefreshTime     = startTime + refreshtime;
nextRefreshTime2    = startTime + refreshtime2;
lastRefreshIterTime = -inf;

% save current (empty) state
state = 'Empty state before first iteration.';
if saveHistory
saveoptimizationresult(paramDefCell, parameterDimVector, optimInfo, [], [], [], [], 0, 0, ...
[], [], startTime, state, DEParams, [], [], objFctName, ...
objFctSettings, objFctParams, 1, 0);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% start main program loop %
%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(repmat('*', 1, textWidth)); % ***************
iterationNr         = 0;
bestval             = minMaxSign * Inf;
displayMeanEvalTime = true;
timeOver            = false;
while ~timeOver && (iterationNr < DEParams.maxiter) && ...
(isempty(VTR) || isnan(VTR) || (bestval * minMaxSign > VTR * minMaxSign))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize or re-initialize population %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
initialization   = (iterationNr == 0);
reinitialization = (valstddev < DEParams.minvalstddev) || (paramstddev < DEParams.minparamstddev) || ...
(noChangeCounter >= DEParams.nochangeiter) || (nofevalCounter > DEParams.nofevaliter);

if initialization

% initialization
if DEParams.useInitParams > 0
% first population member: current parameters
pop(1,:) = initialMem;
istart = 2;
else
% initialize all population members randomly
istart = 1;
end

if nrOfPossibleMembers <= NP && D == 1 && paramDefCell{1,3} ~= 0
% if only one scalar parameter to optimize, set all possible members
% directly (no random initialization)

if DEParams.useInitParams > 0
% get and quantize first member
[ignore, firstMem] = considerparametercontraints(...
objFctParams, paramDefCell, parameterDimVector, pop(1,:)); %#ok
end

% generate equidistant vector
pop = (XVmin(1) : paramDefCell{1,3} :  (XVmax(1) + 0.5*paramDefCell{1,3}))';

if DEParams.useInitParams > 0
% evaluate initial member first, after that those
% that are next to initial member
[ignore, sortIndex] = sort(abs(pop(:,1) - firstMem)); %#ok
pop = pop(sortIndex,:);

% do not use random initialization later
istart = NP+1;
end
pop = pop(1:NP,:);

elseif DEParams.useInitParams == 2
% next population members: current parameter vector plus random noise
pop      = computerandominitialization(1, pop, memIndex, paramDefCell, ...
objFctSettings, parameterDimVector, XVmax, XVmin, validChkHandle);
end
end
initialization = true;

elseif reinitialization

% re-initialization
fprintf('Re-initializing population in iteration %d.\n', iterationNr);
if valstddev < DEParams.minvalstddev
fprintf('Evaluation value std. dev. below threshold (%g < %g).\n', ...
valstddev, DEParams.minvalstddev);
elseif paramstddev < DEParams.minparamstddev
fprintf('Mean parameter std. dev. below threshold (%g < %g).\n', ...
paramstddev, DEParams.minparamstddev);
elseif noChangeCounter >= DEParams.nochangeiter
fprintf('Population did not change in last %d iterations.\n', ...
DEParams.nochangeiter);
elseif nofevalCounter > DEParams.nofevaliter
fprintf('No function evaluations in the last %d iterations.\n', ...
DEParams.nofevaliter);
end
disp(' ');

if minMaxSign > 0
[val, index] = sort( val);
else
[val, index] = sort(-val); % not using 'descend' for backwards compatibility
val = -val;
end
istart = max(1, round(0.2*NP)); % keep best twenty percent of the population members
index = index(1:istart-1);
pop(1:istart-1,:) = pop(index,:);
val(1:istart-1)   = val(index);
reinitialization  = true;
end

if initialization || reinitialization

% initialize population members from istart to NP randomly
pop = computerandominitialization(2, pop, istart:NP, paramDefCell, ...
objFctSettings, parameterDimVector, XVmax, XVmin, validChkHandle);

% feed slave process
if feedSlaveProc
generatefilesforslaveprocess(objFctHandle, objFctParams, ...
objFctSettings, paramDefCell, parameterDimVector, pop, allmem, ...
iterationNr, saveHistory, slaveFileDir, validChkHandle);
end

% evaluate initial population members
for memberNr = 1:NP
if initialization && memberNr == 1
startLoopTime = mbtime;
end

% evaluate member
[val(memberNr), pop(memberNr,:), nfeval, allval, allmem] = ...
computeevaluationvalue(pop(memberNr,:), nfeval, timeOver, ...
objFctParams, paramDefCell, parameterDimVector, objFctSettings, ...
objFctHandle, saveHistory, allval, allmem, iterationNr, ...
memberNr, NP, slaveFileDir, maxMasterEvals, validChkHandle, XVmin, XVmax);

if initialization && (memberNr == 1)

% set initial value
bestval = val(1);
if isnan(bestval)
% set to highest possible value
bestval = minMaxSign * Inf;
elseif isinf(bestval)
if (bestval > 0) && (minMaxSign < 0)
error('Objective function value of initial parameter vector may not be +Inf in maximization problem.');
elseif (bestval < 0) && (minMaxSign > 0)
error('Objective function value of initial parameter vector may not be -Inf in minimization problem.');
end
end
initval = bestval;

% best member so far
bestmem = pop(1,:);
initmem = bestmem;

% display initial member
state = 'Initial parameter set.';
displaybestmember(paramDefCell, parameterDimVector, initval, initmem, ...
allval, initval, initmem, 1, state, optimInfo, sendEmail, playSound);

% save initial state
if saveHistory
saveoptimizationresult(paramDefCell, parameterDimVector, optimInfo, ...
bestval, bestmem, bestval, bestmem, 0, 0, pop, val, startTime, ...
state, DEParams, allval, allmem, objFctName, ...
objFctSettings, objFctParams, 0, 0);
end

% display time needed for first function evaluation
if mbtime - startLoopTime > 1
fprintf('First function evaluation time: %s\n', ...
formattime(mbtime - startLoopTime));
end

else

% check if current member is overall best
if isfinite(val(memberNr)) && ...
(val(memberNr) * minMaxSign < bestval * minMaxSign)
bestval = val(memberNr);
bestmem = pop(memberNr,:);
end

% set current state
if initialization
state = sprintf('In initialization, %d of %d members checked.', ...
memberNr, NP);
else
state = sprintf('In iteration %d, %d of %d members re-initialized.', ...
iterationNr, memberNr, NP);
end
end

% display best member
displaybestmember(paramDefCell, parameterDimVector, bestval, bestmem, ...
allval, initval, initmem, 0, state, optimInfo, sendEmail, playSound);

% save and display current state
if mbtime - nextRefreshTime >= 0
nextRefreshTime = nextRefreshTime + refreshtime * ...
(1 + floor((mbtime - nextRefreshTime) / refreshtime));
disp(state);
if saveHistory
saveoptimizationresult(paramDefCell, parameterDimVector, optimInfo, ...
bestval, bestmem, bestval, bestmem, 0, 0, pop, val, startTime, ...
state, DEParams, allval, allmem, objFctName, ...
objFctSettings, objFctParams, 0, 0);
end
end

% display progress information
if mbtime - nextRefreshTime2 >= 0
nextRefreshTime2 = nextRefreshTime2 + refreshtime2 * ...
(1 + floor((mbtime - nextRefreshTime2) / refreshtime2));
displayprogressinfo(startLoopTime, state, refreshtime3, ...
maxtime, maxclock, timeOver, nfeval, nrOfPossibleMembers, pop, ...
bestval, allval, optimInfo, sendEmail);
end

% display mean function evaluation time
if displayMeanEvalTime && nfeval.overall == 5 && mbtime - startLoopTime > 1
fprintf('Mean function evaluation time after %d runs: %s\n', ...
nfeval.overall, formattime((mbtime - startLoopTime) / nfeval.overall));
displayMeanEvalTime = false;
end

% check time
timeOver = timeovercheck(startTime, maxtime, maxclock, timeOverFileName, timeOver);

end % for memberNr = 1:NP

iterationNr = iterationNr + 1;
end % if initialization || reinitialization

if timeOver
break
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute competing population %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% In function computenewpopulation, you can place your own favorite algorithm!

popold = pop;
popnew = computenewpopulation(pop, bestmem, DEParams);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% check which vectors will survive %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% feed slave process
if feedSlaveProc
generatefilesforslaveprocess(objFctHandle, objFctParams, ...
objFctSettings, paramDefCell, parameterDimVector, popnew, allmem, ...
iterationNr, saveHistory, slaveFileDir, validChkHandle);
end

for memberNr = 1:NP

% check cost of competitor
[tempval, popnew(memberNr,:), nfeval, allval, allmem] = ...
computeevaluationvalue(popnew(memberNr,:), nfeval, timeOver, ...
objFctParams, paramDefCell, parameterDimVector, objFctSettings, ...
objFctHandle, saveHistory, allval, allmem, iterationNr, ...
memberNr, NP, slaveFileDir, maxMasterEvals, validChkHandle);

if isfinite(tempval)
% if competitor is better than value in cost array or old vector was
% not evaluated (tempval = NaN), replace old vector with new one (for
% new iteration) and save value in cost array
if (tempval * minMaxSign < val(memberNr) * minMaxSign) || ~isfinite(val(memberNr))
pop(memberNr,:) = popnew(memberNr,:);
val(memberNr)   = tempval;

% if competitor is better than the best one ever, set new best value
% and set new best parameter vector ever
if (tempval * minMaxSign < bestval * minMaxSign)
bestval = tempval;
bestmem = popnew(memberNr,:);
end
end
else
if isinf(tempval)
if (tempval > 0) && (minMaxSign < 0)
warning('Objective function value of initial parameter vector may not be +Inf in maximization problem.'); %#ok
elseif (tempval < 0) && (minMaxSign > 0)
warning('Objective function value of initial parameter vector may not be -Inf in minimization problem.'); %#ok
end
end
end

% display best member
state = sprintf('In iteration %d, %d of %d competitors checked.', iterationNr, memberNr, NP);
displaybestmember(paramDefCell, parameterDimVector, bestval, bestmem, ...
allval, initval, initmem, 0, state, optimInfo, sendEmail, playSound);

% save and display current state
if mbtime - nextRefreshTime >= 0
nextRefreshTime = nextRefreshTime + refreshtime * ...
(1 + floor((mbtime - nextRefreshTime) / refreshtime));
disp(state);
if saveHistory
saveoptimizationresult(paramDefCell, parameterDimVector, optimInfo, ...
bestval, bestmem, bestvalhist, bestmemhist, valstddevhist, ...
paramstddevhist, pop, val, startTime, state, DEParams, allval, allmem, ...
objFctName, objFctSettings, objFctParams, 0, 0);
end
end

% display progress information
if mbtime - nextRefreshTime2 >= 0
nextRefreshTime2 = nextRefreshTime2 + refreshtime2 * ...
(1 + floor((mbtime - nextRefreshTime2) / refreshtime2));
displayprogressinfo(startLoopTime, state, refreshtime3, maxtime, ...
maxclock, timeOver, nfeval, nrOfPossibleMembers, pop, bestval, ...
allval, optimInfo, sendEmail);
end

% check time
timeOver = timeovercheck(startTime, maxtime, maxclock, timeOverFileName, timeOver);

end % for memberNr = 1:NP

%%%%%%%%%%%%%%%%%%%%
% finish iteration %
%%%%%%%%%%%%%%%%%%%%

% compute cost value and population standard deviation
index = isfinite(val);
if any(index)
valstddev   = std(val(index));
paramstddev = max(std(pop(index,:),0,1)' ./ diff(cell2mat(paramDefCell(:,2)),1,2));
else
valstddev   = inf;
paramstddev = inf;
end

% check if population has changed
if isequal(pop, popold)
noChangeCounter = noChangeCounter + 1;
else
noChangeCounter = 0;
end

% check if any function evaluation was done
if nfeval.last == nfeval.overall
nofevalCounter = nofevalCounter + 1;
else
nofevalCounter = 0;
end

% save history
if saveHistory
bestvalhist    (end+1) = bestval;     %#ok
bestmemhist  (:,end+1) = bestmem';    %#ok
valstddevhist  (end+1) = valstddev;   %#ok
paramstddevhist(end+1) = paramstddev; %#ok

% check if all possible members have been computed
if length(allval) >= nrOfPossibleMembers
timeOver = true;
end
end

% check time
if timeOver
break
end

% display current state
if ((DEParams.refreshiter > 0) && (rem(iterationNr, DEParams.refreshiter) == 0)) && ...
mbtime - lastRefreshIterTime > 10
lastRefreshIterTime = mbtime; % avoid many lines of output if no evaluations were done
fprintf('Iteration %d finished.\n', iterationNr);
end

iterationNr = iterationNr + 1;
end % while (iterationNr < DEParams.maxiter) ...

%%%%%%%%%%%%%%%%%%%%%%%%
% display final result %
%%%%%%%%%%%%%%%%%%%%%%%%
disp(' ');
disp(repmat('*', 1, textWidth)); % ***************

% display why optimization was finished
if iterationNr >= DEParams.maxiter
state = sprintf('''%s'' finished after given maximum number of %d iterations.', ...
optimInfo.title, DEParams.maxiter);
elseif ~isempty(VTR) && (bestval * minMaxSign <= VTR * minMaxSign)
state = sprintf('''%s'' finished after specified cost value of %2.6g was reached.', ...
optimInfo.title, VTR);
elseif length(allval) >= nrOfPossibleMembers
state = sprintf('''%s'' finished after all possible %d members have been tested.', ...
optimInfo.title, nrOfPossibleMembers);
elseif timeOver
if ~isempty(maxtime) && mbtime - startTime > maxtime
state = sprintf('''%s'' finished after given amount of time.', optimInfo.title);
elseif ~isempty(maxclock) && etime(clock, maxclock) > 0
state = sprintf('''%s'' finished at given end time.', optimInfo.title);
elseif ~isempty(timeOverFileName) && ~existfile(timeOverFileName);
state = sprintf('''%s'' finished after ''time over''-file was deleted.', optimInfo.title);
end
end
disp(textwrap2(state, textWidth));
displayprogressinfo(startLoopTime, state, [], maxtime, maxclock, 1, ...
nfeval, nrOfPossibleMembers, pop, bestval, allval, optimInfo, sendEmail);
state = sprintf('Final result after %d iteration(s).', iterationNr);
displaybestmember(paramDefCell, parameterDimVector, bestval, bestmem, allval, ...
initval, initmem, 1, state, optimInfo, sendEmail, playSound);

% save final result
if saveHistory
saveoptimizationresult(paramDefCell, parameterDimVector, optimInfo, ...
bestval, bestmem, bestvalhist, bestmemhist, valstddevhist, paramstddevhist, ...
pop, val, startTime, state, DEParams, allval, allmem, objFctName, ...
objFctSettings, objFctParams, 1, 1);
end
disp(repmat('*', 1, textWidth)); % ***************

% display optimization parameter history
if displayResults && saveHistory
displayoptimizationhistory(paramDefCell, allmem, allval);
end

% remove all remaining slave files
if exist(slaveFileDir, 'dir')
remainingSlaveFiles = findfiles(slaveFileDir, 'iteration_*_member_*_*.mat', 'nonrecursive');
deletewithsemaphores(remainingSlaveFiles);
end

% remove time-over file
if existfile(timeOverFileName)
delete(timeOverFileName);
end

if nargout > 0
% compute parameter set with best member
[bestFctParams, bestmem] = considerparametercontraints(...
objFctParams, paramDefCell, parameterDimVector, bestmem);
varargout = {bestmem', bestval, bestFctParams, iterationNr};
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function timeOver = timeovercheck(startTime, maxtime, maxclock, timeOverFile, timeOver)

persistent lastCheckTime
if isempty(lastCheckTime)
lastCheckTime = mbtime;
end

if timeOver
return
end

curTime = mbtime;
if curTime - lastCheckTime > 1
% only check once a second to save computation time
timeOver = ...
(~isempty(maxtime)      && ((curTime - startTime) > maxtime || maxtime == 0)) || ...
(~isempty(maxclock)     &&  etime(clock, maxclock) > 0) || ...
(~isempty(timeOverFile) && ~existfile(timeOverFile));
lastCheckTime = curTime;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function valid = alwaysvalid(varargin) %#ok

valid = true;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function checkinputs(paramDefCell, objFctParams)

textWidth = 75;

% check dimensions of paramDefCell
if ~iscell(paramDefCell)
error('Input argument paramDefCell has to be a cell array.');
end
if isempty(paramDefCell)
error('Input argument paramDefCell must not be empty.');
end
if size(paramDefCell, 2) < 3 || size(paramDefCell, 2) > 4
error('Input argument paramDefCell has to consist of three or four columns.');
end

% check parameter names
for m = 1:size(paramDefCell, 1)
for n = m+1:size(paramDefCell, 1)
if strcmp(paramDefCell{m,1}, paramDefCell{n,1})
error(['Parameter names in paramDefCell have to be unique. Parameter ', ...
'name %s was found twice.'], paramDefCell{m,1});
end
end
end

% check dimensions of cell contents
for k = 1:size(paramDefCell, 1)
if (k > 1 && (isempty((paramDefCell{k,1})) || ~ischar(paramDefCell{k,1}))) || ...
k == 1 && ~ischar(paramDefCell{k,1})
error(['All cells in the first column of paramDefCell have to contain ', ...
'non-empty strings (except when there is only one row).']);
end
if isempty(paramDefCell{k,2}) || size(paramDefCell{k,2}, 2) ~= 2
error(textwrap2(['All cells in the second column of paramDefCell have to ', ...
'contain matrices with two columns (the parameter limits).'], textWidth));
end
if any(~isfinite(paramDefCell{k,2}))
error(textwrap2(['The parameter limit matrices may not contain Inf or NaN. ', ...
'You have to provide hard parameter bounds in any case, sorry.'], textWidth));
end
if isempty(paramDefCell{k,3}) || size(paramDefCell{k,3}, 2) ~= 1
error(textwrap2(['All cells in the third column of paramDefCell have to ', ...
'contain scalars or column vectors (the parameter quantizations).'], textWidth));
end
if size(paramDefCell{k,2}, 1) ~= size(paramDefCell{k,3}, 1)
error(['All vectors or matrices in the second, third and fourth row ', ...
'of paramDefCell have to have the same number of rows.']);
end
if size(paramDefCell, 2) == 4
if ~isempty(paramDefCell{k,4}) && size(paramDefCell{k,4}, 2) ~= 1
error(['All cells in the fourth column of paramDefCell have to be ', ...
'empty or contain scalars or column vectors (the initial values).']);
end
if size(paramDefCell{k,2}, 1) ~= size(paramDefCell{k,4}, 1)
error(['All vectors or matrices in the second, third and fourth row ', ...
'of paramDefCell have to have the same number of rows.']);
end
end
end

% check dimensions of objFctParams
if ~isempty(objFctParams) && isstruct(objFctParams)
fieldNames = fieldnames(objFctParams);
for k=1:length(fieldNames)
if ~isempty(objFctParams.(fieldNames{k})) && ...
size(objFctParams.(fieldNames{k}), 2) ~= 1
error('Only column vectors are allowed as parameters in objFctParams, sorry.');
end
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function str = getparametername(varargin)

persistent paramDefCell parameterDimVector

if iscell(varargin{1})
% initialization
paramDefCell       = varargin{1};
parameterDimVector = varargin{2};
return
else
% normal operation
parNr    = varargin{1};
nameMode = varargin{2};
end

if strcmp(paramDefCell{1,1}, '_1')
str = sprintf('%d', parNr);
elseif parameterDimVector(parNr) > 1
switch nameMode
case 1
% return for example "bla(2)" for parameter name "bla_2"
str = regexprep(paramDefCell{parNr,1}, '_(\d)+\$', '(\$1)');
case 2
% return for example "bla" for parameter name "bla_2"
str = regexprep(paramDefCell{parNr,1}, '_\d+\$', '');
otherwise
error('Name mode %d unknown.', nameMode);
end
else
str = paramDefCell{parNr,1};
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function pop = computerandominitialization(randMode, pop, memIndex, ...
paramDefCell, objFctSettings, parameterDimVector, XVmax, XVmin, validChkHandle)

if isempty(memIndex)
return
end

switch randMode
case 1
% use first population member plus random noise
baseMem = pop(1,:);
randStdDev1   = 0.1;
case 2
% only use random numbers
baseMem = XVmin;
randStdDev1   = 1;
otherwise
error('Random initialization mode %d unknown.', randMode);
end

% initialize population randomly
D = size(pop, 2);
for n = memIndex
pop(n,:) = baseMem + randStdDev1 * rand(1,D) .* (XVmax - XVmin);
end

% quantize all population vectors
quantPop = pop;
objFctParamsCell = cell(size(pop,1),1);
for n=1:size(pop,1)
[objFctParamsCell{n}, quantPop(n,:)] = considerparametercontraints([], ...
paramDefCell, parameterDimVector, pop(n,:)); %#ok
end

% check for multiple occurences and invalid parameter vectors and recompute
% random vectors
nindex = find(memIndex > 1, 1);
maxNrOfTests = min(1000*D, 10*length(memIndex));
nrOfRecomputations = 0;
for k=1:maxNrOfTests
if nindex == length(memIndex)
break
end
n = memIndex(nindex);

if any(all(abs(quantPop(1:n-1,:) - quantPop(n(ones(n-1,1)), :)) < eps, 2)) || ...
~paramvecvalidity(paramDefCell, objFctSettings, objFctParamsCell{n}, ...
quantPop(n,:), validChkHandle)
% quantized member invalid or already in population, recompute member

% increase random standard deviation
nrOfRecomputations = nrOfRecomputations + 1;

% compute new random parameter vector
randmem = baseMem + randStdDev*rand(1,D).*(XVmax - XVmin);

% consider parameter quantization
[objFctParamsCell{n}, quantPop(n,:), quantmem2] = considerparametercontraints([], ...
paramDefCell, parameterDimVector, randmem); %#ok
pop(n,:) = quantmem2;
else
% quantized member unique, step to next member
nindex = nindex + 1;
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function displayprogressinfo(startLoopTime, state, refreshtime3, ...
maxtime, maxclock, timeOver, nfeval, nrOfPossibleMembers, pop, bestval, ...
allval, optimInfo, sendEmail) %#ok

persistent nextRefreshTime3

% elapsed time
elapsedTime = mbtime - startLoopTime;
str = sprintf('Elapsed time:                   %s\n', formattime(elapsedTime));

% time left
if ~timeOver
timeLeft = inf;
if ~isempty(maxtime) && isfinite(maxtime)
timeLeft = min(timeLeft, maxtime - elapsedTime);
end
if ~isempty(maxclock)
timeLeft = min(timeLeft, -etime(clock, maxclock));
end
if isfinite(timeLeft) && timeLeft > 0
str = [str, sprintf('Time left:                      %s\n', formattime(timeLeft))];
end
end

% function evaluations
str = [str, sprintf('Number of function evaluations: %d\n',   nfeval.overall)];
if nfeval.slave > 0
str = [str, sprintf('Number of slave evaluations:    %d\n',   nfeval.slave)];
end
percentage = round(nfeval.overall / nrOfPossibleMembers * 100);
if percentage > 0
str = [str, sprintf('Percentage of members computed: %d %%\n', percentage)];
end
if nfeval.overall - nfeval.slave > 0
str = [str, sprintf('Mean function evaluation time:  %s\n', ...
formattime(elapsedTime / (nfeval.overall - nfeval.slave)))];
end

if ~isempty(allval)
sameEvaluationValue = length(find(allval == bestval));
if sameEvaluationValue > 1
str = [str, sprintf('Vectors with best value:        %d\n', sameEvaluationValue)];
end
end

disp(' ');
disp(str);

% send E-mail
if sendEmail && ~isempty(refreshtime3)
if isempty(nextRefreshTime3)
nextRefreshTime3 = mbtime + refreshtime3;
elseif mbtime - nextRefreshTime3 >= 0
if ~isempty(refreshtime3)
nextRefreshTime3 = nextRefreshTime3 + refreshtime3 * ...
(1 + floor((mbtime - nextRefreshTime3) / refreshtime3));
end
subject = 'Progress information';
if isfield(optimInfo, 'title')
subject = [subject, sprintf(' (%s, host %s)', optimInfo.title, gethostname)];
else
subject = [subject, sprintf(' (host %s)', gethostname)];
end
sendmailblat(subject, [state, sprintf('\n'), str]);
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function saveoptimizationresult(paramDefCell, parameterDimVector, optimInfo, ...
bestval, bestmem, bestvalhist, bestmemhist, valstddevhist, paramstddevhist, ...
pop, val, startTime, state, DEParams, allval, allmem, objFctName, ...
objFctSettings, objFctParams, forceResultFileNameDisplay, ...

if nargin == 0
resultFileName     = '';
return
end

D = size(paramDefCell, 1);
hostname = gethostname;

% save all interesing values in a structure using more meaningful variable names
optimResult.title = optimInfo.title;
optimInfo = rmfield(optimInfo, 'title');
if ~isempty(fieldnames(optimInfo))
optimResult.info = optimInfo;
end
optimResult.state                    = state;
optimResult.startTime                = datestr(mbdatevec(startTime), 31);
optimResult.startClock               = mbdatevec(startTime);
optimResult.currentTime              = datestr(clock, 31);
optimResult.DEParams                 = DEParams;
optimResult.paramDefCell             = paramDefCell;
optimResult.hostname                 = hostname;
optimResult.bestMember               = [];         % to be overwritten
optimResult.objFctParams             = objFctParams;
optimResult.boundaryValuesReached    = zeros(D,1); % to be overwritten
optimResult.bestEvaluationValue      = bestval;
optimResult.bestMemberHistory        = bestmemhist;
optimResult.bestValueHistory         = bestvalhist;
optimResult.costValueVarianceHistory = valstddevhist;
optimResult.parameterStdDevHistory   = paramstddevhist;
optimResult.currentPopulation        = pop';
optimResult.currentCostValues        = val;
optimResult.allEvaluationValues      = allval;
optimResult.allTestedMembers         = allmem;

% overwrite values in objFctParams with best member
if ~isempty(bestmem)
[ignore, bestmem] = considerparametercontraints([], paramDefCell, ...
parameterDimVector, bestmem); %#ok
optimResult.bestMember = bestmem';
if strcmp(paramDefCell{1,1}, '_1')
optimResult = rmfield(optimResult, 'objFctParams');
else
parNr = 1;
while parNr <= D
index = parNr:parNr+parameterDimVector(parNr)-1;
optimResult.objFctParams.(getparametername(parNr, 2)) = bestmem(index);
parNr = parNr + parameterDimVector(parNr);
end
end
for parNr = 1:D
optimResult.boundaryValuesReached(parNr) = any(bestmem(parNr) == paramDefCell{parNr,2});
end
end
optimResult.objFctSettings = objFctSettings;

% get file number to avoid overwriting old results
if isempty(resultFileName)
fileName = sprintf('%s_lastresultnumber.mat', objFctName);
if exist('./data', 'dir')
fileName = ['data/', fileName];
end

% save current result file number
sem = setfilesemaphore(fileName);
if existfile(fileName)
resultFileNr = mod(resultFileNr, 50) + 1; %#ok
else
resultFileNr = 1;
end
save(fileName, 'resultFileNr');
removefilesemaphore(sem);

% build file name
resultFileName = sprintf('%s_result_%s_%02d.mat', objFctName, hostname, resultFileNr);
if exist('./data', 'dir')
resultFileName = ['data/', resultFileName];
end
forceResultFileNameDisplay = true;
end

if forceResultFileNameDisplay
fprintf('Results are saved in file %s.\n', resultFileName);
end

% save data
sem = setfilesemaphore(resultFileName);
save(resultFileName, 'optimResult');
removefilesemaphore(sem);

if 0 % && ispc
% the following code is deactivated, as it uses a Perl script and an
% external FTP program to upload the current results to a server. Contact
% me if you are interested in this script.
fileUploadPeriod = 15*60; % in seconds
end

try %#ok
% put file to server
% (this is a file access that is not protected by a semaphore, but
% different Matlab processes use different result file names)
system(sprintf('start /B putfiletoserver.pl %s', resultFileName));
end
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function displaybestmember(paramDefCell, parameterDimVector, bestval, bestmem, ...
allval, initval, initmem, forceParameterDisplay, state, optimInfo, sendEmail, playSound)

persistent lastbestmem lastSoundTime lastEmailTime lastSubject lastStr

minTimeBetweenEmails = 30; % avoid sending many E-mails shortly after another

D = size(paramDefCell, 1);

if nargin == 1
% initialize values
lastbestmem   = NaN;
lastEmailTime = -inf;
lastSoundTime = mbtime;
hostname = gethostname;

% get maximum name length for proper display
maxNameLength = 0;
for parNr = 1:D
maxNameLength = max(maxNameLength, length(getparametername(parNr, 1)));
end
return
end

% display current best member if it has changed or if display is forced
if ~any(isnan(lastbestmem)) && any(size(bestmem) ~= size(lastbestmem))
error('Internal error: bestmem and lastbestmem are of different size!');
end
if any(bestmem ~= lastbestmem) || forceParameterDisplay
sendEmailThisTime = 1;
lastbestmem = bestmem;

% display state
str = sprintf('%s\n', state);

% get quantized parameter vector for display
[ignore, bestmem] = considerparametercontraints([], paramDefCell, ...
parameterDimVector, bestmem); %#ok

case 'Markus@Edison'
prefixString = 'param.';
otherwise
prefixString = '';
end
str = [str, sprintf('Best member:\n')];
if all(bestmem == initmem) && isempty(strfind(state, 'Initial'))
str = [str(1:end-2), sprintf(' (same as initial member):\n')];
end
for parNr=1:D
if any(bestmem(parNr) == paramDefCell{parNr, 2});
markStr = '% (boundary value)';
else
markStr = '';
end
if strcmp(paramDefCell{1,1}, '_1')
str = [str, sprintf('%10g; %s\n', bestmem(parNr), markStr)]; %#ok
else
parameterName = getparametername(parNr, 1);
str = [str, sprintf('%s%s = %g; %s\n', prefixString, ...
[parameterName, repmat(' ', 1, maxNameLength - ...
length(parameterName))], bestmem(parNr), markStr)]; %#ok
end
end
if bestval >= 1e-5
str = [str, sprintf('Evaluation value: %2.6f\n', bestval)]; % always print zeros
else
str = [str, sprintf('Evaluation value: %2.6g\n', bestval)]; % use exponential notation for small values
end
if bestval == initval && isempty(strfind(state, 'Initial'))
str = [str(1:end-1), sprintf(' (same as initial value)\n')];
end

if ~isempty(allval)
sameEvaluationValue = length(find(allval == bestval));
if sameEvaluationValue > 1
str = [str, sprintf('Number of vectors with same evaluation value: %d\n', sameEvaluationValue)];
end
end

% display information
disp(str);

% play sound
if playSound && ~isempty(lastSoundTime) && mbtime - lastSoundTime > 60
soundsc(x, fs, bits);
pause(length(x) / fs + 1);
lastSoundTime = mbtime;
end
else
sendEmailThisTime = 0;
end

if ~sendEmail
return
end

if sendEmailThisTime
% build subject and body
if bestval >= 1e-5
formatString = '%2.6f'; % always print zeros
else
formatString = '%2.6g'; % use exponential notation for small values
end
if forceParameterDisplay
if ~isempty(strfind(lower(state), 'initial'))
subject = sprintf(sprintf('Initial eval. value: %s', formatString), bestval);
else
subject = sprintf(sprintf('Best eval. value: %s', formatString), bestval);
end
else
subject = sprintf(sprintf('New best eval. value: %s', formatString), bestval);
end
if isfield(optimInfo, 'title')
subject = [subject, sprintf(' (%s, host %s)', optimInfo.title, gethostname)];
else
subject = [subject, sprintf(' (host %s)', gethostname)];
end

if mbtime - lastEmailTime < minTimeBetweenEmails
% do not send now, save data for sending later
lastSubject       = subject;
lastStr           = str;
sendEmailThisTime = 0;
end

elseif ~isempty(lastSubject) && mbtime - lastEmailTime >= minTimeBetweenEmails
% restore data that was not sent until now
subject           = lastSubject;
str               = lastStr;
sendEmailThisTime = 1;
end

if sendEmailThisTime
sendmailblat(subject, str);
lastEmailTime = mbtime;
lastSubject   = [];
lastStr       = [];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [objFctParams, paramVec, paramVec2] = considerparametercontraints( ...
objFctParams, paramDefCell,  parameterDimVector, paramVec)

paramVec2 = paramVec; % paramVec2: not quantized

for parNr = 1:length(paramVec)
% hard bound constraint
parBounds = paramDefCell{parNr, 2};
paramVec (parNr) = min(max(paramVec(parNr), parBounds(1)), parBounds(2));
paramVec2(parNr) = paramVec(parNr);

% parameter quantization
minQuantization = 1e-14;
parGrid = paramDefCell{parNr, 3};
parGrid = max(parGrid, minQuantization); % quantize with at least some value larger eps, see also checkinputs
paramVec(parNr) = parGrid * round((paramVec(parNr) - parBounds(1)) / ...
parGrid) + parBounds(1);
end

% set parameter vector in objFctParams
% (always use quantized value here!)
if ~strcmp(paramDefCell{1,1}, '_1')
parNr = 1;
while parNr <= length(paramVec)
index = parNr:parNr+parameterDimVector(parNr)-1;
objFctParams.(getparametername(parNr,2)) = paramVec(index)';
parNr = parNr + parameterDimVector(parNr);
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function generatefilesforslaveprocess(objFctHandle, objFctParams, ...
objFctSettings, paramDefCell, parameterDimVector, pop, allmem, ...
iterationNr, saveHistory, slaveFileDir, validChkHandle)

if ~exist(slaveFileDir, 'dir')
return
end

% remove all existing slave files
existingSlaveFiles = findfiles(slaveFileDir, 'iteration_*_member_*_*.mat', 'nonrecursive');
deletewithsemaphores(existingSlaveFiles);

% build slave file name
slaveFileNameTemplate = concatpath(slaveFileDir, ...
sprintf('iteration_%02d_member_XX_parameters.mat', iterationNr));

% generate new slave files
NP = size(pop,1);
if saveHistory
allmem = [allmem, nan(size(pop,2), NP)];
else
allmem = nan(size(pop,2), NP);
end
nrOfCols = size(allmem,2);

for memberNr = NP:-1:1
testmem = pop(memberNr,:);

% get constrained parameter vector
[objFctParams, testmem] = considerparametercontraints(...
objFctParams, paramDefCell, parameterDimVector, testmem); %#ok

if ~paramvecvalidity(paramDefCell, objFctSettings, objFctParams, testmem, validChkHandle)
% parameter vector invalid
continue
end

% check if the current parameter vector was tested before
index = find(all(abs(allmem - repmat(testmem', 1, nrOfCols)) < eps, 1));

if length(index) > 1
disp('Warning: More than one equal test vector in allmem (internal error?).');
elseif ~isempty(index)
continue
end

% save testmem in allmem, so that no two files with the same parameters are saved
allmem(:,nrOfCols-memberNr+1) = testmem';

% build cell array of function arguments
if strcmp(paramDefCell{1,1}, '_1')
% pass parameters as vector
if iscell(objFctSettings)
argumentCell = [objFctSettings, {testmem}];
else
argumentCell = {objFctSettings,  testmem};
end
else
% pass parameters as structure objFctParams
if iscell(objFctSettings)
argumentCell = [objFctSettings, {objFctParams}];
else
argumentCell = {objFctSettings,  objFctParams};
end
end

% save file
memberNrString = sprintf(sprintf('%%0%dd', ceil(log10(NP+1))), memberNr);
slaveFileName = strrep(slaveFileNameTemplate, 'XX', memberNrString);

% save slave file
sem = setfilesemaphore(slaveFileName);
objFctHandle; %#ok
argumentCell; %#ok
save(slaveFileName, 'objFctHandle', 'argumentCell');
removefilesemaphore(sem);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function valid = paramvecvalidity(paramDefCell, objFctSettings, objFctParams, ...
testmem, validChkHandle)

if strcmp(paramDefCell{1,1}, '_1')
% pass parameters as vector
if iscell(objFctSettings)
valid = validChkHandle(objFctSettings{:}, testmem');
else
valid = validChkHandle(objFctSettings,    testmem');
end
else
% pass parameters as structure objFctParams
if iscell(objFctSettings)
valid = validChkHandle(objFctSettings{:}, objFctParams);
else
valid = validChkHandle(objFctSettings,    objFctParams);
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [testval, testmem, nfeval, allval, allmem] = computeevaluationvalue(...
testmem, nfeval, timeOver, objFctParams, paramDefCell, parameterDimVector, ...
objFctSettings, objFctHandle, saveHistory, allval, allmem, iterationNr, ...
memberNr, NP, slaveFileDir, maxMasterEvals, validChkHandle, XVmin, XVmax)

persistent optParVector nrOfMasterEvaluations lastIterationNr

pauseTime    = 0.1;
maxPauseTime = 1;

if nargin == 0
lastIterationNr       = NaN;
nrOfMasterEvaluations = 0;
return
end

% get constrained parameter vector
[objFctParams, testmem, testmem2] = considerparametercontraints(...
objFctParams, paramDefCell, parameterDimVector, testmem);

% check if current parameter vector is valid
if ~paramvecvalidity(paramDefCell, objFctSettings, objFctParams, testmem, validChkHandle)
testval = NaN;

% write non-quantized value into population
testmem = testmem2;
return
end

% check if current parameter vector was tested before
if saveHistory && ~isempty(allmem)
allmemIndex = find(all(abs(allmem - repmat(testmem', 1, size(allmem, 2))) < eps, 1));
if length(allmemIndex) > 1
disp('Warning: More than one equal test vector in allmem (internal error?).');
allmemIndex = allmemIndex(1);
end
else
allmemIndex = [];
end

if ~isempty(allmemIndex)
% parameter vector was tested before
testval      = allval(allmemIndex);
nfeval.saved = nfeval.saved + 1;
else
% parameter vector was not tested before, compute evaluation value or

% build slave file name
if ~exist(slaveFileDir, 'dir')
slaveFileDir = '';
end
if lastIterationNr ~= iterationNr
lastIterationNr = iterationNr;
nrOfMasterEvaluations = 0;
end
if ~isempty(slaveFileDir)
memberNrString = sprintf(sprintf('%%0%dd', ceil(log10(NP+1))), memberNr);
slaveFileName = concatpath(slaveFileDir, sprintf('iteration_%02d_member_%s_result.mat', ...
iterationNr, memberNrString));

% remove slave parameter file if existing
if nrOfMasterEvaluations < maxMasterEvals
slaveParameterFileName = strrep(slaveFileName, 'result', 'parameters');
deletewithsemaphores(slaveParameterFileName);
end
end

if ~isempty(objFctHandle)
testval = [];

% check if the current parameter vector was tested before by slave process
if ~isempty(slaveFileDir)
while 1
sem = setfilesemaphore(slaveFileName);
if existfile(slaveFileName)
delete(slaveFileName);
end
removefilesemaphore(sem);

if ~isempty(testval)
% result of slave process was loaded
nfeval.overall = nfeval.overall + 1;
nfeval.slave   = nfeval.slave   + 1;
break
end
if nrOfMasterEvaluations < maxMasterEvals
% master will evaluate current member
break
else
% Wait until member was evaluated by slave.
% Note: If any of the slave processes is interrupted, this loop
% will never be left!! Limiting the number of master
% evaluations is not safe!!
pause(pauseTime);
pauseTime = min(2*pauseTime, maxPauseTime);
end
end
end

% evaluate objective function
if isempty(testval) && (nrOfMasterEvaluations < maxMasterEvals) && ~timeOver

% run function
if strcmp(paramDefCell{1,1}, '_1')
% pass parameters as vector
if iscell(objFctSettings)
testval = objFctHandle(objFctSettings{:}, testmem');
else
testval = objFctHandle(objFctSettings,    testmem');
end
else
% pass parameters as structure objFctParams
if iscell(objFctSettings)
testval = objFctHandle(objFctSettings{:}, objFctParams);
else
testval = objFctHandle(objFctSettings,    objFctParams);
end
end
nrOfMasterEvaluations = nrOfMasterEvaluations + 1;
nfeval.overall = nfeval.overall + 1;

% check returned value
if isempty(testval)
error('Objective function returned empty evaluation value!');
elseif ~isscalar(testval)
error('Objective function returned vector as evaluation value!');
end
end

else
% compute distance to randomly chosen vector for testing
if isempty(optParVector)
optParVector = XVmin + rand(1, length(XVmin)) .* (XVmax - XVmin);
end
testval = sqrt(sum(testmem - optParVector).^2);
pause(0.2);

end % if ~isempty(objFctHandle)

if ~isempty(testval)
% save tested vector and resulting cost value
if saveHistory
allmem(:,end+1) = testmem';
allval  (end+1) = testval;
end
else
testval = NaN;
end
end % if ~isempty(allmemIndex)

% write non-quantized value into population
testmem = testmem2;

```