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
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`