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 }