differentialevolution.m

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

		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);
maxNOfTests = min(1000, 10*length(memIndex));
nOfRecomputations = 0;
for k=1:maxNOfTests
	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 + nOfRecomputations/maxNOfTests*randStdDevAdd;
		nOfRecomputations = nOfRecomputations + 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, nOfPossibleMembers, 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.local)];
if nfeval.slave > 0
	str = [str, sprintf('Number of slave evaluations:    %d\n',   nfeval.slave)];
end
percentage = round(nfeval.local / nOfPossibleMembers * 100);
if percentage > 0
	str = [str, sprintf('Percentage of members computed: %d %%\n', percentage)];
end
if nfeval.local - nfeval.slave > 0
	str = [str, sprintf('Mean function evaluation time:  %s\n', ...
		formattime(elapsedTime / (nfeval.local-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
	disp(sprintf('Results are saved in file %s.', 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