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