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