View Javadoc

1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  package org.apache.commons.math3.optimization.direct;
19  
20  import java.util.ArrayList;
21  import java.util.Arrays;
22  import java.util.List;
23  
24  import org.apache.commons.math3.analysis.MultivariateFunction;
25  import org.apache.commons.math3.exception.DimensionMismatchException;
26  import org.apache.commons.math3.exception.NotPositiveException;
27  import org.apache.commons.math3.exception.NotStrictlyPositiveException;
28  import org.apache.commons.math3.exception.OutOfRangeException;
29  import org.apache.commons.math3.exception.TooManyEvaluationsException;
30  import org.apache.commons.math3.linear.Array2DRowRealMatrix;
31  import org.apache.commons.math3.linear.EigenDecomposition;
32  import org.apache.commons.math3.linear.MatrixUtils;
33  import org.apache.commons.math3.linear.RealMatrix;
34  import org.apache.commons.math3.optimization.ConvergenceChecker;
35  import org.apache.commons.math3.optimization.OptimizationData;
36  import org.apache.commons.math3.optimization.GoalType;
37  import org.apache.commons.math3.optimization.MultivariateOptimizer;
38  import org.apache.commons.math3.optimization.PointValuePair;
39  import org.apache.commons.math3.optimization.SimpleValueChecker;
40  import org.apache.commons.math3.random.MersenneTwister;
41  import org.apache.commons.math3.random.RandomGenerator;
42  import org.apache.commons.math3.util.MathArrays;
43  
44  /**
45   * <p>An implementation of the active Covariance Matrix Adaptation Evolution Strategy (CMA-ES)
46   * for non-linear, non-convex, non-smooth, global function minimization.
47   * The CMA-Evolution Strategy (CMA-ES) is a reliable stochastic optimization method
48   * which should be applied if derivative-based methods, e.g. quasi-Newton BFGS or
49   * conjugate gradient, fail due to a rugged search landscape (e.g. noise, local
50   * optima, outlier, etc.) of the objective function. Like a
51   * quasi-Newton method, the CMA-ES learns and applies a variable metric
52   * on the underlying search space. Unlike a quasi-Newton method, the
53   * CMA-ES neither estimates nor uses gradients, making it considerably more
54   * reliable in terms of finding a good, or even close to optimal, solution.</p>
55   *
56   * <p>In general, on smooth objective functions the CMA-ES is roughly ten times
57   * slower than BFGS (counting objective function evaluations, no gradients provided).
58   * For up to <math>N=10</math> variables also the derivative-free simplex
59   * direct search method (Nelder and Mead) can be faster, but it is
60   * far less reliable than CMA-ES.</p>
61   *
62   * <p>The CMA-ES is particularly well suited for non-separable
63   * and/or badly conditioned problems. To observe the advantage of CMA compared
64   * to a conventional evolution strategy, it will usually take about
65   * <math>30 N</math> function evaluations. On difficult problems the complete
66   * optimization (a single run) is expected to take <em>roughly</em> between
67   * <math>30 N</math> and <math>300 N<sup>2</sup></math>
68   * function evaluations.</p>
69   *
70   * <p>This implementation is translated and adapted from the Matlab version
71   * of the CMA-ES algorithm as implemented in module {@code cmaes.m} version 3.51.</p>
72   *
73   * For more information, please refer to the following links:
74   * <ul>
75   *  <li><a href="http://www.lri.fr/~hansen/cmaes.m">Matlab code</a></li>
76   *  <li><a href="http://www.lri.fr/~hansen/cmaesintro.html">Introduction to CMA-ES</a></li>
77   *  <li><a href="http://en.wikipedia.org/wiki/CMA-ES">Wikipedia</a></li>
78   * </ul>
79   *
80   * @version $Id: CMAESOptimizer.java 1499808 2013-07-04 17:00:42Z sebb $
81   * @deprecated As of 3.1 (to be removed in 4.0).
82   * @since 3.0
83   */
84  
85  @Deprecated
86  public class CMAESOptimizer
87      extends BaseAbstractMultivariateSimpleBoundsOptimizer<MultivariateFunction>
88      implements MultivariateOptimizer {
89      /** Default value for {@link #checkFeasableCount}: {@value}. */
90      public static final int DEFAULT_CHECKFEASABLECOUNT = 0;
91      /** Default value for {@link #stopFitness}: {@value}. */
92      public static final double DEFAULT_STOPFITNESS = 0;
93      /** Default value for {@link #isActiveCMA}: {@value}. */
94      public static final boolean DEFAULT_ISACTIVECMA = true;
95      /** Default value for {@link #maxIterations}: {@value}. */
96      public static final int DEFAULT_MAXITERATIONS = 30000;
97      /** Default value for {@link #diagonalOnly}: {@value}. */
98      public static final int DEFAULT_DIAGONALONLY = 0;
99      /** 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 }