differentialevolution.m

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

	error('maxiter must be greater than zero.');
end
if DEParams.displayResults && ~DEParams.saveHistory
	disp(sprintf('Warning: Optimization history can not be displayed if not saved.\n'));
	DEParams.displayResults = 0;
end
DEParams.refreshiter = floor(DEParams.refreshiter);

% get parameters
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
	disp(sprintf('Starting optimization at %s.', datestr(clock, 'mmm dd, HH:MM')));
end

% display number of possible parameter vectors
if isfinite(nOfPossibleMembers)
	disp(sprintf('Number of possible parameter vectors: %g', nOfPossibleMembers));
else
	disp(sprintf('Infinite number of possible parameter vectors.'));
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.local    = 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 nOfPossibleMembers <= 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
		disp(sprintf('Re-initializing population in iteration %d.', iterationNr));
		if valstddev < DEParams.minvalstddev
			disp(sprintf('Evaluation value std. dev. below threshold (%g < %g).', ...
				valstddev, DEParams.minvalstddev));
		elseif paramstddev < DEParams.minparamstddev
			disp(sprintf('Mean parameter std. dev. below threshold (%g < %g).', ...
				paramstddev, DEParams.minparamstddev));
		elseif noChangeCounter >= DEParams.nochangeiter
			disp(sprintf('Population did not change in last %d iterations.', ...
				DEParams.nochangeiter));
		elseif nofevalCounter > DEParams.nofevaliter
			disp(sprintf('No function evaluations in the last %d iterations.', ...
				DEParams.nofevaliter));
		end
		disp(' ');

		if minMaxSign > 0
			[val, index] = sort( val); 
		else
			[val, index] = sort(-val); % not using 'descend' for downwards 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