docs: def ask(self, number=None, xmean=None, sigma_fac=1,
              gradf=None, args=()):
          """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, phenotyp?
              `sigma_fac`
                  multiplier for internal sample width (standard
                  deviation)
              `gradf`
                  gradient, ``len(gradf(x)) == len(x)``, if
                  ``gradf is not None`` the third solution in the
                  returned list is "sampled" in supposedly Newton
                  direction ``dot(C, gradf(xmean, *args))``.
              `args`
                  additional arguments passed to gradf
  
          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_geno = self.ask_geno(number, xmean, sigma_fac)
  
          # N,lambda=20,200: overall CPU 7s vs 5s == 40% overhead, even without bounds!
          #                  new data: 11.5s vs 9.5s == 20%
          # TODO: check here, whether this is necessary?
          # return [self.gp.pheno(x, copy=False, into_bounds=self.boundary_handler.repair) 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
          pop_pheno = [self.gp.pheno(x, copy=True, into_bounds=self.boundary_handler.repair) for x in pop_geno]
  
          if gradf is not None:
              # see Hansen (2011), Injecting external solutions into CMA-ES
              if not self.gp.islinear:
                  _print_warning("""
                  using the gradient (option ``gradf``) with a non-linear
                  coordinate-wise transformation (option ``transformation``)
                  has never been tested.""")
                  # TODO: check this out
              def grad_numerical_of_coordinate_map(x, map, epsilon=None):
                  """map is a coordinate-wise independent map, return
                  the estimated diagonal of the Jacobian.
                  """
                  eps = 1e-8 * (1 + abs(x)) if epsilon is None else epsilon
                  return (map(x + eps) - map(x - eps)) / (2 * eps)
              def grad_numerical_sym(x, func, epsilon=None):
                  """return symmetric numerical gradient of func : R^n -> R.
                  """
                  eps = 1e-8 * (1 + abs(x)) if epsilon is None else epsilon
                  grad = np.zeros(len(x))
                  ei = np.zeros(len(x))  # float is 1.6 times faster than int
                  for i in rglen(x):
                      ei[i] = eps[i]
                      grad[i] = (func(x + ei) - func(x - ei)) / (2*eps[i])
                      ei[i] = 0
                  return grad
              try:
                  if self.last_iteration_with_gradient == self.countiter:
                      _print_warning('gradient is used several times in ' +
                              'this iteration', iteration=self.countiter)
                  self.last_iteration_with_gradient = self.countiter
              except AttributeError:
                  pass
              index_for_gradient = min((2, len(pop_pheno)-1))
              xmean = self.mean if xmean is None else xmean
              xpheno = self.gp.pheno(xmean, copy=True,
                                  into_bounds=self.boundary_handler.repair)
              grad_at_mean = gradf(xpheno, *args)
              # lift gradient into geno-space
              if not self.gp.isidentity or (self.boundary_handler is not None
                      and self.boundary_handler.has_bounds()):
                  boundary_repair = None
                  gradpen = 0
                  if isinstance(self.boundary_handler, BoundTransform):
                      boundary_repair = self.boundary_handler.repair
                  elif isinstance(self.boundary_handler, BoundPenalty):
                      fpenalty = lambda x: self.boundary_handler.__call__(
                          x, SolutionDict({tuple(x): {'geno': x}}), self.gp)
                      gradpen = grad_numerical_sym(
                          xmean, fpenalty)
                  elif self.boundary_handler is None or \
                          isinstance(self.boundary_handler, BoundNone):
                      pass
                  else:
                      raise NotImplementedError(
                          "unknown boundary handling method" +
                          str(self.boundary_handler) +
                          " when using gradf")
                  gradgp = grad_numerical_of_coordinate_map(
                      xmean,
                      lambda x: self.gp.pheno(x, copy=True,
                              into_bounds=boundary_repair))
                  grad_at_mean = grad_at_mean * gradgp + gradpen
  
              # TODO: frozen variables brake the code (e.g. at grad of map)
              if len(grad_at_mean) != self.N and self.opts['fixed_variables']:
                  NotImplementedError("""
                  gradient with fixed variables is not yet implemented""")
              v = self.D * dot(self.B.T, self.sigma_vec * grad_at_mean)
              # newton_direction = sv * B * D * D * B^T * sv * gradient = sv * B * D * v
              # v = D^-1 * B^T * sv^-1 * newton_direction = D * B^T * sv * gradient
              q = sum(v**2)
              if q:
                  # Newton direction
                  pop_geno[index_for_gradient] = xmean - self.sigma \
                              * (self.N / q)**0.5 \
                              * (self.sigma_vec * dot(self.B, self.D * v))
              else:
                  pop_geno[index_for_gradient] = xmean
                  _print_warning('gradient zero observed',
                                 iteration=self.countiter)
  
              pop_pheno[index_for_gradient] = self.gp.pheno(
                  pop_geno[index_for_gradient], copy=True,
                  into_bounds=self.boundary_handler.repair)
  
          # insert solutions, this could also (better?) be done in self.gp.pheno
          for i in rglen((pop_geno)):
              self.sent_solutions.insert(pop_pheno[i], geno=pop_geno[i], iteration=self.countiter)
          return pop_pheno