You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by er...@apache.org on 2011/02/17 14:10:00 UTC
svn commit: r1071600 - in /commons/proper/math/trunk/src:
main/java/org/apache/commons/math/optimization/direct/CMAESOptimizer.java
site/xdoc/changes.xml
test/java/org/apache/commons/math/optimization/direct/CMAESOptimizerTest.java
Author: erans
Date: Thu Feb 17 13:10:00 2011
New Revision: 1071600
URL: http://svn.apache.org/viewvc?rev=1071600&view=rev
Log:
MATH-442
Implementation of the CMA-ES optimization algorithm provided by Dietmar Wolz
and Nikolaus Hansen.
Added:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/direct/CMAESOptimizer.java (with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/direct/CMAESOptimizerTest.java (with props)
Modified:
commons/proper/math/trunk/src/site/xdoc/changes.xml
Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/direct/CMAESOptimizer.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/direct/CMAESOptimizer.java?rev=1071600&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/direct/CMAESOptimizer.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/direct/CMAESOptimizer.java Thu Feb 17 13:10:00 2011
@@ -0,0 +1,1317 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math.optimization.direct;
+
+import java.util.ArrayList;
+import java.util.Arrays;
+import java.util.List;
+
+import org.apache.commons.math.analysis.MultivariateRealFunction;
+import org.apache.commons.math.exception.MultiDimensionMismatchException;
+import org.apache.commons.math.exception.NoDataException;
+import org.apache.commons.math.exception.NotPositiveException;
+import org.apache.commons.math.exception.OutOfRangeException;
+import org.apache.commons.math.exception.TooManyEvaluationsException;
+import org.apache.commons.math.linear.Array2DRowRealMatrix;
+import org.apache.commons.math.linear.EigenDecomposition;
+import org.apache.commons.math.linear.EigenDecompositionImpl;
+import org.apache.commons.math.linear.MatrixUtils;
+import org.apache.commons.math.linear.RealMatrix;
+import org.apache.commons.math.optimization.GoalType;
+import org.apache.commons.math.optimization.MultivariateRealOptimizer;
+import org.apache.commons.math.optimization.RealPointValuePair;
+import org.apache.commons.math.random.MersenneTwister;
+import org.apache.commons.math.random.RandomGenerator;
+import org.apache.commons.math.util.MathUtils;
+
+/**
+ * CMA-ES algorithm. This code is translated and adapted from the Matlab version
+ * of this algorithm as implemented in module {@code cmaes.m} version 3.51.
+ *
+ * Implements the active Covariance Matrix Adaptation Evolution Strategy (CMA-ES)
+ * for non-linear, non-convex, non-smooth, global function minimization.
+ * The CMA-Evolution Strategy (CMA-ES) is a reliable stochastic optimization method
+ * which should be applied, if derivative based methods, e.g. quasi-Newton BFGS or
+ * conjugate gradient, fail due to a rugged search landscape (e.g. noise, local
+ * optima, outlier, etc.) of the objective function. Like a
+ * quasi-Newton method the CMA-ES learns and applies a variable metric
+ * of the underlying search space. Unlike a quasi-Newton method the
+ * CMA-ES does neither estimate nor use gradients, making it considerably more
+ * reliable in terms of finding a good, or even close to optimal, solution, finally.
+ *
+ * <p>In general, on smooth objective functions the CMA-ES is roughly ten times
+ * slower than BFGS (counting objective function evaluations, no gradients provided).
+ * For up to <math>N=10</math> variables also the derivative-free simplex
+ * direct search method (Nelder and Mead) can be faster, but it is
+ * far less reliable than CMA-ES.</p>
+ *
+ * <p>The CMA-ES is particularly well suited for non-separable
+ * and/or badly conditioned problems.
+ * To observe the advantage of CMA compared to a conventional
+ * evolution strategy, it will usually take about <math>30 N</math> function
+ * evaluations. On difficult problems the complete
+ * optimization (a single run) is expected to take <em>roughly</em> between
+ * <math>30 N</math> and <math>300 N<sup>2</sup></math>
+ * function evaluations.</p>
+ *
+ * For more information, please refer to the following links:
+ * <ul>
+ * <li><a href="http://www.lri.fr/~hansen/cmaes.m">Matlab code</a></li>
+ * <li><a href="http://www.lri.fr/~hansen/cmaesintro.html">Introduction to CMA-ES</a></li>
+ * <li><a href="http://en.wikipedia.org/wiki/CMA-ES">Wikipedia</a></li>
+ * </ul>
+ *
+ * @version $Revision$ $Date$
+ * @since 3.0
+ */
+
+public class CMAESOptimizer extends
+ BaseAbstractScalarOptimizer<MultivariateRealFunction> implements
+ MultivariateRealOptimizer {
+
+ /** Default value for {@link #checkFeasableCount}: {@value}. */
+ public static final int DEFAULT_CHECKFEASABLECOUNT = 0;
+ /** Default value for {@link #stopfitness}: {@value}. */
+ public static final double DEFAULT_STOPFITNESS = 0;
+ /** Default value for {@link #isActiveCMA}: {@value}. */
+ public static final boolean DEFAULT_ISACTIVECMA = true;
+ /** Default value for {@link #maxIterations}: {@value}. */
+ public static final int DEFAULT_MAXITERATIONS = 30000;
+ /** Default value for {@link #diagonalOnly}: {@value}. */
+ public static final int DEFAULT_DIAGONALONLY = 0;
+ /** Default value for {@link #random}. */
+ public static final RandomGenerator DEFAULT_RANDOMGENERATOR = new MersenneTwister();
+
+ // global search parameters
+ /**
+ * Population size, offspring number. The primary strategy parameter to play
+ * with, which can be increased from its default value. Increasing the
+ * population size improves global search properties in exchange to speed.
+ * Speed decreases, as a rule, at most linearely with increasing population
+ * size. It is advisable to begin with the default small population size.
+ */
+ private int lambda; // population size
+ /**
+ * Covariance update mechanism, default is active CMA. isActiveCMA = true
+ * turns on "active CMA" with a negative update of the covariance matrix and
+ * checks for positive definiteness. OPTS.CMA.active = 2 does not check for
+ * pos. def. and is numerically faster. Active CMA usually speeds up the
+ * adaptation.
+ */
+ private boolean isActiveCMA;
+ /**
+ * Determines how often a new random offspring is generated in case it is
+ * not feasible / beyond the defined limits, default is 0. Only relevant if
+ * boundaries != null.
+ */
+ private int checkFeasableCount;
+ /**
+ * Lower and upper boundaries of the objective variables. boundaries == null
+ * means no boundaries.
+ */
+ private double[][] boundaries;
+ /**
+ * Individual sigma values - initial search volume. inputSigma determines
+ * the initial coordinate wise standard deviations for the search. Setting
+ * SIGMA one third of the initial search region is appropriate.
+ */
+ private double[] inputSigma;
+ /** Number of objective variables/problem dimension */
+ private int dimension;
+ /**
+ * Defines the number of initial iterations, where the covariance matrix
+ * remains diagonal and the algorithm has internally linear time complexity.
+ * diagonalOnly = 1 means keeping the covariance matrix always diagonal and
+ * this setting also exhibits linear space complexity. This can be
+ * particularly useful for dimension > 100.
+ * @see <a href="http://hal.archives-ouvertes.fr/inria-00287367/en">A Simple Modification in CMA-ES</a> .
+ */
+ private int diagonalOnly = 0;
+ /** Number of objective variables/problem dimension */
+ private boolean isMinimize = true;
+ /** Indicates whether statistic data is collected. */
+ private boolean generateStatistics = false;
+
+ // termination criteria
+ /** Maximal number of iterations allowed. */
+ private int maxIterations;
+ /** Limit for fitness value. */
+ private double stopfitness;
+ /** Stop if x-changes larger stopTolUpX. */
+ private double stopTolUpX;
+ /** Stop if x-change smaller stopTolX. */
+ private double stopTolX;
+ /** Stop if fun-changes smaller stopTolFun. */
+ private double stopTolFun;
+ /** Stop if back fun-changes smaller stopTolHistFun. */
+ private double stopTolHistFun;
+
+ // selection strategy parameters
+ /** Number of parents/points for recombination. */
+ private int mu; //
+ /** log(mu + 0.5), stored for efficiency. */
+ private double logMu2;
+ /** Array for weighted recombination. */
+ private RealMatrix weights;
+ /** Variance-effectiveness of sum w_i x_i. */
+ private double mueff; //
+
+ // dynamic strategy parameters and constants
+ /** Overall standard deviation - search volume. */
+ private double sigma;
+ /** Cumulation constant. */
+ private double cc;
+ /** Cumulation constant for step-size. */
+ private double cs;
+ /** Damping for step-size. */
+ private double damps;
+ /** Learning rate for rank-one update. */
+ private double ccov1;
+ /** Learning rate for rank-mu update' */
+ private double ccovmu;
+ /** Expectation of ||N(0,I)|| == norm(randn(N,1)). */
+ private double chiN;
+ /** Learning rate for rank-one update - diagonalOnly */
+ private double ccov1Sep;
+ /** Learning rate for rank-mu update - diagonalOnly */
+ private double ccovmuSep;
+
+ // CMA internal values - updated each generation
+ /** Objective variables. */
+ private RealMatrix xmean;
+ /** Evolution path. */
+ private RealMatrix pc;
+ /** Evolution path for sigma. */
+ private RealMatrix ps;
+ /** Norm of ps, stored for efficiency. */
+ private double normps;
+ /** Coordinate system. */
+ private RealMatrix B;
+ /** Scaling. */
+ private RealMatrix D;
+ /** B*D, stored for efficiency. */
+ private RealMatrix BD;
+ /** Diagonal of sqrt(D), stored for efficiency. */
+ private RealMatrix diagD;
+ /** Covariance matrix. */
+ private RealMatrix C;
+ /** Diagonal of C, used for diagonalOnly. */
+ private RealMatrix diagC;
+ /** Number of iterations already performed. */
+ private int iterations;
+
+ /** History queue of best values. */
+ private double[] fitnessHistory;
+ /** Size of history queue of best values. */
+ private int historySize;
+
+ /** Random generator. */
+ private RandomGenerator random;
+
+ /** History of sigma values. */
+ private List<Double> statisticsSigmaHistory = new ArrayList<Double>();
+ /** History of mean matrix. */
+ private List<RealMatrix> statisticsMeanHistory = new ArrayList<RealMatrix>();
+ /** History of fitness values. */
+ private List<Double> statisticsFitnessHistory = new ArrayList<Double>();
+ /** History of D matrix. */
+ private List<RealMatrix> statisticsDHistory = new ArrayList<RealMatrix>();
+
+ /**
+ * Default constructor, uses default parameters
+ */
+ public CMAESOptimizer() {
+ this(0);
+ }
+
+ /**
+ * @param lambda
+ * Population size.
+ */
+ public CMAESOptimizer(int lambda) {
+ this(lambda, null, null, DEFAULT_MAXITERATIONS, DEFAULT_STOPFITNESS,
+ DEFAULT_ISACTIVECMA, DEFAULT_DIAGONALONLY,
+ DEFAULT_CHECKFEASABLECOUNT, DEFAULT_RANDOMGENERATOR, false);
+ }
+
+ /**
+ * @param lambda
+ * Population size.
+ * @param inputSigma
+ * Initial search volume - sigma of offspring objective
+ * variables.
+ * @param boundaries
+ * Boundaries for objective variables.
+ */
+ public CMAESOptimizer(int lambda, double[] inputSigma,
+ double[][] boundaries) {
+ this(lambda, inputSigma, boundaries, DEFAULT_MAXITERATIONS, DEFAULT_STOPFITNESS,
+ DEFAULT_ISACTIVECMA, DEFAULT_DIAGONALONLY,
+ DEFAULT_CHECKFEASABLECOUNT, DEFAULT_RANDOMGENERATOR, false);
+ }
+
+ /**
+ * @param lambda
+ * Population size.
+ * @param inputSigma
+ * Initial search volume - sigma of offspring objective
+ * variables.
+ * @param boundaries
+ * Boundaries for objective variables.
+ * @param maxIterations
+ * Maximal number of iterations.
+ * @param stopfitness
+ * stop if objective function value < stopfitness.
+ * @param isActiveCMA
+ * Chooses the covariance matrix update method.
+ * @param diagonalOnly
+ * Number of initial iterations, where the covariance matrix
+ * remains diagonal.
+ * @param checkFeasableCount
+ * Determines how often new. random objective variables are
+ * generated in case they are out of bounds.
+ * @param random
+ * Used random generator.
+ * @param generateStatistics
+ * Indicates whether statistic data is collected.
+ */
+ public CMAESOptimizer(int lambda, double[] inputSigma,
+ double[][] boundaries, int maxIterations, double stopfitness,
+ boolean isActiveCMA, int diagonalOnly, int checkFeasableCount,
+ RandomGenerator random, boolean generateStatistics) {
+ this.lambda = lambda;
+ this.inputSigma = inputSigma;
+ this.boundaries = boundaries;
+ this.maxIterations = maxIterations;
+ this.stopfitness = stopfitness;
+ this.isActiveCMA = isActiveCMA;
+ this.diagonalOnly = diagonalOnly;
+ this.checkFeasableCount = checkFeasableCount;
+ this.random = random;
+ this.generateStatistics = generateStatistics;
+ }
+
+ /**
+ * @return History of sigma values.
+ */
+ public List<Double> getStatisticsSigmaHistory() {
+ return statisticsSigmaHistory;
+ }
+
+ /**
+ * @return History of mean matrix.
+ */
+ public List<RealMatrix> getStatisticsMeanHistory() {
+ return statisticsMeanHistory;
+ }
+
+ /**
+ * @return History of fitness values.
+ */
+ public List<Double> getStatisticsFitnessHistory() {
+ return statisticsFitnessHistory;
+ }
+
+ /**
+ * @return History of D matrix.
+ */
+ public List<RealMatrix> getStatisticsDHistory() {
+ return statisticsDHistory;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ protected RealPointValuePair doOptimize() {
+ checkParameters();
+ // -------------------- Initialization --------------------------------
+ isMinimize = getGoalType().equals(GoalType.MINIMIZE);
+ final FitnessFunction fitfun = new FitnessFunction(boundaries,
+ isMinimize);
+ final double[] guess = fitfun.encode(getStartPoint());
+ // number of objective variables/problem dimension
+ dimension = guess.length;
+ initializeCMA(guess);
+ iterations = 0;
+ double bestValue = fitfun.value(guess);
+ push(fitnessHistory, bestValue);
+ RealPointValuePair optimum = new RealPointValuePair(getStartPoint(),
+ isMinimize ? bestValue : -bestValue);
+ RealPointValuePair lastResult = null;
+
+ // -------------------- Generation Loop --------------------------------
+
+ generationLoop:
+ for (iterations = 1; iterations <= maxIterations; iterations++) {
+ // Generate and evaluate lambda offspring
+ RealMatrix arz = randn1(dimension, lambda);
+ RealMatrix arx = zeros(dimension, lambda);
+ double[] fitness = new double[lambda];
+ // generate random offspring
+ for (int k = 0; k < lambda; k++) {
+ RealMatrix arxk = null;
+ for (int i = 0; i < checkFeasableCount+1; i++) {
+ if (diagonalOnly <= 0)
+ arxk = xmean.add(BD.multiply(arz.getColumnMatrix(k))
+ .scalarMultiply(sigma)); // m + sig * Normal(0,C)
+ else
+ arxk = xmean.add(times(diagD,arz.getColumnMatrix(k))
+ .scalarMultiply(sigma));
+ if (i >= checkFeasableCount || fitfun.isFeasible(arxk.getColumn(0)))
+ break;
+ // regenerate random arguments for row
+ arz.setColumn(k, randn(dimension));
+ }
+ copyColumn(arxk, 0, arx, k);
+ try {
+ fitness[k] = fitfun.value(arx.getColumn(k)); // compute fitness
+ } catch (TooManyEvaluationsException e) {
+ break generationLoop;
+ }
+ }
+ // Sort by fitness and compute weighted mean into xmean
+ int[] arindex = sortedIndices(fitness);
+ // Calculate new xmean, this is selection and recombination
+ RealMatrix xold = xmean; // for speed up of Eq. (2) and (3)
+ RealMatrix bestArx = selectColumns(arx, MathUtils.copyOf(arindex, mu));
+ xmean = bestArx.multiply(weights);
+ RealMatrix bestArz = selectColumns(arz, MathUtils.copyOf(arindex, mu));
+ RealMatrix zmean = bestArz.multiply(weights);
+ boolean hsig = updateEvolutionPaths(zmean, xold);
+ if (diagonalOnly <= 0)
+ updateCovariance(hsig, bestArx, arz, arindex, xold);
+ else
+ updateCovarianceDiagonalOnly(hsig, bestArz, xold);
+ // Adapt step size sigma - Eq. (5)
+ sigma *= Math.exp(Math.min(1.0,(normps/chiN - 1.)*cs/damps));
+ double bestFitness = fitness[arindex[0]];
+ double worstFitness = fitness[arindex[arindex.length-1]];
+ if (bestValue > bestFitness) {
+ bestValue = bestFitness;
+ lastResult = optimum;
+ optimum = new RealPointValuePair(
+ fitfun.decode(bestArx.getColumn(0)),
+ isMinimize ? bestFitness : -bestFitness);
+ if (getConvergenceChecker() != null && lastResult != null) {
+ if (getConvergenceChecker().converged(
+ iterations, optimum, lastResult))
+ break generationLoop;
+ }
+ }
+ // handle termination criteria
+ // Break, if fitness is good enough
+ if (stopfitness != 0) { // only if stopfitness is defined
+ if (bestFitness < (isMinimize ? stopfitness : -stopfitness))
+ break generationLoop;
+ }
+ double[] sqrtDiagC = sqrt(diagC).getColumn(0);
+ double[] pcCol = pc.getColumn(0);
+ for (int i = 0; i < dimension; i++) {
+ if (sigma*(Math.max(Math.abs(pcCol[i]), sqrtDiagC[i])) > stopTolX)
+ break;
+ if (i >= dimension-1)
+ break generationLoop;
+ }
+ for (int i = 0; i < dimension; i++)
+ if (sigma*sqrtDiagC[i] > stopTolUpX)
+ break generationLoop;
+ double historyBest = min(fitnessHistory);
+ double historyWorst = max(fitnessHistory);
+ if (iterations > 2 && Math.max(historyWorst, worstFitness) -
+ Math.min(historyBest, bestFitness) < stopTolFun)
+ break generationLoop;
+ if (iterations > fitnessHistory.length &&
+ historyWorst-historyBest < stopTolHistFun)
+ break generationLoop;
+ // condition number of the covariance matrix exceeds 1e14
+ if (max(diagD)/min(diagD) > 1e7)
+ break generationLoop;
+ // user defined termination
+ if (getConvergenceChecker() != null) {
+ RealPointValuePair current =
+ new RealPointValuePair(bestArx.getColumn(0),
+ isMinimize ? bestFitness : -bestFitness);
+ if (lastResult != null &&
+ getConvergenceChecker().converged(
+ iterations, current, lastResult))
+ break generationLoop;
+ lastResult = current;
+ }
+ // Adjust step size in case of equal function values (flat fitness)
+ if (bestValue == fitness[arindex[(int)(0.1+lambda/4.)]])
+ sigma = sigma * Math.exp(0.2+cs/damps);
+ if (iterations > 2 && Math.max(historyWorst, bestFitness) -
+ Math.min(historyBest, bestFitness) == 0)
+ sigma = sigma * Math.exp(0.2+cs/damps);
+ // store best in history
+ push(fitnessHistory,bestFitness);
+ fitfun.setValueRange(worstFitness-bestFitness);
+ if (generateStatistics) {
+ statisticsSigmaHistory.add(sigma);
+ statisticsFitnessHistory.add(bestFitness);
+ statisticsMeanHistory.add(xmean.transpose());
+ statisticsDHistory.add(diagD.transpose().scalarMultiply(1E5));
+ }
+ }
+ return optimum;
+ }
+
+ /**
+ * Checks dimensions and values of boundaries and inputSigma if defined.
+ */
+ private void checkParameters() {
+ double[] init = getStartPoint();
+ if (boundaries != null) {
+ if (boundaries.length != 2)
+ throw new MultiDimensionMismatchException(
+ new Integer[] { boundaries.length },
+ new Integer[] { 2 });
+ if (boundaries[0] == null || boundaries[1] == null)
+ throw new NoDataException();
+ if (boundaries[0].length != init.length)
+ throw new MultiDimensionMismatchException(
+ new Integer[] { boundaries[0].length },
+ new Integer[] { init.length });
+ if (boundaries[1].length != init.length)
+ throw new MultiDimensionMismatchException(
+ new Integer[] { boundaries[1].length },
+ new Integer[] { init.length });
+ for (int i = 0; i < init.length; i++) {
+ if (boundaries[0][i] > init[i] || boundaries[1][i] < init[i])
+ throw new OutOfRangeException(init[i], boundaries[0][i],
+ boundaries[1][i]);
+ }
+ }
+ if (inputSigma != null) {
+ if (inputSigma.length != init.length)
+ throw new MultiDimensionMismatchException(
+ new Integer[] { inputSigma.length },
+ new Integer[] { init.length });
+ for (int i = 0; i < init.length; i++) {
+ if (inputSigma[i] < 0)
+ throw new NotPositiveException(inputSigma[i]);
+ if (boundaries != null) {
+ if (inputSigma[i] > 1.0)
+ throw new OutOfRangeException(inputSigma[i], 0, 1.0);
+ }
+ }
+ }
+ }
+
+ /**
+ * Initialization of the dynamic search parameters
+ *
+ * @param guess
+ * initial guess for the arguments of the fitness function
+ */
+
+ private void initializeCMA(double[] guess) {
+ if (lambda <= 0)
+ lambda = 4 + (int) (3. * Math.log(dimension));
+ // initialize sigma
+ double[][] sigmaArray = new double[guess.length][1];
+ for (int i = 0; i < guess.length; i++)
+ sigmaArray[i][0] = inputSigma != null ? inputSigma[i] : 0.3;
+ RealMatrix insigma = new Array2DRowRealMatrix(sigmaArray, false);
+ sigma = max(insigma); // overall standard deviation
+
+ // initialize termination criteria
+ stopTolUpX = 1e3 * max(insigma);
+ stopTolX = 1e-11 * max(insigma);
+ stopTolFun = 1e-12;
+ stopTolHistFun = 1e-13;
+
+ // initialize selection strategy parameters
+ mu = lambda / 2; // number of parents/points for recombination
+ logMu2 = Math.log(mu + 0.5);
+ weights = log(sequence(1, mu, 1)).scalarMultiply(-1.).scalarAdd(logMu2);
+ double sumw = 0;
+ double sumwq = 0;
+ for (int i = 0; i < mu; i++) {
+ double w = weights.getEntry(i, 0);
+ sumw += w;
+ sumwq += w * w;
+ }
+ weights = weights.scalarMultiply(1. / sumw);
+ mueff = sumw * sumw / sumwq; // variance-effectiveness of sum w_i x_i
+
+ // initialize dynamic strategy parameters and constants
+ cc = (4. + mueff / dimension) /
+ (dimension + 4. + 2. * mueff / dimension);
+ cs = (mueff + 2.) / (dimension + mueff + 3.);
+ damps = (1. + 2. * Math.max(0, Math.sqrt((mueff - 1.) /
+ (dimension + 1.)) - 1.)) *
+ Math.max(0.3, 1. - dimension /
+ (1e-6 + Math.min(maxIterations, getMaxEvaluations() /
+ lambda))) + cs; // minor increment
+ ccov1 = 2. / ((dimension + 1.3) * (dimension + 1.3) + mueff);
+ ccovmu = Math.min(1 - ccov1, 2. * (mueff - 2. + 1. / mueff) /
+ ((dimension + 2.) * (dimension + 2.) + mueff));
+ ccov1Sep = Math.min(1, ccov1 * (dimension + 1.5) / 3.);
+ ccovmuSep = Math.min(1 - ccov1, ccovmu * (dimension + 1.5) / 3.);
+ chiN = Math.sqrt(dimension) *
+ (1. - 1. / (4. * dimension) + 1 / (21. * dimension * dimension));
+ // intialize CMA internal values - updated each generation
+ xmean = MatrixUtils.createColumnRealMatrix(guess); // objective
+ // variables
+ diagD = insigma.scalarMultiply(1. / sigma);
+ diagC = square(diagD);
+ pc = zeros(dimension, 1); // evolution paths for C and sigma
+ ps = zeros(dimension, 1); // B defines the coordinate system
+ normps = norm(ps);
+
+ B = eye(dimension, dimension);
+ D = ones(dimension, 1); // diagonal D defines the scaling
+ BD = times(B, repmat(diagD.transpose(), dimension, 1));
+ C = B.multiply(diag(square(D)).multiply(B.transpose())); // covariance
+ historySize = 10 + (int) (3. * 10. * dimension / lambda);
+ fitnessHistory = new double[historySize]; // history of fitness values
+ for (int i = 0; i < historySize; i++)
+ fitnessHistory[i] = Double.MAX_VALUE;
+ }
+
+ /**
+ * Update of the evolution paths ps and pc
+ *
+ * @param zmean
+ * weighted row matrix of the gaussian random numbers generating
+ * the current offspring
+ * @param xold
+ * xmean matrix of the previous generation
+ * @return hsig flag indicating a small correction
+ */
+ private boolean updateEvolutionPaths(RealMatrix zmean, RealMatrix xold) {
+ ps = ps.scalarMultiply(1. - cs).add(
+ B.multiply(zmean).scalarMultiply(
+ Math.sqrt(cs * (2. - cs) * mueff)));
+ normps = norm(ps);
+ boolean hsig = normps /
+ Math.sqrt(1. - Math.pow(1. - cs, 2. * iterations)) /
+ chiN < 1.4 + 2. / (dimension + 1.);
+ pc = pc.scalarMultiply(1. - cc);
+ if (hsig)
+ pc = pc.add(xmean.subtract(xold).scalarMultiply(
+ Math.sqrt(cc * (2. - cc) * mueff) / sigma));
+ return hsig;
+ }
+
+ /**
+ * Update of the covariance matrix C for diagonalOnly > 0
+ *
+ * @param hsig
+ * flag indicating a small correction
+ * @param bestArz
+ * fitness-sorted matrix of the gaussian random values of the
+ * current offspring
+ * @param xold
+ * xmean matrix of the previous generation
+ */
+ private void updateCovarianceDiagonalOnly(boolean hsig, final RealMatrix bestArz,
+ final RealMatrix xold) {
+ // minor correction if hsig==false
+ double oldFac = hsig ? 0 : ccov1Sep * cc * (2. - cc);
+ oldFac += 1. - ccov1Sep - ccovmuSep;
+ diagC = diagC.scalarMultiply(oldFac) // regard old matrix
+ // plus rank one update
+ .add(square(pc).scalarMultiply(ccov1Sep))
+ // plus rank mu update
+ .add((times(diagC, square(bestArz).multiply(weights)))
+ .scalarMultiply(ccovmuSep));
+ diagD = sqrt(diagC); // replaces eig(C)
+ if (diagonalOnly > 1 && iterations > diagonalOnly) {
+ // full covariance matrix from now on
+ diagonalOnly = 0;
+ B = eye(dimension, dimension);
+ BD = diag(diagD);
+ C = diag(diagC);
+ }
+ }
+
+ /**
+ * Update of the covariance matrix C
+ *
+ * @param hsig
+ * flag indicating a small correction
+ * @param bestArx
+ * fitness-sorted matrix of the argument vectors producing the
+ * current offspring
+ * @param arz
+ * unsorted matrix containing the gaussian random values of the
+ * current offspring
+ * @param arindex
+ * indices indicating the fitness-order of the current offspring
+ * @param xold
+ * xmean matrix of the previous generation
+ */
+ private void updateCovariance(boolean hsig, final RealMatrix bestArx,
+ final RealMatrix arz, final int[] arindex, final RealMatrix xold) {
+ double negccov = 0;
+ if (ccov1 + ccovmu > 0) {
+ RealMatrix arpos = bestArx.subtract(repmat(xold, 1, mu))
+ .scalarMultiply(1. / sigma); // mu difference vectors
+ RealMatrix roneu = pc.multiply(pc.transpose())
+ .scalarMultiply(ccov1); // rank one update
+ // minor correction if hsig==false
+ double oldFac = hsig ? 0 : ccov1 * cc * (2. - cc);
+ oldFac += 1. - ccov1 - ccovmu;
+ if (isActiveCMA) {
+ // Adapt covariance matrix C active CMA
+ negccov = (1. - ccovmu) * 0.25 * mueff /
+ (Math.pow(dimension + 2., 1.5) + 2. * mueff);
+ double negminresidualvariance = 0.66;
+ // keep at least 0.66 in all directions, small popsize are most
+ // critical
+ double negalphaold = 0.5; // where to make up for the variance
+ // loss,
+ // prepare vectors, compute negative updating matrix Cneg
+ int[] arReverseIndex = reverse(arindex);
+ RealMatrix arzneg
+ = selectColumns(arz, MathUtils.copyOf(arReverseIndex, mu));
+ RealMatrix arnorms = sqrt(sumRows(square(arzneg)));
+ int[] idxnorms = sortedIndices(arnorms.getRow(0));
+ RealMatrix arnormsSorted = selectColumns(arnorms, idxnorms);
+ int[] idxReverse = reverse(idxnorms);
+ RealMatrix arnormsReverse = selectColumns(arnorms, idxReverse);
+ arnorms = divide(arnormsReverse, arnormsSorted);
+ int[] idxInv = inverse(idxnorms);
+ RealMatrix arnormsInv = selectColumns(arnorms, idxInv);
+ // check and set learning rate negccov
+ double negcovMax = (1. - negminresidualvariance) /
+ square(arnormsInv).multiply(weights).getEntry(0, 0);
+ if (negccov > negcovMax)
+ negccov = negcovMax;
+ arzneg = times(arzneg, repmat(arnormsInv, dimension, 1));
+ RealMatrix artmp = BD.multiply(arzneg);
+ RealMatrix Cneg = artmp.multiply(diag(weights)).multiply(
+ artmp.transpose());
+ oldFac += negalphaold * negccov;
+ C = C.scalarMultiply(oldFac)
+ // regard old matrix
+ .add(roneu)
+ // plus rank one update
+ .add(arpos.scalarMultiply(
+ // plus rank mu update
+ ccovmu + (1. - negalphaold) * negccov)
+ .multiply(
+ times(repmat(weights, 1, dimension),
+ arpos.transpose())))
+ .subtract(Cneg.scalarMultiply(negccov));
+ } else {
+ // Adapt covariance matrix C - nonactive
+ C = C.scalarMultiply(oldFac) // regard old matrix
+ .add(roneu)
+ // plus rank one update
+ .add(arpos.scalarMultiply(ccovmu) // plus rank mu update
+ .multiply(
+ times(repmat(weights, 1, dimension),
+ arpos.transpose())));
+ }
+ }
+ updateBD(negccov);
+ }
+
+ /**
+ * Update B and D from C
+ *
+ * @param negccov
+ * Negative covariance factor.
+ */
+ private void updateBD(double negccov) {
+ if (ccov1 + ccovmu + negccov > 0 &&
+ (iterations % 1. / (ccov1 + ccovmu + negccov) / dimension / 10.) < 1.) {
+ // to achieve O(N^2)
+ C = triu(C, 0).add(triu(C, 1).transpose());
+ // enforce symmetry to prevent complex numbers
+ EigenDecomposition eig = new EigenDecompositionImpl(C, 1.0);
+ B = eig.getV(); // eigen decomposition, B==normalized eigenvectors
+ D = eig.getD();
+ diagD = diag(D);
+ if (min(diagD) <= 0) {
+ for (int i = 0; i < dimension; i++)
+ if (diagD.getEntry(i, 0) < 0)
+ diagD.setEntry(i, 0, 0.);
+ double tfac = max(diagD) / 1e14;
+ C = C.add(eye(dimension, dimension).scalarMultiply(tfac));
+ diagD = diagD.add(ones(dimension, 1).scalarMultiply(tfac));
+ }
+ if (max(diagD) > 1e14 * min(diagD)) {
+ double tfac = max(diagD) / 1e14 - min(diagD);
+ C = C.add(eye(dimension, dimension).scalarMultiply(tfac));
+ diagD = diagD.add(ones(dimension, 1).scalarMultiply(tfac));
+ }
+ diagC = diag(C);
+ diagD = sqrt(diagD); // D contains standard deviations now
+ BD = times(B, repmat(diagD.transpose(), dimension, 1)); // O(n^2)
+ }
+ }
+
+ /**
+ * Pushes the current best fitness value in a history queue.
+ *
+ * @param vals
+ * the history queue
+ * @param val
+ * current best fitness value
+ */
+ private static void push(double[] vals, double val) {
+ for (int i = vals.length-1; i > 0; i--)
+ vals[i] = vals[i-1];
+ vals[0] = val;
+ }
+
+ /**
+ * Sorts fitness values.
+ *
+ * @param doubles
+ * array of values to be sorted
+ * @return sorted array of indices pointing into doubles
+ */
+ private int[] sortedIndices(final double[] doubles) {
+ DoubleIndex[] dis = new DoubleIndex[doubles.length];
+ for (int i = 0; i < doubles.length; i++)
+ dis[i] = new DoubleIndex(doubles[i], i);
+ Arrays.sort(dis);
+ int[] indices = new int[doubles.length];
+ for (int i = 0; i < doubles.length; i++)
+ indices[i] = dis[i].index;
+ return indices;
+ }
+
+ /**
+ * Used to sort fitness values. Sorting is always in lower value first
+ * order.
+ */
+ private static class DoubleIndex implements Comparable<DoubleIndex> {
+
+ /** Value to compare. */
+ private double value;
+ /** Index into sorted array. */
+ private int index;
+
+ /**
+ * @param value
+ * Value to compare.
+ * @param index
+ * Index into sorted array.
+ */
+ DoubleIndex(double value, int index) {
+ this.value = value;
+ this.index = index;
+ }
+
+ /** {@inheritDoc} */
+ public int compareTo(DoubleIndex o) {
+ return Double.compare(value, o.value);
+ }
+ }
+
+ /**
+ * Normalizes fitness values to the range [0,1]. Adds a penalty to the
+ * fitness value if out of range. The penalty is adjusted by calling
+ * setValueRange().
+ */
+ private class FitnessFunction {
+
+ /** Optional bounds for the objective variables */
+ private double[][] boundaries;
+ /** Determines the penalty for boundary violations */
+ private double valueRange = 1.0;
+ /**
+ * Flag indicating whether the objective variables are forced into their
+ * bounds if defined
+ */
+ private boolean isRepairMode = true;
+ /** Flag indicating the optimization goal. */
+ private boolean isMinimize = true;
+
+ /**
+ * @param boundaries
+ * Bounds for the objective variables.
+ * @param isMinimize
+ * Flag indicating the optimization goal.
+ */
+ private FitnessFunction(final double[][] boundaries,
+ final boolean isMinimize) {
+ this.boundaries = boundaries;
+ this.isMinimize = isMinimize;
+ }
+
+ /**
+ * @param x
+ * Original objective variables.
+ * @return Normalized objective variables.
+ */
+ private double[] encode(final double[] x) {
+ if (boundaries == null)
+ return x;
+ double[] res = new double[x.length];
+ for (int i = 0; i < x.length; i++) {
+ double diff = boundaries[1][i] - boundaries[0][i];
+ res[i] = (x[i] - boundaries[0][i]) / diff;
+ }
+ return res;
+ }
+
+ /**
+ * @param x
+ * Normalized objective variables.
+ * @return Original objective variables.
+ */
+ private double[] decode(final double[] x) {
+ if (boundaries == null)
+ return x;
+ double[] res = new double[x.length];
+ for (int i = 0; i < x.length; i++) {
+ double diff = boundaries[1][i] - boundaries[0][i];
+ res[i] = diff * x[i] + boundaries[0][i];
+ }
+ return res;
+ }
+
+ /**
+ * @param point
+ * Normalized objective variables.
+ * @return Objective value + penalty for violated bounds.
+ */
+ private double value(final double[] point) {
+ double value;
+ if (boundaries != null && isRepairMode) {
+ double[] repaired = repair(point);
+ value = CMAESOptimizer.this
+ .computeObjectiveValue(decode(repaired)) +
+ penalty(point, repaired);
+ } else
+ value = CMAESOptimizer.this
+ .computeObjectiveValue(decode(point));
+ return isMinimize ? value : -value;
+ }
+
+ /**
+ * @param x
+ * Normalized objective variables.
+ * @return True if in bounds
+ */
+ private boolean isFeasible(final double[] x) {
+ if (boundaries == null)
+ return true;
+ for (int i = 0; i < x.length; i++) {
+ if (x[i] < 0)
+ return false;
+ if (x[i] > 1.0)
+ return false;
+ }
+ return true;
+ }
+
+ /**
+ * @param valueRange
+ * Adjusts the penalty computation.
+ */
+ private void setValueRange(double valueRange) {
+ this.valueRange = valueRange;
+ }
+
+ /**
+ * @param x
+ * Normalized objective variables.
+ * @return Repaired objective variables - all in bounds.
+ */
+ private double[] repair(final double[] x) {
+ double[] repaired = new double[x.length];
+ for (int i = 0; i < x.length; i++) {
+ if (x[i] < 0)
+ repaired[i] = 0;
+ else if (x[i] > 1.0)
+ repaired[i] = 1.0;
+ else
+ repaired[i] = x[i];
+ }
+ return repaired;
+ }
+
+ /**
+ * @param x
+ * Normalized objective variables.
+ * @param repaired
+ * Repaired objective variables.
+ * @return Penalty value according to the violation of the bounds.
+ */
+ private double penalty(final double[] x, final double[] repaired) {
+ double penalty = 0;
+ for (int i = 0; i < x.length; i++) {
+ double diff = Math.abs(x[i] - repaired[i]);
+ penalty += diff * valueRange;
+ }
+ return isMinimize ? penalty : -penalty;
+ }
+ }
+
+ // -----Matrix utility functions similar to the Matlab build in functions------
+
+ /**
+ * @param m
+ * Input matrix
+ * @return Matrix representing the element wise logarithm of m.
+ */
+ private static RealMatrix log(final RealMatrix m) {
+ double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
+ for (int r = 0; r < m.getRowDimension(); r++)
+ for (int c = 0; c < m.getColumnDimension(); c++)
+ d[r][c] = Math.log(m.getEntry(r, c));
+ return new Array2DRowRealMatrix(d, false);
+ }
+
+ /**
+ * @param m
+ * Input matrix
+ * @return Matrix representing the element wise square root of m.
+ */
+ private static RealMatrix sqrt(final RealMatrix m) {
+ double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
+ for (int r = 0; r < m.getRowDimension(); r++)
+ for (int c = 0; c < m.getColumnDimension(); c++)
+ d[r][c] = Math.sqrt(m.getEntry(r, c));
+ return new Array2DRowRealMatrix(d, false);
+ }
+
+ /**
+ * @param m Input matrix
+ * @return Matrix representing the element wise square (^2) of m.
+ */
+ private static RealMatrix square(final RealMatrix m) {
+ double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
+ for (int r = 0; r < m.getRowDimension(); r++)
+ for (int c = 0; c < m.getColumnDimension(); c++) {
+ double e = m.getEntry(r, c);
+ d[r][c] = e * e;
+ }
+ return new Array2DRowRealMatrix(d, false);
+ }
+
+ /**
+ * @param m
+ * Input matrix 1.
+ * @param n
+ * Input matrix 2.
+ * @return Matrix where the elements of m and m are element wise multiplied.
+ */
+ private static RealMatrix times(final RealMatrix m, final RealMatrix n) {
+ double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
+ for (int r = 0; r < m.getRowDimension(); r++)
+ for (int c = 0; c < m.getColumnDimension(); c++)
+ d[r][c] = m.getEntry(r, c)*n.getEntry(r, c);
+ return new Array2DRowRealMatrix(d, false);
+ }
+
+ /**
+ * @param m
+ * Input matrix 1.
+ * @param n
+ * Input matrix 2.
+ * @return Matrix where the elements of m and m are element wise divided.
+ */
+ private static RealMatrix divide(final RealMatrix m, final RealMatrix n) {
+ double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
+ for (int r = 0; r < m.getRowDimension(); r++)
+ for (int c = 0; c < m.getColumnDimension(); c++)
+ d[r][c] = m.getEntry(r, c)/n.getEntry(r, c);
+ return new Array2DRowRealMatrix(d, false);
+ }
+
+ /**
+ * @param m Input matrix.
+ * @param cols Columns to select.
+ * @return Matrix representing the selected columns.
+ */
+ private static RealMatrix selectColumns(final RealMatrix m, final int[] cols) {
+ double[][] d = new double[m.getRowDimension()][cols.length];
+ for (int r = 0; r < m.getRowDimension(); r++)
+ for (int c = 0; c < cols.length; c++)
+ d[r][c] = m.getEntry(r, cols[c]);
+ return new Array2DRowRealMatrix(d, false);
+ }
+
+ /**
+ * @param m Input matrix.
+ * @param k diagonal position.
+ * @return Upper triangular part of matrix.
+ */
+ private static RealMatrix triu(final RealMatrix m, int k) {
+ double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
+ for (int r = 0; r < m.getRowDimension(); r++)
+ for (int c = 0; c < m.getColumnDimension(); c++)
+ d[r][c] = r <= c - k ? m.getEntry(r, c) : 0;
+ return new Array2DRowRealMatrix(d, false);
+ }
+
+ /**
+ * @param m
+ * Input matrix.
+ * @return Norm of the matrix.
+ */
+ private static double norm(final RealMatrix m) {
+ double sum = 0;
+ for (int r = 0; r < m.getRowDimension(); r++)
+ for (int c = 0; c < m.getColumnDimension(); c++) {
+ double e = m.getEntry(r, c);
+ sum += e*e;
+ }
+ return Math.sqrt(sum);
+ }
+
+ /**
+ * @param m
+ * Input matrix.
+ * @return Row matrix representing the sums of the rows.
+ */
+ private static RealMatrix sumRows(final RealMatrix m) {
+ double[][] d = new double[1][m.getColumnDimension()];
+ for (int c = 0; c < m.getColumnDimension(); c++) {
+ double sum = 0;
+ for (int r = 0; r < m.getRowDimension(); r++)
+ sum += m.getEntry(r, c);
+ d[0][c] = sum;
+ }
+ return new Array2DRowRealMatrix(d, false);
+ }
+
+ /**
+ * @param m
+ * Input matrix.
+ * @return Diagonal n X n matrix if m is a column matrix, Rolumn matrix
+ * representing the diagonal if m is a nXn matrix.
+ */
+ private static RealMatrix diag(final RealMatrix m) {
+ if (m.getColumnDimension() == 1) {
+ double[][] d = new double[m.getRowDimension()][m.getRowDimension()];
+ for (int i = 0; i < m.getRowDimension(); i++)
+ d[i][i] = m.getEntry(i, 0);
+ return new Array2DRowRealMatrix(d, false);
+ } else {
+ double[][] d = new double[m.getRowDimension()][1];
+ for (int i = 0; i < m.getColumnDimension(); i++)
+ d[i][0] = m.getEntry(i, i);
+ return new Array2DRowRealMatrix(d, false);
+ }
+ }
+
+ /**
+ * Copies a row from m1 to m2.
+ *
+ * @param m1
+ * Source matrix 1.
+ * @param col1
+ * Source column.
+ * @param m2
+ * Target matrix.
+ * @param col2
+ * Target column.
+ */
+ private static void copyColumn(final RealMatrix m1, int col1, RealMatrix m2, int col2) {
+ for (int i = 0; i < m1.getRowDimension(); i++)
+ m2.setEntry(i, col2, m1.getEntry(i, col1));
+ }
+
+ /**
+ * @param n
+ * Number of rows.
+ * @param m
+ * Number of columns.
+ * @return n X m matrix of 1.0-values.
+ */
+ private static RealMatrix ones(int n, int m) {
+ double[][] d = new double[n][m];
+ for (int r = 0; r < n; r++)
+ Arrays.fill(d[r], 1.0);
+ return new Array2DRowRealMatrix(d, false);
+ }
+
+ /**
+ * @param n
+ * Number of rows.
+ * @param m
+ * Number of columns.
+ * @return n X m matrix of 0.0-values, diagonal has values 1.0.
+ */
+ private static RealMatrix eye(int n, int m) {
+ double[][] d = new double[n][m];
+ for (int r = 0; r < n; r++)
+ if (r < m)
+ d[r][r] = 1;
+ return new Array2DRowRealMatrix(d, false);
+ }
+
+ /**
+ * @param n
+ * Number of rows.
+ * @param m
+ * Number of columns.
+ * @return n X m matrix of 0.0-values.
+ */
+ private static RealMatrix zeros(int n, int m) {
+ return new Array2DRowRealMatrix(n, m);
+ }
+
+ /**
+ * @param mat
+ * Input matrix.
+ * @param n
+ * Number of row replicates.
+ * @param m
+ * Number of column replicates.
+ * @return Matrix which replicates the input matrix in both directions.
+ */
+ private static RealMatrix repmat(final RealMatrix mat, int n, int m) {
+ int rd = mat.getRowDimension();
+ int cd = mat.getColumnDimension();
+ double[][] d = new double[n * rd][m * cd];
+ for (int r = 0; r < n * rd; r++)
+ for (int c = 0; c < m * cd; c++)
+ d[r][c] = mat.getEntry(r % rd, c % cd);
+ return new Array2DRowRealMatrix(d, false);
+ }
+
+ /**
+ * @param start
+ * Start value.
+ * @param end
+ * End value.
+ * @param step
+ * Step size.
+ * @return Sequence as column matrix.
+ */
+ private static RealMatrix sequence(double start, double end, double step) {
+ int size = (int) ((end - start) / step + 1);
+ double[][] d = new double[size][1];
+ double value = start;
+ for (int r = 0; r < size; r++) {
+ d[r][0] = value;
+ value += step;
+ }
+ return new Array2DRowRealMatrix(d, false);
+ }
+
+ /**
+ * @param m
+ * Input matrix.
+ * @return Maximum of matrix element values.
+ */
+ private static double max(final RealMatrix m) {
+ double max = -Double.MAX_VALUE;
+ for (int r = 0; r < m.getRowDimension(); r++)
+ for (int c = 0; c < m.getColumnDimension(); c++) {
+ double e = m.getEntry(r, c);
+ if (max < e)
+ max = e;
+ }
+ return max;
+ }
+
+ /**
+ * @param m
+ * Input matrix.
+ * @return Minimum of matrix element values.
+ */
+ private static double min(final RealMatrix m) {
+ double min = Double.MAX_VALUE;
+ for (int r = 0; r < m.getRowDimension(); r++)
+ for (int c = 0; c < m.getColumnDimension(); c++) {
+ double e = m.getEntry(r, c);
+ if (min > e)
+ min = e;
+ }
+ return min;
+ }
+
+ /**
+ * @param m
+ * Input array.
+ * @return Maximum of array values.
+ */
+ private static double max(final double[] m) {
+ double max = -Double.MAX_VALUE;
+ for (int r = 0; r < m.length; r++)
+ if (max < m[r])
+ max = m[r];
+ return max;
+ }
+
+ /**
+ * @param m
+ * Input array.
+ * @return Minimum of array values.
+ */
+ private static double min(final double[] m) {
+ double min = Double.MAX_VALUE;
+ for (int r = 0; r < m.length; r++)
+ if (min > m[r])
+ min = m[r];
+ return min;
+ }
+
+ /**
+ * @param indices
+ * Input index array.
+ * @return Inverse of the mapping defined by indices
+ */
+ private static int[] inverse(final int[] indices) {
+ int[] inverse = new int[indices.length];
+ for (int i = 0; i < indices.length; i++)
+ inverse[indices[i]] = i;
+ return inverse;
+ }
+
+ /**
+ * @param indices
+ * Input index array.
+ * @return Indices in inverse order (last is first)
+ */
+ private static int[] reverse(final int[] indices) {
+ int[] reverse = new int[indices.length];
+ for (int i = 0; i < indices.length; i++)
+ reverse[i] = indices[indices.length - i - 1];
+ return reverse;
+ }
+
+ /**
+ * @param size
+ * Length of random array.
+ * @return Array of gaussian random numbers.
+ */
+ private double[] randn(int size) {
+ double[] randn = new double[size];
+ for (int i = 0; i < size; i++)
+ randn[i] = random.nextGaussian();
+ return randn;
+ }
+
+ /**
+ * @param size
+ * Number of rows.
+ * @param popSize
+ * Population size.
+ * @return 2-dimensional matrix of gaussian random numbers.
+ */
+ private RealMatrix randn1(int size, int popSize) {
+ double[][] d = new double[size][popSize];
+ for (int r = 0; r < size; r++)
+ for (int c = 0; c < popSize; c++)
+ d[r][c] = random.nextGaussian();
+ return new Array2DRowRealMatrix(d, false);
+ }
+}
+
Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/direct/CMAESOptimizer.java
------------------------------------------------------------------------------
svn:eol-style = native
Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=1071600&r1=1071599&r2=1071600&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Thu Feb 17 13:10:00 2011
@@ -52,6 +52,9 @@ The <action> type attribute can be add,u
If the output is not quite correct, check for invisible trailing spaces!
-->
<release version="3.0" date="TBD" description="TBD">
+ <action dev="erans" type="fix" issue="MATH-442" due-to="Dietmar Wolz">
+ Implementation of the CMA-ES optimization algorithm.
+ </action>
<action dev="erans" type="fix" issue="MATH-513">
The interface "ParametricRealFunction" (in package "optimization.fitting") has
been renamed to "ParametricUnivariateRealFunction" and moved to package "analysis".
Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/direct/CMAESOptimizerTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/direct/CMAESOptimizerTest.java?rev=1071600&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/direct/CMAESOptimizerTest.java (added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/direct/CMAESOptimizerTest.java Thu Feb 17 13:10:00 2011
@@ -0,0 +1,667 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math.optimization.direct;
+
+import java.util.Arrays;
+import java.util.Random;
+
+import org.apache.commons.math.MathException;
+import org.apache.commons.math.analysis.MultivariateRealFunction;
+import org.apache.commons.math.exception.MathUserException;
+import org.apache.commons.math.exception.MultiDimensionMismatchException;
+import org.apache.commons.math.exception.NoDataException;
+import org.apache.commons.math.exception.NotPositiveException;
+import org.apache.commons.math.exception.OutOfRangeException;
+import org.apache.commons.math.optimization.GoalType;
+import org.apache.commons.math.optimization.MultivariateRealOptimizer;
+import org.apache.commons.math.optimization.RealPointValuePair;
+import org.apache.commons.math.random.MersenneTwister;
+import org.junit.Assert;
+import org.junit.Test;
+
+/**
+ * Test for {@link CMAESOptimizer}.
+ */
+public class CMAESOptimizerTest {
+
+ static final int DIM = 13;
+ static final int LAMBDA = 4 + (int)(3.*Math.log(DIM));
+
+ @Test(expected = OutOfRangeException.class)
+ public void testInitOutofbounds() throws MathUserException, MathException {
+ double[] startPoint = point(DIM,3);
+ double[] insigma = null;
+ double[][] boundaries = boundaries(DIM,-1,2);
+ RealPointValuePair expected =
+ new RealPointValuePair(point(DIM,1.0),0.0);
+ doTest(new Rosen(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+ 1e-13, 1e-6, 100000, expected);
+ }
+
+ @Test(expected = MultiDimensionMismatchException.class)
+ public void testBoundariesDimensionMismatch() throws MathUserException, MathException {
+ double[] startPoint = point(DIM,0.5);
+ double[] insigma = null;
+ double[][] boundaries = boundaries(DIM+1,-1,2);
+ RealPointValuePair expected =
+ new RealPointValuePair(point(DIM,1.0),0.0);
+ doTest(new Rosen(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+ 1e-13, 1e-6, 100000, expected);
+ }
+
+ @Test(expected = NoDataException.class)
+ public void testBoundariesNoData() throws MathUserException, MathException {
+ double[] startPoint = point(DIM,0.5);
+ double[] insigma = null;
+ double[][] boundaries = boundaries(DIM,-1,2);
+ boundaries[1] = null;
+ RealPointValuePair expected =
+ new RealPointValuePair(point(DIM,1.0),0.0);
+ doTest(new Rosen(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+ 1e-13, 1e-6, 100000, expected);
+ }
+
+ @Test(expected = NotPositiveException.class)
+ public void testInputSigmaNegative() throws MathUserException, MathException {
+ double[] startPoint = point(DIM,0.5);
+ double[] insigma = point(DIM,-0.5);
+ double[][] boundaries = null;
+ RealPointValuePair expected =
+ new RealPointValuePair(point(DIM,1.0),0.0);
+ doTest(new Rosen(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+ 1e-13, 1e-6, 100000, expected);
+ }
+
+ @Test(expected = OutOfRangeException.class)
+ public void testInputSigmaOutOfRange() throws MathUserException, MathException {
+ double[] startPoint = point(DIM,0.5);
+ double[] insigma = point(DIM, 1.1);
+ double[][] boundaries = boundaries(DIM,-1,2);
+ RealPointValuePair expected =
+ new RealPointValuePair(point(DIM,1.0),0.0);
+ doTest(new Rosen(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+ 1e-13, 1e-6, 100000, expected);
+ }
+
+ @Test(expected = MultiDimensionMismatchException.class)
+ public void testInputSigmaDimensionMismatch() throws MathUserException, MathException {
+ double[] startPoint = point(DIM,0.5);
+ double[] insigma = point(DIM+1,-0.5);
+ double[][] boundaries = null;;
+ RealPointValuePair expected =
+ new RealPointValuePair(point(DIM,1.0),0.0);
+ doTest(new Rosen(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+ 1e-13, 1e-6, 100000, expected);
+ }
+
+ @Test
+ public void testRosen() throws MathException {
+ double[] startPoint = point(DIM,0.1);
+ double[] insigma = point(DIM,0.1);
+ double[][] boundaries = null;
+ RealPointValuePair expected =
+ new RealPointValuePair(point(DIM,1.0),0.0);
+ doTest(new Rosen(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+ 1e-13, 1e-6, 100000, expected);
+ doTest(new Rosen(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13,
+ 1e-13, 1e-6, 100000, expected);
+ }
+
+ @Test
+ public void testMaximize() throws MathException {
+ double[] startPoint = point(DIM,1.0);
+ double[] insigma = point(DIM,0.1);
+ double[][] boundaries = null;
+ RealPointValuePair expected =
+ new RealPointValuePair(point(DIM,0.0),1.0);
+ doTest(new MinusElli(), startPoint, insigma, boundaries,
+ GoalType.MAXIMIZE, LAMBDA, true, 0, 1.0-1e-13,
+ 2e-10, 5e-6, 100000, expected);
+ doTest(new MinusElli(), startPoint, insigma, boundaries,
+ GoalType.MAXIMIZE, LAMBDA, false, 0, 1.0-1e-13,
+ 2e-10, 5e-6, 100000, expected);
+ boundaries = boundaries(DIM,-0.3,0.3);
+ startPoint = point(DIM,0.1);
+ doTest(new MinusElli(), startPoint, insigma, boundaries,
+ GoalType.MAXIMIZE, LAMBDA, true, 0, 1.0-1e-13,
+ 2e-10, 5e-6, 100000, expected);
+ }
+
+ @Test
+ public void testEllipse() throws MathException {
+ double[] startPoint = point(DIM,1.0);
+ double[] insigma = point(DIM,0.1);
+ double[][] boundaries = null;
+ RealPointValuePair expected =
+ new RealPointValuePair(point(DIM,0.0),0.0);
+ doTest(new Elli(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+ 1e-13, 1e-6, 100000, expected);
+ doTest(new Elli(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13,
+ 1e-13, 1e-6, 100000, expected);
+ }
+
+ @Test
+ public void testElliRotated() throws MathException {
+ double[] startPoint = point(DIM,1.0);
+ double[] insigma = point(DIM,0.1);
+ double[][] boundaries = null;
+ RealPointValuePair expected =
+ new RealPointValuePair(point(DIM,0.0),0.0);
+ doTest(new ElliRotated(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+ 1e-13, 1e-6, 100000, expected);
+ doTest(new ElliRotated(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13,
+ 1e-13, 1e-6, 100000, expected);
+ }
+
+ @Test
+ public void testCigar() throws MathException {
+ double[] startPoint = point(DIM,1.0);
+ double[] insigma = point(DIM,0.1);
+ double[][] boundaries = null;
+ RealPointValuePair expected =
+ new RealPointValuePair(point(DIM,0.0),0.0);
+ doTest(new Cigar(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+ 1e-13, 1e-6, 200000, expected);
+ doTest(new Cigar(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13,
+ 1e-13, 1e-6, 100000, expected);
+ }
+
+ @Test
+ public void testTwoAxes() throws MathException {
+ double[] startPoint = point(DIM,1.0);
+ double[] insigma = point(DIM,0.1);
+ double[][] boundaries = null;
+ RealPointValuePair expected =
+ new RealPointValuePair(point(DIM,0.0),0.0);
+ doTest(new TwoAxes(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, 2*LAMBDA, true, 0, 1e-13,
+ 1e-13, 1e-6, 200000, expected);
+ doTest(new TwoAxes(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, 2*LAMBDA, false, 0, 1e-13,
+ 1e-8, 1e-3, 200000, expected);
+ }
+
+ @Test
+ public void testCigTab() throws MathException {
+ double[] startPoint = point(DIM,1.0);
+ double[] insigma = point(DIM,0.3);
+ double[][] boundaries = null;
+ RealPointValuePair expected =
+ new RealPointValuePair(point(DIM,0.0),0.0);
+ doTest(new CigTab(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+ 1e-13, 5e-5, 100000, expected);
+ doTest(new CigTab(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13,
+ 1e-13, 5e-5, 100000, expected);
+ }
+
+ @Test
+ public void testSphere() throws MathException {
+ double[] startPoint = point(DIM,1.0);
+ double[] insigma = point(DIM,0.1);
+ double[][] boundaries = null;
+ RealPointValuePair expected =
+ new RealPointValuePair(point(DIM,0.0),0.0);
+ doTest(new Sphere(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+ 1e-13, 1e-6, 100000, expected);
+ doTest(new Sphere(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13,
+ 1e-13, 1e-6, 100000, expected);
+ }
+
+ @Test
+ public void testTablet() throws MathException {
+ double[] startPoint = point(DIM,1.0);
+ double[] insigma = point(DIM,0.1);
+ double[][] boundaries = null;
+ RealPointValuePair expected =
+ new RealPointValuePair(point(DIM,0.0),0.0);
+ doTest(new Tablet(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+ 1e-13, 1e-6, 100000, expected);
+ doTest(new Tablet(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13,
+ 1e-13, 1e-6, 100000, expected);
+ }
+
+ @Test
+ public void testDiffPow() throws MathException {
+ double[] startPoint = point(DIM,1.0);
+ double[] insigma = point(DIM,0.1);
+ double[][] boundaries = null;
+ RealPointValuePair expected =
+ new RealPointValuePair(point(DIM,0.0),0.0);
+ doTest(new DiffPow(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, 10, true, 0, 1e-13,
+ 1e-8, 1e-1, 100000, expected);
+ doTest(new DiffPow(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, 10, false, 0, 1e-13,
+ 1e-8, 2e-1, 100000, expected);
+ }
+
+ @Test
+ public void testSsDiffPow() throws MathException {
+ double[] startPoint = point(DIM,1.0);
+ double[] insigma = point(DIM,0.1);
+ double[][] boundaries = null;
+ RealPointValuePair expected =
+ new RealPointValuePair(point(DIM,0.0),0.0);
+ doTest(new SsDiffPow(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, 10, true, 0, 1e-13,
+ 1e-4, 1e-1, 200000, expected);
+ doTest(new SsDiffPow(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, 10, false, 0, 1e-13,
+ 1e-4, 1e-1, 200000, expected);
+ }
+
+ @Test
+ public void testAckley() throws MathException {
+ double[] startPoint = point(DIM,1.0);
+ double[] insigma = point(DIM,1.0);
+ double[][] boundaries = null;
+ RealPointValuePair expected =
+ new RealPointValuePair(point(DIM,0.0),0.0);
+ doTest(new Ackley(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, 2*LAMBDA, true, 0, 1e-13,
+ 1e-9, 1e-5, 100000, expected);
+ doTest(new Ackley(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, 2*LAMBDA, false, 0, 1e-13,
+ 1e-9, 1e-5, 100000, expected);
+ }
+
+ @Test
+ public void testRastrigin() throws MathException {
+ double[] startPoint = point(DIM,0.1);
+ double[] insigma = point(DIM,0.1);
+ double[][] boundaries = null;
+ RealPointValuePair expected =
+ new RealPointValuePair(point(DIM,0.0),0.0);
+ doTest(new Rastrigin(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, (int)(200*Math.sqrt(DIM)), true, 0, 1e-13,
+ 1e-13, 1e-6, 200000, expected);
+ doTest(new Rastrigin(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, (int)(200*Math.sqrt(DIM)), false, 0, 1e-13,
+ 1e-13, 1e-6, 200000, expected);
+ }
+
+ @Test
+ public void testConstrainedRosen() throws MathException {
+ double[] startPoint = point(DIM,0.1);
+ double[] insigma = point(DIM,0.1);
+ double[][] boundaries = boundaries(DIM,-1,2);
+ RealPointValuePair expected =
+ new RealPointValuePair(point(DIM,1.0),0.0);
+ doTest(new Rosen(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, 2*LAMBDA, true, 0, 1e-13,
+ 1e-13, 1e-6, 100000, expected);
+ doTest(new Rosen(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, 2*LAMBDA, false, 0, 1e-13,
+ 1e-13, 1e-6, 100000, expected);
+ }
+
+ @Test
+ public void testDiagonalRosen() throws MathException {
+ double[] startPoint = point(DIM,0.1);
+ double[] insigma = point(DIM,0.1);
+ double[][] boundaries = null;
+ RealPointValuePair expected =
+ new RealPointValuePair(point(DIM,1.0),0.0);
+ doTest(new Rosen(), startPoint, insigma, boundaries,
+ GoalType.MINIMIZE, LAMBDA, false, 1, 1e-13,
+ 1e-10, 1e-4, 1000000, expected);
+ }
+
+ /**
+ * @param func Function to optimize.
+ * @param startPoint Starting point.
+ * @param inSigma Individual input sigma.
+ * @param boundaries Upper / lower point limit.
+ * @param goal Minimization or maximization.
+ * @param lambda Population size used for offspring.
+ * @param isActive Covariance update mechanism.
+ * @param diagonalOnly Simplified covariance update.
+ * @param stopValue Termination criteria for optimization.
+ * @param fTol Tolerance relative error on the objective function.
+ * @param pointTol Tolerance for checking that the optimum is correct.
+ * @param maxEvaluations Maximum number of evaluations.
+ * @param expected Expected point / value.
+ */
+ private void doTest(MultivariateRealFunction func,
+ double[] startPoint,
+ double[] inSigma,
+ double[][] boundaries,
+ GoalType goal,
+ int lambda,
+ boolean isActive,
+ int diagonalOnly,
+ double stopValue,
+ double fTol,
+ double pointTol,
+ int maxEvaluations,
+ RealPointValuePair expected)
+ throws MathException {
+ int dim = startPoint.length;
+ // test diagonalOnly = 0 - slow but normally fewer feval#
+ MultivariateRealOptimizer optim =
+ new CMAESOptimizer(
+ lambda, inSigma, boundaries, 30000,
+ stopValue, isActive, diagonalOnly, 0, new MersenneTwister(),false);
+ RealPointValuePair result = optim.optimize(maxEvaluations, func, goal, startPoint);
+ Assert.assertEquals(expected.getValue(),
+ result.getValue(), fTol);
+ for (int i = 0; i < dim; i++) {
+ Assert.assertEquals(expected.getPoint()[i],
+ result.getPoint()[i], pointTol);
+ }
+ }
+
+ private static double[] point(int n, double value) {
+ double[] ds = new double[n];
+ Arrays.fill(ds, value);
+ return ds;
+ }
+
+ private static double[][] boundaries(int dim,
+ double lower, double upper) {
+ double[][] boundaries = new double[2][dim];
+ for (int i = 0; i < dim; i++)
+ boundaries[0][i] = lower;
+ for (int i = 0; i < dim; i++)
+ boundaries[1][i] = upper;
+ return boundaries;
+ }
+
+ private static class Sphere implements MultivariateRealFunction {
+
+ public double value(double[] x) {
+ double f = 0;
+ for (int i = 0; i < x.length; ++i)
+ f += x[i] * x[i];
+ return f;
+ }
+ }
+
+ private static class Cigar implements MultivariateRealFunction {
+ private double factor;
+
+ Cigar() {
+ this(1e3);
+ }
+
+ Cigar(double axisratio) {
+ factor = axisratio * axisratio;
+ }
+
+ public double value(double[] x) {
+ double f = x[0] * x[0];
+ for (int i = 1; i < x.length; ++i)
+ f += factor * x[i] * x[i];
+ return f;
+ }
+ }
+
+ private static class Tablet implements MultivariateRealFunction {
+ private double factor;
+
+ Tablet() {
+ this(1e3);
+ }
+
+ Tablet(double axisratio) {
+ factor = axisratio * axisratio;
+ }
+
+ public double value(double[] x) {
+ double f = factor * x[0] * x[0];
+ for (int i = 1; i < x.length; ++i)
+ f += x[i] * x[i];
+ return f;
+ }
+ }
+
+ private static class CigTab implements MultivariateRealFunction {
+ private double factor;
+
+ CigTab() {
+ this(1e4);
+ }
+
+ CigTab(double axisratio) {
+ factor = axisratio;
+ }
+
+ public double value(double[] x) {
+ int end = x.length - 1;
+ double f = x[0] * x[0] / factor + factor * x[end] * x[end];
+ for (int i = 1; i < end; ++i)
+ f += x[i] * x[i];
+ return f;
+ }
+ }
+
+ private static class TwoAxes implements MultivariateRealFunction {
+
+ private double factor;
+
+ TwoAxes() {
+ this(1e6);
+ }
+
+ TwoAxes(double axisratio) {
+ factor = axisratio * axisratio;
+ }
+
+ public double value(double[] x) {
+ double f = 0;
+ for (int i = 0; i < x.length; ++i)
+ f += (i < x.length / 2 ? factor : 1) * x[i] * x[i];
+ return f;
+ }
+ }
+
+ private static class ElliRotated implements MultivariateRealFunction {
+ private Basis B = new Basis();
+ private double factor;
+
+ ElliRotated() {
+ this(1e3);
+ }
+
+ ElliRotated(double axisratio) {
+ factor = axisratio * axisratio;
+ }
+
+ public double value(double[] x) {
+ double f = 0;
+ x = B.Rotate(x);
+ for (int i = 0; i < x.length; ++i)
+ f += Math.pow(factor, i / (x.length - 1.)) * x[i] * x[i];
+ return f;
+ }
+ }
+
+ private static class Elli implements MultivariateRealFunction {
+
+ private double factor;
+
+ Elli() {
+ this(1e3);
+ }
+
+ Elli(double axisratio) {
+ factor = axisratio * axisratio;
+ }
+
+ public double value(double[] x) {
+ double f = 0;
+ for (int i = 0; i < x.length; ++i)
+ f += Math.pow(factor, i / (x.length - 1.)) * x[i] * x[i];
+ return f;
+ }
+ }
+
+ private static class MinusElli implements MultivariateRealFunction {
+
+ public double value(double[] x) {
+ return 1.0-(new Elli().value(x));
+ }
+ }
+
+ private static class DiffPow implements MultivariateRealFunction {
+
+ public double value(double[] x) {
+ double f = 0;
+ for (int i = 0; i < x.length; ++i)
+ f += Math.pow(Math.abs(x[i]), 2. + 10 * (double) i
+ / (x.length - 1.));
+ return f;
+ }
+ }
+
+ private static class SsDiffPow implements MultivariateRealFunction {
+
+ public double value(double[] x) {
+ double f = Math.pow(new DiffPow().value(x), 0.25);
+ return f;
+ }
+ }
+
+ private static class Rosen implements MultivariateRealFunction {
+
+ public double value(double[] x) {
+ double f = 0;
+ for (int i = 0; i < x.length - 1; ++i)
+ f += 1e2 * (x[i] * x[i] - x[i + 1]) * (x[i] * x[i] - x[i + 1])
+ + (x[i] - 1.) * (x[i] - 1.);
+ return f;
+ }
+ }
+
+ private static class Ackley implements MultivariateRealFunction {
+ private double axisratio;
+
+ Ackley(double axra) {
+ axisratio = axra;
+ }
+
+ public Ackley() {
+ this(1);
+ }
+
+ public double value(double[] x) {
+ double f = 0;
+ double res2 = 0;
+ double fac = 0;
+ for (int i = 0; i < x.length; ++i) {
+ fac = Math.pow(axisratio, (i - 1.) / (x.length - 1.));
+ f += fac * fac * x[i] * x[i];
+ res2 += Math.cos(2. * Math.PI * fac * x[i]);
+ }
+ f = (20. - 20. * Math.exp(-0.2 * Math.sqrt(f / x.length))
+ + Math.exp(1.) - Math.exp(res2 / x.length));
+ return f;
+ }
+ }
+
+ private static class Rastrigin implements MultivariateRealFunction {
+
+ private double axisratio;
+ private double amplitude;
+
+ Rastrigin() {
+ this(1, 10);
+ }
+
+ Rastrigin(double axisratio, double amplitude) {
+ this.axisratio = axisratio;
+ this.amplitude = amplitude;
+ }
+
+ public double value(double[] x) {
+ double f = 0;
+ double fac;
+ for (int i = 0; i < x.length; ++i) {
+ fac = Math.pow(axisratio, (i - 1.) / (x.length - 1.));
+ if (i == 0 && x[i] < 0)
+ fac *= 1.;
+ f += fac * fac * x[i] * x[i] + amplitude
+ * (1. - Math.cos(2. * Math.PI * fac * x[i]));
+ }
+ return f;
+ }
+ }
+
+ private static class Basis {
+ double[][] basis;
+ Random rand = new Random(2); // use not always the same basis
+
+ double[] Rotate(double[] x) {
+ GenBasis(x.length);
+ double[] y = new double[x.length];
+ for (int i = 0; i < x.length; ++i) {
+ y[i] = 0;
+ for (int j = 0; j < x.length; ++j)
+ y[i] += basis[i][j] * x[j];
+ }
+ return y;
+ }
+
+ void GenBasis(int DIM) {
+ if (basis != null ? basis.length == DIM : false)
+ return;
+
+ double sp;
+ int i, j, k;
+
+ /* generate orthogonal basis */
+ basis = new double[DIM][DIM];
+ for (i = 0; i < DIM; ++i) {
+ /* sample components gaussian */
+ for (j = 0; j < DIM; ++j)
+ basis[i][j] = rand.nextGaussian();
+ /* substract projection of previous vectors */
+ for (j = i - 1; j >= 0; --j) {
+ for (sp = 0., k = 0; k < DIM; ++k)
+ sp += basis[i][k] * basis[j][k]; /* scalar product */
+ for (k = 0; k < DIM; ++k)
+ basis[i][k] -= sp * basis[j][k]; /* substract */
+ }
+ /* normalize */
+ for (sp = 0., k = 0; k < DIM; ++k)
+ sp += basis[i][k] * basis[i][k]; /* squared norm */
+ for (k = 0; k < DIM; ++k)
+ basis[i][k] /= Math.sqrt(sp);
+ }
+ }
+ }
+}
Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/direct/CMAESOptimizerTest.java
------------------------------------------------------------------------------
svn:eol-style = native