1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.hipparchus.optim.nonlinear.vector.leastsquares;
18
19 import static org.hamcrest.CoreMatchers.is;
20 import static org.hamcrest.CoreMatchers.not;
21 import static org.hamcrest.CoreMatchers.sameInstance;
22
23 import java.io.IOException;
24 import java.util.Arrays;
25
26 import org.hamcrest.MatcherAssert;
27 import org.hipparchus.analysis.MultivariateMatrixFunction;
28 import org.hipparchus.analysis.MultivariateVectorFunction;
29 import org.hipparchus.exception.LocalizedCoreFormats;
30 import org.hipparchus.exception.MathIllegalArgumentException;
31 import org.hipparchus.exception.MathIllegalStateException;
32 import org.hipparchus.geometry.euclidean.twod.Vector2D;
33 import org.hipparchus.linear.Array2DRowRealMatrix;
34 import org.hipparchus.linear.ArrayRealVector;
35 import org.hipparchus.linear.BlockRealMatrix;
36 import org.hipparchus.linear.DiagonalMatrix;
37 import org.hipparchus.linear.MatrixUtils;
38 import org.hipparchus.linear.RealMatrix;
39 import org.hipparchus.linear.RealVector;
40 import org.hipparchus.optim.ConvergenceChecker;
41 import org.hipparchus.optim.SimpleVectorValueChecker;
42 import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresOptimizer.Optimum;
43 import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem.Evaluation;
44 import org.hipparchus.util.FastMath;
45 import org.hipparchus.util.Pair;
46 import org.junit.Assert;
47 import org.junit.Test;
48
49
50
51
52
53
54
55
56
57
58 public abstract class AbstractSequentialLeastSquaresOptimizerAbstractTest {
59
60
61 public static final double TOl = 1e-10;
62
63
64
65
66 protected SequentialGaussNewtonOptimizer optimizer;
67
68 public LeastSquaresBuilder base() {
69 return new LeastSquaresBuilder()
70 .checkerPair(new SimpleVectorValueChecker(1e-6, 1e-6))
71 .maxEvaluations(100)
72 .maxIterations(getMaxIterations());
73 }
74
75 public LeastSquaresBuilder builder(CircleVectorial c) {
76 final double[] weights = new double[c.getN()];
77 Arrays.fill(weights, 1.0);
78 return base()
79 .model(c.getModelFunction(), c.getModelFunctionJacobian())
80 .target(new double[c.getN()])
81 .weight(new DiagonalMatrix(weights));
82 }
83
84 public LeastSquaresBuilder builder(StatisticalReferenceDataset dataset) {
85 StatisticalReferenceDataset.LeastSquaresProblem problem
86 = dataset.getLeastSquaresProblem();
87 final double[] weights = new double[dataset.getNumObservations()];
88 Arrays.fill(weights, 1.0);
89 return base()
90 .model(problem.getModelFunction(), problem.getModelFunctionJacobian())
91 .target(dataset.getData()[1])
92 .weight(new DiagonalMatrix(weights))
93 .start(dataset.getStartingPoint(0));
94 }
95
96 public void fail(LeastSquaresOptimizer optimizer) {
97 Assert.fail("Expected Exception from: " + optimizer.toString());
98 }
99
100
101
102
103
104
105
106 public void assertEquals(double tolerance, RealVector actual, double... expected){
107 for (int i = 0; i < expected.length; i++) {
108 Assert.assertEquals(expected[i], actual.getEntry(i), tolerance);
109 }
110 Assert.assertEquals(expected.length, actual.getDimension());
111 }
112
113
114
115
116
117 public abstract int getMaxIterations();
118
119
120
121
122
123 public abstract void defineOptimizer(Evaluation evaluation);
124
125 @Test
126 public void testGetIterations() {
127 LeastSquaresProblem lsp = base()
128 .target(new double[]{1})
129 .weight(new DiagonalMatrix(new double[]{1}))
130 .start(new double[]{3})
131 .model(new MultivariateJacobianFunction() {
132 public Pair<RealVector, RealMatrix> value(final RealVector point) {
133 return new Pair<RealVector, RealMatrix>(
134 new ArrayRealVector(
135 new double[]{
136 FastMath.pow(point.getEntry(0), 4)
137 },
138 false),
139 new Array2DRowRealMatrix(
140 new double[][]{
141 {0.25 * FastMath.pow(point.getEntry(0), 3)}
142 },
143 false)
144 );
145 }
146 })
147 .build();
148
149 defineOptimizer(null);
150 Optimum optimum = optimizer.optimize(lsp);
151
152
153 Assert.assertTrue(optimum.getIterations() > 0);
154 }
155
156 @Test
157 public void testTrivial() {
158 LinearProblem problem
159 = new LinearProblem(new double[][]{{2}},
160 new double[]{3});
161 LeastSquaresProblem ls = problem.getBuilder().build();
162
163 defineOptimizer(null);
164 Optimum optimum = optimizer.optimize(ls);
165
166 Assert.assertEquals(0, optimum.getRMS(), TOl);
167 assertEquals(TOl, optimum.getPoint(), 1.5);
168 Assert.assertEquals(0.0, optimum.getResiduals().getEntry(0), TOl);
169 }
170
171 @Test
172 public void testQRColumnsPermutation() {
173
174 LinearProblem problem
175 = new LinearProblem(new double[][]{{1, -1}, {0, 2}, {1, -2}},
176 new double[]{4, 6, 1});
177
178 defineOptimizer(null);
179 Optimum optimum = optimizer.optimize(problem.getBuilder().build());
180
181 Assert.assertEquals(0, optimum.getRMS(), TOl);
182 assertEquals(TOl, optimum.getPoint(), 7, 3);
183 assertEquals(TOl, optimum.getResiduals(), 0, 0, 0);
184 }
185
186 @Test
187 public void testNoDependency() {
188 LinearProblem problem = new LinearProblem(new double[][]{
189 {2, 0, 0, 0, 0, 0},
190 {0, 2, 0, 0, 0, 0},
191 {0, 0, 2, 0, 0, 0},
192 {0, 0, 0, 2, 0, 0},
193 {0, 0, 0, 0, 2, 0},
194 {0, 0, 0, 0, 0, 2}
195 }, new double[]{0, 1.1, 2.2, 3.3, 4.4, 5.5});
196
197 defineOptimizer(null);
198 Optimum optimum = optimizer.optimize(problem.getBuilder().build());
199
200 Assert.assertEquals(0, optimum.getRMS(), TOl);
201 for (int i = 0; i < problem.target.length; ++i) {
202 Assert.assertEquals(0.55 * i, optimum.getPoint().getEntry(i), TOl);
203 }
204 }
205
206 @Test
207 public void testOneSet() {
208 LinearProblem problem = new LinearProblem(new double[][]{
209 {1, 0, 0},
210 {-1, 1, 0},
211 {0, -1, 1}
212 }, new double[]{1, 1, 1});
213
214 defineOptimizer(null);
215 Optimum optimum = optimizer.optimize(problem.getBuilder().build());
216
217 Assert.assertEquals(0, optimum.getRMS(), TOl);
218 assertEquals(TOl, optimum.getPoint(), 1, 2, 3);
219 }
220
221 @Test
222 public void testTwoSets() {
223 double epsilon = 1e-7;
224 LinearProblem problem = new LinearProblem(new double[][]{
225 {2, 1, 0, 4, 0, 0},
226 {-4, -2, 3, -7, 0, 0},
227 {4, 1, -2, 8, 0, 0},
228 {0, -3, -12, -1, 0, 0},
229 {0, 0, 0, 0, epsilon, 1},
230 {0, 0, 0, 0, 1, 1}
231 }, new double[]{2, -9, 2, 2, 1 + epsilon * epsilon, 2});
232
233 defineOptimizer(null);
234 Optimum optimum = optimizer.optimize(problem.getBuilder().build());
235
236 Assert.assertEquals(0, optimum.getRMS(), TOl);
237 assertEquals(TOl, optimum.getPoint(), 3, 4, -1, -2, 1 + epsilon, 1 - epsilon);
238 }
239
240 @Test
241 public void testNonInvertible() throws Exception {
242 try {
243 LinearProblem problem = new LinearProblem(new double[][]{
244 {1, 2, -3},
245 {2, 1, 3},
246 {-3, 0, -9}
247 }, new double[]{1, 1, 1});
248
249 defineOptimizer(null);
250
251 optimizer.optimize(problem.getBuilder().build());
252
253 fail(optimizer);
254 } catch (MathIllegalArgumentException miae) {
255
256 } catch (MathIllegalStateException mise) {
257
258 }
259 }
260
261 @Test
262 public void testIllConditioned() {
263 LinearProblem problem1 = new LinearProblem(new double[][]{
264 {10, 7, 8, 7},
265 {7, 5, 6, 5},
266 {8, 6, 10, 9},
267 {7, 5, 9, 10}
268 }, new double[]{32, 23, 33, 31});
269 final double[] start = {0, 1, 2, 3};
270
271 defineOptimizer(null);
272 Optimum optimum = optimizer
273 .optimize(problem1.getBuilder().start(start).build());
274
275 Assert.assertEquals(0, optimum.getRMS(), TOl);
276 assertEquals(TOl, optimum.getPoint(), 1, 1, 1, 1);
277
278 LinearProblem problem2 = new LinearProblem(new double[][]{
279 {10.00, 7.00, 8.10, 7.20},
280 {7.08, 5.04, 6.00, 5.00},
281 {8.00, 5.98, 9.89, 9.00},
282 {6.99, 4.99, 9.00, 9.98}
283 }, new double[]{32, 23, 33, 31});
284
285 optimum = optimizer.optimize(problem2.getBuilder().start(start).build());
286
287 Assert.assertEquals(0, optimum.getRMS(), TOl);
288 assertEquals(1e-8, optimum.getPoint(), -81, 137, -34, 22);
289 }
290
291 @Test
292 public void testMoreEstimatedParametersSimple() {
293 LinearProblem problem = new LinearProblem(new double[][]{
294 {3, 2, 0, 0},
295 {0, 1, -1, 1},
296 {2, 0, 1, 0}
297 }, new double[]{7, 3, 5});
298
299 defineOptimizer(null);
300 Optimum optimum = optimizer
301 .optimize(problem.getBuilder().start(new double[]{7, 6, 5, 4}).build());
302
303 Assert.assertEquals(0, optimum.getRMS(), TOl);
304 }
305
306 @Test
307 public void testMoreEstimatedParametersUnsorted() {
308 LinearProblem problem = new LinearProblem(new double[][]{
309 {1, 1, 0, 0, 0, 0},
310 {0, 0, 1, 1, 1, 0},
311 {0, 0, 0, 0, 1, -1},
312 {0, 0, -1, 1, 0, 1},
313 {0, 0, 0, -1, 1, 0}
314 }, new double[]{3, 12, -1, 7, 1});
315
316 defineOptimizer(null);
317 Optimum optimum = optimizer.optimize(
318 problem.getBuilder().start(new double[]{2, 2, 2, 2, 2, 2}).build());
319
320 Assert.assertEquals(0, optimum.getRMS(), TOl);
321 RealVector point = optimum.getPoint();
322
323
324 Assert.assertEquals(3, point.getEntry(0) + point.getEntry(1), TOl);
325
326 assertEquals(TOl, point.getSubVector(2, 4), 3, 4, 5, 6);
327 }
328
329 @Test
330 public void testRedundantEquations() {
331 LinearProblem problem = new LinearProblem(new double[][]{
332 {1, 1},
333 {1, -1},
334 {1, 3}
335 }, new double[]{3, 1, 5});
336
337 defineOptimizer(null);
338 Optimum optimum = optimizer
339 .optimize(problem.getBuilder().start(new double[]{1, 1}).build());
340
341 Assert.assertEquals(0, optimum.getRMS(), TOl);
342 assertEquals(TOl, optimum.getPoint(), 2, 1);
343 }
344
345 @Test
346 public void testInconsistentEquations() {
347 LinearProblem problem = new LinearProblem(new double[][]{
348 {1, 1},
349 {1, -1},
350 {1, 3}
351 }, new double[]{3, 1, 4});
352
353 defineOptimizer(null);
354 Optimum optimum = optimizer
355 .optimize(problem.getBuilder().start(new double[]{1, 1}).build());
356
357
358 Assert.assertTrue(optimum.getRMS() > 0.1);
359 }
360
361 @Test
362 public void testInconsistentSizes1() {
363 try {
364 LinearProblem problem
365 = new LinearProblem(new double[][]{{1, 0},
366 {0, 1}},
367 new double[]{-1, 1});
368
369 defineOptimizer(null);
370
371 Optimum optimum = optimizer.optimize(problem.getBuilder().build());
372
373 Assert.assertEquals(0, optimum.getRMS(), TOl);
374 assertEquals(TOl, optimum.getPoint(), -1, 1);
375
376
377 optimizer.optimize(
378 problem.getBuilder().weight(new DiagonalMatrix(new double[]{1})).build());
379
380 fail(optimizer);
381 } catch (MathIllegalArgumentException e) {
382 Assert.assertEquals(LocalizedCoreFormats.DIMENSIONS_MISMATCH, e.getSpecifier());
383 }
384 }
385
386 @Test
387 public void testInconsistentSizes2() {
388 try {
389 LinearProblem problem
390 = new LinearProblem(new double[][]{{1, 0}, {0, 1}},
391 new double[]{-1, 1});
392
393 defineOptimizer(null);
394 Optimum optimum = optimizer.optimize(problem.getBuilder().build());
395
396 Assert.assertEquals(0, optimum.getRMS(), TOl);
397 assertEquals(TOl, optimum.getPoint(), -1, 1);
398
399
400 optimizer.optimize(
401 problem.getBuilder()
402 .target(new double[]{1})
403 .weight(new DiagonalMatrix(new double[]{1}))
404 .build()
405 );
406
407 fail(optimizer);
408 } catch (MathIllegalArgumentException e) {
409 Assert.assertEquals(LocalizedCoreFormats.DIMENSIONS_MISMATCH, e.getSpecifier());
410 }
411 }
412
413 @Test
414 public void testSequential() throws Exception {
415
416 CircleVectorial circleAll = new CircleVectorial();
417 CircleVectorial circleFirst = new CircleVectorial();
418 CircleVectorial circleSecond = new CircleVectorial();
419 circleAll.addPoint( 30.0, 68.0);
420 circleFirst.addPoint( 30.0, 68.0);
421 circleAll.addPoint( 50.0, -6.0);
422 circleFirst.addPoint( 50.0, -6.0);
423 circleAll.addPoint(110.0, -20.0);
424 circleFirst.addPoint(110.0, -20.0);
425 circleAll.addPoint( 35.0, 15.0);
426 circleSecond.addPoint( 35.0, 15.0);
427 circleAll.addPoint( 45.0, 97.0);
428 circleSecond.addPoint( 45.0, 97.0);
429
430
431 defineOptimizer(null);
432 Optimum oneRun = optimizer.optimize(builder(circleAll).
433 checkerPair(new SimpleVectorValueChecker(1e-3, 1e-3)).
434 start(new double[] { 98.680, 47.345 }).
435 build());
436 Assert.assertEquals(2, oneRun.getJacobian().getColumnDimension());
437 Vector2D oneShotCenter = new Vector2D(oneRun.getPoint().getEntry(0), oneRun.getPoint().getEntry(1));
438 Assert.assertEquals(96.075901, oneShotCenter.getX(), 1.0e-6);
439 Assert.assertEquals(48.135169, oneShotCenter.getY(), 1.0e-6);
440 Assert.assertEquals(69.960161, circleAll.getRadius(oneShotCenter), 1.0e-6);
441
442
443
444 defineOptimizer(null);
445 Optimum firstRun = optimizer.optimize(builder(circleFirst).
446 checkerPair(new SimpleVectorValueChecker(1e-3, 1e-3)).
447 start(new double[] { 98.680, 47.345 }).
448 build());
449 Assert.assertEquals(2, firstRun.getJacobian().getColumnDimension());
450 Vector2D firstRunCenter = new Vector2D(firstRun.getPoint().getEntry(0), firstRun.getPoint().getEntry(1));
451 Assert.assertEquals(93.650000, firstRunCenter.getX(), 1.0e-6);
452 Assert.assertEquals(45.500000, firstRunCenter.getY(), 1.0e-6);
453 Assert.assertEquals(67.896265, circleAll.getRadius(firstRunCenter), 1.0e-6);
454
455
456
457
458 optimizer = optimizer.withAPrioriData(firstRun.getPoint(), firstRun.getCovariances(1.0e-8));
459 Optimum secondRun = optimizer.optimize(builder(circleSecond).
460 checkerPair(new SimpleVectorValueChecker(1e-3, 1e-3)).
461 start(new double[] { firstRunCenter.getX(), firstRunCenter.getY() }).
462 build());
463 Assert.assertEquals(2, secondRun.getJacobian().getColumnDimension());
464 Vector2D secondRunCenter = new Vector2D(secondRun.getPoint().getEntry(0), secondRun.getPoint().getEntry(1));
465 Assert.assertEquals(97.070437, secondRunCenter.getX(), 1.0e-6);
466 Assert.assertEquals(49.039898, secondRunCenter.getY(), 1.0e-6);
467 Assert.assertEquals(70.789016, circleAll.getRadius(secondRunCenter), 1.0e-6);
468
469 }
470
471 protected void doTestStRD(final StatisticalReferenceDataset dataset,
472 final double errParams, final double errParamsSd) {
473
474 defineOptimizer(null);
475 final Optimum optimum = optimizer.optimize(builder(dataset).build());
476
477 final RealVector actual = optimum.getPoint();
478 for (int i = 0; i < actual.getDimension(); i++) {
479 double expected = dataset.getParameter(i);
480 double delta = FastMath.abs(errParams * expected);
481 Assert.assertEquals(dataset.getName() + ", param #" + i,
482 expected, actual.getEntry(i), delta);
483 }
484 }
485
486 @Test
487 public void testKirby2() throws IOException {
488 doTestStRD(StatisticalReferenceDatasetFactory.createKirby2(), 1E-7, 1E-7);
489 }
490
491 @Test
492 public void testHahn1() throws IOException {
493 doTestStRD(StatisticalReferenceDatasetFactory.createHahn1(), 1E-7, 1E-4);
494 }
495
496 @Test
497 public void testPointCopy() {
498 LinearProblem problem = new LinearProblem(new double[][]{
499 {1, 0, 0},
500 {-1, 1, 0},
501 {0, -1, 1}
502 }, new double[]{1, 1, 1});
503
504 final boolean[] checked = {false};
505
506 final LeastSquaresBuilder builder = problem.getBuilder()
507 .checker(new ConvergenceChecker<Evaluation>() {
508 public boolean converged(int iteration, Evaluation previous, Evaluation current) {
509 MatcherAssert.assertThat(
510 previous.getPoint(),
511 not(sameInstance(current.getPoint())));
512 Assert.assertArrayEquals(new double[3], previous.getPoint().toArray(), 0);
513 Assert.assertArrayEquals(new double[] {1, 2, 3}, current.getPoint().toArray(), TOl);
514 checked[0] = true;
515 return true;
516 }
517 });
518 defineOptimizer(null);
519 optimizer.optimize(builder.build());
520
521 MatcherAssert.assertThat(checked[0], is(true));
522 }
523
524 @Test
525 public void testPointDifferentDim() {
526 LinearProblem problem
527 = new LinearProblem(new double[][]{{2}},
528 new double[]{3});
529
530 LeastSquaresProblem lsp = problem.getBuilder().build();
531 defineOptimizer(new AbstractEvaluation(2){
532 public RealMatrix getJacobian() {
533 return MatrixUtils.createRealMatrix(2, 2);
534 }
535 public RealVector getPoint() {
536 return MatrixUtils.createRealVector(2);
537 }
538 public RealVector getResiduals() {
539 return MatrixUtils.createRealVector(2);
540 }
541 });
542 try {
543 optimizer.optimize(lsp);
544 fail(optimizer);
545 } catch (MathIllegalStateException mise) {
546 Assert.assertEquals(LocalizedCoreFormats.DIMENSIONS_MISMATCH, mise.getSpecifier());
547 }
548 }
549
550 class LinearProblem {
551 private final RealMatrix factors;
552 private final double[] target;
553
554 public LinearProblem(double[][] factors, double[] target) {
555 this.factors = new BlockRealMatrix(factors);
556 this.target = target;
557 }
558
559 public double[] getTarget() {
560 return target;
561 }
562
563 public MultivariateVectorFunction getModelFunction() {
564 return new MultivariateVectorFunction() {
565 public double[] value(double[] params) {
566 return factors.operate(params);
567 }
568 };
569 }
570
571 public MultivariateMatrixFunction getModelFunctionJacobian() {
572 return new MultivariateMatrixFunction() {
573 public double[][] value(double[] params) {
574 return factors.getData();
575 }
576 };
577 }
578
579 public LeastSquaresBuilder getBuilder() {
580 final double[] weights = new double[target.length];
581 Arrays.fill(weights, 1.0);
582 return base()
583 .model(getModelFunction(), getModelFunctionJacobian())
584 .target(target)
585 .weight(new DiagonalMatrix(weights))
586 .start(new double[factors.getColumnDimension()]);
587 }
588 }
589
590 }