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