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.FastMath;
43  import org.apache.commons.math3.util.MathArrays;
44  
45  /**
46   * <p>An implementation of the active Covariance Matrix Adaptation Evolution Strategy (CMA-ES)
47   * for non-linear, non-convex, non-smooth, global function minimization.
48   * The CMA-Evolution Strategy (CMA-ES) is a reliable stochastic optimization method
49   * which should be applied if derivative-based methods, e.g. quasi-Newton BFGS or
50   * conjugate gradient, fail due to a rugged search landscape (e.g. noise, local
51   * optima, outlier, etc.) of the objective function. Like a
52   * quasi-Newton method, the CMA-ES learns and applies a variable metric
53   * on the underlying search space. Unlike a quasi-Newton method, the
54   * CMA-ES neither estimates nor uses gradients, making it considerably more
55   * reliable in terms of finding a good, or even close to optimal, solution.</p>
56   *
57   * <p>In general, on smooth objective functions the CMA-ES is roughly ten times
58   * slower than BFGS (counting objective function evaluations, no gradients provided).
59   * For up to <math>N=10</math> variables also the derivative-free simplex
60   * direct search method (Nelder and Mead) can be faster, but it is
61   * far less reliable than CMA-ES.</p>
62   *
63   * <p>The CMA-ES is particularly well suited for non-separable
64   * and/or badly conditioned problems. To observe the advantage of CMA compared
65   * to a conventional evolution strategy, it will usually take about
66   * <math>30 N</math> function evaluations. On difficult problems the complete
67   * optimization (a single run) is expected to take <em>roughly</em> between
68   * <math>30 N</math> and <math>300 N<sup>2</sup></math>
69   * function evaluations.</p>
70   *
71   * <p>This implementation is translated and adapted from the Matlab version
72   * of the CMA-ES algorithm as implemented in module {@code cmaes.m} version 3.51.</p>
73   *
74   * For more information, please refer to the following links:
75   * <ul>
76   *  <li><a href="http://www.lri.fr/~hansen/cmaes.m">Matlab code</a></li>
77   *  <li><a href="http://www.lri.fr/~hansen/cmaesintro.html">Introduction to CMA-ES</a></li>
78   *  <li><a href="http://en.wikipedia.org/wiki/CMA-ES">Wikipedia</a></li>
79   * </ul>
80   *
81   * @version $Id: CMAESOptimizer.java 1540165 2013-11-08 19:53:58Z tn $
82   * @deprecated As of 3.1 (to be removed in 4.0).
83   * @since 3.0
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 *= FastMath.exp(FastMath.min(1, (normps/chiN - 1) * cs / damps));
563             final double bestFitness = fitness[arindex[0]];
564             final double worstFitness = fitness[arindex[arindex.length - 1]];
565             if (bestValue > bestFitness) {
566                 bestValue = bestFitness;
567                 lastResult = optimum;
568                 optimum = new PointValuePair(fitfun.repair(bestArx.getColumn(0)),
569                                              isMinimize ? bestFitness : -bestFitness);
570                 if (getConvergenceChecker() != null && lastResult != null &&
571                     getConvergenceChecker().converged(iterations, optimum, lastResult)) {
572                     break generationLoop;
573                 }
574             }
575             // handle termination criteria
576             // Break, if fitness is good enough
577             if (stopFitness != 0 && bestFitness < (isMinimize ? stopFitness : -stopFitness)) {
578                 break generationLoop;
579             }
580             final double[] sqrtDiagC = sqrt(diagC).getColumn(0);
581             final double[] pcCol = pc.getColumn(0);
582             for (int i = 0; i < dimension; i++) {
583                 if (sigma * FastMath.max(FastMath.abs(pcCol[i]), sqrtDiagC[i]) > stopTolX) {
584                     break;
585                 }
586                 if (i >= dimension - 1) {
587                     break generationLoop;
588                 }
589             }
590             for (int i = 0; i < dimension; i++) {
591                 if (sigma * sqrtDiagC[i] > stopTolUpX) {
592                     break generationLoop;
593                 }
594             }
595             final double historyBest = min(fitnessHistory);
596             final double historyWorst = max(fitnessHistory);
597             if (iterations > 2 &&
598                 FastMath.max(historyWorst, worstFitness) -
599                 FastMath.min(historyBest, bestFitness) < stopTolFun) {
600                 break generationLoop;
601             }
602             if (iterations > fitnessHistory.length &&
603                 historyWorst-historyBest < stopTolHistFun) {
604                 break generationLoop;
605             }
606             // condition number of the covariance matrix exceeds 1e14
607             if (max(diagD)/min(diagD) > 1e7) {
608                 break generationLoop;
609             }
610             // user defined termination
611             if (getConvergenceChecker() != null) {
612                 final PointValuePair current
613                     = new PointValuePair(bestArx.getColumn(0),
614                                          isMinimize ? bestFitness : -bestFitness);
615                 if (lastResult != null &&
616                     getConvergenceChecker().converged(iterations, current, lastResult)) {
617                     break generationLoop;
618                     }
619                 lastResult = current;
620             }
621             // Adjust step size in case of equal function values (flat fitness)
622             if (bestValue == fitness[arindex[(int)(0.1+lambda/4.)]]) {
623                 sigma *= FastMath.exp(0.2 + cs / damps);
624             }
625             if (iterations > 2 && FastMath.max(historyWorst, bestFitness) -
626                 FastMath.min(historyBest, bestFitness) == 0) {
627                 sigma *= FastMath.exp(0.2 + cs / damps);
628             }
629             // store best in history
630             push(fitnessHistory,bestFitness);
631             fitfun.setValueRange(worstFitness-bestFitness);
632             if (generateStatistics) {
633                 statisticsSigmaHistory.add(sigma);
634                 statisticsFitnessHistory.add(bestFitness);
635                 statisticsMeanHistory.add(xmean.transpose());
636                 statisticsDHistory.add(diagD.transpose().scalarMultiply(1E5));
637             }
638         }
639         return optimum;
640     }
641 
642     /**
643      * Scans the list of (required and optional) optimization data that
644      * characterize the problem.
645      *
646      * @param optData Optimization data. The following data will be looked for:
647      * <ul>
648      *  <li>{@link Sigma}</li>
649      *  <li>{@link PopulationSize}</li>
650      * </ul>
651      */
652     private void parseOptimizationData(OptimizationData... optData) {
653         // The existing values (as set by the previous call) are reused if
654         // not provided in the argument list.
655         for (OptimizationData data : optData) {
656             if (data instanceof Sigma) {
657                 inputSigma = ((Sigma) data).getSigma();
658                 continue;
659             }
660             if (data instanceof PopulationSize) {
661                 lambda = ((PopulationSize) data).getPopulationSize();
662                 continue;
663             }
664         }
665     }
666 
667     /**
668      * Checks dimensions and values of boundaries and inputSigma if defined.
669      */
670     private void checkParameters() {
671         final double[] init = getStartPoint();
672         final double[] lB = getLowerBound();
673         final double[] uB = getUpperBound();
674 
675         if (inputSigma != null) {
676             if (inputSigma.length != init.length) {
677                 throw new DimensionMismatchException(inputSigma.length, init.length);
678             }
679             for (int i = 0; i < init.length; i++) {
680                 if (inputSigma[i] < 0) {
681                     // XXX Remove this block in 4.0 (check performed in "Sigma" class).
682                     throw new NotPositiveException(inputSigma[i]);
683                 }
684                 if (inputSigma[i] > uB[i] - lB[i]) {
685                     throw new OutOfRangeException(inputSigma[i], 0, uB[i] - lB[i]);
686                 }
687             }
688         }
689     }
690 
691     /**
692      * Initialization of the dynamic search parameters
693      *
694      * @param guess Initial guess for the arguments of the fitness function.
695      */
696     private void initializeCMA(double[] guess) {
697         if (lambda <= 0) {
698             // XXX Line below to replace the current one in 4.0 (MATH-879).
699             // throw new NotStrictlyPositiveException(lambda);
700             lambda = 4 + (int) (3 * FastMath.log(dimension));
701         }
702         // initialize sigma
703         final double[][] sigmaArray = new double[guess.length][1];
704         for (int i = 0; i < guess.length; i++) {
705             // XXX Line below to replace the current one in 4.0 (MATH-868).
706             // sigmaArray[i][0] = inputSigma[i];
707             sigmaArray[i][0] = inputSigma == null ? 0.3 : inputSigma[i];
708         }
709         final RealMatrix insigma = new Array2DRowRealMatrix(sigmaArray, false);
710         sigma = max(insigma); // overall standard deviation
711 
712         // initialize termination criteria
713         stopTolUpX = 1e3 * max(insigma);
714         stopTolX = 1e-11 * max(insigma);
715         stopTolFun = 1e-12;
716         stopTolHistFun = 1e-13;
717 
718         // initialize selection strategy parameters
719         mu = lambda / 2; // number of parents/points for recombination
720         logMu2 = FastMath.log(mu + 0.5);
721         weights = log(sequence(1, mu, 1)).scalarMultiply(-1).scalarAdd(logMu2);
722         double sumw = 0;
723         double sumwq = 0;
724         for (int i = 0; i < mu; i++) {
725             double w = weights.getEntry(i, 0);
726             sumw += w;
727             sumwq += w * w;
728         }
729         weights = weights.scalarMultiply(1 / sumw);
730         mueff = sumw * sumw / sumwq; // variance-effectiveness of sum w_i x_i
731 
732         // initialize dynamic strategy parameters and constants
733         cc = (4 + mueff / dimension) /
734                 (dimension + 4 + 2 * mueff / dimension);
735         cs = (mueff + 2) / (dimension + mueff + 3.);
736         damps = (1 + 2 * FastMath.max(0, FastMath.sqrt((mueff - 1) /
737                                                        (dimension + 1)) - 1)) *
738             FastMath.max(0.3,
739                          1 - dimension / (1e-6 + maxIterations)) + cs; // minor increment
740         ccov1 = 2 / ((dimension + 1.3) * (dimension + 1.3) + mueff);
741         ccovmu = FastMath.min(1 - ccov1, 2 * (mueff - 2 + 1 / mueff) /
742                               ((dimension + 2) * (dimension + 2) + mueff));
743         ccov1Sep = FastMath.min(1, ccov1 * (dimension + 1.5) / 3);
744         ccovmuSep = FastMath.min(1 - ccov1, ccovmu * (dimension + 1.5) / 3);
745         chiN = FastMath.sqrt(dimension) *
746             (1 - 1 / ((double) 4 * dimension) + 1 / ((double) 21 * dimension * dimension));
747         // intialize CMA internal values - updated each generation
748         xmean = MatrixUtils.createColumnRealMatrix(guess); // objective variables
749         diagD = insigma.scalarMultiply(1 / sigma);
750         diagC = square(diagD);
751         pc = zeros(dimension, 1); // evolution paths for C and sigma
752         ps = zeros(dimension, 1); // B defines the coordinate system
753         normps = ps.getFrobeniusNorm();
754 
755         B = eye(dimension, dimension);
756         D = ones(dimension, 1); // diagonal D defines the scaling
757         BD = times(B, repmat(diagD.transpose(), dimension, 1));
758         C = B.multiply(diag(square(D)).multiply(B.transpose())); // covariance
759         historySize = 10 + (int) (3 * 10 * dimension / (double) lambda);
760         fitnessHistory = new double[historySize]; // history of fitness values
761         for (int i = 0; i < historySize; i++) {
762             fitnessHistory[i] = Double.MAX_VALUE;
763         }
764     }
765 
766     /**
767      * Update of the evolution paths ps and pc.
768      *
769      * @param zmean Weighted row matrix of the gaussian random numbers generating
770      * the current offspring.
771      * @param xold xmean matrix of the previous generation.
772      * @return hsig flag indicating a small correction.
773      */
774     private boolean updateEvolutionPaths(RealMatrix zmean, RealMatrix xold) {
775         ps = ps.scalarMultiply(1 - cs).add(
776                 B.multiply(zmean).scalarMultiply(FastMath.sqrt(cs * (2 - cs) * mueff)));
777         normps = ps.getFrobeniusNorm();
778         final boolean hsig = normps /
779             FastMath.sqrt(1 - FastMath.pow(1 - cs, 2 * iterations)) /
780             chiN < 1.4 + 2 / ((double) dimension + 1);
781         pc = pc.scalarMultiply(1 - cc);
782         if (hsig) {
783             pc = pc.add(xmean.subtract(xold).scalarMultiply(FastMath.sqrt(cc * (2 - cc) * mueff) / sigma));
784         }
785         return hsig;
786     }
787 
788     /**
789      * Update of the covariance matrix C for diagonalOnly > 0
790      *
791      * @param hsig Flag indicating a small correction.
792      * @param bestArz Fitness-sorted matrix of the gaussian random values of the
793      * current offspring.
794      */
795     private void updateCovarianceDiagonalOnly(boolean hsig,
796                                               final RealMatrix bestArz) {
797         // minor correction if hsig==false
798         double oldFac = hsig ? 0 : ccov1Sep * cc * (2 - cc);
799         oldFac += 1 - ccov1Sep - ccovmuSep;
800         diagC = diagC.scalarMultiply(oldFac) // regard old matrix
801             .add(square(pc).scalarMultiply(ccov1Sep)) // plus rank one update
802             .add((times(diagC, square(bestArz).multiply(weights))) // plus rank mu update
803                  .scalarMultiply(ccovmuSep));
804         diagD = sqrt(diagC); // replaces eig(C)
805         if (diagonalOnly > 1 &&
806             iterations > diagonalOnly) {
807             // full covariance matrix from now on
808             diagonalOnly = 0;
809             B = eye(dimension, dimension);
810             BD = diag(diagD);
811             C = diag(diagC);
812         }
813     }
814 
815     /**
816      * Update of the covariance matrix C.
817      *
818      * @param hsig Flag indicating a small correction.
819      * @param bestArx Fitness-sorted matrix of the argument vectors producing the
820      * current offspring.
821      * @param arz Unsorted matrix containing the gaussian random values of the
822      * current offspring.
823      * @param arindex Indices indicating the fitness-order of the current offspring.
824      * @param xold xmean matrix of the previous generation.
825      */
826     private void updateCovariance(boolean hsig, final RealMatrix bestArx,
827                                   final RealMatrix arz, final int[] arindex,
828                                   final RealMatrix xold) {
829         double negccov = 0;
830         if (ccov1 + ccovmu > 0) {
831             final RealMatrix arpos = bestArx.subtract(repmat(xold, 1, mu))
832                 .scalarMultiply(1 / sigma); // mu difference vectors
833             final RealMatrix roneu = pc.multiply(pc.transpose())
834                 .scalarMultiply(ccov1); // rank one update
835             // minor correction if hsig==false
836             double oldFac = hsig ? 0 : ccov1 * cc * (2 - cc);
837             oldFac += 1 - ccov1 - ccovmu;
838             if (isActiveCMA) {
839                 // Adapt covariance matrix C active CMA
840                 negccov = (1 - ccovmu) * 0.25 * mueff / (FastMath.pow(dimension + 2, 1.5) + 2 * mueff);
841                 // keep at least 0.66 in all directions, small popsize are most
842                 // critical
843                 final double negminresidualvariance = 0.66;
844                 // where to make up for the variance loss
845                 final double negalphaold = 0.5;
846                 // prepare vectors, compute negative updating matrix Cneg
847                 final int[] arReverseIndex = reverse(arindex);
848                 RealMatrix arzneg = selectColumns(arz, MathArrays.copyOf(arReverseIndex, mu));
849                 RealMatrix arnorms = sqrt(sumRows(square(arzneg)));
850                 final int[] idxnorms = sortedIndices(arnorms.getRow(0));
851                 final RealMatrix arnormsSorted = selectColumns(arnorms, idxnorms);
852                 final int[] idxReverse = reverse(idxnorms);
853                 final RealMatrix arnormsReverse = selectColumns(arnorms, idxReverse);
854                 arnorms = divide(arnormsReverse, arnormsSorted);
855                 final int[] idxInv = inverse(idxnorms);
856                 final RealMatrix arnormsInv = selectColumns(arnorms, idxInv);
857                 // check and set learning rate negccov
858                 final double negcovMax = (1 - negminresidualvariance) /
859                     square(arnormsInv).multiply(weights).getEntry(0, 0);
860                 if (negccov > negcovMax) {
861                     negccov = negcovMax;
862                 }
863                 arzneg = times(arzneg, repmat(arnormsInv, dimension, 1));
864                 final RealMatrix artmp = BD.multiply(arzneg);
865                 final RealMatrix Cneg = artmp.multiply(diag(weights)).multiply(artmp.transpose());
866                 oldFac += negalphaold * negccov;
867                 C = C.scalarMultiply(oldFac)
868                     .add(roneu) // regard old matrix
869                     .add(arpos.scalarMultiply( // plus rank one update
870                                               ccovmu + (1 - negalphaold) * negccov) // plus rank mu update
871                          .multiply(times(repmat(weights, 1, dimension),
872                                          arpos.transpose())))
873                     .subtract(Cneg.scalarMultiply(negccov));
874             } else {
875                 // Adapt covariance matrix C - nonactive
876                 C = C.scalarMultiply(oldFac) // regard old matrix
877                     .add(roneu) // plus rank one update
878                     .add(arpos.scalarMultiply(ccovmu) // plus rank mu update
879                          .multiply(times(repmat(weights, 1, dimension),
880                                          arpos.transpose())));
881             }
882         }
883         updateBD(negccov);
884     }
885 
886     /**
887      * Update B and D from C.
888      *
889      * @param negccov Negative covariance factor.
890      */
891     private void updateBD(double negccov) {
892         if (ccov1 + ccovmu + negccov > 0 &&
893             (iterations % 1. / (ccov1 + ccovmu + negccov) / dimension / 10.) < 1) {
894             // to achieve O(N^2)
895             C = triu(C, 0).add(triu(C, 1).transpose());
896             // enforce symmetry to prevent complex numbers
897             final EigenDecomposition eig = new EigenDecomposition(C);
898             B = eig.getV(); // eigen decomposition, B==normalized eigenvectors
899             D = eig.getD();
900             diagD = diag(D);
901             if (min(diagD) <= 0) {
902                 for (int i = 0; i < dimension; i++) {
903                     if (diagD.getEntry(i, 0) < 0) {
904                         diagD.setEntry(i, 0, 0);
905                     }
906                 }
907                 final double tfac = max(diagD) / 1e14;
908                 C = C.add(eye(dimension, dimension).scalarMultiply(tfac));
909                 diagD = diagD.add(ones(dimension, 1).scalarMultiply(tfac));
910             }
911             if (max(diagD) > 1e14 * min(diagD)) {
912                 final double tfac = max(diagD) / 1e14 - min(diagD);
913                 C = C.add(eye(dimension, dimension).scalarMultiply(tfac));
914                 diagD = diagD.add(ones(dimension, 1).scalarMultiply(tfac));
915             }
916             diagC = diag(C);
917             diagD = sqrt(diagD); // D contains standard deviations now
918             BD = times(B, repmat(diagD.transpose(), dimension, 1)); // O(n^2)
919         }
920     }
921 
922     /**
923      * Pushes the current best fitness value in a history queue.
924      *
925      * @param vals History queue.
926      * @param val Current best fitness value.
927      */
928     private static void push(double[] vals, double val) {
929         for (int i = vals.length-1; i > 0; i--) {
930             vals[i] = vals[i-1];
931         }
932         vals[0] = val;
933     }
934 
935     /**
936      * Sorts fitness values.
937      *
938      * @param doubles Array of values to be sorted.
939      * @return a sorted array of indices pointing into doubles.
940      */
941     private int[] sortedIndices(final double[] doubles) {
942         final DoubleIndex[] dis = new DoubleIndex[doubles.length];
943         for (int i = 0; i < doubles.length; i++) {
944             dis[i] = new DoubleIndex(doubles[i], i);
945         }
946         Arrays.sort(dis);
947         final int[] indices = new int[doubles.length];
948         for (int i = 0; i < doubles.length; i++) {
949             indices[i] = dis[i].index;
950         }
951         return indices;
952     }
953 
954     /**
955      * Used to sort fitness values. Sorting is always in lower value first
956      * order.
957      */
958     private static class DoubleIndex implements Comparable<DoubleIndex> {
959         /** Value to compare. */
960         private final double value;
961         /** Index into sorted array. */
962         private final int index;
963 
964         /**
965          * @param value Value to compare.
966          * @param index Index into sorted array.
967          */
968         DoubleIndex(double value, int index) {
969             this.value = value;
970             this.index = index;
971         }
972 
973         /** {@inheritDoc} */
974         public int compareTo(DoubleIndex o) {
975             return Double.compare(value, o.value);
976         }
977 
978         /** {@inheritDoc} */
979         @Override
980         public boolean equals(Object other) {
981 
982             if (this == other) {
983                 return true;
984             }
985 
986             if (other instanceof DoubleIndex) {
987                 return Double.compare(value, ((DoubleIndex) other).value) == 0;
988             }
989 
990             return false;
991         }
992 
993         /** {@inheritDoc} */
994         @Override
995         public int hashCode() {
996             long bits = Double.doubleToLongBits(value);
997             return (int) ((1438542 ^ (bits >>> 32) ^ bits) & 0xffffffff);
998         }
999     }
1000 
1001     /**
1002      * Normalizes fitness values to the range [0,1]. Adds a penalty to the
1003      * fitness value if out of range. The penalty is adjusted by calling
1004      * setValueRange().
1005      */
1006     private class FitnessFunction {
1007         /** Determines the penalty for boundary violations */
1008         private double valueRange;
1009         /**
1010          * Flag indicating whether the objective variables are forced into their
1011          * bounds if defined
1012          */
1013         private final boolean isRepairMode;
1014 
1015         /** Simple constructor.
1016          */
1017         public FitnessFunction() {
1018             valueRange = 1;
1019             isRepairMode = true;
1020         }
1021 
1022         /**
1023          * @param point Normalized objective variables.
1024          * @return the objective value + penalty for violated bounds.
1025          */
1026         public double value(final double[] point) {
1027             double value;
1028             if (isRepairMode) {
1029                 double[] repaired = repair(point);
1030                 value = CMAESOptimizer.this.computeObjectiveValue(repaired) +
1031                     penalty(point, repaired);
1032             } else {
1033                 value = CMAESOptimizer.this.computeObjectiveValue(point);
1034             }
1035             return isMinimize ? value : -value;
1036         }
1037 
1038         /**
1039          * @param x Normalized objective variables.
1040          * @return {@code true} if in bounds.
1041          */
1042         public boolean isFeasible(final double[] x) {
1043             final double[] lB = CMAESOptimizer.this.getLowerBound();
1044             final double[] uB = CMAESOptimizer.this.getUpperBound();
1045 
1046             for (int i = 0; i < x.length; i++) {
1047                 if (x[i] < lB[i]) {
1048                     return false;
1049                 }
1050                 if (x[i] > uB[i]) {
1051                     return false;
1052                 }
1053             }
1054             return true;
1055         }
1056 
1057         /**
1058          * @param valueRange Adjusts the penalty computation.
1059          */
1060         public void setValueRange(double valueRange) {
1061             this.valueRange = valueRange;
1062         }
1063 
1064         /**
1065          * @param x Normalized objective variables.
1066          * @return the repaired (i.e. all in bounds) objective variables.
1067          */
1068         private double[] repair(final double[] x) {
1069             final double[] lB = CMAESOptimizer.this.getLowerBound();
1070             final double[] uB = CMAESOptimizer.this.getUpperBound();
1071 
1072             final double[] repaired = new double[x.length];
1073             for (int i = 0; i < x.length; i++) {
1074                 if (x[i] < lB[i]) {
1075                     repaired[i] = lB[i];
1076                 } else if (x[i] > uB[i]) {
1077                     repaired[i] = uB[i];
1078                 } else {
1079                     repaired[i] = x[i];
1080                 }
1081             }
1082             return repaired;
1083         }
1084 
1085         /**
1086          * @param x Normalized objective variables.
1087          * @param repaired Repaired objective variables.
1088          * @return Penalty value according to the violation of the bounds.
1089          */
1090         private double penalty(final double[] x, final double[] repaired) {
1091             double penalty = 0;
1092             for (int i = 0; i < x.length; i++) {
1093                 double diff = FastMath.abs(x[i] - repaired[i]);
1094                 penalty += diff * valueRange;
1095             }
1096             return isMinimize ? penalty : -penalty;
1097         }
1098     }
1099 
1100     // -----Matrix utility functions similar to the Matlab build in functions------
1101 
1102     /**
1103      * @param m Input matrix
1104      * @return Matrix representing the element-wise logarithm of m.
1105      */
1106     private static RealMatrix log(final RealMatrix m) {
1107         final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
1108         for (int r = 0; r < m.getRowDimension(); r++) {
1109             for (int c = 0; c < m.getColumnDimension(); c++) {
1110                 d[r][c] = FastMath.log(m.getEntry(r, c));
1111             }
1112         }
1113         return new Array2DRowRealMatrix(d, false);
1114     }
1115 
1116     /**
1117      * @param m Input matrix.
1118      * @return Matrix representing the element-wise square root of m.
1119      */
1120     private static RealMatrix sqrt(final RealMatrix m) {
1121         final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
1122         for (int r = 0; r < m.getRowDimension(); r++) {
1123             for (int c = 0; c < m.getColumnDimension(); c++) {
1124                 d[r][c] = FastMath.sqrt(m.getEntry(r, c));
1125             }
1126         }
1127         return new Array2DRowRealMatrix(d, false);
1128     }
1129 
1130     /**
1131      * @param m Input matrix.
1132      * @return Matrix representing the element-wise square of m.
1133      */
1134     private static RealMatrix square(final RealMatrix m) {
1135         final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
1136         for (int r = 0; r < m.getRowDimension(); r++) {
1137             for (int c = 0; c < m.getColumnDimension(); c++) {
1138                 double e = m.getEntry(r, c);
1139                 d[r][c] = e * e;
1140             }
1141         }
1142         return new Array2DRowRealMatrix(d, false);
1143     }
1144 
1145     /**
1146      * @param m Input matrix 1.
1147      * @param n Input matrix 2.
1148      * @return the matrix where the elements of m and n are element-wise multiplied.
1149      */
1150     private static RealMatrix times(final RealMatrix m, final RealMatrix n) {
1151         final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
1152         for (int r = 0; r < m.getRowDimension(); r++) {
1153             for (int c = 0; c < m.getColumnDimension(); c++) {
1154                 d[r][c] = m.getEntry(r, c) * n.getEntry(r, c);
1155             }
1156         }
1157         return new Array2DRowRealMatrix(d, false);
1158     }
1159 
1160     /**
1161      * @param m Input matrix 1.
1162      * @param n Input matrix 2.
1163      * @return Matrix where the elements of m and n are element-wise divided.
1164      */
1165     private static RealMatrix divide(final RealMatrix m, final RealMatrix n) {
1166         final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
1167         for (int r = 0; r < m.getRowDimension(); r++) {
1168             for (int c = 0; c < m.getColumnDimension(); c++) {
1169                 d[r][c] = m.getEntry(r, c) / n.getEntry(r, c);
1170             }
1171         }
1172         return new Array2DRowRealMatrix(d, false);
1173     }
1174 
1175     /**
1176      * @param m Input matrix.
1177      * @param cols Columns to select.
1178      * @return Matrix representing the selected columns.
1179      */
1180     private static RealMatrix selectColumns(final RealMatrix m, final int[] cols) {
1181         final double[][] d = new double[m.getRowDimension()][cols.length];
1182         for (int r = 0; r < m.getRowDimension(); r++) {
1183             for (int c = 0; c < cols.length; c++) {
1184                 d[r][c] = m.getEntry(r, cols[c]);
1185             }
1186         }
1187         return new Array2DRowRealMatrix(d, false);
1188     }
1189 
1190     /**
1191      * @param m Input matrix.
1192      * @param k Diagonal position.
1193      * @return Upper triangular part of matrix.
1194      */
1195     private static RealMatrix triu(final RealMatrix m, int k) {
1196         final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
1197         for (int r = 0; r < m.getRowDimension(); r++) {
1198             for (int c = 0; c < m.getColumnDimension(); c++) {
1199                 d[r][c] = r <= c - k ? m.getEntry(r, c) : 0;
1200             }
1201         }
1202         return new Array2DRowRealMatrix(d, false);
1203     }
1204 
1205     /**
1206      * @param m Input matrix.
1207      * @return Row matrix representing the sums of the rows.
1208      */
1209     private static RealMatrix sumRows(final RealMatrix m) {
1210         final double[][] d = new double[1][m.getColumnDimension()];
1211         for (int c = 0; c < m.getColumnDimension(); c++) {
1212             double sum = 0;
1213             for (int r = 0; r < m.getRowDimension(); r++) {
1214                 sum += m.getEntry(r, c);
1215             }
1216             d[0][c] = sum;
1217         }
1218         return new Array2DRowRealMatrix(d, false);
1219     }
1220 
1221     /**
1222      * @param m Input matrix.
1223      * @return the diagonal n-by-n matrix if m is a column matrix or the column
1224      * matrix representing the diagonal if m is a n-by-n matrix.
1225      */
1226     private static RealMatrix diag(final RealMatrix m) {
1227         if (m.getColumnDimension() == 1) {
1228             final double[][] d = new double[m.getRowDimension()][m.getRowDimension()];
1229             for (int i = 0; i < m.getRowDimension(); i++) {
1230                 d[i][i] = m.getEntry(i, 0);
1231             }
1232             return new Array2DRowRealMatrix(d, false);
1233         } else {
1234             final double[][] d = new double[m.getRowDimension()][1];
1235             for (int i = 0; i < m.getColumnDimension(); i++) {
1236                 d[i][0] = m.getEntry(i, i);
1237             }
1238             return new Array2DRowRealMatrix(d, false);
1239         }
1240     }
1241 
1242     /**
1243      * Copies a column from m1 to m2.
1244      *
1245      * @param m1 Source matrix.
1246      * @param col1 Source column.
1247      * @param m2 Target matrix.
1248      * @param col2 Target column.
1249      */
1250     private static void copyColumn(final RealMatrix m1, int col1,
1251                                    RealMatrix m2, int col2) {
1252         for (int i = 0; i < m1.getRowDimension(); i++) {
1253             m2.setEntry(i, col2, m1.getEntry(i, col1));
1254         }
1255     }
1256 
1257     /**
1258      * @param n Number of rows.
1259      * @param m Number of columns.
1260      * @return n-by-m matrix filled with 1.
1261      */
1262     private static RealMatrix ones(int n, int m) {
1263         final double[][] d = new double[n][m];
1264         for (int r = 0; r < n; r++) {
1265             Arrays.fill(d[r], 1);
1266         }
1267         return new Array2DRowRealMatrix(d, false);
1268     }
1269 
1270     /**
1271      * @param n Number of rows.
1272      * @param m Number of columns.
1273      * @return n-by-m matrix of 0 values out of diagonal, and 1 values on
1274      * the diagonal.
1275      */
1276     private static RealMatrix eye(int n, int m) {
1277         final double[][] d = new double[n][m];
1278         for (int r = 0; r < n; r++) {
1279             if (r < m) {
1280                 d[r][r] = 1;
1281             }
1282         }
1283         return new Array2DRowRealMatrix(d, false);
1284     }
1285 
1286     /**
1287      * @param n Number of rows.
1288      * @param m Number of columns.
1289      * @return n-by-m matrix of zero values.
1290      */
1291     private static RealMatrix zeros(int n, int m) {
1292         return new Array2DRowRealMatrix(n, m);
1293     }
1294 
1295     /**
1296      * @param mat Input matrix.
1297      * @param n Number of row replicates.
1298      * @param m Number of column replicates.
1299      * @return a matrix which replicates the input matrix in both directions.
1300      */
1301     private static RealMatrix repmat(final RealMatrix mat, int n, int m) {
1302         final int rd = mat.getRowDimension();
1303         final int cd = mat.getColumnDimension();
1304         final double[][] d = new double[n * rd][m * cd];
1305         for (int r = 0; r < n * rd; r++) {
1306             for (int c = 0; c < m * cd; c++) {
1307                 d[r][c] = mat.getEntry(r % rd, c % cd);
1308             }
1309         }
1310         return new Array2DRowRealMatrix(d, false);
1311     }
1312 
1313     /**
1314      * @param start Start value.
1315      * @param end End value.
1316      * @param step Step size.
1317      * @return a sequence as column matrix.
1318      */
1319     private static RealMatrix sequence(double start, double end, double step) {
1320         final int size = (int) ((end - start) / step + 1);
1321         final double[][] d = new double[size][1];
1322         double value = start;
1323         for (int r = 0; r < size; r++) {
1324             d[r][0] = value;
1325             value += step;
1326         }
1327         return new Array2DRowRealMatrix(d, false);
1328     }
1329 
1330     /**
1331      * @param m Input matrix.
1332      * @return the maximum of the matrix element values.
1333      */
1334     private static double max(final RealMatrix m) {
1335         double max = -Double.MAX_VALUE;
1336         for (int r = 0; r < m.getRowDimension(); r++) {
1337             for (int c = 0; c < m.getColumnDimension(); c++) {
1338                 double e = m.getEntry(r, c);
1339                 if (max < e) {
1340                     max = e;
1341                 }
1342             }
1343         }
1344         return max;
1345     }
1346 
1347     /**
1348      * @param m Input matrix.
1349      * @return the minimum of the matrix element values.
1350      */
1351     private static double min(final RealMatrix m) {
1352         double min = Double.MAX_VALUE;
1353         for (int r = 0; r < m.getRowDimension(); r++) {
1354             for (int c = 0; c < m.getColumnDimension(); c++) {
1355                 double e = m.getEntry(r, c);
1356                 if (min > e) {
1357                     min = e;
1358                 }
1359             }
1360         }
1361         return min;
1362     }
1363 
1364     /**
1365      * @param m Input array.
1366      * @return the maximum of the array values.
1367      */
1368     private static double max(final double[] m) {
1369         double max = -Double.MAX_VALUE;
1370         for (int r = 0; r < m.length; r++) {
1371             if (max < m[r]) {
1372                 max = m[r];
1373             }
1374         }
1375         return max;
1376     }
1377 
1378     /**
1379      * @param m Input array.
1380      * @return the minimum of the array values.
1381      */
1382     private static double min(final double[] m) {
1383         double min = Double.MAX_VALUE;
1384         for (int r = 0; r < m.length; r++) {
1385             if (min > m[r]) {
1386                 min = m[r];
1387             }
1388         }
1389         return min;
1390     }
1391 
1392     /**
1393      * @param indices Input index array.
1394      * @return the inverse of the mapping defined by indices.
1395      */
1396     private static int[] inverse(final int[] indices) {
1397         final int[] inverse = new int[indices.length];
1398         for (int i = 0; i < indices.length; i++) {
1399             inverse[indices[i]] = i;
1400         }
1401         return inverse;
1402     }
1403 
1404     /**
1405      * @param indices Input index array.
1406      * @return the indices in inverse order (last is first).
1407      */
1408     private static int[] reverse(final int[] indices) {
1409         final int[] reverse = new int[indices.length];
1410         for (int i = 0; i < indices.length; i++) {
1411             reverse[i] = indices[indices.length - i - 1];
1412         }
1413         return reverse;
1414     }
1415 
1416     /**
1417      * @param size Length of random array.
1418      * @return an array of Gaussian random numbers.
1419      */
1420     private double[] randn(int size) {
1421         final double[] randn = new double[size];
1422         for (int i = 0; i < size; i++) {
1423             randn[i] = random.nextGaussian();
1424         }
1425         return randn;
1426     }
1427 
1428     /**
1429      * @param size Number of rows.
1430      * @param popSize Population size.
1431      * @return a 2-dimensional matrix of Gaussian random numbers.
1432      */
1433     private RealMatrix randn1(int size, int popSize) {
1434         final double[][] d = new double[size][popSize];
1435         for (int r = 0; r < size; r++) {
1436             for (int c = 0; c < popSize; c++) {
1437                 d[r][c] = random.nextGaussian();
1438             }
1439         }
1440         return new Array2DRowRealMatrix(d, false);
1441     }
1442 }