barecmaes2

1 #!/usr/bin/env python 2 """Module barecmaes2 implements the CMA-ES without using `numpy`. 3 4 The Covariance Matrix Adaptation Evolution Strategy, CMA-ES, serves for 5 nonlinear function minimization. 6 7 The **main functionality** is implemented in 8 9 - class `CMAES` and 10 - function `fmin()` that is a single-line-usage wrapper around `CMAES`. 11 12 This code has two **purposes**: 13 14 1. it might be used for READING and UNDERSTANDING the basic flow and 15 the details of the CMA-ES *algorithm*. The source code is meant to 16 be read. For short, study the class `CMAES`, in particular its doc 17 string and the code of the method `CMAES.tell()`, where all the real 18 work is done in about 20 lines (see "def tell" in the source). 19 Otherwise, reading from the top is a feasible option. 20 21 2. it might be used when the python module `numpy` is not available. 22 When `numpy` is available, rather use `cma.py` (see http://www.lri 23 .fr/~hansen/cmaes_inmatlab.html#python) to run serious simulations: 24 the latter code has many more lines, but executes much faster 25 (roughly ten times), offers a richer user interface, far better 26 termination options, supposedly quite useful output, and automatic 27 restarts. 28 29 Dependencies: `math.exp`, `math.log` and `random.normalvariate` (modules 30 `matplotlib.pylab` and `sys` are optional). 31 32 Testing: call ``python barecmaes2.py`` at the OS shell. Tested with 33 Python 2.5 (only with removed import of print_function), 2.6, 2.7, 3.2. 34 Small deviations between the versions are expected. 35 36 URL: http://www.lri.fr/~hansen/barecmaes2.py 37 38 Last change: May, 2014, version 2.00 39 40 :Author: Nikolaus Hansen, 2010-2011, hansen[at-symbol]lri.fr 41 :License: This code is released into the public domain (that is, you may 42 use and modify it however you like). 43 44 """ 45 # import future syntax as for Python version >= 3.0 46 from __future__ import division # such that 0 != 1/2 == 0.5 47 from __future__ import print_function # available since 2.6, not needed 48 49 from math import log, exp 50 from random import normalvariate as random_normalvariate 51 # see randn keyword argument to CMAES.__init__ 52 53 # Optional imports, can be out-commented, if not available 54 import sys 55 try: 56 import matplotlib.pylab as pylab 57 pylab.ion() # prevents that execution stops after plotting 58 # for plotting, scitools.easyfiz might be an alternative 59 except ImportError: 60 pylab = None 61 print(' pylab could not be imported ') 62 63 __version__ = '1.10' 64 __author__ = 'Nikolaus Hansen' 65 __docformat__ = 'reStructuredText'

66 67 68 -def fmin(objectivefct, xstart, sigma, args=(), 69 max_eval='1e3*N**2', ftarget=None, 70 verb_disp=100, verb_log=1, verb_plot=100):

71 """non-linear non-convex minimization procedure. 72 The functional interface to CMA-ES. 73 74 Parameters 75 ========== 76 `objectivefct` 77 a function that takes as input a list of floats (like 78 [3.0, 2.2, 1.1]) and returns a single float (a scalar). 79 The objective is to find ``x`` with ``objectivefct(x)`` 80 to be as small as possible. 81 `xstart` 82 list of numbers (like `[3,2,1.2]`), initial solution vector 83 `sigma` 84 float, initial step-size, standard deviation in any coordinate 85 `args` 86 arguments to `objectivefct` 87 `ftarget` 88 float, target function value 89 `max_eval` 90 int or str, maximal number of function evaluations, a string 91 is evaluated with N being the search space dimension 92 `verb_disp` 93 int, display on console every verb_disp iteration, 0 for never 94 `verb_log` 95 int, data logging every verb_log iteration, 0 for never 96 `verb_plot` 97 int, plot logged data every verb_plot iteration, 0 for never 98 99 Returns 100 ======= 101 ``return es.result() + (es.stop(), es, logger)``, that is, 102 ``(xbest, fbest, evalsbest, evals, iterations, xmean,`` 103 `` termination_condition, CMAES_object_instance, data_logger)`` 104 105 Example 106 ======= 107 The following example minimizes the function `Fcts.elli`:: 108 109 >> import barecmaes2 as cma 110 >> res = cma.fmin(cma.Fcts.elli, 10 * [0.5], 0.3, verb_disp=100) 111 evals: ax-ratio max(std) f-value 112 10: 1.0 2.8e-01 198003.585517 113 1000: 8.4 5.5e-02 95.9162313173 114 2000: 40.5 3.6e-02 5.07618122556 115 3000: 149.1 8.5e-03 0.271537247667 116 4000: 302.2 4.2e-03 0.0623570374451 117 5000: 681.1 5.9e-03 0.000485971681802 118 6000: 1146.4 9.5e-06 5.26919100476e-10 119 6510: 1009.1 2.3e-07 3.34128914738e-13 120 termination by {'tolfun': 1e-12} 121 best f-value = 3.34128914738e-13 122 123 >> print(res[0]) 124 [2.1187532328944602e-07, 6.893386424102321e-08, -2.008255256456535e-09, 4.472078873398156e-09, -9.421306741003398e-09, 7.331265238205156e-09, 2.4804701814730273e-10, -6.030651566971234e-10, -6.063921614755129e-10, -1.066906137937511e-10] 125 126 >> print(res[1]) 127 3.34128914738e-13 128 129 >> res[-1].plot() # needs pylab/matplotlib to be installed 130 131 Details 132 ======= 133 By default `fmin()` tries to plot some output. This works only with 134 Python module `matplotlib` being installed. The first and last lines :: 135 136 cma.fmin(cma.Fcts.elli, 10 * [0.5], 0.3, verb_plot=0) 137 logger = cma.CMADataLogger() 138 cma.CMAES(10 * [0.5], 0.3).optimize(cma.Fcts.elli, logger=logger) 139 140 do the same. `fmin()` adds the possibility to plot *during* the 141 execution. 142 143 :See: `CMAES`, `OOOptimizer` 144 145 """ 146 es = CMAES(xstart, sigma, max_eval, ftarget) # new optimizer instance 147 logger = CMAESDataLogger(verb_log).add(es, True) # add data row 148 while not es.stop(): 149 es.optimize(objectivefct, args=args, 150 iterations=verb_plot if verb_plot else None, 151 verb_disp=verb_disp, 152 logger=logger) 153 if verb_plot: # and es.counteval/es.lam % verb_plot == 0: 154 logger.plot() 155 156 return es.result() + (es.stop(), es, logger)

