You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by tn...@apache.org on 2015/02/25 22:49:41 UTC
[13/18] [math] Remove deprecated optimization package.
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);
-
- }
-
-}