differentialevolution.m

This contribution provides functions for finding an optimum parameter set using the evolutionary alg

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.
%		nOfIterations    - 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. 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
%		  'nOfCoefficients', [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 on
%		one processor/computer, function DIFFERENTIALEVOLUTIONSLAVE on one or
%		more other processors/computers. 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
%		evolutional algorithm was taken from function devec3.m.
%
%		Markus Buehren
%		Last modified 26.08.2008 
%
%		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)
		disp(sprintf('Error: Lower bound (%g) of parameter %s is larger than upper bound (%g).', ...
			XVmin(parNr), getparametername__(parNr, 1), XVmax(parNr)));
		errorFound = true;
	end
	if parGridVector(parNr) < 0
		disp(sprintf('Error: Negative quantization values are not allowed (parameter %s).', ...
			getparametername__(parNr, 1)));
		errorFound = true;
	end		
	minQuantizationUser = 1e-12;
	if parGridVector(parNr) > 0 && parGridVector(parNr) < minQuantizationUser
		disp(sprintf('Error: Minimum quantization step size is %g (parameter %s).', ...
			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)
	nOfPossibleMembers = prod(floor((diff(parameterBounds, 1, 2) + 0.5*parGridVector) ./ parGridVector) + 1);
else
	nOfPossibleMembers = inf;
end

% check parameters
DEParams.NP = min(DEParams.NP, nOfPossibleMembers);
if (DEParams.maxiter <= 0)