Merge branch 'issue-738' into develop

e9b8893e · Luc Maisonobe · 6bc05440 · 9c810c2b · e9b8893e · e9b8893e
Commit e9b8893e authored Dec 05, 2020 by Luc Maisonobe
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -21,6 +21,9 @@
  </properties>
  <body>
    <release version="10.3" date="TBD" description="TBD">
+      <action dev="luc" type="add" issue="738">
+        Added user-defined max iteration and convergence criterion in SecularAndHarmonic.
+      </action>
      <action dev="luc" type="add" issue="737">
        Added loading of AGI LeapSecond.dat files.
      </action>

--- a/src/main/java/org/orekit/utils/SecularAndHarmonic.java
+++ b/src/main/java/org/orekit/utils/SecularAndHarmonic.java
@@ -54,6 +54,16 @@ public class SecularAndHarmonic {
    /** Observed points. */
    private List<WeightedObservedPoint> observedPoints;

+    /** RMS for convergence.
+     * @since 10.3
+     */
+    private double convergenceRMS;
+
+    /** Maximum number of iterations.
+     * @since 10.3
+     */
+    private int maxIter;
+
    /** Simple constructor.
     * @param secularDegree degree of polynomial secular part
     * @param pulsations pulsations of harmonic part
@@ -62,6 +72,8 @@ public class SecularAndHarmonic {
        this.secularDegree  = secularDegree;
        this.pulsations     = pulsations.clone();
        this.observedPoints = new ArrayList<WeightedObservedPoint>();
+        this.convergenceRMS = 0.0;
+        this.maxIter        = Integer.MAX_VALUE;
    }

    /** Reset fitting.
@@ -75,6 +87,26 @@ public class SecularAndHarmonic {
        observedPoints.clear();
    }

+    /** Set RMS for convergence.
+     * <p>
+     * The RMS is the square-root of the sum of squared of
+     * the residuals, divided by the number of measurements.
+     * </p>
+     * @param convergenceRMS RMS below which convergence is considered to have been reached
+     * @since 10.3
+     */
+    public void setConvergenceRMS(final double convergenceRMS) {
+        this.convergenceRMS = convergenceRMS;
+    }
+
+    /** Set maximum number of iterations.
+     * @param maxIter maximum number of iterations
+     * @since 10.3
+     */
+    public void setMaxIter(final int maxIter) {
+        this.maxIter = maxIter;
+    }
+
    /** Add a fitting point.
     * @param date date of the point
     * @param osculatingValue osculating value
@@ -130,7 +162,8 @@ public class SecularAndHarmonic {
                // build a new least squares problem set up to fit a secular and harmonic curve to the observed points
                return new LeastSquaresBuilder().
                        maxEvaluations(Integer.MAX_VALUE).
-                        maxIterations(Integer.MAX_VALUE).
+                        maxIterations(maxIter).
+                        checker((iteration, previous, current) -> current.getRMS() <= convergenceRMS).
                        start(fitted).
                        target(target).
                        weight(new DiagonalMatrix(weights)).

--- a/src/test/java/org/orekit/utils/SecularAndHarmonicTest.java
+++ b/src/test/java/org/orekit/utils/SecularAndHarmonicTest.java
@@ -28,6 +28,7 @@ import org.hipparchus.analysis.solvers.BaseUnivariateSolver;
 import org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver;
 import org.hipparchus.analysis.solvers.UnivariateSolverUtils;
 import org.hipparchus.exception.LocalizedCoreFormats;
+import org.hipparchus.exception.MathIllegalStateException;
 import org.hipparchus.exception.MathRuntimeException;
 import org.hipparchus.geometry.euclidean.threed.Vector3D;
 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
@@ -164,6 +165,60 @@ public class SecularAndHarmonicTest {

    }

+    @Test
+    public void testLinear() {
+        SecularAndHarmonic sh = new SecularAndHarmonic(1);
+        sh.setConvergenceRMS(1.0e-6);
+        sh.setMaxIter(5);
+        sh.resetFitting(AbsoluteDate.J2000_EPOCH, 0.0, 0.0);
+        sh.addPoint(sh.getReferenceDate().shiftedBy(1.0), 1.0);
+        sh.addPoint(sh.getReferenceDate().shiftedBy(2.0), 2.0);
+        sh.addPoint(sh.getReferenceDate().shiftedBy(3.0), 3.0);
+        sh.fit();
+        Assert.assertEquals(0.0, sh.getFittedParameters()[0], 1.0e-15);
+        Assert.assertEquals(1.0, sh.getFittedParameters()[1], 1.0e-15);
+    }
+
+    @Test
+    public void testImpossibleConvergence() {
+        SecularAndHarmonic sh = new SecularAndHarmonic(1);
+        sh.setConvergenceRMS(1.0e-6);
+        sh.setMaxIter(5);
+        sh.resetFitting(AbsoluteDate.J2000_EPOCH, 0.0, 0.0);
+        sh.addPoint(sh.getReferenceDate().shiftedBy(1.0), 1.0);
+        sh.addPoint(sh.getReferenceDate().shiftedBy(2.0), 2.0);
+        sh.addPoint(sh.getReferenceDate().shiftedBy(3.0), Double.NaN);
+        try {
+            sh.fit();
+            Assert.fail("an exception should have been thrown");
+        } catch (MathIllegalStateException mise) {
+            Assert.assertEquals(LocalizedCoreFormats.MAX_COUNT_EXCEEDED, mise.getSpecifier());
+        }
+    }
+
+    @Test
+    public void testConvergence() {
+        try {
+            doTestConvergence(0.2, 3);
+            Assert.fail("an exception should have been thrown");
+        } catch (MathIllegalStateException mise) {
+            Assert.assertEquals(LocalizedCoreFormats.MAX_COUNT_EXCEEDED, mise.getSpecifier());
+        }
+        doTestConvergence(0.3, 3);
+    }
+
+    private void doTestConvergence(final double rms, final int maxIter) {
+        SecularAndHarmonic sh = new SecularAndHarmonic(0, 1.0, 2.0, 3.0);
+        sh.setConvergenceRMS(rms);
+        sh.setMaxIter(maxIter);
+        sh.resetFitting(AbsoluteDate.J2000_EPOCH, new double[7]);
+        for (double t = -1.0; t <= 1.0; t += 0.125) {
+            // step function
+            sh.addPoint(sh.getReferenceDate().shiftedBy(t), FastMath.copySign(1.0, t));
+        }
+        sh.fit();
+    }
+
    @Test
    public void testPerfectOsculatingModel() {