docs: def tell(self, solutions, function_values, 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, in particular the two evolution paths, the
          distribution mean, the covariance matrix and a step-size.
  
          Arguments
          ---------
              `solutions`
                  list or array of candidate solution points (of
                  type `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 ``check_points is None``, only solutions that are not generated
                  by `ask()` are possibly clipped (recommended). ``False`` does not clip
                  any solution (not recommended).
                  If ``True``, clips solutions that realize long steps (i.e. also
                  those that are unlikely to be generated with `ask()`). `check_points`
                  can be a list of indices to be checked in solutions.
              `copy`
                  ``solutions`` can be modified in this routine, if ``copy is False``
  
          Details
          -------
          `tell()` updates the parameters of the multivariate
          normal search distribution, 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()`
  
          """
          if self._flgtelldone:
              raise _Error('tell should only be called once per iteration')
  
          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')
  
          if not isscalar(function_values[0]):
              if 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 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, self.sent_solutions, function_values, self.countevals)
  
          flg_diagonal = self.opts['CMA_diagonal'] is True \
                         or self.countiter <= self.opts['CMA_diagonal']
          if not flg_diagonal and len(self.C.shape) == 1:  # C was diagonal ie 1-D
              # enter non-separable phase (no easy return from here)
              self.C = np.diag(self.C)
              if 1 < 3:
                  self.B = np.eye(N)  # identity(N)
                  idx = np.argsort(self.D)
                  self.D = self.D[idx]
                  self.B = self.B[:, idx]
              self._Yneg = 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.boundary_handler.update(function_values, self)(solutions, self.sent_solutions, self.gp)
          # for testing:
          # fit.bndpen = self.boundary_handler.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
          # now get the genotypes
          pop = self.pop_sorted = []  # create pop from input argument solutions
          for k, s in enumerate(solutions):  # use phenotype before Solution.repair()
              if 1 < 3:
                  pop += [self.gp.geno(s,
                              from_bounds=self.boundary_handler.inverse,
                              repair=(self.repair_genotype if check_points not in (False, 0, [], ()) else None),
                              archive=self.sent_solutions)]  # takes genotype from sent_solutions, if available
                  try:
                      self.archive.insert(s, value=self.sent_solutions.pop(s), fitness=function_values[k])
                      # self.sent_solutions.pop(s)
                  except KeyError:
                      pass
          # check that TPA mirrors are available, TODO: move to TPA class?
          if isinstance(self.adapt_sigma, CMAAdaptSigmaTPA) and self.countiter > 3 and not (self.countiter % 3):
              dm = self.mean[0] - self.mean_old[0]
              dx0 = pop[0][0] - self.mean_old[0]
              dx1 = pop[1][0] - self.mean_old[0]
              for i in np.random.randint(1, self.N, 1):
                  try:
                      if not Mh.equals_approximately(
                                  (self.mean[i] - self.mean_old[i])
                                      / (pop[0][i] - self.mean_old[i]),
                                   dm / dx0, 1e-8) or \
                          not Mh.equals_approximately(
                                  (self.mean[i] - self.mean_old[i])
                                      / (pop[1][i] - self.mean_old[i]),
                                  dm / dx1, 1e-8):
                          _print_warning('TPA error with mirrored samples', 'tell',
                                         'CMAEvolutionStrategy', self.countiter)
                  except ZeroDivisionError:
                      _print_warning('zero division encountered in TPA check\n which should be very rare and is likely a bug',
                                     'tell', 'CMAEvolutionStrategy', self.countiter)
  
          try:
              moldold = self.mean_old
          except:
              pass
          self.mean_old = self.mean
          mold = self.mean_old  # just an alias
  
          # check and normalize each x - m
          # check_points is a flag (None is default: check non-known solutions) 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, however is automatic with use_sent_solutions
              try:
                  if len(check_points):
                      idx = check_points
              except:
                  idx = xrange(sp.popsize)
  
              for k in idx:
                  self.repair_genotype(pop[k])
  
          # only arrays can be multiple indexed
          pop = array(pop, copy=False)
  
          # sort pop
          pop = pop[fit.idx]
  
          # prepend best-ever solution to population, in case
          if self.opts['CMA_elitist'] and self.best.f < fit.fit[0]:
              if self.best.x_geno is not None:
                  xp = [self.best.x_geno]
                  # xp = [self.best.xdict['geno']]
                  # xp = [self.gp.geno(self.best.x[:])]  # TODO: remove
                  # print self.mahalanobis_norm(xp[0]-self.mean)
              else:
                  xp = [self.gp.geno(array(self.best.x, copy=True),
                                     self.boundary_handler.inverse,
                                     copy_if_changed=False)]
                  print('genotype for elitist not found')
              self.clip_or_fit_solutions(xp, [0])
              pop = array([xp[0]] + list(pop))
          elif self.opts['CMA_elitist'] == 'initial':  # current solution was better
              self.opts['CMA_elitist'] = False
  
          self.pop_sorted = 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.
          # replaced by repair_geno?
          # simple test case injecting self.mean:
          # self.mean = 1e-4 * self.sigma * np.random.randn(N)
          if 1 < 3:
              cmean = self.sp.cmean
  
          # zzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz
          # get learning rate constants
          cc, c1, cmu = sp.cc, sp.c1, sp.cmu
          if flg_diagonal:
              cc, c1, cmu = sp.cc_sep, sp.c1_sep, sp.cmu_sep
  
          # now the real work can start
  
          hsig = self.adapt_sigma.hsig(self) # ps update must be done here in separable case
  
          # hsig = sum(self.ps**2) / 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
  
          self.pc = (1 - cc) * self.pc + \
                    hsig * (sqrt(cc * (2 - cc) * sp.mueff) / self.sigma / cmean) * \
                    (self.mean - mold) / self.sigma_vec
  
          # 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 flg_diagonal:
                  Y = (pop[0:sp.mu] - mold) / (self.sigma * self.sigma_vec)
                  Y = dot((cmu * sp.weights) * Y.T, Y)  # learning rate integrated
                  if self.sp.neg.cmuexp:
                      tmp = (pop[-sp.neg.mu:] - mold) / (self.sigma * self.sigma_vec)
                      # normalize to constant length (seems preferable in several aspects)
                      for i in xrange(tmp.shape[0]):
                          tmp[i, :] *= N**0.5 / self.mahalanobis_norm(
                                   self.sigma_vec * tmp[i, :]) / self.sigma
                      self._Yneg *= 1 - self.sp.neg.cmuexp  # for some reason necessary?
                      self._Yneg += 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)
  
                  self.C *= 1 - c1a - cmu
                  self.C += np.outer(c1 * self.pc, self.pc) + Y
                  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 * self.sigma_vec)  # 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, this is why adapt_sigma needs to prepare before
                  self.itereigenupdated = self.countiter
  
                  # idx = self._mirror_idx_cov()  # take half of mirrored vectors for negative update
  
          # step-size adaptation, adapt sigma
          # in case of TPA, function_values[0] and [1] must reflect samples colinear to xmean - xmean_old
          self.adapt_sigma.update(self, function_values=function_values)
  
          if self.sigma * min(self.sigma_vec * self.dC**0.5) < self.opts['minstd']:
              self.sigma = self.opts['minstd'] / min(self.sigma_vec * self.dC**0.5)
          if self.sigma * max(self.sigma_vec * self.dC**0.5) > self.opts['maxstd']:
              self.sigma = self.opts['maxstd'] / max(self.sigma_vec * self.dC**0.5)
          # g = self.countiter
          # N = self.N
          # mindx = eval(self.opts['mindx'])
          #  if isinstance(self.opts['mindx'], basestring) else self.opts['mindx']
          if self.sigma * min(self.D) < self.opts['mindx']:  # TODO: sigma_vec is missing here
              self.sigma = self.opts['mindx'] / min(self.D)
  
          if self.sigma > 1e9 * self.sigma0:
              alpha = self.sigma / max(self.D)
              self.multiplyC(alpha)
              self.sigma /= alpha**0.5
              self.opts['tolupsigma'] /= alpha**0.5  # to be compared with sigma
  
          # TODO increase sigma in case of a plateau?
  
          # Uncertainty noise measurement is done on an upper level
  
          # move mean into "feasible preimage", leads to weird behavior on
          # 40-D tablet with bound 0.1, not quite explained (constant
          # dragging is problematic, but why doesn't it settle), still a bug?
          if new_injections:
              self.pop_injection_directions = self.prepare_injection_directions()
          self.pop_sorted = []  # remove this in case pop is still needed
          self._flgtelldone = True