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