157

158 159 -class OOOptimizer(object):

160 """"abstract" base class for an OO optimizer interface. 161 162 """

163 - def __init__(self, xstart, **more_kwargs):

164 """abstract method, ``xstart`` is a mandatory argument """ 165 raise NotImplementedError('method to be implemented in subclass')

166

167 - def ask(self):

168 """abstract method, AKA get, deliver new candidate solution(s), 169 a list of "vectors" """ 170 raise NotImplementedError('method to be implemented in subclass')

171

172 - def tell(self, solutions, function_values):

173 """abstract method, AKA update, prepare for next iteration""" 174 raise NotImplementedError('method to be implemented in subclass') 175

176 - def stop(self):

177 """abstract method, return satisfied termination conditions in 178 a dictionary like ``{'termination reason':value, ...}``, 179 for example ``{'tolfun':1e-12}``, or the empty dictionary ``{}``. 180 The implementation of `stop()` should prevent an infinite loop. 181 182 """ 183 raise NotImplementedError('method to be implemented in subclass') 184

185 - def result(self):

186 """abstract method, return ``(x, f(x), ...)``, that is the 187 minimizer, its function value, ...""" 188 raise NotImplementedError('method to be implemented in subclass') 189

190 - def disp(self, verbosity_modulo=1):

191 """display of some iteration info""" 192 print("method disp of " + str(type(self)) + " is not implemented")

193

194 - def optimize(self, objectivefct, iterations=None, min_iterations=1, 195 args=(), verb_disp=20, logger=None):

196 """iterate at least ``min_iterations`` and at most ``iterations`` 197 over ``OOOptimizer self`` using objective function 198 ``objectivefct``. Prints current information every ``verb_disp``, 199 uses ``OptimDataLogger logger``, and returns ``self``. 200 201 Example 202 ======= 203 :: 204 205 import barecmaes2 as cma 206 es = cma.CMAES(7 * [0.1], 0.5).optimize(cma.Fcts.rosenbrock) 207 print(es.result()[0]) 208 209 """ 210 if iterations is not None and iterations < min_iterations: 211 raise ValueError("iterations was smaller than min_iterations") 212 if iterations == 0: 213 raise ValueError("parameter iterations = 0") 214 215 iteration = 0 216 while not self.stop() or iteration < min_iterations: 217 if iterations and iteration >= iterations: 218 return self 219 iteration += 1 220 221 X = self.ask() # deliver candidate solutions 222 fitvals = [objectivefct(x, *args) for x in X] 223 self.tell(X, fitvals) # all the work is done here 224 225 self.disp(verb_disp) 226 logger.add(self) if logger else None 227 228 logger.add(self, force=True) if logger else None 229 if verb_disp: 230 self.disp(1) 231 print('termination by', self.stop()) 232 print('best f-value =', self.result()[1]) 233 print('solution =', self.result()[0]) 234 235 return self

236

237 -class CMAES(OOOptimizer):

238 """class for non-linear non-convex minimization. The class implements 239 the interface define in `OOOptimizer`, namely the methods 240 `__init__()`, `ask()`, `tell()`, `stop()`, `result()`, and `disp()`. 241 242 Examples 243 -------- 244 All examples minimize the function `elli`, the output is not shown. 245 (A preferred environment to execute all examples is 246 ``ipython -pylab``.) 247 248 First we need to :: 249 250 import barecmaes2 as cma 251 252 The shortest example uses the inherited method 253 `OOOptimizer.optimize()`:: 254 255 res = cma.CMAES(8 * [0.1], 0.5).optimize(cma.Fcts.elli) 256 257 See method `__init__()` for the input parameters to `CMAES`. We might 258 have a look at the result:: 259 260 print(res[0]) # best solution and 261 print(res[1]) # its function value 262 263 `res` is the return value from method `CMAES.results()` appended with 264 `None` (no logger). In order to display more exciting output we rather 265 do :: 266 267 logger = cma.CMAESDataLogger() 268 res = cma.CMAES(9 * [0.5], 0.3).optimize(cma.Fcts.elli, logger) 269 logger.plot() # if matplotlib is available 270 271 Virtually the same example can be written with an explicit loop 272 instead of using `optimize()`. This gives the necessary insight 273 into the `CMAES` class interface and gives entire control over the 274 iteration loop :: 275 276 optim = cma.CMAES(9 * [0.5], 0.3) # calls CMAES.__init__() 277 logger = cma.CMAESDataLogger().register(optim) # logger instance 278 279 # this loop resembles optimize() 280 while not optim.stop(): # iterate 281 X = optim.ask() # get candidate solutions 282 # do whatever needs to be done, however rather don't 283 # change X unless like for example X[2] = optim.ask()[0] 284 f = [cma.Fcts.elli(x) for x in X] # evaluate solutions 285 optim.tell(X, f) # do all the real work 286 optim.disp(20) # display info every 20th iteration 287 logger.add() # log another "data line" 288 289 # final output 290 print('termination by', optim.stop()) 291 print('best f-value =', optim.result()[1]) 292 print('best solution =', optim.result()[0]) 293 294 print('potentially better solution xmean =', optim.result()[5]) 295 print('let\'s check f(xmean) =', cma.Fcts.elli(optim.result()[5])) 296 logger.plot() # if matplotlib is available 297 raw_input('press enter to continue') # prevents closing figures 298 299 A slightly longer example, which also plots within the loop, is 300 the implementation of function `fmin(...)`. 301 302 Details 303 ------- 304 Most of the work is done in the method `tell(...)`. The method 305 `result()` returns more useful output. 306 307 :See: `fmin()`, `OOOptimizer.optimize()` 308 309 """

