1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math.optimization.general;
19
20 import java.awt.geom.Point2D;
21 import java.io.Serializable;
22 import java.util.ArrayList;
23
24 import junit.framework.Test;
25 import junit.framework.TestCase;
26 import junit.framework.TestSuite;
27
28 import org.apache.commons.math.FunctionEvaluationException;
29 import org.apache.commons.math.analysis.DifferentiableMultivariateRealFunction;
30 import org.apache.commons.math.analysis.MultivariateRealFunction;
31 import org.apache.commons.math.analysis.MultivariateVectorialFunction;
32 import org.apache.commons.math.analysis.solvers.BrentSolver;
33 import org.apache.commons.math.linear.BlockRealMatrix;
34 import org.apache.commons.math.linear.RealMatrix;
35 import org.apache.commons.math.optimization.GoalType;
36 import org.apache.commons.math.optimization.OptimizationException;
37 import org.apache.commons.math.optimization.RealPointValuePair;
38 import org.apache.commons.math.optimization.SimpleScalarValueChecker;
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102 public class NonLinearConjugateGradientOptimizerTest
103 extends TestCase {
104
105 public NonLinearConjugateGradientOptimizerTest(String name) {
106 super(name);
107 }
108
109 public void testTrivial() throws FunctionEvaluationException, OptimizationException {
110 LinearProblem problem =
111 new LinearProblem(new double[][] { { 2 } }, new double[] { 3 });
112 NonLinearConjugateGradientOptimizer optimizer =
113 new NonLinearConjugateGradientOptimizer(ConjugateGradientFormula.POLAK_RIBIERE);
114 optimizer.setMaxIterations(100);
115 optimizer.setConvergenceChecker(new SimpleScalarValueChecker(1.0e-6, 1.0e-6));
116 RealPointValuePair optimum =
117 optimizer.optimize(problem, GoalType.MINIMIZE, new double[] { 0 });
118 assertEquals(1.5, optimum.getPoint()[0], 1.0e-10);
119 assertEquals(0.0, optimum.getValue(), 1.0e-10);
120 }
121
122 public void testColumnsPermutation() throws FunctionEvaluationException, OptimizationException {
123
124 LinearProblem problem =
125 new LinearProblem(new double[][] { { 1.0, -1.0 }, { 0.0, 2.0 }, { 1.0, -2.0 } },
126 new double[] { 4.0, 6.0, 1.0 });
127
128 NonLinearConjugateGradientOptimizer optimizer =
129 new NonLinearConjugateGradientOptimizer(ConjugateGradientFormula.POLAK_RIBIERE);
130 optimizer.setMaxIterations(100);
131 optimizer.setConvergenceChecker(new SimpleScalarValueChecker(1.0e-6, 1.0e-6));
132 RealPointValuePair optimum =
133 optimizer.optimize(problem, GoalType.MINIMIZE, new double[] { 0, 0 });
134 assertEquals(7.0, optimum.getPoint()[0], 1.0e-10);
135 assertEquals(3.0, optimum.getPoint()[1], 1.0e-10);
136 assertEquals(0.0, optimum.getValue(), 1.0e-10);
137
138 }
139
140 public void testNoDependency() throws FunctionEvaluationException, OptimizationException {
141 LinearProblem problem = new LinearProblem(new double[][] {
142 { 2, 0, 0, 0, 0, 0 },
143 { 0, 2, 0, 0, 0, 0 },
144 { 0, 0, 2, 0, 0, 0 },
145 { 0, 0, 0, 2, 0, 0 },
146 { 0, 0, 0, 0, 2, 0 },
147 { 0, 0, 0, 0, 0, 2 }
148 }, new double[] { 0.0, 1.1, 2.2, 3.3, 4.4, 5.5 });
149 NonLinearConjugateGradientOptimizer optimizer =
150 new NonLinearConjugateGradientOptimizer(ConjugateGradientFormula.POLAK_RIBIERE);
151 optimizer.setMaxIterations(100);
152 optimizer.setConvergenceChecker(new SimpleScalarValueChecker(1.0e-6, 1.0e-6));
153 RealPointValuePair optimum =
154 optimizer.optimize(problem, GoalType.MINIMIZE, new double[] { 0, 0, 0, 0, 0, 0 });
155 for (int i = 0; i < problem.target.length; ++i) {
156 assertEquals(0.55 * i, optimum.getPoint()[i], 1.0e-10);
157 }
158 }
159
160 public void testOneSet() throws FunctionEvaluationException, OptimizationException {
161
162 LinearProblem problem = new LinearProblem(new double[][] {
163 { 1, 0, 0 },
164 { -1, 1, 0 },
165 { 0, -1, 1 }
166 }, new double[] { 1, 1, 1});
167 NonLinearConjugateGradientOptimizer optimizer =
168 new NonLinearConjugateGradientOptimizer(ConjugateGradientFormula.POLAK_RIBIERE);
169 optimizer.setMaxIterations(100);
170 optimizer.setConvergenceChecker(new SimpleScalarValueChecker(1.0e-6, 1.0e-6));
171 RealPointValuePair optimum =
172 optimizer.optimize(problem, GoalType.MINIMIZE, new double[] { 0, 0, 0 });
173 assertEquals(1.0, optimum.getPoint()[0], 1.0e-10);
174 assertEquals(2.0, optimum.getPoint()[1], 1.0e-10);
175 assertEquals(3.0, optimum.getPoint()[2], 1.0e-10);
176
177 }
178
179 public void testTwoSets() throws FunctionEvaluationException, OptimizationException {
180 final double epsilon = 1.0e-7;
181 LinearProblem problem = new LinearProblem(new double[][] {
182 { 2, 1, 0, 4, 0, 0 },
183 { -4, -2, 3, -7, 0, 0 },
184 { 4, 1, -2, 8, 0, 0 },
185 { 0, -3, -12, -1, 0, 0 },
186 { 0, 0, 0, 0, epsilon, 1 },
187 { 0, 0, 0, 0, 1, 1 }
188 }, new double[] { 2, -9, 2, 2, 1 + epsilon * epsilon, 2});
189
190 NonLinearConjugateGradientOptimizer optimizer =
191 new NonLinearConjugateGradientOptimizer(ConjugateGradientFormula.POLAK_RIBIERE);
192 optimizer.setMaxIterations(100);
193 optimizer.setPreconditioner(new Preconditioner() {
194 public double[] precondition(double[] point, double[] r) {
195 double[] d = r.clone();
196 d[0] /= 72.0;
197 d[1] /= 30.0;
198 d[2] /= 314.0;
199 d[3] /= 260.0;
200 d[4] /= 2 * (1 + epsilon * epsilon);
201 d[5] /= 4.0;
202 return d;
203 }
204 });
205 optimizer.setConvergenceChecker(new SimpleScalarValueChecker(1.0e-13, 1.0e-13));
206
207 RealPointValuePair optimum =
208 optimizer.optimize(problem, GoalType.MINIMIZE, new double[] { 0, 0, 0, 0, 0, 0 });
209 assertEquals( 3.0, optimum.getPoint()[0], 1.0e-10);
210 assertEquals( 4.0, optimum.getPoint()[1], 1.0e-10);
211 assertEquals(-1.0, optimum.getPoint()[2], 1.0e-10);
212 assertEquals(-2.0, optimum.getPoint()[3], 1.0e-10);
213 assertEquals( 1.0 + epsilon, optimum.getPoint()[4], 1.0e-10);
214 assertEquals( 1.0 - epsilon, optimum.getPoint()[5], 1.0e-10);
215
216 }
217
218 public void testNonInversible() throws FunctionEvaluationException, OptimizationException {
219
220 LinearProblem problem = new LinearProblem(new double[][] {
221 { 1, 2, -3 },
222 { 2, 1, 3 },
223 { -3, 0, -9 }
224 }, new double[] { 1, 1, 1 });
225 NonLinearConjugateGradientOptimizer optimizer =
226 new NonLinearConjugateGradientOptimizer(ConjugateGradientFormula.POLAK_RIBIERE);
227 optimizer.setMaxIterations(100);
228 optimizer.setConvergenceChecker(new SimpleScalarValueChecker(1.0e-6, 1.0e-6));
229 RealPointValuePair optimum =
230 optimizer.optimize(problem, GoalType.MINIMIZE, new double[] { 0, 0, 0 });
231 assertTrue(optimum.getValue() > 0.5);
232 }
233
234 public void testIllConditioned() throws FunctionEvaluationException, OptimizationException {
235 LinearProblem problem1 = new LinearProblem(new double[][] {
236 { 10.0, 7.0, 8.0, 7.0 },
237 { 7.0, 5.0, 6.0, 5.0 },
238 { 8.0, 6.0, 10.0, 9.0 },
239 { 7.0, 5.0, 9.0, 10.0 }
240 }, new double[] { 32, 23, 33, 31 });
241 NonLinearConjugateGradientOptimizer optimizer =
242 new NonLinearConjugateGradientOptimizer(ConjugateGradientFormula.POLAK_RIBIERE);
243 optimizer.setMaxIterations(100);
244 optimizer.setConvergenceChecker(new SimpleScalarValueChecker(1.0e-13, 1.0e-13));
245 BrentSolver solver = new BrentSolver();
246 solver.setAbsoluteAccuracy(1.0e-15);
247 solver.setRelativeAccuracy(1.0e-15);
248 optimizer.setLineSearchSolver(solver);
249 RealPointValuePair optimum1 =
250 optimizer.optimize(problem1, GoalType.MINIMIZE, new double[] { 0, 1, 2, 3 });
251 assertEquals(1.0, optimum1.getPoint()[0], 1.0e-5);
252 assertEquals(1.0, optimum1.getPoint()[1], 1.0e-5);
253 assertEquals(1.0, optimum1.getPoint()[2], 1.0e-5);
254 assertEquals(1.0, optimum1.getPoint()[3], 1.0e-5);
255
256 LinearProblem problem2 = new LinearProblem(new double[][] {
257 { 10.00, 7.00, 8.10, 7.20 },
258 { 7.08, 5.04, 6.00, 5.00 },
259 { 8.00, 5.98, 9.89, 9.00 },
260 { 6.99, 4.99, 9.00, 9.98 }
261 }, new double[] { 32, 23, 33, 31 });
262 RealPointValuePair optimum2 =
263 optimizer.optimize(problem2, GoalType.MINIMIZE, new double[] { 0, 1, 2, 3 });
264 assertEquals(-81.0, optimum2.getPoint()[0], 1.0e-1);
265 assertEquals(137.0, optimum2.getPoint()[1], 1.0e-1);
266 assertEquals(-34.0, optimum2.getPoint()[2], 1.0e-1);
267 assertEquals( 22.0, optimum2.getPoint()[3], 1.0e-1);
268
269 }
270
271 public void testMoreEstimatedParametersSimple()
272 throws FunctionEvaluationException, OptimizationException {
273
274 LinearProblem problem = new LinearProblem(new double[][] {
275 { 3.0, 2.0, 0.0, 0.0 },
276 { 0.0, 1.0, -1.0, 1.0 },
277 { 2.0, 0.0, 1.0, 0.0 }
278 }, new double[] { 7.0, 3.0, 5.0 });
279
280 NonLinearConjugateGradientOptimizer optimizer =
281 new NonLinearConjugateGradientOptimizer(ConjugateGradientFormula.POLAK_RIBIERE);
282 optimizer.setMaxIterations(100);
283 optimizer.setConvergenceChecker(new SimpleScalarValueChecker(1.0e-6, 1.0e-6));
284 RealPointValuePair optimum =
285 optimizer.optimize(problem, GoalType.MINIMIZE, new double[] { 7, 6, 5, 4 });
286 assertEquals(0, optimum.getValue(), 1.0e-10);
287
288 }
289
290 public void testMoreEstimatedParametersUnsorted()
291 throws FunctionEvaluationException, OptimizationException {
292 LinearProblem problem = new LinearProblem(new double[][] {
293 { 1.0, 1.0, 0.0, 0.0, 0.0, 0.0 },
294 { 0.0, 0.0, 1.0, 1.0, 1.0, 0.0 },
295 { 0.0, 0.0, 0.0, 0.0, 1.0, -1.0 },
296 { 0.0, 0.0, -1.0, 1.0, 0.0, 1.0 },
297 { 0.0, 0.0, 0.0, -1.0, 1.0, 0.0 }
298 }, new double[] { 3.0, 12.0, -1.0, 7.0, 1.0 });
299 NonLinearConjugateGradientOptimizer optimizer =
300 new NonLinearConjugateGradientOptimizer(ConjugateGradientFormula.POLAK_RIBIERE);
301 optimizer.setMaxIterations(100);
302 optimizer.setConvergenceChecker(new SimpleScalarValueChecker(1.0e-6, 1.0e-6));
303 RealPointValuePair optimum =
304 optimizer.optimize(problem, GoalType.MINIMIZE, new double[] { 2, 2, 2, 2, 2, 2 });
305 assertEquals(0, optimum.getValue(), 1.0e-10);
306 }
307
308 public void testRedundantEquations() throws FunctionEvaluationException, OptimizationException {
309 LinearProblem problem = new LinearProblem(new double[][] {
310 { 1.0, 1.0 },
311 { 1.0, -1.0 },
312 { 1.0, 3.0 }
313 }, new double[] { 3.0, 1.0, 5.0 });
314
315 NonLinearConjugateGradientOptimizer optimizer =
316 new NonLinearConjugateGradientOptimizer(ConjugateGradientFormula.POLAK_RIBIERE);
317 optimizer.setMaxIterations(100);
318 optimizer.setConvergenceChecker(new SimpleScalarValueChecker(1.0e-6, 1.0e-6));
319 RealPointValuePair optimum =
320 optimizer.optimize(problem, GoalType.MINIMIZE, new double[] { 1, 1 });
321 assertEquals(2.0, optimum.getPoint()[0], 1.0e-8);
322 assertEquals(1.0, optimum.getPoint()[1], 1.0e-8);
323
324 }
325
326 public void testInconsistentEquations() throws FunctionEvaluationException, OptimizationException {
327 LinearProblem problem = new LinearProblem(new double[][] {
328 { 1.0, 1.0 },
329 { 1.0, -1.0 },
330 { 1.0, 3.0 }
331 }, new double[] { 3.0, 1.0, 4.0 });
332
333 NonLinearConjugateGradientOptimizer optimizer =
334 new NonLinearConjugateGradientOptimizer(ConjugateGradientFormula.POLAK_RIBIERE);
335 optimizer.setMaxIterations(100);
336 optimizer.setConvergenceChecker(new SimpleScalarValueChecker(1.0e-6, 1.0e-6));
337 RealPointValuePair optimum =
338 optimizer.optimize(problem, GoalType.MINIMIZE, new double[] { 1, 1 });
339 assertTrue(optimum.getValue() > 0.1);
340
341 }
342
343 public void testCircleFitting() throws FunctionEvaluationException, OptimizationException {
344 Circle circle = new Circle();
345 circle.addPoint( 30.0, 68.0);
346 circle.addPoint( 50.0, -6.0);
347 circle.addPoint(110.0, -20.0);
348 circle.addPoint( 35.0, 15.0);
349 circle.addPoint( 45.0, 97.0);
350 NonLinearConjugateGradientOptimizer optimizer =
351 new NonLinearConjugateGradientOptimizer(ConjugateGradientFormula.POLAK_RIBIERE);
352 optimizer.setMaxIterations(100);
353 optimizer.setConvergenceChecker(new SimpleScalarValueChecker(1.0e-30, 1.0e-30));
354 BrentSolver solver = new BrentSolver();
355 solver.setAbsoluteAccuracy(1.0e-13);
356 solver.setRelativeAccuracy(1.0e-15);
357 optimizer.setLineSearchSolver(solver);
358 RealPointValuePair optimum =
359 optimizer.optimize(circle, GoalType.MINIMIZE, new double[] { 98.680, 47.345 });
360 Point2D.Double center = new Point2D.Double(optimum.getPointRef()[0], optimum.getPointRef()[1]);
361 assertEquals(69.960161753, circle.getRadius(center), 1.0e-8);
362 assertEquals(96.075902096, center.x, 1.0e-8);
363 assertEquals(48.135167894, center.y, 1.0e-8);
364 }
365
366 private static class LinearProblem implements DifferentiableMultivariateRealFunction, Serializable {
367
368 private static final long serialVersionUID = 703247177355019415L;
369 final RealMatrix factors;
370 final double[] target;
371 public LinearProblem(double[][] factors, double[] target) {
372 this.factors = new BlockRealMatrix(factors);
373 this.target = target;
374 }
375
376 private double[] gradient(double[] point) {
377 double[] r = factors.operate(point);
378 for (int i = 0; i < r.length; ++i) {
379 r[i] -= target[i];
380 }
381 double[] p = factors.transpose().operate(r);
382 for (int i = 0; i < p.length; ++i) {
383 p[i] *= 2;
384 }
385 return p;
386 }
387
388 public double value(double[] variables) throws FunctionEvaluationException {
389 double[] y = factors.operate(variables);
390 double sum = 0;
391 for (int i = 0; i < y.length; ++i) {
392 double ri = y[i] - target[i];
393 sum += ri * ri;
394 }
395 return sum;
396 }
397
398 public MultivariateVectorialFunction gradient() {
399 return new MultivariateVectorialFunction() {
400 private static final long serialVersionUID = 2621997811350805819L;
401 public double[] value(double[] point) {
402 return gradient(point);
403 }
404 };
405 }
406
407 public MultivariateRealFunction partialDerivative(final int k) {
408 return new MultivariateRealFunction() {
409 private static final long serialVersionUID = -6186178619133562011L;
410 public double value(double[] point) {
411 return gradient(point)[k];
412 }
413 };
414 }
415
416 }
417
418 private static class Circle implements DifferentiableMultivariateRealFunction, Serializable {
419
420 private static final long serialVersionUID = -4711170319243817874L;
421
422 private ArrayList<Point2D.Double> points;
423
424 public Circle() {
425 points = new ArrayList<Point2D.Double>();
426 }
427
428 public void addPoint(double px, double py) {
429 points.add(new Point2D.Double(px, py));
430 }
431
432 public double getRadius(Point2D.Double center) {
433 double r = 0;
434 for (Point2D.Double point : points) {
435 r += point.distance(center);
436 }
437 return r / points.size();
438 }
439
440 private double[] gradient(double[] point) {
441
442
443 Point2D.Double center = new Point2D.Double(point[0], point[1]);
444 double radius = getRadius(center);
445
446
447 double dJdX = 0;
448 double dJdY = 0;
449 for (Point2D.Double pk : points) {
450 double dk = pk.distance(center);
451 dJdX += (center.x - pk.x) * (dk - radius) / dk;
452 dJdY += (center.y - pk.y) * (dk - radius) / dk;
453 }
454 dJdX *= 2;
455 dJdY *= 2;
456
457 return new double[] { dJdX, dJdY };
458
459 }
460
461 public double value(double[] variables)
462 throws IllegalArgumentException, FunctionEvaluationException {
463
464 Point2D.Double center = new Point2D.Double(variables[0], variables[1]);
465 double radius = getRadius(center);
466
467 double sum = 0;
468 for (Point2D.Double point : points) {
469 double di = point.distance(center) - radius;
470 sum += di * di;
471 }
472
473 return sum;
474
475 }
476
477 public MultivariateVectorialFunction gradient() {
478 return new MultivariateVectorialFunction() {
479 private static final long serialVersionUID = 3174909643301201710L;
480 public double[] value(double[] point) {
481 return gradient(point);
482 }
483 };
484 }
485
486 public MultivariateRealFunction partialDerivative(final int k) {
487 return new MultivariateRealFunction() {
488 private static final long serialVersionUID = 3073956364104833888L;
489 public double value(double[] point) {
490 return gradient(point)[k];
491 }
492 };
493 }
494
495 }
496
497 public static Test suite() {
498 return new TestSuite(NonLinearConjugateGradientOptimizerTest.class);
499 }
500
501 }