#!/usr/bin/env python
"""Module cma implements the CMA-ES, Covariance Matrix Adaptation Evolution 
Strategy, a stochastic optimizer for robust non-linear non-convex function 
minimization. 

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. 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-2011

: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 print_function  # for cross-checking, available from python 2.6

__version__ = "0.9.96 $Revision: 3031 $" 
#    $Date: 2011-09-27 23:11:32 +0200 (Tue, 27 Sep 2011) $
#    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 <http://www.gnu.org/licenses/>.
# 

# 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: split tell into a variable transformation part and the "pure" functionality
#       usecase: es.tell(X, [func(es.pheno(x)) for x in X])
# 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}, <allvecs>)
#
# 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   # 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
    matplotlib.use('wxagg')
    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:
# 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

#____________________________________________________________
#____________________________________________________________
#
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 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.f = f if f is not None else np.inf
        self.evals = evals
        self.evalsall = evals
    def update(self, arx, 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 
        
        """
        if arf is not None:  # find failsave minimum
            minidx = np.nanargmin(arf)	    
            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], [arx.f], arx.evals)
        elif minarf < np.inf and (minarf < self.f or self.f is None):
            self.x, self.f = arx[minidx], arf[minidx]
            self.evals = None if not evals else evals - len(arf) + minidx+1
            self.evalsall = evals
        elif evals:
            self.evalsall = evals
    def get(self):
        """return ``(x, f, evals)`` """
        return self.x, self.f, self.evals
    
    
#____________________________________________________________
#____________________________________________________________
#
class Solution(np.ndarray):
    """this class is a stump, but functional and in use: Solution represents 
    a candidate solution (a vector of real numbers), for example returned by 
    `cma.CMAEvolutionStrategy.ask()`. 
    The class inherits from `numpy.ndarray` and can therefore be used as 
    a `numpy` array. Additionally, the solution can be repaired based on domain 
    boundaries. The input arguments to further available methods are subject to 
    future changes.  
    
    Attributes
    ----------
        - unrepaired 
            array, solution before the last call to repair() 
            in case of bounds, or current solution in case bounds are None, 
            or None in case repair was never called
        - isgeno 
            whether the solution represents a genotyp (internal
            representation in CMAEvolutionStrategy)
        - geno 
            array, solution during construction, genotype 
            representation, or None

    TODO: add attributes for the f-value f and the nb of f-evaluations 
    evals? 
    
    Example
    -------
     
        >>> import cma
        >>> x = cma.Solution([0, -2, 3])
        >>> x                         # is a float np.array
        Solution([ 0., -2.,  3.])
        >>> x.repair([0.5, None])     # set lower domain bound to 0.5
        Solution([ 0.5,  0.5,  3. ])
        >>> x                         # check that x was changed
        Solution([ 0.5,  0.5,  3. ])
        >>> x.unrepaired              # original values
        array([ 0., -2.,  3.])
        >>> x.isgeno                  # by default not the internal repr.
        False
        
    """
    def __init__(self, x, copy=False, geno=False, gp=None):
        """return `Solution` instance with value `x`. 
        
        Arguments
        ---------
            `x`
                data vector
            `copy` 
                make a copy of `x`?
            `geno` 
                result represents a genotype?
            `gp` 
                if not `None` and ``geno==False`` a `cma.GenoPheno` object 
                used to transform `x` to the phenotype
                
        """ 
        # f=np.NaN, evals=np.NaN, ldom=None, udom=None):
        # by now __new__(self, x, copy=False, geno=False) has been
        # called implicitly, so there is nothing to do. 
        self.isgeno = geno  # overwrite values from x
        self.geno = None
        if geno is False and gp is not None:  
            self.geno = x  # array(x, copy)
        self.unrepaired = None # array(self, copy=False)
        return
        
    def __new__(cls, data, copy=False, geno=False, gp=None):
        """return a new object instance of type `cls` == `cma.Solution` """ 
        # static classmethod, implicitly called before __init__
        # cls is like the self argument for an instance!?
        # necessary, because numpy.ndarray behaves in part like a
        # built-in class (what a hack)
        
        # initialize self from argument data
        if type(data) == cls and copy is False and gp is None:
            return data
        
        if geno is False and gp is not None:
            self = np.array(gp.pheno(data, copy=True), dtype=np.float_, copy=False).view(type=cls)
            # self = np.asarray(gp.pheno(data, bounds)).view(cls)
            # self = np.ndarray.__new__(Solution, np.asarray(data).shape)
            # self[:] = gp.pheno(data, bounds)
        else:
            self = np.array(data, dtype=np.float_, copy=copy).view(type=cls)
            # self = np.asarray(data).view(cls)
            # self = np.ndarray.__new__(Solution, np.asarray(data).shape)
            # self[:] = np.asarray(data)
        
        # initialize attributes (moved to init)
        if 11 < 3:  # disregard class of data and type(data) is not Solution: 
            self.isgeno = geno  
            self.unrepaired = None
            
        return self  # calls _array_finalize__, if available
        
    #~ def __array_finalize__(self, obj):  # only called if type(obj) is Solution?
        #~ # Purpose: assign additional attributes from initializing obj (copy constructor)
        #~ # This is removed, because it is not very efficient.
        #~ # the effect: any new solution is not isgeno and has no unrepaired version
        #~ if type(obj) is Solution:  # copy attributes
            #~ self.unrepaired = getattr(obj, 'unrepaired', None)
            #~ self.isgeno = getattr(obj, 'isgeno', None)
    
    def get_geno(self, gp):
        """return solution in "internal" representation.
        Argument `gp`: `GenoPheno` class instance (subject to future changes)
        """

        if self.isgeno:
            return self
        # gp.geno copies data
        return Solution(gp.geno(self), copy=False, geno=True) 
        
    def get_pheno(self, gp):
        """return solution in original representation. 
        
        TODO: consider name change
        
        Argument `gp` -- `GenoPheno` class instance (subject to future changes)
        
        """
        if not self.isgeno:
            return self

        # keep genotype, might become out-of-date
        self.geno = array(self)  

        # gp.pheno copies data
        return Solution(gp.pheno(self), copy=False, geno=False)
    
    def repair(self, bounds, keep_unrepaired=True):
        """sets out-of-bounds components of the represented solution 
        on the bounds.  
        
        Arguments
        ---------
            `bounds` 
                can be `None` 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)``
            `keep_unrepaired`    
                make a copy of the original values in case of bounds
        """
        # TODO: CPU(N,lam,iter=20,200,100): 3.3s of 8s for two bounds, 1.8s of 6.5s for one bound
        if self.isgeno:
            raise _Error('genotyp cannot be repaired')
            
        if bounds not in (None, [None, None]):  # solely for effiency 
            # this produces a memory leak (supposably a numpy bug, reported by Daan Sprunken)
            # self.unrepaired = self.view(np.ndarray)  
                    
            self.unrepaired = None
            if keep_unrepaired: 
                # self.unrepaired = self[:]  # make a copy
                self.unrepaired = array(self)  # make a copy, resolves the memory leak
            if bounds[0] is not None:
                if np.isscalar(bounds[0]):
                    for i in xrange(len(self)):
                        self[i] = max(bounds[0], self[i])
                else:
                    for i in xrange(len(self)):
                        self[i] = max(bounds[0][i], self[i])
            if bounds[1] is not None:
                if np.isscalar(bounds[1]):
                    for i in xrange(len(self)):
                        self[i] = min(bounds[1], self[i])
                else:
                    for i in xrange(len(self)):
                        self[i] = min(bounds[1][i], self[i])
            # gp.into_bounds(self, copy_never=True)  # don't copy, because of side effect
        return self  # convenience return

