1
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
46 from __future__ import division
47 from __future__ import print_function
48
49 from math import log, exp
50 from random import normalvariate as random_normalvariate
51
52
53
54 import sys
55 try:
56 import matplotlib.pylab as pylab
57 pylab.ion()
58
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)
147 logger = CMAESDataLogger(verb_log).add(es, True)
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:
154 logger.plot()
155
156 return es.result() + (es.stop(), es, logger)
157
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
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
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
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()
222 fitvals = [objectivefct(x, *args) for x in X]
223 self.tell(X, fitvals)
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,
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
338 N = len(xstart)
339 self.xmean = xstart[:]
340 self.sigma = sigma
341 self.ftarget = ftarget
342 self.max_eval = eval(str(max_eval))
343 self.randn = randn
344
345
346 self.lam = eval(str(popsize))
347 self.mu = int(self.lam / 2)
348 self.weights = [log(self.mu+0.5) - log(i+1) for i in range(self.mu)]
349 self.weights = [w / sum(self.weights) for w in self.weights]
350 self.mueff = sum(self.weights)**2 / sum(w**2 for w in self.weights)
351
352
353 self.cc = (4 + self.mueff/N) / (N+4 + 2 * self.mueff/N)
354 self.cs = (self.mueff + 2) / (N + self.mueff + 5)
355 self.c1 = 2 / ((N + 1.3)**2 + self.mueff)
356 self.cmu = min([1 - self.c1, 2 * (self.mueff - 2 + 1/self.mueff) / ((N + 2)**2 + self.mueff)])
357 self.damps = 2 * self.mueff/self.lam + 0.3 + self.cs
358
359
360 self.pc, self.ps = N * [0], N * [0]
361 self.B = eye(N)
362 self.D = N * [1]
363 self.C = eye(N)
364 self.invsqrtC = eye(N)
365 self.eigeneval = 0
366 self.counteval = 0
367 self.fitvals = []
368
369 self.best = BestSolution()
370
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
389 res['tolx'] = 1e-11
390 return res
391
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
397
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)
402 self.D = [d**0.5 for d in self.D]
403 rg = range(len(self.C))
404 for i in rg:
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
410 res = []
411 for k in range(self.lam):
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
430 self.counteval += len(fitvals)
431 N = len(self.xmean)
432 iN = range(N)
433 xold = self.xmean
434
435
436 self.fitvals, arindex = sorted(fitvals), argsort(fitvals)
437 arx = [arx[arindex[k]] for k in range(self.mu)]
438 del fitvals, arindex
439
440 self.best.update([arx[0]], [self.fitvals[0]], self.counteval)
441
442
443 self.xmean = dot(arx[0:self.mu], self.weights, t=True)
444
445
446
447
448
449 y = minus(self.xmean, xold)
450 z = dot(self.invsqrtC, y)
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]
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]
462 self.pc[i] += c * hsig * y[i]
463
464
465 c1a = self.c1 - (1-hsig**2) * self.c1 * self.cc * (2-self.cc)
466
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
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
476 self.sigma *= exp(min(0.6, (self.cs / self.damps) *
477 (sum(x**2 for x in self.ps)/N - 1)/2))
478
479
480
481
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
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
526 """return ``(x, f, evals)`` """
527 return self.x, self.f, self.evals
528
532 """"abstract" base class for a data logger that can be used with an
533 `OOOptimizer`"""
535 self.optim = None
536 self.dat = None
537
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
548 """display some data trace (not implemented)"""
549 print('method OptimDataLogger.disp() is not implemented in class '
550 + str(type(self)))
551
553 """plot data"""
554 print('method OptimDataLogger.plot() is not implemented in class '
555 + str(type(self)))
556
558 """return logged data in a dictionary"""
559 return self.dat
560
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
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)
581 self.modulo = verb_modulo
582 self.dat = {'eval': [], 'iter': [], 'stds': [], 'D': [],
583 'sig': [], 'fit': [], 'xm': []}
584 self.counter = 0
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
641 if not dat:
642 return
643 try:
644 strpopsize = ' (popsize~' + str(dat['eval'][-2] -
645 dat['eval'][-3]) + ')'
646 except IndexError:
647 strpopsize = ''
648
649
650 subplot(221)
651 hold(False)
652 if dat['fit'][0] is None:
653 dat['fit'][0] = dat['fit'][1]
654
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
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
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
691 subplot(224)
692
693
694
695
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
711
712 -def eye(dimension):
713 """return identity matrix as list of "vectors" """
714 m = [dimension * [0] for i in range(dimension)]
715
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)
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])
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
742 """scalar a times vector b"""
743 return [a * c for c in b]
744
747 """add vectors, return a + b """
748 return [a[i] + b[i] for i in range(len(a))]
749
752 """substract vectors, return a - b"""
753 return [a[i] - b[i] for i in range(len(a))]
754
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
764
765 -class Fcts(object):
766 """versatile collection of test functions"""
767
768 @staticmethod
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
777 """sphere, ``sum(x**2)``, test objective function"""
778 return sum(x[i]**2 for i in range(len(x)))
779
780 @staticmethod
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
794
795
796
797
798
799
800
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
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840 def tred2(n, V, d, e):
841
842
843
844
845
846 num_opt = False
847 if num_opt:
848 import numpy as np
849
850 for j in range(n):
851 d[j] = V[n-1][j]
852
853
854
855 for i in range(n-1, 0, -1):
856
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
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
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
942
943
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
985
986
987 def tql2(n, d, e, V):
988
989
990
991
992
993 num_opt = False
994
995 if not num_opt:
996 for i in range(1, n):
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):
1006
1007
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
1017
1018
1019 if m > l:
1020 iiter = 0
1021 while 1:
1022 iiter += 1
1023
1024
1025
1026 g = d[l]
1027 p = (d[l+1] - g) / (2.0 * e[l])
1028 r = (p**2 + 1)**0.5
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
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
1055 for i in range(m-1, l-1, -1):
1056
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
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
1070
1071 if not num_opt:
1072 for k in range(n):
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:
1077 hh = V.T[i+1].copy()
1078
1079 V.T[i+1] = s * V.T[i] + c * hh
1080 V.T[i] = c * V.T[i] - s * hh
1081
1082
1083
1084 p = -s * s2 * c3 * el1 * e[l] / dl1
1085 e[l] = s * p
1086 d[l] = c * p
1087
1088
1089 if abs(e[l]) <= eps*tst1:
1090 break
1091
1092
1093 d[l] += f
1094 e[l] = 0.0
1095
1096
1097 if 11 < 3:
1098 for i in range(n-1):
1099 k = i
1100 p = d[i]
1101 for j in range(i+1, n):
1102 if d[j] < p:
1103 k = j
1104 p = d[j]
1105
1106 if k != i:
1107 d[k] = d[i]
1108 d[i] = p
1109 for j in range(n):
1110 p = V[j][i]
1111 V[j][i] = V[j][k]
1112 V[j][k] = p
1113
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
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)
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
1155