1function [xmin, ... % minimum search point of last iteration 2 fmin, ... % function value of xmin 3 counteval, ... % number of function evaluations done 4 stopflag, ... % stop criterion reached 5 out, ... % struct with various histories and solutions 6 bestever ... % struct containing overall best solution (for convenience) 7 ] = cmaes( ... 8 fitfun, ... % name of objective/fitness function 9 xstart, ... % objective variables initial point, determines N 10 insigma, ... % initial coordinate wise standard deviation(s) 11 inopts, ... % options struct, see defopts below 12 varargin ) % arguments passed to objective function 13 % cmaes.m, Version 3.56.beta, last change: February, 2012 14 % CMAES implements an Evolution Strategy with Covariance Matrix 15 % Adaptation (CMA-ES) for nonlinear function minimization. For 16 % introductory comments and copyright (GPL) see end of file (type 17 % 'type cmaes'). cmaes.m runs with MATLAB (Windows, Linux) and, 18 % without data logging and plotting, it should run under Octave 19 % (Linux, package octave-forge is needed). 20 % 21 % OPTS = CMAES returns default options. 22 % OPTS = CMAES('defaults') returns default options quietly. 23 % OPTS = CMAES('displayoptions') displays options. 24 % OPTS = CMAES('defaults', OPTS) supplements options OPTS with default 25 % options. 26 % 27 % XMIN = CMAES(FUN, X0, SIGMA[, OPTS]) locates the minimum XMIN of 28 % function FUN starting from column vector X0 with the initial 29 % coordinate wise search standard deviation SIGMA. 30 % 31 % Input arguments: 32 % 33 % FUN is a string function name like 'myfun'. FUN takes as argument a 34 % column vector of size of X0 and returns a scalar. An easy way to 35 % implement a hard non-linear constraint is to return NaN. Then, 36 % this function evaluation is not counted and a newly sampled 37 % point is tried immediately. 38 % 39 % X0 is a column vector, or a matrix, or a string. If X0 is a matrix, 40 % mean(X0, 2) is taken as initial point. If X0 is a string like 41 % '2*rand(10,1)-1', the string is evaluated first. 42 % 43 % SIGMA is a scalar, or a column vector of size(X0,1), or a string 44 % that can be evaluated into one of these. SIGMA determines the 45 % initial coordinate wise standard deviations for the search. 46 % Setting SIGMA one third of the initial search region is 47 % appropriate, e.g., the initial point in [0, 6]^10 and SIGMA=2 48 % means cmaes('myfun', 3*rand(10,1), 2). If SIGMA is missing and 49 % size(X0,2) > 1, SIGMA is set to sqrt(var(X0')'). That is, X0 is 50 % used as a sample for estimating initial mean and variance of the 51 % search distribution. 52 % 53 % OPTS (an optional argument) is a struct holding additional input 54 % options. Valid field names and a short documentation can be 55 % discovered by looking at the default options (type 'cmaes' 56 % without arguments, see above). Empty or missing fields in OPTS 57 % invoke the default value, i.e. OPTS needs not to have all valid 58 % field names. Capitalization does not matter and unambiguous 59 % abbreviations can be used for the field names. If a string is 60 % given where a numerical value is needed, the string is evaluated 61 % by eval, where 'N' expands to the problem dimension 62 % (==size(X0,1)) and 'popsize' to the population size. 63 % 64 % [XMIN, FMIN, COUNTEVAL, STOPFLAG, OUT, BESTEVER] = ... 65 % CMAES(FITFUN, X0, SIGMA) 66 % returns the best (minimal) point XMIN (found in the last 67 % generation); function value FMIN of XMIN; the number of needed 68 % function evaluations COUNTEVAL; a STOPFLAG value as cell array, 69 % where possible entries are 'fitness', 'tolx', 'tolupx', 'tolfun', 70 % 'maxfunevals', 'maxiter', 'stoptoresume', 'manual', 71 % 'warnconditioncov', 'warnnoeffectcoord', 'warnnoeffectaxis', 72 % 'warnequalfunvals', 'warnequalfunvalhist', 'bug' (use 73 % e.g. any(strcmp(STOPFLAG, 'tolx')) or findstr(strcat(STOPFLAG, 74 % 'tolx')) for further processing); a record struct OUT with some 75 % more output, where the struct SOLUTIONS.BESTEVER contains the overall 76 % best evaluated point X with function value F evaluated at evaluation 77 % count EVALS. The last output argument BESTEVER equals 78 % OUT.SOLUTIONS.BESTEVER. Moreover a history of solutions and 79 % parameters is written to files according to the Log-options. 80 % 81 % A regular manual stop can be achieved via the file signals.par. The 82 % program is terminated if the first two non-white sequences in any 83 % line of this file are 'stop' and the value of the LogFilenamePrefix 84 % option (by default 'outcmaes'). Also a run can be skipped. 85 % Given, for example, 'skip outcmaes run 2', skips the second run 86 % if option Restarts is at least 2, and another run will be started. 87 % 88 % To run the code completely silently set Disp, Save, and Log options 89 % to 0. With OPTS.LogModulo > 0 (1 by default) the most important 90 % data are written to ASCII files permitting to investigate the 91 % results (e.g. plot with function plotcmaesdat) even while CMAES is 92 % still running (which can be quite useful on expensive objective 93 % functions). When OPTS.SaveVariables==1 (default) everything is saved 94 % in file OPTS.SaveFilename (default 'variablescmaes.mat') allowing to 95 % resume the search afterwards by using the resume option. 96 % 97 % To find the best ever evaluated point load the variables typing 98 % "es=load('variablescmaes')" and investigate the variable 99 % es.out.solutions.bestever. 100 % 101 % In case of a noisy objective function (uncertainties) set 102 % OPTS.Noise.on = 1. This option interferes presumably with some 103 % termination criteria, because the step-size sigma will presumably 104 % not converge to zero anymore. If CMAES was provided with a 105 % fifth argument (P1 in the below example, which is passed to the 106 % objective function FUN), this argument is multiplied with the 107 % factor given in option Noise.alphaevals, each time the detected 108 % noise exceeds a threshold. This argument can be used within 109 % FUN, for example, as averaging number to reduce the noise level. 110 % 111 % OPTS.DiagonalOnly > 1 defines the number of initial iterations, 112 % where the covariance matrix remains diagonal and the algorithm has 113 % internally linear time complexity. OPTS.DiagonalOnly = 1 means 114 % keeping the covariance matrix always diagonal and this setting 115 % also exhibits linear space complexity. This can be particularly 116 % useful for dimension > 100. The default is OPTS.DiagonalOnly = 0. 117 % 118 % OPTS.CMA.active = 1 turns on "active CMA" with a negative update 119 % of the covariance matrix and checks for positive definiteness. 120 % OPTS.CMA.active = 2 does not check for pos. def. and is numerically 121 % faster. Active CMA usually speeds up the adaptation and might 122 % become a default in near future. 123 % 124 % The primary strategy parameter to play with is OPTS.PopSize, which 125 % can be increased from its default value. Increasing the population 126 % size (by default linked to increasing parent number OPTS.ParentNumber) 127 % improves global search properties in exchange to speed. Speed 128 % decreases, as a rule, at most linearely with increasing population 129 % size. It is advisable to begin with the default small population 130 % size. The options Restarts and IncPopSize can be used for an 131 % automated multistart where the population size is increased by the 132 % factor IncPopSize (two by default) before each restart. X0 (given as 133 % string) is reevaluated for each restart. Stopping options 134 % StopFunEvals, StopIter, MaxFunEvals, and Fitness terminate the 135 % program, all others including MaxIter invoke another restart, where 136 % the iteration counter is reset to zero. 137 % 138 % Examples: 139 % 140 % XMIN = cmaes('myfun', 5*ones(10,1), 1.5); starts the search at 141 % 10D-point 5 and initially searches mainly between 5-3 and 5+3 142 % (+- two standard deviations), but this is not a strict bound. 143 % 'myfun' is a name of a function that returns a scalar from a 10D 144 % column vector. 145 % 146 % opts.LBounds = 0; opts.UBounds = 10; 147 % X=cmaes('myfun', 10*rand(10,1), 5, opts); 148 % search within lower bound of 0 and upper bound of 10. Bounds can 149 % also be given as column vectors. If the optimum is not located 150 % on the boundary, use rather a penalty approach to handle bounds. 151 % 152 % opts=cmaes; opts.StopFitness=1e-10; 153 % X=cmaes('myfun', rand(5,1), 0.5, opts); stops the search, if 154 % the function value is smaller than 1e-10. 155 % 156 % [X, F, E, STOP, OUT] = cmaes('myfun2', 'rand(5,1)', 1, [], P1, P2); 157 % passes two additional parameters to the function MYFUN2. 158 % 159 160% Copyright (C) 2001-2012 Nikolaus Hansen, 161% Copyright (C) 2012-2017 Dynare Team 162% 163% This file is part of Dynare. 164% 165% Dynare is free software: you can redistribute it and/or modify 166% it under the terms of the GNU General Public License as published by 167% the Free Software Foundation, either version 3 of the License, or 168% (at your option) any later version. 169% 170% Dynare is distributed in the hope that it will be useful, 171% but WITHOUT ANY WARRANTY; without even the implied warranty of 172% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 173% GNU General Public License for more details. 174% 175% You should have received a copy of the GNU General Public License 176% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 177 178 179cmaVersion = '3.60.beta'; 180 181% ----------- Set Defaults for Input Parameters and Options ------------- 182% These defaults may be edited for convenience 183% Input Defaults (obsolete, these are obligatory now) 184definput.fitfun = 'felli'; % frosen; fcigar; see end of file for more 185definput.xstart = rand(10,1); % 0.50*ones(10,1); 186definput.sigma = 0.3; 187 188% Options defaults: Stopping criteria % (value of stop flag) 189defopts.StopFitness = '-Inf % stop if f(xmin) < stopfitness, minimization'; 190defopts.MaxFunEvals = 'Inf % maximal number of fevals'; 191defopts.MaxIter = '1e3*(N+5)^2/sqrt(popsize) % maximal number of iterations'; 192defopts.StopFunEvals = 'Inf % stop after resp. evaluation, possibly resume later'; 193defopts.StopIter = 'Inf % stop after resp. iteration, possibly resume later'; 194defopts.TolX = '1e-9*max(insigma) % stop if x-change smaller TolX'; 195defopts.TolUpX = '1e3*max(insigma) % stop if x-changes larger TolUpX'; 196defopts.TolFun = '1e-10 % stop if fun-changes smaller TolFun'; 197defopts.TolHistFun = '1e-11 % stop if back fun-changes smaller TolHistFun'; 198defopts.StopOnStagnation = 'on % stop when fitness stagnates for a long time'; 199defopts.StopOnWarnings = 'yes % ''no''==''off''==0, ''on''==''yes''==1 '; 200defopts.StopOnEqualFunctionValues = '2 + N/3 % number of iterations'; 201 202% Options defaults: Other 203defopts.DiffMaxChange = 'Inf % maximal variable change(s), can be Nx1-vector'; 204defopts.DiffMinChange = '0 % minimal variable change(s), can be Nx1-vector'; 205defopts.WarnOnEqualFunctionValues = ... 206 'yes % ''no''==''off''==0, ''on''==''yes''==1 '; 207defopts.LBounds = '-Inf % lower bounds, scalar or Nx1-vector'; 208defopts.UBounds = 'Inf % upper bounds, scalar or Nx1-vector'; 209defopts.EvalParallel = 'no % objective function FUN accepts NxM matrix, with M>1?'; 210defopts.EvalInitialX = 'yes % evaluation of initial solution'; 211defopts.Restarts = '0 % number of restarts '; 212defopts.IncPopSize = '2 % multiplier for population size before each restart'; 213 214defopts.PopSize = '(4 + floor(3*log(N))) % population size, lambda'; 215defopts.ParentNumber = 'floor(popsize/2) % AKA mu, popsize equals lambda'; 216defopts.RecombinationWeights = 'superlinear decrease % or linear, or equal'; 217defopts.DiagonalOnly = '0*(1+100*N/sqrt(popsize))+(N>=1000) % C is diagonal for given iterations, 1==always'; 218defopts.Noise.on = '0 % uncertainty handling is off by default'; 219defopts.Noise.reevals = '1*ceil(0.05*lambda) % nb. of re-evaluated for uncertainty measurement'; 220defopts.Noise.theta = '0.5 % threshold to invoke uncertainty treatment'; % smaller: more likely to diverge 221defopts.Noise.cum = '0.3 % cumulation constant for uncertainty'; 222defopts.Noise.cutoff = '2*lambda/3 % rank change cutoff for summation'; 223defopts.Noise.alphasigma = '1+2/(N+10) % factor for increasing sigma'; % smaller: slower adaptation 224defopts.Noise.epsilon = '1e-7 % additional relative perturbation before reevaluation'; 225defopts.Noise.minmaxevals = '[1 inf] % min and max value of 2nd arg to fitfun, start value is 5th arg to cmaes'; 226defopts.Noise.alphaevals = '1+2/(N+10) % factor for increasing 2nd arg to fitfun'; 227defopts.Noise.callback = '[] % callback function when uncertainty threshold is exceeded'; 228% defopts.TPA = 0; 229defopts.CMA.cs = '(mueff+2)/(N+mueff+3) % cumulation constant for step-size'; 230%qqq cs = (mueff^0.5)/(N^0.5+mueff^0.5) % the short time horizon version 231defopts.CMA.damps = '1 + 2*max(0,sqrt((mueff-1)/(N+1))-1) + cs % damping for step-size'; 232% defopts.CMA.ccum = '4/(N+4) % cumulation constant for covariance matrix'; 233defopts.CMA.ccum = '(4 + mueff/N) / (N+4 + 2*mueff/N) % cumulation constant for pc'; 234defopts.CMA.ccov1 = '2 / ((N+1.3)^2+mueff) % learning rate for rank-one update'; 235defopts.CMA.ccovmu = '2 * (mueff-2+1/mueff) / ((N+2)^2+mueff) % learning rate for rank-mu update'; 236defopts.CMA.on = 'yes'; 237defopts.CMA.active = '0 % active CMA 1: neg. updates with pos. def. check, 2: neg. updates'; 238 239flg_future_setting = 0; % testing for possible future variant(s) 240if flg_future_setting 241 disp('in the future') 242 243 % damps setting from Brockhoff et al 2010 244 % this damps diverges with popsize 400: 245 % cmaeshtml('benchmarkszero', ones(20,1)*2, 5, o, 15); 246 defopts.CMA.damps = '2*mueff/lambda + 0.3 + cs % damping for step-size'; % cs: for large mueff 247 % how about: 248 % defopts.CMA.damps = '2*mueff/lambda + 0.3 + 2*max(0,sqrt((mueff-1)/(N+1))-1) + cs % damping for step-size'; 249 250 % ccum adjusted for large mueff, better on schefelmult? 251 % TODO: this should also depend on diagonal option!? 252 defopts.CMA.ccum = '(4 + mueff/N) / (N+4 + 2*mueff/N) % cumulation constant for pc'; 253 254 defopts.CMA.active = '1 % active CMA 1: neg. updates with pos. def. check, 2: neg. updates'; 255end 256 257defopts.Resume = 'no % resume former run from SaveFile'; 258defopts.Science = 'on % off==do some additional (minor) problem capturing, NOT IN USE'; 259defopts.ReadSignals = 'on % from file signals.par for termination, yet a stumb'; 260defopts.Seed = 'sum(100*clock) % evaluated if it is a string'; 261defopts.DispFinal = 'on % display messages like initial and final message'; 262defopts.DispModulo = '100 % [0:Inf], disp messages after every i-th iteration'; 263defopts.SaveVariables = 'on % [on|final|off][-v6] save variables to .mat file'; 264defopts.SaveFilename = 'variablescmaes.mat % save all variables, see SaveVariables'; 265defopts.LogModulo = '1 % [0:Inf] if >1 record data less frequently after gen=100'; 266defopts.LogTime = '25 % [0:100] max. percentage of time for recording data'; 267defopts.LogFilenamePrefix = 'outcmaes % files for output data'; 268defopts.LogPlot = 'off % plot while running using output data files'; 269 270%qqqkkk 271%defopts.varopt1 = ''; % 'for temporary and hacking purposes'; 272%defopts.varopt2 = ''; % 'for temporary and hacking purposes'; 273defopts.UserData = 'for saving data/comments associated with the run'; 274defopts.UserDat2 = ''; 'for saving data/comments associated with the run'; 275 276% ---------------------- Handling Input Parameters ---------------------- 277 278if nargin < 1 || isequal(fitfun, 'defaults') % pass default options 279 if nargin < 1 280 disp('Default options returned (type "help cmaes" for help).'); 281 end 282 xmin = defopts; 283 if nargin > 1 % supplement second argument with default options 284 xmin = getoptions(xstart, defopts); 285 end 286 return 287end 288 289if isequal(fitfun, 'displayoptions') 290 names = fieldnames(defopts); 291 for name = names' 292 disp([name{:} repmat(' ', 1, 20-length(name{:})) ': ''' defopts.(name{:}) '''']); 293 end 294 return 295end 296 297input.fitfun = fitfun; % record used input 298if isempty(fitfun) 299 % fitfun = definput.fitfun; 300 % warning(['Objective function not determined, ''' fitfun ''' used']); 301 error(['Objective function not determined']); 302end 303if ~ischar(fitfun) 304 error('first argument FUN must be a string'); 305end 306 307 308if nargin < 2 309 xstart = []; 310end 311 312input.xstart = xstart; 313if isempty(xstart) 314 % xstart = definput.xstart; % objective variables initial point 315 % warning('Initial search point, and problem dimension, not determined'); 316 error('Initial search point, and problem dimension, not determined'); 317end 318 319if nargin < 3 320 insigma = []; 321end 322if isa(insigma, 'struct') 323 error(['Third argument SIGMA must be (or eval to) a scalar '... 324 'or a column vector of size(X0,1)']); 325end 326input.sigma = insigma; 327if isempty(insigma) 328 if all(size(myeval(xstart)) > 1) 329 insigma = std(xstart, 0, 2); 330 if any(insigma == 0) 331 error(['Initial search volume is zero, choose SIGMA or X0 appropriate']); 332 end 333 else 334 % will be captured later 335 % error(['Initial step sizes (SIGMA) not determined']); 336 end 337end 338 339% Compose options opts 340if nargin < 4 || isempty(inopts) % no input options available 341 inopts = []; 342 opts = defopts; 343else 344 opts = getoptions(inopts, defopts); 345end 346i = strfind(opts.SaveFilename, ' '); % remove everything after white space 347if ~isempty(i) 348 opts.SaveFilename = opts.SaveFilename(1:i(1)-1); 349end 350 351%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 352counteval = 0; countevalNaN = 0; 353irun = 0; 354while irun <= myeval(opts.Restarts) % for-loop does not work with resume 355 irun = irun + 1; 356 357 % ------------------------ Initialization ------------------------------- 358 359 % Handle resuming of old run 360 flgresume = myevalbool(opts.Resume); 361 xmean = myeval(xstart); 362 if all(size(xmean) > 1) 363 xmean = mean(xmean, 2); % in case if xstart is a population 364 elseif size(xmean, 2) > 1 365 xmean = xmean'; 366 end 367 if ~flgresume % not resuming a former run 368 % Assign settings from input parameters and options for myeval... 369 N = size(xmean, 1); numberofvariables = N; 370 lambda0 = floor(myeval(opts.PopSize) * myeval(opts.IncPopSize)^(irun-1)); 371 % lambda0 = floor(myeval(opts.PopSize) * 3^floor((irun-1)/2)); 372 popsize = lambda0; 373 lambda = lambda0; 374 insigma = myeval(insigma); 375 if all(size(insigma) == [N 2]) 376 insigma = 0.5 * (insigma(:,2) - insigma(:,1)); 377 end 378 else % flgresume is true, do resume former run 379 tmp = whos('-file', opts.SaveFilename); 380 for i = 1:length(tmp) 381 if strcmp(tmp(i).name, 'localopts') 382 error('Saved variables include variable "localopts", please remove'); 383 end 384 end 385 local.opts = opts; % keep stopping and display options 386 local.varargin = varargin; 387 load(opts.SaveFilename); 388 varargin = local.varargin; 389 flgresume = 1; 390 391 % Overwrite old stopping and display options 392 opts.StopFitness = local.opts.StopFitness; 393 %%opts.MaxFunEvals = local.opts.MaxFunEvals; 394 %%opts.MaxIter = local.opts.MaxIter; 395 opts.StopFunEvals = local.opts.StopFunEvals; 396 opts.StopIter = local.opts.StopIter; 397 opts.TolX = local.opts.TolX; 398 opts.TolUpX = local.opts.TolUpX; 399 opts.TolFun = local.opts.TolFun; 400 opts.TolHistFun = local.opts.TolHistFun; 401 opts.StopOnStagnation = local.opts.StopOnStagnation; 402 opts.StopOnWarnings = local.opts.StopOnWarnings; 403 opts.ReadSignals = local.opts.ReadSignals; 404 opts.DispFinal = local.opts.DispFinal; 405 opts.LogPlot = local.opts.LogPlot; 406 opts.DispModulo = local.opts.DispModulo; 407 opts.SaveVariables = local.opts.SaveVariables; 408 opts.LogModulo = local.opts.LogModulo; 409 opts.LogTime = local.opts.LogTime; 410 clear local; % otherwise local would be overwritten during load 411 end 412 413 %-------------------------------------------------------------- 414 % Evaluate options 415 stopFitness = myeval(opts.StopFitness); 416 stopMaxFunEvals = myeval(opts.MaxFunEvals); 417 stopMaxIter = myeval(opts.MaxIter); 418 stopFunEvals = myeval(opts.StopFunEvals); 419 stopIter = myeval(opts.StopIter); 420 stopTolX = myeval(opts.TolX); 421 stopTolUpX = myeval(opts.TolUpX); 422 stopTolFun = myeval(opts.TolFun); 423 stopTolHistFun = myeval(opts.TolHistFun); 424 stopOnStagnation = myevalbool(opts.StopOnStagnation); 425 stopOnWarnings = myevalbool(opts.StopOnWarnings); 426 flgreadsignals = myevalbool(opts.ReadSignals); 427 flgWarnOnEqualFunctionValues = myevalbool(opts.WarnOnEqualFunctionValues); 428 flgEvalParallel = myevalbool(opts.EvalParallel); 429 stopOnEqualFunctionValues = myeval(opts.StopOnEqualFunctionValues); 430 arrEqualFunvals = zeros(1,10+N); 431 flgDiagonalOnly = myeval(opts.DiagonalOnly); 432 flgActiveCMA = myeval(opts.CMA.active); 433 noiseHandling = myevalbool(opts.Noise.on); 434 noiseMinMaxEvals = myeval(opts.Noise.minmaxevals); 435 noiseAlphaEvals = myeval(opts.Noise.alphaevals); 436 noiseCallback = myeval(opts.Noise.callback); 437 flgdisplay = myevalbool(opts.DispFinal); 438 flgplotting = myevalbool(opts.LogPlot); 439 verbosemodulo = myeval(opts.DispModulo); 440 flgscience = myevalbool(opts.Science); 441 flgsaving = []; 442 strsaving = []; 443 if strfind(opts.SaveVariables, '-v6') 444 i = strfind(opts.SaveVariables, '%'); 445 if isempty(i) || i == 0 || strfind(opts.SaveVariables, '-v6') < i 446 strsaving = '-v6'; 447 flgsaving = 1; 448 flgsavingfinal = 1; 449 end 450 end 451 if strncmp('final', opts.SaveVariables, 5) 452 flgsaving = 0; 453 flgsavingfinal = 1; 454 end 455 if isempty(flgsaving) 456 flgsaving = myevalbool(opts.SaveVariables); 457 flgsavingfinal = flgsaving; 458 end 459 savemodulo = myeval(opts.LogModulo); 460 savetime = myeval(opts.LogTime); 461 462 i = strfind(opts.LogFilenamePrefix, ' '); % remove everything after white space 463 if ~isempty(i) 464 opts.LogFilenamePrefix = opts.LogFilenamePrefix(1:i(1)-1); 465 end 466 467 % TODO here silent option? set disp, save and log options to 0 468 469 %-------------------------------------------------------------- 470 471 if (isfinite(stopFunEvals) || isfinite(stopIter)) && ~flgsaving 472 warning('To resume later the saving option needs to be set'); 473 end 474 475 476 % Do more checking and initialization 477 if flgresume % resume is on 478 time.t0 = clock; 479 if flgdisplay 480 disp([' resumed from ' opts.SaveFilename ]); 481 end 482 if counteval >= stopMaxFunEvals 483 error(['MaxFunEvals exceeded, use StopFunEvals as stopping ' ... 484 'criterion before resume']); 485 end 486 if countiter >= stopMaxIter 487 error(['MaxIter exceeded, use StopIter as stopping criterion ' ... 488 'before resume']); 489 end 490 491 else % flgresume 492 % xmean = mean(myeval(xstart), 2); % evaluate xstart again, because of irun 493 maxdx = myeval(opts.DiffMaxChange); % maximal sensible variable change 494 mindx = myeval(opts.DiffMinChange); % minimal sensible variable change 495 % can both also be defined as Nx1 vectors 496 lbounds = myeval(opts.LBounds); 497 ubounds = myeval(opts.UBounds); 498 if length(lbounds) == 1 499 lbounds = repmat(lbounds, N, 1); 500 end 501 if length(ubounds) == 1 502 ubounds = repmat(ubounds, N, 1); 503 end 504 if isempty(insigma) % last chance to set insigma 505 if all(lbounds > -Inf) && all(ubounds < Inf) 506 if any(lbounds>=ubounds) 507 error('upper bound must be greater than lower bound'); 508 end 509 insigma = 0.3*(ubounds-lbounds); 510 stopTolX = myeval(opts.TolX); % reevaluate these 511 stopTolUpX = myeval(opts.TolUpX); 512 else 513 error(['Initial step sizes (SIGMA) not determined']); 514 end 515 end 516 517 % Check all vector sizes 518 if size(xmean, 2) > 1 || size(xmean,1) ~= N 519 error(['intial search point should be a column vector of size ' ... 520 num2str(N)]); 521 elseif ~(all(size(insigma) == [1 1]) || all(size(insigma) == [N 1])) 522 error(['input parameter SIGMA should be (or eval to) a scalar '... 523 'or a column vector of size ' num2str(N)] ); 524 elseif size(stopTolX, 2) > 1 || ~ismember(size(stopTolX, 1), [1 N]) 525 error(['option TolX should be (or eval to) a scalar '... 526 'or a column vector of size ' num2str(N)] ); 527 elseif size(stopTolUpX, 2) > 1 || ~ismember(size(stopTolUpX, 1), [1 N]) 528 error(['option TolUpX should be (or eval to) a scalar '... 529 'or a column vector of size ' num2str(N)] ); 530 elseif size(maxdx, 2) > 1 || ~ismember(size(maxdx, 1), [1 N]) 531 error(['option DiffMaxChange should be (or eval to) a scalar '... 532 'or a column vector of size ' num2str(N)] ); 533 elseif size(mindx, 2) > 1 || ~ismember(size(mindx, 1), [1 N]) 534 error(['option DiffMinChange should be (or eval to) a scalar '... 535 'or a column vector of size ' num2str(N)] ); 536 elseif size(lbounds, 2) > 1 || ~ismember(size(lbounds, 1), [1 N]) 537 error(['option lbounds should be (or eval to) a scalar '... 538 'or a column vector of size ' num2str(N)] ); 539 elseif size(ubounds, 2) > 1 || ~ismember(size(ubounds, 1), [1 N]) 540 error(['option ubounds should be (or eval to) a scalar '... 541 'or a column vector of size ' num2str(N)] ); 542 end 543 544 % Initialize dynamic internal state parameters 545 if any(insigma <= 0) 546 error(['Initial search volume (SIGMA) must be greater than zero']); 547 end 548 if max(insigma)/min(insigma) > 1e6 549 error(['Initial search volume (SIGMA) badly conditioned']); 550 end 551 sigma = max(insigma); % overall standard deviation 552 pc = zeros(N,1); ps = zeros(N,1); % evolution paths for C and sigma 553 554 if length(insigma) == 1 555 insigma = insigma * ones(N,1) ; 556 end 557 diagD = insigma/max(insigma); % diagonal matrix D defines the scaling 558 diagC = diagD.^2; 559 if flgDiagonalOnly ~= 1 % use at some point full covariance matrix 560 B = eye(N,N); % B defines the coordinate system 561 BD = B.*repmat(diagD',N,1); % B*D for speed up only 562 C = diag(diagC); % covariance matrix == BD*(BD)' 563 end 564 if flgDiagonalOnly 565 B = 1; 566 end 567 568 fitness.hist=NaN*ones(1,10+ceil(3*10*N/lambda)); % history of fitness values 569 fitness.histsel=NaN*ones(1,10+ceil(3*10*N/lambda)); % history of fitness values 570 fitness.histbest=[]; % history of fitness values 571 fitness.histmedian=[]; % history of fitness values 572 573 % Initialize boundary handling 574 bnd.isactive = any(lbounds > -Inf) || any(ubounds < Inf); 575 if bnd.isactive 576 if any(lbounds>ubounds) 577 error('lower bound found to be greater than upper bound'); 578 end 579 [xmean, ti] = xintobounds(xmean, lbounds, ubounds); % just in case 580 if any(ti) 581 warning('Initial point was out of bounds, corrected'); 582 end 583 bnd.weights = zeros(N,1); % weights for bound penalty 584 % scaling is better in axis-parallel case, worse in rotated 585 bnd.flgscale = 0; % scaling will be omitted if zero 586 if bnd.flgscale ~= 0 587 bnd.scale = diagC/mean(diagC); 588 else 589 bnd.scale = ones(N,1); 590 end 591 592 idx = (lbounds > -Inf) | (ubounds < Inf); 593 if length(idx) == 1 594 idx = idx * ones(N,1); 595 end 596 bnd.isbounded = zeros(N,1); 597 bnd.isbounded(find(idx)) = 1; 598 maxdx = min(maxdx, (ubounds - lbounds)/2); 599 if any(sigma*sqrt(diagC) > maxdx) 600 fac = min(maxdx ./ sqrt(diagC))/sigma; 601 sigma = min(maxdx ./ sqrt(diagC)); 602 warning(['Initial SIGMA multiplied by the factor ' num2str(fac) ... 603 ', because it was larger than half' ... 604 ' of one of the boundary intervals']); 605 end 606 idx = (lbounds > -Inf) & (ubounds < Inf); 607 dd = diagC; 608 if any(5*sigma*sqrt(dd(idx)) < ubounds(idx) - lbounds(idx)) 609 warning(['Initial SIGMA is, in at least one coordinate, ' ... 610 'much smaller than the '... 611 'given boundary intervals. For reasonable ' ... 612 'global search performance SIGMA should be ' ... 613 'between 0.2 and 0.5 of the bounded interval in ' ... 614 'each coordinate. If all coordinates have ' ... 615 'lower and upper bounds SIGMA can be empty']); 616 end 617 bnd.dfithist = 1; % delta fit for setting weights 618 bnd.aridxpoints = []; % remember complete outside points 619 bnd.arfitness = []; % and their fitness 620 bnd.validfitval = 0; 621 bnd.iniphase = 1; 622 end 623 624 % ooo initial feval, for output only 625 if irun == 1 626 out.solutions.bestever.x = xmean; 627 out.solutions.bestever.f = Inf; % for simpler comparison below 628 out.solutions.bestever.evals = counteval; 629 bestever = out.solutions.bestever; 630 end 631 if myevalbool(opts.EvalInitialX) 632 fitness.hist(1)=feval(fitfun, xmean, varargin{:}); 633 fitness.histsel(1)=fitness.hist(1); 634 counteval = counteval + 1; 635 if fitness.hist(1) < out.solutions.bestever.f 636 out.solutions.bestever.x = xmean; 637 out.solutions.bestever.f = fitness.hist(1); 638 out.solutions.bestever.evals = counteval; 639 bestever = out.solutions.bestever; 640 end 641 else 642 fitness.hist(1)=NaN; 643 fitness.histsel(1)=NaN; 644 end 645 646 % initialize random number generator 647% $$$ if ischar(opts.Seed) 648% $$$ randn('state', eval(opts.Seed)); % random number generator state 649% $$$ else 650% $$$ randn('state', opts.Seed); 651% $$$ end 652%qqq 653% load(opts.SaveFilename, 'startseed'); 654% randn('state', startseed); 655% disp(['SEED RELOADED FROM ' opts.SaveFilename]); 656% startseed = randn('state'); % for retrieving in saved variables 657 658 % Initialize further constants 659 chiN=N^0.5*(1-1/(4*N)+1/(21*N^2)); % expectation of 660 % ||N(0,I)|| == norm(randn(N,1)) 661 662 countiter = 0; 663 % Initialize records and output 664 if irun == 1 665 time.t0 = clock; 666 667 % TODO: keep also median solution? 668 out.evals = counteval; % should be first entry 669 out.stopflag = {}; 670 671 outiter = 0; 672 673 % Write headers to output data files 674 filenameprefix = opts.LogFilenamePrefix; 675 if savemodulo && savetime 676 filenames = {}; 677 filenames(end+1) = {'axlen'}; 678 filenames(end+1) = {'fit'}; 679 filenames(end+1) = {'stddev'}; 680 filenames(end+1) = {'xmean'}; 681 filenames(end+1) = {'xrecentbest'}; 682 str = [' (startseed=' num2str(startseed(2)) ... 683 ', ' num2str(clock, '%d/%02d/%d %d:%d:%2.2f') ')']; 684 for namecell = filenames(:)' 685 name = namecell{:}; 686 687 [fid, err] = fopen(['./' filenameprefix name '.dat'], 'w'); 688 if fid < 1 % err ~= 0 689 warning(['could not open ' filenameprefix name '.dat']); 690 filenames(find(strcmp(filenames,name))) = []; 691 else 692 % fprintf(fid, '%s\n', ... 693 % ['<CMAES-OUTPUT version="' cmaVersion '">']); 694 % fprintf(fid, [' <NAME>' name '</NAME>\n']); 695 % fprintf(fid, [' <DATE>' date() '</DATE>\n']); 696 % fprintf(fid, ' <PARAMETERS>\n'); 697 % fprintf(fid, [' dimension=' num2str(N) '\n']); 698 % fprintf(fid, ' </PARAMETERS>\n'); 699 % different cases for DATA columns annotations here 700 % fprintf(fid, ' <DATA'); 701 if strcmp(name, 'axlen') 702 fprintf(fid, ['%% columns="iteration, evaluation, sigma, ' ... 703 'max axis length, min axis length, ' ... 704 'all principal axes lengths (sorted square roots ' ... 705 'of eigenvalues of C)"' str]); 706 elseif strcmp(name, 'fit') 707 fprintf(fid, ['%% columns="iteration, evaluation, sigma, axis ratio, bestever,' ... 708 ' best, median, worst fitness function value,' ... 709 ' further objective values of best"' str]); 710 elseif strcmp(name, 'stddev') 711 fprintf(fid, ['%% columns=["iteration, evaluation, sigma, void, void, ' ... 712 'stds==sigma*sqrt(diag(C))"' str]); 713 elseif strcmp(name, 'xmean') 714 fprintf(fid, ['%% columns="iteration, evaluation, void, ' ... 715 'void, void, xmean"' str]); 716 elseif strcmp(name, 'xrecentbest') 717 fprintf(fid, ['%% columns="iteration, evaluation, fitness, ' ... 718 'void, void, xrecentbest"' str]); 719 end 720 fprintf(fid, '\n'); % DATA 721 if strcmp(name, 'xmean') 722 fprintf(fid, '%ld %ld 0 0 0 ', 0, counteval); 723 % fprintf(fid, '%ld %ld 0 0 %e ', countiter, counteval, fmean); 724 %qqq fprintf(fid, msprintf('%e ', genophenotransform(out.genopheno, xmean)) + '\n'); 725 fprintf(fid, '%e ', xmean); 726 fprintf(fid, '\n'); 727 end 728 fclose(fid); 729 clear fid; % preventing 730 end 731 end % for files 732 end % savemodulo 733 end % irun == 1 734 735 end % else flgresume 736 737 % -------------------- Generation Loop -------------------------------- 738 stopflag = {}; 739 while isempty(stopflag) 740 % set internal parameters 741 if countiter == 0 || lambda ~= lambda_last 742 if countiter > 0 && floor(log10(lambda)) ~= floor(log10(lambda_last)) ... 743 && flgdisplay 744 disp([' lambda = ' num2str(lambda)]); 745 lambda_hist(:,end+1) = [countiter+1; lambda]; 746 else 747 lambda_hist = [countiter+1; lambda]; 748 end 749 lambda_last = lambda; 750 % Strategy internal parameter setting: Selection 751 mu = myeval(opts.ParentNumber); % number of parents/points for recombination 752 if strncmp(lower(opts.RecombinationWeights), 'equal', 3) 753 weights = ones(mu,1); % (mu_I,lambda)-CMA-ES 754 elseif strncmp(lower(opts.RecombinationWeights), 'linear', 3) 755 weights = mu+0.5-(1:mu)'; 756 elseif strncmp(lower(opts.RecombinationWeights), 'superlinear', 3) 757 weights = log(mu+0.5)-log(1:mu)'; % muXone array for weighted recombination 758 % qqq mu can be non-integer and 759 % should become ceil(mu-0.5) (minor correction) 760 else 761 error(['Recombination weights to be "' opts.RecombinationWeights ... 762 '" is not implemented']); 763 end 764 mueff=sum(weights)^2/sum(weights.^2); % variance-effective size of mu 765 weights = weights/sum(weights); % normalize recombination weights array 766 if mueff == lambda 767 error(['Combination of values for PopSize, ParentNumber and ' ... 768 ' and RecombinationWeights is not reasonable']); 769 end 770 771 % Strategy internal parameter setting: Adaptation 772 cc = myeval(opts.CMA.ccum); % time constant for cumulation for covariance matrix 773 cs = myeval(opts.CMA.cs); 774 775 % old way TODO: remove this at some point 776 % mucov = mueff; % size of mu used for calculating learning rate ccov 777 % ccov = (1/mucov) * 2/(N+1.41)^2 ... % learning rate for covariance matrix 778 % + (1-1/mucov) * min(1,(2*mucov-1)/((N+2)^2+mucov)); 779 780 % new way 781 if myevalbool(opts.CMA.on) 782 ccov1 = myeval(opts.CMA.ccov1); 783 ccovmu = min(1-ccov1, myeval(opts.CMA.ccovmu)); 784 else 785 ccov1 = 0; 786 ccovmu = 0; 787 end 788 789 % flgDiagonalOnly = -lambda*4*1/ccov; % for ccov==1 it is not needed 790 % 0 : C will never be diagonal anymore 791 % 1 : C will always be diagonal 792 % >1: C is diagonal for first iterations, set to 0 afterwards 793 if flgDiagonalOnly < 1 794 flgDiagonalOnly = 0; 795 end 796 if flgDiagonalOnly 797 ccov1_sep = min(1, ccov1 * (N+1.5) / 3); 798 ccovmu_sep = min(1-ccov1_sep, ccovmu * (N+1.5) / 3); 799 elseif N > 98 && flgdisplay && countiter == 0 800 disp('consider option DiagonalOnly for high-dimensional problems'); 801 end 802 803 % ||ps|| is close to sqrt(mueff/N) for mueff large on linear fitness 804 %damps = ... % damping for step size control, usually close to one 805 % (1 + 2*max(0,sqrt((mueff-1)/(N+1))-1)) ... % limit sigma increase 806 % * max(0.3, ... % reduce damps, if max. iteration number is small 807 % 1 - N/min(stopMaxIter,stopMaxFunEvals/lambda)) + cs; 808 damps = myeval(opts.CMA.damps); 809 if noiseHandling 810 noiseReevals = min(myeval(opts.Noise.reevals), lambda); 811 noiseAlpha = myeval(opts.Noise.alphasigma); 812 noiseEpsilon = myeval(opts.Noise.epsilon); 813 noiseTheta = myeval(opts.Noise.theta); 814 noisecum = myeval(opts.Noise.cum); 815 noiseCutOff = myeval(opts.Noise.cutoff); % arguably of minor relevance 816 else 817 noiseReevals = 0; % more convenient in later coding 818 end 819 820 %qqq hacking of a different parameter setting, e.g. for ccov or damps, 821 % can be done here, but is not necessary anymore, see opts.CMA. 822 % ccov1 = 0.0*ccov1; disp(['CAVE: ccov1=' num2str(ccov1)]); 823 % ccovmu = 0.0*ccovmu; disp(['CAVE: ccovmu=' num2str(ccovmu)]); 824 % damps = inf*damps; disp(['CAVE: damps=' num2str(damps)]); 825 % cc = 1; disp(['CAVE: cc=' num2str(cc)]); 826 827 end 828 829 % Display initial message 830 if countiter == 0 && flgdisplay 831 if mu == 1 832 strw = '100'; 833 elseif mu < 8 834 strw = [sprintf('%.0f', 100*weights(1)) ... 835 sprintf(' %.0f', 100*weights(2:end)')]; 836 else 837 strw = [sprintf('%.2g ', 100*weights(1:2)') ... 838 sprintf('%.2g', 100*weights(3)') '...' ... 839 sprintf(' %.2g', 100*weights(end-1:end)') ']%, ']; 840 end 841 if irun > 1 842 strrun = [', run ' num2str(irun)]; 843 else 844 strrun = ''; 845 end 846 disp([' n=' num2str(N) ': (' num2str(mu) ',' ... 847 num2str(lambda) ')-CMA-ES(w=[' ... 848 strw ']%, ' ... 849 'mu_eff=' num2str(mueff,'%.1f') ... 850 ') on function ' ... 851 (fitfun) strrun]); 852 if flgDiagonalOnly == 1 853 disp(' C is diagonal'); 854 elseif flgDiagonalOnly 855 disp([' C is diagonal for ' num2str(floor(flgDiagonalOnly)) ' iterations']); 856 end 857 end 858 859 flush; 860 861 countiter = countiter + 1; 862 863 % Generate and evaluate lambda offspring 864 865 fitness.raw = repmat(NaN, 1, lambda + noiseReevals); 866 867 % parallel evaluation 868 if flgEvalParallel 869 arz = randn(N,lambda); 870 871 if ~flgDiagonalOnly 872 arx = repmat(xmean, 1, lambda) + sigma * (BD * arz); % Eq. (1) 873 else 874 arx = repmat(xmean, 1, lambda) + repmat(sigma * diagD, 1, lambda) .* arz; 875 end 876 877 if noiseHandling 878 if noiseEpsilon == 0 879 arx = [arx arx(:,1:noiseReevals)]; 880 elseif flgDiagonalOnly 881 arx = [arx arx(:,1:noiseReevals) + ... 882 repmat(noiseEpsilon * sigma * diagD, 1, noiseReevals) ... 883 .* randn(N,noiseReevals)]; 884 else 885 arx = [arx arx(:,1:noiseReevals) + ... 886 noiseEpsilon * sigma * ... 887 (BD * randn(N,noiseReevals))]; 888 end 889 end 890 891 % You may handle constraints here. You may either resample 892 % arz(:,k) and/or multiply it with a factor between -1 and 1 893 % (the latter will decrease the overall step size) and 894 % recalculate arx accordingly. Do not change arx or arz in any 895 % other way. 896 897 if ~bnd.isactive 898 arxvalid = arx; 899 else 900 arxvalid = xintobounds(arx, lbounds, ubounds); 901 end 902 % You may handle constraints here. You may copy and alter 903 % (columns of) arxvalid(:,k) only for the evaluation of the 904 % fitness function. arx and arxvalid should not be changed. 905 fitness.raw = feval(fitfun, arxvalid, varargin{:}); 906 countevalNaN = countevalNaN + sum(isnan(fitness.raw)); 907 counteval = counteval + sum(~isnan(fitness.raw)); 908 end 909 910 % non-parallel evaluation and remaining NaN-values 911 % set also the reevaluated solution to NaN 912 fitness.raw(lambda + find(isnan(fitness.raw(1:noiseReevals)))) = NaN; 913 for k=find(isnan(fitness.raw)) 914 % fitness.raw(k) = NaN; 915 tries = 0; 916 % Resample, until fitness is not NaN 917 while isnan(fitness.raw(k)) 918 if k <= lambda % regular samples (not the re-evaluation-samples) 919 arz(:,k) = randn(N,1); % (re)sample 920 921 if flgDiagonalOnly 922 arx(:,k) = xmean + sigma * diagD .* arz(:,k); % Eq. (1) 923 else 924 arx(:,k) = xmean + sigma * (BD * arz(:,k)); % Eq. (1) 925 end 926 else % re-evaluation solution with index > lambda 927 if flgDiagonalOnly 928 arx(:,k) = arx(:,k-lambda) + (noiseEpsilon * sigma) * diagD .* randn(N,1); 929 else 930 arx(:,k) = arx(:,k-lambda) + (noiseEpsilon * sigma) * (BD * randn(N,1)); 931 end 932 end 933 934 % You may handle constraints here. You may either resample 935 % arz(:,k) and/or multiply it with a factor between -1 and 1 936 % (the latter will decrease the overall step size) and 937 % recalculate arx accordingly. Do not change arx or arz in any 938 % other way. 939 940 if ~bnd.isactive 941 arxvalid(:,k) = arx(:,k); 942 else 943 arxvalid(:,k) = xintobounds(arx(:,k), lbounds, ubounds); 944 end 945 % You may handle constraints here. You may copy and alter 946 % (columns of) arxvalid(:,k) only for the evaluation of the 947 % fitness function. arx should not be changed. 948 fitness.raw(k) = feval(fitfun, arxvalid(:,k), varargin{:}); 949 tries = tries + 1; 950 if isnan(fitness.raw(k)) 951 countevalNaN = countevalNaN + 1; 952 end 953 if mod(tries, 100) == 0 954 warning([num2str(tries) ... 955 ' NaN objective function values at evaluation ' ... 956 num2str(counteval)]); 957 end 958 end 959 counteval = counteval + 1; % retries due to NaN are not counted 960 end 961 962 fitness.sel = fitness.raw; 963 964 % ----- handle boundaries ----- 965 if 1 < 3 && bnd.isactive 966 % Get delta fitness values 967 val = myprctile(fitness.raw, [25 75]); 968 % more precise would be exp(mean(log(diagC))) 969 val = (val(2) - val(1)) / N / mean(diagC) / sigma^2; 970 %val = (myprctile(fitness.raw, 75) - myprctile(fitness.raw, 25)) ... 971 % / N / mean(diagC) / sigma^2; 972 % Catch non-sensible values 973 if ~isfinite(val) 974 warning('Non-finite fitness range'); 975 val = max(bnd.dfithist); 976 elseif val == 0 % happens if all points are out of bounds 977 val = min(bnd.dfithist(bnd.dfithist>0)); % seems not to make sense, given all solutions are out of bounds 978 elseif bnd.validfitval == 0 % flag that first sensible val was found 979 bnd.dfithist = []; 980 bnd.validfitval = 1; 981 end 982 983 % Store delta fitness values 984 if length(bnd.dfithist) < 20+(3*N)/lambda 985 bnd.dfithist = [bnd.dfithist val]; 986 else 987 bnd.dfithist = [bnd.dfithist(2:end) val]; 988 end 989 990 [tx, ti] = xintobounds(xmean, lbounds, ubounds); 991 992 % Set initial weights 993 if bnd.iniphase 994 if any(ti) 995 bnd.weights(find(bnd.isbounded)) = 2.0002 * median(bnd.dfithist); 996 if bnd.flgscale == 0 % scale only initial weights then 997 dd = diagC; 998 idx = find(bnd.isbounded); 999 dd = dd(idx) / mean(dd); % remove mean scaling 1000 bnd.weights(idx) = bnd.weights(idx) ./ dd; 1001 end 1002 if bnd.validfitval && countiter > 2 1003 bnd.iniphase = 0; 1004 end 1005 end 1006 end 1007 1008 % Increase weights 1009 if 1 < 3 && any(ti) % any coordinate of xmean out of bounds 1010 % judge distance of xmean to boundary 1011 tx = xmean - tx; 1012 idx = (ti ~= 0 & abs(tx) > 3*max(1,sqrt(N)/mueff) ... 1013 * sigma*sqrt(diagC)) ; 1014 % only increase if xmean is moving away 1015 idx = idx & (sign(tx) == sign(xmean - xold)); 1016 if ~isempty(idx) % increase 1017 % the factor became 1.2 instead of 1.1, because 1018 % changed from max to min in version 3.52 1019 bnd.weights(idx) = 1.2^(min(1, mueff/10/N)) * bnd.weights(idx); 1020 end 1021 end 1022 1023 % Calculate scaling biased to unity, product is one 1024 if bnd.flgscale ~= 0 1025 bnd.scale = exp(0.9*(log(diagC)-mean(log(diagC)))); 1026 end 1027 1028 % Assigned penalized fitness 1029 bnd.arpenalty = (bnd.weights ./ bnd.scale)' * (arxvalid - arx).^2; 1030 1031 fitness.sel = fitness.raw + bnd.arpenalty; 1032 1033 end % handle boundaries 1034 % ----- end handle boundaries ----- 1035 1036 % compute noise measurement and reduce fitness arrays to size lambda 1037 if noiseHandling 1038 [noiseS] = local_noisemeasurement(fitness.sel(1:lambda), ... 1039 fitness.sel(lambda+(1:noiseReevals)), ... 1040 noiseReevals, noiseTheta, noiseCutOff); 1041 if countiter == 1 % TODO: improve this very rude way of initialization 1042 noiseSS = 0; 1043 noiseN = 0; % counter for mean 1044 end 1045 noiseSS = noiseSS + noisecum * (noiseS - noiseSS); 1046 1047 % noise-handling could be done here, but the original sigma is still needed 1048 % disp([noiseS noiseSS noisecum]) 1049 1050 fitness.rawar12 = fitness.raw; % just documentary 1051 fitness.selar12 = fitness.sel; % just documentary 1052 % qqq refine fitness based on both values 1053 if 11 < 3 % TODO: in case of outliers this mean is counterproductive 1054 % median out of three would be ok 1055 fitness.raw(1:noiseReevals) = ... % not so raw anymore 1056 (fitness.raw(1:noiseReevals) + fitness.raw(lambda+(1:noiseReevals))) / 2; 1057 fitness.sel(1:noiseReevals) = ... 1058 (fitness.sel(1:noiseReevals) + fitness.sel(lambda+(1:noiseReevals))) / 2; 1059 end 1060 fitness.raw = fitness.raw(1:lambda); 1061 fitness.sel = fitness.sel(1:lambda); 1062 end 1063 1064 % Sort by fitness 1065 [fitness.raw, fitness.idx] = sort(fitness.raw); 1066 [fitness.sel, fitness.idxsel] = sort(fitness.sel); % minimization 1067 fitness.hist(2:end) = fitness.hist(1:end-1); % record short history of 1068 fitness.hist(1) = fitness.raw(1); % best fitness values 1069 if length(fitness.histbest) < 120+ceil(30*N/lambda) || ... 1070 (mod(countiter, 5) == 0 && length(fitness.histbest) < 2e4) % 20 percent of 1e5 gen. 1071 fitness.histbest = [fitness.raw(1) fitness.histbest]; % best fitness values 1072 fitness.histmedian = [median(fitness.raw) fitness.histmedian]; % median fitness values 1073 else 1074 fitness.histbest(2:end) = fitness.histbest(1:end-1); 1075 fitness.histmedian(2:end) = fitness.histmedian(1:end-1); 1076 fitness.histbest(1) = fitness.raw(1); % best fitness values 1077 fitness.histmedian(1) = median(fitness.raw); % median fitness values 1078 end 1079 fitness.histsel(2:end) = fitness.histsel(1:end-1); % record short history of 1080 fitness.histsel(1) = fitness.sel(1); % best sel fitness values 1081 1082 % Calculate new xmean, this is selection and recombination 1083 xold = xmean; % for speed up of Eq. (2) and (3) 1084 xmean = arx(:,fitness.idxsel(1:mu))*weights; 1085 zmean = arz(:,fitness.idxsel(1:mu))*weights;%==D^-1*B'*(xmean-xold)/sigma 1086 if mu == 1 1087 fmean = fitness.sel(1); 1088 else 1089 fmean = NaN; % [] does not work in the latter assignment 1090 % fmean = feval(fitfun, xintobounds(xmean, lbounds, ubounds), varargin{:}); 1091 % counteval = counteval + 1; 1092 end 1093 1094 % Cumulation: update evolution paths 1095 ps = (1-cs)*ps + sqrt(cs*(2-cs)*mueff) * (B*zmean); % Eq. (4) 1096 hsig = norm(ps)/sqrt(1-(1-cs)^(2*countiter))/chiN < 1.4 + 2/(N+1); 1097 if flg_future_setting 1098 hsig = sum(ps.^2) / (1-(1-cs)^(2*countiter)) / N < 2 + 4/(N+1); % just simplified 1099 end 1100 % hsig = norm(ps)/sqrt(1-(1-cs)^(2*countiter))/chiN < 1.4 + 2/(N+1); 1101 % hsig = norm(ps)/sqrt(1-(1-cs)^(2*countiter))/chiN < 1.5 + 1/(N-0.5); 1102 % hsig = norm(ps) < 1.5 * sqrt(N); 1103 % hsig = 1; 1104 1105 pc = (1-cc)*pc ... 1106 + hsig*(sqrt(cc*(2-cc)*mueff)/sigma) * (xmean-xold); % Eq. (2) 1107 if hsig == 0 1108 % disp([num2str(countiter) ' ' num2str(counteval) ' pc update stalled']); 1109 end 1110 1111 % Adapt covariance matrix 1112 neg.ccov = 0; % TODO: move parameter setting upwards at some point 1113 if ccov1 + ccovmu > 0 % Eq. (3) 1114 if flgDiagonalOnly % internal linear(?) complexity 1115 diagC = (1-ccov1_sep-ccovmu_sep+(1-hsig)*ccov1_sep*cc*(2-cc)) * diagC ... % regard old matrix 1116 + ccov1_sep * pc.^2 ... % plus rank one update 1117 + ccovmu_sep ... % plus rank mu update 1118 * (diagC .* (arz(:,fitness.idxsel(1:mu)).^2 * weights)); 1119 % * (repmat(diagC,1,mu) .* arz(:,fitness.idxsel(1:mu)).^2 * weights); 1120 diagD = sqrt(diagC); % replaces eig(C) 1121 else 1122 arpos = (arx(:,fitness.idxsel(1:mu))-repmat(xold,1,mu)) / sigma; 1123 % "active" CMA update: negative update, in case controlling pos. definiteness 1124 if flgActiveCMA > 0 1125 % set parameters 1126 neg.mu = mu; 1127 neg.mueff = mueff; 1128 if flgActiveCMA > 10 % flat weights with mu=lambda/2 1129 neg.mu = floor(lambda/2); 1130 neg.mueff = neg.mu; 1131 end 1132 % neg.mu = ceil(min([N, lambda/4, mueff])); neg.mueff = mu; % i.e. neg.mu <= N 1133 % Parameter study: in 3-D lambda=50,100, 10-D lambda=200,400, 30-D lambda=1000,2000 a 1134 % three times larger neg.ccov does not work. 1135 % increasing all ccov rates three times does work (probably because of the factor (1-ccovmu)) 1136 % in 30-D to looks fine 1137 1138 neg.ccov = (1 - ccovmu) * 0.25 * neg.mueff / ((N+2)^1.5 + 2*neg.mueff); 1139 neg.minresidualvariance = 0.66; % keep at least 0.66 in all directions, small popsize are most critical 1140 neg.alphaold = 0.5; % where to make up for the variance loss, 0.5 means no idea what to do 1141 % 1 is slightly more robust and gives a better "guaranty" for pos. def., 1142 % but does it make sense from the learning perspective for large ccovmu? 1143 1144 neg.ccovfinal = neg.ccov; 1145 1146 % prepare vectors, compute negative updating matrix Cneg and checking matrix Ccheck 1147 arzneg = arz(:,fitness.idxsel(lambda:-1:lambda - neg.mu + 1)); 1148 % i-th longest becomes i-th shortest 1149 % TODO: this is not in compliance with the paper Hansen&Ros2010, 1150 % where simply arnorms = arnorms(end:-1:1) ? 1151 [arnorms, idxnorms] = sort(sqrt(sum(arzneg.^2, 1))); 1152 [ignore, idxnorms] = sort(idxnorms); % inverse index 1153 arnormfacs = arnorms(end:-1:1) ./ arnorms; 1154 % arnormfacs = arnorms(randperm(neg.mu)) ./ arnorms; 1155 arnorms = arnorms(end:-1:1); % for the record 1156 if flgActiveCMA < 20 1157 arzneg = arzneg .* repmat(arnormfacs(idxnorms), N, 1); % E x*x' is N 1158 % arzneg = sqrt(N) * arzneg ./ repmat(sqrt(sum(arzneg.^2, 1)), N, 1); % E x*x' is N 1159 end 1160 if flgActiveCMA < 10 && neg.mu == mu % weighted sum 1161 if mod(flgActiveCMA, 10) == 1 % TODO: prevent this with a less tight but more efficient check (see below) 1162 Ccheck = arzneg * diag(weights) * arzneg'; % in order to check the largest EV 1163 end 1164 artmp = BD * arzneg; 1165 Cneg = artmp * diag(weights) * artmp'; 1166 else % simple sum 1167 if mod(flgActiveCMA, 10) == 1 1168 Ccheck = (1/neg.mu) * arzneg*arzneg'; % in order to check largest EV 1169 end 1170 artmp = BD * arzneg; 1171 Cneg = (1/neg.mu) * artmp*artmp'; 1172 1173 end 1174 1175 % check pos.def. and set learning rate neg.ccov accordingly, 1176 % this check makes the original choice of neg.ccov extremly failsafe 1177 % still assuming C == BD*BD', which is only approxim. correct 1178 if mod(flgActiveCMA, 10) == 1 && 1 - neg.ccov * arnorms(idxnorms).^2 * weights < neg.minresidualvariance 1179 % TODO: the simple and cheap way would be to set 1180 % fac = 1 - ccovmu - ccov1 OR 1 - mueff/lambda and 1181 % neg.ccov = fac*(1 - neg.minresidualvariance) / (arnorms(idxnorms).^2 * weights) 1182 % this is the more sophisticated way: 1183 % maxeigenval = eigs(arzneg * arzneg', 1, 'lm', eigsopts); % not faster 1184 maxeigenval = max(eig(Ccheck)); % norm is much slower, because (norm()==max(svd()) 1185 %disp([countiter log10([neg.ccov, maxeigenval, arnorms(idxnorms).^2 * weights, max(arnorms)^2]), ... 1186 % neg.ccov * arnorms(idxnorms).^2 * weights]) 1187 % pause 1188 % remove less than ??34*(1-cmu)%?? of variance in any direction 1189 % 1-ccovmu is the variance left from the old C 1190 neg.ccovfinal = min(neg.ccov, (1-ccovmu)*(1-neg.minresidualvariance)/maxeigenval); 1191 % -ccov1 removed to avoid error message?? 1192 if neg.ccovfinal < neg.ccov 1193 disp(['active CMA at iteration ' num2str(countiter) ... 1194 ': max EV ==', num2str([maxeigenval, neg.ccov, neg.ccovfinal])]); 1195 end 1196 end 1197 % xmean = xold; % the distribution does not degenerate!? 1198 % update C 1199 C = (1-ccov1-ccovmu+neg.alphaold*neg.ccovfinal+(1-hsig)*ccov1*cc*(2-cc)) * C ... % regard old matrix 1200 + ccov1 * pc*pc' ... % plus rank one update 1201 + (ccovmu + (1-neg.alphaold)*neg.ccovfinal) ... % plus rank mu update 1202 * arpos * (repmat(weights,1,N) .* arpos') ... 1203 - neg.ccovfinal * Cneg; % minus rank mu update 1204 else % no active (negative) update 1205 C = (1-ccov1-ccovmu+(1-hsig)*ccov1*cc*(2-cc)) * C ... % regard old matrix 1206 + ccov1 * pc*pc' ... % plus rank one update 1207 + ccovmu ... % plus rank mu update 1208 * arpos * (repmat(weights,1,N) .* arpos'); 1209 % is now O(mu*N^2 + mu*N), was O(mu*N^2 + mu^2*N) when using diag(weights) 1210 % for mu=30*N it is now 10 times faster, overall 3 times faster 1211 end 1212 diagC = diag(C); 1213 end 1214 end 1215 1216 % the following is de-preciated and will be removed in future 1217 % better setting for cc makes this hack obsolete 1218 if 11 < 2 && ~flgscience 1219 % remove momentum in ps, if ps is large and fitness is getting worse. 1220 % this should rarely happen. 1221 % this might very well be counterproductive in dynamic environments 1222 if sum(ps.^2)/N > 1.5 + 10*(2/N)^.5 && ... 1223 fitness.histsel(1) > max(fitness.histsel(2:3)) 1224 ps = ps * sqrt(N*(1+max(0,log(sum(ps.^2)/N))) / sum(ps.^2)); 1225 if flgdisplay 1226 disp(['Momentum in ps removed at [niter neval]=' ... 1227 num2str([countiter counteval]) ']']); 1228 end 1229 end 1230 end 1231 1232 % Adapt sigma 1233 if flg_future_setting % according to a suggestion from Dirk Arnold (2000) 1234 % exp(1) is still not reasonably small enough 1235 sigma = sigma * exp(min(1, (sum(ps.^2)/N - 1)/2 * cs/damps)); % Eq. (5) 1236 else 1237 % exp(1) is still not reasonably small enough 1238 sigma = sigma * exp(min(1, (sqrt(sum(ps.^2))/chiN - 1) * cs/damps)); % Eq. (5) 1239 end 1240 % disp([countiter norm(ps)/chiN]); 1241 1242 if 11 < 3 % testing with optimal step-size 1243 if countiter == 1 1244 disp('*********** sigma set to const * ||x|| ******************'); 1245 end 1246 sigma = 0.04 * mueff * sqrt(sum(xmean.^2)) / N; % 20D,lam=1000:25e3 1247 sigma = 0.3 * mueff * sqrt(sum(xmean.^2)) / N; % 20D,lam=(40,1000):17e3 1248 % 75e3 with def (1.5) 1249 % 35e3 with damps=0.25 1250 end 1251 if 11 < 3 1252 if countiter == 1 1253 disp('*********** xmean set to const ******************'); 1254 end 1255 xmean = ones(N,1); 1256 end 1257 1258 % Update B and D from C 1259 1260 if ~flgDiagonalOnly && (ccov1+ccovmu+neg.ccov) > 0 && mod(countiter, 1/(ccov1+ccovmu+neg.ccov)/N/10) < 1 1261 C=triu(C)+triu(C,1)'; % enforce symmetry to prevent complex numbers 1262 [B,tmp] = eig(C); % eigen decomposition, B==normalized eigenvectors 1263 % effort: approx. 15*N matrix-vector multiplications 1264 diagD = diag(tmp); 1265 1266 if any(~isfinite(diagD)) 1267 clear idx; % prevents error under octave 1268 save(['tmp' opts.SaveFilename]); 1269 error(['function eig returned non-finited eigenvalues, cond(C)=' ... 1270 num2str(cond(C)) ]); 1271 end 1272 if any(any(~isfinite(B))) 1273 clear idx; % prevents error under octave 1274 save(['tmp' opts.SaveFilename]); 1275 error(['function eig returned non-finited eigenvectors, cond(C)=' ... 1276 num2str(cond(C)) ]); 1277 end 1278 1279 % limit condition of C to 1e14 + 1 1280 if min(diagD) <= 0 1281 if stopOnWarnings 1282 stopflag(end+1) = {'warnconditioncov'}; 1283 else 1284 warning(['Iteration ' num2str(countiter) ... 1285 ': Eigenvalue (smaller) zero']); 1286 diagD(diagD<0) = 0; 1287 tmp = max(diagD)/1e14; 1288 C = C + tmp*eye(N,N); diagD = diagD + tmp*ones(N,1); 1289 end 1290 end 1291 if max(diagD) > 1e14*min(diagD) 1292 if stopOnWarnings 1293 stopflag(end+1) = {'warnconditioncov'}; 1294 else 1295 warning(['Iteration ' num2str(countiter) ': condition of C ' ... 1296 'at upper limit' ]); 1297 tmp = max(diagD)/1e14 - min(diagD); 1298 C = C + tmp*eye(N,N); diagD = diagD + tmp*ones(N,1); 1299 end 1300 end 1301 1302 diagC = diag(C); 1303 diagD = sqrt(diagD); % D contains standard deviations now 1304 % diagD = diagD / prod(diagD)^(1/N); C = C / prod(diagD)^(2/N); 1305 BD = B.*repmat(diagD',N,1); % O(n^2) 1306 end % if mod 1307 1308 % Align/rescale order of magnitude of scales of sigma and C for nicer output 1309 % not a very usual case 1310 if 1 < 2 && sigma > 1e10*max(diagD) 1311 fac = sigma / max(diagD); 1312 sigma = sigma/fac; 1313 pc = fac * pc; 1314 diagD = fac * diagD; 1315 if ~flgDiagonalOnly 1316 C = fac^2 * C; % disp(fac); 1317 BD = B.*repmat(diagD',N,1); % O(n^2), but repmat might be inefficient todo? 1318 end 1319 diagC = fac^2 * diagC; 1320 end 1321 1322 if flgDiagonalOnly > 1 && countiter > flgDiagonalOnly 1323 % full covariance matrix from now on 1324 flgDiagonalOnly = 0; 1325 B = eye(N,N); 1326 BD = diag(diagD); 1327 C = diag(diagC); % is better, because correlations are spurious anyway 1328 end 1329 1330 if noiseHandling 1331 if countiter == 1 % assign firstvarargin for noise treatment e.g. as #reevaluations 1332 if ~isempty(varargin) && length(varargin{1}) == 1 && isnumeric(varargin{1}) 1333 if irun == 1 1334 firstvarargin = varargin{1}; 1335 else 1336 varargin{1} = firstvarargin; % reset varargin{1} 1337 end 1338 else 1339 firstvarargin = 0; 1340 end 1341 end 1342 if noiseSS < 0 && noiseMinMaxEvals(2) > noiseMinMaxEvals(1) && firstvarargin 1343 varargin{1} = max(noiseMinMaxEvals(1), varargin{1} / noiseAlphaEvals^(1/4)); % still experimental 1344 elseif noiseSS > 0 1345 if ~isempty(noiseCallback) % to be removed? 1346 res = feval(noiseCallback); % should also work without output argument!? 1347 if ~isempty(res) && res > 1 % TODO: decide for interface of callback 1348 % also a dynamic popsize could be done here 1349 sigma = sigma * noiseAlpha; 1350 end 1351 else 1352 if noiseMinMaxEvals(2) > noiseMinMaxEvals(1) && firstvarargin 1353 varargin{1} = min(noiseMinMaxEvals(2), varargin{1} * noiseAlphaEvals); 1354 end 1355 1356 sigma = sigma * noiseAlpha; 1357 % lambda = ceil(0.1 * sqrt(lambda) + lambda); 1358 % TODO: find smallest increase of lambda with log-linear 1359 % convergence in iterations 1360 end 1361 % qqq experimental: take a mean to estimate the true optimum 1362 noiseN = noiseN + 1; 1363 if noiseN == 1 1364 noiseX = xmean; 1365 else 1366 noiseX = noiseX + (3/noiseN) * (xmean - noiseX); 1367 end 1368 end 1369 end 1370 1371 % ----- numerical error management ----- 1372 % Adjust maximal coordinate axis deviations 1373 if any(sigma*sqrt(diagC) > maxdx) 1374 sigma = min(maxdx ./ sqrt(diagC)); 1375 %warning(['Iteration ' num2str(countiter) ': coordinate axis std ' ... 1376 % 'deviation at upper limit of ' num2str(maxdx)]); 1377 % stopflag(end+1) = {'maxcoorddev'}; 1378 end 1379 % Adjust minimal coordinate axis deviations 1380 if any(sigma*sqrt(diagC) < mindx) 1381 sigma = max(mindx ./ sqrt(diagC)) * exp(0.05+cs/damps); 1382 %warning(['Iteration ' num2str(countiter) ': coordinate axis std ' ... 1383 % 'deviation at lower limit of ' num2str(mindx)]); 1384 % stopflag(end+1) = {'mincoorddev'};; 1385 end 1386 % Adjust too low coordinate axis deviations 1387 if any(xmean == xmean + 0.2*sigma*sqrt(diagC)) 1388 if stopOnWarnings 1389 stopflag(end+1) = {'warnnoeffectcoord'}; 1390 else 1391 warning(['Iteration ' num2str(countiter) ': coordinate axis std ' ... 1392 'deviation too low' ]); 1393 if flgDiagonalOnly 1394 diagC = diagC + (ccov1_sep+ccovmu_sep) * (diagC .* ... 1395 (xmean == xmean + 0.2*sigma*sqrt(diagC))); 1396 else 1397 C = C + (ccov1+ccovmu) * diag(diagC .* ... 1398 (xmean == xmean + 0.2*sigma*sqrt(diagC))); 1399 end 1400 sigma = sigma * exp(0.05+cs/damps); 1401 end 1402 end 1403 % Adjust step size in case of (numerical) precision problem 1404 if flgDiagonalOnly 1405 tmp = 0.1*sigma*diagD; 1406 else 1407 tmp = 0.1*sigma*BD(:,1+floor(mod(countiter,N))); 1408 end 1409 if all(xmean == xmean + tmp) 1410 i = 1+floor(mod(countiter,N)); 1411 if stopOnWarnings 1412 stopflag(end+1) = {'warnnoeffectaxis'}; 1413 else 1414 warning(['Iteration ' num2str(countiter) ... 1415 ': main axis standard deviation ' ... 1416 num2str(sigma*diagD(i)) ' has no effect' ]); 1417 sigma = sigma * exp(0.2+cs/damps); 1418 end 1419 end 1420 % Adjust step size in case of equal function values (flat fitness) 1421 % isequalfuncvalues = 0; 1422 if fitness.sel(1) == fitness.sel(1+ceil(0.1+lambda/4)) 1423 % isequalfuncvalues = 1; 1424 if stopOnEqualFunctionValues 1425 arrEqualFunvals = [countiter arrEqualFunvals(1:end-1)]; 1426 % stop if this happens in more than 33% 1427 if arrEqualFunvals(end) > countiter - 3 * length(arrEqualFunvals) 1428 stopflag(end+1) = {'equalfunvals'}; 1429 end 1430 else 1431 if flgWarnOnEqualFunctionValues 1432 warning(['Iteration ' num2str(countiter) ... 1433 ': equal function values f=' num2str(fitness.sel(1)) ... 1434 ' at maximal main axis sigma ' ... 1435 num2str(sigma*max(diagD))]); 1436 end 1437 sigma = sigma * exp(0.2+cs/damps); 1438 end 1439 end 1440 % Adjust step size in case of equal function values 1441 if countiter > 2 && myrange([fitness.hist fitness.sel(1)]) == 0 1442 if stopOnWarnings 1443 stopflag(end+1) = {'warnequalfunvalhist'}; 1444 else 1445 warning(['Iteration ' num2str(countiter) ... 1446 ': equal function values in history at maximal main ' ... 1447 'axis sigma ' num2str(sigma*max(diagD))]); 1448 sigma = sigma * exp(0.2+cs/damps); 1449 end 1450 end 1451 1452 % ----- end numerical error management ----- 1453 1454 % Keep overall best solution 1455 out.evals = counteval; 1456 out.solutions.evals = counteval; 1457 out.solutions.mean.x = xmean; 1458 out.solutions.mean.f = fmean; 1459 out.solutions.mean.evals = counteval; 1460 out.solutions.recentbest.x = arxvalid(:, fitness.idx(1)); 1461 out.solutions.recentbest.f = fitness.raw(1); 1462 out.solutions.recentbest.evals = counteval + fitness.idx(1) - lambda; 1463 out.solutions.recentworst.x = arxvalid(:, fitness.idx(end)); 1464 out.solutions.recentworst.f = fitness.raw(end); 1465 out.solutions.recentworst.evals = counteval + fitness.idx(end) - lambda; 1466 if fitness.hist(1) < out.solutions.bestever.f 1467 out.solutions.bestever.x = arxvalid(:, fitness.idx(1)); 1468 out.solutions.bestever.f = fitness.hist(1); 1469 out.solutions.bestever.evals = counteval + fitness.idx(1) - lambda; 1470 bestever = out.solutions.bestever; 1471 end 1472 1473 % Set stop flag 1474 if fitness.raw(1) <= stopFitness, stopflag(end+1) = {'fitness'}; end 1475 if counteval >= stopMaxFunEvals, stopflag(end+1) = {'maxfunevals'}; end 1476 if countiter >= stopMaxIter, stopflag(end+1) = {'maxiter'}; end 1477 if all(sigma*(max(abs(pc), sqrt(diagC))) < stopTolX) 1478 stopflag(end+1) = {'tolx'}; 1479 end 1480 if any(sigma*sqrt(diagC) > stopTolUpX) 1481 stopflag(end+1) = {'tolupx'}; 1482 end 1483 if sigma*max(diagD) == 0 % should never happen 1484 stopflag(end+1) = {'bug'}; 1485 end 1486 if countiter > 2 && myrange([fitness.sel fitness.hist]) <= stopTolFun 1487 stopflag(end+1) = {'tolfun'}; 1488 end 1489 if countiter >= length(fitness.hist) && myrange(fitness.hist) <= stopTolHistFun 1490 stopflag(end+1) = {'tolhistfun'}; 1491 end 1492 l = floor(length(fitness.histbest)/3); 1493 if 1 < 2 && stopOnStagnation && ... % leads sometimes early stop on ftablet, fcigtab 1494 countiter > N * (5+100/lambda) && ... 1495 length(fitness.histbest) > 100 && ... 1496 median(fitness.histmedian(1:l)) >= median(fitness.histmedian(end-l:end)) && ... 1497 median(fitness.histbest(1:l)) >= median(fitness.histbest(end-l:end)) 1498 stopflag(end+1) = {'stagnation'}; 1499 end 1500 1501 if counteval >= stopFunEvals || countiter >= stopIter 1502 stopflag(end+1) = {'stoptoresume'}; 1503 if length(stopflag) == 1 && flgsaving == 0 1504 error('To resume later the saving option needs to be set'); 1505 end 1506 end 1507 % read stopping message from file signals.par 1508 if flgreadsignals 1509 fid = fopen('./signals.par', 'rt'); % can be performance critical 1510 while fid > 0 1511 strline = fgetl(fid); %fgets(fid, 300); 1512 if strline < 0 % fgets and fgetl returns -1 at end of file 1513 break 1514 end 1515 % 'stop filename' sets stopflag to manual 1516 str = sscanf(strline, ' %s %s', 2); 1517 if strcmp(str, ['stop' opts.LogFilenamePrefix]) 1518 stopflag(end+1) = {'manual'}; 1519 break 1520 end 1521 % 'skip filename run 3' skips a run, but not the last 1522 str = sscanf(strline, ' %s %s %s', 3); 1523 if strcmp(str, ['skip' opts.LogFilenamePrefix 'run']) 1524 i = strfind(strline, 'run'); 1525 if irun == sscanf(strline(i+3:end), ' %d ', 1) && irun <= myeval(opts.Restarts) 1526 stopflag(end+1) = {'skipped'}; 1527 end 1528 end 1529 end % while, break 1530 if fid > 0 1531 fclose(fid); 1532 clear fid; % prevents strange error under octave 1533 end 1534 end 1535 1536 out.stopflag = stopflag; 1537 1538 % ----- output generation ----- 1539 if verbosemodulo > 0 && isfinite(verbosemodulo) 1540 if countiter == 1 || mod(countiter, 10*verbosemodulo) < 1 1541 disp(['Iterat, #Fevals: Function Value (median,worst) ' ... 1542 '|Axis Ratio|' ... 1543 'idx:Min SD idx:Max SD']); 1544 end 1545 if mod(countiter, verbosemodulo) < 1 ... 1546 || (verbosemodulo > 0 && isfinite(verbosemodulo) && ... 1547 (countiter < 3 || ~isempty(stopflag))) 1548 [minstd, minstdidx] = min(sigma*sqrt(diagC)); 1549 [maxstd, maxstdidx] = max(sigma*sqrt(diagC)); 1550 % format display nicely 1551 disp([repmat(' ',1,4-floor(log10(countiter))) ... 1552 num2str(countiter) ' , ' ... 1553 repmat(' ',1,5-floor(log10(counteval))) ... 1554 num2str(counteval) ' : ' ... 1555 num2str(fitness.hist(1), '%.13e') ... 1556 ' +(' num2str(median(fitness.raw)-fitness.hist(1), '%.0e ') ... 1557 ',' num2str(max(fitness.raw)-fitness.hist(1), '%.0e ') ... 1558 ') | ' ... 1559 num2str(max(diagD)/min(diagD), '%4.2e') ' | ' ... 1560 repmat(' ',1,1-floor(log10(minstdidx))) num2str(minstdidx) ':' ... 1561 num2str(minstd, ' %.1e') ' ' ... 1562 repmat(' ',1,1-floor(log10(maxstdidx))) num2str(maxstdidx) ':' ... 1563 num2str(maxstd, ' %.1e')]); 1564 end 1565 end 1566 1567 % measure time for recording data 1568 if countiter < 3 1569 time.c = 0.05; 1570 time.nonoutput = 0; 1571 time.recording = 0; 1572 time.saving = 0.15; % first saving after 3 seconds of 100 iterations 1573 time.plotting = 0; 1574 elseif countiter > 300 1575 % set backward horizon, must be long enough to cover infrequent plotting etc 1576 % time.c = min(1, time.nonoutput/3 + 1e-9); 1577 time.c = max(1e-5, 0.1/sqrt(countiter)); % mean over all or 1e-5 1578 end 1579 % get average time per iteration 1580 time.t1 = clock; 1581 time.act = max(0,etime(time.t1, time.t0)); 1582 time.nonoutput = (1-time.c) * time.nonoutput ... 1583 + time.c * time.act; 1584 1585 time.recording = (1-time.c) * time.recording; % per iteration 1586 time.saving = (1-time.c) * time.saving; 1587 time.plotting = (1-time.c) * time.plotting; 1588 1589 % record output data, concerning time issues 1590 if savemodulo && savetime && (countiter < 1e2 || ~isempty(stopflag) || ... 1591 countiter >= outiter + savemodulo) 1592 outiter = countiter; 1593 % Save output data to files 1594 for namecell = filenames(:)' 1595 name = namecell{:}; 1596 1597 [fid, err] = fopen(['./' filenameprefix name '.dat'], 'a'); 1598 if fid < 1 % err ~= 0 1599 warning(['could not open ' filenameprefix name '.dat']); 1600 else 1601 if strcmp(name, 'axlen') 1602 fprintf(fid, '%d %d %e %e %e ', countiter, counteval, sigma, ... 1603 max(diagD), min(diagD)); 1604 fprintf(fid, '%e ', sort(diagD)); 1605 fprintf(fid, '\n'); 1606 elseif strcmp(name, 'disp') % TODO 1607 elseif strcmp(name, 'fit') 1608 fprintf(fid, '%ld %ld %e %e %25.18e %25.18e %25.18e %25.18e', ... 1609 countiter, counteval, sigma, max(diagD)/min(diagD), ... 1610 out.solutions.bestever.f, ... 1611 fitness.raw(1), median(fitness.raw), fitness.raw(end)); 1612 if ~isempty(varargin) && length(varargin{1}) == 1 && isnumeric(varargin{1}) && varargin{1} ~= 0 1613 fprintf(fid, ' %f', varargin{1}); 1614 end 1615 fprintf(fid, '\n'); 1616 elseif strcmp(name, 'stddev') 1617 fprintf(fid, '%ld %ld %e 0 0 ', countiter, counteval, sigma); 1618 fprintf(fid, '%e ', sigma*sqrt(diagC)); 1619 fprintf(fid, '\n'); 1620 elseif strcmp(name, 'xmean') 1621 if isnan(fmean) 1622 fprintf(fid, '%ld %ld 0 0 0 ', countiter, counteval); 1623 else 1624 fprintf(fid, '%ld %ld 0 0 %e ', countiter, counteval, fmean); 1625 end 1626 fprintf(fid, '%e ', xmean); 1627 fprintf(fid, '\n'); 1628 elseif strcmp(name, 'xrecentbest') 1629 % TODO: fitness is inconsistent with x-value 1630 fprintf(fid, '%ld %ld %25.18e 0 0 ', countiter, counteval, fitness.raw(1)); 1631 fprintf(fid, '%e ', arx(:,fitness.idx(1))); 1632 fprintf(fid, '\n'); 1633 end 1634 fclose(fid); 1635 end 1636 end 1637 1638 % get average time for recording data 1639 time.t2 = clock; 1640 time.recording = time.recording + time.c * max(0,etime(time.t2, time.t1)); 1641 1642 % plot 1643 if flgplotting && countiter > 1 1644 if countiter == 2 1645 iterplotted = 0; 1646 end 1647 if ~isempty(stopflag) || ... 1648 ((time.nonoutput+time.recording) * (countiter - iterplotted) > 1 && ... 1649 time.plotting < 0.05 * (time.nonoutput+time.recording)) 1650 local_plotcmaesdat(324, filenameprefix); 1651 iterplotted = countiter; 1652 % outplot(out); % outplot defined below 1653 if time.plotting == 0 % disregard opening of the window 1654 time.plotting = time.nonoutput+time.recording; 1655 else 1656 time.plotting = time.plotting + time.c * max(0,etime(clock, time.t2)); 1657 end 1658 end 1659 end 1660 if countiter > 100 + 20 && savemodulo && ... 1661 time.recording * countiter > 0.1 && ... % absolute time larger 0.1 second 1662 time.recording > savetime * (time.nonoutput+time.recording) / 100 1663 savemodulo = floor(1.1 * savemodulo) + 1; 1664 % disp('++savemodulo'); %qqq 1665 end 1666 end % if output 1667 1668 % save everything 1669 time.t3 = clock; 1670 if ~isempty(stopflag) || time.saving < 0.05 * time.nonoutput || countiter == 100 1671 xmin = arxvalid(:, fitness.idx(1)); 1672 fmin = fitness.raw(1); 1673 if flgsaving && countiter > 2 1674 clear idx; % prevents error under octave 1675 % -v6 : non-compressed non-unicode for version 6 and earlier 1676 if ~isempty(strsaving) && ~isoctave 1677 save('-mat', strsaving, opts.SaveFilename); % for inspection and possible restart 1678 else 1679 save('-mat', opts.SaveFilename); % for inspection and possible restart 1680 end 1681 time.saving = time.saving + time.c * max(0,etime(clock, time.t3)); 1682 end 1683 end 1684 time.t0 = clock; 1685 1686 % ----- end output generation ----- 1687 1688 end % while, end generation loop 1689 1690 % -------------------- Final Procedures ------------------------------- 1691 1692 % Evaluate xmean and return best recent point in xmin 1693 fmin = fitness.raw(1); 1694 xmin = arxvalid(:, fitness.idx(1)); % Return best point of last generation. 1695 if length(stopflag) > sum(strcmp(stopflag, 'stoptoresume')) % final stopping 1696 out.solutions.mean.f = ... 1697 feval(fitfun, xintobounds(xmean, lbounds, ubounds), varargin{:}); 1698 counteval = counteval + 1; 1699 out.solutions.mean.evals = counteval; 1700 if out.solutions.mean.f < fitness.raw(1) 1701 fmin = out.solutions.mean.f; 1702 xmin = xintobounds(xmean, lbounds, ubounds); % Return xmean as best point 1703 end 1704 if out.solutions.mean.f < out.solutions.bestever.f 1705 out.solutions.bestever = out.solutions.mean; % Return xmean as bestever point 1706 out.solutions.bestever.x = xintobounds(xmean, lbounds, ubounds); 1707 bestever = out.solutions.bestever; 1708 end 1709 end 1710 1711 % Save everything and display final message 1712 if flgsavingfinal 1713 clear idx; % prevents error under octave 1714 if ~isempty(strsaving) && ~isoctave 1715 save('-mat', strsaving, opts.SaveFilename); % for inspection and possible restart 1716 else 1717 save('-mat', opts.SaveFilename); % for inspection and possible restart 1718 end 1719 message = [' (saved to ' opts.SaveFilename ')']; 1720 else 1721 message = []; 1722 end 1723 1724 if flgdisplay 1725 disp(['#Fevals: f(returned x) | bestever.f | stopflag' ... 1726 message]); 1727 if isoctave 1728 strstop = stopflag(:); 1729 else 1730 strcat(stopflag(:), '.'); 1731 end 1732 strstop = stopflag(:); %strcat(stopflag(:), '.'); 1733 disp([repmat(' ',1,6-floor(log10(counteval))) ... 1734 num2str(counteval, '%6.0f') ': ' num2str(fmin, '%.11e') ' | ' ... 1735 num2str(out.solutions.bestever.f, '%.11e') ' | ' ... 1736 strstop{1:end}]); 1737 if N < 102 1738 disp(['mean solution:' sprintf(' %+.1e', xmean)]); 1739 disp(['std deviation:' sprintf(' %.1e', sigma*sqrt(diagC))]); 1740 disp(sprintf('use plotcmaesdat.m for plotting the output at any time (option LogModulo must not be zero)')); 1741 end 1742 if exist('sfile', 'var') 1743 disp(['Results saved in ' sfile]); 1744 end 1745 end 1746 1747 out.arstopflags{irun} = stopflag; 1748 if any(strcmp(stopflag, 'fitness')) ... 1749 || any(strcmp(stopflag, 'maxfunevals')) ... 1750 || any(strcmp(stopflag, 'stoptoresume')) ... 1751 || any(strcmp(stopflag, 'manual')) 1752 break 1753 end 1754end % while irun <= Restarts 1755 1756% --------------------------------------------------------------- 1757% --------------------------------------------------------------- 1758function [x, idx] = xintobounds(x, lbounds, ubounds) 1759% 1760% x can be a column vector or a matrix consisting of column vectors 1761% 1762if ~isempty(lbounds) 1763 if length(lbounds) == 1 1764 idx = x < lbounds; 1765 x(idx) = lbounds; 1766 else 1767 arbounds = repmat(lbounds, 1, size(x,2)); 1768 idx = x < arbounds; 1769 x(idx) = arbounds(idx); 1770 end 1771else 1772 idx = 0; 1773end 1774if ~isempty(ubounds) 1775 if length(ubounds) == 1 1776 idx2 = x > ubounds; 1777 x(idx2) = ubounds; 1778 else 1779 arbounds = repmat(ubounds, 1, size(x,2)); 1780 idx2 = x > arbounds; 1781 x(idx2) = arbounds(idx2); 1782 end 1783else 1784 idx2 = 0; 1785end 1786idx = idx2-idx; 1787 1788% --------------------------------------------------------------- 1789% --------------------------------------------------------------- 1790function opts=getoptions(inopts, defopts) 1791% OPTS = GETOPTIONS(INOPTS, DEFOPTS) handles an arbitrary number of 1792% optional arguments to a function. The given arguments are collected 1793% in the struct INOPTS. GETOPTIONS matches INOPTS with a default 1794% options struct DEFOPTS and returns the merge OPTS. Empty or missing 1795% fields in INOPTS invoke the default value. Fieldnames in INOPTS can 1796% be abbreviated. 1797% 1798% The returned struct OPTS is first assigned to DEFOPTS. Then any 1799% field value in OPTS is replaced by the respective field value of 1800% INOPTS if (1) the field unambiguously (case-insensitive) matches 1801% with the fieldname in INOPTS (cut down to the length of the INOPTS 1802% fieldname) and (2) the field is not empty. 1803% 1804% Example: 1805% In the source-code of the function that needs optional 1806% arguments, the last argument is the struct of optional 1807% arguments: 1808% 1809% function results = myfunction(mandatory_arg, inopts) 1810% % Define four default options 1811% defopts.PopulationSize = 200; 1812% defopts.ParentNumber = 50; 1813% defopts.MaxIterations = 1e6; 1814% defopts.MaxSigma = 1; 1815% 1816% % merge default options with input options 1817% opts = getoptions(inopts, defopts); 1818% 1819% % Thats it! From now on the values in opts can be used 1820% for i = 1:opts.PopulationSize 1821% % do whatever 1822% if sigma > opts.MaxSigma 1823% % do whatever 1824% end 1825% end 1826% 1827% For calling the function myfunction with default options: 1828% myfunction(argument1, []); 1829% For calling the function myfunction with modified options: 1830% opt.pop = 100; % redefine PopulationSize option 1831% opt.PAR = 10; % redefine ParentNumber option 1832% opt.maxiter = 2; % opt.max=2 is ambiguous and would result in an error 1833% myfunction(argument1, opt); 1834 1835% 1836% 04/07/19: Entries can be structs itself leading to a recursive 1837% call to getoptions. 1838% 1839 1840if nargin < 2 || isempty(defopts) % no default options available 1841 opts=inopts; 1842 return 1843elseif isempty(inopts) % empty inopts invoke default options 1844 opts = defopts; 1845 return 1846elseif ~isstruct(defopts) % handle a single option value 1847 if isempty(inopts) 1848 opts = defopts; 1849 elseif ~isstruct(inopts) 1850 opts = inopts; 1851 else 1852 error('Input options are a struct, while default options are not'); 1853 end 1854 return 1855elseif ~isstruct(inopts) % no valid input options 1856 error('The options need to be a struct or empty'); 1857end 1858 1859opts = defopts; % start from defopts 1860 % if necessary overwrite opts fields by inopts values 1861defnames = fieldnames(defopts); 1862idxmatched = []; % indices of defopts that already matched 1863for name = fieldnames(inopts)' 1864 name = name{1}; % name of i-th inopts-field 1865 if isoctave 1866 for i = 1:size(defnames, 1) 1867 idx(i) = strncmp(lower(defnames(i)), lower(name), length(name)); 1868 end 1869 else 1870 idx = strncmpi(defnames, name, length(name)); 1871 end 1872 if sum(idx) > 1 1873 error(['option "' name '" is not an unambigous abbreviation. ' ... 1874 'Use opts=RMFIELD(opts, ''' name, ... 1875 ''') to remove the field from the struct.']); 1876 end 1877 if sum(idx) == 1 1878 defname = defnames{find(idx)}; 1879 if ismember(find(idx), idxmatched) 1880 error(['input options match more than ones with "' ... 1881 defname '". ' ... 1882 'Use opts=RMFIELD(opts, ''' name, ... 1883 ''') to remove the field from the struct.']); 1884 end 1885 idxmatched = [idxmatched find(idx)]; 1886 val = getfield(inopts, name); 1887 % next line can replace previous line from MATLAB version 6.5.0 on and in octave 1888 % val = inopts.(name); 1889 if isstruct(val) % valid syntax only from version 6.5.0 1890 opts = setfield(opts, defname, ... 1891 getoptions(val, getfield(defopts, defname))); 1892 elseif isstruct(getfield(defopts, defname)) 1893 % next three lines can replace previous three lines from MATLAB 1894 % version 6.5.0 on 1895 % opts.(defname) = ... 1896 % getoptions(val, defopts.(defname)); 1897 % elseif isstruct(defopts.(defname)) 1898 warning(['option "' name '" disregarded (must be struct)']); 1899 elseif ~isempty(val) % empty value: do nothing, i.e. stick to default 1900 opts = setfield(opts, defnames{find(idx)}, val); 1901 % next line can replace previous line from MATLAB version 6.5.0 on 1902 % opts.(defname) = inopts.(name); 1903 end 1904 else 1905 warning(['option "' name '" disregarded (unknown field name)']); 1906 end 1907end 1908 1909% --------------------------------------------------------------- 1910% --------------------------------------------------------------- 1911function res=myeval(s) 1912if ischar(s) 1913 res = evalin('caller', s); 1914else 1915 res = s; 1916end 1917 1918% --------------------------------------------------------------- 1919% --------------------------------------------------------------- 1920function res=myevalbool(s) 1921if ~ischar(s) % s may not and cannot be empty 1922 res = s; 1923else % evaluation string s 1924 if strncmpi(lower(s), 'yes', 3) || strncmpi(s, 'on', 2) ... 1925 || strncmpi(s, 'true', 4) || strncmp(s, '1 ', 2) 1926 res = 1; 1927 elseif strncmpi(s, 'no', 2) || strncmpi(s, 'off', 3) ... 1928 || strncmpi(s, 'false', 5) || strncmp(s, '0 ', 2) 1929 res = 0; 1930 else 1931 try res = evalin('caller', s); catch 1932 error(['String value "' s '" cannot be evaluated']); 1933 end 1934 try res ~= 0; catch 1935 error(['String value "' s '" cannot be evaluated reasonably']); 1936 end 1937 end 1938end 1939 1940 1941% --------------------------------------------------------------- 1942% --------------------------------------------------------------- 1943function res = isoctave 1944% any hack to find out whether we are running octave 1945s = version; 1946res = 0; 1947if exist('fflush', 'builtin') && eval(s(1)) < 7 1948 res = 1; 1949end 1950 1951% --------------------------------------------------------------- 1952% --------------------------------------------------------------- 1953function flush 1954if isoctave 1955 feval('fflush', stdout); 1956end 1957 1958% --------------------------------------------------------------- 1959% --------------------------------------------------------------- 1960% ----- replacements for statistic toolbox functions ------------ 1961% --------------------------------------------------------------- 1962% --------------------------------------------------------------- 1963function res=myrange(x) 1964res = max(x) - min(x); 1965 1966% --------------------------------------------------------------- 1967% --------------------------------------------------------------- 1968function res = myprctile(inar, perc, idx) 1969% 1970% Computes the percentiles in vector perc from vector inar 1971% returns vector with length(res)==length(perc) 1972% idx: optional index-array indicating sorted order 1973% 1974 1975N = length(inar); 1976flgtranspose = 0; 1977 1978% sizes 1979if size(perc,1) > 1 1980 perc = perc'; 1981 flgtranspose = 1; 1982 if size(perc,1) > 1 1983 error('perc must not be a matrix'); 1984 end 1985end 1986if size(inar, 1) > 1 && size(inar,2) > 1 1987 error('data inar must not be a matrix'); 1988end 1989 1990% sort inar 1991if nargin < 3 || isempty(idx) 1992 [sar, idx] = sort(inar); 1993else 1994 sar = inar(idx); 1995end 1996 1997res = []; 1998for p = perc 1999 if p <= 100*(0.5/N) 2000 res(end+1) = sar(1); 2001 elseif p >= 100*((N-0.5)/N) 2002 res(end+1) = sar(N); 2003 else 2004 % find largest index smaller than required percentile 2005 availablepercentiles = 100*((1:N)-0.5)/N; 2006 i = max(find(p > availablepercentiles)); 2007 % interpolate linearly 2008 res(end+1) = sar(i) ... 2009 + (sar(i+1)-sar(i))*(p - availablepercentiles(i)) ... 2010 / (availablepercentiles(i+1) - availablepercentiles(i)); 2011 2012 end 2013end 2014 2015if flgtranspose 2016 res = res'; 2017end 2018 2019 2020% --------------------------------------------------------------- 2021% --------------------------------------------------------------- 2022% --------------------------------------------------------------- 2023% --------------------------------------------------------------- 2024 2025 2026 2027function [s, ranks, rankDelta] = local_noisemeasurement(arf1, arf2, lamreev, theta, cutlimit) 2028% function [s ranks rankDelta] = noisemeasurement(arf1, arf2, lamreev, theta) 2029% 2030% Input: 2031% arf1, arf2 : two arrays of function values. arf1 is of size 1xlambda, 2032% arf2 may be of size 1xlamreev or 1xlambda. The first lamreev values 2033% in arf2 are (re-)evaluations of the respective solutions, i.e. 2034% arf1(1) and arf2(1) are two evaluations of "the first" solution. 2035% lamreev: number of reevaluated individuals in arf2 2036% theta : parameter theta for the rank change limit, between 0 and 1, 2037% typically between 0.2 and 0.7. 2038% cutlimit (optional): output s is computed as a mean of rankchange minus 2039% threshold, where rankchange is <=2*(lambda-1). cutlimit limits 2040% abs(rankchange minus threshold) in this calculation to cutlimit. 2041% cutlimit=1 evaluates basically the sign only. cutlimit=2 could be 2042% the rank change with one solution (both evaluations of it). 2043% 2044% Output: 2045% s : noise measurement, s>0 means the noise measure is above the 2046% acceptance threshold 2047% ranks : 2xlambda array, corresponding to [arf1; arf2], of ranks 2048% of arf1 and arf2 in the set [arf1 arf2], values are in [1:2*lambda] 2049% rankDelta: 1xlambda array of rank movements of arf2 compared to 2050% arf1. rankDelta(i) agrees with the number of values from 2051% the set [arf1 arf2] that lie between arf1(i) and arf2(i). 2052% 2053% Note: equal function values might lead to somewhat spurious results. 2054% For this case a revision is advisable. 2055 2056%%% verify input argument sizes 2057if size(arf1,1) ~= 1 2058 error('arf1 must be an 1xlambda array'); 2059elseif size(arf2,1) ~= 1 2060 error('arf2 must be an 1xsomething array'); 2061elseif size(arf1,2) < size(arf2,2) % not really necessary, but saver 2062 error('arf2 must not be smaller than arf1 in length'); 2063end 2064lam = size(arf1, 2); 2065if size(arf1,2) ~= size(arf2,2) 2066 arf2(end+1:lam) = arf1((size(arf2,2)+1):lam); 2067end 2068if nargin < 5 2069 cutlimit = inf; 2070end 2071 2072%%% capture unusual values 2073if any(diff(arf1) == 0) 2074 % this will presumably interpreted as rank change, because 2075 % sort(ones(...)) returns 1,2,3,... 2076 warning([num2str(sum(diff(arf1)==0)) ' equal function values']); 2077end 2078 2079%%% compute rank changes into rankDelta 2080% compute ranks 2081[ignore, idx] = sort([arf1 arf2]); 2082[ignore, ranks] = sort(idx); 2083ranks = reshape(ranks, lam, 2)'; 2084 2085rankDelta = ranks(1,:) - ranks(2,:) - sign(ranks(1,:) - ranks(2,:)); 2086 2087%%% compute rank change limits using both ranks(1,...) and ranks(2,...) 2088for i = 1:lamreev 2089 sumlim(i) = ... 2090 0.5 * (... 2091 myprctile(abs((1:2*lam-1) - (ranks(1,i) - (ranks(1,i)>ranks(2,i)))), ... 2092 theta*50) ... 2093 + myprctile(abs((1:2*lam-1) - (ranks(2,i) - (ranks(2,i)>ranks(1,i)))), ... 2094 theta*50)); 2095end 2096 2097%%% compute measurement 2098%s = abs(rankDelta(1:lamreev)) - sumlim; % lives roughly in 0..2*lambda 2099 2100% max: 1 rankchange in 2*lambda is always fine 2101s = abs(rankDelta(1:lamreev)) - max(1, sumlim); % lives roughly in 0..2*lambda 2102 2103% cut-off limit 2104idx = abs(s) > cutlimit; 2105s(idx) = sign(s(idx)) * cutlimit; 2106s = mean(s); 2107 2108% --------------------------------------------------------------- 2109% --------------------------------------------------------------- 2110% --------------------------------------------------------------- 2111% --------------------------------------------------------------- 2112% just a "local" copy of plotcmaesdat.m, with manual_mode set to zero 2113function local_plotcmaesdat(figNb, filenameprefix, filenameextension, objectvarname) 2114% PLOTCMAESDAT; 2115% PLOTCMAES(FIGURENUMBER_iBEGIN_iEND, FILENAMEPREFIX, FILENAMEEXTENSION, OBJECTVARNAME); 2116% plots output from CMA-ES, e.g. cmaes.m, Java class CMAEvolutionStrategy... 2117% mod(figNb,100)==1 plots versus iterations. 2118% 2119% PLOTCMAES([101 300]) plots versus iteration, from iteration 300. 2120% PLOTCMAES([100 150 800]) plots versus function evaluations, between iteration 150 and 800. 2121% 2122% Upper left subplot: blue/red: function value of the best solution in the 2123% recent population, cyan: same function value minus best 2124% ever seen function value, green: sigma, red: ratio between 2125% longest and shortest principal axis length which is equivalent 2126% to sqrt(cond(C)). 2127% Upper right plot: time evolution of the distribution mean (default) or 2128% the recent best solution vector. 2129% Lower left: principal axes lengths of the distribution ellipsoid, 2130% equivalent with the sqrt(eig(C)) square root eigenvalues of C. 2131% Lower right: magenta: minimal and maximal "true" standard deviation 2132% (with sigma included) in the coordinates, other colors: sqrt(diag(C)) 2133% of all diagonal elements of C, if C is diagonal they equal to the 2134% lower left. 2135% 2136% Files [FILENAMEPREFIX name FILENAMEEXTENSION] are used, where 2137% name = axlen, OBJECTVARNAME (xmean|xrecentbest), fit, or stddev. 2138% 2139 2140manual_mode = 0; 2141 2142if nargin < 1 || isempty(figNb) 2143 figNb = 325; 2144end 2145if nargin < 2 || isempty(filenameprefix) 2146 filenameprefix = 'outcmaes'; 2147end 2148if nargin < 3 || isempty(filenameextension) 2149 filenameextension = '.dat'; 2150end 2151if nargin < 4 || isempty(objectvarname) 2152 objectvarname = 'xmean'; 2153 objectvarname = 'xrecentbest'; 2154end 2155% load data 2156d.x = load([filenameprefix objectvarname filenameextension]); 2157% d.x = load([filenameprefix 'xmean' filenameextension]); 2158% d.x = load([filenameprefix 'xrecentbest' filenameextension]); 2159d.f = load([filenameprefix 'fit' filenameextension]); 2160d.std = load([filenameprefix 'stddev' filenameextension]); 2161d.D = load([filenameprefix 'axlen' filenameextension]); 2162 2163% interpret entries in figNb for cutting out some data 2164if length(figNb) > 1 2165 iend = inf; 2166 istart = figNb(2); 2167 if length(figNb) > 2 2168 iend = figNb(3); 2169 end 2170 figNb = figNb(1); 2171 d.x = d.x(d.x(:,1) >= istart & d.x(:,1) <= iend, :); 2172 d.f = d.f(d.f(:,1) >= istart & d.f(:,1) <= iend, :); 2173 d.std = d.std(d.std(:,1) >= istart & d.std(:,1) <= iend, :); 2174 d.D = d.D(d.D(:,1) >= istart & d.D(:,1) <= iend, :); 2175end 2176 2177% decide for x-axis 2178iabscissa = 2; % 1== versus iterations, 2==versus fevals 2179if mod(figNb,100) == 1 2180 iabscissa = 1; % a short hack 2181end 2182if iabscissa == 1 2183 xlab ='iterations'; 2184elseif iabscissa == 2 2185 xlab = 'function evaluations'; 2186end 2187 2188if size(d.x, 2) < 1000 2189 minxend = 1.03*d.x(end, iabscissa); 2190else 2191 minxend = 0; 2192end 2193 2194% set up figure window 2195if manual_mode 2196 figure(figNb); % just create and raise the figure window 2197else 2198 if 1 < 3 && evalin('caller', 'iterplotted') == 0 && evalin('caller', 'irun') == 1 2199 figure(figNb); % incomment this, if figure raise in the beginning is not desired 2200 elseif ismember(figNb, findobj('Type', 'figure')) 2201 set(0, 'CurrentFigure', figNb); % prevents raise of existing figure window 2202 else 2203 figure(figNb); 2204 end 2205end 2206 2207% plot fitness etc 2208foffset = 1e-99; 2209dfit = d.f(:,6)-min(d.f(:,6)); 2210[ignore, idxbest] = min(dfit); 2211dfit(dfit<1e-98) = NaN; 2212subplot(2,2,1); hold off; 2213dd = abs(d.f(:,7:8)) + foffset; 2214dd(d.f(:,7:8)==0) = NaN; 2215semilogy(d.f(:,iabscissa), dd, '-k'); hold on; 2216% additional fitness data, for example constraints values 2217if size(d.f,2) > 8 2218 dd = abs(d.f(:,9:end)) + 10*foffset; % a hack 2219 % dd(d.f(:,9:end)==0) = NaN; 2220 semilogy(d.f(:,iabscissa), dd, '-m'); hold on; 2221 if size(d.f,2) > 12 2222 semilogy(d.f(:,iabscissa),abs(d.f(:,[7 8 11 13]))+foffset,'-k'); hold on; 2223 end 2224end 2225 2226idx = find(d.f(:,6)>1e-98); % positive values 2227if ~isempty(idx) % otherwise non-log plot gets hold 2228 semilogy(d.f(idx,iabscissa), d.f(idx,6)+foffset, '.b'); hold on; 2229end 2230idx = find(d.f(:,6) < -1e-98); % negative values 2231if ~isempty(idx) 2232 semilogy(d.f(idx, iabscissa), abs(d.f(idx,6))+foffset,'.r'); hold on; 2233end 2234semilogy(d.f(:,iabscissa),abs(d.f(:,6))+foffset,'-b'); hold on; 2235semilogy(d.f(:,iabscissa),dfit,'-c'); hold on; 2236semilogy(d.f(:,iabscissa),(d.f(:,4)),'-r'); hold on; % AR 2237semilogy(d.std(:,iabscissa), [max(d.std(:,6:end)')' ... 2238 min(d.std(:,6:end)')'], '-m'); % max,min std 2239maxval = max(d.std(end,6:end)); 2240minval = min(d.std(end,6:end)); 2241text(d.std(end,iabscissa), maxval, sprintf('%.0e', maxval)); 2242text(d.std(end,iabscissa), minval, sprintf('%.0e', minval)); 2243 2244semilogy(d.std(:,iabscissa),(d.std(:,3)),'-g'); % sigma 2245 % plot best f 2246semilogy(d.f(idxbest,iabscissa),min(dfit),'*c'); hold on; 2247semilogy(d.f(idxbest,iabscissa),abs(d.f(idxbest,6))+foffset,'*r'); hold on; 2248 2249ax = axis; 2250ax(2) = max(minxend, ax(2)); 2251axis(ax); 2252 2253yannote = 10^(log10(ax(3)) + 0.05*(log10(ax(4))-log10(ax(3)))); 2254text(ax(1), yannote, ... 2255 [ 'f=' num2str(d.f(end,6), '%.15g') ]); 2256 2257title('blue:abs(f), cyan:f-min(f), green:sigma, red:axis ratio'); 2258grid on; 2259 2260subplot(2,2,2); hold off; 2261plot(d.x(:,iabscissa), d.x(:,6:end),'-'); hold on; 2262ax = axis; 2263ax(2) = max(minxend, ax(2)); 2264axis(ax); 2265 2266% add some annotation lines 2267[ignore idx] = sort(d.x(end,6:end)); 2268% choose no more than 25 indices 2269idxs = round(linspace(1, size(d.x,2)-5, min(size(d.x,2)-5, 25))); 2270yy = repmat(NaN, 2, size(d.x,2)-5); 2271yy(1,:) = d.x(end, 6:end); 2272yy(2,idx(idxs)) = linspace(ax(3), ax(4), length(idxs)); 2273plot([d.x(end,iabscissa) ax(2)], yy, '-'); 2274plot(repmat(d.x(end,iabscissa),2), [ax(3) ax(4)], 'k-'); 2275for i = idx(idxs) 2276 text(ax(2), yy(2,i), ... 2277 ['x(' num2str(i) ')=' num2str(yy(1,i))]); 2278end 2279 2280lam = 'NA'; 2281if size(d.x, 1) > 1 && d.x(end, 1) > d.x(end-1, 1) 2282 lam = num2str((d.x(end, 2) - d.x(end-1, 2)) / (d.x(end, 1) - d.x(end-1, 1))); 2283end 2284title(['Object Variables (' num2str(size(d.x, 2)-5) ... 2285 '-D, popsize~' lam ')']);grid on; 2286 2287subplot(2,2,3); hold off; semilogy(d.D(:,iabscissa), d.D(:,6:end), '-'); 2288ax = axis; 2289ax(2) = max(minxend, ax(2)); 2290axis(ax); 2291title('Principal Axes Lengths');grid on; 2292xlabel(xlab); 2293 2294subplot(2,2,4); hold off; 2295% semilogy(d.std(:,iabscissa), d.std(:,6:end), 'k-'); hold on; 2296% remove sigma from stds 2297d.std(:,6:end) = d.std(:,6:end) ./ (d.std(:,3) * ones(1,size(d.std,2)-5)); 2298semilogy(d.std(:,iabscissa), d.std(:,6:end), '-'); hold on; 2299if 11 < 3 % max and min std 2300 semilogy(d.std(:,iabscissa), [d.std(:,3).*max(d.std(:,6:end)')' ... 2301 d.std(:,3).*min(d.std(:,6:end)')'], '-m', 'linewidth', 2); 2302 maxval = max(d.std(end,6:end)); 2303 minval = min(d.std(end,6:end)); 2304 text(d.std(end,iabscissa), d.std(end,3)*maxval, sprintf('max=%.0e', maxval)); 2305 text(d.std(end,iabscissa), d.std(end,3)*minval, sprintf('min=%.0e', minval)); 2306end 2307ax = axis; 2308ax(2) = max(minxend, ax(2)); 2309axis(ax); 2310% add some annotation lines 2311[ignore, idx] = sort(d.std(end,6:end)); 2312% choose no more than 25 indices 2313idxs = round(linspace(1, size(d.x,2)-5, min(size(d.x,2)-5, 25))); 2314yy = repmat(NaN, 2, size(d.std,2)-5); 2315yy(1,:) = d.std(end, 6:end); 2316yy(2,idx(idxs)) = logspace(log10(ax(3)), log10(ax(4)), length(idxs)); 2317semilogy([d.std(end,iabscissa) ax(2)], yy, '-'); 2318semilogy(repmat(d.std(end,iabscissa),2), [ax(3) ax(4)], 'k-'); 2319for i = idx(idxs) 2320 text(ax(2), yy(2,i), [' ' num2str(i)]); 2321end 2322title('Standard Deviations in Coordinates divided by sigma');grid on; 2323xlabel(xlab); 2324 2325if figNb ~= 324 2326 % zoom on; % does not work in Octave 2327end 2328drawnow; 2329 2330% --------------------------------------------------------------- 2331% --------------- TEST OBJECTIVE FUNCTIONS ---------------------- 2332% --------------------------------------------------------------- 2333 2334%%% Unimodal functions 2335 2336function f=fjens1(x) 2337% 2338% use population size about 2*N 2339% 2340f = sum((x>0) .* x.^1, 1); 2341if any(any(x<0)) 2342 idx = sum(x < 0, 1) > 0; 2343 f(idx) = 1e3; 2344 % f = f + 1e3 * sum(x<0, 1); 2345 % f = f + 10 * sum((x<0) .* x.^2, 1); 2346 f(idx) = f(idx) + 1e-3*abs(randn(1,sum(idx))); 2347 % f(idx) = NaN; 2348end 2349 2350function f=fsphere(x) 2351f = sum(x.^2,1); 2352 2353function f=fmax(x) 2354f = max(abs(x), [], 1); 2355 2356function f=fssphere(x) 2357f=sqrt(sum(x.^2, 1)); 2358 2359% lb = -0.512; ub = 512; 2360% xfeas = x; 2361% xfeas(x<lb) = lb; 2362% xfeas(x>ub) = ub; 2363% f=sum(xfeas.^2, 1); 2364% f = f + 1e-9 * sum((xfeas-x).^2); 2365 2366function f=fspherenoise(x, Nevals) 2367if nargin < 2 || isempty(Nevals) 2368 Nevals = 1; 2369end 2370[N,popsi] = size(x); 2371% x = x .* (1 + 0.3e-0 * randn(N, popsi)/(2*N)); % actuator noise 2372fsum = 10.^(0*(0:N-1)/(N-1)) * x.^2; 2373% f = 0*rand(1,1) ... 2374% + fsum ... 2375% + fsum .* (2*randn(1,popsi) ./ randn(1,popsi).^0 / (2*N)) ... 2376% + 1*fsum.^0.9 .* 2*randn(1,popsi) / (2*N); % 2377 2378% f = fsum .* exp(0.1*randn(1,popsi)); 2379f = fsum .* (1 + (10/(N+10)/sqrt(Nevals))*randn(1,popsi)); 2380% f = fsum .* (1 + (0.1/N)*randn(1,popsi)./randn(1,popsi).^1); 2381 2382idx = rand(1,popsi) < 0.0; 2383if sum(idx) > 0 2384 f(idx) = f(idx) + 1e3*exp(randn(1,sum(idx))); 2385end 2386 2387function f=fmixranks(x) 2388N = size(x,1); 2389f=(10.^(0*(0:(N-1))/(N-1))*x.^2).^0.5; 2390if size(x, 2) > 1 % compute ranks, if it is a population 2391 [ignore, idx] = sort(f); 2392 [ignore, ranks] = sort(idx); 2393 k = 9; % number of solutions randomly permuted, lambda/2-1 2394 % works still quite well (two time slower) 2395 for i = k+1:k-0:size(x,2) 2396 idx(i-k+(1:k)) = idx(i-k+randperm(k)); 2397 end 2398 %disp([ranks' f']) 2399 [ignore, ranks] = sort(idx); 2400 %disp([ranks' f']) 2401 %pause 2402 f = ranks+1e-9*randn(1,1); 2403end 2404 2405function f = fsphereoneax(x) 2406f = x(1)^2; 2407f = mean(x)^2; 2408 2409function f=frandsphere(x) 2410N = size(x,1); 2411idx = ceil(N*rand(7,1)); 2412f=sum(x(idx).^2); 2413 2414function f=fspherelb0(x, M) % lbound at zero for 1:M needed 2415if nargin < 2 M = 0; end 2416N = size(x,1); 2417% M active bounds, f_i = 1 for x = 0 2418f = -M + sum((x(1:M) + 1).^2); 2419f = f + sum(x(M+1:N).^2); 2420 2421function f=fspherehull(x) 2422% Patton, Dexter, Goodman, Punch 2423% in -500..500 2424% spherical ridge through zeros(N,1) 2425% worst case start point seems x = 2*100*sqrt(N) 2426% and small step size 2427N = size(x,1); 2428f = norm(x) + (norm(x-100*sqrt(N)) - 100*N)^2; 2429 2430function f=fellilb0(x, idxM, scal) % lbound at zero for 1:M needed 2431N = size(x,1); 2432if nargin < 3 || isempty(scal) 2433 scal = 100; 2434end 2435scale=scal.^((0:N-1)/(N-1)); 2436if nargin < 2 || isempty(idxM) 2437 idxM = 1:N; 2438end 2439%scale(N) = 1e0; 2440% M active bounds 2441xopt = 0.1; 2442x(idxM) = x(idxM) + xopt; 2443f = scale.^2*x.^2; 2444f = f - sum((xopt*scale(idxM)).^2); 2445% f = exp(f) - 1; 2446% f = log10(f+1e-19) + 19; 2447 2448f = f + 1e-19; 2449 2450function f=fcornersphere(x) 2451w = ones(size(x,1)); 2452w(1) = 2.5; w(2)=2.5; 2453idx = x < 0; 2454f = sum(x(idx).^2); 2455idx = x > 0; 2456f = f + 2^2*sum(w(idx).*x(idx).^2); 2457 2458function f=fsectorsphere(x, scal) 2459% 2460% This is deceptive for cumulative sigma control CSA in large dimension: 2461% The strategy (initially) diverges for N=50 and popsize = 150. (Even 2462% for cs==1 this can be observed for larger settings of N and 2463% popsize.) The reason is obvious from the function topology. 2464% Divergence can be avoided by setting boundaries or adding a 2465% penalty for large ||x||. Then, convergence can be observed again. 2466% Conclusion: for popsize>N cumulative sigma control is not completely 2467% reasonable, but I do not know better alternatives. In particular: 2468% TPA takes longer to converge than CSA when the latter still works. 2469% 2470if nargin < 2 || isempty (scal) 2471 scal = 1e3; 2472end 2473f=sum(x.^2,1); 2474idx = x<0; 2475f = f + (scal^2 - 1) * sum((idx.*x).^2,1); 2476if 11 < 3 2477 idxpen = find(f>1e9); 2478 if ~isempty(idxpen) 2479 f(idxpen) = f(idxpen) + 1e8*sum(x(:,idxpen).^2,1); 2480 end 2481end 2482 2483function f=fstepsphere(x, scal) 2484if nargin < 2 || isempty (scal) 2485 scal = 1e0; 2486end 2487N = size(x,1); 2488f=1e-11+sum(scal.^((0:N-1)/(N-1))*floor(x+0.5).^2); 2489f=1e-11+sum(floor(scal.^((0:N-1)/(N-1))'.*x+0.5).^2); 2490% f=1e-11+sum(floor(x+0.5).^2); 2491 2492function f=fstep(x) 2493% in -5.12..5.12 (bounded) 2494N = size(x,1); 2495f=1e-11+6*N+sum(floor(x)); 2496 2497function f=flnorm(x, scal, e) 2498if nargin < 2 || isempty(scal) 2499 scal = 1; 2500end 2501if nargin < 3 || isempty(e) 2502 e = 1; 2503end 2504if e==inf 2505 f = max(abs(x)); 2506else 2507 N = size(x,1); 2508 scale = scal.^((0:N-1)/(N-1))'; 2509 f=sum(abs(scale.*x).^e); 2510end 2511 2512function f=fneumaier3(x) 2513% in -n^2..n^2 2514% x^*-i = i(n+1-i) 2515N = size(x,1); 2516% f = N*(N+4)*(N-1)/6 + sum((x-1).^2) - sum(x(1:N-1).*x(2:N)); 2517f = sum((x-1).^2) - sum(x(1:N-1).*x(2:N)); 2518 2519function f = fmaxmindist(y) 2520% y in [-1,1], y(1:2) is first point on a plane, y(3:4) second etc 2521% points best 2522% 5 1.4142 2523% 8 1.03527618 2524% 10 0.842535997 2525% 20 0.5997 2526pop = size(y,2); 2527N = size(y,1)/2; 2528f = []; 2529for ipop = 1:pop 2530 if any(abs(y(:,ipop)) > 1) 2531 f(ipop) = NaN; 2532 else 2533 x = reshape(y(:,ipop), [2, N]); 2534 f(ipop) = inf; 2535 for i = 1:N 2536 f(ipop) = min(f(ipop), min(sqrt(sum((x(:,[1:i-1 i+1:N]) - repmat(x(:,i), 1, N-1)).^2, 1)))); 2537 end 2538 end 2539end 2540f = -f; 2541 2542function f=fchangingsphere(x) 2543N = size(x,1); 2544global scale_G; global count_G; if isempty(count_G) count_G=-1; end 2545count_G = count_G+1; 2546if mod(count_G,10) == 0 2547 scale_G = 10.^(2*rand(1,N)); 2548end 2549%disp(scale(1)); 2550f = scale_G*x.^2; 2551 2552function f= flogsphere(x) 2553f = 1-exp(-sum(x.^2)); 2554 2555function f= fexpsphere(x) 2556f = exp(sum(x.^2)) - 1; 2557 2558function f=fbaluja(x) 2559% in [-0.16 0.16] 2560y = x(1); 2561for i = 2:length(x) 2562 y(i) = x(i) + y(i-1); 2563end 2564f = 1e5 - 1/(1e-5 + sum(abs(y))); 2565 2566function f=fschwefel(x) 2567f = 0; 2568for i = 1:size(x,1) 2569 f = f+sum(x(1:i))^2; 2570end 2571 2572function f=fcigar(x, ar) 2573if nargin < 2 || isempty(ar) 2574 ar = 1e3; 2575end 2576f = x(1,:).^2 + ar^2*sum(x(2:end,:).^2,1); 2577 2578function f=fcigtab(x) 2579f = x(1,:).^2 + 1e8*x(end,:).^2 + 1e4*sum(x(2:(end-1),:).^2, 1); 2580 2581function f=ftablet(x) 2582f = 1e6*x(1,:).^2 + sum(x(2:end,:).^2, 1); 2583 2584function f=felli(x, lgscal, expon, expon2) 2585% lgscal: log10(axis ratio) 2586% expon: x_i^expon, sphere==2 2587N = size(x,1); if N < 2 error('dimension must be greater one'); end 2588 2589% x = x - repmat(-0.5+(1:N)',1,size(x,2)); % optimum in 1:N 2590if nargin < 2 || isempty(lgscal), lgscal = 3; end 2591if nargin < 3 || isempty(expon), expon = 2; end 2592if nargin < 4 || isempty(expon2), expon2 = 1; end 2593 2594f=((10^(lgscal*expon)).^((0:N-1)/(N-1)) * abs(x).^expon).^(1/expon2); 2595% if rand(1,1) > 0.015 2596% f = NaN; 2597% end 2598% f = f + randn(size(f)); 2599 2600function f=fellitest(x) 2601beta = 0.9; 2602N = size(x,1); 2603f = (1e6.^((0:(N-1))/(N-1))).^beta * (x.^2).^beta; 2604 2605 2606function f=fellii(x, scal) 2607N = size(x,1); if N < 2 error('dimension must be greater one'); end 2608if nargin < 2 2609 scal = 1; 2610end 2611f= (scal*(1:N)).^2 * (x).^2; 2612 2613function f=fellirot(x) 2614N = size(x,1); 2615global ORTHOGONALCOORSYSTEM_G 2616if isempty(ORTHOGONALCOORSYSTEM_G) ... 2617 || length(ORTHOGONALCOORSYSTEM_G) < N ... 2618 || isempty(ORTHOGONALCOORSYSTEM_G{N}) 2619 coordinatesystem(N); 2620end 2621f = felli(ORTHOGONALCOORSYSTEM_G{N}*x); 2622 2623function f=frot(x, fun, varargin) 2624N = size(x,1); 2625global ORTHOGONALCOORSYSTEM_G 2626if isempty(ORTHOGONALCOORSYSTEM_G) ... 2627 || length(ORTHOGONALCOORSYSTEM_G) < N ... 2628 || isempty(ORTHOGONALCOORSYSTEM_G{N}) 2629 coordinatesystem(N); 2630end 2631f = feval(fun, ORTHOGONALCOORSYSTEM_G{N}*x, varargin{:}); 2632 2633function coordinatesystem(N) 2634if nargin < 1 || isempty(N) 2635 arN = 2:30; 2636else 2637 arN = N; 2638end 2639global ORTHOGONALCOORSYSTEM_G 2640ORTHOGONALCOORSYSTEM_G{1} = 1; 2641for N = arN 2642 ar = randn(N,N); 2643 for i = 1:N 2644 for j = 1:i-1 2645 ar(:,i) = ar(:,i) - ar(:,i)'*ar(:,j) * ar(:,j); 2646 end 2647 ar(:,i) = ar(:,i) / norm(ar(:,i)); 2648 end 2649 ORTHOGONALCOORSYSTEM_G{N} = ar; 2650end 2651 2652function f=fplane(x) 2653f=x(1); 2654 2655function f=ftwoaxes(x) 2656f = sum(x(1:floor(end/2),:).^2, 1) + 1e6*sum(x(floor(1+end/2):end,:).^2, 1); 2657 2658function f=fparabR(x) 2659f = -x(1,:) + 100*sum(x(2:end,:).^2,1); 2660 2661function f=fsharpR(x) 2662f = abs(-x(1, :)).^2 + 100 * sqrt(sum(x(2:end,:).^2, 1)); 2663 2664function f=frosen(x) 2665if size(x,1) < 2, error('dimension must be greater one'); end 2666N = size(x,1); 2667popsi = size(x,2); 2668f = 1e2*sum((x(1:end-1,:).^2 - x(2:end,:)).^2,1) + sum((x(1:end-1,:)-1).^2,1); 2669% f = f + f^0.9 .* (2*randn(1,popsi) ./ randn(1,popsi).^0 / (2*N)); 2670 2671function f=frosenlin(x) 2672if size(x,1) < 2 error('dimension must be greater one'); end 2673 2674x_org = x; 2675x(x>30) = 30; 2676x(x<-30) = -30; 2677 2678f = 1e2*sum(-(x(1:end-1,:).^2 - x(2:end,:)),1) + ... 2679 sum((x(1:end-1,:)-1).^2,1); 2680 2681f = f + sum((x-x_org).^2,1); 2682% f(any(abs(x)>30,1)) = NaN; 2683 2684function f=frosenrot(x) 2685N = size(x,1); 2686global ORTHOGONALCOORSYSTEM_G 2687if isempty(ORTHOGONALCOORSYSTEM_G) ... 2688 || length(ORTHOGONALCOORSYSTEM_G) < N ... 2689 || isempty(ORTHOGONALCOORSYSTEM_G{N}) 2690 coordinatesystem(N); 2691end 2692f = frosen(ORTHOGONALCOORSYSTEM_G{N}*x); 2693 2694function f=frosenmodif(x) 2695f = 74 + 100*(x(2)-x(1)^2)^2 + (1-x(1))^2 ... 2696 - 400*exp(-sum((x+1).^2)/2/0.05); 2697 2698function f=fschwefelrosen1(x) 2699% in [-10 10] 2700f=sum((x.^2-x(1)).^2 + (x-1).^2); 2701 2702function f=fschwefelrosen2(x) 2703% in [-10 10] 2704f=sum((x(2:end).^2-x(1)).^2 + (x(2:end)-1).^2); 2705 2706function f=fdiffpow(x) 2707[N, popsi] = size(x); if N < 2 error('dimension must be greater one'); end 2708 2709f = sum(abs(x).^repmat(2+10*(0:N-1)'/(N-1), 1, popsi), 1); 2710f = sqrt(f); 2711 2712function f=fabsprod(x) 2713f = sum(abs(x),1) + prod(abs(x),1); 2714 2715function f=ffloor(x) 2716f = sum(floor(x+0.5).^2,1); 2717 2718function f=fmaxx(x) 2719f = max(abs(x), [], 1); 2720 2721%%% Multimodal functions 2722 2723function f=fbirastrigin(x) 2724% todo: the volume needs to be a constant 2725N = size(x,1); 2726idx = (sum(x, 1) < 0.5*N); % global optimum 2727f = zeros(1,size(x,2)); 2728f(idx) = 10*(N-sum(cos(2*pi*x(:,idx)),1)) + sum(x(:,idx).^2,1); 2729idx = ~idx; 2730f(idx) = 0.1 + 10*(N-sum(cos(2*pi*(x(:,idx)-2)),1)) + sum((x(:,idx)-2).^2,1); 2731 2732function f=fackley(x) 2733% -32.768..32.768 2734% Adding a penalty outside the interval is recommended, 2735% because for large step sizes, fackley imposes like frand 2736% 2737N = size(x,1); 2738f = 20-20*exp(-0.2*sqrt(sum(x.^2)/N)); 2739f = f + (exp(1) - exp(sum(cos(2*pi*x))/N)); 2740% add penalty outside the search interval 2741f = f + sum((x(x>32.768)-32.768).^2) + sum((x(x<-32.768)+32.768).^2); 2742 2743function f = fbohachevsky(x) 2744% -15..15 2745f = sum(x(1:end-1).^2 + 2 * x(2:end).^2 - 0.3 * cos(3*pi*x(1:end-1)) ... 2746 - 0.4 * cos(4*pi*x(2:end)) + 0.7); 2747 2748function f=fconcentric(x) 2749% in +-600 2750s = sum(x.^2); 2751f = s^0.25 * (sin(50*s^0.1)^2 + 1); 2752 2753function f=fgriewank(x) 2754% in [-600 600] 2755[N, P] = size(x); 2756f = 1 - prod(cos(x'./sqrt(1:N))) + sum(x.^2)/4e3; 2757scale = repmat(sqrt(1:N)', 1, P); 2758f = 1 - prod(cos(x./scale), 1) + sum(x.^2, 1)/4e3; 2759% f = f + 1e4*sum(x(abs(x)>5).^2); 2760% if sum(x(abs(x)>5).^2) > 0 2761% f = 1e4 * sum(x(abs(x)>5).^2) + 1e8 * sum(x(x>5)).^2; 2762% end 2763 2764function f=fgriewrosen(x) 2765% F13 or F8F2 2766[N, P] = size(x); 2767scale = repmat(sqrt(1:N)', 1, P); 2768y = [x(2:end,:); x(1,:)]; 2769x = 100 * (x.^2 - y) + (x - 1).^2; % Rosenbrock part 2770f = 1 - prod(cos(x./scale), 1) + sum(x.^2, 1)/4e3; 2771f = sum(1 - cos(x) + x.^2/4e3, 1); 2772 2773function f=fspallpseudorastrigin(x, scal, skewfac, skewstart, amplitude) 2774% by default multi-modal about between -30 and 30 2775if nargin < 5 || isempty(amplitude) 2776 amplitude = 40; 2777end 2778if nargin < 4 || isempty(skewstart) 2779 skewstart = 0; 2780end 2781if nargin < 3 || isempty(skewfac) 2782 skewfac = 1; 2783end 2784if nargin < 2 || isempty(scal) 2785 scal = 1; 2786end 2787N = size(x,1); 2788scale = 1; 2789if N > 1 2790 scale=scal.^((0:N-1)'/(N-1)); 2791end 2792% simple version: 2793% f = amplitude*(N - sum(cos(2*pi*(scale.*x)))) + sum((scale.*x).^2); 2794 2795% skew version: 2796y = repmat(scale, 1, size(x,2)) .* x; 2797idx = find(x > skewstart); 2798if ~isempty(idx) 2799 y(idx) = skewfac*y(idx); 2800end 2801f = amplitude * (0*N-prod(cos((2*pi)^0*y),1)) + 0.05 * sum(y.^2,1) ... 2802 + randn(1,1); 2803 2804function f=frastrigin(x, scal, skewfac, skewstart, amplitude) 2805% by default multi-modal about between -30 and 30 2806if nargin < 5 || isempty(amplitude) 2807 amplitude = 10; 2808end 2809if nargin < 4 || isempty(skewstart) 2810 skewstart = 0; 2811end 2812if nargin < 3 || isempty(skewfac) 2813 skewfac = 1; 2814end 2815if nargin < 2 || isempty(scal) 2816 scal = 1; 2817end 2818N = size(x,1); 2819scale = 1; 2820if N > 1 2821 scale=scal.^((0:N-1)'/(N-1)); 2822end 2823% simple version: 2824% f = amplitude*(N - sum(cos(2*pi*(scale.*x)))) + sum((scale.*x).^2); 2825 2826% skew version: 2827y = repmat(scale, 1, size(x,2)) .* x; 2828idx = find(x > skewstart); 2829% idx = intersect(idx, 2:2:10); 2830if ~isempty(idx) 2831 y(idx) = skewfac*y(idx); 2832end 2833f = amplitude * (N-sum(cos(2*pi*y),1)) + sum(y.^2,1); 2834 2835function f=frastriginmax(x) 2836N = size(x,1); 2837f = (N/20)*807.06580387678 - (10 * (N-sum(cos(2*pi*x),1)) + sum(x.^2,1)); 2838f(any(abs(x) > 5.12)) = 1e2*N; 2839 2840function f = fschaffer(x) 2841% -100..100 2842N = size(x,1); 2843s = x(1:N-1,:).^2 + x(2:N,:).^2; 2844f = sum(s.^0.25 .* (sin(50*s.^0.1).^2+1), 1); 2845 2846function f=fschwefelmult(x) 2847% -500..500 2848% 2849N = size(x,1); 2850f = - sum(x.*sin(sqrt(abs(x))), 1); 2851f = 418.9829*N - 1.27275661e-5*N - sum(x.*sin(sqrt(abs(x))), 1); 2852% penalty term 2853f = f + 1e4*sum((abs(x)>500) .* (abs(x)-500).^2, 1); 2854 2855function f=ftwomax(x) 2856% Boundaries at +/-5 2857N = size(x,1); 2858f = -abs(sum(x)) + 5*N; 2859 2860function f=ftwomaxtwo(x) 2861% Boundaries at +/-10 2862N = size(x,1); 2863f = abs(sum(x)); 2864if f > 30 2865 f = f - 30; 2866end 2867f = -f; 2868 2869function f=frand(x) 2870f=1./(1-rand(1, size(x,2))) - 1; 2871 2872% CHANGES 2873% 12/02/19: "future" setting of ccum, correcting for large mueff, is default now 2874% 11/11/15: bug-fix: max value for ccovmu_sep setting corrected 2875% 10/11/11: (3.52.beta) boundary handling: replace max with min in change 2876% rate formula. Active CMA: check of pos.def. improved. 2877% Plotting: value of lambda appears in the title. 2878% 10/04/03: (3.51.beta) active CMA cleaned up. Equal fitness detection 2879% looks into history now. 2880% 10/03/08: (3.50.beta) "active CMA" revised and bug-fix of ambiguous 2881% option Noise.alpha -> Noise.alphasigma. 2882% 09/10/12: (3.40.beta) a slightly modified version of "active CMA", 2883% that is a negative covariance matrix update, use option 2884% CMA.active. In 10;30;90-D the gain on ftablet is a factor 2885% of 1.6;2.5;4.4 (the scaling improves by sqrt(N)). On 2886% Rosenbrock the gain is about 25%. On sharp ridge the 2887% behavior is improved. Cigar is unchanged. 2888% 09/08/10: local plotcmaesdat remains in backround 2889% 09/08/10: bug-fix in time management for data writing, logtime was not 2890% considered properly (usually not at all). 2891% 09/07/05: V3.24: stagnation termination added 2892% 08/09/27: V3.23: momentum alignment is out-commented and de-preciated 2893% 08/09/25: V3.22: re-alignment of sigma and C was buggy 2894% 08/07/15: V3.20, CMA-parameters are options now. ccov and mucov were replaced 2895% by ccov1 \approx ccov/mucov and ccovmu \approx (1-1/mucov)*ccov 2896% 08/06/30: file name xrecent was change to xrecentbest (compatible with other 2897% versions) 2898% 08/06/29: time stamp added to output files 2899% 08/06/28: bug fixed with resume option, commentary did not work 2900% 08/06/28: V3.10, uncertainty (noise) handling added (re-implemented), according 2901% to reference "A Method for Handling Uncertainty..." from below. 2902% 08/06/28: bug fix: file xrecent was empty 2903% 08/06/01: diagonalonly clean up. >1 means some iterations. 2904% 08/05/05: output is written to file preventing an increasing data 2905% array and ease long runs. 2906% 08/03/27: DiagonalOnly<0 learns for -DiagonalOnly iterations only the 2907% diagonal with a larger learning rate. 2908% 08/03 (2.60): option DiagonalOnly>=1 invokes a time- and space-linear 2909% variant with only diagonal elements of the covariance matrix 2910% updating. This can be useful for large dimensions, say > 100. 2911% 08/02: diag(weights) * ... replaced with repmat(weights,1,N) .* ... 2912% in C update, implies O(mu*N^2) instead of O(mu^2*N + mu*N^2). 2913% 07/09: tolhistfun as termination criterion added, "<" changed to 2914% "<=" also for TolFun to allow for stopping on zero difference. 2915% Name tolfunhist clashes with option tolfun. 2916% 07/07: hsig threshold made slighly smaller for large dimension, 2917% useful for lambda < lambda_default. 2918% 07/06: boundary handling: scaling in the boundary handling 2919% is omitted now, see bnd.flgscale. This seems not to 2920% have a big impact. Using the scaling is worse on rotated 2921% functions, but better on separable ones. 2922% 07/05: boundary handling: weight i is not incremented anymore 2923% if xmean(i) moves towards the feasible space. Increment 2924% factor changed to 1.2 instead of 1.1. 2925% 07/05: boundary handling code simplified not changing the algorithm 2926% 07/04: bug removed for saving in octave 2927% 06/11/10: more testing of outcome of eig, fixed max(D) to max(diag(D)) 2928% 06/10/21: conclusive final bestever assignment in the end 2929% 06/10/21: restart and incpopsize option implemented for restarts 2930% with increasing population size, version 2.50. 2931% 06/09/16: output argument bestever inserted again for convenience and 2932% backward compatibility 2933% 06/08: output argument out and struct out reorganized. 2934% 06/01: Possible parallel evaluation included as option EvalParallel 2935% 05/11: Compatibility to octave implemented, package octave-forge 2936% is needed. 2937% 05/09: Raise of figure and waiting for first plots improved 2938% 05/01: Function coordinatesystem cleaned up. 2939% 05/01: Function prctile, which requires the statistics toolbox, 2940% replaced by myprctile. 2941% 05/01: Option warnonequalfunctionvalues included. 2942% 04/12: Decrease of sigma removed. Problems on fsectorsphere can 2943% be addressed better by adding search space boundaries. 2944% 04/12: Boundary handling simpyfied. 2945% 04/12: Bug when stopping criteria tolx or tolupx are vectors. 2946% 04/11: Three input parameters are obligatory now. 2947% 04/11: Bug in boundary handling removed: Boundary weights can decrease now. 2948% 04/11: Normalization for boundary weights scale changed. 2949% 04/11: VerboseModulo option bug removed. Documentation improved. 2950% 04/11: Condition for increasing boundary weights changed. 2951% 04/10: Decrease of sigma when fitness is getting consistenly 2952% worse. Addresses the problems appearing on fsectorsphere for 2953% large population size. 2954% 04/10: VerboseModulo option included. 2955% 04/10: Bug for condition for increasing boundary weights removed. 2956% 04/07: tolx depends on initial sigma to achieve scale invariance 2957% for this stopping criterion. 2958% 04/06: Objective function value NaN is not counted as function 2959% evaluation and invokes resampling of the search point. 2960% 04/06: Error handling for eigenvalue beeing zero (never happens 2961% with default parameter setting) 2962% 04/05: damps further tuned for large mueff 2963% o Details for stall of pc-adaptation added (variable hsig 2964% introduced). 2965% 04/05: Bug in boundary handling removed: A large initial SIGMA was 2966% corrected not until *after* the first iteration, which could 2967% lead to a complete failure. 2968% 04/05: Call of function range (works with stats toolbox only) 2969% changed to myrange. 2970% 04/04: Parameter cs depends on mueff now and damps \propto sqrt(mueff) 2971% instead of \propto mueff. 2972% o Initial stall to adapt C (flginiphase) is removed and 2973% adaptation of pc is stalled for large norm(ps) instead. 2974% o Returned default options include documentation. 2975% o Resume part reorganized. 2976% 04/03: Stopflag becomes cell-array. 2977 2978% --------------------------------------------------------------- 2979% CMA-ES: Evolution Strategy with Covariance Matrix Adaptation for 2980% nonlinear function minimization. To be used under the terms of the 2981% GNU General Public License (http://www.gnu.org/copyleft/gpl.html). 2982% Author (copyright): Nikolaus Hansen, 2001-2008. 2983% e-mail: nikolaus.hansen AT inria.fr 2984% URL:http://www.bionik.tu-berlin.de/user/niko 2985% References: See below. 2986% --------------------------------------------------------------- 2987% 2988% GENERAL PURPOSE: The CMA-ES (Evolution Strategy with Covariance 2989% Matrix Adaptation) is a robust search method which should be 2990% applied, if derivative based methods, e.g. quasi-Newton BFGS or 2991% conjucate gradient, (supposably) fail due to a rugged search 2992% landscape (e.g. noise, local optima, outlier, etc.). On smooth 2993% landscapes CMA-ES is roughly ten times slower than BFGS. For up to 2994% N=10 variables even the simplex direct search method (Nelder & Mead) 2995% is often faster, but far less robust than CMA-ES. To see the 2996% advantage of the CMA, it will usually take at least 30*N and up to 2997% 300*N function evaluations, where N is the search problem dimension. 2998% On considerably hard problems the complete search (a single run) is 2999% expected to take at least 30*N^2 and up to 300*N^2 function 3000% evaluations. 3001% 3002% SOME MORE COMMENTS: 3003% The adaptation of the covariance matrix (e.g. by the CMA) is 3004% equivalent to a general linear transformation of the problem 3005% coding. Nevertheless every problem specific knowlegde about the best 3006% linear transformation should be exploited before starting the 3007% search. That is, an appropriate a priori transformation should be 3008% applied to the problem. This also makes the identity matrix as 3009% initial covariance matrix the best choice. 3010% 3011% The strategy parameter lambda (population size, opts.PopSize) is the 3012% preferred strategy parameter to play with. If results with the 3013% default strategy are not satisfactory, increase the population 3014% size. (Remark that the crucial parameter mu (opts.ParentNumber) is 3015% increased proportionally to lambda). This will improve the 3016% strategies capability of handling noise and local minima. We 3017% recomment successively increasing lambda by a factor of about three, 3018% starting with initial values between 5 and 20. Casually, population 3019% sizes even beyond 1000+100*N can be sensible. 3020% 3021% 3022% --------------------------------------------------------------- 3023%%% REFERENCES 3024% 3025% The equation numbers refer to 3026% Hansen, N. and S. Kern (2004). Evaluating the CMA Evolution 3027% Strategy on Multimodal Test Functions. Eighth International 3028% Conference on Parallel Problem Solving from Nature PPSN VIII, 3029% Proceedings, pp. 282-291, Berlin: Springer. 3030% (http://www.bionik.tu-berlin.de/user/niko/ppsn2004hansenkern.pdf) 3031% 3032% Further references: 3033% Hansen, N. and A. Ostermeier (2001). Completely Derandomized 3034% Self-Adaptation in Evolution Strategies. Evolutionary Computation, 3035% 9(2), pp. 159-195. 3036% (http://www.bionik.tu-berlin.de/user/niko/cmaartic.pdf). 3037% 3038% Hansen, N., S.D. Mueller and P. Koumoutsakos (2003). Reducing the 3039% Time Complexity of the Derandomized Evolution Strategy with 3040% Covariance Matrix Adaptation (CMA-ES). Evolutionary Computation, 3041% 11(1). (http://mitpress.mit.edu/journals/pdf/evco_11_1_1_0.pdf). 3042% 3043% Ros, R. and N. Hansen (2008). A Simple Modification in CMA-ES 3044% Achieving Linear Time and Space Complexity. To appear in Tenth 3045% International Conference on Parallel Problem Solving from Nature 3046% PPSN X, Proceedings, Berlin: Springer. 3047% 3048% Hansen, N., A.S.P. Niederberger, L. Guzzella and P. Koumoutsakos 3049% (2009?). A Method for Handling Uncertainty in Evolutionary 3050% Optimization with an Application to Feedback Control of 3051% Combustion. To appear in IEEE Transactions on Evolutionary 3052% Computation. 3053