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