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 * https://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 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22 package org.hipparchus.optim.nonlinear.scalar.noderiv;
23
24 import org.hipparchus.exception.LocalizedCoreFormats;
25 import org.hipparchus.exception.MathIllegalArgumentException;
26 import org.hipparchus.exception.MathRuntimeException;
27 import org.hipparchus.optim.ConvergenceChecker;
28 import org.hipparchus.optim.OptimizationData;
29 import org.hipparchus.optim.PointValuePair;
30 import org.hipparchus.optim.nonlinear.scalar.GoalType;
31 import org.hipparchus.optim.nonlinear.scalar.LineSearch;
32 import org.hipparchus.optim.nonlinear.scalar.MultivariateOptimizer;
33 import org.hipparchus.optim.univariate.UnivariatePointValuePair;
34 import org.hipparchus.util.FastMath;
35
36 /**
37 * Powell's algorithm.
38 * This code is translated and adapted from the Python version of this
39 * algorithm (as implemented in module {@code optimize.py} v0.5 of
40 * <em>SciPy</em>).
41 * <br>
42 * The default stopping criterion is based on the differences of the
43 * function value between two successive iterations. It is however possible
44 * to define a custom convergence checker that might terminate the algorithm
45 * earlier.
46 * <br>
47 * Line search is performed by the {@link LineSearch} class.
48 * <br>
49 * Constraints are not supported: the call to
50 * {@link #optimize(OptimizationData...)} optimize} will throw
51 * {@link MathRuntimeException} if bounds are passed to it.
52 * In order to impose simple constraints, the objective function must be
53 * wrapped in an adapter like
54 * {@link org.hipparchus.optim.nonlinear.scalar.MultivariateFunctionMappingAdapter
55 * MultivariateFunctionMappingAdapter} or
56 * {@link org.hipparchus.optim.nonlinear.scalar.MultivariateFunctionPenaltyAdapter
57 * MultivariateFunctionPenaltyAdapter}.
58 *
59 */
60 public class PowellOptimizer
61 extends MultivariateOptimizer {
62 /**
63 * Minimum relative tolerance.
64 */
65 private static final double MIN_RELATIVE_TOLERANCE = 2 * FastMath.ulp(1d);
66 /**
67 * Relative threshold.
68 */
69 private final double relativeThreshold;
70 /**
71 * Absolute threshold.
72 */
73 private final double absoluteThreshold;
74 /**
75 * Line search.
76 */
77 private final LineSearch line;
78
79 /**
80 * This constructor allows to specify a user-defined convergence checker,
81 * in addition to the parameters that control the default convergence
82 * checking procedure.
83 * <br>
84 * The internal line search tolerances are set to the square-root of their
85 * corresponding value in the multivariate optimizer.
86 *
87 * @param rel Relative threshold.
88 * @param abs Absolute threshold.
89 * @param checker Convergence checker.
90 * @throws MathIllegalArgumentException if {@code abs <= 0}.
91 * @throws MathIllegalArgumentException if {@code rel < 2 * FastMath.ulp(1d)}.
92 */
93 public PowellOptimizer(double rel,
94 double abs,
95 ConvergenceChecker<PointValuePair> checker) {
96 this(rel, abs, FastMath.sqrt(rel), FastMath.sqrt(abs), checker);
97 }
98
99 /**
100 * This constructor allows to specify a user-defined convergence checker,
101 * in addition to the parameters that control the default convergence
102 * checking procedure and the line search tolerances.
103 *
104 * @param rel Relative threshold for this optimizer.
105 * @param abs Absolute threshold for this optimizer.
106 * @param lineRel Relative threshold for the internal line search optimizer.
107 * @param lineAbs Absolute threshold for the internal line search optimizer.
108 * @param checker Convergence checker.
109 * @throws MathIllegalArgumentException if {@code abs <= 0}.
110 * @throws MathIllegalArgumentException if {@code rel < 2 * FastMath.ulp(1d)}.
111 */
112 public PowellOptimizer(double rel,
113 double abs,
114 double lineRel,
115 double lineAbs,
116 ConvergenceChecker<PointValuePair> checker) {
117 super(checker);
118
119 if (rel < MIN_RELATIVE_TOLERANCE) {
120 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
121 rel, MIN_RELATIVE_TOLERANCE);
122 }
123 if (abs <= 0) {
124 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL_BOUND_EXCLUDED,
125 abs, 0);
126 }
127 relativeThreshold = rel;
128 absoluteThreshold = abs;
129
130 // Create the line search optimizer.
131 line = new LineSearch(this,
132 lineRel,
133 lineAbs,
134 1d);
135 }
136
137 /**
138 * The parameters control the default convergence checking procedure.
139 * <br>
140 * The internal line search tolerances are set to the square-root of their
141 * corresponding value in the multivariate optimizer.
142 *
143 * @param rel Relative threshold.
144 * @param abs Absolute threshold.
145 * @throws MathIllegalArgumentException if {@code abs <= 0}.
146 * @throws MathIllegalArgumentException if {@code rel < 2 * FastMath.ulp(1d)}.
147 */
148 public PowellOptimizer(double rel,
149 double abs) {
150 this(rel, abs, null);
151 }
152
153 /**
154 * Builds an instance with the default convergence checking procedure.
155 *
156 * @param rel Relative threshold.
157 * @param abs Absolute threshold.
158 * @param lineRel Relative threshold for the internal line search optimizer.
159 * @param lineAbs Absolute threshold for the internal line search optimizer.
160 * @throws MathIllegalArgumentException if {@code abs <= 0}.
161 * @throws MathIllegalArgumentException if {@code rel < 2 * FastMath.ulp(1d)}.
162 */
163 public PowellOptimizer(double rel,
164 double abs,
165 double lineRel,
166 double lineAbs) {
167 this(rel, abs, lineRel, lineAbs, null);
168 }
169
170 /** {@inheritDoc} */
171 @Override
172 protected PointValuePair doOptimize() {
173 checkParameters();
174
175 final GoalType goal = getGoalType();
176 final double[] guess = getStartPoint();
177 final int n = guess.length;
178
179 final double[][] direc = new double[n][n];
180 for (int i = 0; i < n; i++) {
181 direc[i][i] = 1;
182 }
183
184 final ConvergenceChecker<PointValuePair> checker
185 = getConvergenceChecker();
186
187 double[] x = guess;
188 double fVal = computeObjectiveValue(x);
189 double[] x1 = x.clone();
190 while (true) {
191 incrementIterationCount();
192
193 double fX = fVal;
194 double delta = 0;
195 int bigInd = 0;
196
197 for (int i = 0; i < n; i++) {
198 final double[] d = direc[i].clone();
199
200 final double fX2 = fVal;
201
202 final UnivariatePointValuePair optimum = line.search(x, d);
203 fVal = optimum.getValue();
204 final double alphaMin = optimum.getPoint();
205 final double[][] result = newPointAndDirection(x, d, alphaMin);
206 x = result[0];
207
208 if ((fX2 - fVal) > delta) {
209 delta = fX2 - fVal;
210 bigInd = i;
211 }
212 }
213
214 // Default convergence check.
215 boolean stop = 2 * (fX - fVal) <=
216 (relativeThreshold * (FastMath.abs(fX) + FastMath.abs(fVal)) +
217 absoluteThreshold);
218
219 final PointValuePair previous = new PointValuePair(x1, fX);
220 final PointValuePair current = new PointValuePair(x, fVal);
221 if (!stop && checker != null) { // User-defined stopping criteria.
222 stop = checker.converged(getIterations(), previous, current);
223 }
224 if (stop) {
225 if (goal == GoalType.MINIMIZE) {
226 return (fVal < fX) ? current : previous;
227 } else {
228 return (fVal > fX) ? current : previous;
229 }
230 }
231
232 final double[] d = new double[n];
233 final double[] x2 = new double[n];
234 for (int i = 0; i < n; i++) {
235 d[i] = x[i] - x1[i];
236 x2[i] = 2 * x[i] - x1[i];
237 }
238
239 x1 = x.clone();
240 final double fX2 = computeObjectiveValue(x2);
241
242 if (fX > fX2) {
243 double t = 2 * (fX + fX2 - 2 * fVal);
244 double temp = fX - fVal - delta;
245 t *= temp * temp;
246 temp = fX - fX2;
247 t -= delta * temp * temp;
248
249 if (t < 0.0) {
250 final UnivariatePointValuePair optimum = line.search(x, d);
251 fVal = optimum.getValue();
252 final double alphaMin = optimum.getPoint();
253 final double[][] result = newPointAndDirection(x, d, alphaMin);
254 x = result[0];
255
256 final int lastInd = n - 1;
257 direc[bigInd] = direc[lastInd];
258 direc[lastInd] = result[1];
259 }
260 }
261 }
262 }
263
264 /**
265 * Compute a new point (in the original space) and a new direction
266 * vector, resulting from the line search.
267 *
268 * @param p Point used in the line search.
269 * @param d Direction used in the line search.
270 * @param optimum Optimum found by the line search.
271 * @return a 2-element array containing the new point (at index 0) and
272 * the new direction (at index 1).
273 */
274 private double[][] newPointAndDirection(double[] p,
275 double[] d,
276 double optimum) {
277 final int n = p.length;
278 final double[] nP = new double[n];
279 final double[] nD = new double[n];
280 for (int i = 0; i < n; i++) {
281 nD[i] = d[i] * optimum;
282 nP[i] = p[i] + nD[i];
283 }
284
285 final double[][] result = new double[2][];
286 result[0] = nP;
287 result[1] = nD;
288
289 return result;
290 }
291
292 /**
293 * @throws MathRuntimeException if bounds were passed to the
294 * {@link #optimize(OptimizationData...)} optimize} method.
295 */
296 private void checkParameters() {
297 if (getLowerBound() != null ||
298 getUpperBound() != null) {
299 throw new MathRuntimeException(LocalizedCoreFormats.CONSTRAINT);
300 }
301 }
302 }