1 /* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 18 package org.apache.commons.math.optimization.direct; 19 20 import org.apache.commons.math.MaxIterationsExceededException; 21 import org.apache.commons.math.analysis.UnivariateRealFunction; 22 import org.apache.commons.math.FunctionEvaluationException; 23 import org.apache.commons.math.optimization.GoalType; 24 import org.apache.commons.math.optimization.OptimizationException; 25 import org.apache.commons.math.optimization.RealPointValuePair; 26 import org.apache.commons.math.optimization.general.AbstractScalarDifferentiableOptimizer; 27 import org.apache.commons.math.optimization.univariate.AbstractUnivariateRealOptimizer; 28 import org.apache.commons.math.optimization.univariate.BracketFinder; 29 import org.apache.commons.math.optimization.univariate.BrentOptimizer; 30 31 /** 32 * Powell algorithm. 33 * This code is translated and adapted from the Python version of this 34 * algorithm (as implemented in module {@code optimize.py} v0.5 of 35 * <em>SciPy</em>). 36 * 37 * @version $Revision$ $Date$ 38 * @since 2.2 39 */ 40 public class PowellOptimizer 41 extends AbstractScalarDifferentiableOptimizer { 42 /** 43 * Default relative tolerance for line search ({@value}). 44 */ 45 public static final double DEFAULT_LS_RELATIVE_TOLERANCE = 1e-7; 46 /** 47 * Default absolute tolerance for line search ({@value}). 48 */ 49 public static final double DEFAULT_LS_ABSOLUTE_TOLERANCE = 1e-11; 50 /** 51 * Line search. 52 */ 53 private final LineSearch line; 54 55 /** 56 * Constructor with default line search tolerances (see the 57 * {@link #PowellOptimizer(double,double) other constructor}). 58 */ PowellOptimizer()59 public PowellOptimizer() { 60 this(DEFAULT_LS_RELATIVE_TOLERANCE, 61 DEFAULT_LS_ABSOLUTE_TOLERANCE); 62 } 63 64 /** 65 * Constructor with default absolute line search tolerances (see 66 * the {@link #PowellOptimizer(double,double) other constructor}). 67 * 68 * @param lsRelativeTolerance Relative error tolerance for 69 * the line search algorithm ({@link BrentOptimizer}). 70 */ PowellOptimizer(double lsRelativeTolerance)71 public PowellOptimizer(double lsRelativeTolerance) { 72 this(lsRelativeTolerance, 73 DEFAULT_LS_ABSOLUTE_TOLERANCE); 74 } 75 76 /** 77 * @param lsRelativeTolerance Relative error tolerance for 78 * the line search algorithm ({@link BrentOptimizer}). 79 * @param lsAbsoluteTolerance Relative error tolerance for 80 * the line search algorithm ({@link BrentOptimizer}). 81 */ PowellOptimizer(double lsRelativeTolerance, double lsAbsoluteTolerance)82 public PowellOptimizer(double lsRelativeTolerance, 83 double lsAbsoluteTolerance) { 84 line = new LineSearch(lsRelativeTolerance, 85 lsAbsoluteTolerance); 86 } 87 88 /** {@inheritDoc} */ 89 @Override doOptimize()90 protected RealPointValuePair doOptimize() 91 throws FunctionEvaluationException, 92 OptimizationException { 93 final double[] guess = point.clone(); 94 final int n = guess.length; 95 96 final double[][] direc = new double[n][n]; 97 for (int i = 0; i < n; i++) { 98 direc[i][i] = 1; 99 } 100 101 double[] x = guess; 102 double fVal = computeObjectiveValue(x); 103 double[] x1 = x.clone(); 104 while (true) { 105 incrementIterationsCounter(); 106 107 double fX = fVal; 108 double fX2 = 0; 109 double delta = 0; 110 int bigInd = 0; 111 double alphaMin = 0; 112 113 for (int i = 0; i < n; i++) { 114 final double[] d = /* Arrays.*/ copyOf(direc[i], n); // Java 1.5 does not support Arrays.copyOf() 115 116 fX2 = fVal; 117 118 line.search(x, d); 119 fVal = line.getValueAtOptimum(); 120 alphaMin = line.getOptimum(); 121 final double[][] result = newPointAndDirection(x, d, alphaMin); 122 x = result[0]; 123 124 if ((fX2 - fVal) > delta) { 125 delta = fX2 - fVal; 126 bigInd = i; 127 } 128 } 129 130 final RealPointValuePair previous = new RealPointValuePair(x1, fX); 131 final RealPointValuePair current = new RealPointValuePair(x, fVal); 132 if (getConvergenceChecker().converged(getIterations(), previous, current)) { 133 if (goal == GoalType.MINIMIZE) { 134 return (fVal < fX) ? current : previous; 135 } else { 136 return (fVal > fX) ? current : previous; 137 } 138 } 139 140 final double[] d = new double[n]; 141 final double[] x2 = new double[n]; 142 for (int i = 0; i < n; i++) { 143 d[i] = x[i] - x1[i]; 144 x2[i] = 2 * x[i] - x1[i]; 145 } 146 147 x1 = x.clone(); 148 fX2 = computeObjectiveValue(x2); 149 150 if (fX > fX2) { 151 double t = 2 * (fX + fX2 - 2 * fVal); 152 double temp = fX - fVal - delta; 153 t *= temp * temp; 154 temp = fX - fX2; 155 t -= delta * temp * temp; 156 157 if (t < 0.0) { 158 line.search(x, d); 159 fVal = line.getValueAtOptimum(); 160 alphaMin = line.getOptimum(); 161 final double[][] result = newPointAndDirection(x, d, alphaMin); 162 x = result[0]; 163 164 final int lastInd = n - 1; 165 direc[bigInd] = direc[lastInd]; 166 direc[lastInd] = result[1]; 167 } 168 } 169 } 170 } 171 172 /** 173 * Compute a new point (in the original space) and a new direction 174 * vector, resulting from the line search. 175 * The parameters {@code p} and {@code d} will be changed in-place. 176 * 177 * @param p Point used in the line search. 178 * @param d Direction used in the line search. 179 * @param optimum Optimum found by the line search. 180 * @return a 2-element array containing the new point (at index 0) and 181 * the new direction (at index 1). 182 */ newPointAndDirection(double[] p, double[] d, double optimum)183 private double[][] newPointAndDirection(double[] p, 184 double[] d, 185 double optimum) { 186 final int n = p.length; 187 final double[][] result = new double[2][n]; 188 final double[] nP = result[0]; 189 final double[] nD = result[1]; 190 for (int i = 0; i < n; i++) { 191 nD[i] = d[i] * optimum; 192 nP[i] = p[i] + nD[i]; 193 } 194 return result; 195 } 196 197 /** 198 * Class for finding the minimum of the objective function along a given 199 * direction. 200 */ 201 private class LineSearch { 202 /** 203 * Optimizer. 204 */ 205 private final AbstractUnivariateRealOptimizer optim = new BrentOptimizer(); 206 /** 207 * Automatic bracketing. 208 */ 209 private final BracketFinder bracket = new BracketFinder(); 210 /** 211 * Value of the optimum. 212 */ 213 private double optimum = Double.NaN; 214 /** 215 * Value of the objective function at the optimum. 216 */ 217 private double valueAtOptimum = Double.NaN; 218 219 /** 220 * @param relativeTolerance Relative tolerance. 221 * @param absoluteTolerance Absolute tolerance. 222 */ LineSearch(double relativeTolerance, double absoluteTolerance)223 public LineSearch(double relativeTolerance, 224 double absoluteTolerance) { 225 optim.setRelativeAccuracy(relativeTolerance); 226 optim.setAbsoluteAccuracy(absoluteTolerance); 227 } 228 229 /** 230 * Find the minimum of the function {@code f(p + alpha * d)}. 231 * 232 * @param p Starting point. 233 * @param d Search direction. 234 * @throws FunctionEvaluationException if function cannot be evaluated at some test point 235 * @throws OptimizationException if algorithm fails to converge 236 */ search(final double[] p, final double[] d)237 public void search(final double[] p, final double[] d) 238 throws OptimizationException, FunctionEvaluationException { 239 240 // Reset. 241 optimum = Double.NaN; 242 valueAtOptimum = Double.NaN; 243 244 try { 245 final int n = p.length; 246 final UnivariateRealFunction f = new UnivariateRealFunction() { 247 public double value(double alpha) 248 throws FunctionEvaluationException { 249 250 final double[] x = new double[n]; 251 for (int i = 0; i < n; i++) { 252 x[i] = p[i] + alpha * d[i]; 253 } 254 final double obj; 255 obj = computeObjectiveValue(x); 256 return obj; 257 } 258 }; 259 260 bracket.search(f, goal, 0, 1); 261 optimum = optim.optimize(f, goal, 262 bracket.getLo(), 263 bracket.getHi(), 264 bracket.getMid()); 265 valueAtOptimum = optim.getFunctionValue(); 266 } catch (MaxIterationsExceededException e) { 267 throw new OptimizationException(e); 268 } 269 } 270 271 /** 272 * @return the optimum. 273 */ getOptimum()274 public double getOptimum() { 275 return optimum; 276 } 277 /** 278 * @return the value of the function at the optimum. 279 */ getValueAtOptimum()280 public double getValueAtOptimum() { 281 return valueAtOptimum; 282 } 283 } 284 285 /** 286 * Java 1.5 does not support Arrays.copyOf() 287 * 288 * @param source the array to be copied 289 * @param newLen the length of the copy to be returned 290 * @return the copied array, truncated or padded as necessary. 291 */ copyOf(double[] source, int newLen)292 private double[] copyOf(double[] source, int newLen) { 293 double[] output = new double[newLen]; 294 System.arraycopy(source, 0, output, 0, Math.min(source.length, newLen)); 295 return output; 296 } 297 298 } 299