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    }