Code covered by the BSD License  

Highlights from
Differential Evolution

Differential Evolution

by

 

03 Feb 2008 (Updated )

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.
%    emailParams      - Struct with fields serveraddress, fromaddress,
%                       toaddress, and, if needed, username and password.
%                       Parameters are used for sending E-mail
%                       notifications.
%    optimInfo        - Info about current optimization task. Fields 'title'
%                       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\
%    Temp\<UserName>@<HostName>\MATLAB).
%
%    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
%    Last modified 09.05.2013
%
%    See also DIFFERENTIALEVOLUTIONSLAVE, DISPLAYOPTIMIZATIONHISTORY,
%    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;
    objFctHandle            = @rosenbrocksaddle;
    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
      NPAdd = min(round(0.4*NP), NP-istart+1);
      if NPAdd > 0
        memIndex = istart:istart+NPAdd-1;
        istart   = istart+NPAdd+1;
        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;
    randStdDevAdd = 0.9;
  case 2
    % only use random numbers
    baseMem = XVmin;
    randStdDev1   = 1;
    randStdDevAdd = 0;
  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
    randStdDev = randStdDev1 + nrOfRecomputations/maxNrOfTests*randStdDevAdd;
    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, ...
                                forceFileUpload)

persistent resultFileName nextFileUploadTime

if nargin == 0
  nextFileUploadTime = NaN;
  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)
    load(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
  if isnan(nextFileUploadTime)
    nextFileUploadTime = mbtime + fileUploadPeriod;
  end

  if forceFileUpload || mbtime - nextFileUploadTime >= 0
    nextFileUploadTime = mbtime + fileUploadPeriod;
    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
persistent maxNameLength hostname username

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;
  username = getusername;

  % 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

  switch [username '@' hostname]
    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
    [x, fs, bits] = wavread('applause.wav');
    soundsc(x, fs, bits);
    pause(length(x) / fs + 1);
    lastSoundTime = mbtime;
  end
else
  sendEmailThisTime = 0;
end

if ~sendEmail
  return
end

% send E-mail notification
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
  % read from slave
  
  % 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)
          load(slaveFileName, 'testval');
          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;

Contact us