#____________________________________________________________
#____________________________________________________________
#
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): 
        """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 = None  
        
        if bounds is not None:  # make the update stand-alone
            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 __call__(self, x, bounds=None):
        """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 bounds is None:
            bounds = self.bounds
        if bounds is None or bounds == [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:   
            if type(xi) is not Solution or xi.unrepaired is None:
                # 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
                xi = Solution(xi, copy=True).repair(bounds, keep_unrepaired=True) 
            fac = 1  # exp(0.1 * (log(self.scal) - np.mean(self.scal)))
            pen.append(sum(self.gamma * ((xi - xi.unrepaired) / 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)

        ### 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):
    """abstract base class for genotyp-phenotype transformation, 
    to be implemented. 
    
    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(GenoPhenoBase):
    """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 again 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. 

    """
    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):
        """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 bounds[0] is not None:
            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 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, copy=True, bounds=None):
        """maps the genotypic input argument into the phenotypic space, see
        help for class `GenoPheno`"""
        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):
        """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`` might still lead to a copy.
        """
        # TODO: to avoid computational overhead just do
        # if type(y) is Solution and y.geno: 
        #     return y.geno 
        
        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 OptimDataLogger 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 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 `OptimDataLogger` `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 =', Solution([ 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]
            
        """
        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 get 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 mirroring* is implicite 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)``
     
    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()
    <output omitted>
    >>>
    >>> assert es.result()[1] < 1e-9
    >>> assert es.result()[2] < 9000  # by internal termination
    >>> logger.plot()  # plot data
    >>> cma.show()
    >>> print('  *** if execution stalls close the figure window to continue (and check out ipython --pylab) ***')

    Example implementing restarts with increasing popsize (IPOP), output is not displayed:
    
    >>> import cma, numpy as np
    >>> 
    >>> # restart with increasing population size (IPOP)
    >>> bestever = cma.BestSolution()
    >>> for lam in 10 * 2**np.arange(7):  # 10, 20, 40, 80, ..., 10 * 2**6
    ...     es = cma.CMAEvolutionStrategy('6 - 8 * np.random.rand(9)',  # 9-D
    ...                                   5,         # initial std sigma0
    ...                                   {'popsize': lam, 
    ...                                    'verb_append': bestever.evalsall})   # pass options
    ...     logger = cma.CMADataLogger().register(es, append=bestever.evalsall) 
    ...     while not es.stop():
    ...         X = es.ask()    # get list of new solutions
    ...         fit = [cma.fcts.rastrigin(x) for x in X]  # evaluate each solution
    ...         es.tell(X, fit) # besides for termination only the ranking in fit is used
    ... 
    ...         # display some output
    ...         logger.add()  # add a "data point" to the log, writing in files
    ...         es.disp()  # uses option verb_disp with default 100
    ... 
    ...     print('termination:', es.stop())
    ...     cma.pprint(es.best.__dict__)
    ...     
    ...     bestever.update(es.best)
    ... 
    ...     # show a plot    
    ...     logger.plot();
    ...     if bestever.f < 1e-8:  # global optimum was hit
    ...         break
    <output omitted>
    >>> assert es.result()[1] < 1e-8
        
    On the Rastrigin function, usually after five restarts the global optimum 
    is located. 
        
    The final example shows how to resume:
        
    >>> import cma, pickle
    >>> 
    >>> es = cma.CMAEvolutionStrategy(12 * [0.1],  # a new instance, 12-D
    ...                               0.5)         # initial std sigma0
    >>> logger = cma.CMADataLogger().register(es)
    >>> es.optimize(cma.fcts.rosen, logger, iterations=100)
    >>> logger.plot() 
    >>> pickle.dump(es, open('saved-cma-object.pkl', 'w'))
    >>> print('saved')
    >>> del es, logger  # let's start fresh
    >>>
    >>> es = pickle.load(open('saved-cma-object.pkl'))
    >>> print('resumed')
    >>> logger = cma.CMADataLogger(es.opts['verb_filenameprefix']  # use same name
    ...                           ).register(es, True)  # True: append to old log data
    >>> es.optimize(cma.fcts.rosen, logger, verb_disp=200)
    >>> assert es.result()[2] < 15000
    >>> cma.pprint(es.result())
    >>> logger.plot()

    Missing Features
    ================
    Option ``randn`` to pass a random number generator. 
    
    :See: `fmin()`, `Options`, `plot()`, `ask()`, `tell()`, `ask_and_eval()`
            
    """

    # __all__ = ()  # TODO this would be the interface

    #____________________________________________________________
    @property  # read only attribute decorator for a method
    def popsize(self): 
        """number of samples by default returned by` ask()` 
        """
        return self.sp.popsize
        
    # this is not compatible with python2.5:
    #     @popsize.setter
    #     def popsize(self, p):
    #         """popsize cannot be set (this might change in future)
    #         """
    #         raise _Error("popsize cannot be changed (this might change in future)")        
        
    #____________________________________________________________
    #____________________________________________________________
    def stop(self, check=True):
        # this doc string is available via help cma.CMAEvolutionStrategy.stop
        """return a dictionary with the termination status. 
        With ``check==False``, the termination conditions are not checked and 
        the status might not reflect the current situation. 
        """
        return self.stopdict(self if check else None)  # update the stopdict and return a Dict
    
    #____________________________________________________________
    #____________________________________________________________
    def __init__(self, x0, sigma0, inopts = {}):
        """see class `CMAEvolutionStrategy` 

        """
        self.inputargs = dict(locals()) # for the record
        del self.inputargs['self'] # otherwise the instance self has a cyclic reference
        self.inopts = inopts
        opts = Options(inopts).complement()  # Options() == fmin([],[]) == defaultOptions()
        
        if x0 == str(x0):
            x0 = eval(x0)
        self.mean = array(x0)  # should not have column or row, is just 1-D
        if self.mean.ndim == 2:
            print('WARNING: input x0 should be a list or 1-D array, trying to flatten ' +
                    str(self.mean.shape) + '-array')
            if self.mean.shape[0] == 1:
                self.mean = self.mean[0]
            elif self.mean.shape[1] == 1:
                self.mean = array([x[0] for x in self.mean])
        self.x0 = self.mean.copy()  # for the record and output
        if self.mean.ndim != 1:
            raise _Error('x0 must be 1-D array')
        if len(self.mean) <= 0:
            raise _Error('optimization in 1-D is not supported (code was never tested)')
        
        self.N = self.mean.shape[0]
        N = self.N
        self.mean.resize(N) # 1-D array, not really necessary?!
        
        self.sigma0 = sigma0
        if isinstance(sigma0, str):  # TODO: no real need here (do rather in fmin)
            self.sigma0 = eval(sigma0)  # like '1./N' or 'np.random.rand(1)[0]+1e-2'
        if np.size(self.sigma0) != 1 or np.shape(self.sigma0): 
            raise _Error('input argument sigma0 must be (or evaluate to) a scalar')
        self.sigma = self.sigma0
        
        # extract/expand options 
        opts.evalall(locals())  # using only N
        self.opts = opts

        self.randn = opts['randn'] 
        self.gp = GenoPheno(N, opts['scaling_of_variables'], opts['typical_x'], 
            opts['bounds'], opts['fixed_variables'], opts['transformation']) 
        self.boundPenalty = BoundPenalty(self.gp.bounds)
        s = self.gp.geno(self.mean)
        self.mean = self.gp.geno(self.mean, bounds=self.gp.bounds)
        self.N = len(self.mean)  
        N = self.N 
        if (self.mean != s).any():
            print('WARNING: initial solution is out of the domain boundaries:')
            print('  x0   = ' + str(self.inputargs['x0']))
            print('  ldom = ' + str(self.gp.bounds[0]))
            print('  udom = ' + str(self.gp.bounds[1]))
        self.fmean = np.NaN             # TODO name should change? prints nan (OK with matlab&octave)
        self.fmean_noise_free = 0.  # for output only
        
        self.sp = CMAParameters(N, opts)
        self.sp0 = self.sp

        # initialization of state variables
        self.countiter = 0
        self.countevals = max((0, opts['verb_append'])) if type(opts['verb_append']) is not bool else 0
        self.ps = np.zeros(N)
        self.pc = np.zeros(N)

        stds = np.ones(N)
        if np.all(self.opts['CMA_teststds']):  # also 0 would not make sense
            stds = self.opts['CMA_teststds']
            if np.size(stds) != N:
                raise _Error('CMA_teststds option must have dimension = ' + str(N)) 
        if self.opts['CMA_diagonal']:  # is True or > 0
            # linear time and space complexity 
            self.B = array(1) # works fine with np.dot(self.B, anything) and self.B.T
            self.C = stds**2
            self.dC = self.C
        else:
            self.B = np.eye(N) # identity(N), do not from matlib import *, as eye is a matrix there
            # prevent equal eigenvals, a hack for np.linalg:
            self.C = np.diag(stds**2 * exp(1e-6*np.random.randn(N)))
            self.dC = np.diag(self.C)
            self.Zneg = np.zeros((N, N))
        self.D = stds

        self.flgtelldone = True
        self.itereigenupdated = self.countiter
        self.noiseS = 0  # noise "signal"
        self.hsiglist = [] 

        if not opts['seed']:
            np.random.seed() 
            six_decimals = (time.time() - 1e6 * (time.time() // 1e6))
            opts['seed'] = 1e5 * np.random.rand() + six_decimals + 1e5 * (time.time() % 1) 
        opts['seed'] = int(opts['seed'])
        np.random.seed(opts['seed'])
        
        self.best = BestSolution()
        
        out = {}
        out['best'] = self.best
        # out['hsigcount'] = 0
        out['termination'] = {}
        self.out = out
        
        self.const = BlancClass()
        self.const.chiN = N**0.5*(1-1./(4.*N)+1./(21.*N**2)) # expectation of norm(randn(N,1))

        # attribute for stopping criteria in function stop 
        self.stopdict = CMAStopDict()
        self.callbackstop = 0 

        self.fit = BlancClass() 
        self.fit.fit = []   # not really necessary 
        self.fit.hist = []  # short history of best
        self.fit.histbest = []   # long history of best
        self.fit.histmedian = [] # long history of median

        self.more_to_write = []  #[1, 1, 1, 1]  #  N*[1]  # needed when writing takes place before setting

        # say hello
        if opts['verb_disp'] > 0:
            sweighted = '_w' if self.sp.mu > 1 else ''
            smirr = 'mirr%d' % (self.sp.lam_mirr) if self.sp.lam_mirr else ''
            print('(%d' % (self.sp.mu) + sweighted + ',%d' % (self.sp.popsize) + smirr + ')-CMA-ES' +
                  ' (mu_w=%2.1f,w_1=%d%%)' % (self.sp.mueff, int(100*self.sp.weights[0])) +
                  ' in dimension %d (seed=%d, %s)' % (N, opts['seed'], time.asctime())) # + func.__name__ 
            if opts['CMA_diagonal'] and self.sp.CMA_on:
                s = ''
                if opts['CMA_diagonal'] is not True:
                    s = ' for '
                    if opts['CMA_diagonal'] < np.inf:
                        s += str(int(opts['CMA_diagonal']))
                    else:
                        s += str(np.floor(opts['CMA_diagonal']))
                    s += ' iterations'
                    s += ' (1/ccov=' + str(round(1./(self.sp.c1+self.sp.cmu))) + ')'
                print('   Covariance matrix is diagonal' + s) 

    #____________________________________________________________
    #____________________________________________________________
    def ask(self, number=None, xmean=None, sigma_fac=1):
        """get new candidate solutions, sampled from a multi-variate
        normal distribution and transformed to f-representation
        (phenotype) to be evaluated. 
           
        Arguments
        ---------
            `number` 
                number of returned solutions, by default the
                population size ``popsize`` (AKA ``lambda``).
            `xmean` 
                distribution mean 
            `sigma` 
                multiplier for internal sample width (standard
                deviation)
        
        Return
        ------
        A list of N-dimensional candidate solutions to be evaluated
        
        Example
        -------
        >>> import cma
        >>> es = cma.CMAEvolutionStrategy([0,0,0,0], 0.3)
        >>> while not es.stop() and es.best.f > 1e-6:  # my_desired_target_f_value
        ...     X = es.ask()  # get list of new solutions
        ...     fit = [cma.fcts.rosen(x) for x in X]  # call function rosen with each solution
        ...     es.tell(X, fit)  # feed values
        
        :See: `ask_and_eval`, `ask_geno`, `tell`
                
        """
        pop = self.ask_geno(number, xmean, sigma_fac)
        
        # N,lambda=20,200: overall CPU 7s vs 5s == 40% overhead, even without bounds!
        # TODO: check here, whether this is necessary? 
        # return [self.gp.pheno(x, copy=False, bounds=self.gp.bounds) for x in pop]  # probably fine
        # return [Solution(self.gp.pheno(x, copy=False), copy=False) for x in pop]  # here comes the memory leak, now solved  
        return [Solution(self.gp.pheno(x, copy=False), copy=False).repair(self.gp.bounds) for x in pop]  
        
        # TODO: check out this:
        return [Solution(x, gp=self.gp).repair(self.gp.bounds) for x in pop]  
        # TODO: check out this (old):
        # return [Solution(x, gp=self.gp, bounds=self.gp.bounds) for x in pop]
        # and use solution.geno in tell()

    #____________________________________________________________
    #____________________________________________________________
    def ask_geno(self, number=None, xmean=None, sigma_fac=1):
        """get new candidate solutions in genotyp, sampled from a 
        multi-variate normal distribution. 
           
        Arguments are
            `number`
                number of returned solutions, by default the
                population size `popsize` (AKA lambda).
            `xmean` 
                distribution mean 
            `sigma_fac`
                multiplier for internal sample width (standard
                deviation)
        
        `ask_geno` returns a list of N-dimensional candidate solutions 
        in genotyp representation and is called by `ask`. 
        
        :See: `ask`, `ask_and_eval`
                        
        """
        
        if number is None or number < 1:
            number = self.sp.popsize
        if xmean is None:
            xmean = self.mean  

        if self.countiter == 0:
            self.tic = time.clock()  # backward compatible
            self.elapsed_time = ElapsedTime()

        sigma = sigma_fac * self.sigma

        # update parameters for sampling the distribution
        #        fac  0      1      10
        # 150-D cigar: 
        #           50749  50464   50787
        # 200-D elli:               == 6.9
        #                  99900   101160
        #                 100995   103275 == 2% loss
        # 100-D elli:               == 6.9
        #                 363052   369325  < 2% loss
        #                 365075   365755

        # update distribution
        if self.sp.CMA_on and (
                (self.opts['updatecovwait'] is None and 
                 self.countiter >=
                     self.itereigenupdated + 1./(self.sp.c1+self.sp.cmu)/self.N/10
                 ) or
                (self.opts['updatecovwait'] is not None and
                 self.countiter > self.itereigenupdated + self.opts['updatecovwait']
                 )):
            self.updateBD()

        # sample distribution
        if self.flgtelldone:  # could be done in tell()!? 
            self.flgtelldone = False
            self.ary = [] 

        # each row is a solution
        self.ary = np.dot(self.B, (self.D * self.randn((number, self.N))).T).T
        pop = xmean + sigma * self.ary
        self.evaluations_per_f_value = 1

        return pop

    def get_mirror(self, x):
        """return ``pheno(self.mean - (geno(x) - self.mean))``.
        
        TODO: this implementation is yet experimental.
        
        Selectively mirrored sampling improves to a moderate extend but 
        overadditively with active CMA for quite understandable reasons. 
        
        Optimal number of mirrors are suprisingly small: 1,2,3 for maxlam=7,13,20
        however note that 3,6,10 are the respective maximal possible mirrors that
        must be clearly suboptimal.         

        """
        try:
            dx = x.geno - self.mean  
        except:
            dx = self.gp.geno(x, copy=True) - self.mean
            
        dx *= self.sigma * (sum(np.random.randn(self.N)**2) / sum((np.dot(self.B.T, dx) / self.D)**2))**0.5
        return Solution(self.mean - dx, gp=self.gp).repair(self.gp.bounds)
        
    def mirror_penalized(self, f_values, idx):
        """obsolete and subject to removal (TODO), 
        return modified f-values such that for each mirror one becomes worst. 
        
        This function is useless when selective mirroring is applied with no
        more than (lambda-mu)/2 solutions. 
        
        Mirrors are leading and trailing values in ``f_values``.
        
        """
        assert len(f_values) >= 2 * len(idx)
        m = np.max(np.abs(f_values))
        for i in len(idx):
            if f_values[idx[i]] > f_values[-1-i]:
                f_values[idx[i]] += m
            else:
                f_values[-1-i] += m
        return f_values
    
    def mirror_idx_cov(self, f_values, idx1):  # will most likely be removed
        """obsolete and subject to removal (TODO), 
        return indices for negative ("active") update of the covariance matrix
        assuming that ``f_values[idx1[i]]`` and ``f_values[-1-i]`` are
        the corresponding mirrored values 
        
        computes the index of the worse solution sorted by the f-value of the
        better solution. 
        
        TODO: when the actual mirror was rejected, it is better
        to return idx1 instead of idx2. 
        
        Remark: this function might not be necessary at all: if the worst solution
        is the best mirrored, the covariance matrix updates cancel (cave: weights
        and learning rates), which seems what is desirable. If the mirror is bad, 
        as strong negative update is made, again what is desirable. 
        And the fitness--step-length correlation is in part addressed by 
        using flat weights. 
         
        """
        idx2 = np.arange(len(f_values) - 1, len(f_values) - 1 - len(idx1), -1)
        f = []
        for i in xrange(len(idx1)):
            f.append(min((f_values[idx1[i]], f_values[idx2[i]])))
            # idx.append(idx1[i] if f_values[idx1[i]] > f_values[idx2[i]] else idx2[i])
        return idx2[np.argsort(f)][-1::-1]
    
    #____________________________________________________________
    #____________________________________________________________
    # 
    def ask_and_eval(self, func, args=(), number=None, xmean=None, sigma_fac=1, 
                     evaluations=1, aggregation=np.median):
        """samples `number` solutions and evaluates them on `func`, where
        each solution `s` is resampled until ``func(s) not in (numpy.NaN, None)``.  
            
        Arguments
        ---------
            `func` 
                objective function
            `args` 
                additional parameters for `func`
            `number` 
                number of solutions to be sampled, by default
                population size ``popsize`` (AKA lambda)
            `xmean` 
                mean for sampling the solutions, by default ``self.mean``. 
            `sigma_fac` 
                multiplier for sampling width, standard deviation, for example
                to get a small perturbation of solution `xmean`
            `evaluations`
                number of evaluations for each sampled solution
            `aggregation`
                function that aggregates `evaluations` values to
                as single value. 
        
        Return
        ------
        ``(X, fit)``, where
            X -- list of solutions
            fit -- list of respective function values
            
        Details
        -------
        When ``func(x)`` returns `NaN` or `None` a new solution is sampled until
        ``func(x) not in (numpy.NaN, None)``.  The argument to `func` can be 
        freely modified within `func`.
        
        Depending on the ``CMA_mirrors`` option, some solutions are not sampled
        independently but as mirrors of other bad solutions. This is a simple 
        derandomization that can save 10-30% of the evaluations in particular 
        with small populations, for example on the cigar function. 
        
        Example
        -------
        >>> import cma
        >>> x0, sigma0 = 8*[10], 1  # 8-D
        >>> es = cma.CMAEvolutionStrategy(x0, sigma0) 
        >>> while not es.stop():
        ...     X, fit = es.ask_and_eval(cma.fcts.elli)  # handles NaN with resampling
        ...     es.tell(X, fit)  # pass on fitness values
        ...     es.disp(20) # print every 20-th iteration
        >>> print('terminated on ' + str(es.stop()))
        <output omitted>

        A single iteration step can be expressed in one line, such that
        an entire optimization after initialization becomes
        ::
        
            while not es.stop():
                es.tell(*es.ask_and_eval(cma.fcts.elli)) 
        
        """
        # initialize
        popsize = self.sp.popsize  
        if number is not None:
            popsize = number
        selective_mirroring = True
        nmirrors = self.sp.lam_mirr 
        if popsize != self.sp.popsize:
            nmirrors = Mh.sround(popsize * self.sp.lam_mirr / self.sp.popsize)
            # TODO: now selective mirroring might be impaired 
        assert nmirrors <= popsize // 2
        self.mirrors_idx = np.arange(nmirrors)  # might never be used
        self.mirrors_rejected_idx = []  # might never be used
        if xmean is None:
            xmean = self.mean
        
        # do the work
        fit = []  # or np.NaN * np.empty(number)
        X = []
        for k in xrange(int(popsize)):
            nreject = -1
            f = np.NaN
            while f in (np.NaN, None):  # rejection sampling
                nreject += 1
                if k < popsize - nmirrors or nreject:  
                    x = self.ask(1, xmean, sigma_fac)[0] 
                else:  # mirrored sample
                    if k == popsize - nmirrors and selective_mirroring:
                        self.mirrors_idx = np.argsort(fit)[-1:-1-nmirrors:-1]
                    x = self.get_mirror(X[self.mirrors_idx[popsize - 1 - k]])
                if nreject == 1 and k >= popsize - nmirrors:
                    self.mirrors_rejected_idx.append(k)
                f = func(x, *args)
                if f not in (np.NaN, None) and evaluations > 1:
                    f = aggregation([f] + [func(x, *args) for _i in xrange(int(evaluations-1))])
                if nreject + 1 % 1000 == 0:
                    print('  %d solutions rejected (f-value NaN or None) at iteration' % 
                          (nreject, self.countiter))
            fit.append(f)
            X.append(x)
        self.evaluations_per_f_value = int(evaluations)
        return X, fit


    #____________________________________________________________
    def tell(self, solutions, function_values, 
                function_values_reevaluated=None, check_points=None, copy=False):  
        """pass objective function values to prepare for next
        iteration. This core procedure of the CMA-ES algorithm updates 
        all state variables: two evolution paths, the distribution mean, 
        the covariance matrix and a step-size. 

        Arguments
        ---------
            `solutions` 
                list or array of candidate solution points (of
                type `cma.Solution` or `numpy.ndarray`), most presumably before 
                delivered by method `ask()` or `ask_and_eval()`.
            `function_values`
                list or array of objective function values 
                corresponding to the respective points. Beside for termination 
                decisions, only the ranking of values in `function_values`
                is used.
            `check_points`
                if ``True``, allows to savely pass solutions that are 
                not necessarily generated using `ask()`. Might just as well be a 
                list of indices to be checked in solutions. Value ``None`` defaults
                to ``False``. 
            `copy`
                ``solutions`` might be modified, if ``copy is False``
        
        Details
        -------
        `tell()` updates the parameters of the multivariate
        normal search distribtion, namely covariance matrix and
        step-size and updates also the attributes `countiter` and 
        `countevals`. To check the points for consistency is quadratic 
        in the dimension (like sampling points). 
        
        Bugs
        ----
        The effect of changing the solutions delivered by `ask()` depends on whether 
        boundary handling is applied. With boundary handling, modifications are 
        disregarded. This is necessary to apply the default boundary handling that
        uses unrepaired solutions but might change in future.   

        Example
        -------
        ::
        
            import cma
            func = cma.fcts.elli  # choose objective function
            es = cma.CMAEvolutionStrategy(cma.np.random.rand(10), 1)
            while not es.stop():
               X = es.ask()
               es.tell(X, [func(x) for x in X])
            es.result()  # where the result can be found
            
        :See: class `CMAEvolutionStrategy`, `ask()`, `ask_and_eval()`, `fmin()`

        """
    #____________________________________________________________
    # TODO: consider an input argument that flags injected trust-worthy solutions (which means 
    #       that they can be treated "absolut" rather than "relative")
        if self.flgtelldone: 
            raise _Error('tell should only be called once per iteration')
        
        if check_points is None:
            check_points = self.opts['check_points']
            if check_points is None:
                check_points = False

        lam = len(solutions)
        if lam != array(function_values).shape[0]:
            raise _Error('for each candidate solution '
                        + 'a function value must be provided')
        if lam + self.sp.lam_mirr < 3: 
            raise _Error('population size ' + str(lam) + ' is too small when option CMA_mirrors * popsize < 0.5')
#        print function_values
#        if not np.isscalar(function_values[0]):
#            if np.isscalar(function_values[0][0]):
#                if self.countiter <= 1:
#                    print('WARNING: function values are not a list of scalars (further warnings are suppressed)')
#                function_values = [val[0] for val in function_values]
#            else:
#                raise _Error('objective function values must be a list of scalars')
#        
#            
        ### prepare
        N = self.N
        sp = self.sp
        if 11 < 3 and lam != sp.popsize:  # turned off, because mu should stay constant, still not desastrous  
            print('WARNING: population size has changed, recomputing parameters')
            self.sp.set(self.opts, lam)  # not really tested
        if lam < sp.mu:  # rather decrease cmean instead of having mu > lambda//2
            raise _Error('not enough solutions passed to function tell (mu>lambda)')
            
        self.countiter += 1  # >= 1 now
        self.countevals += sp.popsize * self.evaluations_per_f_value
        self.best.update(solutions, function_values, self.countevals)
        
        flgseparable = self.opts['CMA_diagonal'] is True \
                       or self.countiter <= self.opts['CMA_diagonal'] 
        if not flgseparable and len(self.C.shape) == 1:  # C was diagonal ie 1-D
            # enter non-separable phase (no easy return from here)
            self.B = np.eye(N) # identity(N)
            self.C = np.diag(self.C)
            idx = np.argsort(self.D)
            self.D = self.D[idx]
            self.B = self.B[:,idx]
            self.Zneg = np.zeros((N, N))

        ### manage fitness
        fit = self.fit  # make short cut
        
        # CPU for N,lam=20,200: this takes 10s vs 7s 
        fit.bndpen = self.boundPenalty.update(function_values, self)(solutions)
        # for testing: 
        # fit.bndpen = self.boundPenalty.update(function_values, self)([s.unrepaired for s in solutions])
        fit.idx = np.argsort(array(fit.bndpen) + array(function_values))
        fit.fit = array(function_values, copy=False)[fit.idx]

        # update output data TODO: this is obsolete!? However: need communicate current best x-value? 
        # old: out['recent_x'] = self.gp.pheno(pop[0])
        self.out['recent_x'] = array(solutions[fit.idx[0]])  # TODO: change in a data structure(?) and use current as identify
        self.out['recent_f'] = fit.fit[0]

        # fitness histories
        fit.hist.insert(0, fit.fit[0])
        # if len(self.fit.histbest) < 120+30*N/sp.popsize or  # does not help, as tablet in the beginning is the critical counter-case
        if ((self.countiter % 5) == 0):  # 20 percent of 1e5 gen.
            fit.histbest.insert(0, fit.fit[0])
            fit.histmedian.insert(0, np.median(fit.fit) if len(fit.fit<21) 
                                    else fit.fit[self.popsize // 2])
        if len(fit.histbest) > 2e4: # 10 + 30*N/sp.popsize:
            fit.histbest.pop()
            fit.histmedian.pop()
        if len(fit.hist) > 10 + 30*N/sp.popsize:
            fit.hist.pop()

        # TODO: clean up inconsistency when an unrepaired solution is available and used
        pop = []  # create pop from input argument solutions
        for s in solutions:  # use phenotype before Solution.repair() 
            if type(s) is Solution and s.unrepaired is not None:
                pop.append(s.unrepaired)  # TODO: this might not be desirable, if solutions were modified/injected
            else:
                pop.append(s)
            
            if type(pop[-1]) is Solution and pop[-1].isgeno:  # not quite likely
                pass  # need to do nothing
            elif 11 < 3 and type(pop[-1]) is Solution and s.geno is not None:
                pop[-1] = s.geno  # TODO: this might not be desirable, if solutions were modified/injected 
            else:  # most likely case
                # TODO: check whether self.gp.geno is efficient 
                pop[-1] = self.gp.geno(pop[-1], copy=copy) 

        mold = self.mean
        sigma_fac = 1

        # check and normalize each x - m
        # check_points is a flag or an index list
        # should also a number possible (first check_points points)? 
        if check_points not in (None, False, 0, [], ()):  # useful in case of injected solutions and/or adaptive encoding
            try: 
                if len(check_points):
                    idx = check_points
            except:
                idx = xrange(sp.popsize)  
            
            self.clip_or_fit_solutions(pop, idx)
                
        # sort pop  
        if type(pop) is not array: # only arrays can be multiple indexed
            pop = array(pop, copy=False)

        pop = pop[fit.idx] 

        if self.opts['CMA_elitist'] and self.best.f < fit.fit[0]:
            xp = [self.gp.geno(self.best.x[:])]
            # print self.mahalanobisNorm(xp[0]-self.mean)
            self.clip_or_fit_solutions(xp, [0])
            pop = array([xp[0]] + list(pop))

        # compute new mean 
        self.mean = mold + self.sp.cmean * \
                    (sum(sp.weights * pop[0:sp.mu].T, 1) - mold)
        
        
        # check Delta m (this is not default, but could become at some point)
        # CAVE: upper_length=sqrt(2)+2 is too restrictive, test upper_length = sqrt(2*N) thoroughly.
        # simple test case injecting self.mean: 
        # self.mean = 1e-4 * self.sigma * np.random.randn(N)
        if 11 < 3 and self.opts['vv'] and check_points:  # TODO: check_points might be an index-list
            cmean = self.sp.cmean / min(1, (sqrt(self.opts['vv']*N)+2) / ( # abuse of cmean
                (sqrt(self.sp.mueff) / self.sp.cmean) * 
                self.mahalanobisNorm(self.mean - mold))) 
        else:
            cmean = self.sp.cmean

        if 11 < 3:  # plot length of mean - mold
            self.more_to_write = [sqrt(sp.mueff) * 
                sum(((1./self.D) * dot(self.B.T, self.mean - mold))**2)**0.5 / 
                       self.sigma / sqrt(N) / cmean]

        # get learning rate constants
        cc, c1, cmu = sp.cc, sp.c1, sp.cmu 
        if flgseparable:
            cc, c1, cmu = sp.cc_sep, sp.c1_sep, sp.cmu_sep
        
        # now the real work can start 
        
        # evolution paths
        self.ps = (1-sp.cs) * self.ps + \
                  (sqrt(sp.cs*(2-sp.cs)*sp.mueff)  / self.sigma / cmean) * \
                  dot(self.B, (1./self.D) * dot(self.B.T, self.mean - mold))
                  
        # "hsig"
        hsig = sum(self.ps**2) / (1-(1-sp.cs)**(2*self.countiter)) / self.N < 2 + 4./(N+1)
        # adjust missing variance due to hsig, in 4-D with damps=1e99 and sig0 small  
        #       hsig leads to premature convergence of C otherwise
        #hsiga = (1-hsig**2) * c1 * cc * (2-cc)  # to be removed in future
        c1a = c1 - (1-hsig**2) * c1 * cc * (2-cc)  # adjust for variance loss
        
        if 11 < 3:  # diagnostic data
            self.out['hsigcount'] += 1 - hsig
            if not hsig:
                self.hsiglist.append(self.countiter)
        if 11 < 3:  # diagnostic message
            if not hsig:  
                print(str(self.countiter) + ': hsig-stall')
        if 11 < 3:  # for testing purpose
            hsig = 1 # TODO: 
            #       put correction term, but how? 
            if self.countiter == 1:
                print('hsig=1')
        
        self.pc = (1-cc) * self.pc + \
                  hsig * (sqrt(cc*(2-cc)*sp.mueff) / self.sigma / cmean) * \
                  (self.mean - mold) 

        # covariance matrix adaptation/udpate 
        if sp.CMA_on:
            # assert sp.c1 + sp.cmu < sp.mueff / N  # ??
            assert c1 + cmu <= 1
            
            # default full matrix case 
            if not flgseparable:
                Z = (pop[0:sp.mu] - mold) / self.sigma
                Z = dot((cmu * sp.weights) * Z.T, Z)  # learning rate integrated
                if self.sp.neg.cmuexp:
                    tmp = (pop[-sp.neg.mu:] - mold) / self.sigma
                    self.Zneg *= 1 - self.sp.neg.cmuexp  # not necessary?
                    self.Zneg += dot(sp.neg.weights * tmp.T, tmp) - self.C
                    # self.update_exponential(dot(sp.neg.weights * tmp.T, tmp) - 1 * self.C, -1*self.sp.neg.cmuexp)
                        
                if 11 < 3: # ?3 to 5 times slower??
                    Z = np.zeros((N,N))
                    for k in xrange(sp.mu): 
                        z = (pop[k]-mold) 
                        Z += np.outer((cmu * sp.weights[k] / self.sigma**2) * z, z)

                self.C *= 1 - c1a - cmu
                self.C += np.outer(c1 * self.pc, self.pc) + Z
                self.dC = np.diag(self.C)  # for output and termination checking 

            else: # separable/diagonal linear case 
                assert(c1+cmu <= 1)
                Z = np.zeros(N)
                for k in xrange(sp.mu):
                    z = (pop[k]-mold) / self.sigma  # TODO see above
                    Z += sp.weights[k] * z * z  # is 1-D
                self.C = (1-c1a-cmu) * self.C + c1 * self.pc * self.pc + cmu * Z
                # TODO: self.C *= exp(cmuneg * (N - dot(sp.neg.weights,  **2)
                self.dC = self.C
                self.D = sqrt(self.C)  # C is a 1-D array
                self.itereigenupdated = self.countiter
                
                # idx = self.mirror_idx_cov()  # take half of mirrored vectors for negative update
                
        # step-size adaptation, adapt sigma
        self.sigma *= sigma_fac * \
                        np.exp(min(1, (sp.cs/sp.damps) *
                                (sqrt(sum(self.ps**2))/self.const.chiN - 1)))

        if 11 < 3 and self.opts['vv']:
            if self.countiter < 2:
                print('constant sigma applied')
                print(self.opts['vv'])  # N=10,lam=10: 0.8 is optimal
            self.sigma = self.opts['vv'] * self.sp.mueff * sum(self.mean**2)**0.5 / N

        if self.sigma * min(self.dC)**0.5 < self.opts['minstd']:
            self.sigma = self.opts['minstd'] / min(self.dC)**0.5
        if self.sigma * min(self.D) < self.opts['mindx']:
            self.sigma = self.opts['mindx'] / min(self.D)

        # TODO increase sigma in case of a plateau? 

        # Uncertainty noise measurement tobe-inserted

        # output, has moved up, e.g. as part of fmin, TODO to be removed
        if 11 < 3 and self.opts['verb_log'] > 0 and (self.countiter < 4 or
                                          self.countiter % self.opts['verb_log'] == 0):
            # this assumes that two logger with the same name access the same data!
            CMADataLogger(self.opts['verb_filenameprefix']).register(self, append=True).add()
            # self.writeOutput(solutions[fit.idx[0]])
            
        self.flgtelldone = True
    # end tell()
        
    def result(self):
        """return ``(xbest, f(xbest), evaluations_xbest, evaluations, iterations, pheno(xmean), effective_stds)``"""
        # TODO: how about xcurrent? 
        return self.best.get() + (
            self.countevals, self.countiter, self.gp.pheno(self.mean), self.gp.scales * self.sigma * self.dC**0.5)

    def clip_or_fit_solutions(self, pop, idx):
        """make sure that solutions fit to sample distribution, this interface will probably change. 
        
        In particular the frequency of long vectors appearing in pop[idx] - self.mean is limited. 
        
        """
        mold = self.mean
        if 1 < 3:  # hard clip at upper_length
            upper_length = self.N**0.5 + 2 * self.N / (self.N+2)  # should become an Option 
            for k in idx: 
                fac = self.mahalanobisNorm(pop[k] - mold) / upper_length
                
                if fac > 1:
                    pop[k] = (pop[k] - mold) / fac + mold
                    # print self.countiter, k, fac, self.mahalanobisNorm(pop[k] - mold)
                    # adapt also sigma: which are the trust-worthy/injected solutions? 
                elif 11 < 3 and k == fit.idx[0]: 
                    return exp(np.tanh(((upper_length*fac)**2/N-1)/2) / 2)
        else:
            if 'checktail' not in self.__dict__:  # hasattr(self, 'checktail')  
                from check_tail_smooth import CheckTail  # for the time being
                self.checktail = CheckTail()
                print('untested feature checktail is on')
            for k in idx: 
                fac = self.checktail.addchin(self.mahalanobisNorm(pop[k] - mold))
                
                if fac < 1:
                    pop[k] = fac * (pop[k] - mold) + mold
        
        return 1.0  # sigma_fac, not in use
                
            
    #____________________________________________________________
    #____________________________________________________________
    # 
    def updateBD(self): 
        """update internal variables for sampling the distribution with the
        current covariance matrix C. This method is O(N^3), if C is not diagonal.
        
        """
        # itereigenupdated is always up-to-date in the diagonal case 
        # just double check here
        if self.itereigenupdated == self.countiter:
            return

        if self.sp.neg.cmuexp:  # cave: 
            self.update_exponential(self.Zneg, -self.sp.neg.cmuexp)
            # self.C += self.Zpos  # pos update after Zneg would be the correct update, overall: 
            # self.C = self.Zpos + Cs * Mh.expms(-self.sp.neg.cmuexp*Csi*self.Zneg*Csi) * Cs
            self.Zneg = np.zeros((self.N, self.N))

        self.C = (self.C + self.C.T) / 2
        # self.C = np.triu(self.C) + np.triu(self.C,1).T  # should work as well
        # self.D, self.B = eigh(self.C) # hermitian, ie symmetric C is assumed
        if self.opts['CMA_eigenmethod'] == -1:  
            # pygsl
            # easy to install (well, in Windows install gsl binaries first, 
            # set system path to respective libgsl-0.dll (or cp the dll to
            # python\DLLS ?), in unzipped pygsl edit
            # gsl_dist/gsl_site_example.py into gsl_dist/gsl_site.py
            # and run "python setup.py build" and "python setup.py install"
            # in MINGW32)
            if 1 < 3:  # import pygsl on the fly
                try:
                    import pygsl.eigen.eigenvectors  # TODO efficient enough?
                except ImportError:
                    print('WARNING: could not find pygsl.eigen module, either install pygsl \n' +
                          '  or set option CMA_eigenmethod=1 (is much slower), option set to 1')
                    self.opts['CMA_eigenmethod'] = 0  # use 0 if 1 is too slow
                    
                self.D, self.B = pygsl.eigen.eigenvectors(self.C)

        elif self.opts['CMA_eigenmethod'] == 0: 
            # TODO: thoroughly test np.linalg.eigh 
            #       numpy.linalg.eig crashes in 200-D
            #       and EVecs with same EVals are not orthogonal 
            self.D, self.B = np.linalg.eigh(self.C)  # self.B[i] is a row and not an eigenvector
        else:  # is overall two;ten times slower in 10;20-D
            self.D, self.B = Misc.eig(self.C)  # def eig, see below
        
        # assert(sum(self.D-DD) < 1e-6)
        # assert(sum(sum(np.dot(BB, BB.T)-np.eye(self.N))) < 1e-6)
        # assert(sum(sum(np.dot(BB * DD, BB.T) - self.C)) < 1e-6)
        idx = np.argsort(self.D)
        self.D = self.D[idx]
        self.B = self.B[:,idx]  # self.B[i] is a row, columns self.B[:,i] are eigenvectors
        # assert(all(self.B[self.countiter % self.N] == self.B[self.countiter % self.N,:]))
            
        if 11 < 3 and any(abs(sum(self.B[:,0:self.N-1] * self.B[:,1:], 0)) > 1e-6): 
            print('B is not orthogonal')
            print(self.D)
            print(sum(self.B[:,0:self.N-1] * self.B[:,1:], 0)) 
        else:
            # is O(N^3)
            # assert(sum(abs(self.C - np.dot(self.D * self.B,  self.B.T))) < N**2*1e-11) 
            pass
        self.D **= 0.5
        self.itereigenupdated = self.countiter

    def update_exponential(self, Z, eta, BDpair=None):
        """exponential update of C that guarantees positive definiteness, that is, 
        instead of the assignment ``C = C + eta * Z``,  
        C gets C**.5 * exp(eta * C**-.5 * Z * C**-.5) * C**.5. 
        
        Parameter Z should have expectation zero, e.g. sum(w[i] * z[i] * z[i].T) - C
        if E z z.T = C. 
        
        This function conducts two eigendecompositions, assuming that
        B and D are not up to date, unless `BDpair` is given. Given BDpair, 
        B is the eigensystem and D is the vector of sqrt(eigenvalues), one
        eigendecomposition is omitted. 
        
        Reference: Glasmachers et al 2010, Exponential Natural Evolution Strategies
        
        """
        if eta == 0:
            return
        if BDpair:
            B, D = BDpair
        else:
            D, B = np.linalg.eigh(self.C)
            D **= 0.5
        Csi = dot(B, (B / D).T)
        Cs = dot(B, (B * D).T)
        self.C = dot(Cs, dot(Mh.expms(eta * dot(Csi, dot(Z, Csi))), Cs))

    #____________________________________________________________
    #____________________________________________________________
    # 
    def _updateCholesky(self, A, Ainv, p, alpha, beta): 
        """not yet implemented"""
        # BD is A, p is A*Normal(0,I) distributed
        # input is assumed to be numpy arrays
        # Ainv is needed to compute the evolution path
        # this is a stump and is not tested

        raise _Error("not yet implemented")
        # prepare
        alpha = float(alpha)
        beta = float(beta)
        y = np.dot(Ainv, p)
        y_sum = sum(y**2)

        # compute scalars
        tmp = sqrt(1 + beta * y_sum / alpha)
        fac = (sqrt(alpha) / sum(y**2)) * (tmp - 1)
        facinv = (1. / (sqrt(alpha) * sum(y**2))) * (1 - 1. / tmp)

        # update matrices
        A *= sqrt(alpha)
        A += np.outer(fac * p, y)
        Ainv /= sqrt(alpha)
        Ainv -= np.outer(facinv * y, np.dot(y.T, Ainv))

    #____________________________________________________________
    #____________________________________________________________
    def feedForResume(self, X, function_values):
        """Given all "previous" candidate solutions and their respective 
        function values, the state of a `CMAEvolutionStrategy` object 
        can be reconstructed from this history. This is the purpose of 
        function `feedForResume`.  

        Arguments
        ---------
            X
              (all) solution points in chronological order, phenotypic 
              representation. The number of points must be a multiple
              of popsize. 
            function_values 
              respective objective function values

        Details
        -------
        `feedForResume` can be called repeatedly with only parts of 
        the history. The part must have the length of a multiple
        of the population size. 
        `feedForResume` feeds the history in popsize-chunks into `tell`.   
        The state of the random number generator might not be 
        reconstructed, but this would be only relevant for the future.        
        
        Example
        -------
        ::
        
            import cma
            
            # prepare 
            (x0, sigma0) = ... # initial values from previous trial
            X = ... # list of generated solutions from a previous trial
            f = ... # respective list of f-values
            
            # resume
            es = cma.CMAEvolutionStrategy(x0, sigma0)  
            es.feedForResume(X, f)
            
            # continue with func as objective function
            while not es.stop():
               X = es.ask()
               es.tell(X, [func(x) for x in X])
        
        Credits to Dirk Bueche and Fabrice Marchal for the feeding idea.     
        
        :See: class `CMAEvolutionStrategy` for a simple dump/load to resume

        """
        if self.countiter > 0:
            print('WARNING: feed should generally be used with a new object instance')
        if len(X) != len(function_values):
            raise _Error('number of solutions ' + str(len(X)) +
                ' and number function values ' +
                str(len(function_values))+' must not differ')
        popsize = self.sp.popsize
        if (len(X) % popsize) != 0:
            raise _Error('number of solutions ' + str(len(X)) +
                    ' must be a multiple of popsize (lambda) ' + 
                    str(popsize))
        for i in xrange(len(X) / popsize):
            # feed in chunks of size popsize
            self.ask()  # a fake ask, mainly for a conditioned calling of updateBD
                        # and secondary to get possibly the same random state
            self.tell(X[i*popsize:(i+1)*popsize], function_values[i*popsize:(i+1)*popsize])
            
    #____________________________________________________________
    #____________________________________________________________
    def readProperties(self):
        """reads dynamic parameters from property file (not implemented)
        """
        print('not yet implemented')

    #____________________________________________________________
    #____________________________________________________________
    def mahalanobisNorm(self, dx):
        """
        compute the Mahalanobis norm that is induced by the adapted covariance
        matrix C times sigma**2.
         
        Argument
        --------
        A *genotype* difference `dx`. 
        
        Example
        -------
        >>> import cma, numpy
        >>> es = cma.CMAEvolutionStrategy(numpy.ones(10), 1)
        >>> xx = numpy.random.randn(2, 10)
        >>> d = es.mahalanobisNorm(es.gp.geno(xx[0]-xx[1]))  
        
        `d` is the distance "in" the true sample distribution,
        sampled points have a typical distance of ``sqrt(2*es.N)``, 
        where `N` is the dimension. In the example, `d` is the 
        Euclidean distance, because C = I and sigma = 1. 
        
        """
        return sqrt(sum((self.D**-1 * np.dot(self.B.T, dx))**2)) / self.sigma
    
    #____________________________________________________________
    #____________________________________________________________
    def disp_annotation(self):
        """print annotation for `disp()`"""
        print('Iterat #Fevals   function value     axis ratio  sigma   minstd maxstd min:sec')
        sys.stdout.flush()

    #____________________________________________________________
    #____________________________________________________________
    def disp(self, modulo=None):  # TODO: rather assign opt['verb_disp'] as default? 
        """prints some infos according to `disp_annotation()`, if
        ``iteration_counter % modulo == 0``
        
        """
        if modulo is None:
            modulo = self.opts['verb_disp'] 
            
        # console display 
        if modulo:
            if 0.1*(self.countiter-1) % modulo == 0:
                self.disp_annotation()
            if self.countiter > 0 and (self.stop() or self.countiter < 4 
                              or self.countiter % modulo == 0): 
                if self.opts['verb_time']:
                    toc = self.elapsed_time()  
                    stime = str(int(toc//60))+':'+str(round(toc%60,1))
                else:
                    stime = ''
                print(' '.join((repr(self.countiter).rjust(5), 
                                repr(self.countevals).rjust(7), 
                                '%.15e' % (min(self.fit.fit)), 
                                '%4.1e' % (self.D.max()/self.D.min()), 
                                '%6.2e' % self.sigma, 
                                '%6.0e' % (self.sigma * sqrt(min(self.dC))), 
                                '%6.0e' % (self.sigma * sqrt(max(self.dC))), 
                                stime)))
                # if self.countiter < 4:
                sys.stdout.flush()
# end class cma 

#____________________________________________________________
#____________________________________________________________
#
class ElapsedTime(object):
    """32-bit C overflows after int(2**32/1e6) == 4294s about 72 min"""
    def __init__(self):
        self.tic0 = time.clock()
        self.tic = self.tic0
        self.lasttoc = time.clock()
        self.lastdiff = time.clock() - self.lasttoc
        self.time_to_add = 0
        self.messages = 0

    def __call__(self):
        toc = time.clock()
        if toc - self.tic >= self.lasttoc - self.tic: 
            self.lastdiff = toc - self.lasttoc
            self.lasttoc = toc
        else:  # overflow, reset self.tic
            if self.messages < 10:
                self.messages += 1
                print('timing overflow, last difference estimated:')
                print(self.tic0, self.tic, self.lasttoc, toc, toc - self.lasttoc, self.lastdiff)
            
            self.time_to_add += self.lastdiff + self.lasttoc - self.tic
            self.tic = toc  # reset 
            self.lasttoc = toc

        return toc - self.tic + self.time_to_add

#____________________________________________________________
#____________________________________________________________
#
class Options(dict):
    """``Options()`` returns a dictionary with the available options and their
    default values for function fmin and for class CMAEvolutionStrategy.
    
    ``Options(opts)`` returns the subset of recognized options in opts.  
    
    ``Options('pop')`` returns a subset of recognized options that contain
    'pop' in there keyword name, value or description. 
    
    Option values can be "written" in a string and, when passed to fmin
    or CMAEvolutionStrategy, are evaluated using "N" and "popsize" as 
    known values for dimension and population size (sample size, number 
    of new solutions per iteration). All default option values are such
    string. 

    Details
    -------
    All Options are originally defined via the input arguments of 
    `fmin()`.

    Options starting with ``tol`` are termination "tolerances". 

    For `tolstagnation`, the median over the first and the second half 
    of at least `tolstagnation` iterations are compared for both, the
    per-iteration best and per-iteration median function value. 
    Some options are, as mentioned (`restarts`,...), only used with `fmin`.
        
    Example
    -------
    ::
    
        import cma
        cma.Options('tol')  # shortcut for cma.Options().match('tol')
        o = cma.Options()
        o.match('verb')
    
    all options that contain 'tol' in their name or description. 

    :See: `fmin`(), `CMAEvolutionStrategy`, `CMAParameters`
     
    """
    
    # @classmethod # self is the class, not the instance
    # @property 
    # def default(self):
    #     """returns all options with defaults"""
    #     return fmin([],[])
     
    @staticmethod
    def defaults():
        """return a dictionary with default option values and description, 
        calls `fmin([], [])`"""
        return fmin([], [])

    @staticmethod
    def versatileOptions():
        """return list of options that can be changed at any time (not only be 
        initialized), however the list might not be entirely up to date. The 
        string ' #v ' in the default value indicates a 'versatile' option 
        that can be changed any time.
        
        """ 
        return tuple(sorted(i[0] for i in Options.defaults().items() if i[1].find(' #v ') > 0))
    
    def __init__(self, s=None, unchecked=False):
        """return an `Options` instance, either with the default options, 
        if ``s is None``, or with all options whose name or description
        contains `s`, if `s` is a string (case is disregarded), 
        or with entries from dictionary `s` as options, not complemented 
        with default options or settings
        
        Returns: see above.
        
        """
        # if not Options.defaults:  # this is different from self.defaults!!!
        #     Options.defaults = fmin([],[])
        if s is None:
            super(Options, self).__init__(Options.defaults()) 
            # self = Options.defaults() 
        elif type(s) is str:
            super(Options, self).__init__(Options().match(s)) 
            # we could return here
        else:
            super(Options, self).__init__(s) 
            
        if not unchecked: 
            for key in self.keys():
                if key not in Options.defaults():
                    print('Warning in cma.Options.__init__(): invalid key ``' + str(key) + '`` popped')
                    self.pop(key)
        # self.evaluated = False  # would become an option entry
                        
    def init(self, dict_or_str, val=None, warn=True):
        """initialize one or several options.  
        
        Arguments
        ---------
            `dict_or_str`
                a dictionary if ``val is None``, otherwise a key. 
                If `val` is provided `dict_or_str` must be a valid key. 
            `val`
                value for key

        Details
        -------
        Only known keys are accepted. Known keys are in `Options.defaults()`

        """
        #dic = dict_or_key if val is None else {dict_or_key:val} 
        dic = dict_or_str
        if val is not None: 
            dic = {dict_or_str:val}

        for key, val in dic.items():
            if key not in Options.defaults():
                # TODO: find a better solution? 
                if warn:
                    print('Warning in cma.Options.init(): key ' + 
                        str(key) + ' ignored')
            else:
                self[key] = val
                
        return self
            
    def set(self, dic, val=None, warn=True):
        """set can assign versatile options from `Options.versatileOptions()` 
        with a new value, use `init()` for the others. 
        
        Arguments
        ---------
            `dic`
                either a dictionary or a key. In the latter
                case, val must be provided
            `val`
                value for key
            `warn`
                bool, print a warning if the option cannot be changed
                and is therefore omitted
                
        This method will be most probably used with the ``opts`` attribute of 
        a `CMAEvolutionStrategy` instance. 

        """
        if val is not None:  # dic is a key in this case
            dic = {dic:val}  # compose a dictionary
        for key, val in dic.items():
            if key in Options.versatileOptions():
                self[key] = val
            elif warn:
                print('Warning in cma.Options.set(): key ' + str(key) + ' ignored')
        return self  # to allow o = Options(o).set(new)
    
    def complement(self):
        """add all missing options with their default values"""
        
        for key in Options.defaults():
            if key not in self: 
                self[key] = Options.defaults()[key]
        return self
                
    def settable(self):
        """return the subset of those options that are settable at any
        time. 
         
        Settable options are in `versatileOptions()`, but the 
        list might be incomlete. 
        
        """
        return Options([i for i in self.items() 
                                if i[0] in Options.versatileOptions()])

    def __call__(self, key, default=None, loc=None):
        """evaluate and return the value of option `key` on the fly, or
        returns those options whose name or description contains `key`, 
        case disregarded. 
        
        Details
        -------
        Keys that contain `filename` are not evaluated.
        For ``loc==None``, `self` is used as environment
        but this does not define `N`. 
            
        :See: `eval()`, `evalall()`
        
        """
        try:
            val = self[key]
        except:
            return self.match(key)
        
        if loc is None: 
            loc = self  # TODO: this hack is not so useful: popsize could be there, but N is missing
        try: 
            if type(val) is str:
                val = val.split('#')[0].strip()  # remove comments
                if key.find('filename') < 0:
                    val = eval(val, globals(), loc)
            # invoke default
            # TODO: val in ... fails with array type, because it is applied element wise!
            # elif val in (None,(),[],{}) and default is not None:  
            elif val is None and default is not None:  
                val = eval(str(default), globals(), loc)
        except:
            pass  # slighly optimistic: the previous is bug-free
        return val
        
    def eval(self, key, default=None, loc=None):
        """Evaluates and sets the specified option value in 
        environment `loc`. Many options need `N` to be defined in 
        `loc`, some need `popsize`. 

        Details
        -------
        Keys that contain 'filename' are not evaluated.
        For `loc` is None, the self-dict is used as environment
            
        :See: `evalall()`, `__call__`
        
        """
        self[key] = self(key, default, loc)
        return self[key]
        
    def evalall(self, loc=None):
        """Evaluates all option values in environment `loc`. 
        
        :See: `eval()`
        
        """
        # TODO: this needs rather the parameter N instead of loc
        if 'N' in loc.keys():  # TODO: __init__ of CMA can be simplified
            popsize = self('popsize', Options.defaults()['popsize'], loc)
            for k in self.keys():
                self.eval(k, Options.defaults()[k], 
                          {'N':loc['N'], 'popsize':popsize})
        return self
    
    def match(self, s=''):
        """return all options that match, in the name or the description, 
        with string `s`, case is disregarded. 
        
        Example: ``cma.Options().match('verb')`` returns the verbosity options. 
        
        """
        match = s.lower()
        res = {}
        for k in sorted(self):
            s = str(k) + '=\'' + str(self[k]) + '\''
            if match in s.lower():
                res[k] = self[k]
        return Options(res)
        
    def pp(self):
        pprint(self)
        
    def printme(self, linebreak=80):
        for i in sorted(Options.defaults().items()): 
            s = str(i[0]) + "='" + str(i[1]) + "'"
            a = s.split(' ')

            # print s in chunks
            l = ''  # start entire to the left
            while a:
                while a and len(l) + len(a[0]) < linebreak:
                    l += ' ' + a.pop(0)
                print(l)
                l = '        '  # tab for subsequent lines

    #____________________________________________________________
    #____________________________________________________________
    # 
    def timesCroot(self, mat):
        """return C**0.5 times mat, where mat can be a vector or matrix.
        Not functional, because _Croot=C**0.5 is never computed (should be in updateBD)
        """
        if self.opts['CMA_diagonal'] is True \
                       or self.countiter <= self.opts['CMA_diagonal']: 
            res = (self._Croot * mat.T).T
        else:
            res = np.dot(self._Croot, mat)
        return res
    def divCroot(self, mat):
        """return C**-1/2 times mat, where mat can be a vector or matrix"""
        if self.opts['CMA_diagonal'] is True \
                       or self.countiter <= self.opts['CMA_diagonal']: 
            res = (self._Crootinv * mat.T).T
        else:
            res = np.dot(self._Crootinv, mat)
        return res

#____________________________________________________________
#____________________________________________________________
class CMAParameters(object):
    """strategy parameters like population size and learning rates.
    
    Note: 
        contrary to `Options`, `CMAParameters` is not (yet) part of the 
        "user-interface" and subject to future changes (it might become
        a `collections.namedtuple`)
    
    Example
    ------- 
    >>> import cma
    >>> es = cma.CMAEvolutionStrategy(20 * [0.1], 1)
    (6_w,12)-CMA-ES (mu_w=3.7,w_1=40%) in dimension 20 (seed=504519190)  # the seed is "random" by default
    >>>
    >>> type(es.sp)  # sp contains the strategy parameters
    <class 'cma.CMAParameters'>
    >>>
    >>> es.sp.disp()
    {'CMA_on': True,
     'N': 20,
     'c1': 0.004181139918745593,
     'c1_sep': 0.034327992810300939,
     'cc': 0.17176721127681213,
     'cc_sep': 0.25259494835857677,
     'cmean': 1.0,
     'cmu': 0.0085149624979034746,
     'cmu_sep': 0.057796356229390715,
     'cs': 0.21434997799189287,
     'damps': 1.2143499779918929,
     'mu': 6,
     'mu_f': 6.0,
     'mueff': 3.7294589343030671,
     'popsize': 12,
     'rankmualpha': 0.3,
     'weights': array([ 0.40240294,  0.25338908,  0.16622156,  0.10437523,  0.05640348,
            0.01720771])}
    >>>
    >> es.sp == cma.CMAParameters(20, 12, cma.Options().evalall({'N': 20}))
    True
    
    :See: `Options`, `CMAEvolutionStrategy`
    
    """
    def __init__(self, N, opts, ccovfac=1, verbose=True):
        """Compute strategy parameters, mainly depending on
        dimension and population size, by calling `set` 
        
        """
        self.N = N
        if ccovfac == 1:
            ccovfac = opts['CMA_on']  # that's a hack
        self.set(opts, ccovfac=ccovfac, verbose=verbose)
        
    def set(self, opts, popsize=None, ccovfac=1, verbose=True):
        """Compute strategy parameters as a function
        of dimension and population size """

        alpha_cc = 0.0  # TODO: should be come 1.0 in near future
        
        def cone(df, mu, N, alphacov=2.0):
            """rank one update learning rate, ``df`` is disregarded and obsolete, reduce alphacov on noisy problems, say to 0.5"""
            return alphacov / ((N + 1.3)**2 + mu)

        def cmu(df, mu, alphamu=0.0, alphacov=2.0):
            """rank mu learning rate, disregarding the constrant cmu <= 1 - cone"""
            return alphacov * (alphamu + mu - 2 + 1/mu) / ((N + 2)**2 + alphacov * mu / 2)
            
        def conedf(df, mu, N):
            """used for computing separable learning rate"""
            return 1. / (df + 2.*sqrt(df) + float(mu)/N)

        def cmudf(df, mu, alphamu):
            """used for computing separable learning rate"""
            return (alphamu + mu - 2. + 1./mu) / (df + 4.*sqrt(df) + mu/2.)

        sp = self
        N = sp.N
        if popsize:
            opts.evalall({'N':N, 'popsize':popsize})
        else:
            popsize = opts.evalall({'N':N})['popsize']  # the default popsize is computed in Options()
        sp.popsize = popsize 
        if opts['CMA_mirrors'] < 0.5:
            sp.lam_mirr = int(0.5 + opts['CMA_mirrors'] * popsize)
        elif opts['CMA_mirrors'] > 1:
            sp.lam_mirr = int(0.5 + opts['CMA_mirrors'])
        else:
            sp.lam_mirr = int(0.5 + 0.16 * min((popsize, 2 * N + 2)) + 0.29)  # 0.158650... * popsize is optimal
            # lam = arange(2,22)
            # mirr = 0.16 + 0.29/lam
            # print(lam); print([int(0.5 + l) for l in mirr*lam])
            # [ 2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21]
            # [1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 4]
            
        sp.mu_f = sp.popsize / 2.0  # float value of mu
        if opts['CMA_mu'] is not None: 
            sp.mu_f = opts['CMA_mu']
        sp.mu = int(sp.mu_f + 0.499999) # round down for x.5
        # in principle we have mu_opt = popsize/2 + lam_mirr/2, 
        # which means in particular weights should only be negative for q > 0.5+mirr_frac/2
        if sp.mu > sp.popsize - 2 * sp.lam_mirr + 1:
            print("WARNING: pairwise selection is not implemented, therefore " + 
                  " mu = %d > %d = %d - 2*%d + 1 = popsize - 2*mirr + 1 can produce a bias" % (
                    sp.mu, sp.popsize - 2 * sp.lam_mirr + 1, sp.popsize, sp.lam_mirr))
        if sp.lam_mirr > sp.popsize // 2:
            raise _Error("fraction of mirrors in the population as read from option CMA_mirrors cannot be larger 0.5, " +
                         "theoretically optimal is 0.16") 
        sp.weights = log(sp.mu_f + 0.5) - log(1 + np.arange(sp.mu))
        sp.weights /= sum(sp.weights)
        sp.mueff = 1 / sum(sp.weights**2)
        sp.cs = (sp.mueff + 2) / (N + sp.mueff + 3)
        sp.cc = (4 + alpha_cc * sp.mueff / N) / (N + 4 + alpha_cc * 2 * sp.mueff / N)  
        sp.cc_sep = (1 + 1/N + alpha_cc * sp.mueff / N) / (N**0.5 + 1/N + alpha_cc * 2 * sp.mueff / N) # \not\gg\cc
        sp.rankmualpha = opts['CMA_rankmualpha']
        # sp.rankmualpha = _evalOption(opts['CMA_rankmualpha'], 0.3)
        sp.c1 = ccovfac * min(1, sp.popsize/6) * cone((N**2 + N) / 2, sp.mueff, N) # 2. / ((N+1.3)**2 + sp.mucov)
        sp.c1_sep = ccovfac * conedf(N, sp.mueff, N) 
        if 11 < 3:
            sp.c1 = 0.
            print('c1 is zero')
        if opts['CMA_rankmu'] != 0:  # also empty
            sp.cmu = min(1 - sp.c1, ccovfac * cmu((N**2+N)/2, sp.mueff, sp.rankmualpha))
            sp.cmu_sep = min(1 - sp.c1_sep, ccovfac * cmudf(N, sp.mueff, sp.rankmualpha))
        else:
            sp.cmu = sp.cmu_sep = 0

        sp.neg = BlancClass()
        if opts['CMA_active']:
            # in principle we have mu_opt = popsize/2 + lam_mirr/2, 
            # which means in particular weights should only be negative for q > 0.5+mirr_frac/2
            sp.neg.mu_f = popsize - (popsize + sp.lam_mirr) / 2  if popsize > 2 else 1
            sp.neg.weights = log(sp.mu_f + 0.5) - log(1 + np.arange(sp.popsize - int(sp.neg.mu_f), sp.popsize))
            sp.neg.mu = len(sp.neg.weights)  # maybe never useful?
            sp.neg.weights /= sum(sp.neg.weights)
            sp.neg.mueff = 1 / sum(sp.neg.weights**2)
            sp.neg.cmuexp = opts['CMA_activefac'] * 0.25 * sp.neg.mueff / ((N+2)**1.5 + 2 * sp.neg.mueff)
            assert sp.neg.mu >= sp.lam_mirr  # not really necessary
            # sp.neg.minresidualvariance = 0.66  # not it use, keep at least 0.66 in all directions, small popsize is most critical  
        else:
            sp.neg.cmuexp = 0   

        sp.CMA_on = sp.c1 + sp.cmu > 0  
        # print(sp.c1_sep / sp.cc_sep)

        if not opts['CMA_on'] and opts['CMA_on'] not in (None,[],(),''): 
            sp.CMA_on = False
            # sp.c1 = sp.cmu = sp.c1_sep = sp.cmu_sep = 0
        sp.damps = opts['CMA_dampfac'] * (1 - 0*sp.lam_mirr/sp.popsize + 2*max(0,sqrt((sp.mueff-1)/(N+1))-1) + sp.cs) 
        # sp.damps = opts['CMA_dampfac'] * (2 - 0*sp.lam_mirr/sp.popsize) * sp.mueff/sp.popsize + 0.3 + sp.cs  # nicer future setting 
        if 11 < 3:
            sp.damps = 10 * sp.damps  # 1e99 # (1 + 2*max(0,sqrt((sp.mueff-1)/(N+1))-1)) + sp.cs; 
            # sp.damps = 20 # 1. + 20 * sp.cs**-1  # 1e99 # (1 + 2*max(0,sqrt((sp.mueff-1)/(N+1))-1)) + sp.cs; 
            print('damps is %f' % (sp.damps))

        sp.cmean = float(opts['CMA_cmean'])
        # sp.kappa = 1  # 4-D, lam=16, rank1, kappa < 4 does not influence convergence rate
                        # in larger dim it does, 15-D with defaults, kappa=8 factor 2 
        if sp.cmean != 1:
            print('  cmean = %f' % (sp.cmean))
        
        if verbose:
            if not sp.CMA_on:
                print('covariance matrix adaptation turned off')
            if opts['CMA_mu'] != None: 
                print('mu = %f' % (sp.mu_f))

        # return self  # the constructor returns itself
        
    def disp(self):
        pprint(self.__dict__)
            
#____________________________________________________________
#____________________________________________________________
class CMAStopDict(dict):
    """keep and update a termination condition dictionary, which is
    "usually" empty and returned by `CMAEvolutionStrategy.stop()`.
    
    Details
    -------
    This could be a nested class, but nested classes cannot be serialized.  
    
    :See: `stop()`
    
    """
    def __init__(self, d={}):
        update = (type(d) == CMAEvolutionStrategy)
        inherit = (type(d) == CMAStopDict)           
        super(CMAStopDict, self).__init__({} if update else d)
        self._stoplist = d._stoplist if inherit else []    # multiple entries
        self.lastiter = d.lastiter if inherit else 0  # probably not necessary
        if update:
            self._update(d)
        
    def __call__(self, es):
        """update the dictionary"""
        return self._update(es)
 
    def _addstop(self, key, cond, val=None):
        if cond:
            self.stoplist.append(key)  # can have the same key twice
            if key in self.opts.keys():
                val = self.opts[key]
            self[key] = val
    
    def _update(self, es):
        """Test termination criteria and update dictionary. 

        """
        if es.countiter == self.lastiter:
            if es.countiter == 0:
                self.__init__()
                return self
            try:
                if es == self.es:
                    return self
            except: # self.es not yet assigned
                pass
            
        self.lastiter = es.countiter
        self.es = es

        self.stoplist = []
        
        N = es.N
        opts = es.opts
        self.opts = opts  # a hack to get _addstop going
        
        # fitness: generic criterion, user defined w/o default
        self._addstop('ftarget',
                     es.best.f < opts['ftarget'])
        # maxiter, maxfevals: generic criteria
        self._addstop('maxfevals',
                     es.countevals - 1 >= opts['maxfevals'])
        self._addstop('maxiter',
                     es.countiter >= opts['maxiter'])
        # tolx, tolfacupx: generic criteria
        # tolfun, tolfunhist (CEC:tolfun includes hist) 
        self._addstop('tolx',
                     all([es.sigma*xi < opts['tolx'] for xi in es.pc]) and \
                     all([es.sigma*xi < opts['tolx'] for xi in sqrt(es.dC)]))
        self._addstop('tolfacupx',
                     any([es.sigma * sig > es.sigma0 * opts['tolfacupx']
                          for sig in sqrt(es.dC)]))
        self._addstop('tolfun',
                     es.fit.fit[-1] - es.fit.fit[0] < opts['tolfun'] and \
                     max(es.fit.hist) - min(es.fit.hist) < opts['tolfun'])
        self._addstop('tolfunhist',
                     len(es.fit.hist) > 9 and \
                     max(es.fit.hist) - min(es.fit.hist) <  opts['tolfunhist'])

        # worst seen false positive: table N=80,lam=80, getting worse for fevals=35e3 \approx 50 * N**1.5 
        # but the median is not so much getting worse
        # / 5 reflects the sparsity of histbest/median
        # / 2 reflects the left and right part to be compared
        l = int(max(opts['tolstagnation'] / 5. / 2, len(es.fit.histbest) / 10));
        # TODO: the problem in the beginning is only with best ==> ???
        if 11 < 3:  # 
            print(es.countiter, (opts['tolstagnation'], es.countiter > N * (5 + 100 / es.popsize),  
                        len(es.fit.histbest) > 100,
                        np.median(es.fit.histmedian[:l]) >= np.median(es.fit.histmedian[l:2*l]),
                        np.median(es.fit.histbest[:l]) >= np.median(es.fit.histbest[l:2*l])))
        # equality should handle flat fitness
        self._addstop('tolstagnation', # leads sometimes early stop on ftablet, fcigtab, N>=50?
                    1 < 3 and opts['tolstagnation'] and es.countiter > N * (5 + 100 / es.popsize) and  
                    len(es.fit.histbest) > 100 and 2*l < len(es.fit.histbest) and 
                    np.median(es.fit.histmedian[:l]) >= np.median(es.fit.histmedian[l:2*l]) and
                    np.median(es.fit.histbest[:l]) >= np.median(es.fit.histbest[l:2*l]))
        # iiinteger: stagnation termination can prevent to find the optimum
        
        self._addstop('tolupsigma', 
                      es.sigma / es.sigma0 / np.max(es.D) > opts['tolupsigma'])

        if 11 < 3 and 2*l < len(es.fit.histbest):  # TODO: this might go wrong, because the nb of written columns changes
            tmp = np.array((-np.median(es.fit.histmedian[:l]) + np.median(es.fit.histmedian[l:2*l]),
                        -np.median(es.fit.histbest[:l]) + np.median(es.fit.histbest[l:2*l])))
            es.more_to_write = [(10**t if t < 0 else t + 1) for t in tmp] # the latter to get monotonicy

        # TODO: add option stoponequalfunvals and test here...

        if 1 < 3:
            # non-user defined, method specific
            # noeffectaxis (CEC: 0.1sigma), noeffectcoord (CEC:0.2sigma), conditioncov
            self._addstop('noeffectcoord',
                         any([es.mean[i] == es.mean[i] + 0.2*es.sigma*sqrt(es.dC[i])
                              for i in xrange(N)]))
            if opts['CMA_diagonal'] is not True and es.countiter > opts['CMA_diagonal']: 
                self._addstop('noeffectaxis',
                             all([es.mean[i] == es.mean[i] +
                                  0.1*es.sigma * es.D[es.countiter%N]
                                  * es.B[i,es.countiter%N] for i in xrange(N)]))
            self._addstop('conditioncov',
                         es.D[-1] > 1e7 * es.D[0], 1e14)  # TODO 

            self._addstop('callback',
                          es.callbackstop)  # termination_callback
        if len(self):
            self._addstop('flat fitness: please (re)consider how to compute the fitness more elaborate',
                         len(es.fit.hist) > 9 and \
                         max(es.fit.hist) == min(es.fit.hist))

        return self

#_____________________________________________________________________
#_____________________________________________________________________
#
class OptimDataLogger(object):
    """"abstract" base class for a data logger that can be used with an `OOOptimizer`"""
    def add(self, optim=None, more_data=[]):
        """abstract method, add a "data point" from the state of `optim` into the 
        logger, the argument `optim` can be omitted if it was `register()`-ed before"""
        OOOptimizer.abstract()
    def register(self, optim):
        """abstract method, register an optimizer `optim`, only needed if `add()` is 
        called without a value for the `optim` argument"""
        self.optim = optim 
    def disp(self):
        """display some data trace (not implemented)"""
        print('method OptimDataLogger.disp() not implemented, to be done in subclass ' + str(type(self)))
    def plot(self):
        """plot data (not implemented)"""
        print('method OptimDataLogger.plot() is not implemented, to be done in subclass ' + str(type(self)))
    def data(self):
        """return logged data in a dictionary (not implemented)"""
        print('method OptimDataLogger.data() is not implemented, to be done in subclass ' + str(type(self)))
                
#_____________________________________________________________________
#_____________________________________________________________________
#
class CMADataLogger(OptimDataLogger):  # might become a dict at some point
    """data logger for class `CMAEvolutionStrategy`. The logger is
    identified by its name prefix and writes or reads according 
    data files. 
    
    Examples
    ======== 
    ::
    
        import cma
        es = cma.CMAEvolutionStrategy(...)
        data = cma.CMADataLogger().register(es)
        while not es.stop():
            ...
            data.add()  # add can also take an argument

        data.plot() # or a short cut can be used:
        cma.plot()  # plot data from logger with default name
        
        
        data2 = cma.CMADataLogger(another_filename_prefix).load()
        data2.plot()
        data2.disp()

    ::

        import cma
        from pylab import *
        res = cma.fmin(cma.Fcts.sphere, rand(10), 1e-0)
        dat = res[-1]  # the CMADataLogger
        dat.load()  # by "default" data are on disk
        semilogy(dat.f[:,0], dat.f[:,5])  # plot f versus iteration, see file header
        show()
        
    Details
    =======
    After loading data, the logger has the attributes `xmean`, `xrecent`, `std`, `f`, and `D`, 
    corresponding to xmean, xrecentbest, stddev, fit, and axlen filename trails. 

    :See: `disp()`, `plot()`
    
    """
    default_prefix = 'outcmaes'
    names = ('axlen','fit','stddev','xmean','xrecentbest')
    
    def __init__(self, name_prefix=default_prefix, modulo=None):
        """initialize logging of data from a `CMAEvolutionStrategy` instance, 
        default modulo expands to 1 == log with each call
        
        """
        # super(CMAData, self).__init__({'iter':[], 'stds':[], 'D':[], 'sig':[], 'fit':[], 'xm':[]})
        # class properties:
        self.counter = 0  # number of calls of add
        self.modulo = 1 if modulo is None else modulo  # allows calling with None
        self.name_prefix = name_prefix if name_prefix else CMADataLogger.default_prefix 
        if type(self.name_prefix) == CMAEvolutionStrategy:
            self.name_prefix = self.name_prefix.opts.eval('verb_filenameprefix') 

    def register(self, es, append=False, modulo=None):
        """register a `CMAEvolutionStrategy` instance for logging, 
        ``append=True`` appends to previous data logged under the same name, 
        by default previous data are overwritten. 
        
        """
        if type(es) != CMAEvolutionStrategy:
            raise TypeError("only class CMAEvolutionStrategy can be registered for logging")
        self.es = es
        if modulo is not None:
            self.modulo = modulo
        if not append and self.modulo != 0:
            self.initialize()  # write file headers
        return self

    def initialize(self, modulo=None):
        """reset logger, overwrite original files, `modulo`: log only every modulo call"""
        if modulo is not None:
            self.modulo = modulo 
        try:
            es = self.es  # must have been registered
        except AttributeError :
            raise _Error('call register() before init()')
        
        # write headers for output
        fn = self.name_prefix + 'fit.dat'
        strseedtime = 'seed=%d, %s' % (es.opts['seed'], time.asctime())

        try: 
            with open(fn, 'w') as f:
                f.write('% # columns="iteration, evaluation, sigma, axis ratio, ' +
                        'bestever, best, median, worst objective function value, ' +
                        'further objective values of best", ' +
                        strseedtime + 
                        # strftime("%Y/%m/%d %H:%M:%S", localtime()) + # just asctime() would do
                        '\n')
        except (IOError, OSError):
            print('could not open file ' + fn)

        fn = self.name_prefix + 'axlen.dat'
        try: 
            f = open(fn, 'w')
            f.write('%  columns="iteration, evaluation, sigma, max axis length, ' +
                    ' min axis length, all principle axes lengths ' +
                    ' (sorted square roots of eigenvalues of C)", ' +
                    strseedtime + 
                    '\n')
            f.close()
        except (IOError, OSError):
            print('could not open file ' + fn)
        finally:
            f.close()
        fn = self.name_prefix + 'stddev.dat'
        try: 
            f = open(fn, 'w')
            f.write('% # columns=["iteration, evaluation, sigma, void, void, ' +
                    ' stds==sigma*sqrt(diag(C))", ' +
                    strseedtime + 
                    '\n')
            f.close()
        except (IOError, OSError):
            print('could not open file ' + fn)
        finally:
            f.close()

        fn = self.name_prefix + 'xmean.dat'
        try: 
            with open(fn, 'w') as f:
                f.write('% # columns="iteration, evaluation, void, void, void, xmean", ' +
                        strseedtime)
                f.write(' # scaling_of_variables: ')
                if np.size(es.gp.scales) > 1:
                    f.write(' '.join(map(str, es.gp.scales)))
                else:
                    f.write(str(es.gp.scales))
                f.write(', typical_x: ')
                if np.size(es.gp.typical_x) > 1:
                    f.write(' '.join(map(str, es.gp.typical_x)))
                else:
                    f.write(str(es.gp.typical_x))
                f.write('\n')
                f.close()
        except (IOError, OSError):
            print('could not open/write file ' + fn)

        fn = self.name_prefix + 'xrecentbest.dat'
        try: 
            with open(fn, 'w') as f: 
                f.write('% # iter+eval+sigma+0+fitness+xbest, ' + 
                        strseedtime + 
                        '\n')
        except (IOError, OSError):
            print('could not open/write file ' + fn)
            
        return self
    # end def __init__
    
    def load(self, filenameprefix=None):
        """loads data from files written and return a data dictionary, *not*
        a prerequisite for using `plot()` or `disp()`. 
    
        Argument `filenameprefix` is the filename prefix of data to be loaded (five files),
        by default ``'outcmaes'``. 
        
        Return data dictionary with keys `xrecent`, `xmean`, `f`, `D`, `std`
        
        """
        if not filenameprefix:
            filenameprefix = self.name_prefix
        dat = self  # historical 
        dat.xrecent = _fileToMatrix(filenameprefix + 'xrecentbest.dat') 
        dat.xmean = _fileToMatrix(filenameprefix + 'xmean.dat')
        dat.std = _fileToMatrix(filenameprefix + 'stddev' + '.dat')
        # a hack to later write something into the last entry
        for key in ['xmean', 'xrecent', 'std']:
            dat.__dict__[key].append(dat.__dict__[key][-1])
            dat.__dict__[key] = array(dat.__dict__[key])
        dat.f = array(_fileToMatrix(filenameprefix + 'fit.dat'))
        dat.D = array(_fileToMatrix(filenameprefix + 'axlen' + '.dat'))
        return dat

         
    def add(self, es=None, more_data=[], modulo=None): # TODO: find a different way to communicate current x and f
        """append some logging data from `CMAEvolutionStrategy` class instance `es`, 
        if ``number_of_times_called % modulo`` equals to zero, never if ``modulo==0``. 
        
        The sequence ``more_data`` must always have the same length.  
        
        """
        self.counter += 1 
        mod = modulo if modulo is not None else self.modulo
        if mod == 0 or (self.counter > 3 and self.counter % mod):
            return
        if es is None:
            try:
                es = self.es  # must have been registered
            except AttributeError :
                raise _Error('call register() before add()')
        if type(es) is not CMAEvolutionStrategy:
            raise TypeError('<type \'CMAEvolutionStrategy\'> expected, found '
                            + str(type(es)) + ' in CMADataLogger.add')
        
        if 1 < 3:
            try: # TODO: find a more decent interface to store and pass recent_x
                xrecent = es.out['recent_x']
            except:
                if self.counter == 2:  # by now a recent_x should be available
                    print('WARNING: es.out[\'recent_x\'] not found in CMADataLogger.add, count=' 
                          + str(self.counter))
        try: 
            # fit
            if es.countiter > 0:
                fit = es.fit.fit[0]
                if es.fmean_noise_free != 0:
                    fit = es.fmean_noise_free
                fn = self.name_prefix + 'fit.dat'
                with open(fn, 'a') as f:
                    f.write(str(es.countiter) + ' ' 
                            + str(es.countevals) + ' '
                            + str(es.sigma) + ' '
                            + str(es.D.max()/es.D.min()) + ' '
                            + str(es.best.f) + ' '
                            + '%.16e' % fit + ' '
                            + str(es.fit.fit[es.sp.popsize//2]) + ' '
                            + str(es.fit.fit[-1]) + ' '
                            # + str(es.sp.popsize) + ' '
                            # + str(10**es.noiseS) + ' '
                            # + str(es.sp.cmean) + ' '
                            + ' '.join(str(i) for i in es.more_to_write)
                            + ' '.join(str(i) for i in more_data)
                            + '\n')
            # axlen
            fn = self.name_prefix + 'axlen.dat'
            with open(fn, 'a') as f:  # does not rely on reference counting
                f.write(str(es.countiter) + ' ' 
                        + str(es.countevals) + ' '
                        + str(es.sigma) + ' '
                        + str(es.D.max()) + ' '
                        + str(es.D.min()) + ' '
                        + ' '.join(map(str, es.D))
                        + '\n')
            # stddev
            fn = self.name_prefix + 'stddev.dat'
            with open(fn, 'a') as f:
                f.write(str(es.countiter) + ' ' 
                        + str(es.countevals) + ' '
                        + str(es.sigma) + ' '
                        + '0 0 '
                        + ' '.join(map(str, es.sigma*sqrt(es.dC)))
                        + '\n')
            # xmean
            fn = self.name_prefix + 'xmean.dat'
            with open(fn, 'a') as f:
                if es.countevals < es.sp.popsize:
                    f.write('0 0 0 0 0 '
                            + ' '.join(map(str,
                                              # TODO should be optional the phenotyp?
                                              es.gp.geno(es.x0)))  
                            + '\n')
                else:
                    f.write(str(es.countiter) + ' ' 
                            + str(es.countevals) + ' '
                            # + str(es.sigma) + ' '
                            + '0 '
                            + str(es.fmean_noise_free) + ' '
                            + str(es.fmean) + ' '  # TODO: this does not make sense
                            # TODO should be optional the phenotyp?
                            + ' '.join(map(str, es.mean)) 
                            + '\n')
            # xrecent
            fn = self.name_prefix + 'xrecentbest.dat'
            if es.countiter > 0 and xrecent is not None:
                with open(fn, 'a') as f:
                    f.write(str(es.countiter) + ' ' 
                            + str(es.countevals) + ' '
                            + str(es.sigma) + ' '
                            + '0 '
                            + str(es.fit.fit[0]) + ' '
                            + ' '.join(map(str, xrecent))
                            + '\n')
                    
        except (IOError, OSError):
            if es.countiter == 1:
                print('could not open/write file') 

    def closefig(self):
        pylab.close(self.fighandle)
        
    def save(self, nameprefix, switch=False):
        """saves logger data to a different set of files, for 
        ``switch=True`` also the loggers name prefix is switched to 
        the new value
        
        """
        if not nameprefix or type(nameprefix) is not str:
            _Error('filename prefix must be a nonempty string')

        if nameprefix == self.default_prefix:
            _Error('cannot save to default name "' + nameprefix + '...", chose another name')
            
        if nameprefix == self.name_prefix:
            return
        
        for name in CMADataLogger.names:
            open(nameprefix+name+'.dat', 'w').write(open(self.name_prefix+name+'.dat').read())
            
        if switch:
            self.name_prefix = nameprefix
       
    def plot(self, fig=None, iabscissa=1, iteridx=None, plot_mean=True,  # TODO: plot_mean default should be False
             foffset=1e-19, x_opt = None, fontsize=10):
        """
        plot data from a `CMADataLogger` (using the files written by the logger). 
    
        Arguments
        ---------
            `fig`
                figure number, by default 325
            `iabscissa`
                ``0==plot`` versus iteration count,
                ``1==plot`` versus function evaluation number
            `iteridx`
                iteration indices to plot
                
        Return `CMADataLogger` itself. 
        
        Examples
        --------
        ::
        
            import cma
            logger = cma.CMADataLogger()  # with default name
            # try to plot the "default logging" data (e.g. from previous fmin calls)
            logger.plot() # to continue you might need to close the pop-up window 
                          # once and call plot() again. 
                          # This behavior seems to disappear in subsequent 
                          # calls of plot(). Also using ipython with -pylab
                          # option might help. 
            cma.savefig('fig325.png')  # save current figure
            logger.closefig()
    
        Dependencies: matlabplotlib/pylab. 
        
        """
        
        dat = self.load(self.name_prefix)
        
        try: 
            # pylab: prodedural interface for matplotlib
            from  matplotlib.pylab import figure, ioff, ion, subplot, semilogy, hold, plot, grid, \
                 axis, title, text, xlabel, isinteractive, draw, gcf
                 
        except ImportError:
            ImportError('could not find matplotlib.pylab module, function plot() is not available')
            return
            
        if fontsize and pylab.rcParams['font.size'] != fontsize:
            print('global variable pylab.rcParams[\'font.size\'] set (from ' + 
                  str(pylab.rcParams['font.size']) + ') to ' + str(fontsize))
            pylab.rcParams['font.size'] = fontsize  # subtracted in the end, but return can happen inbetween
        
        if fig:
            figure(fig)
        else:
            figure(325)
            # show()  # should not be necessary
        self.fighandle = gcf()  # fighandle.number  

        if iabscissa not in (0,1):
            iabscissa = 1
        interactive_status = isinteractive()
        ioff() # prevents immediate drawing
        
        dat.x = dat.xrecent
        if len(dat.x) < 2:
            print('not enough data to plot')
            return {}

        if plot_mean:
            dat.x = dat.xmean    # this is the genotyp
        if iteridx is not None:
            dat.f = dat.f[np.where(map(lambda x: x in iteridx, dat.f[:,0]))[0],:]
            dat.D = dat.D[np.where(map(lambda x: x in iteridx, dat.D[:,0]))[0],:]
            iteridx.append(dat.x[-1,1])  # last entry is artificial
            dat.x = dat.x[np.where(map(lambda x: x in iteridx, dat.x[:,0]))[0],:]
            dat.std = dat.std[np.where(map(lambda x: x in iteridx, dat.std[:,0]))[0],:]

        if iabscissa == 0:
            xlab = 'iterations' 
        elif iabscissa == 1:
            xlab = 'function evaluations' 
        
        # use fake last entry in x and std for line extension-annotation
        if dat.x.shape[1] < 100:
            minxend = int(1.06*dat.x[-2, iabscissa])
            # write y-values for individual annotation into dat.x
            dat.x[-1, iabscissa] = minxend  # TODO: should be ax[1]
            idx = np.argsort(dat.x[-2,5:])
            idx2 = np.argsort(idx)
            if x_opt is None: 
                dat.x[-1,5+idx] = np.linspace(np.min(dat.x[:,5:]), 
                            np.max(dat.x[:,5:]), dat.x.shape[1]-5)
            else: 
                dat.x[-1,5+idx] = np.logspace(np.log10(np.min(abs(dat.x[:,5:]))), 
                            np.log10(np.max(abs(dat.x[:,5:]))), dat.x.shape[1]-5)
        else:
            minxend = 0

        if len(dat.f) == 0:
            print('nothing to plot')
            return
            
        # not in use anymore, see formatter above
        # xticklocs = np.arange(5) * np.round(minxend/4., -int(np.log10(minxend/4.)))
        
        # dfit(dfit<1e-98) = NaN;
        
        ioff() # turns update off

        # TODO: if abscissa==0 plot in chunks, ie loop over subsets where dat.f[:,0]==countiter is monotonous 

        subplot(2,2,1)
        self.plotdivers(dat, iabscissa, foffset)
        
        # TODO: modularize also the remaining subplots
        subplot(2,2,2)
        hold(False)
        if x_opt is not None:  # TODO: differentate neg and pos?
            semilogy(dat.x[:, iabscissa], abs(dat.x[:,5:]) - x_opt, '-')
        else:
            plot(dat.x[:, iabscissa], dat.x[:,5:],'-')
        hold(True)
        grid(True)
        ax = array(axis())
        # ax[1] = max(minxend, ax[1]) 
        axis(ax)
        ax[1] -= 1e-6 
        if dat.x.shape[1] < 100:
            yy = np.linspace(ax[2]+1e-6, ax[3]-1e-6, dat.x.shape[1]-5)
            #yyl = np.sort(dat.x[-1,5:])
            idx = np.argsort(dat.x[-1,5:])
            idx2 = np.argsort(idx)
            if x_opt is not None:
                semilogy([dat.x[-1, iabscissa], ax[1]], [abs(dat.x[-1,5:]), yy[idx2]], 'k-') # line from last data point
                semilogy(np.dot(dat.x[-2, iabscissa],[1,1]), array([ax[2]+1e-6, ax[3]-1e-6]), 'k-')
            else:
                # plot([dat.x[-1, iabscissa], ax[1]], [dat.x[-1,5:], yy[idx2]], 'k-') # line from last data point
                plot(np.dot(dat.x[-2, iabscissa],[1,1]), array([ax[2]+1e-6, ax[3]-1e-6]), 'k-')
            # plot(array([dat.x[-1, iabscissa], ax[1]]),
            #      reshape(array([dat.x[-1,5:], yy[idx2]]).flatten(), (2,4)), '-k')
            for i in range(len(idx)):
                # TODOqqq: annotate phenotypic value!?
                # text(ax[1], yy[i], 'x(' + str(idx[i]) + ')=' + str(dat.x[-2,5+idx[i]]))
                text(dat.x[-1,iabscissa], dat.x[-1,5+i], 'x(' + str(i) + ')=' + str(dat.x[-2,5+i]))
                
        i = 2  # find smallest i where iteration count differs (in case the same row appears twice)
        while i < len(dat.f) and dat.f[-i][0] == dat.f[-1][0]:
            i += 1
        title('Object Variables (' + ('mean' if plot_mean else 'curr best') + 
                ', ' + str(dat.x.shape[1]-5) + '-D, popsize~' + 
                (str(int((dat.f[-1][1] - dat.f[-i][1]) / (dat.f[-1][0] - dat.f[-i][0]))) 
                    if len(dat.f.T[0]) > 1 and dat.f[-1][0] > dat.f[-i][0] else 'NA') 
                + ')')
        # pylab.xticks(xticklocs)
        
        # Scaling
        subplot(2,2,3)
        hold(False)
        semilogy(dat.D[:, iabscissa], dat.D[:,5:], '-b')
        hold(True)
        grid(True)
        ax = array(axis()) 
        # ax[1] = max(minxend, ax[1]) 
        axis(ax) 
        title('Scaling (All Main Axes)') 
        # pylab.xticks(xticklocs)
        xlabel(xlab) 
        
        # standard deviations
        subplot(2,2,4)
        hold(False)
        # remove sigma from stds (graphs become much better readible)
        dat.std[:,5:] = np.transpose(dat.std[:,5:].T / dat.std[:,2].T) 
        # ax = array(axis())
        # ax[1] = max(minxend, ax[1]) 
        # axis(ax)
        if 1 < 2 and dat.std.shape[1] < 100:
            # use fake last entry in x and std for line extension-annotation
            minxend = int(1.06*dat.x[-2, iabscissa])
            dat.std[-1, iabscissa] = minxend  # TODO: should be ax[1]
            idx = np.argsort(dat.std[-2,5:])
            idx2 = np.argsort(idx)
            dat.std[-1,5+idx] = np.logspace(np.log10(np.min(dat.std[:,5:])), 
                            np.log10(np.max(dat.std[:,5:])), dat.std.shape[1]-5)

            dat.std[-1, iabscissa] = minxend  # TODO: should be ax[1]
            yy = np.logspace(np.log10(ax[2]), np.log10(ax[3]), dat.std.shape[1]-5)
            #yyl = np.sort(dat.std[-1,5:])
            idx = np.argsort(dat.std[-1,5:])
            idx2 = np.argsort(idx)
            # plot(np.dot(dat.std[-2, iabscissa],[1,1]), array([ax[2]+1e-6, ax[3]-1e-6]), 'k-') # vertical separator
            # vertical separator
            plot(np.dot(dat.std[-2, iabscissa],[1,1]), array([np.min(dat.std[-2,5:]), np.max(dat.std[-2,5:])]), 'k-') 
            hold(True)
            # plot([dat.std[-1, iabscissa], ax[1]], [dat.std[-1,5:], yy[idx2]], 'k-') # line from last data point
            for i in xrange(len(idx)):
                # text(ax[1], yy[i], ' '+str(idx[i]))
                text(dat.std[-1, iabscissa], dat.std[-1, 5+i], ' '+str(i))
        semilogy(dat.std[:, iabscissa], dat.std[:,5:], '-') 
        grid(True)
        title('Standard Deviations in All Coordinates')
        # pylab.xticks(xticklocs)
        xlabel(xlab)
        draw()  # does not suffice
        if interactive_status:
            ion()  # turns interactive mode on (again)
            draw()  
        show()
        
        return self
 
    
    #____________________________________________________________
    #____________________________________________________________
    # 
    @staticmethod
    def plotdivers(dat, iabscissa, foffset):
        """helper function for `plot()` that plots all what is
        in the upper left subplot like fitness, sigma, etc. 
        
        Arguments
        ---------
            `iabscissa` in ``(0,1)`` 
                0==versus fevals, 1==versus iteration
            `foffset`
                offset to fitness for log-plot
                
         :See: `plot()` 
        
        """
        from  matplotlib.pylab import semilogy, hold, grid, \
                 axis, title, text
        fontsize = pylab.rcParams['font.size']
        
        hold(False)
    
        dfit = dat.f[:,5]-min(dat.f[:,5]) 
        dfit[dfit<1e-98] = np.NaN
            
        if dat.f.shape[1] > 7:
            # semilogy(dat.f[:, iabscissa], abs(dat.f[:,[6, 7, 10, 12]])+foffset,'-k')
            semilogy(dat.f[:, iabscissa], abs(dat.f[:,[6, 7]])+foffset,'-k')
            hold(True)
    
        # (larger indices): additional fitness data, for example constraints values
        if dat.f.shape[1] > 8:
            # dd = abs(dat.f[:,7:]) + 10*foffset
            # dd = np.where(dat.f[:,7:]==0, np.NaN, dd) # cannot be 
            semilogy(dat.f[:, iabscissa], np.abs(dat.f[:,8:]) + 10*foffset, 'm')
            hold(True)
    
        idx = np.where(dat.f[:,5]>1e-98)[0]  # positive values
        semilogy(dat.f[idx, iabscissa], dat.f[idx,5]+foffset, '.b')
        hold(True)
        grid(True)
        
        idx = np.where(dat.f[:,5] < -1e-98)  # negative values
        semilogy(dat.f[idx, iabscissa], abs(dat.f[idx,5])+foffset,'.r')
        
        semilogy(dat.f[:, iabscissa],abs(dat.f[:,5])+foffset,'-b')
        semilogy(dat.f[:, iabscissa], dfit, '-c')
        
        if 11 < 3:  # delta-fitness as points
            dfit = dat.f[1:, 5] - dat.f[:-1,5]  # should be negative usually
            semilogy(dat.f[1:,iabscissa],  # abs(fit(g) - fit(g-1))
                np.abs(dfit)+foffset, '.c')
            i = dfit > 0
            # print(np.sum(i) / float(len(dat.f[1:,iabscissa])))
            semilogy(dat.f[1:,iabscissa][i],  # abs(fit(g) - fit(g-1))
                np.abs(dfit[i])+foffset, '.r')
    
        # overall minimum
        i = np.argmin(dat.f[:,5])
        semilogy(dat.f[i, iabscissa]*np.ones(2), dat.f[i,5]*np.ones(2), 'rd') 
        # semilogy(dat.f[-1, iabscissa]*np.ones(2), dat.f[-1,4]*np.ones(2), 'rd') 
        
        # AR and sigma
        semilogy(dat.f[:, iabscissa], dat.f[:,3], '-r') # AR
        semilogy(dat.f[:, iabscissa], dat.f[:,2],'-g') # sigma
        semilogy(dat.std[:-1, iabscissa], np.vstack([map(max, dat.std[:-1,5:]), map(min, dat.std[:-1,5:])]).T,
                     '-m', linewidth=2)
        text(dat.std[-2, iabscissa], max(dat.std[-2, 5:]), 'max std', fontsize=fontsize)
        text(dat.std[-2, iabscissa], min(dat.std[-2, 5:]), 'min std', fontsize=fontsize)
        ax = array(axis())
        # ax[1] = max(minxend, ax[1])
        axis(ax)
        text(ax[0], ax[2], # 10**(log10(ax[2])+0.05*(log10(ax[3])-log10(ax[2]))),
             '.f_recent=' + repr(dat.f[-1,5]) )
        
        # title('abs(f) (blue), f-min(f) (cyan), Sigma (green), Axis Ratio (red)')
        title('blue:abs(f), cyan:f-min(f), green:sigma, red:axis ratio', fontsize=fontsize-1)
        # pylab.xticks(xticklocs)

            
    def downsampling(self, factor=10, first=3, switch=True):
        """
        rude downsampling of a `CMADataLogger` data file by `factor`, keeping 
        also the first `first` entries. This function is a stump and subject
        to future changes. 
        
        Arguments
        ---------
           - `factor` -- downsampling factor
           - `first` -- keep first `first` entries
           - `switch` -- switch the new logger name to oldname+'down'
        
        Details
        -------
        ``self.name_prefix+'down'`` files are written
        
        Example
        -------
        ::
        
            import cma
            cma.downsampling()  # takes outcmaes* files
            cma.plot('outcmaesdown')
           
        """
        newprefix = self.name_prefix + 'down'
        for name in CMADataLogger.names:
            f = open(newprefix+name+'.dat','w')
            iline = 0
            cwritten = 0
            for line in open(self.name_prefix+name+'.dat'):
                if iline < first or iline % factor == 0:
                    f.write(line)
                    cwritten += 1
                iline += 1
            f.close()
            print('%d' % (cwritten) + ' lines written in ' + newprefix+name+'.dat')
        if switch:
            self.name_prefix += 'down'
        return self
            
    #____________________________________________________________
    #____________________________________________________________
    # 
    def disp(self, idx=100):  # r_[0:5,1e2:1e9:1e2,-10:0]):
        """displays selected data from (files written by) the class `CMADataLogger`.
        
        Arguments
        ---------
           `idx` 
               indices corresponding to rows in the data file; 
               if idx is a scalar (int), the first two, then every idx-th, 
               and the last three rows are displayed. Too large index values are removed. 
               
        Example
        -------
        >>> import cma, numpy as np
        >>> res = cma.fmin(cma.fcts.elli, 7 * [0.1], 1, verb_disp=1e9)  # generate data
        >>> assert res[1] < 1e-9 
        >>> assert res[2] < 4400
        >>> l = cma.CMADataLogger()  # == res[-1], logger with default name, "points to" above data 
        >>> l.disp([0,-1])  # first and last
        >>> l.disp(20)  # some first/last and every 20-th line
        >>> l.disp(np.r_[0:999999:100, -1]) # every 100-th and last
        >>> l.disp(np.r_[0, -10:0]) # first and ten last
        >>> cma.disp(l.name_prefix, np.r_[0::100, -10:])  # the same as l.disp(...)
        
        Details
        -------
        The data line with the best f-value is displayed as last line. 

        :See: `disp()` 
        
        """
    
        filenameprefix=self.name_prefix
            
        def printdatarow(dat, iter):
            """print data of iteration i"""
            i = np.where(dat.f[:, 0] == iter)[0][0]
            j = np.where(dat.std[:, 0] == iter)[0][0]
            print('%5d' % (int(dat.f[i,0])) + ' %6d' % (int(dat.f[i,1])) + ' %.14e' % (dat.f[i,5]) + 
                  ' %5.1e' % (dat.f[i,3]) + 
                  ' %6.2e' % (max(dat.std[j,5:])) + ' %6.2e' % min(dat.std[j,5:]))
                  
        dat = CMADataLogger(filenameprefix).load()
        ndata = dat.f.shape[0]
    
        # map index to iteration number, is difficult if not all iteration numbers exist
        # idx = idx[np.where(map(lambda x: x in dat.f[:,0], idx))[0]] # TODO: takes pretty long
        # otherwise: 
        if idx is None:
            idx = 100
        if np.isscalar(idx):
            # idx = np.arange(0, ndata, idx)
            if idx:
                idx = np.r_[0, 1, idx:ndata-3:idx, -3:0]
            else:
                idx = np.r_[0, 1, -3:0]
        
        idx = array(idx)
        idx = idx[idx<=ndata]  # TODO: shouldn't this be "<"? 
        idx = idx[-idx<=ndata]
        iters = dat.f[idx, 0]
        idxbest = np.argmin(dat.f[:,5])
        iterbest = dat.f[idxbest, 0]
        
        heading = 'Iterat Nfevals  function value    axis ratio maxstd   minstd'
        print(heading)
        for i in iters:
            printdatarow(dat, i)
        print(heading)
        printdatarow(dat, iterbest)

# end class CMADataLogger    

#____________________________________________________________
#____________________________________________________________
# 
def fmin(func, x0, sigma0=None, args=()
    # the follow string arguments are evaluated, besides the verb_filenameprefix
    , CMA_active='False  # exponential negative update, conducted after the original update'
    , CMA_activefac='1  # learning rate multiplier for active update'
    , CMA_cmean='1  # learning rate for the mean value'
    , CMA_diagonal='0*100*N/sqrt(popsize)  # nb of iterations with diagonal covariance matrix, True for always' # TODO 4/ccov_separable?
    , CMA_eigenmethod='0  # 0=numpy-s eigh, -1=pygsl, otherwise cma.eig (slower)'
    , CMA_elitist='False # elitism impairs global search performance'
    , CMA_mirrors='0  # values <0.5 are interpreted as fraction, values >1 as numbers (rounded), otherwise about 0.16 is used'
    , CMA_mu='None  # parents selection parameter, default is popsize // 2'
    , CMA_on='True  # False or 0 for no adaptation of the covariance matrix'
    , CMA_rankmu='True  # False or 0 for omitting rank-mu update of covariance matrix'
    , CMA_rankmualpha='0.3  # factor of rank-mu update if mu=1, subject to removal, default might change to 0.0'
    , CMA_dampfac='1  #v positive multiplier for step-size damping, 0.3 is close to optimal on the sphere'
    , CMA_teststds='None  # factors for non-isotropic initial distr. mainly for test purpose, see scaling_...'
    , bounds='[None, None]  # lower (=bounds[0]) and upper domain boundaries, each a scalar or a list/vector'
    , check_points='None  # when repairing or injecting solutions, they should be checked (index-list or True)' 
    , eval_initial_x='False  # '
    , fixed_variables='None  # dictionary with index-value pairs like {0:1.1, 2:0.1} that are not optimized'
    , ftarget='-inf  #v target function value, minimization'
    , incpopsize='2  # in fmin(): multiplier for increasing popsize before each restart'
    , maxfevals='inf  #v maximum number of function evaluations'
    , maxiter='long(100 + 50 * (N+3)**2 / sqrt(popsize))  #v maximum number of iterations'
    , mindx='0  #v minimal std in any direction, cave interference with tol*'
    , minstd='0  #v minimal std in any coordinate direction, cave interference with tol*'
    , noise_handling='False  # maximal number of evaluations for noise treatment, only fmin'
    , noise_reevals=' 1.5 + popsize/20  # number of solution to be reevaluated for noise measurement, only fmin'
    , noise_eps='1e-7  # perturbation factor for noise handling reevaluations, only fmin'
    , popsize='4+int(3*log(N))  # population size, AKA lambda, number of new solution per iteration'
    , randn='np.random.standard_normal  #v randn((lam, N)) must return an np.array of shape (lam, N)'
    , restarts='0  # in fmin(): number of restarts' 
    , scaling_of_variables='None  # scale for each variable, sigma0 is interpreted w.r.t. this scale, in that effective_sigma0 = sigma0*scaling. Internally the variables are divided by scaling_of_variables and sigma is unchanged, default is ones(N)'
    , seed='None  # random number seed'
    , termination_callback='None  #v in fmin(): a function returning True for termination, called after each iteration step and could be abused for side effects'
    , tolfacupx='1e3  #v termination when step-size increases by tolfacupx (diverges). That is, the initial step-size was chosen far too small and better solutions were found far away from the initial solution x0'
    , tolupsigma='1e20  #v sigma/sigma0 > tolupsigma * max(sqrt(eivenvals(C))) indicates "creeping behavior" with usually minor improvements'
    , tolfun='1e-11  #v termination criterion: tolerance in function value, quite useful'
    , tolfunhist='1e-12  #v termination criterion: tolerance in function value history'
    , tolstagnation='int(100 * N**1.5 / popsize)  #v termination if no improvement over tolstagnation iterations'
    , tolx='1e-11  #v termination criterion: tolerance in x-changes'
    , transformation='None  # [t0, t1] are two mappings, t0 transforms solutions from CMA-representation to f-representation, t1 is the back transformation, see class GenoPheno'
    , typical_x='None  # used with scaling_of_variables'
    , updatecovwait='None  #v number of iterations without distribution update, name is subject to future changes' # TODO: rename: iterwaitupdatedistribution?
    , verb_append='0  # initial evaluation counter, if append, do not overwrite output files' 
    , verb_disp='100  #v verbosity: display console output every verb_disp iteration'
    , verb_filenameprefix='outcmaes  # output filenames prefix'
    , verb_log='1  #v verbosity: write data to files every verb_log iteration, writing can be time critical on fast to evaluate functions'
    , verb_plot='0  #v in fmin(): plot() is called every verb_plot iteration'
    , verb_time='True  #v output timings on console'
    , vv='0  #? versatile variable for hacking purposes, value found in self.opts[\'vv\']'
     ): 
    """functional interface to the stochastic optimizer CMA-ES 
    for non-convex function minimization. 

    Calling Sequences
    =================
        ``fmin([],[])``
            returns all optional arguments, that is,
            all keyword arguments to fmin with their default values 
            in a dictionary.
        ``fmin(func, x0, sigma0)``
            minimizes `func` starting at `x0` and with standard deviation 
            `sigma0` (step-size)
        ``fmin(func, x0, sigma0, ftarget=1e-5)``
            minimizes `func` up to target function value 1e-5
        ``fmin(func, x0, sigma0, args=('f',), **options)`` 
            minimizes `func` called with an additional argument ``'f'``. 
            `options` is a dictionary with additional keyword arguments, e.g.
            delivered by `Options()`. 
        ``fmin(func, x0, sigma0, **{'ftarget':1e-5, 'popsize':40})`` 
            the same as ``fmin(func, x0, sigma0, ftarget=1e-5, popsize=40)``
        ``fmin(func, esobj, **{'maxfevals': 1e5})``
            uses the `CMAEvolutionStrategy` object instance `esobj` to optimize 
            `func`, similar to `CMAEvolutionStrategy.optimize()`. 

    Arguments
    ========= 
        `func`
            function to be minimized. Called as
            ``func(x,*args)``. `x` is a one-dimensional `numpy.ndarray`. `func`
            can return `numpy.NaN`,
            which is interpreted as outright rejection of solution `x`
            and invokes an immediate resampling and (re-)evaluation
            of a new solution not counting as function evaluation. 
        `x0` 
            list or `numpy.ndarray`, initial guess of minimum solution
            or `cma.CMAEvolutionStrategy` object instance. In this case
            `sigma0` can be omitted.
        `sigma0`
            scalar, initial standard deviation in each coordinate.
            `sigma0` should be about 1/4 of the search domain width where the
            optimum is to be expected. The variables in `func` should be
            scaled such that they presumably have similar sensitivity.
            See also option `scaling_of_variables`. 

    Keyword Arguments
    =================
    All arguments besides `args` and `verb_filenameprefix` are evaluated 
    if they are of type `str`, see class `Options` for details.  
    ::

        args=() -- additional arguments for func, not in `cma.Options()`
        CMA_active='False  # exponential negative update, conducted after the original
                update'
        CMA_activefac='1  # learning rate multiplier for active update'
        CMA_cmean='1  # learning rate for the mean value'
        CMA_dampfac='1  #v positive multiplier for step-size damping, 0.3 is close to
                optimal on the sphere'
        CMA_diagonal='0*100*N/sqrt(popsize)  # nb of iterations with diagonal
                covariance matrix, True for always'
        CMA_eigenmethod='0  # 0=numpy-s eigh, -1=pygsl, otherwise cma.eig (slower)'
        CMA_elitist='False # elitism impairs global search performance'
        CMA_mirrors='0  # values <0.5 are interpreted as fraction, values >1 as numbers
                (rounded), otherwise about 0.16 is used'
        CMA_mu='None  # parents selection parameter, default is popsize // 2'
        CMA_on='True  # False or 0 for no adaptation of the covariance matrix'
        CMA_rankmu='True  # False or 0 for omitting rank-mu update of covariance
                matrix'
        CMA_rankmualpha='0.3  # factor of rank-mu update if mu=1, subject to removal,
                default might change to 0.0'
        CMA_teststds='None  # factors for non-isotropic initial distr. mainly for test
                purpose, see scaling_...'
        bounds='[None, None]  # lower (=bounds[0]) and upper domain boundaries, each a
                scalar or a list/vector'
        check_points='None  # when repairing or injecting solutions, they should be checked 
                (index-list or True)' 
        eval_initial_x='False  # '
        fixed_variables='None  # dictionary with index-value pairs like {0:1.1, 2:0.1}
                that are not optimized'
        ftarget='-inf  #v target function value, minimization'
        incpopsize='2  # in fmin(): multiplier for increasing popsize before each
                restart'
        maxfevals='inf  #v maximum number of function evaluations'
        maxiter='long(1e3*N**2/sqrt(popsize))  #v maximum number of iterations'
        mindx='0  #v minimal std in any direction, cave interference with tol*'
        minstd='0  #v minimal std in any coordinate direction, cave interference with
                tol*'
        noise_eps='1e-7  # perturbation factor for noise handling reevaluations, only
                fmin'
        noise_handling='False  # maximal number of evaluations for noise treatment,
                only fmin'
        noise_reevals=' 1.5 + popsize/20  # number of solution to be reevaluated for
                noise measurement, only fmin'
        popsize='4+int(3*log(N))  # population size, AKA lambda, number of new solution
                per iteration'
        randn='np.random.standard_normal  #v randn((lam, N)) must return an np.array of
                shape (lam, N)'
        restarts='0  # in fmin(): number of restarts'
        scaling_of_variables='None  # scale for each variable, sigma0 is interpreted
                w.r.t. this scale, in that effective_sigma0 = sigma0*scaling.
                Internally the variables are divided by scaling_of_variables and sigma
                is unchanged, default is ones(N)'
        seed='None  # random number seed'
        termination_callback='None  #v in fmin(): a function returning True for
                termination, called after each iteration step and could be abused for
                side effects'
        tolfacupx='1e3  #v termination when step-size increases by tolfacupx
                (diverges). That is, the initial step-size was chosen far too small and
                better solutions were found far away from the initial solution x0'
        tolupsigma='1e20  #v sigma/sigma0 > tolupsigma * max(sqrt(eivenvals(C))) 
                indicates "creeping behavior" with usually minor improvements'
        tolfun='1e-11  #v termination criterion: tolerance in function value, quite
                useful'
        tolfunhist='1e-12  #v termination criterion: tolerance in function value
                history'
        tolstagnation='int(100 * N**1.5 / popsize)  #v termination if no improvement
                over tolstagnation iterations'
        tolx='1e-11  #v termination criterion: tolerance in x-changes'
        transformation='None  # [t0, t1] are two mappings, t0 transforms solutions from
                CMA-representation to f-representation, t1 is the back transformation,
                see class GenoPheno'
        typical_x='None  # used with scaling_of_variables'
        updatecovwait='None  #v number of iterations without distribution update, name
                is subject to future changes'
        verb_append='0  # initial evaluation counter, if append, do not overwrite
                output files'
        verb_disp='100  #v verbosity: display console output every verb_disp iteration'
        verb_filenameprefix='outcmaes  # output filenames prefix'
        verb_log='1  #v verbosity: write data to files every verb_log iteration,
                writing can be time critical on fast to evaluate functions'
        verb_plot='0  #v in fmin(): plot() is called every verb_plot iteration'
        verb_time='True  #v output timings on console'
        vv='0  #? versatile variable for hacking purposes, value found in
                self.opts['vv']'

    Subsets of options can be displayed, for example like ``cma.Options('tol')``, 
    see also class `Options`.  
        
    Return
    ====== 
    Similar to `OOOptimizer.optimize()` and/or `CMAEvolutionStrategy.optimize()`, return the 
    list provided by `CMAEvolutionStrategy.result()` appended with an `OOOptimizer` and an 
    `OptimDataLogger`::
    
        res = optim.result() + (optim.stop(), optim, logger)
    
    where
        - ``res[0]`` (``xopt``) -- best evaluated solution 
        - ``res[1]`` (``fopt``) -- respective function value
        - ``res[2]`` (``evalsopt``) -- respective number of function evaluations
        - ``res[3]`` (``evals``) -- number of overall conducted objective function evaluations
        - ``res[4]`` (``iterations``) -- number of overall conducted iterations
        - ``res[5]`` (``xmean``) -- mean of the final sample distribution
        - ``res[6]`` (``stds``) -- effective stds of the final sample distribution
        - ``res[-3]`` (``stop``) -- termination condition(s) in a dictionary
        - ``res[-2]`` (``cmaes``) -- class `CMAEvolutionStrategy` instance
        - ``res[-1]`` (``logger``) -- class `CMADataLogger` instance
                                                                        
    Details
    =======
    This function is an interface to the class `CMAEvolutionStrategy`. The
    class can be used when full control over the iteration loop of the
    optimizer is desired. 
    
    The noise handling follows closely [Hansen et al 2009, A Method for Handling 
    Uncertainty in Evolutionary Optimization...] in the measurement part, but the 
    implemented treatment is slightly different: for ``noiseS > 0``, ``evaluations`` 
    (time) and sigma are increased by ``alpha``. For ``noiseS < 0``, ``evaluations`` 
    (time) is decreased by ``alpha**(1/4)``. The option ``noise_handling`` switches 
    the uncertainty handling on/off and gives the maximal number of evaluations for 
    a single fitness computation. If ``noise_handling`` is a list, the smallest 
    element defines the minimal number and if the list has three elements, the 
    median value is the start value for ``evaluations``. 
    See also class `NoiseHandler`.  
    
    Examples
    ========
    The following example calls `fmin` optimizing the Rosenbrock function 
    in 10-D with initial solution 0.1 and initial step-size 0.5. The
    options are specified for the usage with the `doctest` module.
    
    >>> import cma
    >>> # cma.Options()  # returns all possible options
    >>> options = {'CMA_diagonal':10, 'seed':1234, 'verb_time':0}
    >>>
    >>> res = cma.fmin(cma.fcts.rosen, [0.1] * 10, 0.5, **options)
    (5_w,10)-CMA-ES (mu_w=3.2,w_1=45%) in dimension 10 (seed=1234)
       Covariance matrix is diagonal for 10 iterations (1/ccov=29.0)
    Iterat #Fevals   function value     axis ratio  sigma   minstd maxstd min:sec
        1      10 1.264232686260072e+02 1.1e+00 4.40e-01  4e-01  4e-01 
        2      20 1.023929748193649e+02 1.1e+00 4.00e-01  4e-01  4e-01 
        3      30 1.214724267489674e+02 1.2e+00 3.70e-01  3e-01  4e-01 
      100    1000 6.366683525319511e+00 6.2e+00 2.49e-02  9e-03  3e-02 
      200    2000 3.347312410388666e+00 1.2e+01 4.52e-02  8e-03  4e-02 
      300    3000 1.027509686232270e+00 1.3e+01 2.85e-02  5e-03  2e-02 
      400    4000 1.279649321170636e-01 2.3e+01 3.53e-02  3e-03  3e-02 
      500    5000 4.302636076186532e-04 4.6e+01 4.78e-03  3e-04  5e-03 
      600    6000 6.943669235595049e-11 5.1e+01 5.41e-06  1e-07  4e-06 
      650    6500 5.557961334063003e-14 5.4e+01 1.88e-07  4e-09  1e-07 
    termination on tolfun : 1e-11
    final/bestever f-value = 5.55796133406e-14 2.62435631419e-14
    mean solution:  [ 1.          1.00000001  1.          1.          
        1.          1.00000001  1.00000002  1.00000003 ...]
    std deviation: [ 3.9193387e-09  3.7792732e-09  4.0062285e-09  4.6605925e-09
        5.4966188e-09   7.4377745e-09   1.3797207e-08   2.6020765e-08 ...]
    >>> 
    >>> print('best solutions fitness = %f' % (res[1]))
    best solutions fitness = 2.62435631419e-14
    >>> assert res[1] < 1e-12
    
    The method ::
        
        cma.plot();  
        
    (based on `matplotlib.pylab`) produces a plot of the run and, if necessary::
        
        cma.show()  
    
    shows the plot in a window. To continue you might need to 
    close the pop-up window. This behavior seems to disappear in 
    subsequent calls of `cma.plot()` and is avoided by using
    `ipython` with `-pylab` option. Finally ::
    
        cma.savefig('myfirstrun')  # savefig from matplotlib.pylab
    
    will save the figure in a png.

    :See: `CMAEvolutionStrategy`, `OOOptimizer.optimize(), `plot()`, `Options`, `scipy.optimize.fmin()` 
    
    """ # style guides say there should be the above empty line
    try: # pass on KeyboardInterrupt
        opts = locals()  # collect all local variables (i.e. arguments) in a dictionary
        del opts['func'] # remove those without a default value 
        del opts['args']
        del opts['x0']      # is not optional, no default available
        del opts['sigma0']  # is not optional for the constructor CMAEvolutionStrategy
        if not func:  # return available options in a dictionary
            return Options(opts, True)  # these opts are by definition valid

        irun = 0
        best = BestSolution()
        # TODO: this is very ugly:
        incpopsize = Options({'incpopsize':incpopsize}).eval('incpopsize') 
        restarts = Options({'restarts':restarts}).eval('restarts') 
        while 1:
            # recover from a CMA object
            if irun == 0 and isinstance(x0, CMAEvolutionStrategy):
                es = x0
                x0 = es.inputargs['x0']  # for the next restarts
                if sigma0 is None or not np.isscalar(array(sigma0)): 
                    sigma0 = es.inputargs['sigma0']  # for the next restarts
                # ignore further input args and keep original options
            else:  # default case
                es = CMAEvolutionStrategy(x0, sigma0, opts) 
                if opts['eval_initial_x']:
                    x = es.gp.pheno(es.mean)
                    es.best.update([x], [func(x, *args)], 1)
                    es.countevals += 1
                
            opts = es.opts  # processed options, unambiguous
                
            append = opts['verb_append'] or es.countiter > 0 or irun > 0
            logger = CMADataLogger(opts['verb_filenameprefix'], opts['verb_log'])
            logger.register(es, append).add()  # initial values, not fitness values
            
            # if es.countiter == 0 and es.opts['verb_log'] > 0 and not es.opts['verb_append']:
            #    logger = CMADataLogger(es.opts['verb_filenameprefix']).register(es)
            #    logger.add()
            # es.writeOutput()  # initial values for sigma etc
                
            noisehandler = NoiseHandler(es.N, opts['noise_handling'], np.median, opts['noise_reevals'], opts['noise_eps'])
            while not es.stop():
                X, fit = es.ask_and_eval(func, args, evaluations=noisehandler.evaluations, 
                                         aggregation=np.median) # treats NaN with resampling
                # TODO: check args and in case use args=(noisehandler.evaluations, )
                
                if 11 < 3 and opts['vv']:  # inject a solution
                    # use option check_point = [0]
                    if 0 * np.random.randn() >= 0:
                        X[0] = 0 + opts['vv'] * es.sigma**0 * np.random.randn(es.N)
                        fit[0] = func(X[0], *args)
                        # print fit[0]

                es.tell(X, fit)  # prepare for next iteration
                if opts['noise_handling']:
                    es.sigma *= noisehandler(X, fit, func, es.ask, args)
                    es.countevals += noisehandler.evaluations_just_done  # TODO: this is a hack, not important though

                es.disp()
                logger.add(more_data=[noisehandler.evaluations, 10**noisehandler.noiseS] if opts['noise_handling'] else [], 
                           modulo=1 if es.stop() and logger.modulo else None)
                if opts['verb_log'] and opts['verb_plot'] and \
                    (es.countiter % max(opts['verb_plot'], opts['verb_log']) == 0 or es.stop()):
                    logger.plot(324, fontsize=10)

                if termination_callback and termination_callback != str(termination_callback): 
                    es.callbackstop = termination_callback(es)

            # end while not es.stop

            mean_pheno = es.gp.pheno(es.mean, bounds=es.gp.bounds)
            fmean = func(mean_pheno, *args)
            es.countevals += 1
            
            es.best.update([mean_pheno], [fmean], es.countevals)
            best.update(es.best)  # in restarted case
            
            # final message
            if opts['verb_disp']:
                for k, v in es.stop().items():
                    print('termination on %s=%s (%s)' % (k, str(v), time.asctime()))
                print('final/bestever f-value = %e %e' % (es.out['recent_f'], best.f))
                if es.N < 9:
                    print('mean solution: ' + str(es.gp.pheno(es.mean)))
                    print('std deviation: ' + str(es.sigma * sqrt(es.dC) * es.gp.scales))
                else:
                    print('mean solution: %s ...]' % (str(es.gp.pheno(es.mean)[:8])[:-1])) 
                    print('std deviations: %s ...]' % (str((es.sigma * sqrt(es.dC) * es.gp.scales)[:8])[:-1]))

            irun += 1
            if irun > restarts or 'ftarget' in es.stopdict or 'maxfunevals' in es.stopdict:
                break
            opts['verb_append'] = es.countevals
            opts['popsize'] = incpopsize * es.sp.popsize # TODO: use rather options? 
            opts['seed'] += 1

        # while irun

        es.out['best'] = best  # TODO: this is a rather suboptimal type for inspection in the shell
        if 1 < 3: 
            return es.result() + (es.stop(), es, logger)
        
        else: # previously: to be removed
            return (best.x.copy(), best.f, es.countevals, 
                    dict((('stopdict', CMAStopDict(es.stopdict)) 
                          ,('mean', es.gp.pheno(es.mean)) 
                          ,('std', es.sigma * sqrt(es.dC) * es.gp.scales)
                          ,('out', es.out)
                          ,('opts', es.opts)  # last state of options 
                          ,('cma', es)
                          ,('inputargs', es.inputargs) 
                          ))
                   )
        # TODO refine output, can #args be flexible?
        # is this well usable as it is now?
    except KeyboardInterrupt:  # Exception, e:
        if opts['verb_disp'] > 0:
            print(' in/outcomment ``raise`` in last line of cma.fmin to prevent/restore KeyboardInterrupt exception')
        raise  # cave: swallowing this exception can silently mess up experiments, if ctrl-C is hit
        
# end fmin()


def plot(name=None, fig=None, abscissa=1, iteridx=None, plot_mean=True,  # TODO: plot_mean default should be False
    foffset=1e-19, x_opt=None, fontsize=None):
    """
    plot data from files written by a `CMADataLogger`, 
    the call ``cma.plot(name, **argsdict)`` is a shortcut for 
    ``cma.CMADataLogger(name).plot(**argsdict)``

    Arguments
    ---------
        `name` 
            name of the logger, filename prefix, None evaluates to 
            the default 'outcmaes'
        `fig`
            filename or figure number, or both as a tuple (any order)
        `abscissa`
            0==plot versus iteration count,
            1==plot versus function evaluation number
        `iteridx`
            iteration indices to plot
    
    Return `None`
          
    Examples
    --------
    ::
    
       cma.plot();  # the optimization might be still 
                    # running in a different shell
       cma.show()  # to continue you might need to close the pop-up window 
                   # once and call cma.plot() again. 
                   # This behavior seems to disappear in subsequent 
                   # calls of cma.plot(). Also using ipython with -pylab
                   # option might help. 
       cma.savefig('fig325.png')
       cma.close()

       cdl = cma.CMADataLogger().downsampling().plot()
       
    Details
    -------
    Data from codes in other languages (C, Java, Matlab, Scilab) have the same 
    format and can be plotted just the same. 
       
    :See: `CMADataLogger`, `CMADataLogger.plot()`

    """
    CMADataLogger(name).plot(fig, abscissa, iteridx, plot_mean, foffset, x_opt, fontsize)

try:
    show = pylab.show
    savefig = pylab.savefig   # we would like to be able to use cma.savefig() etc
    closefig = pylab.close
except:
    def show():
        pass     

def disp(name=None, idx=None):
    """displays selected data from (files written by) the class `CMADataLogger`. 

    The call ``cma.disp(name, idx)`` is a shortcut for ``cma.CMADataLogger(name).disp(idx)``.
    
    Arguments
    ---------
        `name`
            name of the logger, filename prefix, `None` evaluates to 
            the default ``'outcmaes'``
        `idx`   
            indices corresponding to rows in the data file; by
            default the first five, then every 100-th, and the last
            10 rows. Too large index values are removed. 
    
    Examples
    --------
    ::

       import cma, numpy
       # assume some data are available from previous runs
       cma.disp(None,numpy.r_[0,-1])  # first and last
       cma.disp(None,numpy.r_[0:1e9:100,-1]) # every 100-th and last
       cma.disp(idx=numpy.r_[0,-10:0]) # first and ten last
       cma.disp(idx=numpy.r_[0:1e9:1e3,-10:0]) 
    
    :See: `CMADataLogger.disp()`
    
    """
    return CMADataLogger(name if name else 'outcmaes'
                         ).disp(idx)

#____________________________________________________________

def _fileToMatrix(file_name):
    """rudimentary method to read in data from a file"""
    # TODO: np.loadtxt() might be an alternative
    #     try: 
    if 1 < 3: 
        lres = []
        for line in open(file_name, 'r').readlines():
            if len(line) > 0 and line[0] not in ('%', '#'):
                lres.append(map(float, line.split()))
        res = lres
    else:
        fil = open(file_name, 'r')
        fil.readline() # rudimentary, assume one comment line
        lineToRow = lambda line: map(float, line.split())
        res = map(lineToRow, fil.readlines()) 
        fil.close()  # close file could be omitted, reference counting should do during garbage collection, but...

    while res != [] and res[0] == []:  # remove further leading empty lines
        del res[0] 
    return res
    #     except:
    print('could not read file ' + file_name)

#____________________________________________________________
#____________________________________________________________

class NoiseHandler(object):
    """Noise handling according to [Hansen et al 2009, A Method for Handling 
    Uncertainty in Evolutionary Optimization...] 
    
    The interface of this class is yet versatile and subject to changes. 
    
    The attribute ``evaluations`` serves to control the noise via number of 
    evaluations, for example with `ask_and_eval()`, see also parameter 
    ``maxevals`` and compare the example. 
    
    Example
    -------
    >>> import cma, numpy as np
    >>> func = cma.Fcts.noisysphere
    >>> es = cma.CMAEvolutionStrategy(np.ones(10), 1)
    >>> logger = cma.CMADataLogger().register(es)
    >>> nh = cma.NoiseHandler(es.N, maxevals=[1, 30])
    >>> while not es.stop():
    ...     X, fit = es.ask_and_eval(func, evaluations=nh.evaluations) 
    ...     es.tell(X, fit)  # prepare for next iteration
    ...     es.sigma *= nh(X, fit, func, es.ask)  # see method __call__
    ...     es.countevals += nh.evaluations_just_done  # this is a hack, not important though
    ...     logger.add(more_data = [nh.evaluations, nh.noiseS])  # add a data point 
    ...     es.disp()
    ...     # nh.maxevals = ...  it might be useful to start with smaller values and then increase
    >>> print(es.stop())
    >>> print(es.result()[-2])  # take mean value, the best solution is totally off
    >>> assert sum(es.result()[-2]**2) < 1e-9
    >>> print(X[np.argmin(fit)])  # not bad, but probably worse than the mean
    >>> logger.plot()
    
    The noise options of `fmin()` control a `NoiseHandler` instance similar to this
    example. The command ``cma.Options('noise')`` lists in effect the parameters of 
    `__init__` apart from ``aggregate``.  
    
    Details
    -------
    The parameters reevals, theta, c_s, and alpha_t are set differently
    than in the original publication, see method `__init__()`. For a
    very small population size, say popsize <= 5, the measurement
    technique based on rank changes is likely to fail.
    
    Missing Features
    ----------------
    In case no noise is found, ``self.lam_reeval`` should be adaptive
    and get at least as low as 1 (however the possible savings from this
    are rather limited). Another option might be to decide during the
    first call by a quantitative analysis of fitness values whether
    ``lam_reeval`` is set to zero. More generally, an automatic noise
    mode detection might also set the covariance matrix learning rates
    to smaller values.
    
    :See: `fmin()`, `ask_and_eval()`
    
    """
    def __init__(self, N, maxevals=10, aggregate=np.median, reevals=None, epsilon=1e-7):
        """parameters are

            `maxevals`
                maximal value for ``self.evaluations``, where
                ``self.evaluations`` function calls are aggregated for
                noise treatment. With ``maxevals == 0`` the noise
                handler is (temporarily) "switched off". If `maxevals`
                is a list, min value and (for >2 elements) median are
                used to define minimal and initial value of
                ``self.evaluations``. Choosing ``maxevals > 1`` is only
                reasonable, if also the original ``fit`` values (that
                are passed to `__call__`) are computed by aggregation of
                ``self.evaluations`` values (otherwise the values are
                not comparable), as it is done within `fmin()`.
            `aggregate`
                function to aggregate single f-values to a 'fitness', e.g. 
                ``np.median``. 
            `reevals`
                number of solutions to be reevaluated for noise measurement,
                can be a float, by default set to ``1.5 + popsize/20``, 
                zero switches noise handling off. 
            `epsilon`
                multiplier for perturbation of the reevaluated solutions
                 
            :See: `fmin()`, `Options`, `CMAEvolutionStrategy.ask_and_eval()`
            
        """
        self.lam_reeval = reevals  # 2 + popsize/20, see method indices(), originally 2 + popsize/10
        self.epsilon = epsilon
        self.theta = 0.5  # originally 0.2
        self.cum = 0.3  # originally 1, 0.3 allows one disagreement of current point with resulting noiseS
        self.alphasigma = 1 + 2 / (N+10)
        self.alphaevals = 1 + 2 / (N+10)  # originally 1.5
        self.alphaevalsdown = self.alphaevals**-0.25  # originally 1/1.5
        self.evaluations = 1  # to aggregate for a single f-evaluation
        self.minevals = 1
        self.maxevals = np.max(maxevals)
        if hasattr(maxevals, '__contains__'):  # i.e. can deal with ``in``
            if len(maxevals) > 1:
                self.minevals = min(maxevals)
                self.evaluations = self.minevals
            if len(maxevals) > 2:
                self.evaluations = np.median(maxevals)
        self.f_aggregate = aggregate
        self.evaluations_just_done = 0  # actually conducted evals, only for documentation
        self.noiseS = 0
    
    def __call__(self, X, fit, func, ask=None, args=()):
        """proceed with noise measurement, set anew attributes ``evaluations``
        (proposed number of evaluations to "treat" noise) and ``evaluations_just_done`` 
        and return a factor for increasing sigma.
        
        Parameters
        ----------
            `X`
                a list/sequence/vector of solutions
            `fit`
                the respective list of function values
            `func`
                the objective function, ``fit[i]`` corresponds to ``func(X[i], *args)``
            `ask`
                a method to generate a new, slightly disturbed solution. The argument
                is mandatory if ``epsilon`` is not zero, see `__init__()`.
            `args`
                optional additional arguments to `func`
                
        Details
        -------
        Calls the methods ``reeval()``, ``update_measure()`` and ``treat()`` in this order. 
        ``self.evaluations`` is adapted within the method `treat()`. 
                
        """
        self.evaluations_just_done = 0
        if not self.maxevals or self.lam_reeval == 0:
            return 1.0
        res = self.reeval(X, fit, func, ask, args)
        if not len(res):
            return 1.0
        self.update_measure()
        return self.treat()
    
    def get_evaluations(self):
        """return ``self.evaluations``, the number of evalutions to get a single fitness measurement"""
        return self.evaluations
    
    def treat(self):
        """adapt self.evaluations and return ``sigma_fac in (1.0, self.alphasigma)`` 
        
        """
        if self.noiseS > 0:
            self.evaluations = min((self.evaluations * self.alphaevals, self.maxevals))
            return self.alphasigma
        else:
            self.evaluations = max((self.evaluations * self.alphaevalsdown, self.minevals))
            return 1.0

    def reeval(self, X, fit, func, ask, args=()):
        """store two fitness lists, `fit` and ``fitre`` reevaluating some 
        solutions in `X`. 
        ``self.evaluations`` evaluations are done for each reevaluated 
        fitness value. 
        See `__call__()`, where `reeval()` is called. 
        
        """
        self.fit = list(fit)
        self.fitre = list(fit)
        self.idx = self.indices(fit)
        if not len(self.idx):
            return self.idx
        evals = int(self.evaluations) if self.f_aggregate else 1
        fagg = np.median if self.f_aggregate is None else self.f_aggregate
        for i in self.idx:
            if self.epsilon:
                self.fitre[i] = fagg([func(ask(1, X[i], self.epsilon)[0], *args) for _k in xrange(evals)])
            else:
                self.fitre[i] = fagg([func(X[i], *args) for _k in xrange(evals)])
        self.evaluations_just_done = evals * len(self.idx)
        return self.fit, self.fitre, self.idx

    def update_measure(self):
        """updated noise level measure using two fitness lists ``self.fit`` and 
        ``self.fitre``, return ``self.noiseS, all_individual_measures``. 
        
        Assumes that `self.idx` contains the indices where the fitness
        lists differ
        
        """
        lam = len(self.fit)
        idx = np.argsort(self.fit + self.fitre)
        ranks = np.argsort(idx).reshape((2, lam))
        rankDelta = ranks[0] - ranks[1] - np.sign(ranks[0] - ranks[1]);

        # compute rank change limits using both ranks[0] and ranks[1]
        r = np.arange(1, 2 * lam)  # 2 * lam - 2 elements
        limits = [0.5 * (Mh.prctile(np.abs(r - (ranks[0,i] + 1 - (ranks[0,i] > ranks[1,i]))), 
                                      self.theta*50) +
                         Mh.prctile(np.abs(r - (ranks[1,i] + 1 - (ranks[1,i] > ranks[0,i]))), 
                                      self.theta*50))
                    for i in self.idx]
        # compute measurement
        #                               max: 1 rankchange in 2*lambda is always fine
        s = np.abs(rankDelta[self.idx]) - Mh.amax(limits, 1)  # lives roughly in 0..2*lambda
        self.noiseS += self.cum * (np.mean(s) - self.noiseS)
        return self.noiseS, s
        
    def indices(self, fit):
        """return the set of indices to be reevaluted for noise measurement,
        taking the ``lam_reeval`` best from the first ``2 * lam_reeval + 2``
        values. 
        
        Given the first values are the earliest, this is a useful policy also 
        with a time changing objective. 
        
        """
        lam = self.lam_reeval if self.lam_reeval else 2 + len(fit) / 20
        reev = int(lam) + ((lam % 1) > np.random.rand())
        return np.argsort(array(fit, copy=False)[:2 * (reev + 1)])[:reev]
    
#____________________________________________________________
#____________________________________________________________

class Sections(object):
    """plot sections through an objective function. A first 
    rational thing to do, when facing an (expensive) application. 
    By default 6 points in each coordinate are evaluated. 
    This class is still experimental. 
    
    Examples
    --------
    
    >>> import cma, numpy as np
    >>> s = cma.Sections(cma.Fcts.rosen, np.zeros(3)).do(plot=False)
    >>> s.do(plot=False)  # evaluate the same points again, i.e. check for noise
    >>> try:
    ...     s.plot()
    ... except: 
    ...     print('plotting failed: pylab package is missing?')
    
    Details
    -------
    Data are saved after each function call during `do()`. The filename is attribute 
    ``name`` and by default ``str(func)``, see `__init__()`. 
    
    A random (orthogonal) basis can be generated with ``cma.Rotation()(np.eye(3))``. 
    
    The default name is unique in the function name, but it should be unique in all
    parameters of `__init__()` but `plot_cmd` and `load`. 

    ``self.res`` is a dictionary with an entry for each "coordinate" ``i`` and with an 
    entry ``'x'``, the middle point. Each entry ``i`` is again a dictionary with keys 
    being different dx values and the value being a sequence of f-values.
    For example ``self.res[2][0.1] == [0.01, 0.01]``, which is generated using the 
    difference vector ``self.basis[2]`` like 
    ``self.res[2][dx] += func(self.res['x'] + dx * self.basis[2])``.  
                
    :See: `__init__()`
    
    """
    def __init__(self, func, x, args=(), basis=None, name=None, 
                 plot_cmd=pylab.plot if pylab else None, load=True):
        """        
        Parameters
        ----------
            `func`
                objective function
            `x` 
                point in search space, middle point of the sections
            `args`
                arguments passed to `func`
            `basis`
                evaluated points are ``func(x + locations[j] * basis[i]) for i in len(basis) for j in len(locations)``, 
                see `do()` 
            `name`
                filename where to save the result
            `plot_cmd`
                command used to plot the data, typically matplotlib pylabs `plot` or `semilogy`
            `load`
                load previous data from file ``str(func) + '.pkl'``
                
        """
        self.func = func
        self.args = args
        self.x = x
        self.name = name if name else str(func)
        self.plot_cmd = plot_cmd  # or semilogy
        self.basis = np.eye(len(x)) if basis is None else basis 
        
        try:
            self.load()
            if any(self.res['x'] != x):
                self.res = {}
                self.res['x'] = x  # TODO: res['x'] does not look perfect
            else:
                print(self.name + ' loaded')
        except:
            self.res = {}
            self.res['x'] = x
        
    def do(self, repetitions=1, locations=np.arange(-0.5, 0.6, 0.2), plot=True):
        """generates, plots and saves function values ``func(y)``, 
        where ``y`` is 'close' to `x` (see `__init__()`). The data are stored in
        the ``res`` attribute and the class instance is saved in a file 
        with (the weired) name ``str(func)``. 

        Parameters
        ----------
            `repetitions`
                for each point, only for noisy functions is >1 useful. For 
                ``repetitions==0`` only already generated data are plotted. 
            `locations`
                coordinated wise deviations from the middle point given in `__init__`
                
        """
        if not repetitions:
            self.plot()
            return
        
        res = self.res
        for i in range(len(self.basis)): # i-th coordinate
            if i not in res:
                res[i] = {}
            # xx = np.array(self.x)
            # TODO: store res[i]['dx'] = self.basis[i] here?
            for dx in locations:
                xx = self.x + dx * self.basis[i]
                xkey = dx  # xx[i] if (self.basis == np.eye(len(self.basis))).all() else dx  
                if xkey not in res[i]:
                    res[i][xkey] = []
                n = repetitions
                while n > 0:
                    n -= 1
                    res[i][xkey].append(self.func(xx, *self.args))
                    if plot:
                        self.plot()
                    self.save()
        return self
            
    def plot(self, plot_cmd=None, tf=lambda y: y):
        """plot the data we have, return ``self``"""
        if not plot_cmd:
            plot_cmd = self.plot_cmd
        colors = 'bgrcmyk'
        pylab.hold(False)
        res = self.res

        flatx, flatf = self.flattened()
        minf = np.inf
        for i in flatf:
            minf = min((minf, min(flatf[i])))
        addf = 1e-9 - minf  if minf <= 0 else 0
        for i in sorted(res.keys()):  # we plot not all values here
            if type(i) is int:
                color = colors[i % len(colors)]
                arx = sorted(res[i].keys()) 
                plot_cmd(arx, [tf(np.median(res[i][x]) + addf) for x in arx], color + '-')
                pylab.text(arx[-1], tf(np.median(res[i][arx[-1]])), i)
                pylab.hold(True)
                plot_cmd(flatx[i], tf(np.array(flatf[i]) + addf), color + 'o')
        pylab.ylabel('f + ' + str(addf))
        pylab.draw()
        show()
        # raw_input('press return')
        return self
 
    def flattened(self):
        """return flattened data ``(x, f)`` such that for the sweep through 
        coordinate ``i`` we have for data point ``j`` that ``f[i][j] == func(x[i][j])``
        
        """
        flatx = {}
        flatf = {}
        for i in self.res:
            if type(i) is int:
                flatx[i] = []
                flatf[i] = []
                for x in sorted(self.res[i]):
                    for d in sorted(self.res[i][x]):
                        flatx[i].append(x)
                        flatf[i].append(d)
        return flatx, flatf

    def save(self, name=None):
        """save to file"""
        import pickle
        name = name if name else self.name
        fun = self.func
        del self.func  # instance method produces error
        pickle.dump(self, open(name + '.pkl', "wb" ))
        self.func = fun
        return self

    def load(self, name=None):
        """load from file"""
        import pickle
        name = name if name else self.name
        s = pickle.load(open(name + '.pkl', 'r'))
        self.res = s.res  # disregard the class
        return self
#____________________________________________________________
#____________________________________________________________

class _Error(Exception):
    """generic exception of cma module"""
    pass


class Misc(object):
    #____________________________________________________________
    #____________________________________________________________
    # 
    class MathHelperFunctions(object):
        """static convenience math helper functions, if the function name 
        is preceded with an "a", a numpy array is returned
        
        """
        @staticmethod
        def aclamp(x, upper):
            return -Misc.MathHelperFunctions.apos(-x, -upper)
        @staticmethod
        def expms(A):  
            """matrix exponential for a symmetric matrix"""
            # TODO: check that this works reliably for low rank matrices
            # first: symmetrize A
            D, B = np.linalg.eigh(A)
            return np.dot(B, (np.exp(D) * B).T)
        @staticmethod
        def amax(vec, vec_or_scalar):
            return array(Misc.MathHelperFunctions.max(vec, vec_or_scalar))
        @staticmethod
        def max(vec, vec_or_scalar):
            b = vec_or_scalar
            if np.isscalar(b):
                m = [max(x, b) for x in vec]
            else:
                m = [max(vec[i], b[i]) for i in xrange(len(vec))]
            return m
        @staticmethod
        def amin(vec_or_scalar, vec_or_scalar2):
            return array(Misc.MathHelperFunctions.min(vec_or_scalar, vec_or_scalar2))
        @staticmethod
        def min(a, b):
            iss = np.isscalar
            if iss(a) and iss(b):
                return min(a, b)
            if iss(a):
                a, b = b, a
            # now only b can be still a scalar
            if iss(b):
                return [min(x, b) for x in a]
            else:  # two non-scalars must have the same length
                return [min(a[i], b[i]) for i in xrange(len(a))]
        @staticmethod
        def norm(vec, expo=2):
            return sum(vec**expo)**(1/expo)
        @staticmethod
        def apos(x, lower=0):
            """clips argument (scalar or array) from below at lower"""
            if lower == 0:
                return (x > 0) * x
            else: 
                return lower + (x > lower) * (x - lower)
        @staticmethod
        def prctile(data, p_vals=[0, 25, 50, 75, 100], sorted_=False):
            """``prctile(data, 50)`` returns the median, but p_vals can
            also be a sequence.
            
            Provides for small samples better values than matplotlib.mlab.prctile, 
            however also slower. 
            
            """
            ps = [p_vals] if np.isscalar(p_vals) else p_vals
                
            if not sorted_:
                data = sorted(data)
            n = len(data)
            d = []
            for p in ps:
                fi = p * n / 100 - 0.5
                if fi < 0:  # maybe extrapolate?
                    d.append(data[0])
                elif fi > n - 1:
                    d.append(data[-1])
                else:
                    i = int(fi)
                    d.append((i+1 - fi) * data[i] + (fi - i) * data[i+1])
            return d[0] if np.isscalar(p_vals) else d
        @staticmethod
        def sround(nb):  # TODO: to be vectorized
            """return stochastic round: floor(nb) + (rand()<remainder(nb))"""
            return nb // 1 + (np.random.rand(1)[0] < (nb % 1))

    @staticmethod
    def likelihood(x, m=None, Cinv=None, sigma=1, detC=None):
        """return likelihood of x for the normal density N(m, sigma**2 * Cinv**-1)"""
        # testing: MC integrate must be one: mean(p(x_i)) * volume(where x_i are uniformely sampled) 
        # for i in range(3): print mean([cma.likelihood(20*r-10, dim * [0], None, 3) for r in rand(10000,dim)]) * 20**dim
        if m is None:
            dx = x
        else:
            dx = x - m  # array(x) - array(m)
        n = len(x)
        s2pi = (2*np.pi)**(n/2.)
        if Cinv is None:
            return exp(-sum(dx**2) / sigma**2 / 2) / s2pi / sigma**n
        if detC is None:
            detC = 1. / np.linalg.linalg.det(Cinv)
        return  exp(-np.dot(dx, np.dot(Cinv, dx)) / sigma**2 / 2) / s2pi / abs(detC)**0.5 / sigma**n

    @staticmethod
    def loglikelihood(self, x, previous=False):
        """return log-likelihood of `x` regarding the current sample distribution"""
        # testing of original fct: MC integrate must be one: mean(p(x_i)) * volume(where x_i are uniformely sampled) 
        # for i in range(3): print mean([cma.likelihood(20*r-10, dim * [0], None, 3) for r in rand(10000,dim)]) * 20**dim
        # TODO: test this!!
        # c=cma.fmin...
        # c[3]['cma'].loglikelihood(...)
        
        if previous and hasattr(self, 'lastiter'):
            sigma = self.lastiter.sigma
            Crootinv = self.lastiter._Crootinv
            xmean = self.lastiter.mean
            D = self.lastiter.D
        elif previous and self.countiter > 1:
            raise _Error('no previous distribution parameters stored, check options importance_mixing')
        else:
            sigma = self.sigma
            Crootinv = self._Crootinv
            xmean = self.mean
            D = self.D
        
        dx = array(x) - xmean  # array(x) - array(m)
        n = self.N
        logs2pi = n * log(2*np.pi) / 2.
        logdetC = 2 * sum(log(D))
        dx = np.dot(Crootinv, dx)
        res = -sum(dx**2) / sigma**2 / 2 - logs2pi - logdetC/2 - n*log(sigma)
        if 1 < 3: # testing 
            s2pi = (2*np.pi)**(n/2.)
            detC = np.prod(D)**2
            res2 = -sum(dx**2) / sigma**2 / 2 - log(s2pi * abs(detC)**0.5 * sigma**n)
            assert res2 < res + 1e-8 or res2 > res - 1e-8
        return res
    
    #____________________________________________________________
    #____________________________________________________________
    # 
    # C and B are arrays rather than matrices, because they are
    # addressed via B[i][j], matrices can only be addressed via B[i,j] 
    
    # tred2(N, B, diagD, offdiag);
    # tql2(N, diagD, offdiag, B);
    
    
    # Symmetric Householder reduction to tridiagonal form, translated from JAMA package.
    @staticmethod
    def eig(C):
        """eigendecomposition of a symmetric matrix, much slower than 
        `numpy.linalg.eigh`, return ``(EVals, Basis)``, the eigenvalues 
        and an orthonormal basis of the corresponding eigenvectors, where 
        
            ``Basis[i]``
                the i-th row of ``Basis`` 
            columns of ``Basis``, ``[Basis[j][i] for j in range(len(Basis))]``
                the i-th eigenvector with eigenvalue ``EVals[i]``
        
        """
    
    # class eig(object):
    #     def __call__(self, C):
    
    # Householder transformation of a symmetric matrix V into tridiagonal form. 
        # -> n             : dimension
        # -> V             : symmetric nxn-matrix
        # <- V             : orthogonal transformation matrix:
        #                    tridiag matrix == V * V_in * V^t
        # <- d             : diagonal
        # <- e[0..n-1]     : off diagonal (elements 1..n-1)
    
        # Symmetric tridiagonal QL algorithm, iterative 
        # Computes the eigensystem from a tridiagonal matrix in roughtly 3N^3 operations
        # -> n     : Dimension. 
        # -> d     : Diagonale of tridiagonal matrix. 
        # -> e[1..n-1] : off-diagonal, output from Householder
        # -> V     : matrix output von Householder
        # <- d     : eigenvalues
        # <- e     : garbage?
        # <- V     : basis of eigenvectors, according to d
        
    
        #  tred2(N, B, diagD, offdiag); B=C on input
        #  tql2(N, diagD, offdiag, B); 
        
        #  private void tred2 (int n, double V[][], double d[], double e[]) {
        def tred2 (n, V, d, e):
            #  This is derived from the Algol procedures tred2 by
            #  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
            #  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
            #  Fortran subroutine in EISPACK.
    
            num_opt = False  # factor 1.5 in 30-D
    
            for j in range(n): 
                d[j] = V[n-1][j] # d is output argument
    
            # Householder reduction to tridiagonal form.
        
            for i in range(n-1,0,-1):
                # Scale to avoid under/overflow.
                h = 0.0
                if not num_opt:
                    scale = 0.0
                    for k in range(i):
                        scale = scale + abs(d[k])
                else:
                    scale = sum(abs(d[0:i]))
      
                if scale == 0.0:
                    e[i] = d[i-1]
                    for j in range(i):
                        d[j] = V[i-1][j]
                        V[i][j] = 0.0
                        V[j][i] = 0.0
                else:
      
                    # Generate Householder vector.
                    if not num_opt:
                        for k in range(i):
                            d[k] /= scale
                            h += d[k] * d[k]
                    else:
                        d[:i] /= scale
                        h = np.dot(d[:i],d[:i])
       
                    f = d[i-1]
                    g = h**0.5
       
                    if f > 0:
                        g = -g
       
                    e[i] = scale * g
                    h = h - f * g
                    d[i-1] = f - g
                    if not num_opt:
                        for j in range(i):
                            e[j] = 0.0
                    else:
                        e[:i] = 0.0 
           
                    # Apply similarity transformation to remaining columns.
           
                    for j in range(i): 
                        f = d[j]
                        V[j][i] = f
                        g = e[j] + V[j][j] * f
                        if not num_opt:
                            for k in range(j+1, i):
                                g += V[k][j] * d[k]
                                e[k] += V[k][j] * f
                            e[j] = g
                        else:
                            e[j+1:i] += V.T[j][j+1:i] * f
                            e[j] = g + np.dot(V.T[j][j+1:i],d[j+1:i])
       
                    f = 0.0
                    if not num_opt:
                        for j in range(i):
                            e[j] /= h
                            f += e[j] * d[j]
                    else:
                        e[:i] /= h
                        f += np.dot(e[:i],d[:i]) 
       
                    hh = f / (h + h)
                    if not num_opt:
                        for j in range(i): 
                            e[j] -= hh * d[j]
                    else:
                        e[:i] -= hh * d[:i]
                        
                    for j in range(i): 
                        f = d[j]
                        g = e[j]
                        if not num_opt:
                            for k in range(j, i): 
                                V[k][j] -= (f * e[k] + g * d[k])
                        else:
                            V.T[j][j:i] -= (f * e[j:i] + g * d[j:i])
                            
                        d[j] = V[i-1][j]
                        V[i][j] = 0.0
      
                d[i] = h
            # end for i--
        
            # Accumulate transformations.
    
            for i in range(n-1): 
                V[n-1][i] = V[i][i]
                V[i][i] = 1.0
                h = d[i+1]
                if h != 0.0:
                    if not num_opt:
                        for k in range(i+1): 
                            d[k] = V[k][i+1] / h
                    else:
                        d[:i+1] = V.T[i+1][:i+1] / h
                       
                    for j in range(i+1): 
                        if not num_opt:
                            g = 0.0
                            for k in range(i+1): 
                                g += V[k][i+1] * V[k][j]
                            for k in range(i+1): 
                                V[k][j] -= g * d[k]
                        else:
                            g = np.dot(V.T[i+1][0:i+1], V.T[j][0:i+1])
                            V.T[j][:i+1] -= g * d[:i+1]
     
                if not num_opt:
                    for k in range(i+1): 
                        V[k][i+1] = 0.0
                else:
                    V.T[i+1][:i+1] = 0.0
                   
    
            if not num_opt:
                for j in range(n): 
                    d[j] = V[n-1][j]
                    V[n-1][j] = 0.0
            else:
                d[:n] = V[n-1][:n]
                V[n-1][:n] = 0.0
    
            V[n-1][n-1] = 1.0
            e[0] = 0.0
    
    
        # Symmetric tridiagonal QL algorithm, taken from JAMA package.
        # private void tql2 (int n, double d[], double e[], double V[][]) {
        # needs roughly 3N^3 operations
        def tql2 (n, d, e, V):
    
            #  This is derived from the Algol procedures tql2, by
            #  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
            #  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
            #  Fortran subroutine in EISPACK.
        
            num_opt = False  # using vectors from numpy makes it faster
    
            if not num_opt:
                for i in range(1,n): # (int i = 1; i < n; i++): 
                    e[i-1] = e[i]
            else:
                e[0:n-1] = e[1:n]
            e[n-1] = 0.0
        
            f = 0.0
            tst1 = 0.0
            eps = 2.0**-52.0
            for l in range(n): # (int l = 0; l < n; l++) {
    
                # Find small subdiagonal element
         
                tst1 = max(tst1, abs(d[l]) + abs(e[l]))
                m = l
                while m < n: 
                    if abs(e[m]) <= eps*tst1:
                        break
                    m += 1
     
                # If m == l, d[l] is an eigenvalue,
                # otherwise, iterate.
         
                if m > l:
                    iiter = 0
                    while 1: # do {
                        iiter += 1  # (Could check iteration count here.)
         
                        # Compute implicit shift
         
                        g = d[l]
                        p = (d[l+1] - g) / (2.0 * e[l])
                        r = (p**2 + 1)**0.5  # hypot(p,1.0) 
                        if p < 0:
                            r = -r
     
                        d[l] = e[l] / (p + r)
                        d[l+1] = e[l] * (p + r)
                        dl1 = d[l+1]
                        h = g - d[l]
                        if not num_opt: 
                            for i in range(l+2, n):
                                d[i] -= h
                        else:
                            d[l+2:n] -= h
     
                        f = f + h
         
                        # Implicit QL transformation.
         
                        p = d[m]
                        c = 1.0
                        c2 = c
                        c3 = c
                        el1 = e[l+1]
                        s = 0.0
                        s2 = 0.0
     
                        # hh = V.T[0].copy()  # only with num_opt
                        for i in range(m-1, l-1, -1): # (int i = m-1; i >= l; i--) {
                            c3 = c2
                            c2 = c
                            s2 = s
                            g = c * e[i]
                            h = c * p
                            r = (p**2 + e[i]**2)**0.5  # hypot(p,e[i])
                            e[i+1] = s * r
                            s = e[i] / r
                            c = p / r
                            p = c * d[i] - s * g
                            d[i+1] = h + s * (c * g + s * d[i])
         
                            # Accumulate transformation.
         
                            if not num_opt: # overall factor 3 in 30-D
                                for k in range(n): # (int k = 0; k < n; k++) {
                                    h = V[k][i+1]
                                    V[k][i+1] = s * V[k][i] + c * h
                                    V[k][i] = c * V[k][i] - s * h
                            else: # about 20% faster in 10-D
                                hh = V.T[i+1].copy()
                                # hh[:] = V.T[i+1][:]
                                V.T[i+1] = s * V.T[i] + c * hh
                                V.T[i] = c * V.T[i] - s * hh
                                # V.T[i] *= c
                                # V.T[i] -= s * hh
     
                        p = -s * s2 * c3 * el1 * e[l] / dl1
                        e[l] = s * p
                        d[l] = c * p
         
                        # Check for convergence.
                        if abs(e[l]) <= eps*tst1:
                            break
                    # } while (Math.abs(e[l]) > eps*tst1);
     
                d[l] = d[l] + f
                e[l] = 0.0
           
    
            # Sort eigenvalues and corresponding vectors.
            if 11 < 3:
                for i in range(n-1): # (int i = 0; i < n-1; i++) {
                    k = i
                    p = d[i]
                    for j in range(i+1, n): # (int j = i+1; j < n; j++) {
                        if d[j] < p: # NH find smallest k>i
                            k = j
                            p = d[j]
    
                    if k != i: 
                        d[k] = d[i] # swap k and i 
                        d[i] = p   
                        for j in range(n): # (int j = 0; j < n; j++) {
                            p = V[j][i]
                            V[j][i] = V[j][k]
                            V[j][k] = p
        # tql2
    
        N = len(C[0])
        if 11 < 3:
            V = np.array([x[:] for x in C])  # copy each "row"
            N = V[0].size
            d = np.zeros(N)
            e = np.zeros(N)
        else:
            V = [[x[i] for i in xrange(N)] for x in C]  # copy each "row"
            d = N * [0.]
            e = N * [0.]
            
        tred2(N, V, d, e)
        tql2(N, d, e, V)
        return (array(d), array(V))

Mh = Misc.MathHelperFunctions
 
#____________________________________________________________
#____________________________________________________________
# 
def pprint(to_be_printed):
    """nicely formated print"""
    try:
        import pprint as pp
        # generate an instance PrettyPrinter
        # pp.PrettyPrinter().pprint(to_be_printed)
        pp.pprint(to_be_printed) 
    except ImportError:
        print('could not use pprint module, will apply regular print')
        print(to_be_printed)

#____________________________________________________________
#____________________________________________________________
# 
class Rotation(object):
    """Rotation class that implements an orthogonal linear transformation, 
    one for each dimension. Used to implement non-separable test functions. 

    Example:
    
    >>> import cma, numpy as np
    >>> R = cma.Rotation()
    >>> R2 = cma.Rotation() # another rotation
    >>> x = np.array((1,2,3))
    >>> print(R(R(x), inverse=1))
    [ 1.  2.  3.]
       
    """
    dicMatrices = {}  # store matrix if necessary, for each dimension
    def __init__(self):
        self.dicMatrices = {} # otherwise there might be shared bases which is probably not what we want 
    def __call__(self, x, inverse=False): # function when calling an object
        """Rotates the input array `x` with a fixed rotation matrix
           (``self.dicMatrices['str(len(x))']``)
        """
        N = x.shape[0]  # can be an array or matrix, TODO: accept also a list of arrays?
        if not self.dicMatrices.has_key(str(N)): # create new N-basis for once and all
            B = np.random.randn(N, N)
            for i in xrange(N):
                for j in xrange(0, i):
                    B[i] -= np.dot(B[i], B[j]) * B[j]
                B[i] /= sum(B[i]**2)**0.5
            self.dicMatrices[str(N)] = B
        if inverse: 
            return np.dot(self.dicMatrices[str(N)].T, x)  # compute rotation
        else: 
            return np.dot(self.dicMatrices[str(N)], x)  # compute rotation

# Use rotate(x) to rotate x 
rotate = Rotation()  # should not be capitalized? 

#____________________________________________________________
#____________________________________________________________
# 
class FitnessFunctions(object):
    """ versatile container for test objective functions """

    def __init__(self):
        self.counter = 0  # number of calls or any other practical use
    
    def rot(self, x, fun, rot=1, args=()):
        """returns ``fun(rotation(x), *args)``, ie. `fun` applied to a rotated argument"""
        if len(np.shape(array(x))) > 1:  # parallelized
            res = []
            for x in x:
                res.append(self.rot(x, fun, rot, args))
            return res
    
        if rot:
            return fun(rotate(x, *args))
        else:
            return fun(x)

    def somenan(self, x, fun, p=0.1):
        """returns sometimes np.NaN, otherwise fun(x)"""
        if np.random.rand(1) < p:
            return np.NaN
        else:
            return fun(x)
            
    def rand(self, x): 
        """Random test objective function"""
        return np.random.random(1)[0]
    
    def linear(self, x):
        return -x[0]
    
    def sphere(self, x):
        """Sphere (squared norm) test objective function"""
        # return np.random.rand(1)[0]**0 * sum(x**2) + 1 * np.random.rand(1)[0]
        return sum((x+0)**2) 

    def noisysphere(self, x, noise=5.0):
        return sum((x+0)**2) * (1 + noise * np.random.randn() / len(x))

    def spherew(self, x):
        """Sphere (squared norm) with sum x_i = 1 test objective function"""
        # return np.random.rand(1)[0]**0 * sum(x**2) + 1 * np.random.rand(1)[0]
        # s = sum(abs(x))
        # return sum((x/s+0)**2) - 1/len(x) 
        # return sum((x/s)**2) - 1/len(x) 
        return -0.01*x[0] + abs(x[0])**-2 * sum(x[1:]**2)

    def partsphere(self, x):
        """Sphere (squared norm) test objective function"""
        self.counter += 1
        # return np.random.rand(1)[0]**0 * sum(x**2) + 1 * np.random.rand(1)[0]
        dim = len(x)
        x = array([x[i % dim] for i in range(2*dim)])
        N = 8
        i = self.counter % dim
        #f = sum(x[i:i + N]**2)
        f = sum(x[np.random.randint(dim, size=N)]**2)
        return f
    
    def sectorsphere(self, x):
        """asymmetric Sphere (squared norm) test objective function"""
        return sum(x**2) + (1e6-1) * sum(x[x,0]**2)
    
    def cornersphere(self, x):
        """Sphere (squared norm) test objective function"""
        if any(x < 0):
            return np.NaN
        return sum(x**2)

    def normalSkew(self, f):
        N = np.random.randn(1)[0]**2
        if N < 1:
            N = f * N  # diminish blow up lower part
        return N
    
    def noiseC(self, x, func=sphere, fac=10, expon=0.8):
        f = func(self, x)
        N = np.random.randn(1)[0]/np.random.randn(1)[0]
        return max(1e-19, f + (float(fac)/len(x)) * f**expon * N)
    
    def noise(self, x, func=sphere, fac=10, expon=1):
        f = func(self, x)
        #R = np.random.randn(1)[0]
        R = np.log10(f) + expon * abs(10-np.log10(f)) * np.random.rand(1)[0]
        # sig = float(fac)/float(len(x))
        # R = log(f) + 0.5*log(f) * random.randn(1)[0] 
        # return max(1e-19, f + sig * (f**np.log10(f)) * np.exp(R)) 
        # return max(1e-19, f * np.exp(sig * N / f**expon)) 
        # return max(1e-19, f * normalSkew(f**expon)**sig)
        return f + 10**R  # == f + f**(1+0.5*RN)
    
    def cigar(self, x, rot=0):
        """Cigar test objective function"""
        if rot:
            x = rotate(x)  
        x = [x] if np.isscalar(x[0]) else x  # scalar into list
        f = [x[0]**2 + 1e6 * sum(x[1:]**2) for x in x]
        return f if len(f) > 1 else f[0]  # 1-element-list into scalar
        
    def tablet(self, x, rot=0):
        """Tablet test objective function"""
        if rot:
            x = rotate(x)  
        x = [x] if np.isscalar(x[0]) else x  # scalar into list
        f = [1e6*x[0]**2 + sum(x[1:]**2) for x in x]
        return f if len(f) > 1 else f[0]  # 1-element-list into scalar

    def cigtab(self, y):
        """Cigtab test objective function"""
        X = [y] if np.isscalar(y[0]) else y
        f = [1e-4 * x[0]**2 + 1e4 * x[1]**2 + sum(x[2:]**2) for x in X]
        return f if len(f) > 1 else f[0]
        
    def twoaxes(self, y):
        """Cigtab test objective function"""
        X = [y] if np.isscalar(y[0]) else y
        N2 = len(X[0]) // 2
        f = [1e6 * sum(x[0:N2]**2) + sum(x[N2:]**2) for x in X]
        return f if len(f) > 1 else f[0]
     
    def ellirot(self, x):
        return fcts.elli(array(x), 1)
       
    def elli(self, x, rot=0, xoffset=0, cond=1e6, actuator_noise=0.0, both=False):
        """Ellipsoid test objective function"""
        if not np.isscalar(x[0]):  # parallel evaluation
            return [self.elli(xi, rot) for xi in x]  # could save 20% overall
        if rot:
            x = rotate(x)  
        N = len(x)
        if actuator_noise:  
            x = x + actuator_noise * np.random.randn(N)
            
        ftrue = sum(cond**(np.arange(N)/(N-1.))*(x+xoffset)**2)

        alpha = 0.49 + 1./N 
        beta = 1      
        felli = np.random.rand(1)[0]**beta * ftrue * \
                max(1, (10.**9 / (ftrue+1e-99))**(alpha*np.random.rand(1)[0]))
        # felli = ftrue + 1*np.random.randn(1)[0] / (1e-30 +
        #                                           np.abs(np.random.randn(1)[0]))**0
        if both:
            return (felli, ftrue)
        else:
            return ftrue # + np.random.randn()
            # return felli
        
    def elliconstraint(self, x, cfac = 1e8, tough=True, cond=1e6):  
        """ellipsoid test objective function with "constraints" """
        N = len(x)
        f = sum(cond**(np.arange(N)[-1::-1]/(N-1)) * x**2)
        cvals = (x[0] + 1, 
                 x[0] + 1 + 100*x[1], 
                 x[0] + 1 - 100*x[1])
        if tough:
            f += cfac * sum(max(0,c) for c in cvals)
        else: 
            f += cfac * sum(max(0,c+1e-3)**2 for c in cvals)
        return f 

    def rosen(self, x):  
        """Rosenbrock test objective function"""
        x = [x] if np.isscalar(x[0]) else x  # scalar into list
        f = [sum(100.*(x[:-1]**2-x[1:])**2 + (1.-x[:-1])**2) for x in x]
        return f if len(f) > 1 else f[0]  # 1-element-list into scalar
        
    def diffpow(self, x, rot=0):
        """Diffpow test objective function"""
        N = len(x)
        if rot:
            x = rotate(x)  
        return sum(np.abs(x)**(2.+4.*np.arange(N)/(N-1.)))**0.5
        
    def ridge(self, x, expo=2):
        x = [x] if np.isscalar(x[0]) else x  # scalar into list
        f = [x[0] + 100*np.sum(x[1:]**2)**(expo/2.) for x in x]
        return f if len(f) > 1 else f[0]  # 1-element-list into scalar
        
    def flat(self,x):
        return 1 
        return 1 if np.random.rand(1) < 0.9 else 1.1
        return np.random.randint(1,30)
    def griewank(self, x):
        # was in [-600 600]
        x = (600./5) * x
        return 1 - np.prod(np.cos(x/sqrt(1.+np.arange(len(x))))) + sum(x**2)/4e3

    def rastrigin(self, x):
        """Rastrigin test objective function"""
        if not np.isscalar(x[0]):
            N = len(x[0])
            return [10*N + sum(xi**2 - 10*np.cos(2*np.pi*xi)) for xi in x]
            # return 10*N + sum(x**2 - 10*np.cos(2*np.pi*x), axis=1)
        N = len(x)
        return 10*N + sum(x**2 - 10*np.cos(2*np.pi*x))
        
    def schwefelmult(self, x, pen_fac = 1e4):
        """multimodal Schwefel function with domain -500..500"""
        y = [x] if np.isscalar(x[0]) else x
        N = len(y[0]) 
        f = array([418.9829*N - 1.27275661e-5*N - sum(x * np.sin(np.abs(x)**0.5))
                + pen_fac * sum((abs(x) > 500) * (abs(x) - 500)**2) for x in y])
        return f if len(f) > 1 else f[0]
        
    def optprob(self, x):
        n = np.arange(len(x)) + 1
        f = n * x * (1-x)**(n-1)
        return sum(1-f)
    
    def lincon(self, x, theta=0.01):
        """ridge like linear function with one linear constraint"""
        if x[0] < 0:
            return np.NaN
        return theta * x[1] + x[0]
        
fcts = FitnessFunctions()
Fcts = fcts  # for cross compatibility, as if the functions were static members of class Fcts
 
#____________________________________________
#____________________________________________________________

def _test(module=None):  # None is fine when called from inside the module
    import doctest
    print(doctest.testmod(module))  # this is pretty coool!

def process_test(stream=None):
    """ """
    import fileinput
    s1 = ""
    s2 = ""
    s3 = ""
    state = 0
    for line in fileinput.input(stream):  # takes argv as file or stdin
        if 1 < 3:
            s3 += line
            if state < -1 and line.startswith('***'):
                print(s3)
            if line.startswith('***'):
                s3 = ""
            
        if state == -1:  # found a failed example line
            s1 += '\n\n*** Failed Example:' + line
            s2 += '\n\n\n'   # line
            # state = 0  # wait for 'Expected:' line 
            
        if line.startswith('Expected:'):
            state = 1
            continue
        elif line.startswith('Got:'):
            state = 2
            continue
        elif line.startswith('***'):  # marks end of failed example
            state = 0
        elif line.startswith('Failed example:'):
            state = -1
        elif line.startswith('Exception raised'):
            state = -2

        # in effect more else:
        if state == 1:
            s1 += line + ''
        if state == 2:
            s2 += line + ''

#____________________________________________________________
#____________________________________________________________
#
def main(argv=None):
    """provides a versatile interface to do some testing from the command line.
    
    Usage: python cma.py [options | func dim sig0 [optkey optval][optkey optval]...]
    Use option -t or --test to run the doctest, -t -v to get (much) verbosity
    and -t -q to run it quietly with output only in case of errors. 
    Use options --fcts and --doc for more infos or start ipython --pylab. 
    
    Examples 
    --------
    A single run: python cma.py elli 10 1 
    
    Run doctest: python cma.py --test --quiet

    """
    if argv is None:
        argv = sys.argv  # should have better been sys.argv[1:]

    # uncomment for unit test
    # _test() 
    # handle input arguments, getopt might be helpful ;-)
    if len(argv) >= 1:  # function and help
        if len(argv) == 1 or argv[1].startswith('-h') or argv[1].startswith('--help'):
            print(main.__doc__)
            fun = None
        elif argv[1].startswith('-t') or argv[1].startswith('--test'):
            import doctest
            if len(argv) > 2 and (argv[2].startswith('--qu') or argv[2].startswith('-q')):
                print('doctest for cma.py: launching (it might be necessary to close a few pop up windows to finish)')
                fn = '__cma_doctest__.txt'
                stdout = sys.stdout
                try:
                    with open(fn, 'w') as f:
                        sys.stdout = f
                        doctest.testmod(report=True)  # this is quite cool!
                finally:
                    sys.stdout = stdout
                process_test(fn)
                print('doctest for cma.py: finished (no other output should be seen after launching)')
            else:
                print('doctest for cma.py: due to different platforms and python versions')
                print('and in some cases due to a missing unique random seed')
                print('many examples will "fail". This is OK, if they give a similar')
                print('to the expected result and if no exception occurs. ')
                # if argv[1][2] == 'v':
                doctest.testmod(report=True)  # this is quite cool!
            return
        elif argv[1] == '--doc':
            print(__doc__)
            print(CMAEvolutionStrategy.__doc__)
            print(fmin.__doc__)
            fun = None
        elif argv[1] == '--fcts':
            print('List of valid function names:')
            print([d for d in dir(fcts) if not d.startswith('_')])
            fun = None
        elif len(argv) > 3:
            fun = eval('fcts.' + argv[1])
        else:
            print('try -h option')
            fun = None
            
    if fun is not None:
        
        if len(argv) > 2:  # dimension
            x0 = np.ones(eval(argv[2]))
        if len(argv) > 3:  # sigma
            sig0 = eval(argv[3])

        opts = {}
        for i in xrange(5, len(argv), 2):
            opts[argv[i-1]] = eval(argv[i])

        # run fmin
        if fun is not None:
            tic = time.time()
            fmin(fun, x0, sig0, **opts)  # ftarget=1e-9, tolfacupx=1e9, verb_log=10)
            # plot()
            # print ' best function value ', res[2]['es'].best[1]
            print('elapsed time [s]: + %.2f', round(time.time() - tic, 2))
    
    elif not len(argv):
        fmin(fcts.elli, np.ones(6)*0.1, 0.1, ftarget=1e-9)


#____________________________________________________________
#____________________________________________________________
#
# mainly for testing purpose
# executed when called from an OS shell
if __name__ == "__main__":
    # for i in range(1000):  # how to find the memory leak
    #     main(["cma.py", "rastrigin", "10", "5", "popsize", "200", "maxfevals", "24999", "verb_log", "0"])
    main()