310 - def __init__(self, xstart, sigma, # mandatory 311 max_eval='1e3*N**2', 312 ftarget=None, 313 popsize="4 + int(3 * log(N))", 314 randn=random_normalvariate):

315 """Initialize` CMAES` object instance, the first two arguments are 316 mandatory. 317 318 Parameters 319 ---------- 320 `xstart` 321 ``list`` of numbers (like ``[3, 2, 1.2]``), initial 322 solution vector 323 `sigma` 324 ``float``, initial step-size (standard deviation in each 325 coordinate) 326 `max_eval` 327 ``int`` or ``str``, maximal number of function 328 evaluations, a string is evaluated with ``N`` being the 329 search space dimension 330 `ftarget` 331 `float`, target function value 332 `randn` 333 normal random number generator, by default 334 ``random.normalvariate`` 335 336 """ 337 # process input parameters 338 N = len(xstart) # number of objective variables/problem dimension 339 self.xmean = xstart[:] # initial point, distribution mean, a copy 340 self.sigma = sigma 341 self.ftarget = ftarget # stop if fitness < ftarget 342 self.max_eval = eval(str(max_eval)) # eval a string 343 self.randn = randn 344 345 # Strategy parameter setting: Selection 346 self.lam = eval(str(popsize)) # population size, offspring number 347 self.mu = int(self.lam / 2) # number of parents/points for recombination 348 self.weights = [log(self.mu+0.5) - log(i+1) for i in range(self.mu)] # recombination weights 349 self.weights = [w / sum(self.weights) for w in self.weights] # normalize recombination weights array 350 self.mueff = sum(self.weights)**2 / sum(w**2 for w in self.weights) # variance-effectiveness of sum w_i x_i 351 352 # Strategy parameter setting: Adaptation 353 self.cc = (4 + self.mueff/N) / (N+4 + 2 * self.mueff/N) # time constant for cumulation for C 354 self.cs = (self.mueff + 2) / (N + self.mueff + 5) # t-const for cumulation for sigma control 355 self.c1 = 2 / ((N + 1.3)**2 + self.mueff) # learning rate for rank-one update of C 356 self.cmu = min([1 - self.c1, 2 * (self.mueff - 2 + 1/self.mueff) / ((N + 2)**2 + self.mueff)]) # and for rank-mu update 357 self.damps = 2 * self.mueff/self.lam + 0.3 + self.cs # damping for sigma, usually close to 1 358 359 # Initialize dynamic (internal) state variables 360 self.pc, self.ps = N * [0], N * [0] # evolution paths for C,sigma 361 self.B = eye(N) # B defines the coordinate system 362 self.D = N * [1] # diagonal D defines the scaling 363 self.C = eye(N) # covariance matrix 364 self.invsqrtC = eye(N) # C^-1/2 365 self.eigeneval = 0 # tracking the update of B and D 366 self.counteval = 0 367 self.fitvals = [] # for bookkeeping output and termination 368 369 self.best = BestSolution()

370

371 - def stop(self):

372 """return satisfied termination conditions in a dictionary like 373 {'termination reason':value, ...}, for example {'tolfun':1e-12}, 374 or the empty dict {}""" 375 res = {} 376 if self.counteval > 0: 377 if self.counteval >= self.max_eval: 378 res['evals'] = self.max_eval 379 if self.ftarget is not None and len(self.fitvals) > 0 \ 380 and self.fitvals[0] <= self.ftarget: 381 res['ftarget'] = self.ftarget 382 if max(self.D) > 1e7 * min(self.D): 383 res['condition'] = 1e7 384 if len(self.fitvals) > 1 \ 385 and self.fitvals[-1] - self.fitvals[0] < 1e-12: 386 res['tolfun'] = 1e-12 387 if self.sigma * max(self.D) < 1e-11: 388 # remark: max(D) >= max(diag(C))**0.5 389 res['tolx'] = 1e-11 390 return res

391

392 - def ask(self):

393 """return a list of lambda candidate solutions according to 394 m + sig * Normal(0,C) = m + sig * B * D * Normal(0,I)""" 395 396 # Eigendecomposition: first update B, D and invsqrtC from C 397 # postpone in case to achieve O(N**2) 398 if self.counteval - self.eigeneval > \ 399 self.lam/(self.c1+self.cmu)/len(self.C)/10: 400 self.eigeneval = self.counteval 401 self.D, self.B = eig(self.C) # eigen decomposition, B==normalized eigenvectors, O(N**3) 402 self.D = [d**0.5 for d in self.D] # D contains standard deviations now 403 rg = range(len(self.C)) 404 for i in rg: # compute invsqrtC = C**(-1/2) = B D**(-1/2) B' 405 for j in rg: 406 self.invsqrtC[i][j] = sum(self.B[i][k] * self.B[j][k] 407 / self.D[k] for k in rg) 408 409 # lam vectors x_k = m + sigma * B * D * randn_k(N) 410 res = [] 411 for k in range(self.lam): # repeat lam times 412 z = [d * self.randn(0, 1) for d in self.D] 413 res.append(plus(self.xmean, dot1(self.sigma, dot(self.B, z)))) 414 return res

415

416 - def tell(self, arx, fitvals):

417 """update the evolution paths and the distribution parameters m, 418 sigma, and C within CMA-ES. 419 420 Parameters 421 ---------- 422 `arx` 423 a list of solutions, presumably from `ask()` 424 `fitvals` 425 the corresponding objective function values 426 427 """ 428 429 # bookkeeping, preparation 430 self.counteval += len(fitvals) # slightly artificial to do here 431 N = len(self.xmean) # convenience short cuts 432 iN = range(N) 433 xold = self.xmean 434 435 # Sort by fitness and compute weighted mean into xmean 436 self.fitvals, arindex = sorted(fitvals), argsort(fitvals) # min 437 arx = [arx[arindex[k]] for k in range(self.mu)] # sorted arx 438 del fitvals, arindex # to prevent misuse 439 # self.fitvals is kept for termination and display only 440 self.best.update([arx[0]], [self.fitvals[0]], self.counteval) 441 442 # xmean = [x_1=best, x_2, ..., x_mu] * weights 443 self.xmean = dot(arx[0:self.mu], self.weights, t=True) 444 # recombination, new mean value 445 # == [sum(self.weights[k] * arx[k][i] for k in range(self.mu)) 446 # for i in iN] 447 448 # Cumulation: update evolution paths 449 y = minus(self.xmean, xold) 450 z = dot(self.invsqrtC, y) # == C**(-1/2) * (xnew - xold) 451 452 c = (self.cs * (2-self.cs) * self.mueff)**0.5 / self.sigma 453 for i in iN: 454 self.ps[i] -= self.cs * self.ps[i] # exponential decay on ps 455 self.ps[i] += c * z[i] 456 hsig = (sum(x**2 for x in self.ps) 457 / (1-(1-self.cs)**(2*self.counteval/self.lam)) / N 458 < 2 + 4./(N+1)) 459 c = (self.cc * (2-self.cc) * self.mueff)**0.5 / self.sigma 460 for i in iN: 461 self.pc[i] -= self.cc * self.pc[i] # exponential decay on pc 462 self.pc[i] += c * hsig * y[i] 463 464 # Adapt covariance matrix C 465 c1a = self.c1 - (1-hsig**2) * self.c1 * self.cc * (2-self.cc) 466 # for a minor adjustment to the variance loss by hsig 467 for i in iN: 468 for j in iN: 469 Cmuij = sum(self.weights[k] * (arx[k][i]-xold[i]) 470 * (arx[k][j]-xold[j]) for k in range(self.mu) 471 ) / self.sigma**2 # rank-mu update 472 self.C[i][j] += (-c1a - self.cmu) * self.C[i][j] \ 473 + self.c1 * self.pc[i]*self.pc[j] + self.cmu * Cmuij 474 475 # Adapt step-size sigma with factor <= exp(0.6) \approx 1.82 476 self.sigma *= exp(min(0.6, (self.cs / self.damps) * 477 (sum(x**2 for x in self.ps)/N - 1)/2))

478 # chiN = (1-1./(4.*N)+1./(21.*N**2)) 479 # self.sigma *= exp(min(0.6, (self.cs/self.damps) * 480 # ((sum(x**2 for x in self.ps) / N)**0.5 / chiN - 1))) 481

482 - def result(self):

483 """return (xbest, f(xbest), evaluations_xbest, evaluations, 484 iterations, xmean) 485 486 """ 487 return self.best.get() + (self.counteval, 488 int(self.counteval / self.lam), 489 self.xmean)

490

491 - def disp(self, verb_modulo=1):

492 """display some iteration info""" 493 iteration = self.counteval / self.lam 494 495 if iteration == 1 or iteration % (10*verb_modulo) == 0: 496 print('evals: ax-ratio max(std) f-value') 497 if iteration <= 2 or iteration % verb_modulo == 0: 498 max_std = max([self.C[i][i] for i in range(len(self.C))])**0.5 499 print(repr(self.counteval).rjust(5) + ': ' + 500 ' %6.1f %8.1e ' % (max(self.D) / min(self.D), 501 self.sigma * max_std) + 502 str(self.fitvals[0])) 503 sys.stdout.flush() 504 505 return None

506

507 508 # ----------------------------------------------- 509 -class BestSolution(object):

510 """container to keep track of the best solution seen"""

511 - def __init__(self, x=None, f=None, evals=None):

512 """take `x`, `f`, and `evals` to initialize the best solution. 513 The better solutions have smaller `f`-values. """ 514 self.x, self.f, self.evals = x, f, evals

515

516 - def update(self, arx, arf, evals=None):

517 """initialize the best solution with `x`, `f`, and `evals`. 518 Better solutions have smaller `f`-values.""" 519 if self.f is None or min(arf) < self.f: 520 i = arf.index(min(arf)) 521 self.x, self.f = arx[i], arf[i] 522 self.evals = None if not evals else evals - len(arf) + i + 1 523 return self

524

525 - def get(self):

526 """return ``(x, f, evals)`` """ 527 return self.x, self.f, self.evals

528

529 530 # ----------------------------------------------- 531 -class OptimDataLogger(object):

532 """"abstract" base class for a data logger that can be used with an 533 `OOOptimizer`"""

534 - def __init__(self):

535 self.optim = None 536 self.dat = None

537

538 - def register(self, optim):

539 self.optim = optim 540 return self

541

542 - def add(self, optim=None, force=False):

543 """abstract method, add a "data point" from the state of optim 544 into the logger""" 545 raise NotImplementedError('to be implemented in a derived class')

546

547 - def disp(self):

548 """display some data trace (not implemented)""" 549 print('method OptimDataLogger.disp() is not implemented in class ' 550 + str(type(self)))

551

552 - def plot(self):

553 """plot data""" 554 print('method OptimDataLogger.plot() is not implemented in class ' 555 + str(type(self)))

556

557 - def data(self):

558 """return logged data in a dictionary""" 559 return self.dat

560

561 562 # ----------------------------------------------- 563 -class CMAESDataLogger(OptimDataLogger):

564 """data logger for class CMAES, that can store and plot data. 565 566 An even more useful logger would write the data to files, to 567 prevent memory overload. 568 569 """ 570 571 plotted = 0 572 "plot count for all instances" 573

574 - def __init__(self, verb_modulo=1):

575 """cma is the `OOOptimizer` class instance that is to be logged, 576 `verb_modulo` controls whether logging takes place for each call 577 to the method `add()` 578 579 """ 580 OptimDataLogger.__init__(self) # not necessary 581 self.modulo = verb_modulo 582 self.dat = {'eval': [], 'iter': [], 'stds': [], 'D': [], 583 'sig': [], 'fit': [], 'xm': []} 584 self.counter = 0 # number of calls of add

585

586 - def add(self, cma=None, force=False):

587 """append some logging data from CMAES class instance cma, 588 if ``number_of_times_called modulo verb_modulo`` equals zero""" 589 cma = self.optim if cma is None else cma 590 if type(cma) is not CMAES: 591 raise TypeError('logged object must be a CMAES instance') 592 dat = self.dat 593 self.counter += 1 594 if force and self.counter == 1: 595 self.counter = 0 596 if (self.modulo 597 and 598 (len(dat['eval']) == 0 or cma.counteval != dat['eval'][-1]) 599 and 600 (self.counter < 4 or force or 601 int(self.counter) % self.modulo == 0)): 602 dat['eval'].append(cma.counteval) 603 dat['iter'].append(cma.counteval/cma.lam) 604 dat['stds'].append([cma.C[i][i]**0.5 605 for i in range(len(cma.C))]) 606 dat['D'].append(sorted(cma.D)) 607 dat['sig'].append(cma.sigma) 608 dat['fit'].append(cma.fitvals[0] if hasattr(cma, 'fitvals') 609 and cma.fitvals 610 else None) 611 dat['xm'].append([x for x in cma.xmean]) 612 return self

613

614 - def plot(self, fig_number=322):

615 """plot the stored data in figure fig_number. 616 617 Dependencies: matlabplotlib/pylab. 618 619 Example 620 ======= 621 :: 622 623 >> import barecmaes2 as cma 624 >> es = cma.CMAES(3 * [0.1], 1) 625 >> logger = cma.CMAESDataLogger().register(es) 626 >> while not es.stop(): 627 >> X = es.ask() 628 >> es.tell(X, [bc.Fcts.elli(x) for x in X]) 629 >> logger.add() 630 >> logger.plot() 631 632 """ 633 if pylab is None: 634 return None 635 pylab.figure(fig_number) 636 637 from pylab import text, hold, plot, ylabel, grid, semilogy, \ 638 xlabel, draw, show, subplot 639 640 dat = self.dat # dictionary with entries as given in __init__ 641 if not dat: 642 return 643 try: # a hack to get the presumable population size lambda 644 strpopsize = ' (popsize~' + str(dat['eval'][-2] - 645 dat['eval'][-3]) + ')' 646 except IndexError: 647 strpopsize = '' 648 649 # plot fit, Delta fit, sigma 650 subplot(221) 651 hold(False) 652 if dat['fit'][0] is None: 653 dat['fit'][0] = dat['fit'][1] 654 # should be reverted later, but let's be lazy 655 assert dat['fit'].count(None) == 0 656 dmin = min(dat['fit']) 657 i = dat['fit'].index(dmin) 658 dat['fit'][i] = max(dat['fit']) 659 dmin2 = min(dat['fit']) 660 dat['fit'][i] = dmin 661 semilogy(dat['iter'], [d - dmin + 1e-19 if d >= dmin2 else 662 dmin2 - dmin for d in dat['fit']], 'c', 663 linewidth=1) 664 hold(True) 665 semilogy(dat['iter'], [abs(d) for d in dat['fit']], 'b') 666 semilogy(dat['iter'][i], abs(dmin), 'r*') 667 semilogy(dat['iter'], dat['sig'], 'g') 668 ylabel('f-value, Delta-f-value, sigma') 669 grid(True) 670 671 # plot xmean 672 subplot(222) 673 hold(False) 674 plot(dat['iter'], dat['xm']) 675 hold(True) 676 for i in range(len(dat['xm'][-1])): 677 text(dat['iter'][0], dat['xm'][0][i], str(i)) 678 text(dat['iter'][-1], dat['xm'][-1][i], str(i)) 679 ylabel('mean solution', ha='center') 680 grid(True) 681 682 # plot D 683 subplot(223) 684 hold(False) 685 semilogy(dat['iter'], dat['D'], 'm') 686 xlabel('iterations' + strpopsize) 687 ylabel('axes lengths') 688 grid(True) 689 690 # plot stds 691 subplot(224) 692 # if len(gcf().axes) > 1: 693 # sca(pylab.gcf().axes[1]) 694 # else: 695 # twinx() 696 hold(False) 697 semilogy(dat['iter'], dat['stds']) 698 for i in range(len(dat['stds'][-1])): 699 text(dat['iter'][-1], dat['stds'][-1][i], str(i)) 700 ylabel('coordinate stds disregarding sigma', ha='center') 701 grid(True) 702 xlabel('iterations' + strpopsize) 703 sys.stdout.flush() 704 draw() 705 show() 706 CMAESDataLogger.plotted += 1

707

708 709 #_____________________________________________________________________ 710 #_______________________ Helper Functions _____________________________ 711 # 712 -def eye(dimension):

713 """return identity matrix as list of "vectors" """ 714 m = [dimension * [0] for i in range(dimension)] 715 # m = N * [N * [0]] fails because it gives N times the same reference 716 for i in range(dimension): 717 m[i][i] = 1 718 return m

719

720 721 -def dot(A, b, t=False):

722 """ usual dot product of "matrix" A with "vector" b, 723 where A[i] is the i-th row of A. With t=True, A transposed is used. 724 725 :rtype : "vector" (list of float) 726 """ 727 n = len(b) 728 if not t: 729 m = len(A) # number of rows, like printed by pprint 730 v = m * [0] 731 for i in range(m): 732 v[i] = sum(b[j] * A[i][j] for j in range(n)) 733 else: 734 m = len(A[0]) # number of columns 735 v = m * [0] 736 for i in range(m): 737 v[i] = sum(b[j] * A[j][i] for j in range(n)) 738 return v

739

740 741 -def dot1(a, b):

742 """scalar a times vector b""" 743 return [a * c for c in b]

744

745 746 -def plus(a, b):

747 """add vectors, return a + b """ 748 return [a[i] + b[i] for i in range(len(a))]

749

750 751 -def minus(a, b):

752 """substract vectors, return a - b""" 753 return [a[i] - b[i] for i in range(len(a))]

754

755 756 -def argsort(a):

757 """return index list to get a in order, ie 758 a[argsort(a)[i]] == sorted(a)[i]""" 759 return [a.index(val) for val in sorted(a)]

760

761 762 #_____________________________________________________________________ 763 #_________________ Fitness (Objective) Functions _____________________ 764 765 -class Fcts(object): # instead of a submodule

766 """versatile collection of test functions""" 767 768 @staticmethod # syntax available since 2.4

769 - def elli(x):

770 """ellipsoid test objective function""" 771 n = len(x) 772 aratio = 1e3 773 return sum(x[i]**2 * aratio**(2.*i/(n-1)) for i in range(n))

774 775 @staticmethod

776 - def sphere(x):

777 """sphere, ``sum(x**2)``, test objective function""" 778 return sum(x[i]**2 for i in range(len(x)))

779 780 @staticmethod

781 - def rosenbrock(x):

782 """Rosenbrock test objective function""" 783 n = len(x) 784 if n < 2: 785 raise ValueError('dimension must be greater one') 786 return sum(100 * (x[i]**2 - x[i+1])**2 + (x[i] - 1)**2 for i 787 in range(n-1))

788

789 790 #____________________________________________________________ 791 #____________________________________________________________ 792 # 793 # C and B are arrays rather than matrices, because they are 794 # addressed via B[i][j], matrices can only be addressed via B[i,j] 795 796 # tred2(N, B, diagD, offdiag); 797 # tql2(N, diagD, offdiag, B); 798 799 # Symmetric Householder reduction to tridiagonal form, translated from 800 # JAMA package. 801 802 -def eig(C):

803 """eigendecomposition of a symmetric matrix, much slower than 804 numpy.linalg.eigh, return the eigenvalues and an orthonormal basis 805 of the corresponding eigenvectors, ``(EVals, Basis)``, where 806 807 ``Basis[i]`` 808 the i-th row of ``Basis`` and columns of ``Basis``, 809 ``[Basis[j][i] for j in range(len(Basis))]`` 810 the i-th eigenvector with eigenvalue ``EVals[i]`` 811 812 """ 813 # class eig(object): 814 # def __call__(self, C): 815 816 # Householder transformation of a symmetric matrix V into tridiagonal 817 # form. 818 # -> n : dimension 819 # -> V : symmetric nxn-matrix 820 # <- V : orthogonal transformation matrix: 821 # tridiag matrix == V * V_in * V^t 822 # <- d : diagonal 823 # <- e[0..n-1] : off diagonal (elements 1..n-1) 824 825 # Symmetric tridiagonal QL algorithm, iterative 826 # Computes the eigensystem from a tridiagonal matrix in roughtly 3N^3 827 # operations 828 # -> n : Dimension. 829 # -> d : Diagonale of tridiagonal matrix. 830 # -> e[1..n-1] : off-diagonal, output from Householder 831 # -> V : matrix output von Householder 832 # <- d : eigenvalues 833 # <- e : garbage? 834 # <- V : basis of eigenvectors, according to d 835 836 # tred2(N, B, diagD, offdiag); B=C on input 837 # tql2(N, diagD, offdiag, B); 838 839 # private void tred2 (int n, double V[][], double d[], double e[]) { 840 def tred2(n, V, d, e): 841 # This is derived from the Algol procedures tred2 by 842 # Bowdler, Martin, Reinsch, and Wilkinson, Handbook for 843 # Auto. Comp., Vol.ii-Linear Algebra, and the corresponding 844 # Fortran subroutine in EISPACK. 845 846 num_opt = False # factor 1.5 in 30-D 847 if num_opt: 848 import numpy as np 849 850 for j in range(n): 851 d[j] = V[n-1][j] # d is output argument 852 853 # Householder reduction to tridiagonal form. 854 855 for i in range(n-1, 0, -1): 856 # Scale to avoid under/overflow. 857 h = 0.0 858 if not num_opt: 859 scale = 0.0 860 for k in range(i): 861 scale = scale + abs(d[k]) 862 else: 863 scale = sum(abs(d[0:i])) 864 865 if scale == 0.0: 866 e[i] = d[i-1] 867 for j in range(i): 868 d[j] = V[i-1][j] 869 V[i][j] = 0.0 870 V[j][i] = 0.0 871 else: 872 873 # Generate Householder vector. 874 if not num_opt: 875 for k in range(i): 876 d[k] /= scale 877 h += d[k] * d[k] 878 else: 879 d[:i] /= scale 880 h = np.dot(d[:i], d[:i]) 881 882 f = d[i-1] 883 g = h**0.5 884 885 if f > 0: 886 g = -g 887 888 e[i] = scale * g 889 h -= f * g 890 d[i-1] = f - g 891 if not num_opt: 892 for j in range(i): 893 e[j] = 0.0 894 else: 895 e[:i] = 0.0 896 897 # Apply similarity transformation to remaining columns. 898 899 for j in range(i): 900 f = d[j] 901 V[j][i] = f 902 g = e[j] + V[j][j] * f 903 if not num_opt: 904 for k in range(j+1, i): 905 g += V[k][j] * d[k] 906 e[k] += V[k][j] * f 907 e[j] = g 908 else: 909 e[j+1:i] += V.T[j][j+1:i] * f 910 e[j] = g + np.dot(V.T[j][j+1:i], d[j+1:i]) 911 912 f = 0.0 913 if not num_opt: 914 for j in range(i): 915 e[j] /= h 916 f += e[j] * d[j] 917 else: 918 e[:i] /= h 919 f += np.dot(e[:i], d[:i]) 920 921 hh = f / (h + h) 922 if not num_opt: 923 for j in range(i): 924 e[j] -= hh * d[j] 925 else: 926 e[:i] -= hh * d[:i] 927 928 for j in range(i): 929 f = d[j] 930 g = e[j] 931 if not num_opt: 932 for k in range(j, i): 933 V[k][j] -= (f * e[k] + g * d[k]) 934 else: 935 V.T[j][j:i] -= (f * e[j:i] + g * d[j:i]) 936 937 d[j] = V[i-1][j] 938 V[i][j] = 0.0 939 940 d[i] = h 941 # end for i-- 942 943 # Accumulate transformations. 944 945 for i in range(n-1): 946 V[n-1][i] = V[i][i] 947 V[i][i] = 1.0 948 h = d[i+1] 949 if h != 0.0: 950 if not num_opt: 951 for k in range(i+1): 952 d[k] = V[k][i+1] / h 953 else: 954 d[:i+1] = V.T[i+1][:i+1] / h 955 956 for j in range(i+1): 957 if not num_opt: 958 g = 0.0 959 for k in range(i+1): 960 g += V[k][i+1] * V[k][j] 961 for k in range(i+1): 962 V[k][j] -= g * d[k] 963 else: 964 g = np.dot(V.T[i+1][0:i+1], V.T[j][0:i+1]) 965 V.T[j][:i+1] -= g * d[:i+1] 966 967 if not num_opt: 968 for k in range(i+1): 969 V[k][i+1] = 0.0 970 else: 971 V.T[i+1][:i+1] = 0.0 972 973 if not num_opt: 974 for j in range(n): 975 d[j] = V[n-1][j] 976 V[n-1][j] = 0.0 977 else: 978 d[:n] = V[n-1][:n] 979 V[n-1][:n] = 0.0 980 981 V[n-1][n-1] = 1.0 982 e[0] = 0.0

983 984 # Symmetric tridiagonal QL algorithm, taken from JAMA package. 985 # private void tql2 (int n, double d[], double e[], double V[][]) { 986 # needs roughly 3N^3 operations 987 def tql2(n, d, e, V): 988 # This is derived from the Algol procedures tql2, by 989 # Bowdler, Martin, Reinsch, and Wilkinson, Handbook for 990 # Auto. Comp., Vol.ii-Linear Algebra, and the corresponding 991 # Fortran subroutine in EISPACK. 992 993 num_opt = False # using vectors from numpy make it faster 994 995 if not num_opt: 996 for i in range(1, n): # (int i = 1; i < n; i++): 997 e[i-1] = e[i] 998 else: 999 e[0:n-1] = e[1:n] 1000 e[n-1] = 0.0 1001 1002 f = 0.0 1003 tst1 = 0.0 1004 eps = 2.0**-52.0 1005 for l in range(n): # (int l = 0; l < n; l++) { 1006 1007 # Find small subdiagonal element 1008 1009 tst1 = max(tst1, abs(d[l]) + abs(e[l])) 1010 m = l 1011 while m < n: 1012 if abs(e[m]) <= eps*tst1: 1013 break 1014 m += 1 1015 1016 # If m == l, d[l] is an eigenvalue, 1017 # otherwise, iterate. 1018 1019 if m > l: 1020 iiter = 0 1021 while 1: # do { 1022 iiter += 1 # (Could check iteration count here.) 1023 1024 # Compute implicit shift 1025 1026 g = d[l] 1027 p = (d[l+1] - g) / (2.0 * e[l]) 1028 r = (p**2 + 1)**0.5 # hypot(p, 1.0) 1029 if p < 0: 1030 r = -r 1031 1032 d[l] = e[l] / (p + r) 1033 d[l+1] = e[l] * (p + r) 1034 dl1 = d[l+1] 1035 h = g - d[l] 1036 if not num_opt: 1037 for i in range(l+2, n): 1038 d[i] -= h 1039 else: 1040 d[l+2:n] -= h 1041 1042 f = f + h 1043 1044 # Implicit QL transformation. 1045 1046 p = d[m] 1047 c = 1.0 1048 c2 = c 1049 c3 = c 1050 el1 = e[l+1] 1051 s = 0.0 1052 s2 = 0.0 1053 1054 # hh = V.T[0].copy() # only with num_opt 1055 for i in range(m-1, l-1, -1): 1056 # (int i = m-1; i >= l; i--) { 1057 c3 = c2 1058 c2 = c 1059 s2 = s 1060 g = c * e[i] 1061 h = c * p 1062 r = (p**2 + e[i]**2)**0.5 # hypot(p,e[i]) 1063 e[i+1] = s * r 1064 s = e[i] / r 1065 c = p / r 1066 p = c * d[i] - s * g 1067 d[i+1] = h + s * (c * g + s * d[i]) 1068 1069 # Accumulate transformation. 1070 1071 if not num_opt: # overall factor 3 in 30-D 1072 for k in range(n): # (int k = 0; k < n; k++){ 1073 h = V[k][i+1] 1074 V[k][i+1] = s * V[k][i] + c * h 1075 V[k][i] = c * V[k][i] - s * h 1076 else: # about 20% faster in 10-D 1077 hh = V.T[i+1].copy() 1078 # hh[:] = V.T[i+1][:] 1079 V.T[i+1] = s * V.T[i] + c * hh 1080 V.T[i] = c * V.T[i] - s * hh 1081 # V.T[i] *= c 1082 # V.T[i] -= s * hh 1083 1084 p = -s * s2 * c3 * el1 * e[l] / dl1 1085 e[l] = s * p 1086 d[l] = c * p 1087 1088 # Check for convergence. 1089 if abs(e[l]) <= eps*tst1: 1090 break 1091 # } while (Math.abs(e[l]) > eps*tst1); 1092 1093 d[l] += f 1094 e[l] = 0.0 1095 1096 # Sort eigenvalues and corresponding vectors. 1097 if 11 < 3: 1098 for i in range(n-1): # (int i = 0; i < n-1; i++) { 1099 k = i 1100 p = d[i] 1101 for j in range(i+1, n): # (int j = i+1; j < n; j++) { 1102 if d[j] < p: # NH find smallest k>i 1103 k = j 1104 p = d[j] 1105 1106 if k != i: 1107 d[k] = d[i] # swap k and i 1108 d[i] = p 1109 for j in range(n): # (int j = 0; j < n; j++) { 1110 p = V[j][i] 1111 V[j][i] = V[j][k] 1112 V[j][k] = p 1113 # tql2 1114 N = len(C[0]) 1115 V = [C[i][:] for i in range(N)] 1116 d = N * [0] 1117 e = N * [0] 1118 tred2(N, V, d, e) 1119 tql2(N, d, e, V) 1120 return d, V 1121

1122 1123 -def _test():

1124 """ 1125 >>> import barecmaes2 as cma 1126 >>> import random 1127 >>> random.seed(5) 1128 >>> x = cma.fmin(cma.Fcts.rosenbrock, 4 * [0.5], 0.5, verb_plot=0) 1129 evals: ax-ratio max(std) f-value 1130 8: 1.0 4.3e-01 8.17705253283 1131 16: 1.1 4.2e-01 66.2003009423 1132 800: 19.1 3.8e-02 0.0423501042891 1133 1600: 46.1 4.7e-06 1.14801950628e-11 1134 1736: 59.2 5.4e-07 1.26488753189e-13 1135 termination by {'tolfun': 1e-12} 1136 best f-value = 4.93281796387e-14 1137 solution = [0.9999999878867273, 0.9999999602211, 0.9999999323618144, 0.9999998579200512] 1138 1139 """ 1140 1141 import doctest 1142 print('launching doctest') 1143 doctest.testmod(report=True) # module test 1144 print("""done (ideally no line between launching and done was printed, 1145 however slight deviations are possible)""")

1146 1147 #_____________________________________________________________________ 1148 #_____________________________________________________________________ 1149 # 1150 if __name__ == "__main__": 1151 1152 _test() 1153 1154 # fmin(Fcts.rosenbrock, 10 * [0.5], 0.5) 1155

Source Code for Module barecmaes2