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 86d5499ccdbbc94470f5045e2c1b0b07a3b37b06..78f406bd052e42911ffd080795b5090306347aad 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 4660fba74c0586128b61c1c9dcb84a9d8c1cbe73..a93136ac06221d8c73473a4d957964f18b88d154 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 adef7fa7ec77bd45f1f0c2108653d8934dcd4524..2feba848f73bb31702fb8b6fc5e706876af8d41b 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 81a11b83d94c0633043bd956772e55e6ad7518d7..26691ecf5adea1ebad2b47b32aa6b9df9376484e 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 e295da84a6e48c29629b757df0678a7ead5a3365..1cbd11ed9acea14023d504dea62ec82bbd89b3a2 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 6b6b8f2be9a3e148cc5d7ec3910fe88238e2855a..e8031f79b02d0f74d5283454fb9cdad684a4289f 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 1acbd1cb5ccae0b85db64954a6a7ff7514515ef4..03006374887772c62f10ecad823a8e1a77bd22f5 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 e0d9f914125690b1ebd131774adffbc766078b20..609dacaefabdb36d88a7d3c3fc2e17c1af0bb3aa 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 09e8f49e62400ed1da052f95d5fc87f71386437f..b6ae58f83a221f924615dc7c86d20775a0731216 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 696acc7d5c5e07b626827e206086915b629b31ec..48a3eaefa83eabc6fc1ba98127d2b091a5de00b3 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 7352c89a906fd24d4c62fd19754e2c5ae59815f2..ce01c5bed557b29b69995ed41b7ff28024b87ac0 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 90a11e2080bec61cd1f39994f8d1a5a0d87e739f..7b2112c39ec3461e8a73727b4a7414fbe6b2b89b 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 d3dbfddc360858b3524a63c6cb8cad761aa66613..12394e04f9dba05618ea2c5983cb8125121c7e27 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 bedb295e633f80f093421ae3012b69d03c604be7..c6ad042e2869ef3bdf0b90d0ff2b5976d03bcba9 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 9db9c86f2781b0f4bf54e15ce6caf371a625dfe0..54528e60a6c0c8a4e79a1f5e7c33b0773464cac4 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 1f020713127df122a8f94addca45a960147c76a2..1495fd70229ac13c3ab63a8bf180d66819be8b30 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 0000000000000000000000000000000000000000..62b7e83c2208457d92ca0ebd32bf3cbd28a4b8b3 --- /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); + } + + } + +} +