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.  */

  18. package org.apache.commons.math3.optim.nonlinear.scalar.noderiv;

  20. import java.util.ArrayList;
  21. import java.util.Arrays;
  22. import java.util.List;

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

  42. /**
  43.  * An implementation of the active Covariance Matrix Adaptation Evolution Strategy (CMA-ES)
  44.  * for non-linear, non-convex, non-smooth, global function minimization.
  45.  * <p>
  46.  * The CMA-Evolution Strategy (CMA-ES) is a reliable stochastic optimization method
  47.  * which should be applied if derivative-based methods, e.g. quasi-Newton BFGS or
  48.  * conjugate gradient, fail due to a rugged search landscape (e.g. noise, local
  49.  * optima, outlier, etc.) of the objective function. Like a
  50.  * quasi-Newton method, the CMA-ES learns and applies a variable metric
  51.  * on the underlying search space. Unlike a quasi-Newton method, the
  52.  * CMA-ES neither estimates nor uses gradients, making it considerably more
  53.  * reliable in terms of finding a good, or even close to optimal, solution.
  54.  * <p>
  55.  * In general, on smooth objective functions the CMA-ES is roughly ten times
  56.  * slower than BFGS (counting objective function evaluations, no gradients provided).
  57.  * For up to <math>N=10</math> variables also the derivative-free simplex
  58.  * direct search method (Nelder and Mead) can be faster, but it is
  59.  * far less reliable than CMA-ES.
  60.  * <p>
  61.  * The CMA-ES is particularly well suited for non-separable
  62.  * and/or badly conditioned problems. To observe the advantage of CMA compared
  63.  * to a conventional evolution strategy, it will usually take about
  64.  * <math>30 N</math> function evaluations. On difficult problems the complete
  65.  * optimization (a single run) is expected to take <em>roughly</em> between
  66.  * <math>30 N</math> and <math>300 N<sup>2</sup></math>
  67.  * function evaluations.
  68.  * <p>
  69.  * This implementation is translated and adapted from the Matlab version
  70.  * of the CMA-ES algorithm as implemented in module {@code cmaes.m} version 3.51.
  71.  * <p>
  72.  * For more information, please refer to the following links:
  73.  * <ul>
  74.  *  <li><a href="http://www.lri.fr/~hansen/cmaes.m">Matlab code</a></li>
  75.  *  <li><a href="http://www.lri.fr/~hansen/cmaesintro.html">Introduction to CMA-ES</a></li>
  76.  *  <li><a href="http://en.wikipedia.org/wiki/CMA-ES">Wikipedia</a></li>
  77.  * </ul>
  78.  *
  79.  * @since 3.0
  80.  */
  81. public class CMAESOptimizer
  82.     extends MultivariateOptimizer {
  83.     // global search parameters
  84.     /**
  85.      * Population size, offspring number. The primary strategy parameter to play
  86.      * with, which can be increased from its default value. Increasing the
  87.      * population size improves global search properties in exchange to speed.
  88.      * Speed decreases, as a rule, at most linearly with increasing population
  89.      * size. It is advisable to begin with the default small population size.
  90.      */
  91.     private int lambda; // population size
  92.     /**
  93.      * Covariance update mechanism, default is active CMA. isActiveCMA = true
  94.      * turns on "active CMA" with a negative update of the covariance matrix and
  95.      * checks for positive definiteness. OPTS.CMA.active = 2 does not check for
  96.      * pos. def. and is numerically faster. Active CMA usually speeds up the
  97.      * adaptation.
  98.      */
  99.     private final boolean isActiveCMA;
  100.     /**
  101.      * Determines how often a new random offspring is generated in case it is
  102.      * not feasible / beyond the defined limits, default is 0.
  103.      */
  104.     private final int checkFeasableCount;
  105.     /**
  106.      * @see Sigma
  107.      */
  108.     private double[] inputSigma;
  109.     /** Number of objective variables/problem dimension */
  110.     private int dimension;
  111.     /**
  112.      * Defines the number of initial iterations, where the covariance matrix
  113.      * remains diagonal and the algorithm has internally linear time complexity.
  114.      * diagonalOnly = 1 means keeping the covariance matrix always diagonal and
  115.      * this setting also exhibits linear space complexity. This can be
  116.      * particularly useful for dimension > 100.
  117.      * @see <a href="http://hal.archives-ouvertes.fr/inria-00287367/en">A Simple Modification in CMA-ES</a>
  118.      */
  119.     private int diagonalOnly;
  120.     /** Number of objective variables/problem dimension */
  121.     private boolean isMinimize = true;
  122.     /** Indicates whether statistic data is collected. */
  123.     private final boolean generateStatistics;

  125.     // termination criteria
  126.     /** Maximal number of iterations allowed. */
  127.     private final int maxIterations;
  128.     /** Limit for fitness value. */
  129.     private final double stopFitness;
  130.     /** Stop if x-changes larger stopTolUpX. */
  131.     private double stopTolUpX;
  132.     /** Stop if x-change smaller stopTolX. */
  133.     private double stopTolX;
  134.     /** Stop if fun-changes smaller stopTolFun. */
  135.     private double stopTolFun;
  136.     /** Stop if back fun-changes smaller stopTolHistFun. */
  137.     private double stopTolHistFun;

  139.     // selection strategy parameters
  140.     /** Number of parents/points for recombination. */
  141.     private int mu; //
  142.     /** log(mu + 0.5), stored for efficiency. */
  143.     private double logMu2;
  144.     /** Array for weighted recombination. */
  145.     private RealMatrix weights;
  146.     /** Variance-effectiveness of sum w_i x_i. */
  147.     private double mueff; //

  149.     // dynamic strategy parameters and constants
  150.     /** Overall standard deviation - search volume. */
  151.     private double sigma;
  152.     /** Cumulation constant. */
  153.     private double cc;
  154.     /** Cumulation constant for step-size. */
  155.     private double cs;
  156.     /** Damping for step-size. */
  157.     private double damps;
  158.     /** Learning rate for rank-one update. */
  159.     private double ccov1;
  160.     /** Learning rate for rank-mu update' */
  161.     private double ccovmu;
  162.     /** Expectation of ||N(0,I)|| == norm(randn(N,1)). */
  163.     private double chiN;
  164.     /** Learning rate for rank-one update - diagonalOnly */
  165.     private double ccov1Sep;
  166.     /** Learning rate for rank-mu update - diagonalOnly */
  167.     private double ccovmuSep;

  169.     // CMA internal values - updated each generation
  170.     /** Objective variables. */
  171.     private RealMatrix xmean;
  172.     /** Evolution path. */
  173.     private RealMatrix pc;
  174.     /** Evolution path for sigma. */
  175.     private RealMatrix ps;
  176.     /** Norm of ps, stored for efficiency. */
  177.     private double normps;
  178.     /** Coordinate system. */
  179.     private RealMatrix B;
  180.     /** Scaling. */
  181.     private RealMatrix D;
  182.     /** B*D, stored for efficiency. */
  183.     private RealMatrix BD;
  184.     /** Diagonal of sqrt(D), stored for efficiency. */
  185.     private RealMatrix diagD;
  186.     /** Covariance matrix. */
  187.     private RealMatrix C;
  188.     /** Diagonal of C, used for diagonalOnly. */
  189.     private RealMatrix diagC;
  190.     /** Number of iterations already performed. */
  191.     private int iterations;

  193.     /** History queue of best values. */
  194.     private double[] fitnessHistory;
  195.     /** Size of history queue of best values. */
  196.     private int historySize;

  198.     /** Random generator. */
  199.     private final RandomGenerator random;

  201.     /** History of sigma values. */
  202.     private final List<Double> statisticsSigmaHistory = new ArrayList<Double>();
  203.     /** History of mean matrix. */
  204.     private final List<RealMatrix> statisticsMeanHistory = new ArrayList<RealMatrix>();
  205.     /** History of fitness values. */
  206.     private final List<Double> statisticsFitnessHistory = new ArrayList<Double>();
  207.     /** History of D matrix. */
  208.     private final List<RealMatrix> statisticsDHistory = new ArrayList<RealMatrix>();

  210.     /**
  211.      * @param maxIterations Maximal number of iterations.
  212.      * @param stopFitness Whether to stop if objective function value is smaller than
  213.      * {@code stopFitness}.
  214.      * @param isActiveCMA Chooses the covariance matrix update method.
  215.      * @param diagonalOnly Number of initial iterations, where the covariance matrix
  216.      * remains diagonal.
  217.      * @param checkFeasableCount Determines how often new random objective variables are
  218.      * generated in case they are out of bounds.
  219.      * @param random Random generator.
  220.      * @param generateStatistics Whether statistic data is collected.
  221.      * @param checker Convergence checker.
  222.      *
  223.      * @since 3.1
  224.      */
  225.     public CMAESOptimizer(int maxIterations,
  226.                           double stopFitness,
  227.                           boolean isActiveCMA,
  228.                           int diagonalOnly,
  229.                           int checkFeasableCount,
  230.                           RandomGenerator random,
  231.                           boolean generateStatistics,
  232.                           ConvergenceChecker<PointValuePair> checker) {
  233.         super(checker);
  234.         this.maxIterations = maxIterations;
  235.         this.stopFitness = stopFitness;
  236.         this.isActiveCMA = isActiveCMA;
  237.         this.diagonalOnly = diagonalOnly;
  238.         this.checkFeasableCount = checkFeasableCount;
  239.         this.random = random;
  240.         this.generateStatistics = generateStatistics;
  241.     }

  243.     /**
  244.      * @return History of sigma values.
  245.      */
  246.     public List<Double> getStatisticsSigmaHistory() {
  247.         return statisticsSigmaHistory;
  248.     }

  250.     /**
  251.      * @return History of mean matrix.
  252.      */
  253.     public List<RealMatrix> getStatisticsMeanHistory() {
  254.         return statisticsMeanHistory;
  255.     }

  257.     /**
  258.      * @return History of fitness values.
  259.      */
  260.     public List<Double> getStatisticsFitnessHistory() {
  261.         return statisticsFitnessHistory;
  262.     }

  264.     /**
  265.      * @return History of D matrix.
  266.      */
  267.     public List<RealMatrix> getStatisticsDHistory() {
  268.         return statisticsDHistory;
  269.     }

  271.     /**
  272.      * Input sigma values.
  273.      * They define the initial coordinate-wise standard deviations for
  274.      * sampling new search points around the initial guess.
  275.      * It is suggested to set them to the estimated distance from the
  276.      * initial to the desired optimum.
  277.      * Small values induce the search to be more local (and very small
  278.      * values are more likely to find a local optimum close to the initial
  279.      * guess).
  280.      * Too small values might however lead to early termination.
  281.      */
  282.     public static class Sigma implements OptimizationData {
  283.         /** Sigma values. */
  284.         private final double[] sigma;

  286.         /**
  287.          * @param s Sigma values.
  288.          * @throws NotPositiveException if any of the array entries is smaller
  289.          * than zero.
  290.          */
  291.         public Sigma(double[] s)
  292.             throws NotPositiveException {
  293.             for (int i = 0; i < s.length; i++) {
  294.                 if (s[i] < 0) {
  295.                     throw new NotPositiveException(s[i]);
  296.                 }
  297.             }

  299.             sigma = s.clone();
  300.         }

  302.         /**
  303.          * @return the sigma values.
  304.          */
  305.         public double[] getSigma() {
  306.             return sigma.clone();
  307.         }
  308.     }

  310.     /**
  311.      * Population size.
  312.      * The number of offspring is the primary strategy parameter.
  313.      * In the absence of better clues, a good default could be an
  314.      * integer close to {@code 4 + 3 ln(n)}, where {@code n} is the
  315.      * number of optimized parameters.
  316.      * Increasing the population size improves global search properties
  317.      * at the expense of speed (which in general decreases at most
  318.      * linearly with increasing population size).
  319.      */
  320.     public static class PopulationSize implements OptimizationData {
  321.         /** Population size. */
  322.         private final int lambda;

  324.         /**
  325.          * @param size Population size.
  326.          * @throws NotStrictlyPositiveException if {@code size <= 0}.
  327.          */
  328.         public PopulationSize(int size)
  329.             throws NotStrictlyPositiveException {
  330.             if (size <= 0) {
  331.                 throw new NotStrictlyPositiveException(size);
  332.             }
  333.             lambda = size;
  334.         }

  336.         /**
  337.          * @return the population size.
  338.          */
  339.         public int getPopulationSize() {
  340.             return lambda;
  341.         }
  342.     }

  344.     /**
  345.      * {@inheritDoc}
  346.      *
  347.      * @param optData Optimization data. In addition to those documented in
  348.      * {@link MultivariateOptimizer#parseOptimizationData(OptimizationData[])
  349.      * MultivariateOptimizer}, this method will register the following data:
  350.      * <ul>
  351.      *  <li>{@link Sigma}</li>
  352.      *  <li>{@link PopulationSize}</li>
  353.      * </ul>
  354.      * @return {@inheritDoc}
  355.      * @throws TooManyEvaluationsException if the maximal number of
  356.      * evaluations is exceeded.
  357.      * @throws DimensionMismatchException if the initial guess, target, and weight
  358.      * arguments have inconsistent dimensions.
  359.      */
  360.     @Override
  361.     public PointValuePair optimize(OptimizationData... optData)
  362.         throws TooManyEvaluationsException,
  363.                DimensionMismatchException {
  364.         // Set up base class and perform computation.
  365.         return super.optimize(optData);
  366.     }

  368.     /** {@inheritDoc} */
  369.     @Override
  370.     protected PointValuePair doOptimize() {
  371.          // -------------------- Initialization --------------------------------
  372.         isMinimize = getGoalType().equals(GoalType.MINIMIZE);
  373.         final FitnessFunction fitfun = new FitnessFunction();
  374.         final double[] guess = getStartPoint();
  375.         // number of objective variables/problem dimension
  376.         dimension = guess.length;
  377.         initializeCMA(guess);
  378.         iterations = 0;
  379.         ValuePenaltyPair valuePenalty = fitfun.value(guess);
  380.         double bestValue = valuePenalty.value+valuePenalty.penalty;
  381.         push(fitnessHistory, bestValue);
  382.         PointValuePair optimum
  383.             = new PointValuePair(getStartPoint(),
  384.                                  isMinimize ? bestValue : -bestValue);
  385.         PointValuePair lastResult = null;

  387.         // -------------------- Generation Loop --------------------------------

  389.         generationLoop:
  390.         for (iterations = 1; iterations <= maxIterations; iterations++) {
  391.             incrementIterationCount();

  393.             // Generate and evaluate lambda offspring
  394.             final RealMatrix arz = randn1(dimension, lambda);
  395.             final RealMatrix arx = zeros(dimension, lambda);
  396.             final double[] fitness = new double[lambda];
  397.             final ValuePenaltyPair[] valuePenaltyPairs = new ValuePenaltyPair[lambda];
  398.             // generate random offspring
  399.             for (int k = 0; k < lambda; k++) {
  400.                 RealMatrix arxk = null;
  401.                 for (int i = 0; i < checkFeasableCount + 1; i++) {
  402.                     if (diagonalOnly <= 0) {
  403.                         arxk = xmean.add(BD.multiply(arz.getColumnMatrix(k))
  404.                                          .scalarMultiply(sigma)); // m + sig * Normal(0,C)
  405.                     } else {
  406.                         arxk = xmean.add(times(diagD,arz.getColumnMatrix(k))
  407.                                          .scalarMultiply(sigma));
  408.                     }
  409.                     if (i >= checkFeasableCount ||
  410.                         fitfun.isFeasible(arxk.getColumn(0))) {
  411.                         break;
  412.                     }
  413.                     // regenerate random arguments for row
  414.                     arz.setColumn(k, randn(dimension));
  415.                 }
  416.                 copyColumn(arxk, 0, arx, k);
  417.                 try {
  418.                     valuePenaltyPairs[k] = fitfun.value(arx.getColumn(k)); // compute fitness
  419.                 } catch (TooManyEvaluationsException e) {
  420.                     break generationLoop;
  421.                 }
  422.             }

  424.             // Compute fitnesses by adding value and penalty after scaling by value range.
  425.             double valueRange = valueRange(valuePenaltyPairs);
  426.             for (int iValue=0;iValue<valuePenaltyPairs.length;iValue++) {
  427.                  fitness[iValue] = valuePenaltyPairs[iValue].value + valuePenaltyPairs[iValue].penalty*valueRange;
  428.             }

  430.             // Sort by fitness and compute weighted mean into xmean
  431.             final int[] arindex = sortedIndices(fitness);
  432.             // Calculate new xmean, this is selection and recombination
  433.             final RealMatrix xold = xmean; // for speed up of Eq. (2) and (3)
  434.             final RealMatrix bestArx = selectColumns(arx, MathArrays.copyOf(arindex, mu));
  435.             xmean = bestArx.multiply(weights);
  436.             final RealMatrix bestArz = selectColumns(arz, MathArrays.copyOf(arindex, mu));
  437.             final RealMatrix zmean = bestArz.multiply(weights);
  438.             final boolean hsig = updateEvolutionPaths(zmean, xold);
  439.             if (diagonalOnly <= 0) {
  440.                 updateCovariance(hsig, bestArx, arz, arindex, xold);
  441.             } else {
  442.                 updateCovarianceDiagonalOnly(hsig, bestArz);
  443.             }
  444.             // Adapt step size sigma - Eq. (5)
  445.             sigma *= FastMath.exp(FastMath.min(1, (normps/chiN - 1) * cs / damps));
  446.             final double bestFitness = fitness[arindex[0]];
  447.             final double worstFitness = fitness[arindex[arindex.length - 1]];
  448.             if (bestValue > bestFitness) {
  449.                 bestValue = bestFitness;
  450.                 lastResult = optimum;
  451.                 optimum = new PointValuePair(fitfun.repair(bestArx.getColumn(0)),
  452.                                              isMinimize ? bestFitness : -bestFitness);
  453.                 if (getConvergenceChecker() != null && lastResult != null &&
  454.                     getConvergenceChecker().converged(iterations, optimum, lastResult)) {
  455.                     break generationLoop;
  456.                 }
  457.             }
  458.             // handle termination criteria
  459.             // Break, if fitness is good enough
  460.             if (stopFitness != 0 && bestFitness < (isMinimize ? stopFitness : -stopFitness)) {
  461.                 break generationLoop;
  462.             }
  463.             final double[] sqrtDiagC = sqrt(diagC).getColumn(0);
  464.             final double[] pcCol = pc.getColumn(0);
  465.             for (int i = 0; i < dimension; i++) {
  466.                 if (sigma * FastMath.max(FastMath.abs(pcCol[i]), sqrtDiagC[i]) > stopTolX) {
  467.                     break;
  468.                 }
  469.                 if (i >= dimension - 1) {
  470.                     break generationLoop;
  471.                 }
  472.             }
  473.             for (int i = 0; i < dimension; i++) {
  474.                 if (sigma * sqrtDiagC[i] > stopTolUpX) {
  475.                     break generationLoop;
  476.                 }
  477.             }
  478.             final double historyBest = min(fitnessHistory);
  479.             final double historyWorst = max(fitnessHistory);
  480.             if (iterations > 2 &&
  481.                 FastMath.max(historyWorst, worstFitness) -
  482.                 FastMath.min(historyBest, bestFitness) < stopTolFun) {
  483.                 break generationLoop;
  484.             }
  485.             if (iterations > fitnessHistory.length &&
  486.                 historyWorst - historyBest < stopTolHistFun) {
  487.                 break generationLoop;
  488.             }
  489.             // condition number of the covariance matrix exceeds 1e14
  490.             if (max(diagD) / min(diagD) > 1e7) {
  491.                 break generationLoop;
  492.             }
  493.             // user defined termination
  494.             if (getConvergenceChecker() != null) {
  495.                 final PointValuePair current
  496.                     = new PointValuePair(bestArx.getColumn(0),
  497.                                          isMinimize ? bestFitness : -bestFitness);
  498.                 if (lastResult != null &&
  499.                     getConvergenceChecker().converged(iterations, current, lastResult)) {
  500.                     break generationLoop;
  501.                     }
  502.                 lastResult = current;
  503.             }
  504.             // Adjust step size in case of equal function values (flat fitness)
  505.             if (bestValue == fitness[arindex[(int)(0.1+lambda/4.)]]) {
  506.                 sigma *= FastMath.exp(0.2 + cs / damps);
  507.             }
  508.             if (iterations > 2 && FastMath.max(historyWorst, bestFitness) -
  509.                 FastMath.min(historyBest, bestFitness) == 0) {
  510.                 sigma *= FastMath.exp(0.2 + cs / damps);
  511.             }
  512.             // store best in history
  513.             push(fitnessHistory,bestFitness);
  514.             if (generateStatistics) {
  515.                 statisticsSigmaHistory.add(sigma);
  516.                 statisticsFitnessHistory.add(bestFitness);
  517.                 statisticsMeanHistory.add(xmean.transpose());
  518.                 statisticsDHistory.add(diagD.transpose().scalarMultiply(1E5));
  519.             }
  520.         }
  521.         return optimum;
  522.     }

  524.     /**
  525.      * Scans the list of (required and optional) optimization data that
  526.      * characterize the problem.
  527.      *
  528.      * @param optData Optimization data. The following data will be looked for:
  529.      * <ul>
  530.      *  <li>{@link Sigma}</li>
  531.      *  <li>{@link PopulationSize}</li>
  532.      * </ul>
  533.      */
  534.     @Override
  535.     protected void parseOptimizationData(OptimizationData... optData) {
  536.         // Allow base class to register its own data.
  537.         super.parseOptimizationData(optData);

  539.         // The existing values (as set by the previous call) are reused if
  540.         // not provided in the argument list.
  541.         for (OptimizationData data : optData) {
  542.             if (data instanceof Sigma) {
  543.                 inputSigma = ((Sigma) data).getSigma();
  544.                 continue;
  545.             }
  546.             if (data instanceof PopulationSize) {
  547.                 lambda = ((PopulationSize) data).getPopulationSize();
  548.                 continue;
  549.             }
  550.         }

  552.         checkParameters();
  553.     }

  555.     /**
  556.      * Checks dimensions and values of boundaries and inputSigma if defined.
  557.      */
  558.     private void checkParameters() {
  559.         final double[] init = getStartPoint();
  560.         final double[] lB = getLowerBound();
  561.         final double[] uB = getUpperBound();

  563.         if (inputSigma != null) {
  564.             if (inputSigma.length != init.length) {
  565.                 throw new DimensionMismatchException(inputSigma.length, init.length);
  566.             }
  567.             for (int i = 0; i < init.length; i++) {
  568.                 if (inputSigma[i] > uB[i] - lB[i]) {
  569.                     throw new OutOfRangeException(inputSigma[i], 0, uB[i] - lB[i]);
  570.                 }
  571.             }
  572.         }
  573.     }

  575.     /**
  576.      * Initialization of the dynamic search parameters
  577.      *
  578.      * @param guess Initial guess for the arguments of the fitness function.
  579.      */
  580.     private void initializeCMA(double[] guess) {
  581.         if (lambda <= 0) {
  582.             throw new NotStrictlyPositiveException(lambda);
  583.         }
  584.         // initialize sigma
  585.         final double[][] sigmaArray = new double[guess.length][1];
  586.         for (int i = 0; i < guess.length; i++) {
  587.             sigmaArray[i][0] = inputSigma[i];
  588.         }
  589.         final RealMatrix insigma = new Array2DRowRealMatrix(sigmaArray, false);
  590.         sigma = max(insigma); // overall standard deviation

  592.         // initialize termination criteria
  593.         stopTolUpX = 1e3 * max(insigma);
  594.         stopTolX = 1e-11 * max(insigma);
  595.         stopTolFun = 1e-12;
  596.         stopTolHistFun = 1e-13;

  598.         // initialize selection strategy parameters
  599.         mu = lambda / 2; // number of parents/points for recombination
  600.         logMu2 = FastMath.log(mu + 0.5);
  601.         weights = log(sequence(1, mu, 1)).scalarMultiply(-1).scalarAdd(logMu2);
  602.         double sumw = 0;
  603.         double sumwq = 0;
  604.         for (int i = 0; i < mu; i++) {
  605.             double w = weights.getEntry(i, 0);
  606.             sumw += w;
  607.             sumwq += w * w;
  608.         }
  609.         weights = weights.scalarMultiply(1 / sumw);
  610.         mueff = sumw * sumw / sumwq; // variance-effectiveness of sum w_i x_i

  612.         // initialize dynamic strategy parameters and constants
  613.         cc = (4 + mueff / dimension) /
  614.                 (dimension + 4 + 2 * mueff / dimension);
  615.         cs = (mueff + 2) / (dimension + mueff + 3.);
  616.         damps = (1 + 2 * FastMath.max(0, FastMath.sqrt((mueff - 1) /
  617.                                                        (dimension + 1)) - 1)) *
  618.             FastMath.max(0.3,
  619.                          1 - dimension / (1e-6 + maxIterations)) + cs; // minor increment
  620.         ccov1 = 2 / ((dimension + 1.3) * (dimension + 1.3) + mueff);
  621.         ccovmu = FastMath.min(1 - ccov1, 2 * (mueff - 2 + 1 / mueff) /
  622.                               ((dimension + 2) * (dimension + 2) + mueff));
  623.         ccov1Sep = FastMath.min(1, ccov1 * (dimension + 1.5) / 3);
  624.         ccovmuSep = FastMath.min(1 - ccov1, ccovmu * (dimension + 1.5) / 3);
  625.         chiN = FastMath.sqrt(dimension) *
  626.                 (1 - 1 / ((double) 4 * dimension) + 1 / ((double) 21 * dimension * dimension));
  627.         // intialize CMA internal values - updated each generation
  628.         xmean = MatrixUtils.createColumnRealMatrix(guess); // objective variables
  629.         diagD = insigma.scalarMultiply(1 / sigma);
  630.         diagC = square(diagD);
  631.         pc = zeros(dimension, 1); // evolution paths for C and sigma
  632.         ps = zeros(dimension, 1); // B defines the coordinate system
  633.         normps = ps.getFrobeniusNorm();

  635.         B = eye(dimension, dimension);
  636.         D = ones(dimension, 1); // diagonal D defines the scaling
  637.         BD = times(B, repmat(diagD.transpose(), dimension, 1));
  638.         C = B.multiply(diag(square(D)).multiply(B.transpose())); // covariance
  639.         historySize = 10 + (int) (3 * 10 * dimension / (double) lambda);
  640.         fitnessHistory = new double[historySize]; // history of fitness values
  641.         for (int i = 0; i < historySize; i++) {
  642.             fitnessHistory[i] = Double.MAX_VALUE;
  643.         }
  644.     }

  646.     /**
  647.      * Update of the evolution paths ps and pc.
  648.      *
  649.      * @param zmean Weighted row matrix of the gaussian random numbers generating
  650.      * the current offspring.
  651.      * @param xold xmean matrix of the previous generation.
  652.      * @return hsig flag indicating a small correction.
  653.      */
  654.     private boolean updateEvolutionPaths(RealMatrix zmean, RealMatrix xold) {
  655.         ps = ps.scalarMultiply(1 - cs).add(
  656.                 B.multiply(zmean).scalarMultiply(
  657.                         FastMath.sqrt(cs * (2 - cs) * mueff)));
  658.         normps = ps.getFrobeniusNorm();
  659.         final boolean hsig = normps /
  660.             FastMath.sqrt(1 - FastMath.pow(1 - cs, 2 * iterations)) /
  661.             chiN < 1.4 + 2 / ((double) dimension + 1);
  662.         pc = pc.scalarMultiply(1 - cc);
  663.         if (hsig) {
  664.             pc = pc.add(xmean.subtract(xold).scalarMultiply(FastMath.sqrt(cc * (2 - cc) * mueff) / sigma));
  665.         }
  666.         return hsig;
  667.     }

  669.     /**
  670.      * Update of the covariance matrix C for diagonalOnly > 0
  671.      *
  672.      * @param hsig Flag indicating a small correction.
  673.      * @param bestArz Fitness-sorted matrix of the gaussian random values of the
  674.      * current offspring.
  675.      */
  676.     private void updateCovarianceDiagonalOnly(boolean hsig,
  677.                                               final RealMatrix bestArz) {
  678.         // minor correction if hsig==false
  679.         double oldFac = hsig ? 0 : ccov1Sep * cc * (2 - cc);
  680.         oldFac += 1 - ccov1Sep - ccovmuSep;
  681.         diagC = diagC.scalarMultiply(oldFac) // regard old matrix
  682.             .add(square(pc).scalarMultiply(ccov1Sep)) // plus rank one update
  683.             .add((times(diagC, square(bestArz).multiply(weights))) // plus rank mu update
  684.                  .scalarMultiply(ccovmuSep));
  685.         diagD = sqrt(diagC); // replaces eig(C)
  686.         if (diagonalOnly > 1 &&
  687.             iterations > diagonalOnly) {
  688.             // full covariance matrix from now on
  689.             diagonalOnly = 0;
  690.             B = eye(dimension, dimension);
  691.             BD = diag(diagD);
  692.             C = diag(diagC);
  693.         }
  694.     }

  696.     /**
  697.      * Update of the covariance matrix C.
  698.      *
  699.      * @param hsig Flag indicating a small correction.
  700.      * @param bestArx Fitness-sorted matrix of the argument vectors producing the
  701.      * current offspring.
  702.      * @param arz Unsorted matrix containing the gaussian random values of the
  703.      * current offspring.
  704.      * @param arindex Indices indicating the fitness-order of the current offspring.
  705.      * @param xold xmean matrix of the previous generation.
  706.      */
  707.     private void updateCovariance(boolean hsig, final RealMatrix bestArx,
  708.                                   final RealMatrix arz, final int[] arindex,
  709.                                   final RealMatrix xold) {
  710.         double negccov = 0;
  711.         if (ccov1 + ccovmu > 0) {
  712.             final RealMatrix arpos = bestArx.subtract(repmat(xold, 1, mu))
  713.                 .scalarMultiply(1 / sigma); // mu difference vectors
  714.             final RealMatrix roneu = pc.multiply(pc.transpose())
  715.                 .scalarMultiply(ccov1); // rank one update
  716.             // minor correction if hsig==false
  717.             double oldFac = hsig ? 0 : ccov1 * cc * (2 - cc);
  718.             oldFac += 1 - ccov1 - ccovmu;
  719.             if (isActiveCMA) {
  720.                 // Adapt covariance matrix C active CMA
  721.                 negccov = (1 - ccovmu) * 0.25 * mueff /
  722.                     (FastMath.pow(dimension + 2, 1.5) + 2 * mueff);
  723.                 // keep at least 0.66 in all directions, small popsize are most
  724.                 // critical
  725.                 final double negminresidualvariance = 0.66;
  726.                 // where to make up for the variance loss
  727.                 final double negalphaold = 0.5;
  728.                 // prepare vectors, compute negative updating matrix Cneg
  729.                 final int[] arReverseIndex = reverse(arindex);
  730.                 RealMatrix arzneg = selectColumns(arz, MathArrays.copyOf(arReverseIndex, mu));
  731.                 RealMatrix arnorms = sqrt(sumRows(square(arzneg)));
  732.                 final int[] idxnorms = sortedIndices(arnorms.getRow(0));
  733.                 final RealMatrix arnormsSorted = selectColumns(arnorms, idxnorms);
  734.                 final int[] idxReverse = reverse(idxnorms);
  735.                 final RealMatrix arnormsReverse = selectColumns(arnorms, idxReverse);
  736.                 arnorms = divide(arnormsReverse, arnormsSorted);
  737.                 final int[] idxInv = inverse(idxnorms);
  738.                 final RealMatrix arnormsInv = selectColumns(arnorms, idxInv);
  739.                 // check and set learning rate negccov
  740.                 final double negcovMax = (1 - negminresidualvariance) /
  741.                     square(arnormsInv).multiply(weights).getEntry(0, 0);
  742.                 if (negccov > negcovMax) {
  743.                     negccov = negcovMax;
  744.                 }
  745.                 arzneg = times(arzneg, repmat(arnormsInv, dimension, 1));
  746.                 final RealMatrix artmp = BD.multiply(arzneg);
  747.                 final RealMatrix Cneg = artmp.multiply(diag(weights)).multiply(artmp.transpose());
  748.                 oldFac += negalphaold * negccov;
  749.                 C = C.scalarMultiply(oldFac)
  750.                     .add(roneu) // regard old matrix
  751.                     .add(arpos.scalarMultiply( // plus rank one update
  752.                                               ccovmu + (1 - negalphaold) * negccov) // plus rank mu update
  753.                          .multiply(times(repmat(weights, 1, dimension),
  754.                                          arpos.transpose())))
  755.                     .subtract(Cneg.scalarMultiply(negccov));
  756.             } else {
  757.                 // Adapt covariance matrix C - nonactive
  758.                 C = C.scalarMultiply(oldFac) // regard old matrix
  759.                     .add(roneu) // plus rank one update
  760.                     .add(arpos.scalarMultiply(ccovmu) // plus rank mu update
  761.                          .multiply(times(repmat(weights, 1, dimension),
  762.                                          arpos.transpose())));
  763.             }
  764.         }
  765.         updateBD(negccov);
  766.     }

  768.     /**
  769.      * Update B and D from C.
  770.      *
  771.      * @param negccov Negative covariance factor.
  772.      */
  773.     private void updateBD(double negccov) {
  774.         if (ccov1 + ccovmu + negccov > 0 &&
  775.             (iterations % 1. / (ccov1 + ccovmu + negccov) / dimension / 10.) < 1) {
  776.             // to achieve O(N^2)
  777.             C = triu(C, 0).add(triu(C, 1).transpose());
  778.             // enforce symmetry to prevent complex numbers
  779.             final EigenDecomposition eig = new EigenDecomposition(C);
  780.             B = eig.getV(); // eigen decomposition, B==normalized eigenvectors
  781.             D = eig.getD();
  782.             diagD = diag(D);
  783.             if (min(diagD) <= 0) {
  784.                 for (int i = 0; i < dimension; i++) {
  785.                     if (diagD.getEntry(i, 0) < 0) {
  786.                         diagD.setEntry(i, 0, 0);
  787.                     }
  788.                 }
  789.                 final double tfac = max(diagD) / 1e14;
  790.                 C = C.add(eye(dimension, dimension).scalarMultiply(tfac));
  791.                 diagD = diagD.add(ones(dimension, 1).scalarMultiply(tfac));
  792.             }
  793.             if (max(diagD) > 1e14 * min(diagD)) {
  794.                 final double tfac = max(diagD) / 1e14 - min(diagD);
  795.                 C = C.add(eye(dimension, dimension).scalarMultiply(tfac));
  796.                 diagD = diagD.add(ones(dimension, 1).scalarMultiply(tfac));
  797.             }
  798.             diagC = diag(C);
  799.             diagD = sqrt(diagD); // D contains standard deviations now
  800.             BD = times(B, repmat(diagD.transpose(), dimension, 1)); // O(n^2)
  801.         }
  802.     }

  804.     /**
  805.      * Pushes the current best fitness value in a history queue.
  806.      *
  807.      * @param vals History queue.
  808.      * @param val Current best fitness value.
  809.      */
  810.     private static void push(double[] vals, double val) {
  811.         for (int i = vals.length-1; i > 0; i--) {
  812.             vals[i] = vals[i-1];
  813.         }
  814.         vals[0] = val;
  815.     }

  817.     /**
  818.      * Sorts fitness values.
  819.      *
  820.      * @param doubles Array of values to be sorted.
  821.      * @return a sorted array of indices pointing into doubles.
  822.      */
  823.     private int[] sortedIndices(final double[] doubles) {
  824.         final DoubleIndex[] dis = new DoubleIndex[doubles.length];
  825.         for (int i = 0; i < doubles.length; i++) {
  826.             dis[i] = new DoubleIndex(doubles[i], i);
  827.         }
  828.         Arrays.sort(dis);
  829.         final int[] indices = new int[doubles.length];
  830.         for (int i = 0; i < doubles.length; i++) {
  831.             indices[i] = dis[i].index;
  832.         }
  833.         return indices;
  834.     }
  835.    /**
  836.      * Get range of values.
  837.      *
  838.      * @param vpPairs Array of valuePenaltyPairs to get range from.
  839.      * @return a double equal to maximum value minus minimum value.
  840.      */
  841.     private double valueRange(final ValuePenaltyPair[] vpPairs) {
  842.         double max = Double.NEGATIVE_INFINITY;
  843.         double min = Double.MAX_VALUE;
  844.         for (ValuePenaltyPair vpPair:vpPairs) {
  845.             if (vpPair.value > max) {
  846.                 max = vpPair.value;
  847.             }
  848.             if (vpPair.value < min) {
  849.                 min = vpPair.value;
  850.             }
  851.         }
  852.         return max-min;
  853.     }

  855.     /**
  856.      * Used to sort fitness values. Sorting is always in lower value first
  857.      * order.
  858.      */
  859.     private static class DoubleIndex implements Comparable<DoubleIndex> {
  860.         /** Value to compare. */
  861.         private final double value;
  862.         /** Index into sorted array. */
  863.         private final int index;

  865.         /**
  866.          * @param value Value to compare.
  867.          * @param index Index into sorted array.
  868.          */
  869.         DoubleIndex(double value, int index) {
  870.             this.value = value;
  871.             this.index = index;
  872.         }

  874.         /** {@inheritDoc} */
  875.         public int compareTo(DoubleIndex o) {
  876.             return Double.compare(value, o.value);
  877.         }

  879.         /** {@inheritDoc} */
  880.         @Override
  881.         public boolean equals(Object other) {

  883.             if (this == other) {
  884.                 return true;
  885.             }

  887.             if (other instanceof DoubleIndex) {
  888.                 return Double.compare(value, ((DoubleIndex) other).value) == 0;
  889.             }

  891.             return false;
  892.         }

  894.         /** {@inheritDoc} */
  895.         @Override
  896.         public int hashCode() {
  897.             long bits = Double.doubleToLongBits(value);
  898.             return (int) ((1438542 ^ (bits >>> 32) ^ bits) & 0xffffffff);
  899.         }
  900.     }
  901.     /**
  902.      * Stores the value and penalty (for repair of out of bounds point).
  903.      */
  904.     private static class ValuePenaltyPair {
  905.         /** Objective function value. */
  906.         private double value;
  907.         /** Penalty value for repair of out out of bounds points. */
  908.         private double penalty;

  910.         /**
  911.          * @param value Function value.
  912.          * @param penalty Out-of-bounds penalty.
  913.         */
  914.         public ValuePenaltyPair(final double value, final double penalty) {
  915.             this.value   = value;
  916.             this.penalty = penalty;
  917.         }
  918.     }


  921.     /**
  922.      * Normalizes fitness values to the range [0,1]. Adds a penalty to the
  923.      * fitness value if out of range.
  924.      */
  925.     private class FitnessFunction {
  926.         /**
  927.          * Flag indicating whether the objective variables are forced into their
  928.          * bounds if defined
  929.          */
  930.         private final boolean isRepairMode;

  932.         /** Simple constructor.
  933.          */
  934.         public FitnessFunction() {
  935.             isRepairMode = true;
  936.         }

  938.         /**
  939.          * @param point Normalized objective variables.
  940.          * @return the objective value + penalty for violated bounds.
  941.          */
  942.         public ValuePenaltyPair value(final double[] point) {
  943.             double value;
  944.             double penalty=0.0;
  945.             if (isRepairMode) {
  946.                 double[] repaired = repair(point);
  947.                 value = CMAESOptimizer.this.computeObjectiveValue(repaired);
  948.                 penalty =  penalty(point, repaired);
  949.             } else {
  950.                 value = CMAESOptimizer.this.computeObjectiveValue(point);
  951.             }
  952.             value = isMinimize ? value : -value;
  953.             penalty = isMinimize ? penalty : -penalty;
  954.             return new ValuePenaltyPair(value,penalty);
  955.         }

  957.         /**
  958.          * @param x Normalized objective variables.
  959.          * @return {@code true} if in bounds.
  960.          */
  961.         public boolean isFeasible(final double[] x) {
  962.             final double[] lB = CMAESOptimizer.this.getLowerBound();
  963.             final double[] uB = CMAESOptimizer.this.getUpperBound();

  965.             for (int i = 0; i < x.length; i++) {
  966.                 if (x[i] < lB[i]) {
  967.                     return false;
  968.                 }
  969.                 if (x[i] > uB[i]) {
  970.                     return false;
  971.                 }
  972.             }
  973.             return true;
  974.         }

  976.         /**
  977.          * @param x Normalized objective variables.
  978.          * @return the repaired (i.e. all in bounds) objective variables.
  979.          */
  980.         private double[] repair(final double[] x) {
  981.             final double[] lB = CMAESOptimizer.this.getLowerBound();
  982.             final double[] uB = CMAESOptimizer.this.getUpperBound();

  984.             final double[] repaired = new double[x.length];
  985.             for (int i = 0; i < x.length; i++) {
  986.                 if (x[i] < lB[i]) {
  987.                     repaired[i] = lB[i];
  988.                 } else if (x[i] > uB[i]) {
  989.                     repaired[i] = uB[i];
  990.                 } else {
  991.                     repaired[i] = x[i];
  992.                 }
  993.             }
  994.             return repaired;
  995.         }

  997.         /**
  998.          * @param x Normalized objective variables.
  999.          * @param repaired Repaired objective variables.
  1000.          * @return Penalty value according to the violation of the bounds.
  1001.          */
  1002.         private double penalty(final double[] x, final double[] repaired) {
  1003.             double penalty = 0;
  1004.             for (int i = 0; i < x.length; i++) {
  1005.                 double diff = FastMath.abs(x[i] - repaired[i]);
  1006.                 penalty += diff;
  1007.             }
  1008.             return isMinimize ? penalty : -penalty;
  1009.         }
  1010.     }

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

  1014.     /**
  1015.      * @param m Input matrix
  1016.      * @return Matrix representing the element-wise logarithm of m.
  1017.      */
  1018.     private static RealMatrix log(final RealMatrix m) {
  1019.         final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
  1020.         for (int r = 0; r < m.getRowDimension(); r++) {
  1021.             for (int c = 0; c < m.getColumnDimension(); c++) {
  1022.                 d[r][c] = FastMath.log(m.getEntry(r, c));
  1023.             }
  1024.         }
  1025.         return new Array2DRowRealMatrix(d, false);
  1026.     }

  1028.     /**
  1029.      * @param m Input matrix.
  1030.      * @return Matrix representing the element-wise square root of m.
  1031.      */
  1032.     private static RealMatrix sqrt(final RealMatrix m) {
  1033.         final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
  1034.         for (int r = 0; r < m.getRowDimension(); r++) {
  1035.             for (int c = 0; c < m.getColumnDimension(); c++) {
  1036.                 d[r][c] = FastMath.sqrt(m.getEntry(r, c));
  1037.             }
  1038.         }
  1039.         return new Array2DRowRealMatrix(d, false);
  1040.     }

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

  1057.     /**
  1058.      * @param m Input matrix 1.
  1059.      * @param n Input matrix 2.
  1060.      * @return the matrix where the elements of m and n are element-wise multiplied.
  1061.      */
  1062.     private static RealMatrix times(final RealMatrix m, final RealMatrix n) {
  1063.         final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
  1064.         for (int r = 0; r < m.getRowDimension(); r++) {
  1065.             for (int c = 0; c < m.getColumnDimension(); c++) {
  1066.                 d[r][c] = m.getEntry(r, c) * n.getEntry(r, c);
  1067.             }
  1068.         }
  1069.         return new Array2DRowRealMatrix(d, false);
  1070.     }

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

  1087.     /**
  1088.      * @param m Input matrix.
  1089.      * @param cols Columns to select.
  1090.      * @return Matrix representing the selected columns.
  1091.      */
  1092.     private static RealMatrix selectColumns(final RealMatrix m, final int[] cols) {
  1093.         final double[][] d = new double[m.getRowDimension()][cols.length];
  1094.         for (int r = 0; r < m.getRowDimension(); r++) {
  1095.             for (int c = 0; c < cols.length; c++) {
  1096.                 d[r][c] = m.getEntry(r, cols[c]);
  1097.             }
  1098.         }
  1099.         return new Array2DRowRealMatrix(d, false);
  1100.     }

  1102.     /**
  1103.      * @param m Input matrix.
  1104.      * @param k Diagonal position.
  1105.      * @return Upper triangular part of matrix.
  1106.      */
  1107.     private static RealMatrix triu(final RealMatrix m, int k) {
  1108.         final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
  1109.         for (int r = 0; r < m.getRowDimension(); r++) {
  1110.             for (int c = 0; c < m.getColumnDimension(); c++) {
  1111.                 d[r][c] = r <= c - k ? m.getEntry(r, c) : 0;
  1112.             }
  1113.         }
  1114.         return new Array2DRowRealMatrix(d, false);
  1115.     }

  1117.     /**
  1118.      * @param m Input matrix.
  1119.      * @return Row matrix representing the sums of the rows.
  1120.      */
  1121.     private static RealMatrix sumRows(final RealMatrix m) {
  1122.         final double[][] d = new double[1][m.getColumnDimension()];
  1123.         for (int c = 0; c < m.getColumnDimension(); c++) {
  1124.             double sum = 0;
  1125.             for (int r = 0; r < m.getRowDimension(); r++) {
  1126.                 sum += m.getEntry(r, c);
  1127.             }
  1128.             d[0][c] = sum;
  1129.         }
  1130.         return new Array2DRowRealMatrix(d, false);
  1131.     }

  1133.     /**
  1134.      * @param m Input matrix.
  1135.      * @return the diagonal n-by-n matrix if m is a column matrix or the column
  1136.      * matrix representing the diagonal if m is a n-by-n matrix.
  1137.      */
  1138.     private static RealMatrix diag(final RealMatrix m) {
  1139.         if (m.getColumnDimension() == 1) {
  1140.             final double[][] d = new double[m.getRowDimension()][m.getRowDimension()];
  1141.             for (int i = 0; i < m.getRowDimension(); i++) {
  1142.                 d[i][i] = m.getEntry(i, 0);
  1143.             }
  1144.             return new Array2DRowRealMatrix(d, false);
  1145.         } else {
  1146.             final double[][] d = new double[m.getRowDimension()][1];
  1147.             for (int i = 0; i < m.getColumnDimension(); i++) {
  1148.                 d[i][0] = m.getEntry(i, i);
  1149.             }
  1150.             return new Array2DRowRealMatrix(d, false);
  1151.         }
  1152.     }

  1154.     /**
  1155.      * Copies a column from m1 to m2.
  1156.      *
  1157.      * @param m1 Source matrix.
  1158.      * @param col1 Source column.
  1159.      * @param m2 Target matrix.
  1160.      * @param col2 Target column.
  1161.      */
  1162.     private static void copyColumn(final RealMatrix m1, int col1,
  1163.                                    RealMatrix m2, int col2) {
  1164.         for (int i = 0; i < m1.getRowDimension(); i++) {
  1165.             m2.setEntry(i, col2, m1.getEntry(i, col1));
  1166.         }
  1167.     }

  1169.     /**
  1170.      * @param n Number of rows.
  1171.      * @param m Number of columns.
  1172.      * @return n-by-m matrix filled with 1.
  1173.      */
  1174.     private static RealMatrix ones(int n, int m) {
  1175.         final double[][] d = new double[n][m];
  1176.         for (int r = 0; r < n; r++) {
  1177.             Arrays.fill(d[r], 1);
  1178.         }
  1179.         return new Array2DRowRealMatrix(d, false);
  1180.     }

  1182.     /**
  1183.      * @param n Number of rows.
  1184.      * @param m Number of columns.
  1185.      * @return n-by-m matrix of 0 values out of diagonal, and 1 values on
  1186.      * the diagonal.
  1187.      */
  1188.     private static RealMatrix eye(int n, int m) {
  1189.         final double[][] d = new double[n][m];
  1190.         for (int r = 0; r < n; r++) {
  1191.             if (r < m) {
  1192.                 d[r][r] = 1;
  1193.             }
  1194.         }
  1195.         return new Array2DRowRealMatrix(d, false);
  1196.     }

  1198.     /**
  1199.      * @param n Number of rows.
  1200.      * @param m Number of columns.
  1201.      * @return n-by-m matrix of zero values.
  1202.      */
  1203.     private static RealMatrix zeros(int n, int m) {
  1204.         return new Array2DRowRealMatrix(n, m);
  1205.     }

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

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

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

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

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

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

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

  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.     }

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

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