#!/usr/bin/env python
"""Module cma implements the CMA-ES, Covariance Matrix Adaptation Evolution
Strategy, a stochastic optimizer for robust non-linear non-convex
derivative-free function minimization for Python versions 2.6 and 2.7
(in Python 2.5 the collections import does not provide the class MutableMapping).
CMA-ES searches for a minimizer (a solution x in R**n) of an
objective function f (cost function), such that f(x) is
minimal. Regarding f, only function values for candidate solutions
need to be available, gradients are not necessary. Even less
restrictive, only a passably reliable ranking of the candidate
solutions in each iteration is necessary, the function values
itself do not matter. Some termination criteria however depend
on actual f-values.
Two interfaces are provided:
- function `fmin(func, x0, sigma0,...)`
runs a complete minimization
of the objective function func with CMA-ES.
- class `CMAEvolutionStrategy`
allows for minimization such that the
control of the iteration loop remains with the user.
Used packages:
- unavoidable: `numpy` (see `barecmaes2.py` if `numpy` is not
available),
- avoidable with small changes: `time`, `sys`
- optional: `matplotlib.pylab` (for `plot` etc., highly
recommended), `pprint` (pretty print), `pickle` (in class
`Sections`), `doctest`, `inspect`, `pygsl` (never by default)
Testing
-------
The code can be tested on a given system. Typing::
python cma.py --test --quiet
or in the Python shell ``ipython -pylab``::
run cma.py --test --quiet
runs ``doctest.testmod(cma)`` showing only exceptions (and not the
tests that fail due to small differences in the output) and should
run without complaints in about under two minutes. On some systems,
the pop up windows must be closed manually to continue and finish
the test.
Example
-------
::
import cma
help(cma) # "this" help message, use cma? in ipython
help(cma.fmin)
help(cma.CMAEvolutionStrategy)
help(cma.Options)
cma.Options('tol') # display 'tolerance' termination options
cma.Options('verb') # display verbosity options
res = cma.fmin(cma.Fcts.tablet, 15 * [1], 1)
res[0] # best evaluated solution
res[5] # mean solution, presumably better with noise
:See: `fmin()`, `Options`, `CMAEvolutionStrategy`
:Author: Nikolaus Hansen, 2008-2012
:License: GPL 2 and 3
"""
from __future__ import division # future is >= 3.0, this code has been used with 2.5 & 2.6 & 2.7
from __future__ import with_statement # only necessary for python 2.5 and not in heavy use
# from __future__ import collections.MutableMapping # does not exist in future
# from __future__ import print_function # for cross-checking, available from python 2.6
__version__ = "0.91.00 $Revision: 3103 $"
# $Date: 2012-01-31 15:14:46 +0100 (Tue, 31 Jan 2012) $
# bash: svn propset svn:keywords 'Date Revision' cma.py
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, version 2 or 3.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see .
#
# for testing:
# pyflakes cma.py # finds bugs by static analysis
# pychecker --limit 60 cma.py # also executes, gives 60 warnings (all checked)
# python cma.py -t -quiet # executes implemented tests based on doctest
# to create a html documentation file:
# pydoc -w cma # edit the header (remove local pointers)
# epydoc cma.py # comes close to javadoc but does not find the
# # links of function references etc
# doxygen needs @package cma as first line in the module docstring
# some things like class attributes are not interpreted correctly
# sphinx: doc style of doc.python.org, could not make it work
# TODO: add a default logger in CMAEvolutionStrategy, see fmin() and optimize() first
# tell() should probably not add data. Should optimize() handle even an after_iteration_handler?
# TODO: CMAEvolutionStrategy(ones(10), 1).optimize(cma.fcts.elli) # should work like fmin
# one problem: the data logger is not default and seemingly cannot be attached in one line
# TODO: check combination of boundary handling and transformation: penalty must be computed
# on gp.pheno(x_geno, bounds=None), but without bounds, check/remove usage of .geno everywhere
# TODO: check whether all new solutions are put into self.sent_solutions
# TODO: separate initialize==reset_state from __init__
# TODO: introduce Zpos == diffC which makes the code more consistent and the active update "exact"
# TODO: split tell into a variable transformation part and the "pure" functionality
# usecase: es.tell_geno(X, [func(es.pheno(x)) for x in X])
# genotypic repair is not part of tell_geno
# TODO: read settable "options" from a (properties) file, see myproperties.py
#
# typical parameters in scipy.optimize: disp, xtol, ftol, maxiter, maxfun, callback=None
# maxfev, diag (A sequency of N positive entries that serve as
# scale factors for the variables.)
# full_output -- non-zero to return all optional outputs.
# If xtol < 0.0, xtol is set to sqrt(machine_precision)
# 'infot -- a dictionary of optional outputs with the keys:
# 'nfev': the number of function calls...
#
# see eg fmin_powell
# typical returns
# x, f, dictionary d
# (xopt, {fopt, gopt, Hopt, func_calls, grad_calls, warnflag}, )
#
# TODO: keep best ten solutions
# TODO: implement constraints handling
# TODO: option full_output -- non-zero to return all optional outputs.
# TODO: extend function unitdoctest, or use unittest?
# TODO: implement equal-fitness termination, covered by stagnation?
# TODO: apply style guide: no capitalizations!?
# TODO: check and test dispdata()
# TODO: eigh(): thorough testing would not hurt
#
# TODO (later): implement readSignals from a file like properties file (to be called after tell())
import sys, time, collections # not really essential
import numpy as np # arange, cos, size, eye, inf, dot, floor, outer, zeros, linalg.eigh, sort, argsort, random, ones,...
from numpy import array, dot, log, sqrt, sum, exp, inf # to access the built-in sum fct: __builtins__.sum or del sum removes the imported sum and recovers the shadowed
try:
import matplotlib.pylab as pylab # also: use ipython -pylab
except:
pylab = None
print(' Could not import matplotlib.pylab, therefore ``cma.plot()`` etc. is not available')
__docformat__ = "reStructuredText" # this hides some comments entirely?
sys.py3kwarning = True # TODO: out-comment from version 2.6
# why not package math?
# TODO: check scitools.easyviz and how big the adaptation would be
# changes:
# 12/01/30: class Solution and more old stuff removed r3101
# 12/01/29: class Solution is depreciated, GenoPheno and SolutionDict do the job (v0.91.00, r3100)
# 12/01/06: CMA_eigenmethod option now takes a function (integer still works)
# 11/09/30: flat fitness termination checks also history length
# 11/09/30: elitist option (using method clip_or_fit_solutions)
# 11/09/xx: method clip_or_fit_solutions for check_points option for all sorts of
# injected or modified solutions and even reliable adaptive encoding
# 11/08/19: fixed: scaling and typical_x type clashes 1 vs array(1) vs ones(dim) vs dim * [1]
# 11/07/25: fixed: fmin wrote first and last line even with verb_log==0
# fixed: method settableOptionsList, also renamed to versatileOptions
# default seed depends on time now
# 11/07/xx (0.9.92): added: active CMA, selective mirrored sampling, noise/uncertainty handling
# fixed: output argument ordering in fmin, print now only used as function
# removed: parallel option in fmin
# 11/07/01: another try to get rid of the memory leak by replacing self.unrepaired = self[:]
# 11/07/01: major clean-up and reworking of abstract base classes and of the documentation,
# also the return value of fmin changed and attribute stop is now a method.
# 11/04/22: bug-fix: option fixed_variables in combination with scaling
# 11/04/21: stopdict is not a copy anymore
# 11/04/15: option fixed_variables implemented
# 11/03/23: bug-fix boundary update was computed even without boundaries
# 11/03/12: bug-fix of variable annotation in plots
# 11/02/05: work around a memory leak in numpy
# 11/02/05: plotting routines improved
# 10/10/17: cleaning up, now version 0.9.30
# 10/10/17: bug-fix: return values of fmin now use phenotyp (relevant
# if input scaling_of_variables is given)
# 08/10/01: option evalparallel introduced,
# bug-fix for scaling being a vector
# 08/09/26: option CMAseparable becomes CMA_diagonal
# 08/10/18: some names change, test functions go into a class
# 08/10/24: more refactorizing
# 10/03/09: upper bound exp(min(1,...)) for step-size control
# TODO: this would define the visible interface
# __all__ = ['fmin', 'CMAEvolutionStrategy', 'plot', ...]
#
# emptysets = ('', (), [], {}) # array([]) does not work but also np.size(.) == 0
# "x in emptysets" cannot be well replaced by "not x"
# which is also True for array([]) and None, but also for 0 and False, and False for NaN
use_sent_solutions = True # 5-30% CPU slower, particularly for large lambda, will be mandatory soon
#____________________________________________________________
#____________________________________________________________
#
def unitdoctest():
"""is used to describe test cases and might in future become helpful
as an experimental tutorial as well. The main testing feature at the
moment is by doctest with ``cma._test()`` or conveniently by
``python cma.py --test``. Unfortunately, depending on the
system, the results will slightly differ and many "failed" test cases
might be reported. This is prevented with the --quiet option.
A simple first overall test:
>>> import cma
>>> res = cma.fmin(cma.fcts.elli, 3*[1], 1, CMA_diagonal=2, seed=1, verb_time=0)
(3_w,7)-CMA-ES (mu_w=2.3,w_1=58%) in dimension 3 (seed=1)
Covariance matrix is diagonal for 2 iterations (1/ccov=7.0)
Iterat #Fevals function value axis ratio sigma minstd maxstd min:sec
1 7 1.453161670768570e+04 1.2e+00 1.08e+00 1e+00 1e+00
2 14 3.281197961927601e+04 1.3e+00 1.22e+00 1e+00 2e+00
3 21 1.082851071704020e+04 1.3e+00 1.24e+00 1e+00 2e+00
100 700 8.544042012075362e+00 1.4e+02 3.18e-01 1e-03 2e-01
200 1400 5.691152415221861e-12 1.0e+03 3.82e-05 1e-09 1e-06
220 1540 3.890107746209078e-15 9.5e+02 4.56e-06 8e-11 7e-08
termination on tolfun : 1e-11
final/bestever f-value = 3.89010774621e-15 2.52273602735e-15
mean solution: [ -4.63614606e-08 -3.42761465e-10 1.59957987e-11]
std deviation: [ 6.96066282e-08 2.28704425e-09 7.63875911e-11]
Test on the Rosenbrock function with 3 restarts. The first trial only
finds the local optimum, which happens in about 20% of the cases.
>>> import cma
>>> res = cma.fmin(cma.fcts.rosen, 4*[-1],1, ftarget=1e-6, restarts=3, verb_time=0, verb_disp=500, seed=3)
(4_w,8)-CMA-ES (mu_w=2.6,w_1=52%) in dimension 4 (seed=3)
Iterat #Fevals function value axis ratio sigma minstd maxstd min:sec
1 8 4.875315645656848e+01 1.0e+00 8.43e-01 8e-01 8e-01
2 16 1.662319948123120e+02 1.1e+00 7.67e-01 7e-01 8e-01
3 24 6.747063604799602e+01 1.2e+00 7.08e-01 6e-01 7e-01
184 1472 3.701428610430019e+00 4.3e+01 9.41e-07 3e-08 5e-08
termination on tolfun : 1e-11
final/bestever f-value = 3.70142861043 3.70142861043
mean solution: [-0.77565922 0.61309336 0.38206284 0.14597202]
std deviation: [ 2.54211502e-08 3.88803698e-08 4.74481641e-08 3.64398108e-08]
(8_w,16)-CMA-ES (mu_w=4.8,w_1=32%) in dimension 4 (seed=4)
Iterat #Fevals function value axis ratio sigma minstd maxstd min:sec
1 1489 2.011376859371495e+02 1.0e+00 8.90e-01 8e-01 9e-01
2 1505 4.157106647905128e+01 1.1e+00 8.02e-01 7e-01 7e-01
3 1521 3.548184889359060e+01 1.1e+00 1.02e+00 8e-01 1e+00
111 3249 6.831867555502181e-07 5.1e+01 2.62e-02 2e-04 2e-03
termination on ftarget : 1e-06
final/bestever f-value = 6.8318675555e-07 1.18576673231e-07
mean solution: [ 0.99997004 0.99993938 0.99984868 0.99969505]
std deviation: [ 0.00018973 0.00038006 0.00076479 0.00151402]
>>> assert res[1] <= 1e-6
Notice the different termination conditions. Termination on the target
function value ftarget prevents further restarts.
Test of scaling_of_variables option
>>> import cma
>>> opts = cma.Options()
>>> opts['seed'] = 456
>>> opts['verb_disp'] = 0
>>> opts['CMA_active'] = 1
>>> # rescaling of third variable: for searching in roughly
>>> # x0 plus/minus 1e3*sigma0 (instead of plus/minus sigma0)
>>> opts.scaling_of_variables = [1, 1, 1e3, 1]
>>> res = cma.fmin(cma.fcts.rosen, 4 * [0.1], 0.1, **opts)
termination on tolfun : 1e-11
final/bestever f-value = 2.68096173031e-14 1.09714829146e-14
mean solution: [ 1.00000001 1.00000002 1.00000004 1.00000007]
std deviation: [ 3.00466854e-08 5.88400826e-08 1.18482371e-07 2.34837383e-07]
The printed std deviations reflect the actual true value (not the one
in the internal representation which would be different).
:See: cma.main(), cma._test()
"""
pass
#____________________________________________________________
#____________________________________________________________
#
class BlancClass(object):
"""blanc container class for having a collection of attributes"""
#_____________________________________________________________________
#_____________________________________________________________________
#
class DerivedDictBase(collections.MutableMapping):
"""for conveniently adding features to a dictionary. The actual
dictionary is in ``self.data``. Copy-paste
and modify setitem, getitem, and delitem, if necessary"""
def __init__(self, *args, **kwargs):
collections.MutableMapping.__init__(self)
# super(SolutionDict, self).__init__() # the same
self.data = dict(*args, **kwargs)
def __len__(self):
return len(self.data)
def __contains__(self, value):
return value in self.data
def __iter__(self):
return iter(self.data)
def __setitem__(self, key, value):
"""defines self[key] = value"""
self.data[key] = value
def __getitem__(self, key):
"""defines self[key]"""
return self.data[key]
def __delitem__(self, key):
del self.data[key]
class SolutionDict(DerivedDictBase):
"""dictionary with computation of an hash key for the inserted solutions and
stack of previously inserted same solutions.
Each entry is meant to store additional information related to the solution.
d = SolutionDict()
x = array([1,2,4])
d[x] = {'x': x, 'iteration': 1}
d.get(x) == d[x] if d.key(x) in d.keys() else d.get(x) is None
TODO: not yet tested
TODO: data_with_same_key behaves like a stack (see setitem and delitem), but rather should behave like a queue?!
A queue is less consistent with the operation self[key] = ..., if self.data_with_same_key[key] is not empty.
"""
def __init__(self, *args, **kwargs):
DerivedDictBase.__init__(self, *args, **kwargs)
self.data_with_same_key = {}
def key(self, x):
try:
return tuple(x)
except TypeError:
return x
def __setitem__(self, key, value):
"""defines self[key] = value"""
key = self.key(key)
if key in self.data_with_same_key:
self.data_with_same_key[key] += [self.data[key]]
elif key in self.data:
self.data_with_same_key[key] = [self.data[key]]
self.data[key] = value
def __getitem__(self, key):
"""defines self[key]"""
return self.data[self.key(key)]
def __delitem__(self, key):
"""remove only most current key-entry"""
key = self.key(key)
if key in self.data_with_same_key:
if len(self.data_with_same_key[key]) == 1:
self.data[key] = self.data_with_same_key.pop(key)[0]
else:
self.data[key] = self.data_with_same_key[key].pop(-1)
else:
del self.data[key]
def truncate(self, max_len, min_iter):
if len(self) > max_len:
for k in self.keys():
if self[k]['iteration'] < min_iter:
del self[k] # only deletes one item with k as key, should delete all?
class SolutionDictOld(dict):
"""depreciated, SolutionDict should do, to be removed after SolutionDict
has been successfully applied.
dictionary with computation of an hash key for the inserted solutions and
stack of previously inserted same solutions.
Each entry is meant to store additional information related to the solution.
Methods ``pop`` and ``get`` are modified accordingly.
d = SolutionDict()
x = array([1,2,4])
d.insert(x, {'x': x, 'iteration': 1})
d.get(x) == d[d.key(x)] if d.key(x) in d.keys() else d.get(x) is None
TODO: not yet tested
TODO: behaves like a stack (see _pop_derived), but rather should behave like a queue?!
A queue is less consistent with the operation self[key] = ..., if self.more[key] is not empty.
"""
def __init__(self):
self.more = {} # previously inserted same solutions
self._pop_base = self.pop
self.pop = self._pop_derived
self._get_base = self.get
self.get = self._get_derived
def key(self, x):
"""compute the hash key of ``x``"""
return tuple(x)
def insert(self, x, datadict):
key = self.key(x)
if key in self.more:
self.more[key] += [self[key]]
elif key in self:
self.more[key] = [self[key]]
self[key] = datadict
def _get_derived(self, x, default=None):
return self._get_base(self.key(x), default)
def _pop_derived(self, x):
key = self.key(x)
res = self[key]
if key in self.more:
if len(self.more[key]) == 1:
self[key] = self.more.pop(key)[0]
else:
self[key] = self.more[key].pop(-1)
return res
class BestSolution(object):
"""container to keep track of the best solution seen"""
def __init__(self, x=None, f=np.inf, evals=None):
"""initialize the best solution with `x`, `f`, and `evals`.
Better solutions have smaller `f`-values.
"""
self.x = x
self.x_geno = None
self.f = f if f is not None and f is not np.nan else np.inf
self.evals = evals
self.evalsall = evals
self.last = BlancClass()
self.last.x = x
self.last.f = f
def update(self, arx, xarchive=None, arf=None, evals=None):
"""checks for better solutions in list `arx`, based on the smallest
corresponding value in `arf`, alternatively, `update` may be called
with a `BestSolution` instance like ``update(another_best_solution)``
in which case the better solution becomes the current best.
`xarchive` is used to retrieve the genotype of a solution.
"""
if arf is not None: # find failsave minimum
minidx = np.nanargmin(arf)
if minidx is np.nan:
return
minarf = arf[minidx]
# minarf = reduce(lambda x, y: y if y and y is not np.nan and y < x else x, arf, np.inf)
if type(arx) == BestSolution:
self.evalsall = max((self.evalsall, arx.evalsall))
if arx.f is not None and arx.f < np.inf:
self.update([arx.x], xarchive, [arx.f], arx.evals)
return self
elif minarf < np.inf and (minarf < self.f or self.f is None):
self.x, self.f = arx[minidx], arf[minidx]
self.x_geno = xarchive[self.x]['geno'] if xarchive is not None else None
self.evals = None if not evals else evals - len(arf) + minidx+1
self.evalsall = evals
elif evals:
self.evalsall = evals
self.last.x = arx[minidx]
self.last.f = minarf
def get(self):
"""return ``(x, f, evals)`` """
return self.x, self.f, self.evals, self.x_geno
#____________________________________________________________
#____________________________________________________________
#
class BoundPenalty(object):
"""Computes the boundary penalty. Must be updated each iteration,
using the `update` method.
Details
-------
The penalty computes like ``sum(w[i] * (x[i]-xfeas[i])**2)``,
where `xfeas` is the closest feasible (in-bounds) solution from `x`.
The weight `w[i]` should be updated during each iteration using
the update method.
This class uses `GenoPheno.into_bounds` in method `update` to access
domain boundary values and repair. This inconsistency might be
removed in future.
"""
def __init__(self, bounds=None):
"""Argument bounds can be `None` or ``bounds[0]`` and ``bounds[1]``
are lower and upper domain boundaries, each is either `None` or
a scalar or a list or array of appropriate size.
"""
##
# bounds attribute reminds the domain boundary values
self.bounds = bounds
self.gamma = 1 # a very crude assumption
self.weights_initialized = False # gamma becomes a vector after initialization
self.hist = [] # delta-f history
def has_bounds(self):
"""return True, if any variable is bounded"""
bounds = self.bounds
if bounds in (None, [None, None]):
return False
for i in xrange(bounds[0]):
if bounds[0][i] is not None and bound[0][i] > -np.inf:
return True
for i in xrange(bounds[1]):
if bounds[1][i] is not None and bound[1][i] < np.inf:
return True
return False
def repair(self, x, bounds=None, copy=False, copy_always=False):
"""sets out-of-bounds components of ``x`` on the bounds.
Arguments
---------
`bounds`
can be `None`, in which case the "default" bounds are used,
or ``[lb, ub]``, where `lb` and `ub`
represent lower and upper domain bounds respectively that
can be `None` or a scalar or a list or array of length ``len(self)``
code is more or less copy-paste from Solution.repair, but never tested
"""
# TODO (old data): CPU(N,lam,iter=20,200,100): 3.3s of 8s for two bounds, 1.8s of 6.5s for one bound
# TODO: test whether np.max([bounds[0], x], axis=0) etc is speed relevant
if bounds is None:
bounds = self.bounds
if copy_always:
x_out = array(x, copy=True)
if bounds not in (None, [None, None], (None, None)): # solely for effiency
x_out = array(x, copy=True) if copy and not copy_always else x
if bounds[0] is not None:
if np.isscalar(bounds[0]):
for i in xrange(len(x)):
x_out[i] = max([bounds[0], x[i]])
else:
for i in xrange(len(x)):
if bounds[0][i] is not None:
x_out[i] = max([bounds[0][i], x[i]])
if bounds[1] is not None:
if np.isscalar(bounds[1]):
for i in xrange(len(x)):
x_out[i] = min([bounds[1], x[i]])
else:
for i in xrange(len(x)):
if bounds[1][i] is not None:
x_out[i] = min([bounds[1][i], x[i]])
return x_out # convenience return
#____________________________________________________________
#
def __call__(self, x, archive, gp):
"""returns the boundary violation penalty for `x` ,where `x` is a
single solution or a list or array of solutions.
If `bounds` is not `None`, the values in `bounds` are used, see `__init__`"""
if x in (None, (), []):
return x
if gp.bounds in (None, [None, None], (None, None)):
return 0.0 if np.isscalar(x[0]) else [0.0] * len(x) # no penalty
x_is_single_vector = np.isscalar(x[0])
x = [x] if x_is_single_vector else x
pen = []
for xi in x:
# CAVE: this does not work with already repaired values!!
# CPU(N,lam,iter=20,200,100)?: 3s of 10s, array(xi): 1s (check again)
# remark: one deep copy can be prevented by xold = xi first
xpheno = gp.pheno(archive[xi]['geno'])
xinbounds = gp.into_bounds(xpheno)
fac = 1 # exp(0.1 * (log(self.scal) - np.mean(self.scal)))
pen.append(sum(self.gamma * ((xinbounds - xpheno) / fac)**2) / len(xi))
return pen[0] if x_is_single_vector else pen
#____________________________________________________________
#
def feasible_ratio(self, solutions):
"""counts for each coordinate the number of feasible values in
``solutions`` and returns an array of length ``len(solutions[0])``
with the ratios.
`solutions` is a list or array of repaired `Solution` instances
"""
count = np.zeros(len(solutions[0]))
for x in solutions:
count += x.unrepaired == x
return count / float(len(solutions))
#____________________________________________________________
#
def update(self, function_values, es, bounds=None):
"""updates the weights for computing a boundary penalty.
Arguments
---------
`function_values`
all function values of recent population of solutions
`es`
`CMAEvolutionStrategy` object instance, in particular the
method `into_bounds` of the attribute `gp` of type `GenoPheno`
is used.
`bounds`
not (yet) in use other than for ``bounds == [None, None]`` nothing
is updated.
Reference: Hansen et al 2009, A Method for Handling Uncertainty...
IEEE TEC, with addendum at http://www.lri.fr/~hansen/TEC2009online.pdf
"""
if bounds is None:
bounds = self.bounds
if bounds is None or (bounds[0] is None and bounds[1] is None): # no bounds ==> no penalty
return self # len(function_values) * [0.0] # case without voilations
N = es.N
### prepare
# compute varis = sigma**2 * C_ii
varis = es.sigma**2 * array(N * [es.C] if np.isscalar(es.C) else ( # scalar case
es.C if np.isscalar(es.C[0]) else # diagonal matrix case
[es.C[i][i] for i in xrange(N)])) # full matrix case
dmean = (es.mean - es.gp.into_bounds(es.mean)) / varis**0.5
### Store/update a history of delta fitness value
fvals = sorted(function_values)
l = 1 + len(fvals)
val = fvals[3*l // 4] - fvals[l // 4] # exact interquartile range apart interpolation
val = val / np.mean(varis) # new: val is normalized with sigma of the same iteration
# insert val in history
if np.isfinite(val) and val > 0:
self.hist.insert(0, val)
elif val == inf and len(self.hist) > 1:
self.hist.insert(0, max(self.hist))
else:
pass # ignore 0 or nan values
if len(self.hist) > 20 + (3*N) / es.popsize:
self.hist.pop()
### prepare
dfit = np.median(self.hist) # median interquartile range
damp = min(1, es.sp.mueff/10./N)
### set/update weights
# Throw initialization error
if len(self.hist) == 0:
raise _Error('wrongful initialization, no feasible solution sampled. ' +
'Reasons can be mistakenly set bounds (lower bound not smaller than upper bound) or a too large initial sigma0 or... ' +
'See description of argument func in help(cma.fmin) or an example handling infeasible solutions in help(cma.CMAEvolutionStrategy). ')
# initialize weights
if (dmean.any() and (not self.weights_initialized or es.countiter == 2)): # TODO
self.gamma = array(N * [2*dfit])
self.weights_initialized = True
# update weights gamma
if self.weights_initialized:
edist = array(abs(dmean) - 3 * max(1, N**0.5/es.sp.mueff))
if 1 < 3: # this is better, around a factor of two
# increase single weights possibly with a faster rate than they can decrease
# value unit of edst is std dev, 3==random walk of 9 steps
self.gamma *= exp((edist>0) * np.tanh(edist/3) / 2.)**damp
# decrease all weights up to the same level to avoid single extremely small weights
# use a constant factor for pseudo-keeping invariance
self.gamma[self.gamma > 5 * dfit] *= exp(-1./3)**damp
# self.gamma[idx] *= exp(5*dfit/self.gamma[idx] - 1)**(damp/3)
elif 1 < 3 and (edist>0).any(): # previous method
# CAVE: min was max in TEC 2009
self.gamma[edist>0] *= 1.1**min(1, es.sp.mueff/10./N)
# max fails on cigtab(N=12,bounds=[0.1,None]):
# self.gamma[edist>0] *= 1.1**max(1, es.sp.mueff/10./N) # this was a bug!?
# self.gamma *= exp((edist>0) * np.tanh(edist))**min(1, es.sp.mueff/10./N)
else: # alternative version, but not better
solutions = es.pop # this has not been checked
r = self.feasible_ratio(solutions) # has to be the averaged over N iterations
self.gamma *= exp(np.max([N*[0], 0.3 - r], axis=0))**min(1, es.sp.mueff/10/N)
es.more_to_write = self.gamma if self.weights_initialized else np.ones(N)
### return penalty
# es.more_to_write = self.gamma if not np.isscalar(self.gamma) else N*[1]
return self # bound penalty values
#____________________________________________________________
#____________________________________________________________
#
class GenoPhenoBase(object):
"""depreciated, abstract base class for genotyp-phenotype transformation,
to be implemented.
See (and rather use) option ``transformation`` of ``fmin`` or ``CMAEvolutionStrategy``.
Example
-------
::
import cma
class Mygpt(cma.GenoPhenoBase):
def pheno(self, x):
return x # identity for the time being
gpt = Mygpt()
optim = cma.CMAEvolutionStrategy(...)
while not optim.stop():
X = optim.ask()
f = [func(gpt.pheno(x)) for x in X]
optim.tell(X, f)
In case of a repair, we might pass the repaired solution into `tell()`
(with check_points being True).
TODO: check usecases in `CMAEvolutionStrategy` and implement option GenoPhenoBase
"""
def pheno(self, x):
raise NotImplementedError()
return x
#____________________________________________________________
#____________________________________________________________
#
class GenoPheno(object):
"""Genotype-phenotype transformation.
Method `pheno` provides the transformation from geno- to phenotype,
that is from the internal representation to the representation used
in the objective function. Method `geno` provides the "inverse" pheno-
to genotype transformation. The geno-phenotype transformation comprises,
in this order:
- insert fixed variables (with the phenotypic and therefore quite
possibly "wrong" values)
- affine linear transformation (scaling and shift)
- user-defined transformation
- projection into feasible domain (boundaries)
- assign fixed variables their original phenotypic value
By default all transformations are the identity. The boundary
transformation is only applied, if the boundaries are given as argument to
the method `pheno` or `geno` respectively.
``geno`` is not really necessary and might disappear in future.
"""
def __init__(self, dim, scaling=None, typical_x=None, bounds=None, fixed_values=None, tf=None):
"""return `GenoPheno` instance with fixed dimension `dim`.
Keyword Arguments
-----------------
`scaling`
the diagonal of a scaling transformation matrix, multipliers
in the genotyp-phenotyp transformation, see `typical_x`
`typical_x`
``pheno = scaling*geno + typical_x``
`bounds` (obsolete, might disappear)
list with two elements,
lower and upper bounds both can be a scalar or a "vector"
of length dim or `None`. Without effect, as `bounds` must
be given as argument to `pheno()`.
`fixed_values`
a dictionary of variable indices and values, like ``{0:2.0, 2:1.1}``,
that a not subject to change, negative indices are dropped
(they act like incommenting the index), values are phenotypic
values.
`tf`
list of two user-defined transformation functions, or `None`.
``tf[0]`` is a function that transforms the internal representation
as used by the optimizer into a solution as used by the
objective function. ``tf[1]`` does the back-transformation.
For example ::
tf_0 = lambda x: [xi**2 for xi in x]
tf_1 = lambda x: [abs(xi)**0.5 fox xi in x]
or "equivalently" without the `lambda` construct ::
def tf_0(x):
return [xi**2 for xi in x]
def tf_1(x):
return [abs(xi)**0.5 fox xi in x]
``tf=[tf_0, tf_1]`` is a reasonable way to guaranty that only positive
values are used in the objective function.
Details
-------
In future it might be possible to omit `tf_1`
"""
self.N = dim
self.bounds = bounds
self.fixed_values = fixed_values
if tf is not None:
self.tf_pheno = tf[0]
self.tf_geno = tf[1] # TODO: should not necessarily be needed
# r = np.random.randn(dim)
# assert all(tf[0](tf[1](r)) - r < 1e-7)
# r = np.random.randn(dim)
# assert all(tf[0](tf[1](r)) - r > -1e-7)
print("WARNING in class GenoPheno: user defined transformations have not been tested thoroughly")
else:
self.tf_geno = None
self.tf_pheno = None
if fixed_values:
if type(fixed_values) is not dict:
raise _Error("fixed_values must be a dictionary {index:value,...}")
if max(fixed_values.keys()) >= dim:
raise _Error("max(fixed_values.keys()) = " + str(max(fixed_values.keys())) +
" >= dim=N=" + str(dim) + " is not a feasible index")
# convenience commenting functionality: drop negative keys
for k in fixed_values.keys():
if k < 0:
fixed_values.pop(k)
if bounds:
if len(bounds) != 2:
raise _Error('len(bounds) must be 2 for lower and upper bounds')
for i in (0,1):
if bounds[i] is not None:
bounds[i] = array(dim * [bounds[i]] if np.isscalar(bounds[i]) else
[b for b in bounds[i]])
def vec_is_default(vec, default_val=0):
"""None or [None] are also recognized as default"""
try:
if len(vec) == 1:
vec = vec[0] # [None] becomes None and is always default
else:
return False
except TypeError:
pass # vec is a scalar
if vec is None or vec == array(None) or vec == default_val:
return True
return False
self.scales = array(scaling)
if vec_is_default(self.scales, 1):
self.scales = 1 # CAVE: 1 is not array(1)
elif self.scales.shape is not () and len(self.scales) != self.N:
raise _Error('len(scales) == ' + str(len(self.scales)) +
' does not match dimension N == ' + str(self.N))
self.typical_x = array(typical_x)
if vec_is_default(self.typical_x, 0):
self.typical_x = 0
elif self.typical_x.shape is not () and len(self.typical_x) != self.N:
raise _Error('len(typical_x) == ' + str(len(self.typical_x)) +
' does not match dimension N == ' + str(self.N))
if (self.scales is 1 and
self.typical_x is 0 and
self.bounds in (None, [None, None]) and
self.fixed_values is None and
self.tf_pheno is None):
self.isidentity = True
else:
self.isidentity = False
def into_bounds(self, y, bounds=None, copy_never=False, copy_always=False):
"""Argument `y` is a phenotypic vector,
return `y` put into boundaries, as a copy iff ``y != into_bounds(y)``.
Note: this code is duplicated in `Solution.repair` and might
disappear in future.
"""
bounds = bounds if bounds is not None else self.bounds
if bounds in (None, [None, None]):
return y if not copy_always else array(y, copy=True)
if bounds[0] is not None:
if len(bounds[0]) not in (1, len(y)):
raise ValueError('len(bounds[0]) = ' + str(len(bounds[0])) +
' and len of initial solution (' + str(len(y)) + ') disagree')
if copy_never: # is rather slower
for i in xrange(len(y)):
y[i] = max(bounds[0][i], y[i])
else:
y = np.max([bounds[0], y], axis=0)
if bounds[1] is not None:
if len(bounds[1]) not in (1, len(y)):
raise ValueError('len(bounds[1]) = ' + str(len(bounds[1])) +
' and initial solution (' + str(len(y)) + ') disagree')
if copy_never:
for i in xrange(len(y)):
y[i] = min(bounds[1][i], y[i])
else:
y = np.min([bounds[1], y], axis=0)
return y
def pheno(self, x, bounds=None, copy=True, copy_always=False):
"""maps the genotypic input argument into the phenotypic space,
boundaries are only applied if argument ``bounds is not None``, see
help for class `GenoPheno`
"""
if copy_always and not copy:
raise ValueError('arguments copy_always=' + str(copy_always) +
' and copy=' + str(copy) + ' have inconsistent values')
if self.isidentity and bounds in (None, [None, None], (None, None)):
return x if not copy_always else array(x, copy=copy_always)
if self.fixed_values is None:
y = array(x, copy=copy) # make a copy, in case
else: # expand with fixed values
y = list(x) # is a copy
for i in sorted(self.fixed_values.keys()):
y.insert(i, self.fixed_values[i])
y = array(y, copy=False)
if self.scales is not 1: # just for efficiency
y *= self.scales
if self.typical_x is not 0:
y += self.typical_x
if self.tf_pheno is not None:
y = array(self.tf_pheno(y), copy=False)
if bounds is not None:
y = self.into_bounds(y, bounds)
if self.fixed_values is not None:
for i, k in self.fixed_values.items():
y[i] = k
return y
def geno(self, y, bounds=None, copy=True, copy_always=False, archive=None):
"""maps the phenotypic input argument into the genotypic space.
If `bounds` are given, first `y` is projected into the feasible
domain. In this case ``copy==False`` leads to a copy.
by default a copy is made only to prevent to modify ``y``
method geno is only needed if external solutions are injected
(geno(initial_solution) is depreciated and will disappear)
TODO: arg copy=True should become copy_never=False
"""
if archive is not None and bounds is not None:
try:
return archive[y]['geno']
except:
pass
x = array(y, copy=(copy and not self.isidentity) or copy_always)
# bounds = self.bounds if bounds is None else bounds
if bounds is not None: # map phenotyp into bounds first
x = self.into_bounds(x, bounds)
if self.isidentity:
return x
# user-defined transformation
if self.tf_geno is not None:
x = array(self.tf_geno(x), copy=False)
# affine-linear transformation: shift and scaling
if self.typical_x is not 0:
x -= self.typical_x
if self.scales is not 1: # just for efficiency
x /= self.scales
# kick out fixed_values
if self.fixed_values is not None:
# keeping the transformed values does not help much
# therefore it is omitted
if 1 < 3:
keys = sorted(self.fixed_values.keys())
x = array([x[i] for i in range(len(x)) if i not in keys], copy=False)
else: # TODO: is this more efficient?
x = list(x)
for key in sorted(self.fixed_values.keys(), reverse=True):
x.remove(key)
x = array(x, copy=False)
return x
#____________________________________________________________
#____________________________________________________________
# check out built-in package abc: class ABCMeta, abstractmethod, abstractproperty...
# see http://docs.python.org/whatsnew/2.6.html PEP 3119 abstract base classes
#
class OOOptimizer(object):
""""abstract" base class for an OO optimizer interface with methods
`__init__`, `ask`, `tell`, `stop`, `result`, and `optimize`. Only
`optimize` is fully implemented in this base class.
Examples
--------
All examples minimize the function `elli`, the output is not shown.
(A preferred environment to execute all examples is ``ipython -pylab``.)
First we need ::
from cma import CMAEvolutionStrategy, CMADataLogger # CMAEvolutionStrategy derives from the OOOptimizer class
elli = lambda x: sum(1e3**((i-1.)/(len(x)-1.)*x[i])**2 for i in range(len(x)))
The shortest example uses the inherited method `OOOptimizer.optimize()`::
res = CMAEvolutionStrategy(8 * [0.1], 0.5).optimize(elli)
The input parameters to `CMAEvolutionStrategy` are specific to this
inherited class. The remaining functionality is based on interface
defined by `OOOptimizer`. We might have a look at the result::
print(res[0]) # best solution and
print(res[1]) # its function value
`res` is the return value from method
`CMAEvolutionStrategy.result()` appended with `None` (no logger).
In order to display more exciting output we rather do ::
logger = CMADataLogger() # derives from the abstract BaseDataLogger class
res = CMAEvolutionStrategy(9 * [0.5], 0.3).optimize(elli, logger)
logger.plot() # if matplotlib is available, logger == res[-1]
or even shorter ::
res = CMAEvolutionStrategy(9 * [0.5], 0.3).optimize(elli, CMADataLogger())
res[-1].plot() # if matplotlib is available
Virtually the same example can be written with an explicit loop
instead of using `optimize()`. This gives the necessary insight into
the `OOOptimizer` class interface and gives entire control over the
iteration loop::
optim = CMAEvolutionStrategy(9 * [0.5], 0.3) # a new CMAEvolutionStrategy instance calling CMAEvolutionStrategy.__init__()
logger = CMADataLogger(optim) # get a logger instance
# this loop resembles optimize()
while not optim.stop(): # iterate
X = optim.ask() # get candidate solutions
f = [elli(x) for x in X] # evaluate solutions
# maybe do something else that needs to be done
optim.tell(X, f) # do all the real work: prepare for next iteration
optim.disp(20) # display info every 20th iteration
logger.add() # log another "data line"
# final output
print('termination by', optim.stop())
print('best f-value =', optim.result()[1])
print('best solution =', optim.result()[0])
logger.plot() # if matplotlib is available
raw_input('press enter to continue') # prevents exiting and closing figures
Details
-------
Most of the work is done in the method `tell(...)`. The method `result()` returns
more useful output.
"""
@staticmethod
def abstract():
"""marks a method as abstract, ie to be implemented by a subclass"""
import inspect
caller = inspect.getouterframes(inspect.currentframe())[1][3]
raise NotImplementedError('method ' + caller + '() must be implemented in subclass')
def __init__(self, xstart, **more_args):
"""abstract method, ``xstart`` is a mandatory argument"""
OOOptimizer.abstract()
def initialize(self):
"""(re-)set to the initial state"""
OOOptimizer.abstract()
def ask(self):
"""abstract method, AKA "get", deliver new candidate solution(s), a list of "vectors"
"""
OOOptimizer.abstract()
def tell(self, solutions, function_values):
"""abstract method, AKA "update", prepare for next iteration"""
OOOptimizer.abstract()
def stop(self):
"""abstract method, return satisfied termination conditions in a dictionary like
``{'termination reason': value, ...}``, for example ``{'tolfun': 1e-12}``, or the empty
dictionary ``{}``. The implementation of `stop()` should prevent an infinite loop.
"""
OOOptimizer.abstract()
def disp(self, modulo=None):
"""abstract method, display some iteration infos if ``self.iteration_counter % modulo == 0``"""
OOOptimizer.abstract()
def result(self):
"""abstract method, return ``(x, f(x), ...)``, that is, the minimizer, its function value, ..."""
OOOptimizer.abstract()
def optimize(self, objectivefct, logger=None, verb_disp=20, iterations=None):
"""find minimizer of `objectivefct` by iterating over `OOOptimizer` `self`
with verbosity `verb_disp`, using `BaseDataLogger` `logger` with at
most `iterations` iterations. ::
return self.result() + (self.stop(), self, logger)
Example
-------
>>> import cma
>>> res = cma.CMAEvolutionStrategy(7 * [0.1], 0.5).optimize(cma.fcts.rosen, None, 100)
(4_w,9)-CMA-ES (mu_w=2.8,w_1=49%) in dimension 7 (seed=630721393)
Iterat #Fevals function value axis ratio sigma minstd maxstd min:sec
1 9 3.163954777181882e+01 1.0e+00 4.12e-01 4e-01 4e-01 0:0.0
2 18 3.299006223906629e+01 1.0e+00 3.60e-01 3e-01 4e-01 0:0.0
3 27 1.389129389866704e+01 1.1e+00 3.18e-01 3e-01 3e-01 0:0.0
100 900 2.494847340045985e+00 8.6e+00 5.03e-02 2e-02 5e-02 0:0.3
200 1800 3.428234862999135e-01 1.7e+01 3.77e-02 6e-03 3e-02 0:0.5
300 2700 3.216640032470860e-04 5.6e+01 6.62e-03 4e-04 9e-03 0:0.8
400 3600 6.155215286199821e-12 6.6e+01 7.44e-06 1e-07 4e-06 0:1.1
438 3942 1.187372505161762e-14 6.0e+01 3.27e-07 4e-09 9e-08 0:1.2
438 3942 1.187372505161762e-14 6.0e+01 3.27e-07 4e-09 9e-08 0:1.2
('termination by', {'tolfun': 1e-11})
('best f-value =', 1.1189867885201275e-14)
('solution =', array([ 1. , 1. , 1. , 0.99999999, 0.99999998,
0.99999996, 0.99999992]))
>>> print(res[0])
[ 1. 1. 1. 0.99999999 0.99999998 0.99999996
0.99999992]
"""
if logger is None and hasattr(self, 'logger'):
logger = self.logger
citer = 0
while not self.stop():
if iterations is not None and citer >= iterations:
return self.result()
citer += 1
X = self.ask() # deliver candidate solutions
fitvals = [objectivefct(x) for x in X]
self.tell(X, fitvals) # all the work is done here
self.disp(verb_disp)
logger.add(self) if logger else None
logger.add(self, modulo=bool(logger.modulo)) if logger else None
if verb_disp:
self.disp(1)
print('termination by', self.stop())
print('best f-value =', self.result()[1])
print('solution =', self.result()[0])
return self.result() + (self.stop(), self, logger)
#____________________________________________________________
#____________________________________________________________
#
class CMAEvolutionStrategy(OOOptimizer):
"""CMA-ES stochastic optimizer class with ask-and-tell interface.
See `fmin` for the one-line-call functional interface.
Calling sequence
================
``optim = CMAEvolutionStrategy(x0, sigma0, opts)``
returns a class instance.
Arguments
---------
`x0`
initial solution, starting point.
`sigma0`
initial standard deviation. The problem
variables should have been scaled, such that a single
standard deviation on all variables is useful and the
optimum is expected to lie within about `x0` +- ``3*sigma0``.
See also options `scaling_of_variables`.
Often one wants to check for solutions close to the initial
point. This allows for an easier check for consistency of
the objective function and its interfacing with the optimizer.
In this case, a much smaller `sigma0` is advisable.
`opts`
options, a dictionary with optional settings,
see class `Options`.
Main interface / usage
======================
The ask-and-tell interface is inherited from the generic `OOOptimizer`
interface for iterative optimization algorithms (see there). With ::
optim = CMAEvolutionStrategy(8 * [0.5], 0.2)
an object instance is generated. In each iteration ::
solutions = optim.ask()
is used to ask for new candidate solutions (possibly several times) and ::
optim.tell(solutions, func_values)
passes the respective function values to `optim`. Instead of `ask()`,
the class `CMAEvolutionStrategy` also provides ::
(solutions, func_values) = optim.ask_and_eval(objective_func)
Therefore, after initialization, an entire optimization can be written
in two lines like ::
while not optim.stop():
optim.tell(*optim.ask_and_eval(objective_func))
Without the freedom of executing additional lines within the iteration,
the same reads in a single line as ::
optim.optimize(objective_func)
Besides for termination criteria, in CMA-ES only
the ranks of the `func_values` are relevant.
Attributes and Properties
=========================
- `inputargs` -- passed input arguments
- `inopts` -- passed options
- `opts` -- actually used options, some of them can be changed any
time, see class `Options`
- `popsize` -- population size lambda, number of candidate solutions
returned by `ask()`
Details
=======
The following two enhancements are turned off by default.
**Active CMA** is implemented with option ``CMA_active`` and conducts
an update of the covariance matrix with negative weights. The
exponential update is implemented, where from a mathematical
viewpoint positive definiteness is guarantied. The update is applied
after the default update and only before the covariance matrix is
decomposed, which limits the additional computational burden to be
at most a factor of three (typically smaller). A typical speed up
factor (number of f-evaluations) is between 1.1 and two.
References: Jastrebski and Arnold, CEC 2006, Glasmachers et al, GECCO 2010.
**Selective mirroring** is implemented with option ``CMA_mirrors`` in
the method ``get_mirror()``. Only the method `ask_and_eval()` will
then sample selectively mirrored vectors. In selective mirroring, only
the worst solutions are mirrored. With the default small number of mirrors,
*pairwise selection* (where at most one of the two mirrors contribute to the
update of the distribution mean) is implicitely guarantied under selective
mirroring and therefore not explicitly implemented.
References: Brockhoff et al, PPSN 2010, Auger et al, GECCO 2011.
Examples
========
Super-short example, with output shown:
>>> import cma
>>> # construct an object instance in 4-D, sigma0=1
>>> es = cma.CMAEvolutionStrategy(4 * [1], 1, {'seed':234})
(4_w,8)-CMA-ES (mu_w=2.6,w_1=52%) in dimension 4 (seed=234)
>>>
>>> # iterate until termination
>>> while not es.stop():
... X = es.ask()
... es.tell(X, [cma.fcts.elli(x) for x in X])
... es.disp() # by default sparse, see option verb_disp
Iterat #Fevals function value axis ratio sigma minstd maxstd min:sec
1 8 2.093015112685775e+04 1.0e+00 9.27e-01 9e-01 9e-01 0:0.0
2 16 4.964814235917688e+04 1.1e+00 9.54e-01 9e-01 1e+00 0:0.0
3 24 2.876682459926845e+05 1.2e+00 1.02e+00 9e-01 1e+00 0:0.0
100 800 6.809045875281943e-01 1.3e+02 1.41e-02 1e-04 1e-02 0:0.2
200 1600 2.473662150861846e-10 8.0e+02 3.08e-05 1e-08 8e-06 0:0.5
233 1864 2.766344961865341e-14 8.6e+02 7.99e-07 8e-11 7e-08 0:0.6
>>>
>>> cma.pprint(es.result())
(Solution([ -1.98546755e-09, -1.10214235e-09, 6.43822409e-11,
-1.68621326e-11]),
4.5119610261406537e-16,
1666,
1672,
209,
array([ -9.13545269e-09, -1.45520541e-09, -6.47755631e-11,
-1.00643523e-11]),
array([ 3.20258681e-08, 3.15614974e-09, 2.75282215e-10,
3.27482983e-11]))
>>>
>>> # help(es.result) shows
result(self) method of cma.CMAEvolutionStrategy instance
return ``(xbest, f(xbest), evaluations_xbest, evaluations, iterations, pheno(xmean), effective_stds)``
Using the multiprocessing module, we can evaluate the function in parallel with a simple
modification of the example ::
import multiprocessing
# prepare es = ...
pool = multiprocessing.Pool(es.popsize)
while not es.stop():
X = es.ask()
es.tell(X, pool.map_async(cma.fcts.elli, X))
Example with a data logger, lower bounds (at zero) and handling infeasible solutions:
>>> import cma
>>> import numpy as np
>>> es = cma.CMAEvolutionStrategy(10 * [0.2], 0.5, {'bounds': [0, np.inf]})
>>> logger = cma.CMADataLogger().register(es)
>>> while not es.stop():
... fit, X = [], []
... while len(X) < es.popsize:
... curr_fit = np.NaN
... while curr_fit is np.NaN:
... x = es.ask(1)[0]
... curr_fit = cma.fcts.somenan(x, cma.fcts.elli) # might return np.NaN
... X.append(x)
... fit.append(curr_fit)
... es.tell(X, fit)
... logger.add()
... es.disp()