From 64172568144a6f4d5bdd280e9550473838f6b5b8 Mon Sep 17 00:00:00 2001
From: Luc Maisonobe <luc@orekit.org>
Date: Tue, 8 Jan 2019 18:20:26 +0100
Subject: [PATCH] Fixed scaling error in ParameterFunction differentiation.
Fixes #509
---
.../RangeIonosphericDelayModifier.java | 2 +-
.../RangeRateIonosphericDelayModifier.java | 2 +-
.../RangeRateTroposphericDelayModifier.java | 2 +-
.../RangeTroposphericDelayModifier.java | 2 +-
...rnAroundRangeIonosphericDelayModifier.java | 2 +-
...nAroundRangeTroposphericDelayModifier.java | 2 +-
.../org/orekit/utils/Differentiation.java | 64 +++++++++++--------
src/site/xdoc/changes.xml | 3 +
.../measurements/AngularAzElTest.java | 2 +-
.../measurements/AngularRaDecTest.java | 2 +-
.../estimation/measurements/PhaseTest.java | 2 +-
.../estimation/measurements/Range2Test.java | 2 +-
.../measurements/RangeRateTest.java | 6 +-
.../estimation/measurements/RangeTest.java | 2 +-
.../TurnAroundRangeAnalyticTest.java | 2 +-
.../measurements/TurnAroundRangeTest.java | 2 +-
.../org/orekit/utils/DifferentiationTest.java | 60 +++++++++++++++++
17 files changed, 115 insertions(+), 44 deletions(-)
create mode 100644 src/test/java/org/orekit/utils/DifferentiationTest.java
diff --git a/src/main/java/org/orekit/estimation/measurements/modifiers/RangeIonosphericDelayModifier.java b/src/main/java/org/orekit/estimation/measurements/modifiers/RangeIonosphericDelayModifier.java
index 86d5499cc..78f406bd0 100644
--- a/src/main/java/org/orekit/estimation/measurements/modifiers/RangeIonosphericDelayModifier.java
+++ b/src/main/java/org/orekit/estimation/measurements/modifiers/RangeIonosphericDelayModifier.java
@@ -140,7 +140,7 @@ public class RangeIonosphericDelayModifier implements EstimationModifier<Range>
};
final ParameterFunction rangeErrorDerivative =
- Differentiation.differentiate(rangeError, driver, 3, 10.0);
+ Differentiation.differentiate(rangeError, 3, 10.0 * driver.getScale());
return rangeErrorDerivative.value(driver);
diff --git a/src/main/java/org/orekit/estimation/measurements/modifiers/RangeRateIonosphericDelayModifier.java b/src/main/java/org/orekit/estimation/measurements/modifiers/RangeRateIonosphericDelayModifier.java
index 4660fba74..a93136ac0 100644
--- a/src/main/java/org/orekit/estimation/measurements/modifiers/RangeRateIonosphericDelayModifier.java
+++ b/src/main/java/org/orekit/estimation/measurements/modifiers/RangeRateIonosphericDelayModifier.java
@@ -175,7 +175,7 @@ public class RangeRateIonosphericDelayModifier implements EstimationModifier<Ran
};
final ParameterFunction rangeErrorDerivative =
- Differentiation.differentiate(rangeError, driver, 3, 10.0);
+ Differentiation.differentiate(rangeError, 3, 10.0 * driver.getScale());
return rangeErrorDerivative.value(driver);
diff --git a/src/main/java/org/orekit/estimation/measurements/modifiers/RangeRateTroposphericDelayModifier.java b/src/main/java/org/orekit/estimation/measurements/modifiers/RangeRateTroposphericDelayModifier.java
index adef7fa7e..2feba848f 100644
--- a/src/main/java/org/orekit/estimation/measurements/modifiers/RangeRateTroposphericDelayModifier.java
+++ b/src/main/java/org/orekit/estimation/measurements/modifiers/RangeRateTroposphericDelayModifier.java
@@ -174,7 +174,7 @@ public class RangeRateTroposphericDelayModifier implements EstimationModifier<Ra
};
final ParameterFunction rangeErrorDerivative =
- Differentiation.differentiate(rangeError, driver, 3, 10.0);
+ Differentiation.differentiate(rangeError, 3, 10.0 * driver.getScale());
return rangeErrorDerivative.value(driver);
diff --git a/src/main/java/org/orekit/estimation/measurements/modifiers/RangeTroposphericDelayModifier.java b/src/main/java/org/orekit/estimation/measurements/modifiers/RangeTroposphericDelayModifier.java
index 81a11b83d..26691ecf5 100644
--- a/src/main/java/org/orekit/estimation/measurements/modifiers/RangeTroposphericDelayModifier.java
+++ b/src/main/java/org/orekit/estimation/measurements/modifiers/RangeTroposphericDelayModifier.java
@@ -143,7 +143,7 @@ public class RangeTroposphericDelayModifier implements EstimationModifier<Range>
};
final ParameterFunction rangeErrorDerivative =
- Differentiation.differentiate(rangeError, driver, 3, 10.0);
+ Differentiation.differentiate(rangeError, 3, 10.0 * driver.getScale());
return rangeErrorDerivative.value(driver);
diff --git a/src/main/java/org/orekit/estimation/measurements/modifiers/TurnAroundRangeIonosphericDelayModifier.java b/src/main/java/org/orekit/estimation/measurements/modifiers/TurnAroundRangeIonosphericDelayModifier.java
index e295da84a..1cbd11ed9 100644
--- a/src/main/java/org/orekit/estimation/measurements/modifiers/TurnAroundRangeIonosphericDelayModifier.java
+++ b/src/main/java/org/orekit/estimation/measurements/modifiers/TurnAroundRangeIonosphericDelayModifier.java
@@ -138,7 +138,7 @@ public class TurnAroundRangeIonosphericDelayModifier implements EstimationModifi
};
final ParameterFunction rangeErrorDerivative =
- Differentiation.differentiate(rangeError, driver, 3, 10.0);
+ Differentiation.differentiate(rangeError, 3, 10.0 * driver.getScale());
return rangeErrorDerivative.value(driver);
diff --git a/src/main/java/org/orekit/estimation/measurements/modifiers/TurnAroundRangeTroposphericDelayModifier.java b/src/main/java/org/orekit/estimation/measurements/modifiers/TurnAroundRangeTroposphericDelayModifier.java
index 6b6b8f2be..e8031f79b 100644
--- a/src/main/java/org/orekit/estimation/measurements/modifiers/TurnAroundRangeTroposphericDelayModifier.java
+++ b/src/main/java/org/orekit/estimation/measurements/modifiers/TurnAroundRangeTroposphericDelayModifier.java
@@ -139,7 +139,7 @@ public class TurnAroundRangeTroposphericDelayModifier implements EstimationModif
}
};
- final ParameterFunction rangeErrorDerivative = Differentiation.differentiate(rangeError, driver, 3, 10.0);
+ final ParameterFunction rangeErrorDerivative = Differentiation.differentiate(rangeError, 3, 10.0 * driver.getScale());
return rangeErrorDerivative.value(driver);
diff --git a/src/main/java/org/orekit/utils/Differentiation.java b/src/main/java/org/orekit/utils/Differentiation.java
index 1acbd1cb5..030063748 100644
--- a/src/main/java/org/orekit/utils/Differentiation.java
+++ b/src/main/java/org/orekit/utils/Differentiation.java
@@ -21,11 +21,8 @@ import org.hipparchus.analysis.UnivariateVectorFunction;
import org.hipparchus.analysis.differentiation.DSFactory;
import org.hipparchus.analysis.differentiation.DerivativeStructure;
import org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator;
-import org.hipparchus.analysis.differentiation.UnivariateDifferentiableFunction;
import org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction;
import org.orekit.attitudes.AttitudeProvider;
-import org.orekit.errors.OrekitException;
-import org.orekit.errors.OrekitMessages;
import org.orekit.orbits.Orbit;
import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngle;
@@ -50,41 +47,52 @@ public class Differentiation {
* @param function function to differentiate
* @param driver driver for the parameter
* @param nbPoints number of points used for finite differences
- * @param step step for finite differences
+ * @param normalizedStep step for finite differences, in <em>normalized</em> units
* @return scalar function evaluating to the derivative of the original function
+ * @deprecated as of 9.3, replaced by {@link #differentiate(ParameterFunction, int, double)}
+ */
+ @Deprecated
+ public static ParameterFunction differentiate(final ParameterFunction function, final ParameterDriver driver,
+ final int nbPoints, final double normalizedStep) {
+ return differentiate(function, nbPoints, normalizedStep * driver.getScale());
+ }
+
+ /** Differentiate a scalar function using finite differences.
+ * @param function function to differentiate
+ * @param nbPoints number of points used for finite differences
+ * @param step step for finite differences, in <em>physical</em> units
+ * @return scalar function evaluating to the derivative of the original function
+ * @since 9.3
*/
public static ParameterFunction differentiate(final ParameterFunction function,
- final ParameterDriver driver,
final int nbPoints, final double step) {
- final UnivariateFunction uf = new UnivariateFunction() {
- /** {@inheritDoc} */
- @Override
- public double value(final double normalizedValue) {
- final double saved = driver.getNormalizedValue();
- driver.setNormalizedValue(normalizedValue);
- final double functionValue = function.value(driver);
- driver.setNormalizedValue(saved);
- return functionValue;
- }
- };
+ return new ParameterFunction() {
- final FiniteDifferencesDifferentiator differentiator =
- new FiniteDifferencesDifferentiator(nbPoints, step);
- final UnivariateDifferentiableFunction differentiated =
- differentiator.differentiate(uf);
+ /** Finite differences differentiator to use. */
+ private final FiniteDifferencesDifferentiator differentiator =
+ new FiniteDifferencesDifferentiator(nbPoints, step);
- return new ParameterFunction() {
/** {@inheritDoc} */
@Override
- public double value(final ParameterDriver parameterDriver) {
- if (!parameterDriver.getName().equals(driver.getName())) {
- throw new OrekitException(OrekitMessages.UNSUPPORTED_PARAMETER_NAME,
- parameterDriver.getName(), driver.getName());
- }
- final DerivativeStructure dsParam = FACTORY.variable(0, parameterDriver.getNormalizedValue());
- final DerivativeStructure dsValue = differentiated.value(dsParam);
+ public double value(final ParameterDriver driver) {
+
+ final UnivariateFunction uf = new UnivariateFunction() {
+ /** {@inheritDoc} */
+ @Override
+ public double value(final double value) {
+ final double saved = driver.getValue();
+ driver.setValue(value);
+ final double functionValue = function.value(driver);
+ driver.setValue(saved);
+ return functionValue;
+ }
+ };
+
+ final DerivativeStructure dsParam = FACTORY.variable(0, driver.getValue());
+ final DerivativeStructure dsValue = differentiator.differentiate(uf).value(dsParam);
return dsValue.getPartialDerivative(1);
+
}
};
diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml
index e0d9f9141..609dacaef 100644
--- a/src/site/xdoc/changes.xml
+++ b/src/site/xdoc/changes.xml
@@ -21,6 +21,9 @@
</properties>
<body>
<release version="TBD" date="TBD" description="TBD">
+ <action dev="luc" type="fix" issue="509">
+ Fixed scaling error in ParameterFunction differentiation.
+ </action>
<action dev="luc" type="fix" issue="508">
Fixed inconsistency leading to inaccuracies in conversions from AbsoluteDate to FieldAbsoluteDate.
</action>
diff --git a/src/test/java/org/orekit/estimation/measurements/AngularAzElTest.java b/src/test/java/org/orekit/estimation/measurements/AngularAzElTest.java
index 09e8f49e6..b6ae58f83 100644
--- a/src/test/java/org/orekit/estimation/measurements/AngularAzElTest.java
+++ b/src/test/java/org/orekit/estimation/measurements/AngularAzElTest.java
@@ -254,7 +254,7 @@ public class AngularAzElTest {
public double value(final ParameterDriver parameterDriver) {
return measurement.estimate(0, 0, new SpacecraftState[] { state }).getEstimatedValue()[k];
}
- }, drivers[i], 3, 50.0);
+ }, 3, 50.0 * drivers[i].getScale());
final double ref = dMkdP.value(drivers[i]);
if (ref > 1.e-12) {
diff --git a/src/test/java/org/orekit/estimation/measurements/AngularRaDecTest.java b/src/test/java/org/orekit/estimation/measurements/AngularRaDecTest.java
index 696acc7d5..48a3eaefa 100644
--- a/src/test/java/org/orekit/estimation/measurements/AngularRaDecTest.java
+++ b/src/test/java/org/orekit/estimation/measurements/AngularRaDecTest.java
@@ -254,7 +254,7 @@ public class AngularRaDecTest {
public double value(final ParameterDriver parameterDriver) {
return measurement.estimate(0, 0, new SpacecraftState[] { state }).getEstimatedValue()[k];
}
- }, drivers[i], 3, 50.0);
+ }, 3, 50.0 * drivers[i].getScale());
final double ref = dMkdP.value(drivers[i]);
if (ref > 1.e-12) {
diff --git a/src/test/java/org/orekit/estimation/measurements/PhaseTest.java b/src/test/java/org/orekit/estimation/measurements/PhaseTest.java
index 7352c89a9..ce01c5bed 100644
--- a/src/test/java/org/orekit/estimation/measurements/PhaseTest.java
+++ b/src/test/java/org/orekit/estimation/measurements/PhaseTest.java
@@ -479,7 +479,7 @@ public class PhaseTest {
public double value(final ParameterDriver parameterDriver) {
return measurement.estimate(0, 0, new SpacecraftState[] { state }).getEstimatedValue()[0];
}
- }, drivers[i], 3, 20.0);
+ }, 3, 20.0 * drivers[i].getScale());
final double ref = dMkdP.value(drivers[i]);
if (printResults) {
diff --git a/src/test/java/org/orekit/estimation/measurements/Range2Test.java b/src/test/java/org/orekit/estimation/measurements/Range2Test.java
index 90a11e208..7b2112c39 100644
--- a/src/test/java/org/orekit/estimation/measurements/Range2Test.java
+++ b/src/test/java/org/orekit/estimation/measurements/Range2Test.java
@@ -522,7 +522,7 @@ public class Range2Test {
public double value(final ParameterDriver parameterDriver) {
return measurement.estimate(0, 0, new SpacecraftState[] { state }).getEstimatedValue()[0];
}
- }, drivers[i], 3, 20.0);
+ }, 3, 20.0 * drivers[i].getScale());
final double ref = dMkdP.value(drivers[i]);
if (printResults) {
diff --git a/src/test/java/org/orekit/estimation/measurements/RangeRateTest.java b/src/test/java/org/orekit/estimation/measurements/RangeRateTest.java
index d3dbfddc3..12394e04f 100644
--- a/src/test/java/org/orekit/estimation/measurements/RangeRateTest.java
+++ b/src/test/java/org/orekit/estimation/measurements/RangeRateTest.java
@@ -312,7 +312,7 @@ public class RangeRateTest {
public double value(final ParameterDriver parameterDriver) {
return measurement.estimate(0, 0, new SpacecraftState[] { state }).getEstimatedValue()[0];
}
- }, drivers[i], 3, 20.0);
+ }, 3, 20.0 * drivers[i].getScale());
final double ref = dMkdP.value(drivers[i]);
maxRelativeError = FastMath.max(maxRelativeError, FastMath.abs((ref - gradient[0]) / ref));
}
@@ -382,7 +382,7 @@ public class RangeRateTest {
public double value(final ParameterDriver parameterDriver) {
return measurement.estimate(0, 0, new SpacecraftState[] { state }).getEstimatedValue()[0];
}
- }, drivers[i], 3, 20.0);
+ }, 3, 20.0 * drivers[i].getScale());
final double ref = dMkdP.value(drivers[i]);
maxRelativeError = FastMath.max(maxRelativeError, FastMath.abs((ref - gradient[0]) / ref));
}
@@ -516,7 +516,7 @@ public class RangeRateTest {
public double value(final ParameterDriver parameterDriver) {
return measurement.estimate(0, 0, new SpacecraftState[] { state }).getEstimatedValue()[0];
}
- }, drivers[i], 3, 20.0);
+ }, 3, 20.0 * drivers[i].getScale());
final double ref = dMkdP.value(drivers[i]);
maxRelativeError = FastMath.max(maxRelativeError, FastMath.abs((ref - gradient[0]) / ref));
}
diff --git a/src/test/java/org/orekit/estimation/measurements/RangeTest.java b/src/test/java/org/orekit/estimation/measurements/RangeTest.java
index bedb295e6..c6ad042e2 100644
--- a/src/test/java/org/orekit/estimation/measurements/RangeTest.java
+++ b/src/test/java/org/orekit/estimation/measurements/RangeTest.java
@@ -513,7 +513,7 @@ public class RangeTest {
public double value(final ParameterDriver parameterDriver) {
return measurement.estimate(0, 0, new SpacecraftState[] { state }).getEstimatedValue()[0];
}
- }, drivers[i], 3, 20.0);
+ }, 3, 20.0 * drivers[i].getScale());
final double ref = dMkdP.value(drivers[i]);
if (printResults) {
diff --git a/src/test/java/org/orekit/estimation/measurements/TurnAroundRangeAnalyticTest.java b/src/test/java/org/orekit/estimation/measurements/TurnAroundRangeAnalyticTest.java
index 9db9c86f2..54528e60a 100644
--- a/src/test/java/org/orekit/estimation/measurements/TurnAroundRangeAnalyticTest.java
+++ b/src/test/java/org/orekit/estimation/measurements/TurnAroundRangeAnalyticTest.java
@@ -591,7 +591,7 @@ public class TurnAroundRangeAnalyticTest {
public double value(final ParameterDriver parameterDriver) {
return measurement.estimate(0, 0, new SpacecraftState[] { state }).getEstimatedValue()[0];
}
- }, drivers[i], 3, 20.0);
+ }, 3, 20.0 * drivers[i].getScale());
ref = dMkdP.value(drivers[i]);
} else {
// Compute a reference value using TurnAroundRange function
diff --git a/src/test/java/org/orekit/estimation/measurements/TurnAroundRangeTest.java b/src/test/java/org/orekit/estimation/measurements/TurnAroundRangeTest.java
index 1f0207131..1495fd702 100644
--- a/src/test/java/org/orekit/estimation/measurements/TurnAroundRangeTest.java
+++ b/src/test/java/org/orekit/estimation/measurements/TurnAroundRangeTest.java
@@ -502,7 +502,7 @@ public class TurnAroundRangeTest {
public double value(final ParameterDriver parameterDriver) {
return measurement.estimate(0, 0, new SpacecraftState[] { state }).getEstimatedValue()[0];
}
- }, drivers[i], 3, 20.0);
+ }, 3, 20.0 * drivers[i].getScale());
final double ref = dMkdP.value(drivers[i]);
// Deltas
diff --git a/src/test/java/org/orekit/utils/DifferentiationTest.java b/src/test/java/org/orekit/utils/DifferentiationTest.java
new file mode 100644
index 000000000..62b7e83c2
--- /dev/null
+++ b/src/test/java/org/orekit/utils/DifferentiationTest.java
@@ -0,0 +1,60 @@
+/* Copyright 2002-2019 CS Systèmes d'Information
+ * Licensed to CS Systèmes d'Information (CS) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * CS licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.orekit.utils;
+
+
+import org.hipparchus.util.FastMath;
+import org.hipparchus.util.Precision;
+import org.junit.Assert;
+import org.junit.Test;
+
+public class DifferentiationTest {
+
+ @Test
+ public void testScaleOne() {
+ // with this step, computation is exact in IEEE754
+ doTestScale(1.0, FastMath.pow(1.0, -3), Precision.SAFE_MIN);
+ }
+
+ @Test
+ public void testScalePowerOfTwoStepRepresentableNumber() {
+ // with this step, computation is exact in IEEE754
+ doTestScale(FastMath.scalb(1.0, -10), FastMath.pow(1.0, -7), Precision.SAFE_MIN);
+ }
+
+ @Test
+ public void testScalePowerOfTwoStepNonRepresentableNumber() {
+ // with this step, computation has numerical noise
+ doTestScale(FastMath.scalb(1.0, -10), 0.007, 1.7e-12);
+ }
+
+ private void doTestScale(final double scale, final double step, final double tolerance) {
+ ParameterDriver driver = new ParameterDriver("", -100.0, scale,
+ Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY);
+ ParameterFunction f0 = d -> 3 * d.getValue() * d.getValue() - 2 * d.getValue();
+ ParameterFunction f1Diff = Differentiation.differentiate(f0, 4, step);
+ ParameterFunction f1Ref = d -> 6 * d.getValue() - 2;
+
+ for (double x = -3.0; x < 3.0; x += 0.125) {
+ driver.setValue(x);
+ Assert.assertEquals(f1Ref.value(driver), f1Diff.value(driver), tolerance);
+ }
+
+ }
+
+}
+
--
GitLab