diff --git a/src/main/java/org/apache/commons/math/optimization/direct/CMAESOptimizer.java b/src/main/java/org/apache/commons/math/optimization/direct/CMAESOptimizer.java deleted file mode 100644 index 3f2a764cf..000000000 --- a/src/main/java/org/apache/commons/math/optimization/direct/CMAESOptimizer.java +++ /dev/null @@ -1,1317 +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.math.optimization.direct; - -import java.util.ArrayList; -import java.util.Arrays; -import java.util.List; - -import org.apache.commons.math.analysis.MultivariateRealFunction; -import org.apache.commons.math.exception.MultiDimensionMismatchException; -import org.apache.commons.math.exception.NoDataException; -import org.apache.commons.math.exception.NotPositiveException; -import org.apache.commons.math.exception.OutOfRangeException; -import org.apache.commons.math.exception.TooManyEvaluationsException; -import org.apache.commons.math.linear.Array2DRowRealMatrix; -import org.apache.commons.math.linear.EigenDecomposition; -import org.apache.commons.math.linear.EigenDecompositionImpl; -import org.apache.commons.math.linear.MatrixUtils; -import org.apache.commons.math.linear.RealMatrix; -import org.apache.commons.math.optimization.GoalType; -import org.apache.commons.math.optimization.MultivariateRealOptimizer; -import org.apache.commons.math.optimization.RealPointValuePair; -import org.apache.commons.math.random.MersenneTwister; -import org.apache.commons.math.random.RandomGenerator; -import org.apache.commons.math.util.MathUtils; - -/** - * CMA-ES algorithm. This code is translated and adapted from the Matlab version - * of this algorithm as implemented in module {@code cmaes.m} version 3.51. - * - * Implements the active Covariance Matrix Adaptation Evolution Strategy (CMA-ES) - * for non-linear, non-convex, non-smooth, global function minimization. - * The CMA-Evolution Strategy (CMA-ES) is a reliable stochastic optimization method - * which should be applied, if derivative based methods, e.g. quasi-Newton BFGS or - * conjugate gradient, fail due to a rugged search landscape (e.g. noise, local - * optima, outlier, etc.) of the objective function. Like a - * quasi-Newton method the CMA-ES learns and applies a variable metric - * of the underlying search space. Unlike a quasi-Newton method the - * CMA-ES does neither estimate nor use gradients, making it considerably more - * reliable in terms of finding a good, or even close to optimal, solution, finally. - * - *

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 N=10 variables also the derivative-free simplex - * direct search method (Nelder and Mead) can be faster, but it is - * far less reliable than CMA-ES.

- * - *

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 30 N function - * evaluations. On difficult problems the complete - * optimization (a single run) is expected to take roughly between - * 30 N and 300 N2 - * function evaluations.

- * - * For more information, please refer to the following links: - * - * - * @version $Revision$ $Date$ - * @since 3.0 - */ - -public class CMAESOptimizer extends - BaseAbstractScalarOptimizer implements - MultivariateRealOptimizer { - - /** Default value for {@link #checkFeasableCount}: {@value}. */ - public static final int DEFAULT_CHECKFEASABLECOUNT = 0; - /** Default value for {@link #stopfitness}: {@value}. */ - public static final double DEFAULT_STOPFITNESS = 0; - /** Default value for {@link #isActiveCMA}: {@value}. */ - public static final boolean DEFAULT_ISACTIVECMA = true; - /** Default value for {@link #maxIterations}: {@value}. */ - public static final int DEFAULT_MAXITERATIONS = 30000; - /** Default value for {@link #diagonalOnly}: {@value}. */ - public static final int DEFAULT_DIAGONALONLY = 0; - /** Default value for {@link #random}. */ - public static final RandomGenerator DEFAULT_RANDOMGENERATOR = new MersenneTwister(); - - // global search parameters - /** - * Population size, offspring number. The primary strategy parameter to play - * with, which can be increased from its default value. Increasing the - * population size improves global search properties in exchange to speed. - * Speed decreases, as a rule, at most linearely with increasing population - * size. It is advisable to begin with the default small population size. - */ - private int lambda; // population size - /** - * Covariance update mechanism, default is active CMA. isActiveCMA = true - * turns on "active CMA" with a negative update of the covariance matrix and - * checks for positive definiteness. OPTS.CMA.active = 2 does not check for - * pos. def. and is numerically faster. Active CMA usually speeds up the - * adaptation. - */ - private boolean isActiveCMA; - /** - * Determines how often a new random offspring is generated in case it is - * not feasible / beyond the defined limits, default is 0. Only relevant if - * boundaries != null. - */ - private int checkFeasableCount; - /** - * Lower and upper boundaries of the objective variables. boundaries == null - * means no boundaries. - */ - private double[][] boundaries; - /** - * Individual sigma values - initial search volume. inputSigma determines - * the initial coordinate wise standard deviations for the search. Setting - * SIGMA one third of the initial search region is appropriate. - */ - private double[] inputSigma; - /** Number of objective variables/problem dimension */ - private int dimension; - /** - * Defines the number of initial iterations, where the covariance matrix - * remains diagonal and the algorithm has internally linear time complexity. - * diagonalOnly = 1 means keeping the covariance matrix always diagonal and - * this setting also exhibits linear space complexity. This can be - * particularly useful for dimension > 100. - * @see A Simple Modification in CMA-ES . - */ - 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 statisticsSigmaHistory = new ArrayList(); - /** History of mean matrix. */ - private List statisticsMeanHistory = new ArrayList(); - /** History of fitness values. */ - private List statisticsFitnessHistory = new ArrayList(); - /** History of D matrix. */ - private List statisticsDHistory = new ArrayList(); - - /** - * Default constructor, uses default parameters - */ - public CMAESOptimizer() { - this(0); - } - - /** - * @param lambda - * Population size. - */ - public CMAESOptimizer(int lambda) { - this(lambda, null, null, DEFAULT_MAXITERATIONS, DEFAULT_STOPFITNESS, - DEFAULT_ISACTIVECMA, DEFAULT_DIAGONALONLY, - DEFAULT_CHECKFEASABLECOUNT, DEFAULT_RANDOMGENERATOR, false); - } - - /** - * @param lambda - * Population size. - * @param inputSigma - * Initial search volume - sigma of offspring objective - * variables. - * @param boundaries - * Boundaries for objective variables. - */ - public CMAESOptimizer(int lambda, double[] inputSigma, - double[][] boundaries) { - this(lambda, inputSigma, boundaries, DEFAULT_MAXITERATIONS, DEFAULT_STOPFITNESS, - DEFAULT_ISACTIVECMA, DEFAULT_DIAGONALONLY, - DEFAULT_CHECKFEASABLECOUNT, DEFAULT_RANDOMGENERATOR, false); - } - - /** - * @param lambda - * Population size. - * @param inputSigma - * Initial search volume - sigma of offspring objective - * variables. - * @param boundaries - * Boundaries for objective variables. - * @param maxIterations - * Maximal number of iterations. - * @param stopfitness - * stop if objective function value < stopfitness. - * @param isActiveCMA - * Chooses the covariance matrix update method. - * @param diagonalOnly - * Number of initial iterations, where the covariance matrix - * remains diagonal. - * @param checkFeasableCount - * Determines how often new. random objective variables are - * generated in case they are out of bounds. - * @param random - * Used random generator. - * @param generateStatistics - * Indicates whether statistic data is collected. - */ - public CMAESOptimizer(int lambda, double[] inputSigma, - double[][] boundaries, int maxIterations, double stopfitness, - boolean isActiveCMA, int diagonalOnly, int checkFeasableCount, - RandomGenerator random, boolean generateStatistics) { - this.lambda = lambda; - this.inputSigma = inputSigma; - this.boundaries = boundaries; - this.maxIterations = maxIterations; - this.stopfitness = stopfitness; - this.isActiveCMA = isActiveCMA; - this.diagonalOnly = diagonalOnly; - this.checkFeasableCount = checkFeasableCount; - this.random = random; - this.generateStatistics = generateStatistics; - } - - /** - * @return History of sigma values. - */ - public List getStatisticsSigmaHistory() { - return statisticsSigmaHistory; - } - - /** - * @return History of mean matrix. - */ - public List getStatisticsMeanHistory() { - return statisticsMeanHistory; - } - - /** - * @return History of fitness values. - */ - public List getStatisticsFitnessHistory() { - return statisticsFitnessHistory; - } - - /** - * @return History of D matrix. - */ - public List getStatisticsDHistory() { - return statisticsDHistory; - } - - /** {@inheritDoc} */ - @Override - protected RealPointValuePair doOptimize() { - checkParameters(); - // -------------------- Initialization -------------------------------- - isMinimize = getGoalType().equals(GoalType.MINIMIZE); - final FitnessFunction fitfun = new FitnessFunction(boundaries, - isMinimize); - final double[] guess = fitfun.encode(getStartPoint()); - // number of objective variables/problem dimension - dimension = guess.length; - initializeCMA(guess); - iterations = 0; - double bestValue = fitfun.value(guess); - push(fitnessHistory, bestValue); - RealPointValuePair optimum = new RealPointValuePair(getStartPoint(), - isMinimize ? bestValue : -bestValue); - RealPointValuePair lastResult = null; - - // -------------------- Generation Loop -------------------------------- - - generationLoop: - for (iterations = 1; iterations <= maxIterations; iterations++) { - // Generate and evaluate lambda offspring - RealMatrix arz = randn1(dimension, lambda); - RealMatrix arx = zeros(dimension, lambda); - double[] fitness = new double[lambda]; - // generate random offspring - for (int k = 0; k < lambda; k++) { - RealMatrix arxk = null; - for (int i = 0; i < checkFeasableCount+1; i++) { - if (diagonalOnly <= 0) - arxk = xmean.add(BD.multiply(arz.getColumnMatrix(k)) - .scalarMultiply(sigma)); // m + sig * Normal(0,C) - else - arxk = xmean.add(times(diagD,arz.getColumnMatrix(k)) - .scalarMultiply(sigma)); - if (i >= checkFeasableCount || fitfun.isFeasible(arxk.getColumn(0))) - break; - // regenerate random arguments for row - arz.setColumn(k, randn(dimension)); - } - copyColumn(arxk, 0, arx, k); - try { - fitness[k] = fitfun.value(arx.getColumn(k)); // compute fitness - } catch (TooManyEvaluationsException e) { - break generationLoop; - } - } - // Sort by fitness and compute weighted mean into xmean - int[] arindex = sortedIndices(fitness); - // Calculate new xmean, this is selection and recombination - RealMatrix xold = xmean; // for speed up of Eq. (2) and (3) - RealMatrix bestArx = selectColumns(arx, MathUtils.copyOf(arindex, mu)); - xmean = bestArx.multiply(weights); - RealMatrix bestArz = selectColumns(arz, MathUtils.copyOf(arindex, mu)); - RealMatrix zmean = bestArz.multiply(weights); - boolean hsig = updateEvolutionPaths(zmean, xold); - if (diagonalOnly <= 0) - updateCovariance(hsig, bestArx, arz, arindex, xold); - else - updateCovarianceDiagonalOnly(hsig, bestArz, xold); - // Adapt step size sigma - Eq. (5) - sigma *= Math.exp(Math.min(1.0,(normps/chiN - 1.)*cs/damps)); - double bestFitness = fitness[arindex[0]]; - double worstFitness = fitness[arindex[arindex.length-1]]; - if (bestValue > bestFitness) { - bestValue = bestFitness; - lastResult = optimum; - optimum = new RealPointValuePair( - fitfun.decode(bestArx.getColumn(0)), - isMinimize ? bestFitness : -bestFitness); - if (getConvergenceChecker() != null && lastResult != null) { - if (getConvergenceChecker().converged( - iterations, optimum, lastResult)) - break generationLoop; - } - } - // handle termination criteria - // Break, if fitness is good enough - if (stopfitness != 0) { // only if stopfitness is defined - if (bestFitness < (isMinimize ? stopfitness : -stopfitness)) - break generationLoop; - } - double[] sqrtDiagC = sqrt(diagC).getColumn(0); - double[] pcCol = pc.getColumn(0); - for (int i = 0; i < dimension; i++) { - if (sigma*(Math.max(Math.abs(pcCol[i]), sqrtDiagC[i])) > stopTolX) - break; - if (i >= dimension-1) - break generationLoop; - } - for (int i = 0; i < dimension; i++) - if (sigma*sqrtDiagC[i] > stopTolUpX) - break generationLoop; - double historyBest = min(fitnessHistory); - double historyWorst = max(fitnessHistory); - if (iterations > 2 && Math.max(historyWorst, worstFitness) - - Math.min(historyBest, bestFitness) < stopTolFun) - break generationLoop; - if (iterations > fitnessHistory.length && - historyWorst-historyBest < stopTolHistFun) - break generationLoop; - // condition number of the covariance matrix exceeds 1e14 - if (max(diagD)/min(diagD) > 1e7) - break generationLoop; - // user defined termination - if (getConvergenceChecker() != null) { - RealPointValuePair current = - new RealPointValuePair(bestArx.getColumn(0), - isMinimize ? bestFitness : -bestFitness); - if (lastResult != null && - getConvergenceChecker().converged( - iterations, current, lastResult)) - break generationLoop; - lastResult = current; - } - // Adjust step size in case of equal function values (flat fitness) - if (bestValue == fitness[arindex[(int)(0.1+lambda/4.)]]) - sigma = sigma * Math.exp(0.2+cs/damps); - if (iterations > 2 && Math.max(historyWorst, bestFitness) - - Math.min(historyBest, bestFitness) == 0) - sigma = sigma * Math.exp(0.2+cs/damps); - // store best in history - push(fitnessHistory,bestFitness); - fitfun.setValueRange(worstFitness-bestFitness); - if (generateStatistics) { - statisticsSigmaHistory.add(sigma); - statisticsFitnessHistory.add(bestFitness); - statisticsMeanHistory.add(xmean.transpose()); - statisticsDHistory.add(diagD.transpose().scalarMultiply(1E5)); - } - } - return optimum; - } - - /** - * Checks dimensions and values of boundaries and inputSigma if defined. - */ - private void checkParameters() { - double[] init = getStartPoint(); - if (boundaries != null) { - if (boundaries.length != 2) - throw new MultiDimensionMismatchException( - new Integer[] { boundaries.length }, - new Integer[] { 2 }); - if (boundaries[0] == null || boundaries[1] == null) - throw new NoDataException(); - if (boundaries[0].length != init.length) - throw new MultiDimensionMismatchException( - new Integer[] { boundaries[0].length }, - new Integer[] { init.length }); - if (boundaries[1].length != init.length) - throw new MultiDimensionMismatchException( - new Integer[] { boundaries[1].length }, - new Integer[] { init.length }); - for (int i = 0; i < init.length; i++) { - if (boundaries[0][i] > init[i] || boundaries[1][i] < init[i]) - throw new OutOfRangeException(init[i], boundaries[0][i], - boundaries[1][i]); - } - } - if (inputSigma != null) { - if (inputSigma.length != init.length) - throw new MultiDimensionMismatchException( - new Integer[] { inputSigma.length }, - new Integer[] { init.length }); - for (int i = 0; i < init.length; i++) { - if (inputSigma[i] < 0) - throw new NotPositiveException(inputSigma[i]); - if (boundaries != null) { - if (inputSigma[i] > 1.0) - throw new OutOfRangeException(inputSigma[i], 0, 1.0); - } - } - } - } - - /** - * Initialization of the dynamic search parameters - * - * @param guess - * initial guess for the arguments of the fitness function - */ - - private void initializeCMA(double[] guess) { - if (lambda <= 0) - lambda = 4 + (int) (3. * Math.log(dimension)); - // initialize sigma - double[][] sigmaArray = new double[guess.length][1]; - for (int i = 0; i < guess.length; i++) - sigmaArray[i][0] = inputSigma != null ? inputSigma[i] : 0.3; - RealMatrix insigma = new Array2DRowRealMatrix(sigmaArray, false); - sigma = max(insigma); // overall standard deviation - - // initialize termination criteria - stopTolUpX = 1e3 * max(insigma); - stopTolX = 1e-11 * max(insigma); - stopTolFun = 1e-12; - stopTolHistFun = 1e-13; - - // initialize selection strategy parameters - mu = lambda / 2; // number of parents/points for recombination - logMu2 = Math.log(mu + 0.5); - weights = log(sequence(1, mu, 1)).scalarMultiply(-1.).scalarAdd(logMu2); - double sumw = 0; - double sumwq = 0; - for (int i = 0; i < mu; i++) { - double w = weights.getEntry(i, 0); - sumw += w; - sumwq += w * w; - } - weights = weights.scalarMultiply(1. / sumw); - mueff = sumw * sumw / sumwq; // variance-effectiveness of sum w_i x_i - - // initialize dynamic strategy parameters and constants - cc = (4. + mueff / dimension) / - (dimension + 4. + 2. * mueff / dimension); - cs = (mueff + 2.) / (dimension + mueff + 3.); - damps = (1. + 2. * Math.max(0, Math.sqrt((mueff - 1.) / - (dimension + 1.)) - 1.)) * - Math.max(0.3, 1. - dimension / - (1e-6 + Math.min(maxIterations, getMaxEvaluations() / - lambda))) + cs; // minor increment - ccov1 = 2. / ((dimension + 1.3) * (dimension + 1.3) + mueff); - ccovmu = Math.min(1 - ccov1, 2. * (mueff - 2. + 1. / mueff) / - ((dimension + 2.) * (dimension + 2.) + mueff)); - ccov1Sep = Math.min(1, ccov1 * (dimension + 1.5) / 3.); - ccovmuSep = Math.min(1 - ccov1, ccovmu * (dimension + 1.5) / 3.); - chiN = Math.sqrt(dimension) * - (1. - 1. / (4. * dimension) + 1 / (21. * dimension * dimension)); - // intialize CMA internal values - updated each generation - xmean = MatrixUtils.createColumnRealMatrix(guess); // objective - // variables - diagD = insigma.scalarMultiply(1. / sigma); - diagC = square(diagD); - pc = zeros(dimension, 1); // evolution paths for C and sigma - ps = zeros(dimension, 1); // B defines the coordinate system - normps = norm(ps); - - B = eye(dimension, dimension); - D = ones(dimension, 1); // diagonal D defines the scaling - BD = times(B, repmat(diagD.transpose(), dimension, 1)); - C = B.multiply(diag(square(D)).multiply(B.transpose())); // covariance - historySize = 10 + (int) (3. * 10. * dimension / lambda); - fitnessHistory = new double[historySize]; // history of fitness values - for (int i = 0; i < historySize; i++) - fitnessHistory[i] = Double.MAX_VALUE; - } - - /** - * Update of the evolution paths ps and pc - * - * @param zmean - * weighted row matrix of the gaussian random numbers generating - * the current offspring - * @param xold - * xmean matrix of the previous generation - * @return hsig flag indicating a small correction - */ - private boolean updateEvolutionPaths(RealMatrix zmean, RealMatrix xold) { - ps = ps.scalarMultiply(1. - cs).add( - B.multiply(zmean).scalarMultiply( - Math.sqrt(cs * (2. - cs) * mueff))); - normps = norm(ps); - boolean hsig = normps / - Math.sqrt(1. - Math.pow(1. - cs, 2. * iterations)) / - chiN < 1.4 + 2. / (dimension + 1.); - pc = pc.scalarMultiply(1. - cc); - if (hsig) - pc = pc.add(xmean.subtract(xold).scalarMultiply( - Math.sqrt(cc * (2. - cc) * mueff) / sigma)); - return hsig; - } - - /** - * Update of the covariance matrix C for diagonalOnly > 0 - * - * @param hsig - * flag indicating a small correction - * @param bestArz - * fitness-sorted matrix of the gaussian random values of the - * current offspring - * @param xold - * xmean matrix of the previous generation - */ - private void updateCovarianceDiagonalOnly(boolean hsig, final RealMatrix bestArz, - final RealMatrix xold) { - // minor correction if hsig==false - double oldFac = hsig ? 0 : ccov1Sep * cc * (2. - cc); - oldFac += 1. - ccov1Sep - ccovmuSep; - diagC = diagC.scalarMultiply(oldFac) // regard old matrix - // plus rank one update - .add(square(pc).scalarMultiply(ccov1Sep)) - // plus rank mu update - .add((times(diagC, square(bestArz).multiply(weights))) - .scalarMultiply(ccovmuSep)); - diagD = sqrt(diagC); // replaces eig(C) - if (diagonalOnly > 1 && iterations > diagonalOnly) { - // full covariance matrix from now on - diagonalOnly = 0; - B = eye(dimension, dimension); - BD = diag(diagD); - C = diag(diagC); - } - } - - /** - * Update of the covariance matrix C - * - * @param hsig - * flag indicating a small correction - * @param bestArx - * fitness-sorted matrix of the argument vectors producing the - * current offspring - * @param arz - * unsorted matrix containing the gaussian random values of the - * current offspring - * @param arindex - * indices indicating the fitness-order of the current offspring - * @param xold - * xmean matrix of the previous generation - */ - private void updateCovariance(boolean hsig, final RealMatrix bestArx, - final RealMatrix arz, final int[] arindex, final RealMatrix xold) { - double negccov = 0; - if (ccov1 + ccovmu > 0) { - RealMatrix arpos = bestArx.subtract(repmat(xold, 1, mu)) - .scalarMultiply(1. / sigma); // mu difference vectors - RealMatrix roneu = pc.multiply(pc.transpose()) - .scalarMultiply(ccov1); // rank one update - // minor correction if hsig==false - double oldFac = hsig ? 0 : ccov1 * cc * (2. - cc); - oldFac += 1. - ccov1 - ccovmu; - if (isActiveCMA) { - // Adapt covariance matrix C active CMA - negccov = (1. - ccovmu) * 0.25 * mueff / - (Math.pow(dimension + 2., 1.5) + 2. * mueff); - double negminresidualvariance = 0.66; - // keep at least 0.66 in all directions, small popsize are most - // critical - double negalphaold = 0.5; // where to make up for the variance - // loss, - // prepare vectors, compute negative updating matrix Cneg - int[] arReverseIndex = reverse(arindex); - RealMatrix arzneg - = selectColumns(arz, MathUtils.copyOf(arReverseIndex, mu)); - RealMatrix arnorms = sqrt(sumRows(square(arzneg))); - int[] idxnorms = sortedIndices(arnorms.getRow(0)); - RealMatrix arnormsSorted = selectColumns(arnorms, idxnorms); - int[] idxReverse = reverse(idxnorms); - RealMatrix arnormsReverse = selectColumns(arnorms, idxReverse); - arnorms = divide(arnormsReverse, arnormsSorted); - int[] idxInv = inverse(idxnorms); - RealMatrix arnormsInv = selectColumns(arnorms, idxInv); - // check and set learning rate negccov - double negcovMax = (1. - negminresidualvariance) / - square(arnormsInv).multiply(weights).getEntry(0, 0); - if (negccov > negcovMax) - negccov = negcovMax; - arzneg = times(arzneg, repmat(arnormsInv, dimension, 1)); - RealMatrix artmp = BD.multiply(arzneg); - RealMatrix Cneg = artmp.multiply(diag(weights)).multiply( - artmp.transpose()); - oldFac += negalphaold * negccov; - C = C.scalarMultiply(oldFac) - // regard old matrix - .add(roneu) - // plus rank one update - .add(arpos.scalarMultiply( - // plus rank mu update - ccovmu + (1. - negalphaold) * negccov) - .multiply( - times(repmat(weights, 1, dimension), - arpos.transpose()))) - .subtract(Cneg.scalarMultiply(negccov)); - } else { - // Adapt covariance matrix C - nonactive - C = C.scalarMultiply(oldFac) // regard old matrix - .add(roneu) - // plus rank one update - .add(arpos.scalarMultiply(ccovmu) // plus rank mu update - .multiply( - times(repmat(weights, 1, dimension), - arpos.transpose()))); - } - } - updateBD(negccov); - } - - /** - * Update B and D from C - * - * @param negccov - * Negative covariance factor. - */ - private void updateBD(double negccov) { - if (ccov1 + ccovmu + negccov > 0 && - (iterations % 1. / (ccov1 + ccovmu + negccov) / dimension / 10.) < 1.) { - // to achieve O(N^2) - C = triu(C, 0).add(triu(C, 1).transpose()); - // enforce symmetry to prevent complex numbers - EigenDecomposition eig = new EigenDecompositionImpl(C, 1.0); - B = eig.getV(); // eigen decomposition, B==normalized eigenvectors - D = eig.getD(); - diagD = diag(D); - if (min(diagD) <= 0) { - for (int i = 0; i < dimension; i++) - if (diagD.getEntry(i, 0) < 0) - diagD.setEntry(i, 0, 0.); - double tfac = max(diagD) / 1e14; - C = C.add(eye(dimension, dimension).scalarMultiply(tfac)); - diagD = diagD.add(ones(dimension, 1).scalarMultiply(tfac)); - } - if (max(diagD) > 1e14 * min(diagD)) { - double tfac = max(diagD) / 1e14 - min(diagD); - C = C.add(eye(dimension, dimension).scalarMultiply(tfac)); - diagD = diagD.add(ones(dimension, 1).scalarMultiply(tfac)); - } - diagC = diag(C); - diagD = sqrt(diagD); // D contains standard deviations now - BD = times(B, repmat(diagD.transpose(), dimension, 1)); // O(n^2) - } - } - - /** - * Pushes the current best fitness value in a history queue. - * - * @param vals - * the history queue - * @param val - * current best fitness value - */ - private static void push(double[] vals, double val) { - for (int i = vals.length-1; i > 0; i--) - vals[i] = vals[i-1]; - vals[0] = val; - } - - /** - * Sorts fitness values. - * - * @param doubles - * array of values to be sorted - * @return sorted array of indices pointing into doubles - */ - private int[] sortedIndices(final double[] doubles) { - DoubleIndex[] dis = new DoubleIndex[doubles.length]; - for (int i = 0; i < doubles.length; i++) - dis[i] = new DoubleIndex(doubles[i], i); - Arrays.sort(dis); - int[] indices = new int[doubles.length]; - for (int i = 0; i < doubles.length; i++) - indices[i] = dis[i].index; - return indices; - } - - /** - * Used to sort fitness values. Sorting is always in lower value first - * order. - */ - private static class DoubleIndex implements Comparable { - - /** Value to compare. */ - private double value; - /** Index into sorted array. */ - private int index; - - /** - * @param value - * Value to compare. - * @param index - * Index into sorted array. - */ - DoubleIndex(double value, int index) { - this.value = value; - this.index = index; - } - - /** {@inheritDoc} */ - public int compareTo(DoubleIndex o) { - return Double.compare(value, o.value); - } - } - - /** - * Normalizes fitness values to the range [0,1]. Adds a penalty to the - * fitness value if out of range. The penalty is adjusted by calling - * setValueRange(). - */ - private class FitnessFunction { - - /** Optional bounds for the objective variables */ - private double[][] boundaries; - /** Determines the penalty for boundary violations */ - private double valueRange = 1.0; - /** - * Flag indicating whether the objective variables are forced into their - * bounds if defined - */ - private boolean isRepairMode = true; - /** Flag indicating the optimization goal. */ - private boolean isMinimize = true; - - /** - * @param boundaries - * Bounds for the objective variables. - * @param isMinimize - * Flag indicating the optimization goal. - */ - private FitnessFunction(final double[][] boundaries, - final boolean isMinimize) { - this.boundaries = boundaries; - this.isMinimize = isMinimize; - } - - /** - * @param x - * Original objective variables. - * @return Normalized objective variables. - */ - private double[] encode(final double[] x) { - if (boundaries == null) - return x; - double[] res = new double[x.length]; - for (int i = 0; i < x.length; i++) { - double diff = boundaries[1][i] - boundaries[0][i]; - res[i] = (x[i] - boundaries[0][i]) / diff; - } - return res; - } - - /** - * @param x - * Normalized objective variables. - * @return Original objective variables. - */ - private double[] decode(final double[] x) { - if (boundaries == null) - return x; - double[] res = new double[x.length]; - for (int i = 0; i < x.length; i++) { - double diff = boundaries[1][i] - boundaries[0][i]; - res[i] = diff * x[i] + boundaries[0][i]; - } - return res; - } - - /** - * @param point - * Normalized objective variables. - * @return Objective value + penalty for violated bounds. - */ - private double value(final double[] point) { - double value; - if (boundaries != null && isRepairMode) { - double[] repaired = repair(point); - value = CMAESOptimizer.this - .computeObjectiveValue(decode(repaired)) + - penalty(point, repaired); - } else - value = CMAESOptimizer.this - .computeObjectiveValue(decode(point)); - return isMinimize ? value : -value; - } - - /** - * @param x - * Normalized objective variables. - * @return True if in bounds - */ - private boolean isFeasible(final double[] x) { - if (boundaries == null) - return true; - for (int i = 0; i < x.length; i++) { - if (x[i] < 0) - return false; - if (x[i] > 1.0) - return false; - } - return true; - } - - /** - * @param valueRange - * Adjusts the penalty computation. - */ - private void setValueRange(double valueRange) { - this.valueRange = valueRange; - } - - /** - * @param x - * Normalized objective variables. - * @return Repaired objective variables - all in bounds. - */ - private double[] repair(final double[] x) { - double[] repaired = new double[x.length]; - for (int i = 0; i < x.length; i++) { - if (x[i] < 0) - repaired[i] = 0; - else if (x[i] > 1.0) - repaired[i] = 1.0; - else - repaired[i] = x[i]; - } - return repaired; - } - - /** - * @param x - * Normalized objective variables. - * @param repaired - * Repaired objective variables. - * @return Penalty value according to the violation of the bounds. - */ - private double penalty(final double[] x, final double[] repaired) { - double penalty = 0; - for (int i = 0; i < x.length; i++) { - double diff = Math.abs(x[i] - repaired[i]); - penalty += diff * valueRange; - } - return isMinimize ? penalty : -penalty; - } - } - - // -----Matrix utility functions similar to the Matlab build in functions------ - - /** - * @param m - * Input matrix - * @return Matrix representing the element wise logarithm of m. - */ - private static RealMatrix log(final RealMatrix m) { - double[][] d = new double[m.getRowDimension()][m.getColumnDimension()]; - for (int r = 0; r < m.getRowDimension(); r++) - for (int c = 0; c < m.getColumnDimension(); c++) - d[r][c] = Math.log(m.getEntry(r, c)); - return new Array2DRowRealMatrix(d, false); - } - - /** - * @param m - * Input matrix - * @return Matrix representing the element wise square root of m. - */ - private static RealMatrix sqrt(final RealMatrix m) { - double[][] d = new double[m.getRowDimension()][m.getColumnDimension()]; - for (int r = 0; r < m.getRowDimension(); r++) - for (int c = 0; c < m.getColumnDimension(); c++) - d[r][c] = Math.sqrt(m.getEntry(r, c)); - return new Array2DRowRealMatrix(d, false); - } - - /** - * @param m Input matrix - * @return Matrix representing the element wise square (^2) of m. - */ - private static RealMatrix square(final RealMatrix m) { - double[][] d = new double[m.getRowDimension()][m.getColumnDimension()]; - for (int r = 0; r < m.getRowDimension(); r++) - for (int c = 0; c < m.getColumnDimension(); c++) { - double e = m.getEntry(r, c); - d[r][c] = e * e; - } - return new Array2DRowRealMatrix(d, false); - } - - /** - * @param m - * Input matrix 1. - * @param n - * Input matrix 2. - * @return Matrix where the elements of m and m are element wise multiplied. - */ - private static RealMatrix times(final RealMatrix m, final RealMatrix n) { - double[][] d = new double[m.getRowDimension()][m.getColumnDimension()]; - for (int r = 0; r < m.getRowDimension(); r++) - for (int c = 0; c < m.getColumnDimension(); c++) - d[r][c] = m.getEntry(r, c)*n.getEntry(r, c); - return new Array2DRowRealMatrix(d, false); - } - - /** - * @param m - * Input matrix 1. - * @param n - * Input matrix 2. - * @return Matrix where the elements of m and m are element wise divided. - */ - private static RealMatrix divide(final RealMatrix m, final RealMatrix n) { - double[][] d = new double[m.getRowDimension()][m.getColumnDimension()]; - for (int r = 0; r < m.getRowDimension(); r++) - for (int c = 0; c < m.getColumnDimension(); c++) - d[r][c] = m.getEntry(r, c)/n.getEntry(r, c); - return new Array2DRowRealMatrix(d, false); - } - - /** - * @param m Input matrix. - * @param cols Columns to select. - * @return Matrix representing the selected columns. - */ - private static RealMatrix selectColumns(final RealMatrix m, final int[] cols) { - double[][] d = new double[m.getRowDimension()][cols.length]; - for (int r = 0; r < m.getRowDimension(); r++) - for (int c = 0; c < cols.length; c++) - d[r][c] = m.getEntry(r, cols[c]); - return new Array2DRowRealMatrix(d, false); - } - - /** - * @param m Input matrix. - * @param k diagonal position. - * @return Upper triangular part of matrix. - */ - private static RealMatrix triu(final RealMatrix m, int k) { - double[][] d = new double[m.getRowDimension()][m.getColumnDimension()]; - for (int r = 0; r < m.getRowDimension(); r++) - for (int c = 0; c < m.getColumnDimension(); c++) - d[r][c] = r <= c - k ? m.getEntry(r, c) : 0; - return new Array2DRowRealMatrix(d, false); - } - - /** - * @param m - * Input matrix. - * @return Norm of the matrix. - */ - private static double norm(final RealMatrix m) { - double sum = 0; - for (int r = 0; r < m.getRowDimension(); r++) - for (int c = 0; c < m.getColumnDimension(); c++) { - double e = m.getEntry(r, c); - sum += e*e; - } - return Math.sqrt(sum); - } - - /** - * @param m - * Input matrix. - * @return Row matrix representing the sums of the rows. - */ - private static RealMatrix sumRows(final RealMatrix m) { - double[][] d = new double[1][m.getColumnDimension()]; - for (int c = 0; c < m.getColumnDimension(); c++) { - double sum = 0; - for (int r = 0; r < m.getRowDimension(); r++) - sum += m.getEntry(r, c); - d[0][c] = sum; - } - return new Array2DRowRealMatrix(d, false); - } - - /** - * @param m - * Input matrix. - * @return Diagonal n X n matrix if m is a column matrix, Rolumn matrix - * representing the diagonal if m is a nXn matrix. - */ - private static RealMatrix diag(final RealMatrix m) { - if (m.getColumnDimension() == 1) { - double[][] d = new double[m.getRowDimension()][m.getRowDimension()]; - for (int i = 0; i < m.getRowDimension(); i++) - d[i][i] = m.getEntry(i, 0); - return new Array2DRowRealMatrix(d, false); - } else { - double[][] d = new double[m.getRowDimension()][1]; - for (int i = 0; i < m.getColumnDimension(); i++) - d[i][0] = m.getEntry(i, i); - return new Array2DRowRealMatrix(d, false); - } - } - - /** - * Copies a row from m1 to m2. - * - * @param m1 - * Source matrix 1. - * @param col1 - * Source column. - * @param m2 - * Target matrix. - * @param col2 - * Target column. - */ - private static void copyColumn(final RealMatrix m1, int col1, RealMatrix m2, int col2) { - for (int i = 0; i < m1.getRowDimension(); i++) - m2.setEntry(i, col2, m1.getEntry(i, col1)); - } - - /** - * @param n - * Number of rows. - * @param m - * Number of columns. - * @return n X m matrix of 1.0-values. - */ - private static RealMatrix ones(int n, int m) { - double[][] d = new double[n][m]; - for (int r = 0; r < n; r++) - Arrays.fill(d[r], 1.0); - return new Array2DRowRealMatrix(d, false); - } - - /** - * @param n - * Number of rows. - * @param m - * Number of columns. - * @return n X m matrix of 0.0-values, diagonal has values 1.0. - */ - private static RealMatrix eye(int n, int m) { - double[][] d = new double[n][m]; - for (int r = 0; r < n; r++) - if (r < m) - d[r][r] = 1; - return new Array2DRowRealMatrix(d, false); - } - - /** - * @param n - * Number of rows. - * @param m - * Number of columns. - * @return n X m matrix of 0.0-values. - */ - private static RealMatrix zeros(int n, int m) { - return new Array2DRowRealMatrix(n, m); - } - - /** - * @param mat - * Input matrix. - * @param n - * Number of row replicates. - * @param m - * Number of column replicates. - * @return Matrix which replicates the input matrix in both directions. - */ - private static RealMatrix repmat(final RealMatrix mat, int n, int m) { - int rd = mat.getRowDimension(); - int cd = mat.getColumnDimension(); - double[][] d = new double[n * rd][m * cd]; - for (int r = 0; r < n * rd; r++) - for (int c = 0; c < m * cd; c++) - d[r][c] = mat.getEntry(r % rd, c % cd); - return new Array2DRowRealMatrix(d, false); - } - - /** - * @param start - * Start value. - * @param end - * End value. - * @param step - * Step size. - * @return Sequence as column matrix. - */ - private static RealMatrix sequence(double start, double end, double step) { - int size = (int) ((end - start) / step + 1); - double[][] d = new double[size][1]; - double value = start; - for (int r = 0; r < size; r++) { - d[r][0] = value; - value += step; - } - return new Array2DRowRealMatrix(d, false); - } - - /** - * @param m - * Input matrix. - * @return Maximum of matrix element values. - */ - private static double max(final RealMatrix m) { - double max = -Double.MAX_VALUE; - for (int r = 0; r < m.getRowDimension(); r++) - for (int c = 0; c < m.getColumnDimension(); c++) { - double e = m.getEntry(r, c); - if (max < e) - max = e; - } - return max; - } - - /** - * @param m - * Input matrix. - * @return Minimum of matrix element values. - */ - private static double min(final RealMatrix m) { - double min = Double.MAX_VALUE; - for (int r = 0; r < m.getRowDimension(); r++) - for (int c = 0; c < m.getColumnDimension(); c++) { - double e = m.getEntry(r, c); - if (min > e) - min = e; - } - return min; - } - - /** - * @param m - * Input array. - * @return Maximum of array values. - */ - private static double max(final double[] m) { - double max = -Double.MAX_VALUE; - for (int r = 0; r < m.length; r++) - if (max < m[r]) - max = m[r]; - return max; - } - - /** - * @param m - * Input array. - * @return Minimum of array values. - */ - private static double min(final double[] m) { - double min = Double.MAX_VALUE; - for (int r = 0; r < m.length; r++) - if (min > m[r]) - min = m[r]; - return min; - } - - /** - * @param indices - * Input index array. - * @return Inverse of the mapping defined by indices - */ - private static int[] inverse(final int[] indices) { - int[] inverse = new int[indices.length]; - for (int i = 0; i < indices.length; i++) - inverse[indices[i]] = i; - return inverse; - } - - /** - * @param indices - * Input index array. - * @return Indices in inverse order (last is first) - */ - private static int[] reverse(final int[] indices) { - int[] reverse = new int[indices.length]; - for (int i = 0; i < indices.length; i++) - reverse[i] = indices[indices.length - i - 1]; - return reverse; - } - - /** - * @param size - * Length of random array. - * @return Array of gaussian random numbers. - */ - private double[] randn(int size) { - double[] randn = new double[size]; - for (int i = 0; i < size; i++) - randn[i] = random.nextGaussian(); - return randn; - } - - /** - * @param size - * Number of rows. - * @param popSize - * Population size. - * @return 2-dimensional matrix of gaussian random numbers. - */ - private RealMatrix randn1(int size, int popSize) { - double[][] d = new double[size][popSize]; - for (int r = 0; r < size; r++) - for (int c = 0; c < popSize; c++) - d[r][c] = random.nextGaussian(); - return new Array2DRowRealMatrix(d, false); - } -} - diff --git a/src/main/java/org/apache/commons/math/util/MathUtils.java b/src/main/java/org/apache/commons/math/util/MathUtils.java index 0a56b9945..59e4b1c79 100644 --- a/src/main/java/org/apache/commons/math/util/MathUtils.java +++ b/src/main/java/org/apache/commons/math/util/MathUtils.java @@ -2220,7 +2220,10 @@ public final class MathUtils { * @return the copied array. */ public static int[] copyOf(int[] source) { - return copyOf(source, source.length); + final int len = source.length; + final int[] output = new int[len]; + System.arraycopy(source, 0, output, 0, len); + return output; } /** @@ -2230,36 +2233,9 @@ public final class MathUtils { * @return the copied array. */ public static double[] copyOf(double[] source) { - return copyOf(source, source.length); - } - - /** - * Creates a copy of the {@code source} array. - * - * @param source Array to be copied. - * @param len Number of entries to copy. If smaller then the source - * length, the copy will be truncated, if larger it will padded with - * zeroes. - * @return the copied array. - */ - public static int[] copyOf(int[] source, int len) { - final int[] output = new int[len]; - System.arraycopy(source, 0, output, 0, FastMath.min(len, source.length)); - return output; - } - - /** - * Creates a copy of the {@code source} array. - * - * @param source Array to be copied. - * @param len Number of entries to copy. If smaller then the source - * length, the copy will be truncated, if larger it will padded with - * zeroes. - * @return the copied array. - */ - public static double[] copyOf(double[] source, int len) { + final int len = source.length; final double[] output = new double[len]; - System.arraycopy(source, 0, output, 0, FastMath.min(len, source.length)); + System.arraycopy(source, 0, output, 0, len); return output; } } diff --git a/src/test/java/org/apache/commons/math/optimization/direct/CMAESOptimizerTest.java b/src/test/java/org/apache/commons/math/optimization/direct/CMAESOptimizerTest.java deleted file mode 100644 index f493fdc32..000000000 --- a/src/test/java/org/apache/commons/math/optimization/direct/CMAESOptimizerTest.java +++ /dev/null @@ -1,667 +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.math.optimization.direct; - -import java.util.Arrays; -import java.util.Random; - -import org.apache.commons.math.MathException; -import org.apache.commons.math.analysis.MultivariateRealFunction; -import org.apache.commons.math.exception.MathUserException; -import org.apache.commons.math.exception.MultiDimensionMismatchException; -import org.apache.commons.math.exception.NoDataException; -import org.apache.commons.math.exception.NotPositiveException; -import org.apache.commons.math.exception.OutOfRangeException; -import org.apache.commons.math.optimization.GoalType; -import org.apache.commons.math.optimization.MultivariateRealOptimizer; -import org.apache.commons.math.optimization.RealPointValuePair; -import org.apache.commons.math.random.MersenneTwister; -import org.junit.Assert; -import org.junit.Test; - -/** - * Test for {@link CMAESOptimizer}. - */ -public class CMAESOptimizerTest { - - static final int DIM = 13; - static final int LAMBDA = 4 + (int)(3.*Math.log(DIM)); - - @Test(expected = OutOfRangeException.class) - public void testInitOutofbounds() throws MathUserException, MathException { - double[] startPoint = point(DIM,3); - double[] insigma = null; - double[][] boundaries = boundaries(DIM,-1,2); - RealPointValuePair expected = - new RealPointValuePair(point(DIM,1.0),0.0); - doTest(new Rosen(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13, - 1e-13, 1e-6, 100000, expected); - } - - @Test(expected = MultiDimensionMismatchException.class) - public void testBoundariesDimensionMismatch() throws MathUserException, MathException { - double[] startPoint = point(DIM,0.5); - double[] insigma = null; - double[][] boundaries = boundaries(DIM+1,-1,2); - RealPointValuePair expected = - new RealPointValuePair(point(DIM,1.0),0.0); - doTest(new Rosen(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13, - 1e-13, 1e-6, 100000, expected); - } - - @Test(expected = NoDataException.class) - public void testBoundariesNoData() throws MathUserException, MathException { - double[] startPoint = point(DIM,0.5); - double[] insigma = null; - double[][] boundaries = boundaries(DIM,-1,2); - boundaries[1] = null; - RealPointValuePair expected = - new RealPointValuePair(point(DIM,1.0),0.0); - doTest(new Rosen(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13, - 1e-13, 1e-6, 100000, expected); - } - - @Test(expected = NotPositiveException.class) - public void testInputSigmaNegative() throws MathUserException, MathException { - double[] startPoint = point(DIM,0.5); - double[] insigma = point(DIM,-0.5); - double[][] boundaries = null; - RealPointValuePair expected = - new RealPointValuePair(point(DIM,1.0),0.0); - doTest(new Rosen(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13, - 1e-13, 1e-6, 100000, expected); - } - - @Test(expected = OutOfRangeException.class) - public void testInputSigmaOutOfRange() throws MathUserException, MathException { - double[] startPoint = point(DIM,0.5); - double[] insigma = point(DIM, 1.1); - double[][] boundaries = boundaries(DIM,-1,2); - RealPointValuePair expected = - new RealPointValuePair(point(DIM,1.0),0.0); - doTest(new Rosen(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13, - 1e-13, 1e-6, 100000, expected); - } - - @Test(expected = MultiDimensionMismatchException.class) - public void testInputSigmaDimensionMismatch() throws MathUserException, MathException { - double[] startPoint = point(DIM,0.5); - double[] insigma = point(DIM+1,-0.5); - double[][] boundaries = null;; - RealPointValuePair expected = - new RealPointValuePair(point(DIM,1.0),0.0); - doTest(new Rosen(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13, - 1e-13, 1e-6, 100000, expected); - } - - @Test - public void testRosen() throws MathException { - double[] startPoint = point(DIM,0.1); - double[] insigma = point(DIM,0.1); - double[][] boundaries = null; - RealPointValuePair expected = - new RealPointValuePair(point(DIM,1.0),0.0); - doTest(new Rosen(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13, - 1e-13, 1e-6, 100000, expected); - doTest(new Rosen(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13, - 1e-13, 1e-6, 100000, expected); - } - - @Test - public void testMaximize() throws MathException { - double[] startPoint = point(DIM,1.0); - double[] insigma = point(DIM,0.1); - double[][] boundaries = null; - RealPointValuePair expected = - new RealPointValuePair(point(DIM,0.0),1.0); - doTest(new MinusElli(), startPoint, insigma, boundaries, - GoalType.MAXIMIZE, LAMBDA, true, 0, 1.0-1e-13, - 2e-10, 5e-6, 100000, expected); - doTest(new MinusElli(), startPoint, insigma, boundaries, - GoalType.MAXIMIZE, LAMBDA, false, 0, 1.0-1e-13, - 2e-10, 5e-6, 100000, expected); - boundaries = boundaries(DIM,-0.3,0.3); - startPoint = point(DIM,0.1); - doTest(new MinusElli(), startPoint, insigma, boundaries, - GoalType.MAXIMIZE, LAMBDA, true, 0, 1.0-1e-13, - 2e-10, 5e-6, 100000, expected); - } - - @Test - public void testEllipse() throws MathException { - double[] startPoint = point(DIM,1.0); - double[] insigma = point(DIM,0.1); - double[][] boundaries = null; - RealPointValuePair expected = - new RealPointValuePair(point(DIM,0.0),0.0); - doTest(new Elli(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13, - 1e-13, 1e-6, 100000, expected); - doTest(new Elli(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13, - 1e-13, 1e-6, 100000, expected); - } - - @Test - public void testElliRotated() throws MathException { - double[] startPoint = point(DIM,1.0); - double[] insigma = point(DIM,0.1); - double[][] boundaries = null; - RealPointValuePair expected = - new RealPointValuePair(point(DIM,0.0),0.0); - doTest(new ElliRotated(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13, - 1e-13, 1e-6, 100000, expected); - doTest(new ElliRotated(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13, - 1e-13, 1e-6, 100000, expected); - } - - @Test - public void testCigar() throws MathException { - double[] startPoint = point(DIM,1.0); - double[] insigma = point(DIM,0.1); - double[][] boundaries = null; - RealPointValuePair expected = - new RealPointValuePair(point(DIM,0.0),0.0); - doTest(new Cigar(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13, - 1e-13, 1e-6, 200000, expected); - doTest(new Cigar(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13, - 1e-13, 1e-6, 100000, expected); - } - - @Test - public void testTwoAxes() throws MathException { - double[] startPoint = point(DIM,1.0); - double[] insigma = point(DIM,0.1); - double[][] boundaries = null; - RealPointValuePair expected = - new RealPointValuePair(point(DIM,0.0),0.0); - doTest(new TwoAxes(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, 2*LAMBDA, true, 0, 1e-13, - 1e-13, 1e-6, 200000, expected); - doTest(new TwoAxes(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, 2*LAMBDA, false, 0, 1e-13, - 1e-8, 1e-3, 200000, expected); - } - - @Test - public void testCigTab() throws MathException { - double[] startPoint = point(DIM,1.0); - double[] insigma = point(DIM,0.3); - double[][] boundaries = null; - RealPointValuePair expected = - new RealPointValuePair(point(DIM,0.0),0.0); - doTest(new CigTab(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13, - 1e-13, 5e-5, 100000, expected); - doTest(new CigTab(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13, - 1e-13, 5e-5, 100000, expected); - } - - @Test - public void testSphere() throws MathException { - double[] startPoint = point(DIM,1.0); - double[] insigma = point(DIM,0.1); - double[][] boundaries = null; - RealPointValuePair expected = - new RealPointValuePair(point(DIM,0.0),0.0); - doTest(new Sphere(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13, - 1e-13, 1e-6, 100000, expected); - doTest(new Sphere(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13, - 1e-13, 1e-6, 100000, expected); - } - - @Test - public void testTablet() throws MathException { - double[] startPoint = point(DIM,1.0); - double[] insigma = point(DIM,0.1); - double[][] boundaries = null; - RealPointValuePair expected = - new RealPointValuePair(point(DIM,0.0),0.0); - doTest(new Tablet(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13, - 1e-13, 1e-6, 100000, expected); - doTest(new Tablet(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13, - 1e-13, 1e-6, 100000, expected); - } - - @Test - public void testDiffPow() throws MathException { - double[] startPoint = point(DIM,1.0); - double[] insigma = point(DIM,0.1); - double[][] boundaries = null; - RealPointValuePair expected = - new RealPointValuePair(point(DIM,0.0),0.0); - doTest(new DiffPow(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, 10, true, 0, 1e-13, - 1e-8, 1e-1, 100000, expected); - doTest(new DiffPow(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, 10, false, 0, 1e-13, - 1e-8, 2e-1, 100000, expected); - } - - @Test - public void testSsDiffPow() throws MathException { - double[] startPoint = point(DIM,1.0); - double[] insigma = point(DIM,0.1); - double[][] boundaries = null; - RealPointValuePair expected = - new RealPointValuePair(point(DIM,0.0),0.0); - doTest(new SsDiffPow(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, 10, true, 0, 1e-13, - 1e-4, 1e-1, 200000, expected); - doTest(new SsDiffPow(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, 10, false, 0, 1e-13, - 1e-4, 1e-1, 200000, expected); - } - - @Test - public void testAckley() throws MathException { - double[] startPoint = point(DIM,1.0); - double[] insigma = point(DIM,1.0); - double[][] boundaries = null; - RealPointValuePair expected = - new RealPointValuePair(point(DIM,0.0),0.0); - doTest(new Ackley(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, 2*LAMBDA, true, 0, 1e-13, - 1e-9, 1e-5, 100000, expected); - doTest(new Ackley(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, 2*LAMBDA, false, 0, 1e-13, - 1e-9, 1e-5, 100000, expected); - } - - @Test - public void testRastrigin() throws MathException { - double[] startPoint = point(DIM,0.1); - double[] insigma = point(DIM,0.1); - double[][] boundaries = null; - RealPointValuePair expected = - new RealPointValuePair(point(DIM,0.0),0.0); - doTest(new Rastrigin(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, (int)(200*Math.sqrt(DIM)), true, 0, 1e-13, - 1e-13, 1e-6, 200000, expected); - doTest(new Rastrigin(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, (int)(200*Math.sqrt(DIM)), false, 0, 1e-13, - 1e-13, 1e-6, 200000, expected); - } - - @Test - public void testConstrainedRosen() throws MathException { - double[] startPoint = point(DIM,0.1); - double[] insigma = point(DIM,0.1); - double[][] boundaries = boundaries(DIM,-1,2); - RealPointValuePair expected = - new RealPointValuePair(point(DIM,1.0),0.0); - doTest(new Rosen(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, 2*LAMBDA, true, 0, 1e-13, - 1e-13, 1e-6, 100000, expected); - doTest(new Rosen(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, 2*LAMBDA, false, 0, 1e-13, - 1e-13, 1e-6, 100000, expected); - } - - @Test - public void testDiagonalRosen() throws MathException { - double[] startPoint = point(DIM,0.1); - double[] insigma = point(DIM,0.1); - double[][] boundaries = null; - RealPointValuePair expected = - new RealPointValuePair(point(DIM,1.0),0.0); - doTest(new Rosen(), startPoint, insigma, boundaries, - GoalType.MINIMIZE, LAMBDA, false, 1, 1e-13, - 1e-10, 1e-4, 1000000, expected); - } - - /** - * @param func Function to optimize. - * @param startPoint Starting point. - * @param inSigma Individual input sigma. - * @param boundaries Upper / lower point limit. - * @param goal Minimization or maximization. - * @param lambda Population size used for offspring. - * @param isActive Covariance update mechanism. - * @param diagonalOnly Simplified covariance update. - * @param stopValue Termination criteria for optimization. - * @param fTol Tolerance relative error on the objective function. - * @param pointTol Tolerance for checking that the optimum is correct. - * @param maxEvaluations Maximum number of evaluations. - * @param expected Expected point / value. - */ - private void doTest(MultivariateRealFunction func, - double[] startPoint, - double[] inSigma, - double[][] boundaries, - GoalType goal, - int lambda, - boolean isActive, - int diagonalOnly, - double stopValue, - double fTol, - double pointTol, - int maxEvaluations, - RealPointValuePair expected) - throws MathException { - int dim = startPoint.length; - // test diagonalOnly = 0 - slow but normally fewer feval# - MultivariateRealOptimizer optim = - new CMAESOptimizer( - lambda, inSigma, boundaries, 30000, - stopValue, isActive, diagonalOnly, 0, new MersenneTwister(),false); - RealPointValuePair result = optim.optimize(maxEvaluations, func, goal, startPoint); - Assert.assertEquals(expected.getValue(), - result.getValue(), fTol); - for (int i = 0; i < dim; i++) { - Assert.assertEquals(expected.getPoint()[i], - result.getPoint()[i], pointTol); - } - } - - private static double[] point(int n, double value) { - double[] ds = new double[n]; - Arrays.fill(ds, value); - return ds; - } - - private static double[][] boundaries(int dim, - double lower, double upper) { - double[][] boundaries = new double[2][dim]; - for (int i = 0; i < dim; i++) - boundaries[0][i] = lower; - for (int i = 0; i < dim; i++) - boundaries[1][i] = upper; - return boundaries; - } - - private static class Sphere implements MultivariateRealFunction { - - public double value(double[] x) { - double f = 0; - for (int i = 0; i < x.length; ++i) - f += x[i] * x[i]; - return f; - } - } - - private static class Cigar implements MultivariateRealFunction { - private double factor; - - Cigar() { - this(1e3); - } - - Cigar(double axisratio) { - factor = axisratio * axisratio; - } - - public double value(double[] x) { - double f = x[0] * x[0]; - for (int i = 1; i < x.length; ++i) - f += factor * x[i] * x[i]; - return f; - } - } - - private static class Tablet implements MultivariateRealFunction { - private double factor; - - Tablet() { - this(1e3); - } - - Tablet(double axisratio) { - factor = axisratio * axisratio; - } - - public double value(double[] x) { - double f = factor * x[0] * x[0]; - for (int i = 1; i < x.length; ++i) - f += x[i] * x[i]; - return f; - } - } - - private static class CigTab implements MultivariateRealFunction { - private double factor; - - CigTab() { - this(1e4); - } - - CigTab(double axisratio) { - factor = axisratio; - } - - public double value(double[] x) { - int end = x.length - 1; - double f = x[0] * x[0] / factor + factor * x[end] * x[end]; - for (int i = 1; i < end; ++i) - f += x[i] * x[i]; - return f; - } - } - - private static class TwoAxes implements MultivariateRealFunction { - - private double factor; - - TwoAxes() { - this(1e6); - } - - TwoAxes(double axisratio) { - factor = axisratio * axisratio; - } - - public double value(double[] x) { - double f = 0; - for (int i = 0; i < x.length; ++i) - f += (i < x.length / 2 ? factor : 1) * x[i] * x[i]; - return f; - } - } - - private static class ElliRotated implements MultivariateRealFunction { - private Basis B = new Basis(); - private double factor; - - ElliRotated() { - this(1e3); - } - - ElliRotated(double axisratio) { - factor = axisratio * axisratio; - } - - public double value(double[] x) { - double f = 0; - x = B.Rotate(x); - for (int i = 0; i < x.length; ++i) - f += Math.pow(factor, i / (x.length - 1.)) * x[i] * x[i]; - return f; - } - } - - private static class Elli implements MultivariateRealFunction { - - private double factor; - - Elli() { - this(1e3); - } - - Elli(double axisratio) { - factor = axisratio * axisratio; - } - - public double value(double[] x) { - double f = 0; - for (int i = 0; i < x.length; ++i) - f += Math.pow(factor, i / (x.length - 1.)) * x[i] * x[i]; - return f; - } - } - - private static class MinusElli implements MultivariateRealFunction { - - public double value(double[] x) { - return 1.0-(new Elli().value(x)); - } - } - - private static class DiffPow implements MultivariateRealFunction { - - public double value(double[] x) { - double f = 0; - for (int i = 0; i < x.length; ++i) - f += Math.pow(Math.abs(x[i]), 2. + 10 * (double) i - / (x.length - 1.)); - return f; - } - } - - private static class SsDiffPow implements MultivariateRealFunction { - - public double value(double[] x) { - double f = Math.pow(new DiffPow().value(x), 0.25); - return f; - } - } - - private static class Rosen implements MultivariateRealFunction { - - public double value(double[] x) { - double f = 0; - for (int i = 0; i < x.length - 1; ++i) - f += 1e2 * (x[i] * x[i] - x[i + 1]) * (x[i] * x[i] - x[i + 1]) - + (x[i] - 1.) * (x[i] - 1.); - return f; - } - } - - private static class Ackley implements MultivariateRealFunction { - private double axisratio; - - Ackley(double axra) { - axisratio = axra; - } - - public Ackley() { - this(1); - } - - public double value(double[] x) { - double f = 0; - double res2 = 0; - double fac = 0; - for (int i = 0; i < x.length; ++i) { - fac = Math.pow(axisratio, (i - 1.) / (x.length - 1.)); - f += fac * fac * x[i] * x[i]; - res2 += Math.cos(2. * Math.PI * fac * x[i]); - } - f = (20. - 20. * Math.exp(-0.2 * Math.sqrt(f / x.length)) - + Math.exp(1.) - Math.exp(res2 / x.length)); - return f; - } - } - - private static class Rastrigin implements MultivariateRealFunction { - - private double axisratio; - private double amplitude; - - Rastrigin() { - this(1, 10); - } - - Rastrigin(double axisratio, double amplitude) { - this.axisratio = axisratio; - this.amplitude = amplitude; - } - - public double value(double[] x) { - double f = 0; - double fac; - for (int i = 0; i < x.length; ++i) { - fac = Math.pow(axisratio, (i - 1.) / (x.length - 1.)); - if (i == 0 && x[i] < 0) - fac *= 1.; - f += fac * fac * x[i] * x[i] + amplitude - * (1. - Math.cos(2. * Math.PI * fac * x[i])); - } - return f; - } - } - - private static class Basis { - double[][] basis; - Random rand = new Random(2); // use not always the same basis - - double[] Rotate(double[] x) { - GenBasis(x.length); - double[] y = new double[x.length]; - for (int i = 0; i < x.length; ++i) { - y[i] = 0; - for (int j = 0; j < x.length; ++j) - y[i] += basis[i][j] * x[j]; - } - return y; - } - - void GenBasis(int DIM) { - if (basis != null ? basis.length == DIM : false) - return; - - double sp; - int i, j, k; - - /* generate orthogonal basis */ - basis = new double[DIM][DIM]; - for (i = 0; i < DIM; ++i) { - /* sample components gaussian */ - for (j = 0; j < DIM; ++j) - basis[i][j] = rand.nextGaussian(); - /* substract projection of previous vectors */ - for (j = i - 1; j >= 0; --j) { - for (sp = 0., k = 0; k < DIM; ++k) - sp += basis[i][k] * basis[j][k]; /* scalar product */ - for (k = 0; k < DIM; ++k) - basis[i][k] -= sp * basis[j][k]; /* substract */ - } - /* normalize */ - for (sp = 0., k = 0; k < DIM; ++k) - sp += basis[i][k] * basis[i][k]; /* squared norm */ - for (k = 0; k < DIM; ++k) - basis[i][k] /= Math.sqrt(sp); - } - } - } -} diff --git a/src/test/java/org/apache/commons/math/util/MathUtilsTest.java b/src/test/java/org/apache/commons/math/util/MathUtilsTest.java index 5b49797e7..c5787770e 100644 --- a/src/test/java/org/apache/commons/math/util/MathUtilsTest.java +++ b/src/test/java/org/apache/commons/math/util/MathUtilsTest.java @@ -1557,35 +1557,6 @@ public final class MathUtilsTest extends TestCase { } } - public void testCopyOfInt2() { - final int[] source = { Integer.MIN_VALUE, - -1, 0, 1, 3, 113, 4769, - Integer.MAX_VALUE }; - final int offset = 3; - final int[] dest = MathUtils.copyOf(source, source.length - offset); - - assertEquals(dest.length, source.length - offset); - for (int i = 0; i < source.length - offset; i++) { - assertEquals(source[i], dest[i]); - } - } - - public void testCopyOfInt3() { - final int[] source = { Integer.MIN_VALUE, - -1, 0, 1, 3, 113, 4769, - Integer.MAX_VALUE }; - final int offset = 3; - final int[] dest = MathUtils.copyOf(source, source.length + offset); - - assertEquals(dest.length, source.length + offset); - for (int i = 0; i < source.length; i++) { - assertEquals(source[i], dest[i]); - } - for (int i = source.length; i < source.length + offset; i++) { - assertEquals(0, dest[i], 0); - } - } - public void testCopyOfDouble() { final double[] source = { Double.NEGATIVE_INFINITY, -Double.MAX_VALUE, @@ -1602,43 +1573,4 @@ public final class MathUtilsTest extends TestCase { assertEquals(source[i], dest[i], 0); } } - - public void testCopyOfDouble2() { - final double[] source = { Double.NEGATIVE_INFINITY, - -Double.MAX_VALUE, - -1, 0, - Double.MIN_VALUE, - Math.ulp(1d), - 1, 3, 113, 4769, - Double.MAX_VALUE, - Double.POSITIVE_INFINITY }; - final int offset = 3; - final double[] dest = MathUtils.copyOf(source, source.length - offset); - - assertEquals(dest.length, source.length - offset); - for (int i = 0; i < source.length - offset; i++) { - assertEquals(source[i], dest[i], 0); - } - } - - public void testCopyOfDouble3() { - final double[] source = { Double.NEGATIVE_INFINITY, - -Double.MAX_VALUE, - -1, 0, - Double.MIN_VALUE, - Math.ulp(1d), - 1, 3, 113, 4769, - Double.MAX_VALUE, - Double.POSITIVE_INFINITY }; - final int offset = 3; - final double[] dest = MathUtils.copyOf(source, source.length + offset); - - assertEquals(dest.length, source.length + offset); - for (int i = 0; i < source.length; i++) { - assertEquals(source[i], dest[i], 0); - } - for (int i = source.length; i < source.length + offset; i++) { - assertEquals(0, dest[i], 0); - } - } }