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