    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017
    018    package org.apache.commons.math3.optim.nonlinear.scalar.noderiv;
    019
    020    import java.util.ArrayList;
    021    import java.util.Arrays;
    022    import java.util.List;
    023    import org.apache.commons.math3.exception.DimensionMismatchException;
    024    import org.apache.commons.math3.exception.NotPositiveException;
    025    import org.apache.commons.math3.exception.NotStrictlyPositiveException;
    026    import org.apache.commons.math3.exception.OutOfRangeException;
    027    import org.apache.commons.math3.exception.TooManyEvaluationsException;
    028    import org.apache.commons.math3.linear.Array2DRowRealMatrix;
    029    import org.apache.commons.math3.linear.EigenDecomposition;
    030    import org.apache.commons.math3.linear.MatrixUtils;
    031    import org.apache.commons.math3.linear.RealMatrix;
    032    import org.apache.commons.math3.optim.ConvergenceChecker;
    033    import org.apache.commons.math3.optim.OptimizationData;
    034    import org.apache.commons.math3.optim.nonlinear.scalar.GoalType;
    035    import org.apache.commons.math3.optim.PointValuePair;
    036    import org.apache.commons.math3.optim.nonlinear.scalar.MultivariateOptimizer;
    037    import org.apache.commons.math3.random.RandomGenerator;
    038    import org.apache.commons.math3.util.MathArrays;
    039
    040    /**
    041     * <p>An implementation of the active Covariance Matrix Adaptation Evolution Strategy (CMA-ES)
    042     * for non-linear, non-convex, non-smooth, global function minimization.
    043     * The CMA-Evolution Strategy (CMA-ES) is a reliable stochastic optimization method
    044     * which should be applied if derivative-based methods, e.g. quasi-Newton BFGS or
    045     * conjugate gradient, fail due to a rugged search landscape (e.g. noise, local
    046     * optima, outlier, etc.) of the objective function. Like a
    047     * quasi-Newton method, the CMA-ES learns and applies a variable metric
    048     * on the underlying search space. Unlike a quasi-Newton method, the
    049     * CMA-ES neither estimates nor uses gradients, making it considerably more
    050     * reliable in terms of finding a good, or even close to optimal, solution.</p>
    051     *
    052     * <p>In general, on smooth objective functions the CMA-ES is roughly ten times
    053     * slower than BFGS (counting objective function evaluations, no gradients provided).
    054     * For up to <math>N=10</math> variables also the derivative-free simplex
    055     * direct search method (Nelder and Mead) can be faster, but it is
    056     * far less reliable than CMA-ES.</p>
    057     *
    058     * <p>The CMA-ES is particularly well suited for non-separable
    059     * and/or badly conditioned problems. To observe the advantage of CMA compared
    060     * to a conventional evolution strategy, it will usually take about
    061     * <math>30 N</math> function evaluations. On difficult problems the complete
    062     * optimization (a single run) is expected to take <em>roughly</em> between
    063     * <math>30 N</math> and <math>300 N<sup>2</sup></math>
    064     * function evaluations.</p>
    065     *
    066     * <p>This implementation is translated and adapted from the Matlab version
    067     * of the CMA-ES algorithm as implemented in module {@code cmaes.m} version 3.51.</p>
    068     *
    069     * For more information, please refer to the following links:
    070     * <ul>
    071     *  <li><a href="http://www.lri.fr/~hansen/cmaes.m">Matlab code</a></li>
    072     *  <li><a href="http://www.lri.fr/~hansen/cmaesintro.html">Introduction to CMA-ES</a></li>
    073     *  <li><a href="http://en.wikipedia.org/wiki/CMA-ES">Wikipedia</a></li>
    074     * </ul>
    075     *
    076     * @version $Id: CMAESOptimizer.java 1400108 2012-10-19 14:20:16Z erans $
    077     * @since 3.0
    078     */
    079    public class CMAESOptimizer
    080        extends MultivariateOptimizer {
    081        // global search parameters
    082        /**
    083         * Population size, offspring number. The primary strategy parameter to play
    084         * with, which can be increased from its default value. Increasing the
    085         * population size improves global search properties in exchange to speed.
    086         * Speed decreases, as a rule, at most linearly with increasing population
    087         * size. It is advisable to begin with the default small population size.
    088         */
    089        private int lambda; // population size
    090        /**
    091         * Covariance update mechanism, default is active CMA. isActiveCMA = true
    092         * turns on "active CMA" with a negative update of the covariance matrix and
    093         * checks for positive definiteness. OPTS.CMA.active = 2 does not check for
    094         * pos. def. and is numerically faster. Active CMA usually speeds up the
    095         * adaptation.
    096         */
    097        private final boolean isActiveCMA;
    098        /**
    099         * Determines how often a new random offspring is generated in case it is
    100         * not feasible / beyond the defined limits, default is 0.
    101         */
    102        private final int checkFeasableCount;
    103        /**
    104         * @see Sigma
    105         */
    106        private double[] inputSigma;
    107        /** Number of objective variables/problem dimension */
    108        private int dimension;
    109        /**
    110         * Defines the number of initial iterations, where the covariance matrix
    111         * remains diagonal and the algorithm has internally linear time complexity.
    112         * diagonalOnly = 1 means keeping the covariance matrix always diagonal and
    113         * this setting also exhibits linear space complexity. This can be
    114         * particularly useful for dimension > 100.
    115         * @see <a href="http://hal.archives-ouvertes.fr/inria-00287367/en">A Simple Modification in CMA-ES</a>
    116         */
    117        private int diagonalOnly;
    118        /** Number of objective variables/problem dimension */
    119        private boolean isMinimize = true;
    120        /** Indicates whether statistic data is collected. */
    121        private final boolean generateStatistics;
    122
    123        // termination criteria
    124        /** Maximal number of iterations allowed. */
    125        private final int maxIterations;
    126        /** Limit for fitness value. */
    127        private final double stopFitness;
    128        /** Stop if x-changes larger stopTolUpX. */
    129        private double stopTolUpX;
    130        /** Stop if x-change smaller stopTolX. */
    131        private double stopTolX;
    132        /** Stop if fun-changes smaller stopTolFun. */
    133        private double stopTolFun;
    134        /** Stop if back fun-changes smaller stopTolHistFun. */
    135        private double stopTolHistFun;
    136
    137        // selection strategy parameters
    138        /** Number of parents/points for recombination. */
    139        private int mu; //
    140        /** log(mu + 0.5), stored for efficiency. */
    141        private double logMu2;
    142        /** Array for weighted recombination. */
    143        private RealMatrix weights;
    144        /** Variance-effectiveness of sum w_i x_i. */
    145        private double mueff; //
    146
    147        // dynamic strategy parameters and constants
    148        /** Overall standard deviation - search volume. */
    149        private double sigma;
    150        /** Cumulation constant. */
    151        private double cc;
    152        /** Cumulation constant for step-size. */
    153        private double cs;
    154        /** Damping for step-size. */
    155        private double damps;
    156        /** Learning rate for rank-one update. */
    157        private double ccov1;
    158        /** Learning rate for rank-mu update' */
    159        private double ccovmu;
    160        /** Expectation of ||N(0,I)|| == norm(randn(N,1)). */
    161        private double chiN;
    162        /** Learning rate for rank-one update - diagonalOnly */
    163        private double ccov1Sep;
    164        /** Learning rate for rank-mu update - diagonalOnly */
    165        private double ccovmuSep;
    166
    167        // CMA internal values - updated each generation
    168        /** Objective variables. */
    169        private RealMatrix xmean;
    170        /** Evolution path. */
    171        private RealMatrix pc;
    172        /** Evolution path for sigma. */
    173        private RealMatrix ps;
    174        /** Norm of ps, stored for efficiency. */
    175        private double normps;
    176        /** Coordinate system. */
    177        private RealMatrix B;
    178        /** Scaling. */
    179        private RealMatrix D;
    180        /** B*D, stored for efficiency. */
    181        private RealMatrix BD;
    182        /** Diagonal of sqrt(D), stored for efficiency. */
    183        private RealMatrix diagD;
    184        /** Covariance matrix. */
    185        private RealMatrix C;
    186        /** Diagonal of C, used for diagonalOnly. */
    187        private RealMatrix diagC;
    188        /** Number of iterations already performed. */
    189        private int iterations;
    190
    191        /** History queue of best values. */
    192        private double[] fitnessHistory;
    193        /** Size of history queue of best values. */
    194        private int historySize;
    195
    196        /** Random generator. */
    197        private final RandomGenerator random;
    198
    199        /** History of sigma values. */
    200        private final List<Double> statisticsSigmaHistory = new ArrayList<Double>();
    201        /** History of mean matrix. */
    202        private final List<RealMatrix> statisticsMeanHistory = new ArrayList<RealMatrix>();
    203        /** History of fitness values. */
    204        private final List<Double> statisticsFitnessHistory = new ArrayList<Double>();
    205        /** History of D matrix. */
    206        private final List<RealMatrix> statisticsDHistory = new ArrayList<RealMatrix>();
    207
    208        /**
    209         * @param maxIterations Maximal number of iterations.
    210         * @param stopFitness Whether to stop if objective function value is smaller than
    211         * {@code stopFitness}.
    212         * @param isActiveCMA Chooses the covariance matrix update method.
    213         * @param diagonalOnly Number of initial iterations, where the covariance matrix
    214         * remains diagonal.
    215         * @param checkFeasableCount Determines how often new random objective variables are
    216         * generated in case they are out of bounds.
    217         * @param random Random generator.
    218         * @param generateStatistics Whether statistic data is collected.
    219         * @param checker Convergence checker.
    220         *
    221         * @since 3.1
    222         */
    223        public CMAESOptimizer(int maxIterations,
    224                              double stopFitness,
    225                              boolean isActiveCMA,
    226                              int diagonalOnly,
    227                              int checkFeasableCount,
    228                              RandomGenerator random,
    229                              boolean generateStatistics,
    230                              ConvergenceChecker<PointValuePair> checker) {
    231            super(checker);
    232            this.maxIterations = maxIterations;
    233            this.stopFitness = stopFitness;
    234            this.isActiveCMA = isActiveCMA;
    235            this.diagonalOnly = diagonalOnly;
    236            this.checkFeasableCount = checkFeasableCount;
    237            this.random = random;
    238            this.generateStatistics = generateStatistics;
    239        }
    240
    241        /**
    242         * @return History of sigma values.
    243         */
    244        public List<Double> getStatisticsSigmaHistory() {
    245            return statisticsSigmaHistory;
    246        }
    247
    248        /**
    249         * @return History of mean matrix.
    250         */
    251        public List<RealMatrix> getStatisticsMeanHistory() {
    252            return statisticsMeanHistory;
    253        }
    254
    255        /**
    256         * @return History of fitness values.
    257         */
    258        public List<Double> getStatisticsFitnessHistory() {
    259            return statisticsFitnessHistory;
    260        }
    261
    262        /**
    263         * @return History of D matrix.
    264         */
    265        public List<RealMatrix> getStatisticsDHistory() {
    266            return statisticsDHistory;
    267        }
    268
    269        /**
    270         * Input sigma values.
    271         * They define the initial coordinate-wise standard deviations for
    272         * sampling new search points around the initial guess.
    273         * It is suggested to set them to the estimated distance from the
    274         * initial to the desired optimum.
    275         * Small values induce the search to be more local (and very small
    276         * values are more likely to find a local optimum close to the initial
    277         * guess).
    278         * Too small values might however lead to early termination.
    279         */
    280        public static class Sigma implements OptimizationData {
    281            /** Sigma values. */
    282            private final double[] sigma;
    283
    284            /**
    285             * @param s Sigma values.
    286             * @throws NotPositiveException if any of the array entries is smaller
    287             * than zero.
    288             */
    289            public Sigma(double[] s)
    290                throws NotPositiveException {
    291                for (int i = 0; i < s.length; i++) {
    292                    if (s[i] < 0) {
    293                        throw new NotPositiveException(s[i]);
    294                    }
    295                }
    296
    297                sigma = s.clone();
    298            }
    299
    300            /**
    301             * @return the sigma values.
    302             */
    303            public double[] getSigma() {
    304                return sigma.clone();
    305            }
    306        }
    307
    308        /**
    309         * Population size.
    310         * The number of offspring is the primary strategy parameter.
    311         * In the absence of better clues, a good default could be an
    312         * integer close to {@code 4 + 3 ln(n)}, where {@code n} is the
    313         * number of optimized parameters.
    314         * Increasing the population size improves global search properties
    315         * at the expense of speed (which in general decreases at most
    316         * linearly with increasing population size).
    317         */
    318        public static class PopulationSize implements OptimizationData {
    319            /** Population size. */
    320            private final int lambda;
    321
    322            /**
    323             * @param size Population size.
    324             * @throws NotStrictlyPositiveException if {@code size <= 0}.
    325             */
    326            public PopulationSize(int size)
    327                throws NotStrictlyPositiveException {
    328                if (size <= 0) {
    329                    throw new NotStrictlyPositiveException(size);
    330                }
    331                lambda = size;
    332            }
    333
    334            /**
    335             * @return the population size.
    336             */
    337            public int getPopulationSize() {
    338                return lambda;
    339            }
    340        }
    341
    342        /**
    343         * {@inheritDoc}
    344         *
    345         * @param optData Optimization data. The following data will be looked for:
    346         * <ul>
    347         *  <li>{@link org.apache.commons.math3.optim.MaxEval}</li>
    348         *  <li>{@link org.apache.commons.math3.optim.InitialGuess}</li>
    349         *  <li>{@link org.apache.commons.math3.optim.SimpleBounds}</li>
    350         *  <li>{@link org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunction}</li>
    351         *  <li>{@link Sigma}</li>
    352         *  <li>{@link PopulationSize}</li>
    353         * </ul>
    354         * @return {@inheritDoc}
    355         * @throws TooManyEvaluationsException if the maximal number of
    356         * evaluations is exceeded.
    357         * @throws DimensionMismatchException if the initial guess, target, and weight
    358         * arguments have inconsistent dimensions.
    359         */
    360        @Override
    361        public PointValuePair optimize(OptimizationData... optData)
    362            throws TooManyEvaluationsException,
    363                   DimensionMismatchException {
    364            // Retrieve settings.
    365            parseOptimizationData(optData);
    366            // Set up base class and perform computation.
    367            return super.optimize(optData);
    368        }
    369
    370        /** {@inheritDoc} */
    371        @Override
    372        protected PointValuePair doOptimize() {
    373            checkParameters();
    374             // -------------------- Initialization --------------------------------
    375            isMinimize = getGoalType().equals(GoalType.MINIMIZE);
    376            final FitnessFunction fitfun = new FitnessFunction();
    377            final double[] guess = getStartPoint();
    378            // number of objective variables/problem dimension
    379            dimension = guess.length;
    380            initializeCMA(guess);
    381            iterations = 0;
    382            double bestValue = fitfun.value(guess);
    383            push(fitnessHistory, bestValue);
    384            PointValuePair optimum
    385                = new PointValuePair(getStartPoint(),
    386                                     isMinimize ? bestValue : -bestValue);
    387            PointValuePair lastResult = null;
    388
    389            // -------------------- Generation Loop --------------------------------
    390
    391            generationLoop:
    392            for (iterations = 1; iterations <= maxIterations; iterations++) {
    393                // Generate and evaluate lambda offspring
    394                final RealMatrix arz = randn1(dimension, lambda);
    395                final RealMatrix arx = zeros(dimension, lambda);
    396                final double[] fitness = new double[lambda];
    397                // generate random offspring
    398                for (int k = 0; k < lambda; k++) {
    399                    RealMatrix arxk = null;
    400                    for (int i = 0; i < checkFeasableCount + 1; i++) {
    401                        if (diagonalOnly <= 0) {
    402                            arxk = xmean.add(BD.multiply(arz.getColumnMatrix(k))
    403                                             .scalarMultiply(sigma)); // m + sig * Normal(0,C)
    404                        } else {
    405                            arxk = xmean.add(times(diagD,arz.getColumnMatrix(k))
    406                                             .scalarMultiply(sigma));
    407                        }
    408                        if (i >= checkFeasableCount ||
    409                            fitfun.isFeasible(arxk.getColumn(0))) {
    410                            break;
    411                        }
    412                        // regenerate random arguments for row
    413                        arz.setColumn(k, randn(dimension));
    414                    }
    415                    copyColumn(arxk, 0, arx, k);
    416                    try {
    417                        fitness[k] = fitfun.value(arx.getColumn(k)); // compute fitness
    418                    } catch (TooManyEvaluationsException e) {
    419                        break generationLoop;
    420                    }
    421                }
    422                // Sort by fitness and compute weighted mean into xmean
    423                final int[] arindex = sortedIndices(fitness);
    424                // Calculate new xmean, this is selection and recombination
    425                final RealMatrix xold = xmean; // for speed up of Eq. (2) and (3)
    426                final RealMatrix bestArx = selectColumns(arx, MathArrays.copyOf(arindex, mu));
    427                xmean = bestArx.multiply(weights);
    428                final RealMatrix bestArz = selectColumns(arz, MathArrays.copyOf(arindex, mu));
    429                final RealMatrix zmean = bestArz.multiply(weights);
    430                final boolean hsig = updateEvolutionPaths(zmean, xold);
    431                if (diagonalOnly <= 0) {
    432                    updateCovariance(hsig, bestArx, arz, arindex, xold);
    433                } else {
    434                    updateCovarianceDiagonalOnly(hsig, bestArz);
    435                }
    436                // Adapt step size sigma - Eq. (5)
    437                sigma *= Math.exp(Math.min(1, (normps/chiN - 1) * cs / damps));
    438                final double bestFitness = fitness[arindex[0]];
    439                final double worstFitness = fitness[arindex[arindex.length - 1]];
    440                if (bestValue > bestFitness) {
    441                    bestValue = bestFitness;
    442                    lastResult = optimum;
    443                    optimum = new PointValuePair(fitfun.repair(bestArx.getColumn(0)),
    444                                                 isMinimize ? bestFitness : -bestFitness);
    445                    if (getConvergenceChecker() != null &&
    446                        lastResult != null) {
    447                        if (getConvergenceChecker().converged(iterations, optimum, lastResult)) {
    448                            break generationLoop;
    449                        }
    450                    }
    451                }
    452                // handle termination criteria
    453                // Break, if fitness is good enough
    454                if (stopFitness != 0) { // only if stopFitness is defined
    455                    if (bestFitness < (isMinimize ? stopFitness : -stopFitness)) {
    456                        break generationLoop;
    457                    }
    458                }
    459                final double[] sqrtDiagC = sqrt(diagC).getColumn(0);
    460                final double[] pcCol = pc.getColumn(0);
    461                for (int i = 0; i < dimension; i++) {
    462                    if (sigma * Math.max(Math.abs(pcCol[i]), sqrtDiagC[i]) > stopTolX) {
    463                        break;
    464                    }
    465                    if (i >= dimension - 1) {
    466                        break generationLoop;
    467                    }
    468                }
    469                for (int i = 0; i < dimension; i++) {
    470                    if (sigma * sqrtDiagC[i] > stopTolUpX) {
    471                        break generationLoop;
    472                    }
    473                }
    474                final double historyBest = min(fitnessHistory);
    475                final double historyWorst = max(fitnessHistory);
    476                if (iterations > 2 &&
    477                    Math.max(historyWorst, worstFitness) -
    478                    Math.min(historyBest, bestFitness) < stopTolFun) {
    479                    break generationLoop;
    480                }
    481                if (iterations > fitnessHistory.length &&
    482                    historyWorst - historyBest < stopTolHistFun) {
    483                    break generationLoop;
    484                }
    485                // condition number of the covariance matrix exceeds 1e14
    486                if (max(diagD) / min(diagD) > 1e7) {
    487                    break generationLoop;
    488                }
    489                // user defined termination
    490                if (getConvergenceChecker() != null) {
    491                    final PointValuePair current
    492                        = new PointValuePair(bestArx.getColumn(0),
    493                                             isMinimize ? bestFitness : -bestFitness);
    494                    if (lastResult != null &&
    495                        getConvergenceChecker().converged(iterations, current, lastResult)) {
    496                        break generationLoop;
    497                        }
    498                    lastResult = current;
    499                }
    500                // Adjust step size in case of equal function values (flat fitness)
    501                if (bestValue == fitness[arindex[(int)(0.1+lambda/4.)]]) {
    502                    sigma = sigma * Math.exp(0.2 + cs / damps);
    503                }
    504                if (iterations > 2 && Math.max(historyWorst, bestFitness) -
    505                    Math.min(historyBest, bestFitness) == 0) {
    506                    sigma = sigma * Math.exp(0.2 + cs / damps);
    507                }
    508                // store best in history
    509                push(fitnessHistory,bestFitness);
    510                fitfun.setValueRange(worstFitness-bestFitness);
    511                if (generateStatistics) {
    512                    statisticsSigmaHistory.add(sigma);
    513                    statisticsFitnessHistory.add(bestFitness);
    514                    statisticsMeanHistory.add(xmean.transpose());
    515                    statisticsDHistory.add(diagD.transpose().scalarMultiply(1E5));
    516                }
    517            }
    518            return optimum;
    519        }
    520
    521        /**
    522         * Scans the list of (required and optional) optimization data that
    523         * characterize the problem.
    524         *
    525         * @param optData Optimization data. The following data will be looked for:
    526         * <ul>
    527         *  <li>{@link Sigma}</li>
    528         *  <li>{@link PopulationSize}</li>
    529         * </ul>
    530         */
    531        private void parseOptimizationData(OptimizationData... optData) {
    532            // The existing values (as set by the previous call) are reused if
    533            // not provided in the argument list.
    534            for (OptimizationData data : optData) {
    535                if (data instanceof Sigma) {
    536                    inputSigma = ((Sigma) data).getSigma();
    537                    continue;
    538                }
    539                if (data instanceof PopulationSize) {
    540                    lambda = ((PopulationSize) data).getPopulationSize();
    541                    continue;
    542                }
    543            }
    544        }
    545
    546        /**
    547         * Checks dimensions and values of boundaries and inputSigma if defined.
    548         */
    549        private void checkParameters() {
    550            final double[] init = getStartPoint();
    551            final double[] lB = getLowerBound();
    552            final double[] uB = getUpperBound();
    553
    554            if (inputSigma != null) {
    555                if (inputSigma.length != init.length) {
    556                    throw new DimensionMismatchException(inputSigma.length, init.length);
    557                }
    558                for (int i = 0; i < init.length; i++) {
    559                    if (inputSigma[i] > uB[i] - lB[i]) {
    560                        throw new OutOfRangeException(inputSigma[i], 0, uB[i] - lB[i]);
    561                    }
    562                }
    563            }
    564        }
    565
    566        /**
    567         * Initialization of the dynamic search parameters
    568         *
    569         * @param guess Initial guess for the arguments of the fitness function.
    570         */
    571        private void initializeCMA(double[] guess) {
    572            if (lambda <= 0) {
    573                throw new NotStrictlyPositiveException(lambda);
    574            }
    575            // initialize sigma
    576            final double[][] sigmaArray = new double[guess.length][1];
    577            for (int i = 0; i < guess.length; i++) {
    578                sigmaArray[i][0] = inputSigma[i];
    579            }
    580            final RealMatrix insigma = new Array2DRowRealMatrix(sigmaArray, false);
    581            sigma = max(insigma); // overall standard deviation
    582
    583            // initialize termination criteria
    584            stopTolUpX = 1e3 * max(insigma);
    585            stopTolX = 1e-11 * max(insigma);
    586            stopTolFun = 1e-12;
    587            stopTolHistFun = 1e-13;
    588
    589            // initialize selection strategy parameters
    590            mu = lambda / 2; // number of parents/points for recombination
    591            logMu2 = Math.log(mu + 0.5);
    592            weights = log(sequence(1, mu, 1)).scalarMultiply(-1).scalarAdd(logMu2);
    593            double sumw = 0;
    594            double sumwq = 0;
    595            for (int i = 0; i < mu; i++) {
    596                double w = weights.getEntry(i, 0);
    597                sumw += w;
    598                sumwq += w * w;
    599            }
    600            weights = weights.scalarMultiply(1 / sumw);
    601            mueff = sumw * sumw / sumwq; // variance-effectiveness of sum w_i x_i
    602
    603            // initialize dynamic strategy parameters and constants
    604            cc = (4 + mueff / dimension) /
    605                    (dimension + 4 + 2 * mueff / dimension);
    606            cs = (mueff + 2) / (dimension + mueff + 3.);
    607            damps = (1 + 2 * Math.max(0, Math.sqrt((mueff - 1) /
    608                                                   (dimension + 1)) - 1)) *
    609                Math.max(0.3,
    610                         1 - dimension / (1e-6 + maxIterations)) + cs; // minor increment
    611            ccov1 = 2 / ((dimension + 1.3) * (dimension + 1.3) + mueff);
    612            ccovmu = Math.min(1 - ccov1, 2 * (mueff - 2 + 1 / mueff) /
    613                              ((dimension + 2) * (dimension + 2) + mueff));
    614            ccov1Sep = Math.min(1, ccov1 * (dimension + 1.5) / 3);
    615            ccovmuSep = Math.min(1 - ccov1, ccovmu * (dimension + 1.5) / 3);
    616            chiN = Math.sqrt(dimension) *
    617                (1 - 1 / ((double) 4 * dimension) + 1 / ((double) 21 * dimension * dimension));
    618            // intialize CMA internal values - updated each generation
    619            xmean = MatrixUtils.createColumnRealMatrix(guess); // objective variables
    620            diagD = insigma.scalarMultiply(1 / sigma);
    621            diagC = square(diagD);
    622            pc = zeros(dimension, 1); // evolution paths for C and sigma
    623            ps = zeros(dimension, 1); // B defines the coordinate system
    624            normps = ps.getFrobeniusNorm();
    625
    626            B = eye(dimension, dimension);
    627            D = ones(dimension, 1); // diagonal D defines the scaling
    628            BD = times(B, repmat(diagD.transpose(), dimension, 1));
    629            C = B.multiply(diag(square(D)).multiply(B.transpose())); // covariance
    630            historySize = 10 + (int) (3 * 10 * dimension / (double) lambda);
    631            fitnessHistory = new double[historySize]; // history of fitness values
    632            for (int i = 0; i < historySize; i++) {
    633                fitnessHistory[i] = Double.MAX_VALUE;
    634            }
    635        }
    636
    637        /**
    638         * Update of the evolution paths ps and pc.
    639         *
    640         * @param zmean Weighted row matrix of the gaussian random numbers generating
    641         * the current offspring.
    642         * @param xold xmean matrix of the previous generation.
    643         * @return hsig flag indicating a small correction.
    644         */
    645        private boolean updateEvolutionPaths(RealMatrix zmean, RealMatrix xold) {
    646            ps = ps.scalarMultiply(1 - cs).add(
    647                    B.multiply(zmean).scalarMultiply(
    648                            Math.sqrt(cs * (2 - cs) * mueff)));
    649            normps = ps.getFrobeniusNorm();
    650            final boolean hsig = normps /
    651                Math.sqrt(1 - Math.pow(1 - cs, 2 * iterations)) /
    652                chiN < 1.4 + 2 / ((double) dimension + 1);
    653            pc = pc.scalarMultiply(1 - cc);
    654            if (hsig) {
    655                pc = pc.add(xmean.subtract(xold).scalarMultiply(Math.sqrt(cc * (2 - cc) * mueff) / sigma));
    656            }
    657            return hsig;
    658        }
    659
    660        /**
    661         * Update of the covariance matrix C for diagonalOnly > 0
    662         *
    663         * @param hsig Flag indicating a small correction.
    664         * @param bestArz Fitness-sorted matrix of the gaussian random values of the
    665         * current offspring.
    666         */
    667        private void updateCovarianceDiagonalOnly(boolean hsig,
    668                                                  final RealMatrix bestArz) {
    669            // minor correction if hsig==false
    670            double oldFac = hsig ? 0 : ccov1Sep * cc * (2 - cc);
    671            oldFac += 1 - ccov1Sep - ccovmuSep;
    672            diagC = diagC.scalarMultiply(oldFac) // regard old matrix
    673                .add(square(pc).scalarMultiply(ccov1Sep)) // plus rank one update
    674                .add((times(diagC, square(bestArz).multiply(weights))) // plus rank mu update
    675                     .scalarMultiply(ccovmuSep));
    676            diagD = sqrt(diagC); // replaces eig(C)
    677            if (diagonalOnly > 1 &&
    678                iterations > diagonalOnly) {
    679                // full covariance matrix from now on
    680                diagonalOnly = 0;
    681                B = eye(dimension, dimension);
    682                BD = diag(diagD);
    683                C = diag(diagC);
    684            }
    685        }
    686
    687        /**
    688         * Update of the covariance matrix C.
    689         *
    690         * @param hsig Flag indicating a small correction.
    691         * @param bestArx Fitness-sorted matrix of the argument vectors producing the
    692         * current offspring.
    693         * @param arz Unsorted matrix containing the gaussian random values of the
    694         * current offspring.
    695         * @param arindex Indices indicating the fitness-order of the current offspring.
    696         * @param xold xmean matrix of the previous generation.
    697         */
    698        private void updateCovariance(boolean hsig, final RealMatrix bestArx,
    699                                      final RealMatrix arz, final int[] arindex,
    700                                      final RealMatrix xold) {
    701            double negccov = 0;
    702            if (ccov1 + ccovmu > 0) {
    703                final RealMatrix arpos = bestArx.subtract(repmat(xold, 1, mu))
    704                    .scalarMultiply(1 / sigma); // mu difference vectors
    705                final RealMatrix roneu = pc.multiply(pc.transpose())
    706                    .scalarMultiply(ccov1); // rank one update
    707                // minor correction if hsig==false
    708                double oldFac = hsig ? 0 : ccov1 * cc * (2 - cc);
    709                oldFac += 1 - ccov1 - ccovmu;
    710                if (isActiveCMA) {
    711                    // Adapt covariance matrix C active CMA
    712                    negccov = (1 - ccovmu) * 0.25 * mueff /
    713                        (Math.pow(dimension + 2, 1.5) + 2 * mueff);
    714                    // keep at least 0.66 in all directions, small popsize are most
    715                    // critical
    716                    final double negminresidualvariance = 0.66;
    717                    // where to make up for the variance loss
    718                    final double negalphaold = 0.5;
    719                    // prepare vectors, compute negative updating matrix Cneg
    720                    final int[] arReverseIndex = reverse(arindex);
    721                    RealMatrix arzneg = selectColumns(arz, MathArrays.copyOf(arReverseIndex, mu));
    722                    RealMatrix arnorms = sqrt(sumRows(square(arzneg)));
    723                    final int[] idxnorms = sortedIndices(arnorms.getRow(0));
    724                    final RealMatrix arnormsSorted = selectColumns(arnorms, idxnorms);
    725                    final int[] idxReverse = reverse(idxnorms);
    726                    final RealMatrix arnormsReverse = selectColumns(arnorms, idxReverse);
    727                    arnorms = divide(arnormsReverse, arnormsSorted);
    728                    final int[] idxInv = inverse(idxnorms);
    729                    final RealMatrix arnormsInv = selectColumns(arnorms, idxInv);
    730                    // check and set learning rate negccov
    731                    final double negcovMax = (1 - negminresidualvariance) /
    732                        square(arnormsInv).multiply(weights).getEntry(0, 0);
    733                    if (negccov > negcovMax) {
    734                        negccov = negcovMax;
    735                    }
    736                    arzneg = times(arzneg, repmat(arnormsInv, dimension, 1));
    737                    final RealMatrix artmp = BD.multiply(arzneg);
    738                    final RealMatrix Cneg = artmp.multiply(diag(weights)).multiply(artmp.transpose());
    739                    oldFac += negalphaold * negccov;
    740                    C = C.scalarMultiply(oldFac)
    741                        .add(roneu) // regard old matrix
    742                        .add(arpos.scalarMultiply( // plus rank one update
    743                                                  ccovmu + (1 - negalphaold) * negccov) // plus rank mu update
    744                             .multiply(times(repmat(weights, 1, dimension),
    745                                             arpos.transpose())))
    746                        .subtract(Cneg.scalarMultiply(negccov));
    747                } else {
    748                    // Adapt covariance matrix C - nonactive
    749                    C = C.scalarMultiply(oldFac) // regard old matrix
    750                        .add(roneu) // plus rank one update
    751                        .add(arpos.scalarMultiply(ccovmu) // plus rank mu update
    752                             .multiply(times(repmat(weights, 1, dimension),
    753                                             arpos.transpose())));
    754                }
    755            }
    756            updateBD(negccov);
    757        }
    758
    759        /**
    760         * Update B and D from C.
    761         *
    762         * @param negccov Negative covariance factor.
    763         */
    764        private void updateBD(double negccov) {
    765            if (ccov1 + ccovmu + negccov > 0 &&
    766                (iterations % 1. / (ccov1 + ccovmu + negccov) / dimension / 10.) < 1) {
    767                // to achieve O(N^2)
    768                C = triu(C, 0).add(triu(C, 1).transpose());
    769                // enforce symmetry to prevent complex numbers
    770                final EigenDecomposition eig = new EigenDecomposition(C);
    771                B = eig.getV(); // eigen decomposition, B==normalized eigenvectors
    772                D = eig.getD();
    773                diagD = diag(D);
    774                if (min(diagD) <= 0) {
    775                    for (int i = 0; i < dimension; i++) {
    776                        if (diagD.getEntry(i, 0) < 0) {
    777                            diagD.setEntry(i, 0, 0);
    778                        }
    779                    }
    780                    final double tfac = max(diagD) / 1e14;
    781                    C = C.add(eye(dimension, dimension).scalarMultiply(tfac));
    782                    diagD = diagD.add(ones(dimension, 1).scalarMultiply(tfac));
    783                }
    784                if (max(diagD) > 1e14 * min(diagD)) {
    785                    final double tfac = max(diagD) / 1e14 - min(diagD);
    786                    C = C.add(eye(dimension, dimension).scalarMultiply(tfac));
    787                    diagD = diagD.add(ones(dimension, 1).scalarMultiply(tfac));
    788                }
    789                diagC = diag(C);
    790                diagD = sqrt(diagD); // D contains standard deviations now
    791                BD = times(B, repmat(diagD.transpose(), dimension, 1)); // O(n^2)
    792            }
    793        }
    794
    795        /**
    796         * Pushes the current best fitness value in a history queue.
    797         *
    798         * @param vals History queue.
    799         * @param val Current best fitness value.
    800         */
    801        private static void push(double[] vals, double val) {
    802            for (int i = vals.length-1; i > 0; i--) {
    803                vals[i] = vals[i-1];
    804            }
    805            vals[0] = val;
    806        }
    807
    808        /**
    809         * Sorts fitness values.
    810         *
    811         * @param doubles Array of values to be sorted.
    812         * @return a sorted array of indices pointing into doubles.
    813         */
    814        private int[] sortedIndices(final double[] doubles) {
    815            final DoubleIndex[] dis = new DoubleIndex[doubles.length];
    816            for (int i = 0; i < doubles.length; i++) {
    817                dis[i] = new DoubleIndex(doubles[i], i);
    818            }
    819            Arrays.sort(dis);
    820            final int[] indices = new int[doubles.length];
    821            for (int i = 0; i < doubles.length; i++) {
    822                indices[i] = dis[i].index;
    823            }
    824            return indices;
    825        }
    826
    827        /**
    828         * Used to sort fitness values. Sorting is always in lower value first
    829         * order.
    830         */
    831        private static class DoubleIndex implements Comparable<DoubleIndex> {
    832            /** Value to compare. */
    833            private final double value;
    834            /** Index into sorted array. */
    835            private final int index;
    836
    837            /**
    838             * @param value Value to compare.
    839             * @param index Index into sorted array.
    840             */
    841            DoubleIndex(double value, int index) {
    842                this.value = value;
    843                this.index = index;
    844            }
    845
    846            /** {@inheritDoc} */
    847            public int compareTo(DoubleIndex o) {
    848                return Double.compare(value, o.value);
    849            }
    850
    851            /** {@inheritDoc} */
    852            @Override
    853            public boolean equals(Object other) {
    854
    855                if (this == other) {
    856                    return true;
    857                }
    858
    859                if (other instanceof DoubleIndex) {
    860                    return Double.compare(value, ((DoubleIndex) other).value) == 0;
    861                }
    862
    863                return false;
    864            }
    865
    866            /** {@inheritDoc} */
    867            @Override
    868            public int hashCode() {
    869                long bits = Double.doubleToLongBits(value);
    870                return (int) ((1438542 ^ (bits >>> 32) ^ bits) & 0xffffffff);
    871            }
    872        }
    873
    874        /**
    875         * Normalizes fitness values to the range [0,1]. Adds a penalty to the
    876         * fitness value if out of range. The penalty is adjusted by calling
    877         * setValueRange().
    878         */
    879        private class FitnessFunction {
    880            /** Determines the penalty for boundary violations */
    881            private double valueRange;
    882            /**
    883             * Flag indicating whether the objective variables are forced into their
    884             * bounds if defined
    885             */
    886            private final boolean isRepairMode;
    887
    888            /** Simple constructor.
    889             */
    890            public FitnessFunction() {
    891                valueRange = 1;
    892                isRepairMode = true;
    893            }
    894
    895            /**
    896             * @param point Normalized objective variables.
    897             * @return the objective value + penalty for violated bounds.
    898             */
    899            public double value(final double[] point) {
    900                double value;
    901                if (isRepairMode) {
    902                    double[] repaired = repair(point);
    903                    value = CMAESOptimizer.this.computeObjectiveValue(repaired) +
    904                        penalty(point, repaired);
    905                } else {
    906                    value = CMAESOptimizer.this.computeObjectiveValue(point);
    907                }
    908                return isMinimize ? value : -value;
    909            }
    910
    911            /**
    912             * @param x Normalized objective variables.
    913             * @return {@code true} if in bounds.
    914             */
    915            public boolean isFeasible(final double[] x) {
    916                final double[] lB = CMAESOptimizer.this.getLowerBound();
    917                final double[] uB = CMAESOptimizer.this.getUpperBound();
    918
    919                for (int i = 0; i < x.length; i++) {
    920                    if (x[i] < lB[i]) {
    921                        return false;
    922                    }
    923                    if (x[i] > uB[i]) {
    924                        return false;
    925                    }
    926                }
    927                return true;
    928            }
    929
    930            /**
    931             * @param valueRange Adjusts the penalty computation.
    932             */
    933            public void setValueRange(double valueRange) {
    934                this.valueRange = valueRange;
    935            }
    936
    937            /**
    938             * @param x Normalized objective variables.
    939             * @return the repaired (i.e. all in bounds) objective variables.
    940             */
    941            private double[] repair(final double[] x) {
    942                final double[] lB = CMAESOptimizer.this.getLowerBound();
    943                final double[] uB = CMAESOptimizer.this.getUpperBound();
    944
    945                final double[] repaired = new double[x.length];
    946                for (int i = 0; i < x.length; i++) {
    947                    if (x[i] < lB[i]) {
    948                        repaired[i] = lB[i];
    949                    } else if (x[i] > uB[i]) {
    950                        repaired[i] = uB[i];
    951                    } else {
    952                        repaired[i] = x[i];
    953                    }
    954                }
    955                return repaired;
    956            }
    957
    958            /**
    959             * @param x Normalized objective variables.
    960             * @param repaired Repaired objective variables.
    961             * @return Penalty value according to the violation of the bounds.
    962             */
    963            private double penalty(final double[] x, final double[] repaired) {
    964                double penalty = 0;
    965                for (int i = 0; i < x.length; i++) {
    966                    double diff = Math.abs(x[i] - repaired[i]);
    967                    penalty += diff * valueRange;
    968                }
    969                return isMinimize ? penalty : -penalty;
    970            }
    971        }
    972
    973        // -----Matrix utility functions similar to the Matlab build in functions------
    974
    975        /**
    976         * @param m Input matrix
    977         * @return Matrix representing the element-wise logarithm of m.
    978         */
    979        private static RealMatrix log(final RealMatrix m) {
    980            final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
    981            for (int r = 0; r < m.getRowDimension(); r++) {
    982                for (int c = 0; c < m.getColumnDimension(); c++) {
    983                    d[r][c] = Math.log(m.getEntry(r, c));
    984                }
    985            }
    986            return new Array2DRowRealMatrix(d, false);
    987        }
    988
    989        /**
    990         * @param m Input matrix.
    991         * @return Matrix representing the element-wise square root of m.
    992         */
    993        private static RealMatrix sqrt(final RealMatrix m) {
    994            final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
    995            for (int r = 0; r < m.getRowDimension(); r++) {
    996                for (int c = 0; c < m.getColumnDimension(); c++) {
    997                    d[r][c] = Math.sqrt(m.getEntry(r, c));
    998                }
    999            }
    1000            return new Array2DRowRealMatrix(d, false);
    1001        }
    1002
    1003        /**
    1004         * @param m Input matrix.
    1005         * @return Matrix representing the element-wise square of m.
    1006         */
    1007        private static RealMatrix square(final RealMatrix m) {
    1008            final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
    1009            for (int r = 0; r < m.getRowDimension(); r++) {
    1010                for (int c = 0; c < m.getColumnDimension(); c++) {
    1011                    double e = m.getEntry(r, c);
    1012                    d[r][c] = e * e;
    1013                }
    1014            }
    1015            return new Array2DRowRealMatrix(d, false);
    1016        }
    1017
    1018        /**
    1019         * @param m Input matrix 1.
    1020         * @param n Input matrix 2.
    1021         * @return the matrix where the elements of m and n are element-wise multiplied.
    1022         */
    1023        private static RealMatrix times(final RealMatrix m, final RealMatrix n) {
    1024            final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
    1025            for (int r = 0; r < m.getRowDimension(); r++) {
    1026                for (int c = 0; c < m.getColumnDimension(); c++) {
    1027                    d[r][c] = m.getEntry(r, c) * n.getEntry(r, c);
    1028                }
    1029            }
    1030            return new Array2DRowRealMatrix(d, false);
    1031        }
    1032
    1033        /**
    1034         * @param m Input matrix 1.
    1035         * @param n Input matrix 2.
    1036         * @return Matrix where the elements of m and n are element-wise divided.
    1037         */
    1038        private static RealMatrix divide(final RealMatrix m, final RealMatrix n) {
    1039            final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
    1040            for (int r = 0; r < m.getRowDimension(); r++) {
    1041                for (int c = 0; c < m.getColumnDimension(); c++) {
    1042                    d[r][c] = m.getEntry(r, c) / n.getEntry(r, c);
    1043                }
    1044            }
    1045            return new Array2DRowRealMatrix(d, false);
    1046        }
    1047
    1048        /**
    1049         * @param m Input matrix.
    1050         * @param cols Columns to select.
    1051         * @return Matrix representing the selected columns.
    1052         */
    1053        private static RealMatrix selectColumns(final RealMatrix m, final int[] cols) {
    1054            final double[][] d = new double[m.getRowDimension()][cols.length];
    1055            for (int r = 0; r < m.getRowDimension(); r++) {
    1056                for (int c = 0; c < cols.length; c++) {
    1057                    d[r][c] = m.getEntry(r, cols[c]);
    1058                }
    1059            }
    1060            return new Array2DRowRealMatrix(d, false);
    1061        }
    1062
    1063        /**
    1064         * @param m Input matrix.
    1065         * @param k Diagonal position.
    1066         * @return Upper triangular part of matrix.
    1067         */
    1068        private static RealMatrix triu(final RealMatrix m, int k) {
    1069            final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
    1070            for (int r = 0; r < m.getRowDimension(); r++) {
    1071                for (int c = 0; c < m.getColumnDimension(); c++) {
    1072                    d[r][c] = r <= c - k ? m.getEntry(r, c) : 0;
    1073                }
    1074            }
    1075            return new Array2DRowRealMatrix(d, false);
    1076        }
    1077
    1078        /**
    1079         * @param m Input matrix.
    1080         * @return Row matrix representing the sums of the rows.
    1081         */
    1082        private static RealMatrix sumRows(final RealMatrix m) {
    1083            final double[][] d = new double[1][m.getColumnDimension()];
    1084            for (int c = 0; c < m.getColumnDimension(); c++) {
    1085                double sum = 0;
    1086                for (int r = 0; r < m.getRowDimension(); r++) {
    1087                    sum += m.getEntry(r, c);
    1088                }
    1089                d[0][c] = sum;
    1090            }
    1091            return new Array2DRowRealMatrix(d, false);
    1092        }
    1093
    1094        /**
    1095         * @param m Input matrix.
    1096         * @return the diagonal n-by-n matrix if m is a column matrix or the column
    1097         * matrix representing the diagonal if m is a n-by-n matrix.
    1098         */
    1099        private static RealMatrix diag(final RealMatrix m) {
    1100            if (m.getColumnDimension() == 1) {
    1101                final double[][] d = new double[m.getRowDimension()][m.getRowDimension()];
    1102                for (int i = 0; i < m.getRowDimension(); i++) {
    1103                    d[i][i] = m.getEntry(i, 0);
    1104                }
    1105                return new Array2DRowRealMatrix(d, false);
    1106            } else {
    1107                final double[][] d = new double[m.getRowDimension()][1];
    1108                for (int i = 0; i < m.getColumnDimension(); i++) {
    1109                    d[i][0] = m.getEntry(i, i);
    1110                }
    1111                return new Array2DRowRealMatrix(d, false);
    1112            }
    1113        }
    1114
    1115        /**
    1116         * Copies a column from m1 to m2.
    1117         *
    1118         * @param m1 Source matrix.
    1119         * @param col1 Source column.
    1120         * @param m2 Target matrix.
    1121         * @param col2 Target column.
    1122         */
    1123        private static void copyColumn(final RealMatrix m1, int col1,
    1124                                       RealMatrix m2, int col2) {
    1125            for (int i = 0; i < m1.getRowDimension(); i++) {
    1126                m2.setEntry(i, col2, m1.getEntry(i, col1));
    1127            }
    1128        }
    1129
    1130        /**
    1131         * @param n Number of rows.
    1132         * @param m Number of columns.
    1133         * @return n-by-m matrix filled with 1.
    1134         */
    1135        private static RealMatrix ones(int n, int m) {
    1136            final double[][] d = new double[n][m];
    1137            for (int r = 0; r < n; r++) {
    1138                Arrays.fill(d[r], 1);
    1139            }
    1140            return new Array2DRowRealMatrix(d, false);
    1141        }
    1142
    1143        /**
    1144         * @param n Number of rows.
    1145         * @param m Number of columns.
    1146         * @return n-by-m matrix of 0 values out of diagonal, and 1 values on
    1147         * the diagonal.
    1148         */
    1149        private static RealMatrix eye(int n, int m) {
    1150            final double[][] d = new double[n][m];
    1151            for (int r = 0; r < n; r++) {
    1152                if (r < m) {
    1153                    d[r][r] = 1;
    1154                }
    1155            }
    1156            return new Array2DRowRealMatrix(d, false);
    1157        }
    1158
    1159        /**
    1160         * @param n Number of rows.
    1161         * @param m Number of columns.
    1162         * @return n-by-m matrix of zero values.
    1163         */
    1164        private static RealMatrix zeros(int n, int m) {
    1165            return new Array2DRowRealMatrix(n, m);
    1166        }
    1167
    1168        /**
    1169         * @param mat Input matrix.
    1170         * @param n Number of row replicates.
    1171         * @param m Number of column replicates.
    1172         * @return a matrix which replicates the input matrix in both directions.
    1173         */
    1174        private static RealMatrix repmat(final RealMatrix mat, int n, int m) {
    1175            final int rd = mat.getRowDimension();
    1176            final int cd = mat.getColumnDimension();
    1177            final double[][] d = new double[n * rd][m * cd];
    1178            for (int r = 0; r < n * rd; r++) {
    1179                for (int c = 0; c < m * cd; c++) {
    1180                    d[r][c] = mat.getEntry(r % rd, c % cd);
    1181                }
    1182            }
    1183            return new Array2DRowRealMatrix(d, false);
    1184        }
    1185
    1186        /**
    1187         * @param start Start value.
    1188         * @param end End value.
    1189         * @param step Step size.
    1190         * @return a sequence as column matrix.
    1191         */
    1192        private static RealMatrix sequence(double start, double end, double step) {
    1193            final int size = (int) ((end - start) / step + 1);
    1194            final double[][] d = new double[size][1];
    1195            double value = start;
    1196            for (int r = 0; r < size; r++) {
    1197                d[r][0] = value;
    1198                value += step;
    1199            }
    1200            return new Array2DRowRealMatrix(d, false);
    1201        }
    1202
    1203        /**
    1204         * @param m Input matrix.
    1205         * @return the maximum of the matrix element values.
    1206         */
    1207        private static double max(final RealMatrix m) {
    1208            double max = -Double.MAX_VALUE;
    1209            for (int r = 0; r < m.getRowDimension(); r++) {
    1210                for (int c = 0; c < m.getColumnDimension(); c++) {
    1211                    double e = m.getEntry(r, c);
    1212                    if (max < e) {
    1213                        max = e;
    1214                    }
    1215                }
    1216            }
    1217            return max;
    1218        }
    1219
    1220        /**
    1221         * @param m Input matrix.
    1222         * @return the minimum of the matrix element values.
    1223         */
    1224        private static double min(final RealMatrix m) {
    1225            double min = Double.MAX_VALUE;
    1226            for (int r = 0; r < m.getRowDimension(); r++) {
    1227                for (int c = 0; c < m.getColumnDimension(); c++) {
    1228                    double e = m.getEntry(r, c);
    1229                    if (min > e) {
    1230                        min = e;
    1231                    }
    1232                }
    1233            }
    1234            return min;
    1235        }
    1236
    1237        /**
    1238         * @param m Input array.
    1239         * @return the maximum of the array values.
    1240         */
    1241        private static double max(final double[] m) {
    1242            double max = -Double.MAX_VALUE;
    1243            for (int r = 0; r < m.length; r++) {
    1244                if (max < m[r]) {
    1245                    max = m[r];
    1246                }
    1247            }
    1248            return max;
    1249        }
    1250
    1251        /**
    1252         * @param m Input array.
    1253         * @return the minimum of the array values.
    1254         */
    1255        private static double min(final double[] m) {
    1256            double min = Double.MAX_VALUE;
    1257            for (int r = 0; r < m.length; r++) {
    1258                if (min > m[r]) {
    1259                    min = m[r];
    1260                }
    1261            }
    1262            return min;
    1263        }
    1264
    1265        /**
    1266         * @param indices Input index array.
    1267         * @return the inverse of the mapping defined by indices.
    1268         */
    1269        private static int[] inverse(final int[] indices) {
    1270            final int[] inverse = new int[indices.length];
    1271            for (int i = 0; i < indices.length; i++) {
    1272                inverse[indices[i]] = i;
    1273            }
    1274            return inverse;
    1275        }
    1276
    1277        /**
    1278         * @param indices Input index array.
    1279         * @return the indices in inverse order (last is first).
    1280         */
    1281        private static int[] reverse(final int[] indices) {
    1282            final int[] reverse = new int[indices.length];
    1283            for (int i = 0; i < indices.length; i++) {
    1284                reverse[i] = indices[indices.length - i - 1];
    1285            }
    1286            return reverse;
    1287        }
    1288
    1289        /**
    1290         * @param size Length of random array.
    1291         * @return an array of Gaussian random numbers.
    1292         */
    1293        private double[] randn(int size) {
    1294            final double[] randn = new double[size];
    1295            for (int i = 0; i < size; i++) {
    1296                randn[i] = random.nextGaussian();
    1297            }
    1298            return randn;
    1299        }
    1300
    1301        /**
    1302         * @param size Number of rows.
    1303         * @param popSize Population size.
    1304         * @return a 2-dimensional matrix of Gaussian random numbers.
    1305         */
    1306        private RealMatrix randn1(int size, int popSize) {
    1307            final double[][] d = new double[size][popSize];
    1308            for (int r = 0; r < size; r++) {
    1309                for (int c = 0; c < popSize; c++) {
    1310                    d[r][c] = random.nextGaussian();
    1311                }
    1312            }
    1313            return new Array2DRowRealMatrix(d, false);
    1314        }
    1315    }