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.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 1499808 2013-07-04 17:00:42Z sebb $
081 * @deprecated As of 3.1 (to be removed in 4.0).
082 * @since 3.0
083 */
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 *= Math.exp(Math.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 * Math.max(Math.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                Math.max(historyWorst, worstFitness) -
599                Math.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 = sigma * Math.exp(0.2 + cs / damps);
624            }
625            if (iterations > 2 && Math.max(historyWorst, bestFitness) -
626                Math.min(historyBest, bestFitness) == 0) {
627                sigma = sigma * Math.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 * Math.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 = Math.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 * Math.max(0, Math.sqrt((mueff - 1) /
737                                               (dimension + 1)) - 1)) *
738            Math.max(0.3,
739                     1 - dimension / (1e-6 + maxIterations)) + cs; // minor increment
740        ccov1 = 2 / ((dimension + 1.3) * (dimension + 1.3) + mueff);
741        ccovmu = Math.min(1 - ccov1, 2 * (mueff - 2 + 1 / mueff) /
742                          ((dimension + 2) * (dimension + 2) + mueff));
743        ccov1Sep = Math.min(1, ccov1 * (dimension + 1.5) / 3);
744        ccovmuSep = Math.min(1 - ccov1, ccovmu * (dimension + 1.5) / 3);
745        chiN = Math.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(
777                        Math.sqrt(cs * (2 - cs) * mueff)));
778        normps = ps.getFrobeniusNorm();
779        final boolean hsig = normps /
780            Math.sqrt(1 - Math.pow(1 - cs, 2 * iterations)) /
781            chiN < 1.4 + 2 / ((double) dimension + 1);
782        pc = pc.scalarMultiply(1 - cc);
783        if (hsig) {
784            pc = pc.add(xmean.subtract(xold).scalarMultiply(Math.sqrt(cc * (2 - cc) * mueff) / sigma));
785        }
786        return hsig;
787    }
788
789    /**
790     * Update of the covariance matrix C for diagonalOnly > 0
791     *
792     * @param hsig Flag indicating a small correction.
793     * @param bestArz Fitness-sorted matrix of the gaussian random values of the
794     * current offspring.
795     */
796    private void updateCovarianceDiagonalOnly(boolean hsig,
797                                              final RealMatrix bestArz) {
798        // minor correction if hsig==false
799        double oldFac = hsig ? 0 : ccov1Sep * cc * (2 - cc);
800        oldFac += 1 - ccov1Sep - ccovmuSep;
801        diagC = diagC.scalarMultiply(oldFac) // regard old matrix
802            .add(square(pc).scalarMultiply(ccov1Sep)) // plus rank one update
803            .add((times(diagC, square(bestArz).multiply(weights))) // plus rank mu update
804                 .scalarMultiply(ccovmuSep));
805        diagD = sqrt(diagC); // replaces eig(C)
806        if (diagonalOnly > 1 &&
807            iterations > diagonalOnly) {
808            // full covariance matrix from now on
809            diagonalOnly = 0;
810            B = eye(dimension, dimension);
811            BD = diag(diagD);
812            C = diag(diagC);
813        }
814    }
815
816    /**
817     * Update of the covariance matrix C.
818     *
819     * @param hsig Flag indicating a small correction.
820     * @param bestArx Fitness-sorted matrix of the argument vectors producing the
821     * current offspring.
822     * @param arz Unsorted matrix containing the gaussian random values of the
823     * current offspring.
824     * @param arindex Indices indicating the fitness-order of the current offspring.
825     * @param xold xmean matrix of the previous generation.
826     */
827    private void updateCovariance(boolean hsig, final RealMatrix bestArx,
828                                  final RealMatrix arz, final int[] arindex,
829                                  final RealMatrix xold) {
830        double negccov = 0;
831        if (ccov1 + ccovmu > 0) {
832            final RealMatrix arpos = bestArx.subtract(repmat(xold, 1, mu))
833                .scalarMultiply(1 / sigma); // mu difference vectors
834            final RealMatrix roneu = pc.multiply(pc.transpose())
835                .scalarMultiply(ccov1); // rank one update
836            // minor correction if hsig==false
837            double oldFac = hsig ? 0 : ccov1 * cc * (2 - cc);
838            oldFac += 1 - ccov1 - ccovmu;
839            if (isActiveCMA) {
840                // Adapt covariance matrix C active CMA
841                negccov = (1 - ccovmu) * 0.25 * mueff /
842                    (Math.pow(dimension + 2, 1.5) + 2 * mueff);
843                // keep at least 0.66 in all directions, small popsize are most
844                // critical
845                final double negminresidualvariance = 0.66;
846                // where to make up for the variance loss
847                final double negalphaold = 0.5;
848                // prepare vectors, compute negative updating matrix Cneg
849                final int[] arReverseIndex = reverse(arindex);
850                RealMatrix arzneg = selectColumns(arz, MathArrays.copyOf(arReverseIndex, mu));
851                RealMatrix arnorms = sqrt(sumRows(square(arzneg)));
852                final int[] idxnorms = sortedIndices(arnorms.getRow(0));
853                final RealMatrix arnormsSorted = selectColumns(arnorms, idxnorms);
854                final int[] idxReverse = reverse(idxnorms);
855                final RealMatrix arnormsReverse = selectColumns(arnorms, idxReverse);
856                arnorms = divide(arnormsReverse, arnormsSorted);
857                final int[] idxInv = inverse(idxnorms);
858                final RealMatrix arnormsInv = selectColumns(arnorms, idxInv);
859                // check and set learning rate negccov
860                final double negcovMax = (1 - negminresidualvariance) /
861                    square(arnormsInv).multiply(weights).getEntry(0, 0);
862                if (negccov > negcovMax) {
863                    negccov = negcovMax;
864                }
865                arzneg = times(arzneg, repmat(arnormsInv, dimension, 1));
866                final RealMatrix artmp = BD.multiply(arzneg);
867                final RealMatrix Cneg = artmp.multiply(diag(weights)).multiply(artmp.transpose());
868                oldFac += negalphaold * negccov;
869                C = C.scalarMultiply(oldFac)
870                    .add(roneu) // regard old matrix
871                    .add(arpos.scalarMultiply( // plus rank one update
872                                              ccovmu + (1 - negalphaold) * negccov) // plus rank mu update
873                         .multiply(times(repmat(weights, 1, dimension),
874                                         arpos.transpose())))
875                    .subtract(Cneg.scalarMultiply(negccov));
876            } else {
877                // Adapt covariance matrix C - nonactive
878                C = C.scalarMultiply(oldFac) // regard old matrix
879                    .add(roneu) // plus rank one update
880                    .add(arpos.scalarMultiply(ccovmu) // plus rank mu update
881                         .multiply(times(repmat(weights, 1, dimension),
882                                         arpos.transpose())));
883            }
884        }
885        updateBD(negccov);
886    }
887
888    /**
889     * Update B and D from C.
890     *
891     * @param negccov Negative covariance factor.
892     */
893    private void updateBD(double negccov) {
894        if (ccov1 + ccovmu + negccov > 0 &&
895            (iterations % 1. / (ccov1 + ccovmu + negccov) / dimension / 10.) < 1) {
896            // to achieve O(N^2)
897            C = triu(C, 0).add(triu(C, 1).transpose());
898            // enforce symmetry to prevent complex numbers
899            final EigenDecomposition eig = new EigenDecomposition(C);
900            B = eig.getV(); // eigen decomposition, B==normalized eigenvectors
901            D = eig.getD();
902            diagD = diag(D);
903            if (min(diagD) <= 0) {
904                for (int i = 0; i < dimension; i++) {
905                    if (diagD.getEntry(i, 0) < 0) {
906                        diagD.setEntry(i, 0, 0);
907                    }
908                }
909                final double tfac = max(diagD) / 1e14;
910                C = C.add(eye(dimension, dimension).scalarMultiply(tfac));
911                diagD = diagD.add(ones(dimension, 1).scalarMultiply(tfac));
912            }
913            if (max(diagD) > 1e14 * min(diagD)) {
914                final double tfac = max(diagD) / 1e14 - min(diagD);
915                C = C.add(eye(dimension, dimension).scalarMultiply(tfac));
916                diagD = diagD.add(ones(dimension, 1).scalarMultiply(tfac));
917            }
918            diagC = diag(C);
919            diagD = sqrt(diagD); // D contains standard deviations now
920            BD = times(B, repmat(diagD.transpose(), dimension, 1)); // O(n^2)
921        }
922    }
923
924    /**
925     * Pushes the current best fitness value in a history queue.
926     *
927     * @param vals History queue.
928     * @param val Current best fitness value.
929     */
930    private static void push(double[] vals, double val) {
931        for (int i = vals.length-1; i > 0; i--) {
932            vals[i] = vals[i-1];
933        }
934        vals[0] = val;
935    }
936
937    /**
938     * Sorts fitness values.
939     *
940     * @param doubles Array of values to be sorted.
941     * @return a sorted array of indices pointing into doubles.
942     */
943    private int[] sortedIndices(final double[] doubles) {
944        final DoubleIndex[] dis = new DoubleIndex[doubles.length];
945        for (int i = 0; i < doubles.length; i++) {
946            dis[i] = new DoubleIndex(doubles[i], i);
947        }
948        Arrays.sort(dis);
949        final int[] indices = new int[doubles.length];
950        for (int i = 0; i < doubles.length; i++) {
951            indices[i] = dis[i].index;
952        }
953        return indices;
954    }
955
956    /**
957     * Used to sort fitness values. Sorting is always in lower value first
958     * order.
959     */
960    private static class DoubleIndex implements Comparable<DoubleIndex> {
961        /** Value to compare. */
962        private final double value;
963        /** Index into sorted array. */
964        private final int index;
965
966        /**
967         * @param value Value to compare.
968         * @param index Index into sorted array.
969         */
970        DoubleIndex(double value, int index) {
971            this.value = value;
972            this.index = index;
973        }
974
975        /** {@inheritDoc} */
976        public int compareTo(DoubleIndex o) {
977            return Double.compare(value, o.value);
978        }
979
980        /** {@inheritDoc} */
981        @Override
982        public boolean equals(Object other) {
983
984            if (this == other) {
985                return true;
986            }
987
988            if (other instanceof DoubleIndex) {
989                return Double.compare(value, ((DoubleIndex) other).value) == 0;
990            }
991
992            return false;
993        }
994
995        /** {@inheritDoc} */
996        @Override
997        public int hashCode() {
998            long bits = Double.doubleToLongBits(value);
999            return (int) ((1438542 ^ (bits >>> 32) ^ bits) & 0xffffffff);
1000        }
1001    }
1002
1003    /**
1004     * Normalizes fitness values to the range [0,1]. Adds a penalty to the
1005     * fitness value if out of range. The penalty is adjusted by calling
1006     * setValueRange().
1007     */
1008    private class FitnessFunction {
1009        /** Determines the penalty for boundary violations */
1010        private double valueRange;
1011        /**
1012         * Flag indicating whether the objective variables are forced into their
1013         * bounds if defined
1014         */
1015        private final boolean isRepairMode;
1016
1017        /** Simple constructor.
1018         */
1019        public FitnessFunction() {
1020            valueRange = 1;
1021            isRepairMode = true;
1022        }
1023
1024        /**
1025         * @param point Normalized objective variables.
1026         * @return the objective value + penalty for violated bounds.
1027         */
1028        public double value(final double[] point) {
1029            double value;
1030            if (isRepairMode) {
1031                double[] repaired = repair(point);
1032                value = CMAESOptimizer.this.computeObjectiveValue(repaired) +
1033                    penalty(point, repaired);
1034            } else {
1035                value = CMAESOptimizer.this.computeObjectiveValue(point);
1036            }
1037            return isMinimize ? value : -value;
1038        }
1039
1040        /**
1041         * @param x Normalized objective variables.
1042         * @return {@code true} if in bounds.
1043         */
1044        public boolean isFeasible(final double[] x) {
1045            final double[] lB = CMAESOptimizer.this.getLowerBound();
1046            final double[] uB = CMAESOptimizer.this.getUpperBound();
1047
1048            for (int i = 0; i < x.length; i++) {
1049                if (x[i] < lB[i]) {
1050                    return false;
1051                }
1052                if (x[i] > uB[i]) {
1053                    return false;
1054                }
1055            }
1056            return true;
1057        }
1058
1059        /**
1060         * @param valueRange Adjusts the penalty computation.
1061         */
1062        public void setValueRange(double valueRange) {
1063            this.valueRange = valueRange;
1064        }
1065
1066        /**
1067         * @param x Normalized objective variables.
1068         * @return the repaired (i.e. all in bounds) objective variables.
1069         */
1070        private double[] repair(final double[] x) {
1071            final double[] lB = CMAESOptimizer.this.getLowerBound();
1072            final double[] uB = CMAESOptimizer.this.getUpperBound();
1073
1074            final double[] repaired = new double[x.length];
1075            for (int i = 0; i < x.length; i++) {
1076                if (x[i] < lB[i]) {
1077                    repaired[i] = lB[i];
1078                } else if (x[i] > uB[i]) {
1079                    repaired[i] = uB[i];
1080                } else {
1081                    repaired[i] = x[i];
1082                }
1083            }
1084            return repaired;
1085        }
1086
1087        /**
1088         * @param x Normalized objective variables.
1089         * @param repaired Repaired objective variables.
1090         * @return Penalty value according to the violation of the bounds.
1091         */
1092        private double penalty(final double[] x, final double[] repaired) {
1093            double penalty = 0;
1094            for (int i = 0; i < x.length; i++) {
1095                double diff = Math.abs(x[i] - repaired[i]);
1096                penalty += diff * valueRange;
1097            }
1098            return isMinimize ? penalty : -penalty;
1099        }
1100    }
1101
1102    // -----Matrix utility functions similar to the Matlab build in functions------
1103
1104    /**
1105     * @param m Input matrix
1106     * @return Matrix representing the element-wise logarithm of m.
1107     */
1108    private static RealMatrix log(final RealMatrix m) {
1109        final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
1110        for (int r = 0; r < m.getRowDimension(); r++) {
1111            for (int c = 0; c < m.getColumnDimension(); c++) {
1112                d[r][c] = Math.log(m.getEntry(r, c));
1113            }
1114        }
1115        return new Array2DRowRealMatrix(d, false);
1116    }
1117
1118    /**
1119     * @param m Input matrix.
1120     * @return Matrix representing the element-wise square root of m.
1121     */
1122    private static RealMatrix sqrt(final RealMatrix m) {
1123        final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
1124        for (int r = 0; r < m.getRowDimension(); r++) {
1125            for (int c = 0; c < m.getColumnDimension(); c++) {
1126                d[r][c] = Math.sqrt(m.getEntry(r, c));
1127            }
1128        }
1129        return new Array2DRowRealMatrix(d, false);
1130    }
1131
1132    /**
1133     * @param m Input matrix.
1134     * @return Matrix representing the element-wise square of m.
1135     */
1136    private static RealMatrix square(final RealMatrix m) {
1137        final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
1138        for (int r = 0; r < m.getRowDimension(); r++) {
1139            for (int c = 0; c < m.getColumnDimension(); c++) {
1140                double e = m.getEntry(r, c);
1141                d[r][c] = e * e;
1142            }
1143        }
1144        return new Array2DRowRealMatrix(d, false);
1145    }
1146
1147    /**
1148     * @param m Input matrix 1.
1149     * @param n Input matrix 2.
1150     * @return the matrix where the elements of m and n are element-wise multiplied.
1151     */
1152    private static RealMatrix times(final RealMatrix m, final RealMatrix n) {
1153        final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
1154        for (int r = 0; r < m.getRowDimension(); r++) {
1155            for (int c = 0; c < m.getColumnDimension(); c++) {
1156                d[r][c] = m.getEntry(r, c) * n.getEntry(r, c);
1157            }
1158        }
1159        return new Array2DRowRealMatrix(d, false);
1160    }
1161
1162    /**
1163     * @param m Input matrix 1.
1164     * @param n Input matrix 2.
1165     * @return Matrix where the elements of m and n are element-wise divided.
1166     */
1167    private static RealMatrix divide(final RealMatrix m, final RealMatrix n) {
1168        final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
1169        for (int r = 0; r < m.getRowDimension(); r++) {
1170            for (int c = 0; c < m.getColumnDimension(); c++) {
1171                d[r][c] = m.getEntry(r, c) / n.getEntry(r, c);
1172            }
1173        }
1174        return new Array2DRowRealMatrix(d, false);
1175    }
1176
1177    /**
1178     * @param m Input matrix.
1179     * @param cols Columns to select.
1180     * @return Matrix representing the selected columns.
1181     */
1182    private static RealMatrix selectColumns(final RealMatrix m, final int[] cols) {
1183        final double[][] d = new double[m.getRowDimension()][cols.length];
1184        for (int r = 0; r < m.getRowDimension(); r++) {
1185            for (int c = 0; c < cols.length; c++) {
1186                d[r][c] = m.getEntry(r, cols[c]);
1187            }
1188        }
1189        return new Array2DRowRealMatrix(d, false);
1190    }
1191
1192    /**
1193     * @param m Input matrix.
1194     * @param k Diagonal position.
1195     * @return Upper triangular part of matrix.
1196     */
1197    private static RealMatrix triu(final RealMatrix m, int k) {
1198        final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
1199        for (int r = 0; r < m.getRowDimension(); r++) {
1200            for (int c = 0; c < m.getColumnDimension(); c++) {
1201                d[r][c] = r <= c - k ? m.getEntry(r, c) : 0;
1202            }
1203        }
1204        return new Array2DRowRealMatrix(d, false);
1205    }
1206
1207    /**
1208     * @param m Input matrix.
1209     * @return Row matrix representing the sums of the rows.
1210     */
1211    private static RealMatrix sumRows(final RealMatrix m) {
1212        final double[][] d = new double[1][m.getColumnDimension()];
1213        for (int c = 0; c < m.getColumnDimension(); c++) {
1214            double sum = 0;
1215            for (int r = 0; r < m.getRowDimension(); r++) {
1216                sum += m.getEntry(r, c);
1217            }
1218            d[0][c] = sum;
1219        }
1220        return new Array2DRowRealMatrix(d, false);
1221    }
1222
1223    /**
1224     * @param m Input matrix.
1225     * @return the diagonal n-by-n matrix if m is a column matrix or the column
1226     * matrix representing the diagonal if m is a n-by-n matrix.
1227     */
1228    private static RealMatrix diag(final RealMatrix m) {
1229        if (m.getColumnDimension() == 1) {
1230            final double[][] d = new double[m.getRowDimension()][m.getRowDimension()];
1231            for (int i = 0; i < m.getRowDimension(); i++) {
1232                d[i][i] = m.getEntry(i, 0);
1233            }
1234            return new Array2DRowRealMatrix(d, false);
1235        } else {
1236            final double[][] d = new double[m.getRowDimension()][1];
1237            for (int i = 0; i < m.getColumnDimension(); i++) {
1238                d[i][0] = m.getEntry(i, i);
1239            }
1240            return new Array2DRowRealMatrix(d, false);
1241        }
1242    }
1243
1244    /**
1245     * Copies a column from m1 to m2.
1246     *
1247     * @param m1 Source matrix.
1248     * @param col1 Source column.
1249     * @param m2 Target matrix.
1250     * @param col2 Target column.
1251     */
1252    private static void copyColumn(final RealMatrix m1, int col1,
1253                                   RealMatrix m2, int col2) {
1254        for (int i = 0; i < m1.getRowDimension(); i++) {
1255            m2.setEntry(i, col2, m1.getEntry(i, col1));
1256        }
1257    }
1258
1259    /**
1260     * @param n Number of rows.
1261     * @param m Number of columns.
1262     * @return n-by-m matrix filled with 1.
1263     */
1264    private static RealMatrix ones(int n, int m) {
1265        final double[][] d = new double[n][m];
1266        for (int r = 0; r < n; r++) {
1267            Arrays.fill(d[r], 1);
1268        }
1269        return new Array2DRowRealMatrix(d, false);
1270    }
1271
1272    /**
1273     * @param n Number of rows.
1274     * @param m Number of columns.
1275     * @return n-by-m matrix of 0 values out of diagonal, and 1 values on
1276     * the diagonal.
1277     */
1278    private static RealMatrix eye(int n, int m) {
1279        final double[][] d = new double[n][m];
1280        for (int r = 0; r < n; r++) {
1281            if (r < m) {
1282                d[r][r] = 1;
1283            }
1284        }
1285        return new Array2DRowRealMatrix(d, false);
1286    }
1287
1288    /**
1289     * @param n Number of rows.
1290     * @param m Number of columns.
1291     * @return n-by-m matrix of zero values.
1292     */
1293    private static RealMatrix zeros(int n, int m) {
1294        return new Array2DRowRealMatrix(n, m);
1295    }
1296
1297    /**
1298     * @param mat Input matrix.
1299     * @param n Number of row replicates.
1300     * @param m Number of column replicates.
1301     * @return a matrix which replicates the input matrix in both directions.
1302     */
1303    private static RealMatrix repmat(final RealMatrix mat, int n, int m) {
1304        final int rd = mat.getRowDimension();
1305        final int cd = mat.getColumnDimension();
1306        final double[][] d = new double[n * rd][m * cd];
1307        for (int r = 0; r < n * rd; r++) {
1308            for (int c = 0; c < m * cd; c++) {
1309                d[r][c] = mat.getEntry(r % rd, c % cd);
1310            }
1311        }
1312        return new Array2DRowRealMatrix(d, false);
1313    }
1314
1315    /**
1316     * @param start Start value.
1317     * @param end End value.
1318     * @param step Step size.
1319     * @return a sequence as column matrix.
1320     */
1321    private static RealMatrix sequence(double start, double end, double step) {
1322        final int size = (int) ((end - start) / step + 1);
1323        final double[][] d = new double[size][1];
1324        double value = start;
1325        for (int r = 0; r < size; r++) {
1326            d[r][0] = value;
1327            value += step;
1328        }
1329        return new Array2DRowRealMatrix(d, false);
1330    }
1331
1332    /**
1333     * @param m Input matrix.
1334     * @return the maximum of the matrix element values.
1335     */
1336    private static double max(final RealMatrix m) {
1337        double max = -Double.MAX_VALUE;
1338        for (int r = 0; r < m.getRowDimension(); r++) {
1339            for (int c = 0; c < m.getColumnDimension(); c++) {
1340                double e = m.getEntry(r, c);
1341                if (max < e) {
1342                    max = e;
1343                }
1344            }
1345        }
1346        return max;
1347    }
1348
1349    /**
1350     * @param m Input matrix.
1351     * @return the minimum of the matrix element values.
1352     */
1353    private static double min(final RealMatrix m) {
1354        double min = Double.MAX_VALUE;
1355        for (int r = 0; r < m.getRowDimension(); r++) {
1356            for (int c = 0; c < m.getColumnDimension(); c++) {
1357                double e = m.getEntry(r, c);
1358                if (min > e) {
1359                    min = e;
1360                }
1361            }
1362        }
1363        return min;
1364    }
1365
1366    /**
1367     * @param m Input array.
1368     * @return the maximum of the array values.
1369     */
1370    private static double max(final double[] m) {
1371        double max = -Double.MAX_VALUE;
1372        for (int r = 0; r < m.length; r++) {
1373            if (max < m[r]) {
1374                max = m[r];
1375            }
1376        }
1377        return max;
1378    }
1379
1380    /**
1381     * @param m Input array.
1382     * @return the minimum of the array values.
1383     */
1384    private static double min(final double[] m) {
1385        double min = Double.MAX_VALUE;
1386        for (int r = 0; r < m.length; r++) {
1387            if (min > m[r]) {
1388                min = m[r];
1389            }
1390        }
1391        return min;
1392    }
1393
1394    /**
1395     * @param indices Input index array.
1396     * @return the inverse of the mapping defined by indices.
1397     */
1398    private static int[] inverse(final int[] indices) {
1399        final int[] inverse = new int[indices.length];
1400        for (int i = 0; i < indices.length; i++) {
1401            inverse[indices[i]] = i;
1402        }
1403        return inverse;
1404    }
1405
1406    /**
1407     * @param indices Input index array.
1408     * @return the indices in inverse order (last is first).
1409     */
1410    private static int[] reverse(final int[] indices) {
1411        final int[] reverse = new int[indices.length];
1412        for (int i = 0; i < indices.length; i++) {
1413            reverse[i] = indices[indices.length - i - 1];
1414        }
1415        return reverse;
1416    }
1417
1418    /**
1419     * @param size Length of random array.
1420     * @return an array of Gaussian random numbers.
1421     */
1422    private double[] randn(int size) {
1423        final double[] randn = new double[size];
1424        for (int i = 0; i < size; i++) {
1425            randn[i] = random.nextGaussian();
1426        }
1427        return randn;
1428    }
1429
1430    /**
1431     * @param size Number of rows.
1432     * @param popSize Population size.
1433     * @return a 2-dimensional matrix of Gaussian random numbers.
1434     */
1435    private RealMatrix randn1(int size, int popSize) {
1436        final double[][] d = new double[size][popSize];
1437        for (int r = 0; r < size; r++) {
1438            for (int c = 0; c < popSize; c++) {
1439                d[r][c] = random.nextGaussian();
1440            }
1441        }
1442        return new Array2DRowRealMatrix(d, false);
1443    }
1444}