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