1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math4.legacy.fitting.leastsquares;
19
20 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
21 import org.apache.commons.math4.legacy.exception.TooManyEvaluationsException;
22 import org.apache.commons.math4.legacy.fitting.leastsquares.LeastSquaresOptimizer.Optimum;
23 import org.apache.commons.math4.legacy.fitting.leastsquares.LeastSquaresProblem.Evaluation;
24 import org.apache.commons.math4.legacy.linear.DiagonalMatrix;
25 import org.apache.commons.math4.legacy.linear.RealMatrix;
26 import org.apache.commons.math4.legacy.linear.RealVector;
27 import org.apache.commons.math4.legacy.linear.SingularMatrixException;
28 import org.apache.commons.math4.legacy.optim.ConvergenceChecker;
29 import org.apache.commons.math4.core.jdkmath.JdkMath;
30 import org.apache.commons.numbers.core.Precision;
31 import org.junit.Assert;
32 import org.junit.Test;
33
34
35
36
37
38
39
40
41
42 public class LevenbergMarquardtOptimizerTest
43 extends AbstractLeastSquaresOptimizerAbstractTest{
44
45 public LeastSquaresBuilder builder(BevingtonProblem problem){
46 return base()
47 .model(problem.getModelFunction(), problem.getModelFunctionJacobian());
48 }
49
50 public LeastSquaresBuilder builder(CircleProblem problem){
51 return base()
52 .model(problem.getModelFunction(), problem.getModelFunctionJacobian())
53 .target(problem.target())
54 .weight(new DiagonalMatrix(problem.weight()));
55 }
56
57 @Override
58 public int getMaxIterations() {
59 return 25;
60 }
61
62 @Override
63 public LeastSquaresOptimizer getOptimizer() {
64 return new LevenbergMarquardtOptimizer();
65 }
66
67 @Override
68 @Test
69 public void testNonInvertible() {
70 try{
71
72
73
74
75 LinearProblem problem = new LinearProblem(new double[][] {
76 { 1, 2, -3 },
77 { 2, 1, 3 },
78 { -3, 0, -9 }
79 }, new double[] { 1, 1, 1 });
80
81 final Optimum optimum = optimizer.optimize(
82 problem.getBuilder().maxIterations(20).build());
83
84
85 Assert.assertTrue(JdkMath.sqrt(problem.getTarget().length) * optimum.getRMS() > 0.6);
86
87 optimum.getCovariances(1.5e-14);
88
89 fail(optimizer);
90 }catch (SingularMatrixException e){
91
92 }
93 }
94
95 @Test
96 public void testControlParameters() {
97 CircleVectorial circle = new CircleVectorial();
98 circle.addPoint( 30.0, 68.0);
99 circle.addPoint( 50.0, -6.0);
100 circle.addPoint(110.0, -20.0);
101 circle.addPoint( 35.0, 15.0);
102 circle.addPoint( 45.0, 97.0);
103 checkEstimate(
104 circle, 0.1, 10, 1.0e-14, 1.0e-16, 1.0e-10, false);
105 checkEstimate(
106 circle, 0.1, 10, 1.0e-15, 1.0e-17, 1.0e-10, false);
107 checkEstimate(
108 circle, 0.1, 5, 1.0e-15, 1.0e-16, 1.0e-10, true);
109 circle.addPoint(300, -300);
110
111
112 checkEstimate(
113 circle, 0.1, 20, 1.0e-18, 1.0e-16, 1.0e-10, false);
114 }
115
116 private void checkEstimate(CircleVectorial circle,
117 double initialStepBoundFactor, int maxCostEval,
118 double costRelativeTolerance, double parRelativeTolerance,
119 double orthoTolerance, boolean shouldFail) {
120 try {
121 final LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer()
122 .withInitialStepBoundFactor(initialStepBoundFactor)
123 .withCostRelativeTolerance(costRelativeTolerance)
124 .withParameterRelativeTolerance(parRelativeTolerance)
125 .withOrthoTolerance(orthoTolerance)
126 .withRankingThreshold(Precision.SAFE_MIN);
127
128 final LeastSquaresProblem problem = builder(circle)
129 .maxEvaluations(maxCostEval)
130 .maxIterations(100)
131 .start(new double[] { 98.680, 47.345 })
132 .build();
133
134 optimizer.optimize(problem);
135
136 Assert.assertTrue(!shouldFail);
137
138 } catch (DimensionMismatchException ee) {
139 Assert.assertTrue(shouldFail);
140 } catch (TooManyEvaluationsException ee) {
141 Assert.assertTrue(shouldFail);
142 }
143 }
144
145
146
147
148
149
150
151
152 @Test
153 public void testBevington() {
154 final double[][] dataPoints = {
155
156 { 15, 30, 45, 60, 75, 90, 105, 120, 135, 150,
157 165, 180, 195, 210, 225, 240, 255, 270, 285, 300,
158 315, 330, 345, 360, 375, 390, 405, 420, 435, 450,
159 465, 480, 495, 510, 525, 540, 555, 570, 585, 600,
160 615, 630, 645, 660, 675, 690, 705, 720, 735, 750,
161 765, 780, 795, 810, 825, 840, 855, 870, 885, },
162
163 { 775, 479, 380, 302, 185, 157, 137, 119, 110, 89,
164 74, 61, 66, 68, 48, 54, 51, 46, 55, 29,
165 28, 37, 49, 26, 35, 29, 31, 24, 25, 35,
166 24, 30, 26, 28, 21, 18, 20, 27, 17, 17,
167 14, 17, 24, 11, 22, 17, 12, 10, 13, 16,
168 9, 9, 14, 21, 17, 13, 12, 18, 10, },
169 };
170 final double[] start = {10, 900, 80, 27, 225};
171
172 final BevingtonProblem problem = new BevingtonProblem();
173
174 final int len = dataPoints[0].length;
175 final double[] weights = new double[len];
176 for (int i = 0; i < len; i++) {
177 problem.addPoint(dataPoints[0][i],
178 dataPoints[1][i]);
179
180 weights[i] = 1 / dataPoints[1][i];
181 }
182
183 final Optimum optimum = optimizer.optimize(
184 builder(problem)
185 .target(dataPoints[1])
186 .weight(new DiagonalMatrix(weights))
187 .start(start)
188 .maxIterations(20)
189 .build()
190 );
191
192 final RealVector solution = optimum.getPoint();
193 final double[] expectedSolution = { 10.4, 958.3, 131.4, 33.9, 205.0 };
194
195 final RealMatrix covarMatrix = optimum.getCovariances(1e-14);
196 final double[][] expectedCovarMatrix = {
197 { 3.38, -3.69, 27.98, -2.34, -49.24 },
198 { -3.69, 2492.26, 81.89, -69.21, -8.9 },
199 { 27.98, 81.89, 468.99, -44.22, -615.44 },
200 { -2.34, -69.21, -44.22, 6.39, 53.80 },
201 { -49.24, -8.9, -615.44, 53.8, 929.45 }
202 };
203
204 final int numParams = expectedSolution.length;
205
206
207 for (int i = 0; i < numParams; i++) {
208 final double error = JdkMath.sqrt(expectedCovarMatrix[i][i]);
209 Assert.assertEquals("Parameter " + i, expectedSolution[i], solution.getEntry(i), error);
210 }
211
212
213
214 for (int i = 0; i < numParams; i++) {
215 for (int j = 0; j < numParams; j++) {
216 Assert.assertEquals("Covariance matrix [" + i + "][" + j + "]",
217 expectedCovarMatrix[i][j],
218 covarMatrix.getEntry(i, j),
219 JdkMath.abs(0.1 * expectedCovarMatrix[i][j]));
220 }
221 }
222
223
224 final double chi2 = optimum.getChiSquare();
225 final double cost = optimum.getCost();
226 final double rms = optimum.getRMS();
227 final double reducedChi2 = optimum.getReducedChiSquare(start.length);
228
229
230
231 final double expectedChi2 = 66.07852350839286;
232 final double expectedReducedChi2 = 1.2014277001525975;
233 final double expectedCost = 8.128869755900439;
234 final double expectedRms = 1.0582887010256337;
235
236 final double tol = 1e-14;
237 Assert.assertEquals(expectedChi2, chi2, tol);
238 Assert.assertEquals(expectedReducedChi2, reducedChi2, tol);
239 Assert.assertEquals(expectedCost, cost, tol);
240 Assert.assertEquals(expectedRms, rms, tol);
241 }
242
243 @Test
244 public void testCircleFitting2() {
245 final double xCenter = 123.456;
246 final double yCenter = 654.321;
247 final double xSigma = 10;
248 final double ySigma = 15;
249 final double radius = 111.111;
250
251 final RandomCirclePointGenerator factory
252 = new RandomCirclePointGenerator(xCenter, yCenter, radius,
253 xSigma, ySigma);
254 final CircleProblem circle = new CircleProblem(xSigma, ySigma);
255
256 final int numPoints = 10;
257 factory.samples(numPoints).forEach(circle::addPoint);
258
259
260 final double[] init = { 118, 659, 115 };
261
262 final Optimum optimum = optimizer.optimize(
263 builder(circle).maxIterations(50).start(init).build());
264
265 final double[] paramFound = optimum.getPoint().toArray();
266
267
268 final double[] asymptoticStandardErrorFound = optimum.getSigma(1e-14).toArray();
269
270
271 Assert.assertEquals("Delta=" + 2 * asymptoticStandardErrorFound[0], xCenter, paramFound[0], 2 * asymptoticStandardErrorFound[0]);
272 Assert.assertEquals("Delta=" + 2 * asymptoticStandardErrorFound[1], yCenter, paramFound[1], 2 * asymptoticStandardErrorFound[1]);
273 Assert.assertEquals("Delta=" + 2 * asymptoticStandardErrorFound[2], radius, paramFound[2], asymptoticStandardErrorFound[2]);
274 }
275
276 @Test
277 public void testParameterValidator() {
278
279 final double xCenter = 123.456;
280 final double yCenter = 654.321;
281 final double xSigma = 10;
282 final double ySigma = 15;
283 final double radius = 111.111;
284 final RandomCirclePointGenerator factory
285 = new RandomCirclePointGenerator(xCenter, yCenter, radius,
286 xSigma, ySigma);
287 final CircleProblem circle = new CircleProblem(xSigma, ySigma);
288
289 final int numPoints = 10;
290 factory.samples(numPoints).forEach(circle::addPoint);
291
292
293 final double[] init = { 118, 659, 115 };
294
295 final Optimum optimum = optimizer.optimize(
296 builder(circle).maxIterations(50).start(init).build());
297
298 final int numEval = optimum.getEvaluations();
299 Assert.assertTrue(numEval > 1);
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314 final ParameterValidator cheatValidator
315 = new ParameterValidator() {
316 @Override
317 public RealVector validate(RealVector params) {
318
319 final RealVector direction = optimum.getPoint().subtract(params);
320 return params.add(direction.mapMultiply(0.75));
321 }
322 };
323
324 final Optimum cheatOptimum
325 = optimizer.optimize(builder(circle).maxIterations(50).start(init).parameterValidator(cheatValidator).build());
326 final int cheatNumEval = cheatOptimum.getEvaluations();
327 Assert.assertTrue("n=" + numEval + " nc=" + cheatNumEval, cheatNumEval < numEval);
328
329 }
330
331 @Test
332 public void testEvaluationCount() {
333
334 LeastSquaresProblem lsp = new LinearProblem(new double[][] {{1}}, new double[] {1})
335 .getBuilder()
336 .checker(new ConvergenceChecker<Evaluation>() {
337 @Override
338 public boolean converged(int iteration, Evaluation previous, Evaluation current) {
339 return true;
340 }
341 })
342 .build();
343
344
345 Optimum optimum = optimizer.optimize(lsp);
346
347
348
349 Assert.assertEquals(1, optimum.getIterations());
350 Assert.assertEquals(2, optimum.getEvaluations());
351 }
352 }