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