1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv;
18
19 import java.util.Arrays;
20 import org.apache.commons.math4.legacy.analysis.MultivariateFunction;
21 import org.apache.commons.math4.legacy.exception.MathUnsupportedOperationException;
22 import org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException;
23 import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
24 import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
25 import org.apache.commons.math4.legacy.optim.ConvergenceChecker;
26 import org.apache.commons.math4.legacy.optim.PointValuePair;
27 import org.apache.commons.math4.legacy.optim.nonlinear.scalar.GoalType;
28 import org.apache.commons.math4.legacy.optim.nonlinear.scalar.MultivariateOptimizer;
29 import org.apache.commons.math4.legacy.optim.univariate.UnivariatePointValuePair;
30 import org.apache.commons.math4.core.jdkmath.JdkMath;
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55 public class PowellOptimizer
56 extends MultivariateOptimizer {
57
58
59
60 private static final double MIN_RELATIVE_TOLERANCE = 2 * JdkMath.ulp(1d);
61
62 private final double relativeThreshold;
63
64 private final double absoluteThreshold;
65
66 private final double lineSearchRelativeThreshold;
67
68 private final double lineSearchAbsoluteThreshold;
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84 public PowellOptimizer(double rel,
85 double abs,
86 ConvergenceChecker<PointValuePair> checker) {
87 this(rel, abs, JdkMath.sqrt(rel), JdkMath.sqrt(abs), checker);
88 }
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103 public PowellOptimizer(double rel,
104 double abs,
105 double lineRel,
106 double lineAbs,
107 ConvergenceChecker<PointValuePair> checker) {
108 super(checker);
109
110 if (rel < MIN_RELATIVE_TOLERANCE) {
111 throw new NumberIsTooSmallException(rel, MIN_RELATIVE_TOLERANCE, true);
112 }
113 if (abs <= 0) {
114 throw new NotStrictlyPositiveException(abs);
115 }
116
117 relativeThreshold = rel;
118 absoluteThreshold = abs;
119 lineSearchRelativeThreshold = lineRel;
120 lineSearchAbsoluteThreshold = lineAbs;
121 }
122
123
124
125
126
127
128
129
130
131
132
133
134 public PowellOptimizer(double rel,
135 double abs) {
136 this(rel, abs, null);
137 }
138
139
140
141
142
143
144
145
146
147
148
149 public PowellOptimizer(double rel,
150 double abs,
151 double lineRel,
152 double lineAbs) {
153 this(rel, abs, lineRel, lineAbs, null);
154 }
155
156
157 @Override
158 protected PointValuePair doOptimize() {
159 checkParameters();
160
161
162 createLineSearch();
163
164 final GoalType goal = getGoalType();
165 final double[] guess = getStartPoint();
166 final MultivariateFunction func = getObjectiveFunction();
167 final int n = guess.length;
168
169 final double[][] direc = new double[n][n];
170 for (int i = 0; i < n; i++) {
171 direc[i][i] = 1;
172 }
173
174 final ConvergenceChecker<PointValuePair> checker
175 = getConvergenceChecker();
176
177 double[] x = guess;
178 double fVal = func.value(x);
179 double[] x1 = x.clone();
180 while (true) {
181 incrementIterationCount();
182
183 double fX = fVal;
184 double fX2 = 0;
185 double delta = 0;
186 int bigInd = 0;
187 double alphaMin = 0;
188
189 for (int i = 0; i < n; i++) {
190 final double[] d = Arrays.copyOf(direc[i], direc[i].length);
191
192 fX2 = fVal;
193
194 final UnivariatePointValuePair optimum = lineSearch(x, d);
195 fVal = optimum.getValue();
196 alphaMin = optimum.getPoint();
197 final double[][] result = newPointAndDirection(x, d, alphaMin);
198 x = result[0];
199
200 if ((fX2 - fVal) > delta) {
201 delta = fX2 - fVal;
202 bigInd = i;
203 }
204 }
205
206
207 boolean stop = 2 * (fX - fVal) <=
208 (relativeThreshold * (JdkMath.abs(fX) + JdkMath.abs(fVal)) +
209 absoluteThreshold);
210
211 final PointValuePair previous = new PointValuePair(x1, fX);
212 final PointValuePair current = new PointValuePair(x, fVal);
213 if (!stop && checker != null) {
214 stop = checker.converged(getIterations(), previous, current);
215 }
216 if (stop) {
217 if (goal == GoalType.MINIMIZE) {
218 return (fVal < fX) ? current : previous;
219 } else {
220 return (fVal > fX) ? current : previous;
221 }
222 }
223
224 final double[] d = new double[n];
225 final double[] x2 = new double[n];
226 for (int i = 0; i < n; i++) {
227 d[i] = x[i] - x1[i];
228 x2[i] = 2 * x[i] - x1[i];
229 }
230
231 x1 = x.clone();
232 fX2 = func.value(x2);
233
234 if (fX > fX2) {
235 double t = 2 * (fX + fX2 - 2 * fVal);
236 double temp = fX - fVal - delta;
237 t *= temp * temp;
238 temp = fX - fX2;
239 t -= delta * temp * temp;
240
241 if (t < 0.0) {
242 final UnivariatePointValuePair optimum = lineSearch(x, d);
243 fVal = optimum.getValue();
244 alphaMin = optimum.getPoint();
245 final double[][] result = newPointAndDirection(x, d, alphaMin);
246 x = result[0];
247
248 final int lastInd = n - 1;
249 direc[bigInd] = direc[lastInd];
250 direc[lastInd] = result[1];
251 }
252 }
253 }
254 }
255
256
257
258
259
260
261
262
263
264
265
266 private double[][] newPointAndDirection(double[] p,
267 double[] d,
268 double optimum) {
269 final int n = p.length;
270 final double[] nP = new double[n];
271 final double[] nD = new double[n];
272 for (int i = 0; i < n; i++) {
273 nD[i] = d[i] * optimum;
274 nP[i] = p[i] + nD[i];
275 }
276
277 final double[][] result = new double[2][];
278 result[0] = nP;
279 result[1] = nD;
280
281 return result;
282 }
283
284
285
286
287
288 private void checkParameters() {
289 if (getLowerBound() != null ||
290 getUpperBound() != null) {
291 throw new MathUnsupportedOperationException(LocalizedFormats.CONSTRAINT);
292 }
293 }
294 }