CMAESOptimizer.java

  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. package org.apache.commons.math3.optimization.direct;

  18. import java.util.ArrayList;
  19. import java.util.Arrays;
  20. import java.util.List;

  21. import org.apache.commons.math3.analysis.MultivariateFunction;
  22. import org.apache.commons.math3.exception.DimensionMismatchException;
  23. import org.apache.commons.math3.exception.NotPositiveException;
  24. import org.apache.commons.math3.exception.NotStrictlyPositiveException;
  25. import org.apache.commons.math3.exception.OutOfRangeException;
  26. import org.apache.commons.math3.exception.TooManyEvaluationsException;
  27. import org.apache.commons.math3.linear.Array2DRowRealMatrix;
  28. import org.apache.commons.math3.linear.EigenDecomposition;
  29. import org.apache.commons.math3.linear.MatrixUtils;
  30. import org.apache.commons.math3.linear.RealMatrix;
  31. import org.apache.commons.math3.optimization.ConvergenceChecker;
  32. import org.apache.commons.math3.optimization.OptimizationData;
  33. import org.apache.commons.math3.optimization.GoalType;
  34. import org.apache.commons.math3.optimization.MultivariateOptimizer;
  35. import org.apache.commons.math3.optimization.PointValuePair;
  36. import org.apache.commons.math3.optimization.SimpleValueChecker;
  37. import org.apache.commons.math3.random.MersenneTwister;
  38. import org.apache.commons.math3.random.RandomGenerator;
  39. import org.apache.commons.math3.util.FastMath;
  40. import org.apache.commons.math3.util.MathArrays;

  41. /**
  42.  * <p>An implementation of the active Covariance Matrix Adaptation Evolution Strategy (CMA-ES)
  43.  * for non-linear, non-convex, non-smooth, global function minimization.
  44.  * The CMA-Evolution Strategy (CMA-ES) is a reliable stochastic optimization method
  45.  * which should be applied if derivative-based methods, e.g. quasi-Newton BFGS or
  46.  * conjugate gradient, fail due to a rugged search landscape (e.g. noise, local
  47.  * optima, outlier, etc.) of the objective function. Like a
  48.  * quasi-Newton method, the CMA-ES learns and applies a variable metric
  49.  * on the underlying search space. Unlike a quasi-Newton method, the
  50.  * CMA-ES neither estimates nor uses gradients, making it considerably more
  51.  * reliable in terms of finding a good, or even close to optimal, solution.</p>
  52.  *
  53.  * <p>In general, on smooth objective functions the CMA-ES is roughly ten times
  54.  * slower than BFGS (counting objective function evaluations, no gradients provided).
  55.  * For up to <math>N=10</math> variables also the derivative-free simplex
  56.  * direct search method (Nelder and Mead) can be faster, but it is
  57.  * far less reliable than CMA-ES.</p>
  58.  *
  59.  * <p>The CMA-ES is particularly well suited for non-separable
  60.  * and/or badly conditioned problems. To observe the advantage of CMA compared
  61.  * to a conventional evolution strategy, it will usually take about
  62.  * <math>30 N</math> function evaluations. On difficult problems the complete
  63.  * optimization (a single run) is expected to take <em>roughly</em> between
  64.  * <math>30 N</math> and <math>300 N<sup>2</sup></math>
  65.  * function evaluations.</p>
  66.  *
  67.  * <p>This implementation is translated and adapted from the Matlab version
  68.  * of the CMA-ES algorithm as implemented in module {@code cmaes.m} version 3.51.</p>
  69.  *
  70.  * For more information, please refer to the following links:
  71.  * <ul>
  72.  *  <li><a href="http://www.lri.fr/~hansen/cmaes.m">Matlab code</a></li>
  73.  *  <li><a href="http://www.lri.fr/~hansen/cmaesintro.html">Introduction to CMA-ES</a></li>
  74.  *  <li><a href="http://en.wikipedia.org/wiki/CMA-ES">Wikipedia</a></li>
  75.  * </ul>
  76.  *
  77.  * @deprecated As of 3.1 (to be removed in 4.0).
  78.  * @since 3.0
  79.  */
  80. @Deprecated
  81. public class CMAESOptimizer
  82.     extends BaseAbstractMultivariateSimpleBoundsOptimizer<MultivariateFunction>
  83.     implements MultivariateOptimizer {
  84.     /** Default value for {@link #checkFeasableCount}: {@value}. */
  85.     public static final int DEFAULT_CHECKFEASABLECOUNT = 0;
  86.     /** Default value for {@link #stopFitness}: {@value}. */
  87.     public static final double DEFAULT_STOPFITNESS = 0;
  88.     /** Default value for {@link #isActiveCMA}: {@value}. */
  89.     public static final boolean DEFAULT_ISACTIVECMA = true;
  90.     /** Default value for {@link #maxIterations}: {@value}. */
  91.     public static final int DEFAULT_MAXITERATIONS = 30000;
  92.     /** Default value for {@link #diagonalOnly}: {@value}. */
  93.     public static final int DEFAULT_DIAGONALONLY = 0;
  94.     /** Default value for {@link #random}. */
  95.     public static final RandomGenerator DEFAULT_RANDOMGENERATOR = new MersenneTwister();

  96.     // global search parameters
  97.     /**
  98.      * Population size, offspring number. The primary strategy parameter to play
  99.      * with, which can be increased from its default value. Increasing the
  100.      * population size improves global search properties in exchange to speed.
  101.      * Speed decreases, as a rule, at most linearly with increasing population
  102.      * size. It is advisable to begin with the default small population size.
  103.      */
  104.     private int lambda; // population size
  105.     /**
  106.      * Covariance update mechanism, default is active CMA. isActiveCMA = true
  107.      * turns on "active CMA" with a negative update of the covariance matrix and
  108.      * checks for positive definiteness. OPTS.CMA.active = 2 does not check for
  109.      * pos. def. and is numerically faster. Active CMA usually speeds up the
  110.      * adaptation.
  111.      */
  112.     private boolean isActiveCMA;
  113.     /**
  114.      * Determines how often a new random offspring is generated in case it is
  115.      * not feasible / beyond the defined limits, default is 0.
  116.      */
  117.     private int checkFeasableCount;
  118.     /**
  119.      * @see Sigma
  120.      */
  121.     private double[] inputSigma;
  122.     /** Number of objective variables/problem dimension */
  123.     private int dimension;
  124.     /**
  125.      * Defines the number of initial iterations, where the covariance matrix
  126.      * remains diagonal and the algorithm has internally linear time complexity.
  127.      * diagonalOnly = 1 means keeping the covariance matrix always diagonal and
  128.      * this setting also exhibits linear space complexity. This can be
  129.      * particularly useful for dimension > 100.
  130.      * @see <a href="http://hal.archives-ouvertes.fr/inria-00287367/en">A Simple Modification in CMA-ES</a>
  131.      */
  132.     private int diagonalOnly = 0;
  133.     /** Number of objective variables/problem dimension */
  134.     private boolean isMinimize = true;
  135.     /** Indicates whether statistic data is collected. */
  136.     private boolean generateStatistics = false;

  137.     // termination criteria
  138.     /** Maximal number of iterations allowed. */
  139.     private int maxIterations;
  140.     /** Limit for fitness value. */
  141.     private double stopFitness;
  142.     /** Stop if x-changes larger stopTolUpX. */
  143.     private double stopTolUpX;
  144.     /** Stop if x-change smaller stopTolX. */
  145.     private double stopTolX;
  146.     /** Stop if fun-changes smaller stopTolFun. */
  147.     private double stopTolFun;
  148.     /** Stop if back fun-changes smaller stopTolHistFun. */
  149.     private double stopTolHistFun;

  150.     // selection strategy parameters
  151.     /** Number of parents/points for recombination. */
  152.     private int mu; //
  153.     /** log(mu + 0.5), stored for efficiency. */
  154.     private double logMu2;
  155.     /** Array for weighted recombination. */
  156.     private RealMatrix weights;
  157.     /** Variance-effectiveness of sum w_i x_i. */
  158.     private double mueff; //

  159.     // dynamic strategy parameters and constants
  160.     /** Overall standard deviation - search volume. */
  161.     private double sigma;
  162.     /** Cumulation constant. */
  163.     private double cc;
  164.     /** Cumulation constant for step-size. */
  165.     private double cs;
  166.     /** Damping for step-size. */
  167.     private double damps;
  168.     /** Learning rate for rank-one update. */
  169.     private double ccov1;
  170.     /** Learning rate for rank-mu update' */
  171.     private double ccovmu;
  172.     /** Expectation of ||N(0,I)|| == norm(randn(N,1)). */
  173.     private double chiN;
  174.     /** Learning rate for rank-one update - diagonalOnly */
  175.     private double ccov1Sep;
  176.     /** Learning rate for rank-mu update - diagonalOnly */
  177.     private double ccovmuSep;

  178.     // CMA internal values - updated each generation
  179.     /** Objective variables. */
  180.     private RealMatrix xmean;
  181.     /** Evolution path. */
  182.     private RealMatrix pc;
  183.     /** Evolution path for sigma. */
  184.     private RealMatrix ps;
  185.     /** Norm of ps, stored for efficiency. */
  186.     private double normps;
  187.     /** Coordinate system. */
  188.     private RealMatrix B;
  189.     /** Scaling. */
  190.     private RealMatrix D;
  191.     /** B*D, stored for efficiency. */
  192.     private RealMatrix BD;
  193.     /** Diagonal of sqrt(D), stored for efficiency. */
  194.     private RealMatrix diagD;
  195.     /** Covariance matrix. */
  196.     private RealMatrix C;
  197.     /** Diagonal of C, used for diagonalOnly. */
  198.     private RealMatrix diagC;
  199.     /** Number of iterations already performed. */
  200.     private int iterations;

  201.     /** History queue of best values. */
  202.     private double[] fitnessHistory;
  203.     /** Size of history queue of best values. */
  204.     private int historySize;

  205.     /** Random generator. */
  206.     private RandomGenerator random;

  207.     /** History of sigma values. */
  208.     private List<Double> statisticsSigmaHistory = new ArrayList<Double>();
  209.     /** History of mean matrix. */
  210.     private List<RealMatrix> statisticsMeanHistory = new ArrayList<RealMatrix>();
  211.     /** History of fitness values. */
  212.     private List<Double> statisticsFitnessHistory = new ArrayList<Double>();
  213.     /** History of D matrix. */
  214.     private List<RealMatrix> statisticsDHistory = new ArrayList<RealMatrix>();

  215.     /**
  216.      * Default constructor, uses default parameters
  217.      *
  218.      * @deprecated As of version 3.1: Parameter {@code lambda} must be
  219.      * passed with the call to {@link #optimize(int,MultivariateFunction,GoalType,OptimizationData[])
  220.      * optimize} (whereas in the current code it is set to an undocumented value).
  221.      */
  222.     @Deprecated
  223.     public CMAESOptimizer() {
  224.         this(0);
  225.     }

  226.     /**
  227.      * @param lambda Population size.
  228.      * @deprecated As of version 3.1: Parameter {@code lambda} must be
  229.      * passed with the call to {@link #optimize(int,MultivariateFunction,GoalType,OptimizationData[])
  230.      * optimize} (whereas in the current code it is set to an undocumented value)..
  231.      */
  232.     @Deprecated
  233.     public CMAESOptimizer(int lambda) {
  234.         this(lambda, null, DEFAULT_MAXITERATIONS, DEFAULT_STOPFITNESS,
  235.              DEFAULT_ISACTIVECMA, DEFAULT_DIAGONALONLY,
  236.              DEFAULT_CHECKFEASABLECOUNT, DEFAULT_RANDOMGENERATOR,
  237.              false, null);
  238.     }

  239.     /**
  240.      * @param lambda Population size.
  241.      * @param inputSigma Initial standard deviations to sample new points
  242.      * around the initial guess.
  243.      * @deprecated As of version 3.1: Parameters {@code lambda} and {@code inputSigma} must be
  244.      * passed with the call to {@link #optimize(int,MultivariateFunction,GoalType,OptimizationData[])
  245.      * optimize}.
  246.      */
  247.     @Deprecated
  248.     public CMAESOptimizer(int lambda, double[] inputSigma) {
  249.         this(lambda, inputSigma, DEFAULT_MAXITERATIONS, DEFAULT_STOPFITNESS,
  250.              DEFAULT_ISACTIVECMA, DEFAULT_DIAGONALONLY,
  251.              DEFAULT_CHECKFEASABLECOUNT, DEFAULT_RANDOMGENERATOR, false);
  252.     }

  253.     /**
  254.      * @param lambda Population size.
  255.      * @param inputSigma Initial standard deviations to sample new points
  256.      * around the initial guess.
  257.      * @param maxIterations Maximal number of iterations.
  258.      * @param stopFitness Whether to stop if objective function value is smaller than
  259.      * {@code stopFitness}.
  260.      * @param isActiveCMA Chooses the covariance matrix update method.
  261.      * @param diagonalOnly Number of initial iterations, where the covariance matrix
  262.      * remains diagonal.
  263.      * @param checkFeasableCount Determines how often new random objective variables are
  264.      * generated in case they are out of bounds.
  265.      * @param random Random generator.
  266.      * @param generateStatistics Whether statistic data is collected.
  267.      * @deprecated See {@link SimpleValueChecker#SimpleValueChecker()}
  268.      */
  269.     @Deprecated
  270.     public CMAESOptimizer(int lambda, double[] inputSigma,
  271.                           int maxIterations, double stopFitness,
  272.                           boolean isActiveCMA, int diagonalOnly, int checkFeasableCount,
  273.                           RandomGenerator random, boolean generateStatistics) {
  274.         this(lambda, inputSigma, maxIterations, stopFitness, isActiveCMA,
  275.              diagonalOnly, checkFeasableCount, random, generateStatistics,
  276.              new SimpleValueChecker());
  277.     }

  278.     /**
  279.      * @param lambda Population size.
  280.      * @param inputSigma Initial standard deviations to sample new points
  281.      * around the initial guess.
  282.      * @param maxIterations Maximal number of iterations.
  283.      * @param stopFitness Whether to stop if objective function value is smaller than
  284.      * {@code stopFitness}.
  285.      * @param isActiveCMA Chooses the covariance matrix update method.
  286.      * @param diagonalOnly Number of initial iterations, where the covariance matrix
  287.      * remains diagonal.
  288.      * @param checkFeasableCount Determines how often new random objective variables are
  289.      * generated in case they are out of bounds.
  290.      * @param random Random generator.
  291.      * @param generateStatistics Whether statistic data is collected.
  292.      * @param checker Convergence checker.
  293.      * @deprecated As of version 3.1: Parameters {@code lambda} and {@code inputSigma} must be
  294.      * passed with the call to {@link #optimize(int,MultivariateFunction,GoalType,OptimizationData[])
  295.      * optimize}.
  296.      */
  297.     @Deprecated
  298.     public CMAESOptimizer(int lambda, double[] inputSigma,
  299.                           int maxIterations, double stopFitness,
  300.                           boolean isActiveCMA, int diagonalOnly, int checkFeasableCount,
  301.                           RandomGenerator random, boolean generateStatistics,
  302.                           ConvergenceChecker<PointValuePair> checker) {
  303.         super(checker);
  304.         this.lambda = lambda;
  305.         this.inputSigma = inputSigma == null ? null : (double[]) inputSigma.clone();
  306.         this.maxIterations = maxIterations;
  307.         this.stopFitness = stopFitness;
  308.         this.isActiveCMA = isActiveCMA;
  309.         this.diagonalOnly = diagonalOnly;
  310.         this.checkFeasableCount = checkFeasableCount;
  311.         this.random = random;
  312.         this.generateStatistics = generateStatistics;
  313.     }

  314.     /**
  315.      * @param maxIterations Maximal number of iterations.
  316.      * @param stopFitness Whether to stop if objective function value is smaller than
  317.      * {@code stopFitness}.
  318.      * @param isActiveCMA Chooses the covariance matrix update method.
  319.      * @param diagonalOnly Number of initial iterations, where the covariance matrix
  320.      * remains diagonal.
  321.      * @param checkFeasableCount Determines how often new random objective variables are
  322.      * generated in case they are out of bounds.
  323.      * @param random Random generator.
  324.      * @param generateStatistics Whether statistic data is collected.
  325.      * @param checker Convergence checker.
  326.      *
  327.      * @since 3.1
  328.      */
  329.     public CMAESOptimizer(int maxIterations,
  330.                           double stopFitness,
  331.                           boolean isActiveCMA,
  332.                           int diagonalOnly,
  333.                           int checkFeasableCount,
  334.                           RandomGenerator random,
  335.                           boolean generateStatistics,
  336.                           ConvergenceChecker<PointValuePair> checker) {
  337.         super(checker);
  338.         this.maxIterations = maxIterations;
  339.         this.stopFitness = stopFitness;
  340.         this.isActiveCMA = isActiveCMA;
  341.         this.diagonalOnly = diagonalOnly;
  342.         this.checkFeasableCount = checkFeasableCount;
  343.         this.random = random;
  344.         this.generateStatistics = generateStatistics;
  345.     }

  346.     /**
  347.      * @return History of sigma values.
  348.      */
  349.     public List<Double> getStatisticsSigmaHistory() {
  350.         return statisticsSigmaHistory;
  351.     }

  352.     /**
  353.      * @return History of mean matrix.
  354.      */
  355.     public List<RealMatrix> getStatisticsMeanHistory() {
  356.         return statisticsMeanHistory;
  357.     }

  358.     /**
  359.      * @return History of fitness values.
  360.      */
  361.     public List<Double> getStatisticsFitnessHistory() {
  362.         return statisticsFitnessHistory;
  363.     }

  364.     /**
  365.      * @return History of D matrix.
  366.      */
  367.     public List<RealMatrix> getStatisticsDHistory() {
  368.         return statisticsDHistory;
  369.     }

  370.     /**
  371.      * Input sigma values.
  372.      * They define the initial coordinate-wise standard deviations for
  373.      * sampling new search points around the initial guess.
  374.      * It is suggested to set them to the estimated distance from the
  375.      * initial to the desired optimum.
  376.      * Small values induce the search to be more local (and very small
  377.      * values are more likely to find a local optimum close to the initial
  378.      * guess).
  379.      * Too small values might however lead to early termination.
  380.      * @since 3.1
  381.      */
  382.     public static class Sigma implements OptimizationData {
  383.         /** Sigma values. */
  384.         private final double[] sigma;

  385.         /**
  386.          * @param s Sigma values.
  387.          * @throws NotPositiveException if any of the array entries is smaller
  388.          * than zero.
  389.          */
  390.         public Sigma(double[] s)
  391.             throws NotPositiveException {
  392.             for (int i = 0; i < s.length; i++) {
  393.                 if (s[i] < 0) {
  394.                     throw new NotPositiveException(s[i]);
  395.                 }
  396.             }

  397.             sigma = s.clone();
  398.         }

  399.         /**
  400.          * @return the sigma values.
  401.          */
  402.         public double[] getSigma() {
  403.             return sigma.clone();
  404.         }
  405.     }

  406.     /**
  407.      * Population size.
  408.      * The number of offspring is the primary strategy parameter.
  409.      * In the absence of better clues, a good default could be an
  410.      * integer close to {@code 4 + 3 ln(n)}, where {@code n} is the
  411.      * number of optimized parameters.
  412.      * Increasing the population size improves global search properties
  413.      * at the expense of speed (which in general decreases at most
  414.      * linearly with increasing population size).
  415.      * @since 3.1
  416.      */
  417.     public static class PopulationSize implements OptimizationData {
  418.         /** Population size. */
  419.         private final int lambda;

  420.         /**
  421.          * @param size Population size.
  422.          * @throws NotStrictlyPositiveException if {@code size <= 0}.
  423.          */
  424.         public PopulationSize(int size)
  425.             throws NotStrictlyPositiveException {
  426.             if (size <= 0) {
  427.                 throw new NotStrictlyPositiveException(size);
  428.             }
  429.             lambda = size;
  430.         }

  431.         /**
  432.          * @return the population size.
  433.          */
  434.         public int getPopulationSize() {
  435.             return lambda;
  436.         }
  437.     }

  438.     /**
  439.      * Optimize an objective function.
  440.      *
  441.      * @param maxEval Allowed number of evaluations of the objective function.
  442.      * @param f Objective function.
  443.      * @param goalType Optimization type.
  444.      * @param optData Optimization data. The following data will be looked for:
  445.      * <ul>
  446.      *  <li>{@link org.apache.commons.math3.optimization.InitialGuess InitialGuess}</li>
  447.      *  <li>{@link Sigma}</li>
  448.      *  <li>{@link PopulationSize}</li>
  449.      * </ul>
  450.      * @return the point/value pair giving the optimal value for objective
  451.      * function.
  452.      */
  453.     @Override
  454.     protected PointValuePair optimizeInternal(int maxEval, MultivariateFunction f,
  455.                                               GoalType goalType,
  456.                                               OptimizationData... optData) {
  457.         // Scan "optData" for the input specific to this optimizer.
  458.         parseOptimizationData(optData);

  459.         // The parent's method will retrieve the common parameters from
  460.         // "optData" and call "doOptimize".
  461.         return super.optimizeInternal(maxEval, f, goalType, optData);
  462.     }

  463.     /** {@inheritDoc} */
  464.     @Override
  465.     protected PointValuePair doOptimize() {
  466.         checkParameters();
  467.          // -------------------- Initialization --------------------------------
  468.         isMinimize = getGoalType().equals(GoalType.MINIMIZE);
  469.         final FitnessFunction fitfun = new FitnessFunction();
  470.         final double[] guess = getStartPoint();
  471.         // number of objective variables/problem dimension
  472.         dimension = guess.length;
  473.         initializeCMA(guess);
  474.         iterations = 0;
  475.         double bestValue = fitfun.value(guess);
  476.         push(fitnessHistory, bestValue);
  477.         PointValuePair optimum = new PointValuePair(getStartPoint(),
  478.                 isMinimize ? bestValue : -bestValue);
  479.         PointValuePair lastResult = null;

  480.         // -------------------- Generation Loop --------------------------------

  481.         generationLoop:
  482.         for (iterations = 1; iterations <= maxIterations; iterations++) {
  483.             // Generate and evaluate lambda offspring
  484.             final RealMatrix arz = randn1(dimension, lambda);
  485.             final RealMatrix arx = zeros(dimension, lambda);
  486.             final double[] fitness = new double[lambda];
  487.             // generate random offspring
  488.             for (int k = 0; k < lambda; k++) {
  489.                 RealMatrix arxk = null;
  490.                 for (int i = 0; i < checkFeasableCount + 1; i++) {
  491.                     if (diagonalOnly <= 0) {
  492.                         arxk = xmean.add(BD.multiply(arz.getColumnMatrix(k))
  493.                                          .scalarMultiply(sigma)); // m + sig * Normal(0,C)
  494.                     } else {
  495.                         arxk = xmean.add(times(diagD,arz.getColumnMatrix(k))
  496.                                          .scalarMultiply(sigma));
  497.                     }
  498.                     if (i >= checkFeasableCount ||
  499.                         fitfun.isFeasible(arxk.getColumn(0))) {
  500.                         break;
  501.                     }
  502.                     // regenerate random arguments for row
  503.                     arz.setColumn(k, randn(dimension));
  504.                 }
  505.                 copyColumn(arxk, 0, arx, k);
  506.                 try {
  507.                     fitness[k] = fitfun.value(arx.getColumn(k)); // compute fitness
  508.                 } catch (TooManyEvaluationsException e) {
  509.                     break generationLoop;
  510.                 }
  511.             }
  512.             // Sort by fitness and compute weighted mean into xmean
  513.             final int[] arindex = sortedIndices(fitness);
  514.             // Calculate new xmean, this is selection and recombination
  515.             final RealMatrix xold = xmean; // for speed up of Eq. (2) and (3)
  516.             final RealMatrix bestArx = selectColumns(arx, MathArrays.copyOf(arindex, mu));
  517.             xmean = bestArx.multiply(weights);
  518.             final RealMatrix bestArz = selectColumns(arz, MathArrays.copyOf(arindex, mu));
  519.             final RealMatrix zmean = bestArz.multiply(weights);
  520.             final boolean hsig = updateEvolutionPaths(zmean, xold);
  521.             if (diagonalOnly <= 0) {
  522.                 updateCovariance(hsig, bestArx, arz, arindex, xold);
  523.             } else {
  524.                 updateCovarianceDiagonalOnly(hsig, bestArz);
  525.             }
  526.             // Adapt step size sigma - Eq. (5)
  527.             sigma *= FastMath.exp(FastMath.min(1, (normps/chiN - 1) * cs / damps));
  528.             final double bestFitness = fitness[arindex[0]];
  529.             final double worstFitness = fitness[arindex[arindex.length - 1]];
  530.             if (bestValue > bestFitness) {
  531.                 bestValue = bestFitness;
  532.                 lastResult = optimum;
  533.                 optimum = new PointValuePair(fitfun.repair(bestArx.getColumn(0)),
  534.                                              isMinimize ? bestFitness : -bestFitness);
  535.                 if (getConvergenceChecker() != null && lastResult != null &&
  536.                     getConvergenceChecker().converged(iterations, optimum, lastResult)) {
  537.                     break generationLoop;
  538.                 }
  539.             }
  540.             // handle termination criteria
  541.             // Break, if fitness is good enough
  542.             if (stopFitness != 0 && bestFitness < (isMinimize ? stopFitness : -stopFitness)) {
  543.                 break generationLoop;
  544.             }
  545.             final double[] sqrtDiagC = sqrt(diagC).getColumn(0);
  546.             final double[] pcCol = pc.getColumn(0);
  547.             for (int i = 0; i < dimension; i++) {
  548.                 if (sigma * FastMath.max(FastMath.abs(pcCol[i]), sqrtDiagC[i]) > stopTolX) {
  549.                     break;
  550.                 }
  551.                 if (i >= dimension - 1) {
  552.                     break generationLoop;
  553.                 }
  554.             }
  555.             for (int i = 0; i < dimension; i++) {
  556.                 if (sigma * sqrtDiagC[i] > stopTolUpX) {
  557.                     break generationLoop;
  558.                 }
  559.             }
  560.             final double historyBest = min(fitnessHistory);
  561.             final double historyWorst = max(fitnessHistory);
  562.             if (iterations > 2 &&
  563.                 FastMath.max(historyWorst, worstFitness) -
  564.                 FastMath.min(historyBest, bestFitness) < stopTolFun) {
  565.                 break generationLoop;
  566.             }
  567.             if (iterations > fitnessHistory.length &&
  568.                 historyWorst-historyBest < stopTolHistFun) {
  569.                 break generationLoop;
  570.             }
  571.             // condition number of the covariance matrix exceeds 1e14
  572.             if (max(diagD)/min(diagD) > 1e7) {
  573.                 break generationLoop;
  574.             }
  575.             // user defined termination
  576.             if (getConvergenceChecker() != null) {
  577.                 final PointValuePair current
  578.                     = new PointValuePair(bestArx.getColumn(0),
  579.                                          isMinimize ? bestFitness : -bestFitness);
  580.                 if (lastResult != null &&
  581.                     getConvergenceChecker().converged(iterations, current, lastResult)) {
  582.                     break generationLoop;
  583.                     }
  584.                 lastResult = current;
  585.             }
  586.             // Adjust step size in case of equal function values (flat fitness)
  587.             if (bestValue == fitness[arindex[(int)(0.1+lambda/4.)]]) {
  588.                 sigma *= FastMath.exp(0.2 + cs / damps);
  589.             }
  590.             if (iterations > 2 && FastMath.max(historyWorst, bestFitness) -
  591.                 FastMath.min(historyBest, bestFitness) == 0) {
  592.                 sigma *= FastMath.exp(0.2 + cs / damps);
  593.             }
  594.             // store best in history
  595.             push(fitnessHistory,bestFitness);
  596.             fitfun.setValueRange(worstFitness-bestFitness);
  597.             if (generateStatistics) {
  598.                 statisticsSigmaHistory.add(sigma);
  599.                 statisticsFitnessHistory.add(bestFitness);
  600.                 statisticsMeanHistory.add(xmean.transpose());
  601.                 statisticsDHistory.add(diagD.transpose().scalarMultiply(1E5));
  602.             }
  603.         }
  604.         return optimum;
  605.     }

  606.     /**
  607.      * Scans the list of (required and optional) optimization data that
  608.      * characterize the problem.
  609.      *
  610.      * @param optData Optimization data. The following data will be looked for:
  611.      * <ul>
  612.      *  <li>{@link Sigma}</li>
  613.      *  <li>{@link PopulationSize}</li>
  614.      * </ul>
  615.      */
  616.     private void parseOptimizationData(OptimizationData... optData) {
  617.         // The existing values (as set by the previous call) are reused if
  618.         // not provided in the argument list.
  619.         for (OptimizationData data : optData) {
  620.             if (data instanceof Sigma) {
  621.                 inputSigma = ((Sigma) data).getSigma();
  622.                 continue;
  623.             }
  624.             if (data instanceof PopulationSize) {
  625.                 lambda = ((PopulationSize) data).getPopulationSize();
  626.                 continue;
  627.             }
  628.         }
  629.     }

  630.     /**
  631.      * Checks dimensions and values of boundaries and inputSigma if defined.
  632.      */
  633.     private void checkParameters() {
  634.         final double[] init = getStartPoint();
  635.         final double[] lB = getLowerBound();
  636.         final double[] uB = getUpperBound();

  637.         if (inputSigma != null) {
  638.             if (inputSigma.length != init.length) {
  639.                 throw new DimensionMismatchException(inputSigma.length, init.length);
  640.             }
  641.             for (int i = 0; i < init.length; i++) {
  642.                 if (inputSigma[i] < 0) {
  643.                     // XXX Remove this block in 4.0 (check performed in "Sigma" class).
  644.                     throw new NotPositiveException(inputSigma[i]);
  645.                 }
  646.                 if (inputSigma[i] > uB[i] - lB[i]) {
  647.                     throw new OutOfRangeException(inputSigma[i], 0, uB[i] - lB[i]);
  648.                 }
  649.             }
  650.         }
  651.     }

  652.     /**
  653.      * Initialization of the dynamic search parameters
  654.      *
  655.      * @param guess Initial guess for the arguments of the fitness function.
  656.      */
  657.     private void initializeCMA(double[] guess) {
  658.         if (lambda <= 0) {
  659.             // XXX Line below to replace the current one in 4.0 (MATH-879).
  660.             // throw new NotStrictlyPositiveException(lambda);
  661.             lambda = 4 + (int) (3 * FastMath.log(dimension));
  662.         }
  663.         // initialize sigma
  664.         final double[][] sigmaArray = new double[guess.length][1];
  665.         for (int i = 0; i < guess.length; i++) {
  666.             // XXX Line below to replace the current one in 4.0 (MATH-868).
  667.             // sigmaArray[i][0] = inputSigma[i];
  668.             sigmaArray[i][0] = inputSigma == null ? 0.3 : inputSigma[i];
  669.         }
  670.         final RealMatrix insigma = new Array2DRowRealMatrix(sigmaArray, false);
  671.         sigma = max(insigma); // overall standard deviation

  672.         // initialize termination criteria
  673.         stopTolUpX = 1e3 * max(insigma);
  674.         stopTolX = 1e-11 * max(insigma);
  675.         stopTolFun = 1e-12;
  676.         stopTolHistFun = 1e-13;

  677.         // initialize selection strategy parameters
  678.         mu = lambda / 2; // number of parents/points for recombination
  679.         logMu2 = FastMath.log(mu + 0.5);
  680.         weights = log(sequence(1, mu, 1)).scalarMultiply(-1).scalarAdd(logMu2);
  681.         double sumw = 0;
  682.         double sumwq = 0;
  683.         for (int i = 0; i < mu; i++) {
  684.             double w = weights.getEntry(i, 0);
  685.             sumw += w;
  686.             sumwq += w * w;
  687.         }
  688.         weights = weights.scalarMultiply(1 / sumw);
  689.         mueff = sumw * sumw / sumwq; // variance-effectiveness of sum w_i x_i

  690.         // initialize dynamic strategy parameters and constants
  691.         cc = (4 + mueff / dimension) /
  692.                 (dimension + 4 + 2 * mueff / dimension);
  693.         cs = (mueff + 2) / (dimension + mueff + 3.);
  694.         damps = (1 + 2 * FastMath.max(0, FastMath.sqrt((mueff - 1) /
  695.                                                        (dimension + 1)) - 1)) *
  696.             FastMath.max(0.3,
  697.                          1 - dimension / (1e-6 + maxIterations)) + cs; // minor increment
  698.         ccov1 = 2 / ((dimension + 1.3) * (dimension + 1.3) + mueff);
  699.         ccovmu = FastMath.min(1 - ccov1, 2 * (mueff - 2 + 1 / mueff) /
  700.                               ((dimension + 2) * (dimension + 2) + mueff));
  701.         ccov1Sep = FastMath.min(1, ccov1 * (dimension + 1.5) / 3);
  702.         ccovmuSep = FastMath.min(1 - ccov1, ccovmu * (dimension + 1.5) / 3);
  703.         chiN = FastMath.sqrt(dimension) *
  704.             (1 - 1 / ((double) 4 * dimension) + 1 / ((double) 21 * dimension * dimension));
  705.         // intialize CMA internal values - updated each generation
  706.         xmean = MatrixUtils.createColumnRealMatrix(guess); // objective variables
  707.         diagD = insigma.scalarMultiply(1 / sigma);
  708.         diagC = square(diagD);
  709.         pc = zeros(dimension, 1); // evolution paths for C and sigma
  710.         ps = zeros(dimension, 1); // B defines the coordinate system
  711.         normps = ps.getFrobeniusNorm();

  712.         B = eye(dimension, dimension);
  713.         D = ones(dimension, 1); // diagonal D defines the scaling
  714.         BD = times(B, repmat(diagD.transpose(), dimension, 1));
  715.         C = B.multiply(diag(square(D)).multiply(B.transpose())); // covariance
  716.         historySize = 10 + (int) (3 * 10 * dimension / (double) lambda);
  717.         fitnessHistory = new double[historySize]; // history of fitness values
  718.         for (int i = 0; i < historySize; i++) {
  719.             fitnessHistory[i] = Double.MAX_VALUE;
  720.         }
  721.     }

  722.     /**
  723.      * Update of the evolution paths ps and pc.
  724.      *
  725.      * @param zmean Weighted row matrix of the gaussian random numbers generating
  726.      * the current offspring.
  727.      * @param xold xmean matrix of the previous generation.
  728.      * @return hsig flag indicating a small correction.
  729.      */
  730.     private boolean updateEvolutionPaths(RealMatrix zmean, RealMatrix xold) {
  731.         ps = ps.scalarMultiply(1 - cs).add(
  732.                 B.multiply(zmean).scalarMultiply(FastMath.sqrt(cs * (2 - cs) * mueff)));
  733.         normps = ps.getFrobeniusNorm();
  734.         final boolean hsig = normps /
  735.             FastMath.sqrt(1 - FastMath.pow(1 - cs, 2 * iterations)) /
  736.             chiN < 1.4 + 2 / ((double) dimension + 1);
  737.         pc = pc.scalarMultiply(1 - cc);
  738.         if (hsig) {
  739.             pc = pc.add(xmean.subtract(xold).scalarMultiply(FastMath.sqrt(cc * (2 - cc) * mueff) / sigma));
  740.         }
  741.         return hsig;
  742.     }

  743.     /**
  744.      * Update of the covariance matrix C for diagonalOnly > 0
  745.      *
  746.      * @param hsig Flag indicating a small correction.
  747.      * @param bestArz Fitness-sorted matrix of the gaussian random values of the
  748.      * current offspring.
  749.      */
  750.     private void updateCovarianceDiagonalOnly(boolean hsig,
  751.                                               final RealMatrix bestArz) {
  752.         // minor correction if hsig==false
  753.         double oldFac = hsig ? 0 : ccov1Sep * cc * (2 - cc);
  754.         oldFac += 1 - ccov1Sep - ccovmuSep;
  755.         diagC = diagC.scalarMultiply(oldFac) // regard old matrix
  756.             .add(square(pc).scalarMultiply(ccov1Sep)) // plus rank one update
  757.             .add((times(diagC, square(bestArz).multiply(weights))) // plus rank mu update
  758.                  .scalarMultiply(ccovmuSep));
  759.         diagD = sqrt(diagC); // replaces eig(C)
  760.         if (diagonalOnly > 1 &&
  761.             iterations > diagonalOnly) {
  762.             // full covariance matrix from now on
  763.             diagonalOnly = 0;
  764.             B = eye(dimension, dimension);
  765.             BD = diag(diagD);
  766.             C = diag(diagC);
  767.         }
  768.     }

  769.     /**
  770.      * Update of the covariance matrix C.
  771.      *
  772.      * @param hsig Flag indicating a small correction.
  773.      * @param bestArx Fitness-sorted matrix of the argument vectors producing the
  774.      * current offspring.
  775.      * @param arz Unsorted matrix containing the gaussian random values of the
  776.      * current offspring.
  777.      * @param arindex Indices indicating the fitness-order of the current offspring.
  778.      * @param xold xmean matrix of the previous generation.
  779.      */
  780.     private void updateCovariance(boolean hsig, final RealMatrix bestArx,
  781.                                   final RealMatrix arz, final int[] arindex,
  782.                                   final RealMatrix xold) {
  783.         double negccov = 0;
  784.         if (ccov1 + ccovmu > 0) {
  785.             final RealMatrix arpos = bestArx.subtract(repmat(xold, 1, mu))
  786.                 .scalarMultiply(1 / sigma); // mu difference vectors
  787.             final RealMatrix roneu = pc.multiply(pc.transpose())
  788.                 .scalarMultiply(ccov1); // rank one update
  789.             // minor correction if hsig==false
  790.             double oldFac = hsig ? 0 : ccov1 * cc * (2 - cc);
  791.             oldFac += 1 - ccov1 - ccovmu;
  792.             if (isActiveCMA) {
  793.                 // Adapt covariance matrix C active CMA
  794.                 negccov = (1 - ccovmu) * 0.25 * mueff / (FastMath.pow(dimension + 2, 1.5) + 2 * mueff);
  795.                 // keep at least 0.66 in all directions, small popsize are most
  796.                 // critical
  797.                 final double negminresidualvariance = 0.66;
  798.                 // where to make up for the variance loss
  799.                 final double negalphaold = 0.5;
  800.                 // prepare vectors, compute negative updating matrix Cneg
  801.                 final int[] arReverseIndex = reverse(arindex);
  802.                 RealMatrix arzneg = selectColumns(arz, MathArrays.copyOf(arReverseIndex, mu));
  803.                 RealMatrix arnorms = sqrt(sumRows(square(arzneg)));
  804.                 final int[] idxnorms = sortedIndices(arnorms.getRow(0));
  805.                 final RealMatrix arnormsSorted = selectColumns(arnorms, idxnorms);
  806.                 final int[] idxReverse = reverse(idxnorms);
  807.                 final RealMatrix arnormsReverse = selectColumns(arnorms, idxReverse);
  808.                 arnorms = divide(arnormsReverse, arnormsSorted);
  809.                 final int[] idxInv = inverse(idxnorms);
  810.                 final RealMatrix arnormsInv = selectColumns(arnorms, idxInv);
  811.                 // check and set learning rate negccov
  812.                 final double negcovMax = (1 - negminresidualvariance) /
  813.                     square(arnormsInv).multiply(weights).getEntry(0, 0);
  814.                 if (negccov > negcovMax) {
  815.                     negccov = negcovMax;
  816.                 }
  817.                 arzneg = times(arzneg, repmat(arnormsInv, dimension, 1));
  818.                 final RealMatrix artmp = BD.multiply(arzneg);
  819.                 final RealMatrix Cneg = artmp.multiply(diag(weights)).multiply(artmp.transpose());
  820.                 oldFac += negalphaold * negccov;
  821.                 C = C.scalarMultiply(oldFac)
  822.                     .add(roneu) // regard old matrix
  823.                     .add(arpos.scalarMultiply( // plus rank one update
  824.                                               ccovmu + (1 - negalphaold) * negccov) // plus rank mu update
  825.                          .multiply(times(repmat(weights, 1, dimension),
  826.                                          arpos.transpose())))
  827.                     .subtract(Cneg.scalarMultiply(negccov));
  828.             } else {
  829.                 // Adapt covariance matrix C - nonactive
  830.                 C = C.scalarMultiply(oldFac) // regard old matrix
  831.                     .add(roneu) // plus rank one update
  832.                     .add(arpos.scalarMultiply(ccovmu) // plus rank mu update
  833.                          .multiply(times(repmat(weights, 1, dimension),
  834.                                          arpos.transpose())));
  835.             }
  836.         }
  837.         updateBD(negccov);
  838.     }

  839.     /**
  840.      * Update B and D from C.
  841.      *
  842.      * @param negccov Negative covariance factor.
  843.      */
  844.     private void updateBD(double negccov) {
  845.         if (ccov1 + ccovmu + negccov > 0 &&
  846.             (iterations % 1. / (ccov1 + ccovmu + negccov) / dimension / 10.) < 1) {
  847.             // to achieve O(N^2)
  848.             C = triu(C, 0).add(triu(C, 1).transpose());
  849.             // enforce symmetry to prevent complex numbers
  850.             final EigenDecomposition eig = new EigenDecomposition(C);
  851.             B = eig.getV(); // eigen decomposition, B==normalized eigenvectors
  852.             D = eig.getD();
  853.             diagD = diag(D);
  854.             if (min(diagD) <= 0) {
  855.                 for (int i = 0; i < dimension; i++) {
  856.                     if (diagD.getEntry(i, 0) < 0) {
  857.                         diagD.setEntry(i, 0, 0);
  858.                     }
  859.                 }
  860.                 final double tfac = max(diagD) / 1e14;
  861.                 C = C.add(eye(dimension, dimension).scalarMultiply(tfac));
  862.                 diagD = diagD.add(ones(dimension, 1).scalarMultiply(tfac));
  863.             }
  864.             if (max(diagD) > 1e14 * min(diagD)) {
  865.                 final double tfac = max(diagD) / 1e14 - min(diagD);
  866.                 C = C.add(eye(dimension, dimension).scalarMultiply(tfac));
  867.                 diagD = diagD.add(ones(dimension, 1).scalarMultiply(tfac));
  868.             }
  869.             diagC = diag(C);
  870.             diagD = sqrt(diagD); // D contains standard deviations now
  871.             BD = times(B, repmat(diagD.transpose(), dimension, 1)); // O(n^2)
  872.         }
  873.     }

  874.     /**
  875.      * Pushes the current best fitness value in a history queue.
  876.      *
  877.      * @param vals History queue.
  878.      * @param val Current best fitness value.
  879.      */
  880.     private static void push(double[] vals, double val) {
  881.         for (int i = vals.length-1; i > 0; i--) {
  882.             vals[i] = vals[i-1];
  883.         }
  884.         vals[0] = val;
  885.     }

  886.     /**
  887.      * Sorts fitness values.
  888.      *
  889.      * @param doubles Array of values to be sorted.
  890.      * @return a sorted array of indices pointing into doubles.
  891.      */
  892.     private int[] sortedIndices(final double[] doubles) {
  893.         final DoubleIndex[] dis = new DoubleIndex[doubles.length];
  894.         for (int i = 0; i < doubles.length; i++) {
  895.             dis[i] = new DoubleIndex(doubles[i], i);
  896.         }
  897.         Arrays.sort(dis);
  898.         final int[] indices = new int[doubles.length];
  899.         for (int i = 0; i < doubles.length; i++) {
  900.             indices[i] = dis[i].index;
  901.         }
  902.         return indices;
  903.     }

  904.     /**
  905.      * Used to sort fitness values. Sorting is always in lower value first
  906.      * order.
  907.      */
  908.     private static class DoubleIndex implements Comparable<DoubleIndex> {
  909.         /** Value to compare. */
  910.         private final double value;
  911.         /** Index into sorted array. */
  912.         private final int index;

  913.         /**
  914.          * @param value Value to compare.
  915.          * @param index Index into sorted array.
  916.          */
  917.         DoubleIndex(double value, int index) {
  918.             this.value = value;
  919.             this.index = index;
  920.         }

  921.         /** {@inheritDoc} */
  922.         public int compareTo(DoubleIndex o) {
  923.             return Double.compare(value, o.value);
  924.         }

  925.         /** {@inheritDoc} */
  926.         @Override
  927.         public boolean equals(Object other) {

  928.             if (this == other) {
  929.                 return true;
  930.             }

  931.             if (other instanceof DoubleIndex) {
  932.                 return Double.compare(value, ((DoubleIndex) other).value) == 0;
  933.             }

  934.             return false;
  935.         }

  936.         /** {@inheritDoc} */
  937.         @Override
  938.         public int hashCode() {
  939.             long bits = Double.doubleToLongBits(value);
  940.             return (int) ((1438542 ^ (bits >>> 32) ^ bits) & 0xffffffff);
  941.         }
  942.     }

  943.     /**
  944.      * Normalizes fitness values to the range [0,1]. Adds a penalty to the
  945.      * fitness value if out of range. The penalty is adjusted by calling
  946.      * setValueRange().
  947.      */
  948.     private class FitnessFunction {
  949.         /** Determines the penalty for boundary violations */
  950.         private double valueRange;
  951.         /**
  952.          * Flag indicating whether the objective variables are forced into their
  953.          * bounds if defined
  954.          */
  955.         private final boolean isRepairMode;

  956.         /** Simple constructor.
  957.          */
  958.         public FitnessFunction() {
  959.             valueRange = 1;
  960.             isRepairMode = true;
  961.         }

  962.         /**
  963.          * @param point Normalized objective variables.
  964.          * @return the objective value + penalty for violated bounds.
  965.          */
  966.         public double value(final double[] point) {
  967.             double value;
  968.             if (isRepairMode) {
  969.                 double[] repaired = repair(point);
  970.                 value = CMAESOptimizer.this.computeObjectiveValue(repaired) +
  971.                     penalty(point, repaired);
  972.             } else {
  973.                 value = CMAESOptimizer.this.computeObjectiveValue(point);
  974.             }
  975.             return isMinimize ? value : -value;
  976.         }

  977.         /**
  978.          * @param x Normalized objective variables.
  979.          * @return {@code true} if in bounds.
  980.          */
  981.         public boolean isFeasible(final double[] x) {
  982.             final double[] lB = CMAESOptimizer.this.getLowerBound();
  983.             final double[] uB = CMAESOptimizer.this.getUpperBound();

  984.             for (int i = 0; i < x.length; i++) {
  985.                 if (x[i] < lB[i]) {
  986.                     return false;
  987.                 }
  988.                 if (x[i] > uB[i]) {
  989.                     return false;
  990.                 }
  991.             }
  992.             return true;
  993.         }

  994.         /**
  995.          * @param valueRange Adjusts the penalty computation.
  996.          */
  997.         public void setValueRange(double valueRange) {
  998.             this.valueRange = valueRange;
  999.         }

  1000.         /**
  1001.          * @param x Normalized objective variables.
  1002.          * @return the repaired (i.e. all in bounds) objective variables.
  1003.          */
  1004.         private double[] repair(final double[] x) {
  1005.             final double[] lB = CMAESOptimizer.this.getLowerBound();
  1006.             final double[] uB = CMAESOptimizer.this.getUpperBound();

  1007.             final double[] repaired = new double[x.length];
  1008.             for (int i = 0; i < x.length; i++) {
  1009.                 if (x[i] < lB[i]) {
  1010.                     repaired[i] = lB[i];
  1011.                 } else if (x[i] > uB[i]) {
  1012.                     repaired[i] = uB[i];
  1013.                 } else {
  1014.                     repaired[i] = x[i];
  1015.                 }
  1016.             }
  1017.             return repaired;
  1018.         }

  1019.         /**
  1020.          * @param x Normalized objective variables.
  1021.          * @param repaired Repaired objective variables.
  1022.          * @return Penalty value according to the violation of the bounds.
  1023.          */
  1024.         private double penalty(final double[] x, final double[] repaired) {
  1025.             double penalty = 0;
  1026.             for (int i = 0; i < x.length; i++) {
  1027.                 double diff = FastMath.abs(x[i] - repaired[i]);
  1028.                 penalty += diff * valueRange;
  1029.             }
  1030.             return isMinimize ? penalty : -penalty;
  1031.         }
  1032.     }

  1033.     // -----Matrix utility functions similar to the Matlab build in functions------

  1034.     /**
  1035.      * @param m Input matrix
  1036.      * @return Matrix representing the element-wise logarithm of m.
  1037.      */
  1038.     private static RealMatrix log(final RealMatrix m) {
  1039.         final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
  1040.         for (int r = 0; r < m.getRowDimension(); r++) {
  1041.             for (int c = 0; c < m.getColumnDimension(); c++) {
  1042.                 d[r][c] = FastMath.log(m.getEntry(r, c));
  1043.             }
  1044.         }
  1045.         return new Array2DRowRealMatrix(d, false);
  1046.     }

  1047.     /**
  1048.      * @param m Input matrix.
  1049.      * @return Matrix representing the element-wise square root of m.
  1050.      */
  1051.     private static RealMatrix sqrt(final RealMatrix m) {
  1052.         final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
  1053.         for (int r = 0; r < m.getRowDimension(); r++) {
  1054.             for (int c = 0; c < m.getColumnDimension(); c++) {
  1055.                 d[r][c] = FastMath.sqrt(m.getEntry(r, c));
  1056.             }
  1057.         }
  1058.         return new Array2DRowRealMatrix(d, false);
  1059.     }

  1060.     /**
  1061.      * @param m Input matrix.
  1062.      * @return Matrix representing the element-wise square of m.
  1063.      */
  1064.     private static RealMatrix square(final RealMatrix m) {
  1065.         final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
  1066.         for (int r = 0; r < m.getRowDimension(); r++) {
  1067.             for (int c = 0; c < m.getColumnDimension(); c++) {
  1068.                 double e = m.getEntry(r, c);
  1069.                 d[r][c] = e * e;
  1070.             }
  1071.         }
  1072.         return new Array2DRowRealMatrix(d, false);
  1073.     }

  1074.     /**
  1075.      * @param m Input matrix 1.
  1076.      * @param n Input matrix 2.
  1077.      * @return the matrix where the elements of m and n are element-wise multiplied.
  1078.      */
  1079.     private static RealMatrix times(final RealMatrix m, final RealMatrix n) {
  1080.         final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
  1081.         for (int r = 0; r < m.getRowDimension(); r++) {
  1082.             for (int c = 0; c < m.getColumnDimension(); c++) {
  1083.                 d[r][c] = m.getEntry(r, c) * n.getEntry(r, c);
  1084.             }
  1085.         }
  1086.         return new Array2DRowRealMatrix(d, false);
  1087.     }

  1088.     /**
  1089.      * @param m Input matrix 1.
  1090.      * @param n Input matrix 2.
  1091.      * @return Matrix where the elements of m and n are element-wise divided.
  1092.      */
  1093.     private static RealMatrix divide(final RealMatrix m, final RealMatrix n) {
  1094.         final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
  1095.         for (int r = 0; r < m.getRowDimension(); r++) {
  1096.             for (int c = 0; c < m.getColumnDimension(); c++) {
  1097.                 d[r][c] = m.getEntry(r, c) / n.getEntry(r, c);
  1098.             }
  1099.         }
  1100.         return new Array2DRowRealMatrix(d, false);
  1101.     }

  1102.     /**
  1103.      * @param m Input matrix.
  1104.      * @param cols Columns to select.
  1105.      * @return Matrix representing the selected columns.
  1106.      */
  1107.     private static RealMatrix selectColumns(final RealMatrix m, final int[] cols) {
  1108.         final double[][] d = new double[m.getRowDimension()][cols.length];
  1109.         for (int r = 0; r < m.getRowDimension(); r++) {
  1110.             for (int c = 0; c < cols.length; c++) {
  1111.                 d[r][c] = m.getEntry(r, cols[c]);
  1112.             }
  1113.         }
  1114.         return new Array2DRowRealMatrix(d, false);
  1115.     }

  1116.     /**
  1117.      * @param m Input matrix.
  1118.      * @param k Diagonal position.
  1119.      * @return Upper triangular part of matrix.
  1120.      */
  1121.     private static RealMatrix triu(final RealMatrix m, int k) {
  1122.         final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
  1123.         for (int r = 0; r < m.getRowDimension(); r++) {
  1124.             for (int c = 0; c < m.getColumnDimension(); c++) {
  1125.                 d[r][c] = r <= c - k ? m.getEntry(r, c) : 0;
  1126.             }
  1127.         }
  1128.         return new Array2DRowRealMatrix(d, false);
  1129.     }

  1130.     /**
  1131.      * @param m Input matrix.
  1132.      * @return Row matrix representing the sums of the rows.
  1133.      */
  1134.     private static RealMatrix sumRows(final RealMatrix m) {
  1135.         final double[][] d = new double[1][m.getColumnDimension()];
  1136.         for (int c = 0; c < m.getColumnDimension(); c++) {
  1137.             double sum = 0;
  1138.             for (int r = 0; r < m.getRowDimension(); r++) {
  1139.                 sum += m.getEntry(r, c);
  1140.             }
  1141.             d[0][c] = sum;
  1142.         }
  1143.         return new Array2DRowRealMatrix(d, false);
  1144.     }

  1145.     /**
  1146.      * @param m Input matrix.
  1147.      * @return the diagonal n-by-n matrix if m is a column matrix or the column
  1148.      * matrix representing the diagonal if m is a n-by-n matrix.
  1149.      */
  1150.     private static RealMatrix diag(final RealMatrix m) {
  1151.         if (m.getColumnDimension() == 1) {
  1152.             final double[][] d = new double[m.getRowDimension()][m.getRowDimension()];
  1153.             for (int i = 0; i < m.getRowDimension(); i++) {
  1154.                 d[i][i] = m.getEntry(i, 0);
  1155.             }
  1156.             return new Array2DRowRealMatrix(d, false);
  1157.         } else {
  1158.             final double[][] d = new double[m.getRowDimension()][1];
  1159.             for (int i = 0; i < m.getColumnDimension(); i++) {
  1160.                 d[i][0] = m.getEntry(i, i);
  1161.             }
  1162.             return new Array2DRowRealMatrix(d, false);
  1163.         }
  1164.     }

  1165.     /**
  1166.      * Copies a column from m1 to m2.
  1167.      *
  1168.      * @param m1 Source matrix.
  1169.      * @param col1 Source column.
  1170.      * @param m2 Target matrix.
  1171.      * @param col2 Target column.
  1172.      */
  1173.     private static void copyColumn(final RealMatrix m1, int col1,
  1174.                                    RealMatrix m2, int col2) {
  1175.         for (int i = 0; i < m1.getRowDimension(); i++) {
  1176.             m2.setEntry(i, col2, m1.getEntry(i, col1));
  1177.         }
  1178.     }

  1179.     /**
  1180.      * @param n Number of rows.
  1181.      * @param m Number of columns.
  1182.      * @return n-by-m matrix filled with 1.
  1183.      */
  1184.     private static RealMatrix ones(int n, int m) {
  1185.         final double[][] d = new double[n][m];
  1186.         for (int r = 0; r < n; r++) {
  1187.             Arrays.fill(d[r], 1);
  1188.         }
  1189.         return new Array2DRowRealMatrix(d, false);
  1190.     }

  1191.     /**
  1192.      * @param n Number of rows.
  1193.      * @param m Number of columns.
  1194.      * @return n-by-m matrix of 0 values out of diagonal, and 1 values on
  1195.      * the diagonal.
  1196.      */
  1197.     private static RealMatrix eye(int n, int m) {
  1198.         final double[][] d = new double[n][m];
  1199.         for (int r = 0; r < n; r++) {
  1200.             if (r < m) {
  1201.                 d[r][r] = 1;
  1202.             }
  1203.         }
  1204.         return new Array2DRowRealMatrix(d, false);
  1205.     }

  1206.     /**
  1207.      * @param n Number of rows.
  1208.      * @param m Number of columns.
  1209.      * @return n-by-m matrix of zero values.
  1210.      */
  1211.     private static RealMatrix zeros(int n, int m) {
  1212.         return new Array2DRowRealMatrix(n, m);
  1213.     }

  1214.     /**
  1215.      * @param mat Input matrix.
  1216.      * @param n Number of row replicates.
  1217.      * @param m Number of column replicates.
  1218.      * @return a matrix which replicates the input matrix in both directions.
  1219.      */
  1220.     private static RealMatrix repmat(final RealMatrix mat, int n, int m) {
  1221.         final int rd = mat.getRowDimension();
  1222.         final int cd = mat.getColumnDimension();
  1223.         final double[][] d = new double[n * rd][m * cd];
  1224.         for (int r = 0; r < n * rd; r++) {
  1225.             for (int c = 0; c < m * cd; c++) {
  1226.                 d[r][c] = mat.getEntry(r % rd, c % cd);
  1227.             }
  1228.         }
  1229.         return new Array2DRowRealMatrix(d, false);
  1230.     }

  1231.     /**
  1232.      * @param start Start value.
  1233.      * @param end End value.
  1234.      * @param step Step size.
  1235.      * @return a sequence as column matrix.
  1236.      */
  1237.     private static RealMatrix sequence(double start, double end, double step) {
  1238.         final int size = (int) ((end - start) / step + 1);
  1239.         final double[][] d = new double[size][1];
  1240.         double value = start;
  1241.         for (int r = 0; r < size; r++) {
  1242.             d[r][0] = value;
  1243.             value += step;
  1244.         }
  1245.         return new Array2DRowRealMatrix(d, false);
  1246.     }

  1247.     /**
  1248.      * @param m Input matrix.
  1249.      * @return the maximum of the matrix element values.
  1250.      */
  1251.     private static double max(final RealMatrix m) {
  1252.         double max = -Double.MAX_VALUE;
  1253.         for (int r = 0; r < m.getRowDimension(); r++) {
  1254.             for (int c = 0; c < m.getColumnDimension(); c++) {
  1255.                 double e = m.getEntry(r, c);
  1256.                 if (max < e) {
  1257.                     max = e;
  1258.                 }
  1259.             }
  1260.         }
  1261.         return max;
  1262.     }

  1263.     /**
  1264.      * @param m Input matrix.
  1265.      * @return the minimum of the matrix element values.
  1266.      */
  1267.     private static double min(final RealMatrix m) {
  1268.         double min = Double.MAX_VALUE;
  1269.         for (int r = 0; r < m.getRowDimension(); r++) {
  1270.             for (int c = 0; c < m.getColumnDimension(); c++) {
  1271.                 double e = m.getEntry(r, c);
  1272.                 if (min > e) {
  1273.                     min = e;
  1274.                 }
  1275.             }
  1276.         }
  1277.         return min;
  1278.     }

  1279.     /**
  1280.      * @param m Input array.
  1281.      * @return the maximum of the array values.
  1282.      */
  1283.     private static double max(final double[] m) {
  1284.         double max = -Double.MAX_VALUE;
  1285.         for (int r = 0; r < m.length; r++) {
  1286.             if (max < m[r]) {
  1287.                 max = m[r];
  1288.             }
  1289.         }
  1290.         return max;
  1291.     }

  1292.     /**
  1293.      * @param m Input array.
  1294.      * @return the minimum of the array values.
  1295.      */
  1296.     private static double min(final double[] m) {
  1297.         double min = Double.MAX_VALUE;
  1298.         for (int r = 0; r < m.length; r++) {
  1299.             if (min > m[r]) {
  1300.                 min = m[r];
  1301.             }
  1302.         }
  1303.         return min;
  1304.     }

  1305.     /**
  1306.      * @param indices Input index array.
  1307.      * @return the inverse of the mapping defined by indices.
  1308.      */
  1309.     private static int[] inverse(final int[] indices) {
  1310.         final int[] inverse = new int[indices.length];
  1311.         for (int i = 0; i < indices.length; i++) {
  1312.             inverse[indices[i]] = i;
  1313.         }
  1314.         return inverse;
  1315.     }

  1316.     /**
  1317.      * @param indices Input index array.
  1318.      * @return the indices in inverse order (last is first).
  1319.      */
  1320.     private static int[] reverse(final int[] indices) {
  1321.         final int[] reverse = new int[indices.length];
  1322.         for (int i = 0; i < indices.length; i++) {
  1323.             reverse[i] = indices[indices.length - i - 1];
  1324.         }
  1325.         return reverse;
  1326.     }

  1327.     /**
  1328.      * @param size Length of random array.
  1329.      * @return an array of Gaussian random numbers.
  1330.      */
  1331.     private double[] randn(int size) {
  1332.         final double[] randn = new double[size];
  1333.         for (int i = 0; i < size; i++) {
  1334.             randn[i] = random.nextGaussian();
  1335.         }
  1336.         return randn;
  1337.     }

  1338.     /**
  1339.      * @param size Number of rows.
  1340.      * @param popSize Population size.
  1341.      * @return a 2-dimensional matrix of Gaussian random numbers.
  1342.      */
  1343.     private RealMatrix randn1(int size, int popSize) {
  1344.         final double[][] d = new double[size][popSize];
  1345.         for (int r = 0; r < size; r++) {
  1346.             for (int c = 0; c < popSize; c++) {
  1347.                 d[r][c] = random.nextGaussian();
  1348.             }
  1349.         }
  1350.         return new Array2DRowRealMatrix(d, false);
  1351.     }
  1352. }