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