001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018 package org.apache.commons.math.optimization.direct;
019
020 import org.apache.commons.math.MaxIterationsExceededException;
021 import org.apache.commons.math.analysis.UnivariateRealFunction;
022 import org.apache.commons.math.FunctionEvaluationException;
023 import org.apache.commons.math.optimization.GoalType;
024 import org.apache.commons.math.optimization.OptimizationException;
025 import org.apache.commons.math.optimization.RealPointValuePair;
026 import org.apache.commons.math.optimization.general.AbstractScalarDifferentiableOptimizer;
027 import org.apache.commons.math.optimization.univariate.AbstractUnivariateRealOptimizer;
028 import org.apache.commons.math.optimization.univariate.BracketFinder;
029 import org.apache.commons.math.optimization.univariate.BrentOptimizer;
030
031 /**
032 * Powell algorithm.
033 * This code is translated and adapted from the Python version of this
034 * algorithm (as implemented in module {@code optimize.py} v0.5 of
035 * <em>SciPy</em>).
036 *
037 * @version $Revision$ $Date$
038 * @since 2.2
039 */
040 public class PowellOptimizer
041 extends AbstractScalarDifferentiableOptimizer {
042 /**
043 * Default relative tolerance for line search ({@value}).
044 */
045 public static final double DEFAULT_LS_RELATIVE_TOLERANCE = 1e-7;
046 /**
047 * Default absolute tolerance for line search ({@value}).
048 */
049 public static final double DEFAULT_LS_ABSOLUTE_TOLERANCE = 1e-11;
050 /**
051 * Line search.
052 */
053 private final LineSearch line;
054
055 /**
056 * Constructor with default line search tolerances (see the
057 * {@link #PowellOptimizer(double,double) other constructor}).
058 */
059 public PowellOptimizer() {
060 this(DEFAULT_LS_RELATIVE_TOLERANCE,
061 DEFAULT_LS_ABSOLUTE_TOLERANCE);
062 }
063
064 /**
065 * Constructor with default absolute line search tolerances (see
066 * the {@link #PowellOptimizer(double,double) other constructor}).
067 *
068 * @param lsRelativeTolerance Relative error tolerance for
069 * the line search algorithm ({@link BrentOptimizer}).
070 */
071 public PowellOptimizer(double lsRelativeTolerance) {
072 this(lsRelativeTolerance,
073 DEFAULT_LS_ABSOLUTE_TOLERANCE);
074 }
075
076 /**
077 * @param lsRelativeTolerance Relative error tolerance for
078 * the line search algorithm ({@link BrentOptimizer}).
079 * @param lsAbsoluteTolerance Relative error tolerance for
080 * the line search algorithm ({@link BrentOptimizer}).
081 */
082 public PowellOptimizer(double lsRelativeTolerance,
083 double lsAbsoluteTolerance) {
084 line = new LineSearch(lsRelativeTolerance,
085 lsAbsoluteTolerance);
086 }
087
088 /** {@inheritDoc} */
089 @Override
090 protected RealPointValuePair doOptimize()
091 throws FunctionEvaluationException,
092 OptimizationException {
093 final double[] guess = point.clone();
094 final int n = guess.length;
095
096 final double[][] direc = new double[n][n];
097 for (int i = 0; i < n; i++) {
098 direc[i][i] = 1;
099 }
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 */
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 */
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 */
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 */
274 public double getOptimum() {
275 return optimum;
276 }
277 /**
278 * @return the value of the function at the optimum.
279 */
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 */
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 }