http://git-wip-us.apache.org/repos/asf/commons-math/blob/b4669aad/src/main/java/org/apache/commons/math4/optimization/direct/CMAESOptimizer.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/optimization/direct/CMAESOptimizer.java b/src/main/java/org/apache/commons/math4/optimization/direct/CMAESOptimizer.java deleted file mode 100644 index 17d84af..0000000 --- a/src/main/java/org/apache/commons/math4/optimization/direct/CMAESOptimizer.java +++ /dev/null @@ -1,1441 +0,0 @@ -/* - * 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.math4.optimization.direct; - -import java.util.ArrayList; -import java.util.Arrays; -import java.util.List; - -import org.apache.commons.math4.analysis.MultivariateFunction; -import org.apache.commons.math4.exception.DimensionMismatchException; -import org.apache.commons.math4.exception.NotPositiveException; -import org.apache.commons.math4.exception.NotStrictlyPositiveException; -import org.apache.commons.math4.exception.OutOfRangeException; -import org.apache.commons.math4.exception.TooManyEvaluationsException; -import org.apache.commons.math4.linear.Array2DRowRealMatrix; -import org.apache.commons.math4.linear.EigenDecomposition; -import org.apache.commons.math4.linear.MatrixUtils; -import org.apache.commons.math4.linear.RealMatrix; -import org.apache.commons.math4.optimization.ConvergenceChecker; -import org.apache.commons.math4.optimization.GoalType; -import org.apache.commons.math4.optimization.MultivariateOptimizer; -import org.apache.commons.math4.optimization.OptimizationData; -import org.apache.commons.math4.optimization.PointValuePair; -import org.apache.commons.math4.optimization.SimpleValueChecker; -import org.apache.commons.math4.random.MersenneTwister; -import org.apache.commons.math4.random.RandomGenerator; -import org.apache.commons.math4.util.FastMath; -import org.apache.commons.math4.util.MathArrays; - -/** - * <p>An implementation of 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 - * on the underlying search space. Unlike a quasi-Newton method, the - * CMA-ES neither estimates nor uses gradients, making it considerably more - * reliable in terms of finding a good, or even close to optimal, solution.</p> - * - * <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> - * - * <p>This implementation is translated and adapted from the Matlab version - * of the CMA-ES algorithm as implemented in module {@code cmaes.m} version 3.51.</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> - * - * @deprecated As of 3.1 (to be removed in 4.0). - * @since 3.0 - */ -@Deprecated -public class CMAESOptimizer - extends BaseAbstractMultivariateSimpleBoundsOptimizer<MultivariateFunction> - implements MultivariateOptimizer { - /** 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 linearly 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. - */ - private int checkFeasableCount; - /** - * @see Sigma - */ - 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 - * - * @deprecated As of version 3.1: Parameter {@code lambda} must be - * passed with the call to {@link #optimize(int,MultivariateFunction,GoalType,OptimizationData[]) - * optimize} (whereas in the current code it is set to an undocumented value). - */ - @Deprecated - public CMAESOptimizer() { - this(0); - } - - /** - * @param lambda Population size. - * @deprecated As of version 3.1: Parameter {@code lambda} must be - * passed with the call to {@link #optimize(int,MultivariateFunction,GoalType,OptimizationData[]) - * optimize} (whereas in the current code it is set to an undocumented value).. - */ - @Deprecated - public CMAESOptimizer(int lambda) { - this(lambda, null, DEFAULT_MAXITERATIONS, DEFAULT_STOPFITNESS, - DEFAULT_ISACTIVECMA, DEFAULT_DIAGONALONLY, - DEFAULT_CHECKFEASABLECOUNT, DEFAULT_RANDOMGENERATOR, - false, null); - } - - /** - * @param lambda Population size. - * @param inputSigma Initial standard deviations to sample new points - * around the initial guess. - * @deprecated As of version 3.1: Parameters {@code lambda} and {@code inputSigma} must be - * passed with the call to {@link #optimize(int,MultivariateFunction,GoalType,OptimizationData[]) - * optimize}. - */ - @Deprecated - public CMAESOptimizer(int lambda, double[] inputSigma) { - this(lambda, inputSigma, DEFAULT_MAXITERATIONS, DEFAULT_STOPFITNESS, - DEFAULT_ISACTIVECMA, DEFAULT_DIAGONALONLY, - DEFAULT_CHECKFEASABLECOUNT, DEFAULT_RANDOMGENERATOR, false); - } - - /** - * @param lambda Population size. - * @param inputSigma Initial standard deviations to sample new points - * around the initial guess. - * @param maxIterations Maximal number of iterations. - * @param stopFitness Whether to stop if objective function value is smaller than - * {@code 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 Random generator. - * @param generateStatistics Whether statistic data is collected. - * @deprecated See {@link SimpleValueChecker#SimpleValueChecker()} - */ - @Deprecated - public CMAESOptimizer(int lambda, double[] inputSigma, - int maxIterations, double stopFitness, - boolean isActiveCMA, int diagonalOnly, int checkFeasableCount, - RandomGenerator random, boolean generateStatistics) { - this(lambda, inputSigma, maxIterations, stopFitness, isActiveCMA, - diagonalOnly, checkFeasableCount, random, generateStatistics, - new SimpleValueChecker()); - } - - /** - * @param lambda Population size. - * @param inputSigma Initial standard deviations to sample new points - * around the initial guess. - * @param maxIterations Maximal number of iterations. - * @param stopFitness Whether to stop if objective function value is smaller than - * {@code 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 Random generator. - * @param generateStatistics Whether statistic data is collected. - * @param checker Convergence checker. - * @deprecated As of version 3.1: Parameters {@code lambda} and {@code inputSigma} must be - * passed with the call to {@link #optimize(int,MultivariateFunction,GoalType,OptimizationData[]) - * optimize}. - */ - @Deprecated - public CMAESOptimizer(int lambda, double[] inputSigma, - int maxIterations, double stopFitness, - boolean isActiveCMA, int diagonalOnly, int checkFeasableCount, - RandomGenerator random, boolean generateStatistics, - ConvergenceChecker<PointValuePair> checker) { - super(checker); - this.lambda = lambda; - this.inputSigma = inputSigma == null ? null : (double[]) inputSigma.clone(); - this.maxIterations = maxIterations; - this.stopFitness = stopFitness; - this.isActiveCMA = isActiveCMA; - this.diagonalOnly = diagonalOnly; - this.checkFeasableCount = checkFeasableCount; - this.random = random; - this.generateStatistics = generateStatistics; - } - - /** - * @param maxIterations Maximal number of iterations. - * @param stopFitness Whether to stop if objective function value is smaller than - * {@code 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 Random generator. - * @param generateStatistics Whether statistic data is collected. - * @param checker Convergence checker. - * - * @since 3.1 - */ - public CMAESOptimizer(int maxIterations, - double stopFitness, - boolean isActiveCMA, - int diagonalOnly, - int checkFeasableCount, - RandomGenerator random, - boolean generateStatistics, - ConvergenceChecker<PointValuePair> checker) { - super(checker); - 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; - } - - /** - * Input sigma values. - * They define the initial coordinate-wise standard deviations for - * sampling new search points around the initial guess. - * It is suggested to set them to the estimated distance from the - * initial to the desired optimum. - * Small values induce the search to be more local (and very small - * values are more likely to find a local optimum close to the initial - * guess). - * Too small values might however lead to early termination. - * @since 3.1 - */ - public static class Sigma implements OptimizationData { - /** Sigma values. */ - private final double[] sigma; - - /** - * @param s Sigma values. - * @throws NotPositiveException if any of the array entries is smaller - * than zero. - */ - public Sigma(double[] s) - throws NotPositiveException { - for (int i = 0; i < s.length; i++) { - if (s[i] < 0) { - throw new NotPositiveException(s[i]); - } - } - - sigma = s.clone(); - } - - /** - * @return the sigma values. - */ - public double[] getSigma() { - return sigma.clone(); - } - } - - /** - * Population size. - * The number of offspring is the primary strategy parameter. - * In the absence of better clues, a good default could be an - * integer close to {@code 4 + 3 ln(n)}, where {@code n} is the - * number of optimized parameters. - * Increasing the population size improves global search properties - * at the expense of speed (which in general decreases at most - * linearly with increasing population size). - * @since 3.1 - */ - public static class PopulationSize implements OptimizationData { - /** Population size. */ - private final int lambda; - - /** - * @param size Population size. - * @throws NotStrictlyPositiveException if {@code size <= 0}. - */ - public PopulationSize(int size) - throws NotStrictlyPositiveException { - if (size <= 0) { - throw new NotStrictlyPositiveException(size); - } - lambda = size; - } - - /** - * @return the population size. - */ - public int getPopulationSize() { - return lambda; - } - } - - /** - * Optimize an objective function. - * - * @param maxEval Allowed number of evaluations of the objective function. - * @param f Objective function. - * @param goalType Optimization type. - * @param optData Optimization data. The following data will be looked for: - * <ul> - * <li>{@link org.apache.commons.math4.optimization.InitialGuess InitialGuess}</li> - * <li>{@link Sigma}</li> - * <li>{@link PopulationSize}</li> - * </ul> - * @return the point/value pair giving the optimal value for objective - * function. - */ - @Override - protected PointValuePair optimizeInternal(int maxEval, MultivariateFunction f, - GoalType goalType, - OptimizationData... optData) { - // Scan "optData" for the input specific to this optimizer. - parseOptimizationData(optData); - - // The parent's method will retrieve the common parameters from - // "optData" and call "doOptimize". - return super.optimizeInternal(maxEval, f, goalType, optData); - } - - /** {@inheritDoc} */ - @Override - protected PointValuePair doOptimize() { - checkParameters(); - // -------------------- Initialization -------------------------------- - isMinimize = getGoalType().equals(GoalType.MINIMIZE); - final FitnessFunction fitfun = new FitnessFunction(); - final double[] guess = getStartPoint(); - // number of objective variables/problem dimension - dimension = guess.length; - initializeCMA(guess); - iterations = 0; - double bestValue = fitfun.value(guess); - push(fitnessHistory, bestValue); - PointValuePair optimum = new PointValuePair(getStartPoint(), - isMinimize ? bestValue : -bestValue); - PointValuePair lastResult = null; - - // -------------------- Generation Loop -------------------------------- - - generationLoop: - for (iterations = 1; iterations <= maxIterations; iterations++) { - // Generate and evaluate lambda offspring - final RealMatrix arz = randn1(dimension, lambda); - final RealMatrix arx = zeros(dimension, lambda); - final 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 - final int[] arindex = sortedIndices(fitness); - // Calculate new xmean, this is selection and recombination - final RealMatrix xold = xmean; // for speed up of Eq. (2) and (3) - final RealMatrix bestArx = selectColumns(arx, MathArrays.copyOf(arindex, mu)); - xmean = bestArx.multiply(weights); - final RealMatrix bestArz = selectColumns(arz, MathArrays.copyOf(arindex, mu)); - final RealMatrix zmean = bestArz.multiply(weights); - final boolean hsig = updateEvolutionPaths(zmean, xold); - if (diagonalOnly <= 0) { - updateCovariance(hsig, bestArx, arz, arindex, xold); - } else { - updateCovarianceDiagonalOnly(hsig, bestArz); - } - // Adapt step size sigma - Eq. (5) - sigma *= FastMath.exp(FastMath.min(1, (normps/chiN - 1) * cs / damps)); - final double bestFitness = fitness[arindex[0]]; - final double worstFitness = fitness[arindex[arindex.length - 1]]; - if (bestValue > bestFitness) { - bestValue = bestFitness; - lastResult = optimum; - optimum = new PointValuePair(fitfun.repair(bestArx.getColumn(0)), - isMinimize ? bestFitness : -bestFitness); - if (getConvergenceChecker() != null && lastResult != null && - getConvergenceChecker().converged(iterations, optimum, lastResult)) { - break generationLoop; - } - } - // handle termination criteria - // Break, if fitness is good enough - if (stopFitness != 0 && bestFitness < (isMinimize ? stopFitness : -stopFitness)) { - break generationLoop; - } - final double[] sqrtDiagC = sqrt(diagC).getColumn(0); - final double[] pcCol = pc.getColumn(0); - for (int i = 0; i < dimension; i++) { - if (sigma * FastMath.max(FastMath.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; - } - } - final double historyBest = min(fitnessHistory); - final double historyWorst = max(fitnessHistory); - if (iterations > 2 && - FastMath.max(historyWorst, worstFitness) - - FastMath.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) { - final PointValuePair current - = new PointValuePair(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 *= FastMath.exp(0.2 + cs / damps); - } - if (iterations > 2 && FastMath.max(historyWorst, bestFitness) - - FastMath.min(historyBest, bestFitness) == 0) { - sigma *= FastMath.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; - } - - /** - * Scans the list of (required and optional) optimization data that - * characterize the problem. - * - * @param optData Optimization data. The following data will be looked for: - * <ul> - * <li>{@link Sigma}</li> - * <li>{@link PopulationSize}</li> - * </ul> - */ - private void parseOptimizationData(OptimizationData... optData) { - // The existing values (as set by the previous call) are reused if - // not provided in the argument list. - for (OptimizationData data : optData) { - if (data instanceof Sigma) { - inputSigma = ((Sigma) data).getSigma(); - continue; - } - if (data instanceof PopulationSize) { - lambda = ((PopulationSize) data).getPopulationSize(); - continue; - } - } - } - - /** - * Checks dimensions and values of boundaries and inputSigma if defined. - */ - private void checkParameters() { - final double[] init = getStartPoint(); - final double[] lB = getLowerBound(); - final double[] uB = getUpperBound(); - - if (inputSigma != null) { - if (inputSigma.length != init.length) { - throw new DimensionMismatchException(inputSigma.length, init.length); - } - for (int i = 0; i < init.length; i++) { - if (inputSigma[i] < 0) { - // XXX Remove this block in 4.0 (check performed in "Sigma" class). - throw new NotPositiveException(inputSigma[i]); - } - if (inputSigma[i] > uB[i] - lB[i]) { - throw new OutOfRangeException(inputSigma[i], 0, uB[i] - lB[i]); - } - } - } - } - - /** - * 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) { - // XXX Line below to replace the current one in 4.0 (MATH-879). - // throw new NotStrictlyPositiveException(lambda); - lambda = 4 + (int) (3 * FastMath.log(dimension)); - } - // initialize sigma - final double[][] sigmaArray = new double[guess.length][1]; - for (int i = 0; i < guess.length; i++) { - // XXX Line below to replace the current one in 4.0 (MATH-868). - // sigmaArray[i][0] = inputSigma[i]; - sigmaArray[i][0] = inputSigma == null ? 0.3 : inputSigma[i]; - } - final 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 = FastMath.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 * FastMath.max(0, FastMath.sqrt((mueff - 1) / - (dimension + 1)) - 1)) * - FastMath.max(0.3, - 1 - dimension / (1e-6 + maxIterations)) + cs; // minor increment - ccov1 = 2 / ((dimension + 1.3) * (dimension + 1.3) + mueff); - ccovmu = FastMath.min(1 - ccov1, 2 * (mueff - 2 + 1 / mueff) / - ((dimension + 2) * (dimension + 2) + mueff)); - ccov1Sep = FastMath.min(1, ccov1 * (dimension + 1.5) / 3); - ccovmuSep = FastMath.min(1 - ccov1, ccovmu * (dimension + 1.5) / 3); - chiN = FastMath.sqrt(dimension) * - (1 - 1 / ((double) 4 * dimension) + 1 / ((double) 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 = ps.getFrobeniusNorm(); - - 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 / (double) 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(FastMath.sqrt(cs * (2 - cs) * mueff))); - normps = ps.getFrobeniusNorm(); - final boolean hsig = normps / - FastMath.sqrt(1 - FastMath.pow(1 - cs, 2 * iterations)) / - chiN < 1.4 + 2 / ((double) dimension + 1); - pc = pc.scalarMultiply(1 - cc); - if (hsig) { - pc = pc.add(xmean.subtract(xold).scalarMultiply(FastMath.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. - */ - private void updateCovarianceDiagonalOnly(boolean hsig, - final RealMatrix bestArz) { - // minor correction if hsig==false - double oldFac = hsig ? 0 : ccov1Sep * cc * (2 - cc); - oldFac += 1 - ccov1Sep - ccovmuSep; - diagC = diagC.scalarMultiply(oldFac) // regard old matrix - .add(square(pc).scalarMultiply(ccov1Sep)) // plus rank one update - .add((times(diagC, square(bestArz).multiply(weights))) // plus rank mu update - .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) { - final RealMatrix arpos = bestArx.subtract(repmat(xold, 1, mu)) - .scalarMultiply(1 / sigma); // mu difference vectors - final 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 / (FastMath.pow(dimension + 2, 1.5) + 2 * mueff); - // keep at least 0.66 in all directions, small popsize are most - // critical - final double negminresidualvariance = 0.66; - // where to make up for the variance loss - final double negalphaold = 0.5; - // prepare vectors, compute negative updating matrix Cneg - final int[] arReverseIndex = reverse(arindex); - RealMatrix arzneg = selectColumns(arz, MathArrays.copyOf(arReverseIndex, mu)); - RealMatrix arnorms = sqrt(sumRows(square(arzneg))); - final int[] idxnorms = sortedIndices(arnorms.getRow(0)); - final RealMatrix arnormsSorted = selectColumns(arnorms, idxnorms); - final int[] idxReverse = reverse(idxnorms); - final RealMatrix arnormsReverse = selectColumns(arnorms, idxReverse); - arnorms = divide(arnormsReverse, arnormsSorted); - final int[] idxInv = inverse(idxnorms); - final RealMatrix arnormsInv = selectColumns(arnorms, idxInv); - // check and set learning rate negccov - final double negcovMax = (1 - negminresidualvariance) / - square(arnormsInv).multiply(weights).getEntry(0, 0); - if (negccov > negcovMax) { - negccov = negcovMax; - } - arzneg = times(arzneg, repmat(arnormsInv, dimension, 1)); - final RealMatrix artmp = BD.multiply(arzneg); - final RealMatrix Cneg = artmp.multiply(diag(weights)).multiply(artmp.transpose()); - oldFac += negalphaold * negccov; - C = C.scalarMultiply(oldFac) - .add(roneu) // regard old matrix - .add(arpos.scalarMultiply( // plus rank one update - ccovmu + (1 - negalphaold) * negccov) // plus rank mu update - .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 - final EigenDecomposition eig = new EigenDecomposition(C); - 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); - } - } - final 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)) { - final 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 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 a sorted array of indices pointing into doubles. - */ - private int[] sortedIndices(final double[] doubles) { - final DoubleIndex[] dis = new DoubleIndex[doubles.length]; - for (int i = 0; i < doubles.length; i++) { - dis[i] = new DoubleIndex(doubles[i], i); - } - Arrays.sort(dis); - final 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 final double value; - /** Index into sorted array. */ - private final 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); - } - - /** {@inheritDoc} */ - @Override - public boolean equals(Object other) { - - if (this == other) { - return true; - } - - if (other instanceof DoubleIndex) { - return Double.compare(value, ((DoubleIndex) other).value) == 0; - } - - return false; - } - - /** {@inheritDoc} */ - @Override - public int hashCode() { - long bits = Double.doubleToLongBits(value); - return (int) ((1438542 ^ (bits >>> 32) ^ bits) & 0xffffffff); - } - } - - /** - * 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 { - /** Determines the penalty for boundary violations */ - private double valueRange; - /** - * Flag indicating whether the objective variables are forced into their - * bounds if defined - */ - private final boolean isRepairMode; - - /** Simple constructor. - */ - public FitnessFunction() { - valueRange = 1; - isRepairMode = true; - } - - /** - * @param point Normalized objective variables. - * @return the objective value + penalty for violated bounds. - */ - public double value(final double[] point) { - double value; - if (isRepairMode) { - double[] repaired = repair(point); - value = CMAESOptimizer.this.computeObjectiveValue(repaired) + - penalty(point, repaired); - } else { - value = CMAESOptimizer.this.computeObjectiveValue(point); - } - return isMinimize ? value : -value; - } - - /** - * @param x Normalized objective variables. - * @return {@code true} if in bounds. - */ - public boolean isFeasible(final double[] x) { - final double[] lB = CMAESOptimizer.this.getLowerBound(); - final double[] uB = CMAESOptimizer.this.getUpperBound(); - - for (int i = 0; i < x.length; i++) { - if (x[i] < lB[i]) { - return false; - } - if (x[i] > uB[i]) { - return false; - } - } - return true; - } - - /** - * @param valueRange Adjusts the penalty computation. - */ - public void setValueRange(double valueRange) { - this.valueRange = valueRange; - } - - /** - * @param x Normalized objective variables. - * @return the repaired (i.e. all in bounds) objective variables. - */ - private double[] repair(final double[] x) { - final double[] lB = CMAESOptimizer.this.getLowerBound(); - final double[] uB = CMAESOptimizer.this.getUpperBound(); - - final double[] repaired = new double[x.length]; - for (int i = 0; i < x.length; i++) { - if (x[i] < lB[i]) { - repaired[i] = lB[i]; - } else if (x[i] > uB[i]) { - repaired[i] = uB[i]; - } 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 = FastMath.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) { - final 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] = FastMath.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) { - final 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] = FastMath.sqrt(m.getEntry(r, c)); - } - } - return new Array2DRowRealMatrix(d, false); - } - - /** - * @param m Input matrix. - * @return Matrix representing the element-wise square of m. - */ - private static RealMatrix square(final RealMatrix m) { - final 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 the matrix where the elements of m and n are element-wise multiplied. - */ - private static RealMatrix times(final RealMatrix m, final RealMatrix n) { - final 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 n are element-wise divided. - */ - private static RealMatrix divide(final RealMatrix m, final RealMatrix n) { - final 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) { - final 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) { - final 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 Row matrix representing the sums of the rows. - */ - private static RealMatrix sumRows(final RealMatrix m) { - final 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 the diagonal n-by-n matrix if m is a column matrix or the column - * matrix representing the diagonal if m is a n-by-n matrix. - */ - private static RealMatrix diag(final RealMatrix m) { - if (m.getColumnDimension() == 1) { - final 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 { - final 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 column from m1 to m2. - * - * @param m1 Source matrix. - * @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-by-m matrix filled with 1. - */ - private static RealMatrix ones(int n, int m) { - final double[][] d = new double[n][m]; - for (int r = 0; r < n; r++) { - Arrays.fill(d[r], 1); - } - return new Array2DRowRealMatrix(d, false); - } - - /** - * @param n Number of rows. - * @param m Number of columns. - * @return n-by-m matrix of 0 values out of diagonal, and 1 values on - * the diagonal. - */ - private static RealMatrix eye(int n, int m) { - final 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-by-m matrix of zero 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 a matrix which replicates the input matrix in both directions. - */ - private static RealMatrix repmat(final RealMatrix mat, int n, int m) { - final int rd = mat.getRowDimension(); - final int cd = mat.getColumnDimension(); - final 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 a sequence as column matrix. - */ - private static RealMatrix sequence(double start, double end, double step) { - final int size = (int) ((end - start) / step + 1); - final 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 the maximum of the 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 the minimum of the 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 the maximum of the 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 the minimum of the 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 the inverse of the mapping defined by indices. - */ - private static int[] inverse(final int[] indices) { - final 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 the indices in inverse order (last is first). - */ - private static int[] reverse(final int[] indices) { - final 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 an array of Gaussian random numbers. - */ - private double[] randn(int size) { - final 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 a 2-dimensional matrix of Gaussian random numbers. - */ - private RealMatrix randn1(int size, int popSize) { - final 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); - } -}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/b4669aad/src/main/java/org/apache/commons/math4/optimization/direct/MultiDirectionalSimplex.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/optimization/direct/MultiDirectionalSimplex.java b/src/main/java/org/apache/commons/math4/optimization/direct/MultiDirectionalSimplex.java deleted file mode 100644 index cdc0bab..0000000 --- a/src/main/java/org/apache/commons/math4/optimization/direct/MultiDirectionalSimplex.java +++ /dev/null @@ -1,218 +0,0 @@ -/* - * 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.math4.optimization.direct; - -import java.util.Comparator; - -import org.apache.commons.math4.analysis.MultivariateFunction; -import org.apache.commons.math4.optimization.PointValuePair; - -/** - * This class implements the multi-directional direct search method. - * - * @deprecated As of 3.1 (to be removed in 4.0). - * @since 3.0 - */ -@Deprecated -public class MultiDirectionalSimplex extends AbstractSimplex { - /** Default value for {@link #khi}: {@value}. */ - private static final double DEFAULT_KHI = 2; - /** Default value for {@link #gamma}: {@value}. */ - private static final double DEFAULT_GAMMA = 0.5; - /** Expansion coefficient. */ - private final double khi; - /** Contraction coefficient. */ - private final double gamma; - - /** - * Build a multi-directional simplex with default coefficients. - * The default values are 2.0 for khi and 0.5 for gamma. - * - * @param n Dimension of the simplex. - */ - public MultiDirectionalSimplex(final int n) { - this(n, 1d); - } - - /** - * Build a multi-directional simplex with default coefficients. - * The default values are 2.0 for khi and 0.5 for gamma. - * - * @param n Dimension of the simplex. - * @param sideLength Length of the sides of the default (hypercube) - * simplex. See {@link AbstractSimplex#AbstractSimplex(int,double)}. - */ - public MultiDirectionalSimplex(final int n, double sideLength) { - this(n, sideLength, DEFAULT_KHI, DEFAULT_GAMMA); - } - - /** - * Build a multi-directional simplex with specified coefficients. - * - * @param n Dimension of the simplex. See - * {@link AbstractSimplex#AbstractSimplex(int,double)}. - * @param khi Expansion coefficient. - * @param gamma Contraction coefficient. - */ - public MultiDirectionalSimplex(final int n, - final double khi, final double gamma) { - this(n, 1d, khi, gamma); - } - - /** - * Build a multi-directional simplex with specified coefficients. - * - * @param n Dimension of the simplex. See - * {@link AbstractSimplex#AbstractSimplex(int,double)}. - * @param sideLength Length of the sides of the default (hypercube) - * simplex. See {@link AbstractSimplex#AbstractSimplex(int,double)}. - * @param khi Expansion coefficient. - * @param gamma Contraction coefficient. - */ - public MultiDirectionalSimplex(final int n, double sideLength, - final double khi, final double gamma) { - super(n, sideLength); - - this.khi = khi; - this.gamma = gamma; - } - - /** - * Build a multi-directional simplex with default coefficients. - * The default values are 2.0 for khi and 0.5 for gamma. - * - * @param steps Steps along the canonical axes representing box edges. - * They may be negative but not zero. See - */ - public MultiDirectionalSimplex(final double[] steps) { - this(steps, DEFAULT_KHI, DEFAULT_GAMMA); - } - - /** - * Build a multi-directional simplex with specified coefficients. - * - * @param steps Steps along the canonical axes representing box edges. - * They may be negative but not zero. See - * {@link AbstractSimplex#AbstractSimplex(double[])}. - * @param khi Expansion coefficient. - * @param gamma Contraction coefficient. - */ - public MultiDirectionalSimplex(final double[] steps, - final double khi, final double gamma) { - super(steps); - - this.khi = khi; - this.gamma = gamma; - } - - /** - * Build a multi-directional simplex with default coefficients. - * The default values are 2.0 for khi and 0.5 for gamma. - * - * @param referenceSimplex Reference simplex. See - * {@link AbstractSimplex#AbstractSimplex(double[][])}. - */ - public MultiDirectionalSimplex(final double[][] referenceSimplex) { - this(referenceSimplex, DEFAULT_KHI, DEFAULT_GAMMA); - } - - /** - * Build a multi-directional simplex with specified coefficients. - * - * @param referenceSimplex Reference simplex. See - * {@link AbstractSimplex#AbstractSimplex(double[][])}. - * @param khi Expansion coefficient. - * @param gamma Contraction coefficient. - * @throws org.apache.commons.math4.exception.NotStrictlyPositiveException - * if the reference simplex does not contain at least one point. - * @throws org.apache.commons.math4.exception.DimensionMismatchException - * if there is a dimension mismatch in the reference simplex. - */ - public MultiDirectionalSimplex(final double[][] referenceSimplex, - final double khi, final double gamma) { - super(referenceSimplex); - - this.khi = khi; - this.gamma = gamma; - } - - /** {@inheritDoc} */ - @Override - public void iterate(final MultivariateFunction evaluationFunction, - final Comparator<PointValuePair> comparator) { - // Save the original simplex. - final PointValuePair[] original = getPoints(); - final PointValuePair best = original[0]; - - // Perform a reflection step. - final PointValuePair reflected = evaluateNewSimplex(evaluationFunction, - original, 1, comparator); - if (comparator.compare(reflected, best) < 0) { - // Compute the expanded simplex. - final PointValuePair[] reflectedSimplex = getPoints(); - final PointValuePair expanded = evaluateNewSimplex(evaluationFunction, - original, khi, comparator); - if (comparator.compare(reflected, expanded) <= 0) { - // Keep the reflected simplex. - setPoints(reflectedSimplex); - } - // Keep the expanded simplex. - return; - } - - // Compute the contracted simplex. - evaluateNewSimplex(evaluationFunction, original, gamma, comparator); - - } - - /** - * Compute and evaluate a new simplex. - * - * @param evaluationFunction Evaluation function. - * @param original Original simplex (to be preserved). - * @param coeff Linear coefficient. - * @param comparator Comparator to use to sort simplex vertices from best - * to poorest. - * @return the best point in the transformed simplex. - * @throws org.apache.commons.math4.exception.TooManyEvaluationsException - * if the maximal number of evaluations is exceeded. - */ - private PointValuePair evaluateNewSimplex(final MultivariateFunction evaluationFunction, - final PointValuePair[] original, - final double coeff, - final Comparator<PointValuePair> comparator) { - final double[] xSmallest = original[0].getPointRef(); - // Perform a linear transformation on all the simplex points, - // except the first one. - setPoint(0, original[0]); - final int dim = getDimension(); - for (int i = 1; i < getSize(); i++) { - final double[] xOriginal = original[i].getPointRef(); - final double[] xTransformed = new double[dim]; - for (int j = 0; j < dim; j++) { - xTransformed[j] = xSmallest[j] + coeff * (xSmallest[j] - xOriginal[j]); - } - setPoint(i, new PointValuePair(xTransformed, Double.NaN, false)); - } - - // Evaluate the simplex. - evaluate(evaluationFunction, comparator); - - return getPoint(0); - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/b4669aad/src/main/java/org/apache/commons/math4/optimization/direct/MultivariateFunctionMappingAdapter.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/optimization/direct/MultivariateFunctionMappingAdapter.java b/src/main/java/org/apache/commons/math4/optimization/direct/MultivariateFunctionMappingAdapter.java deleted file mode 100644 index d246ed4..0000000 --- a/src/main/java/org/apache/commons/math4/optimization/direct/MultivariateFunctionMappingAdapter.java +++ /dev/null @@ -1,301 +0,0 @@ -/* - * 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.math4.optimization.direct; - -import org.apache.commons.math4.analysis.MultivariateFunction; -import org.apache.commons.math4.analysis.UnivariateFunction; -import org.apache.commons.math4.analysis.function.Logit; -import org.apache.commons.math4.analysis.function.Sigmoid; -import org.apache.commons.math4.exception.DimensionMismatchException; -import org.apache.commons.math4.exception.NumberIsTooSmallException; -import org.apache.commons.math4.util.FastMath; -import org.apache.commons.math4.util.MathUtils; - -/** - * <p>Adapter for mapping bounded {@link MultivariateFunction} to unbounded ones.</p> - * - * <p> - * This adapter can be used to wrap functions subject to simple bounds on - * parameters so they can be used by optimizers that do <em>not</em> directly - * support simple bounds. - * </p> - * <p> - * The principle is that the user function that will be wrapped will see its - * parameters bounded as required, i.e when its {@code value} method is called - * with argument array {@code point}, the elements array will fulfill requirement - * {@code lower[i] <= point[i] <= upper[i]} for all i. Some of the components - * may be unbounded or bounded only on one side if the corresponding bound is - * set to an infinite value. The optimizer will not manage the user function by - * itself, but it will handle this adapter and it is this adapter that will take - * care the bounds are fulfilled. The adapter {@link #value(double[])} method will - * be called by the optimizer with unbound parameters, and the adapter will map - * the unbounded value to the bounded range using appropriate functions like - * {@link Sigmoid} for double bounded elements for example. - * </p> - * <p> - * As the optimizer sees only unbounded parameters, it should be noted that the - * start point or simplex expected by the optimizer should be unbounded, so the - * user is responsible for converting his bounded point to unbounded by calling - * {@link #boundedToUnbounded(double[])} before providing them to the optimizer. - * For the same reason, the point returned by the {@link - * org.apache.commons.math4.optimization.BaseMultivariateOptimizer#optimize(int, - * MultivariateFunction, org.apache.commons.math4.optimization.GoalType, double[])} - * method is unbounded. So to convert this point to bounded, users must call - * {@link #unboundedToBounded(double[])} by themselves!</p> - * <p> - * This adapter is only a poor man solution to simple bounds optimization constraints - * that can be used with simple optimizers like {@link SimplexOptimizer} with {@link - * NelderMeadSimplex} or {@link MultiDirectionalSimplex}. A better solution is to use - * an optimizer that directly supports simple bounds like {@link CMAESOptimizer} or - * {@link BOBYQAOptimizer}. One caveat of this poor man solution is that behavior near - * the bounds may be numerically unstable as bounds are mapped from infinite values. - * Another caveat is that convergence values are evaluated by the optimizer with respect - * to unbounded variables, so there will be scales differences when converted to bounded - * variables. - * </p> - * - * @see MultivariateFunctionPenaltyAdapter - * - * @deprecated As of 3.1 (to be removed in 4.0). - * @since 3.0 - */ - -@Deprecated -public class MultivariateFunctionMappingAdapter implements MultivariateFunction { - - /** Underlying bounded function. */ - private final MultivariateFunction bounded; - - /** Mapping functions. */ - private final Mapper[] mappers; - - /** Simple constructor. - * @param bounded bounded function - * @param lower lower bounds for each element of the input parameters array - * (some elements may be set to {@code Double.NEGATIVE_INFINITY} for - * unbounded values) - * @param upper upper bounds for each element of the input parameters array - * (some elements may be set to {@code Double.POSITIVE_INFINITY} for - * unbounded values) - * @exception DimensionMismatchException if lower and upper bounds are not - * consistent, either according to dimension or to values - */ - public MultivariateFunctionMappingAdapter(final MultivariateFunction bounded, - final double[] lower, final double[] upper) { - - // safety checks - MathUtils.checkNotNull(lower); - MathUtils.checkNotNull(upper); - if (lower.length != upper.length) { - throw new DimensionMismatchException(lower.length, upper.length); - } - for (int i = 0; i < lower.length; ++i) { - // note the following test is written in such a way it also fails for NaN - if (!(upper[i] >= lower[i])) { - throw new NumberIsTooSmallException(upper[i], lower[i], true); - } - } - - this.bounded = bounded; - this.mappers = new Mapper[lower.length]; - for (int i = 0; i < mappers.length; ++i) { - if (Double.isInfinite(lower[i])) { - if (Double.isInfinite(upper[i])) { - // element is unbounded, no transformation is needed - mappers[i] = new NoBoundsMapper(); - } else { - // element is simple-bounded on the upper side - mappers[i] = new UpperBoundMapper(upper[i]); - } - } else { - if (Double.isInfinite(upper[i])) { - // element is simple-bounded on the lower side - mappers[i] = new LowerBoundMapper(lower[i]); - } else { - // element is double-bounded - mappers[i] = new LowerUpperBoundMapper(lower[i], upper[i]); - } - } - } - - } - - /** Map an array from unbounded to bounded. - * @param point unbounded value - * @return bounded value - */ - public double[] unboundedToBounded(double[] point) { - - // map unbounded input point to bounded point - final double[] mapped = new double[mappers.length]; - for (int i = 0; i < mappers.length; ++i) { - mapped[i] = mappers[i].unboundedToBounded(point[i]); - } - - return mapped; - - } - - /** Map an array from bounded to unbounded. - * @param point bounded value - * @return unbounded value - */ - public double[] boundedToUnbounded(double[] point) { - - // map bounded input point to unbounded point - final double[] mapped = new double[mappers.length]; - for (int i = 0; i < mappers.length; ++i) { - mapped[i] = mappers[i].boundedToUnbounded(point[i]); - } - - return mapped; - - } - - /** Compute the underlying function value from an unbounded point. - * <p> - * This method simply bounds the unbounded point using the mappings - * set up at construction and calls the underlying function using - * the bounded point. - * </p> - * @param point unbounded value - * @return underlying function value - * @see #unboundedToBounded(double[]) - */ - public double value(double[] point) { - return bounded.value(unboundedToBounded(point)); - } - - /** Mapping interface. */ - private interface Mapper { - - /** Map a value from unbounded to bounded. - * @param y unbounded value - * @return bounded value - */ - double unboundedToBounded(double y); - - /** Map a value from bounded to unbounded. - * @param x bounded value - * @return unbounded value - */ - double boundedToUnbounded(double x); - - } - - /** Local class for no bounds mapping. */ - private static class NoBoundsMapper implements Mapper { - - /** Simple constructor. - */ - public NoBoundsMapper() { - } - - /** {@inheritDoc} */ - public double unboundedToBounded(final double y) { - return y; - } - - /** {@inheritDoc} */ - public double boundedToUnbounded(final double x) { - return x; - } - - } - - /** Local class for lower bounds mapping. */ - private static class LowerBoundMapper implements Mapper { - - /** Low bound. */ - private final double lower; - - /** Simple constructor. - * @param lower lower bound - */ - public LowerBoundMapper(final double lower) { - this.lower = lower; - } - - /** {@inheritDoc} */ - public double unboundedToBounded(final double y) { - return lower + FastMath.exp(y); - } - - /** {@inheritDoc} */ - public double boundedToUnbounded(final double x) { - return FastMath.log(x - lower); - } - - } - - /** Local class for upper bounds mapping. */ - private static class UpperBoundMapper implements Mapper { - - /** Upper bound. */ - private final double upper; - - /** Simple constructor. - * @param upper upper bound - */ - public UpperBoundMapper(final double upper) { - this.upper = upper; - } - - /** {@inheritDoc} */ - public double unboundedToBounded(final double y) { - return upper - FastMath.exp(-y); - } - - /** {@inheritDoc} */ - public double boundedToUnbounded(final double x) { - return -FastMath.log(upper - x); - } - - } - - /** Local class for lower and bounds mapping. */ - private static class LowerUpperBoundMapper implements Mapper { - - /** Function from unbounded to bounded. */ - private final UnivariateFunction boundingFunction; - - /** Function from bounded to unbounded. */ - private final UnivariateFunction unboundingFunction; - - /** Simple constructor. - * @param lower lower bound - * @param upper upper bound - */ - public LowerUpperBoundMapper(final double lower, final double upper) { - boundingFunction = new Sigmoid(lower, upper); - unboundingFunction = new Logit(lower, upper); - } - - /** {@inheritDoc} */ - public double unboundedToBounded(final double y) { - return boundingFunction.value(y); - } - - /** {@inheritDoc} */ - public double boundedToUnbounded(final double x) { - return unboundingFunction.value(x); - } - - } - -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/b4669aad/src/main/java/org/apache/commons/math4/optimization/direct/MultivariateFunctionPenaltyAdapter.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/optimization/direct/MultivariateFunctionPenaltyAdapter.java b/src/main/java/org/apache/commons/math4/optimization/direct/MultivariateFunctionPenaltyAdapter.java deleted file mode 100644 index 113ebc8..0000000 --- a/src/main/java/org/apache/commons/math4/optimization/direct/MultivariateFunctionPenaltyAdapter.java +++ /dev/null @@ -1,190 +0,0 @@ -/* - * 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.math4.optimization.direct; - -import org.apache.commons.math4.analysis.MultivariateFunction; -import org.apache.commons.math4.exception.DimensionMismatchException; -import org.apache.commons.math4.exception.NumberIsTooSmallException; -import org.apache.commons.math4.util.FastMath; -import org.apache.commons.math4.util.MathUtils; - -/** - * <p>Adapter extending bounded {@link MultivariateFunction} to an unbouded - * domain using a penalty function.</p> - * - * <p> - * This adapter can be used to wrap functions subject to simple bounds on - * parameters so they can be used by optimizers that do <em>not</em> directly - * support simple bounds. - * </p> - * <p> - * The principle is that the user function that will be wrapped will see its - * parameters bounded as required, i.e when its {@code value} method is called - * with argument array {@code point}, the elements array will fulfill requirement - * {@code lower[i] <= point[i] <= upper[i]} for all i. Some of the components - * may be unbounded or bounded only on one side if the corresponding bound is - * set to an infinite value. The optimizer will not manage the user function by - * itself, but it will handle this adapter and it is this adapter that will take - * care the bounds are fulfilled. The adapter {@link #value(double[])} method will - * be called by the optimizer with unbound parameters, and the adapter will check - * if the parameters is within range or not. If it is in range, then the underlying - * user function will be called, and if it is not the value of a penalty function - * will be returned instead. - * </p> - * <p> - * This adapter is only a poor man solution to simple bounds optimization constraints - * that can be used with simple optimizers like {@link SimplexOptimizer} with {@link - * NelderMeadSimplex} or {@link MultiDirectionalSimplex}. A better solution is to use - * an optimizer that directly supports simple bounds like {@link CMAESOptimizer} or - * {@link BOBYQAOptimizer}. One caveat of this poor man solution is that if start point - * or start simplex is completely outside of the allowed range, only the penalty function - * is used, and the optimizer may converge without ever entering the range. - * </p> - * - * @see MultivariateFunctionMappingAdapter - * - * @deprecated As of 3.1 (to be removed in 4.0). - * @since 3.0 - */ - -@Deprecated -public class MultivariateFunctionPenaltyAdapter implements MultivariateFunction { - - /** Underlying bounded function. */ - private final MultivariateFunction bounded; - - /** Lower bounds. */ - private final double[] lower; - - /** Upper bounds. */ - private final double[] upper; - - /** Penalty offset. */ - private final double offset; - - /** Penalty scales. */ - private final double[] scale; - - /** Simple constructor. - * <p> - * When the optimizer provided points are out of range, the value of the - * penalty function will be used instead of the value of the underlying - * function. In order for this penalty to be effective in rejecting this - * point during the optimization process, the penalty function value should - * be defined with care. This value is computed as: - * <pre> - * penalty(point) = offset + ∑<sub>i</sub>[scale[i] * √|point[i]-boundary[i]|] - * </pre> - * where indices i correspond to all the components that violates their boundaries. - * </p> - * <p> - * So when attempting a function minimization, offset should be larger than - * the maximum expected value of the underlying function and scale components - * should all be positive. When attempting a function maximization, offset - * should be lesser than the minimum expected value of the underlying function - * and scale components should all be negative. - * minimization, and lesser than the minimum expected value of the underlying - * function when attempting maximization. - * </p> - * <p> - * These choices for the penalty function have two properties. First, all out - * of range points will return a function value that is worse than the value - * returned by any in range point. Second, the penalty is worse for large - * boundaries violation than for small violations, so the optimizer has an hint - * about the direction in which it should search for acceptable points. - * </p> - * @param bounded bounded function - * @param lower lower bounds for each element of the input parameters array - * (some elements may be set to {@code Double.NEGATIVE_INFINITY} for - * unbounded values) - * @param upper upper bounds for each element of the input parameters array - * (some elements may be set to {@code Double.POSITIVE_INFINITY} for - * unbounded values) - * @param offset base offset of the penalty function - * @param scale scale of the penalty function - * @exception DimensionMismatchException if lower bounds, upper bounds and - * scales are not consistent, either according to dimension or to bounadary - * values - */ - public MultivariateFunctionPenaltyAdapter(final MultivariateFunction bounded, - final double[] lower, final double[] upper, - final double offset, final double[] scale) { - - // safety checks - MathUtils.checkNotNull(lower); - MathUtils.checkNotNull(upper); - MathUtils.checkNotNull(scale); - if (lower.length != upper.length) { - throw new DimensionMismatchException(lower.length, upper.length); - } - if (lower.length != scale.length) { - throw new DimensionMismatchException(lower.length, scale.length); - } - for (int i = 0; i < lower.length; ++i) { - // note the following test is written in such a way it also fails for NaN - if (!(upper[i] >= lower[i])) { - throw new NumberIsTooSmallException(upper[i], lower[i], true); - } - } - - this.bounded = bounded; - this.lower = lower.clone(); - this.upper = upper.clone(); - this.offset = offset; - this.scale = scale.clone(); - - } - - /** Compute the underlying function value from an unbounded point. - * <p> - * This method simply returns the value of the underlying function - * if the unbounded point already fulfills the bounds, and compute - * a replacement value using the offset and scale if bounds are - * violated, without calling the function at all. - * </p> - * @param point unbounded point - * @return either underlying function value or penalty function value - */ - public double value(double[] point) { - - for (int i = 0; i < scale.length; ++i) { - if ((point[i] < lower[i]) || (point[i] > upper[i])) { - // bound violation starting at this component - double sum = 0; - for (int j = i; j < scale.length; ++j) { - final double overshoot; - if (point[j] < lower[j]) { - overshoot = scale[j] * (lower[j] - point[j]); - } else if (point[j] > upper[j]) { - overshoot = scale[j] * (point[j] - upper[j]); - } else { - overshoot = 0; - } - sum += FastMath.sqrt(overshoot); - } - return offset + sum; - } - } - - // all boundaries are fulfilled, we are in the expected - // domain of the underlying function - return bounded.value(point); - - } - -}
