This is an automated email from the ASF dual-hosted git repository. desruisseaux pushed a commit to branch geoapi-4.0 in repository https://gitbox.apache.org/repos/asf/sis.git
commit b3d55b4f4eda56fe7a6c3471f275188877bd001d Author: Martin Desruisseaux <[email protected]> AuthorDate: Sat Jan 16 00:38:45 2021 +0100 Add an `iterationReachedPrecisionLimit` flag in the test for deciding more accurately when assertion tolerance threshold should be relaxed. --- .../sis/referencing/GeodesicsOnEllipsoid.java | 3 ++ .../sis/referencing/GeodesicsOnEllipsoidTest.java | 54 ++++++++++++++++++++-- .../sis/referencing/GeodeticCalculatorTest.java | 49 ++++++++++++++++---- 3 files changed, 92 insertions(+), 14 deletions(-) diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodesicsOnEllipsoid.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodesicsOnEllipsoid.java index 6d36f7d..52fd130 100644 --- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodesicsOnEllipsoid.java +++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodesicsOnEllipsoid.java @@ -796,6 +796,9 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator { */ if (α1 == (α1 -= dα1)) { moreRefinements = 0; + if (STORE_LOCAL_VARIABLES) { + store("dα₁ ≪ α₁", dα1); // Flag for `iterationReachedPrecisionLimit` in tests. + } } if (STORE_LOCAL_VARIABLES) { // For comparing values against Karney table 5 and 6. final double I1_σ2; diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodesicsOnEllipsoidTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodesicsOnEllipsoidTest.java index 500ab82..c39d1ec 100644 --- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodesicsOnEllipsoidTest.java +++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodesicsOnEllipsoidTest.java @@ -42,16 +42,17 @@ import static org.apache.sis.internal.metadata.ReferencingServices.NAUTICAL_MILE * * @author Matthieu Bastianelli (Geomatys) * @author Martin Desruisseaux (Geomatys) - * @version 1.0 + * @version 1.1 * @since 1.0 * @module */ @DependsOn(GeodeticCalculatorTest.class) public final strictfp class GeodesicsOnEllipsoidTest extends GeodeticCalculatorTest { /** - * The instance to be tested. + * The {@link GeodesicsOnEllipsoid} instance to be tested. + * A specialized type is used for tracking locale variables. */ - private GeodesicsOnEllipsoid testedEarth; + private Calculator testedEarth; /** * Values of local variables in {@link GeodesicsOnEllipsoid} methods. If values for the same key are added @@ -92,10 +93,12 @@ public final strictfp class GeodesicsOnEllipsoidTest extends GeodeticCalculatorT private void createTracked() { localVariables = new HashMap<>(); testedEarth = new Calculator(HardCodedCRS.WGS84) { + /** Replaces a computed value by the value given in Karney table. */ @Override double computedToGiven(final double α1) { return (abs(TRUNCATED_α1 - toDegrees(α1)) < 1E-3) ? toRadians(TRUNCATED_α1) : α1; } + /** Invoked when {@link GeodesicsOnEllipsoid} computed an intermediate value. */ @Override void store(final String name, final double value) { super.store(name, value); if (verifyConsistency) { @@ -119,6 +122,15 @@ public final strictfp class GeodesicsOnEllipsoidTest extends GeodeticCalculatorT * by {@link MathFunctions#polynomialRoots(double...)}. */ private static class Calculator extends GeodesicsOnEllipsoid { + /** + * {@code true} if iteration stopped before to reach the desired accuracy because of limitation + * in {@code double} precision. This field must be reset to {@code false} before any new point. + * + * @see #iterationReachedPrecisionLimit() + * @see #clear() + */ + private boolean iterationReachedPrecisionLimit; + /** Values needed for computation of μ. */ private double x, y; @@ -129,6 +141,9 @@ public final strictfp class GeodesicsOnEllipsoidTest extends GeodeticCalculatorT /** Invoked when {@link GeodesicsOnEllipsoid} computed an intermediate value. */ @Override void store(final String name, final double value) { + if (name.equals("dα₁ ≪ α₁")) { + iterationReachedPrecisionLimit = true; + } if (name.length() == 1) { switch (name.charAt(0)) { case 'x': x = value; break; @@ -154,6 +169,21 @@ public final strictfp class GeodesicsOnEllipsoidTest extends GeodeticCalculatorT } /** + * Clears the tested {@link GeodeticCalculator} before to test a new point. + * This is invoked by parent class between two tests using the same calculator. + * The intent is to make sure that data from previous test are not mixed with current test. + */ + @Override + void clear() { + if (localVariables != null) { + localVariables.clear(); + } + testedEarth.x = Double.NaN; + testedEarth.y = Double.NaN; + testedEarth.iterationReachedPrecisionLimit = false; + } + + /** * Asserts that variable of the given name is equal to the given value. This is used for comparing an * intermediate value computed by Apache SIS against the value published in a table of Karney (2013) * <cite>Algorithms for geodesics</cite>. @@ -493,6 +523,22 @@ public final strictfp class GeodesicsOnEllipsoidTest extends GeodeticCalculatorT } /** + * Returns {@code true} if iteration stopped before to reach the desired accuracy because of limitation in + * {@code double} precision. This problem may happen in the {@link GeodesicsOnEllipsoid#computeDistance()} + * method when {@literal dα₁ ≪ α₁}. If locale variable storage is enabled, this situation is flagged by the + * {@code "dα₁ ≪ α₁"} key. Otherwise we conservatively assume that this situation occurred. + */ + @Override + boolean iterationReachedPrecisionLimit() { + if (GeodesicsOnEllipsoid.STORE_LOCAL_VARIABLES) { + return testedEarth.iterationReachedPrecisionLimit; + } else { + // Conservative value in absence of information. + return true; + } + } + + /** * Tells whether failure to compute geodesic for the given data should cause the test case to fail. * This is invoked by {@link #compareAgainstDataset()}. * @@ -505,7 +551,7 @@ public final strictfp class GeodesicsOnEllipsoidTest extends GeodeticCalculatorT final double φ2 = expected[COLUMN_φ2]; final double Δλ = expected[COLUMN_λ2] - expected[COLUMN_λ1]; if (Δλ > 90 && max(abs(φ1), abs(φ2)) < 2E-4) { - return false; // Ignore equatorial case. + return true; // Ignore equatorial case in previous version (not anymore). } if (Δλ > 179 && abs(φ1 + φ2) < 0.002) { return false; // Ignore antipodal case. diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java index 96933cd..1df1bb6 100644 --- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java +++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java @@ -120,6 +120,15 @@ public strictfp class GeodeticCalculatorTest extends TestCase { } /** + * Clears test data before to test a new point. The default implementation does nothing. + * This is overridden when specialized {@link GeodesicsOnEllipsoid} subclasses are used + * for tracking local variables, for making sure that we do not mix the variables of two + * different test cases. + */ + void clear() { + } + + /** * Tests some simple azimuth directions. The expected directions are approximately North, East, * South and West, but not exactly because of Earth curvature. The test verify merely that the * azimuths are approximately correct. @@ -150,6 +159,7 @@ public strictfp class GeodeticCalculatorTest extends TestCase { c.setEndGeographicPoint (-20, 30); assertEquals( 180, c.getStartingAzimuth(), tolerance); c.setEndGeographicPoint (-90, 30); assertEquals( 180, c.getStartingAzimuth(), tolerance); + clear(); c.setStartGeographicPoint( 90, 0); c.setEndGeographicPoint ( 20, 20); assertEquals( 160, c.getStartingAzimuth(), tolerance); c.setEndGeographicPoint ( 20, -20); assertEquals(-160, c.getStartingAzimuth(), tolerance); @@ -367,6 +377,7 @@ public strictfp class GeodeticCalculatorTest extends TestCase { } } } + clear(); } if (errors != null) { out.println("Distance between points on Bézier curve and points on geodesic."); @@ -456,12 +467,14 @@ public strictfp class GeodeticCalculatorTest extends TestCase { final double cosφ1 = abs(cos(toRadians(expected[COLUMN_φ1]))); // For adjusting longitude tolerance. final double cosφ2 = abs(cos(toRadians(expected[COLUMN_φ2]))); double linearTolerance, latitudeTolerance, longitudeTolerance, azimuthTolerance; + final boolean isLongArcOnEquator; if (isSphere) { /* * When spherical formulas are used instead than ellipsoidal formulas, an error up to 1% is expected * in distance calculations (source: Wikipedia). A yet larger error is observed for azimuth values, * especially near poles and between antipodal points. Following are empirical thresholds. */ + isLongArcOnEquator = false; // Not a problem in spherical case. linearTolerance = expected[COLUMN_Δs] * 0.01; latitudeTolerance = toDegrees(linearTolerance / c.semiMajorAxis); longitudeTolerance = expected[COLUMN_φ2] > 89.5 ? 180 : latitudeTolerance / cosφ2; @@ -492,13 +505,12 @@ public strictfp class GeodeticCalculatorTest extends TestCase { } } /* - * If the start and end points are both within about 3 meters from equator and their - * distance is more than 10,000 km, relax the tolerance to 2 millimetres instead of 1. + * If the start and end points are both within about 3 meters from equator and + * their distance is more than 10,000 km, we may need to relax the tolerance. * This is an empirical adjustment based on test failures observation. */ - if (Math.min(cosφ1, cosφ2) >= 0.9999999999999 && expected[COLUMN_Δs] > 1E+7) { - linearTolerance *= 2; - } + isLongArcOnEquator = Math.min(cosφ1, cosφ2) >= 0.9999999999999 && expected[COLUMN_Δs] > 1E+7; + assertEquals("Consistency with accuracy reported in Javadoc.", 0.001, linearTolerance, 0.0005); } /* * Set input values, compute then verify results. The azimuth tolerance is divided by cos(φ). @@ -523,7 +535,9 @@ public strictfp class GeodeticCalculatorTest extends TestCase { assertEquals("φ₂", expected[COLUMN_φ2], end.getOrdinate(0), latitudeTolerance); assertEquals("λ₂", expected[COLUMN_λ2], end.getOrdinate(1), longitudeTolerance); assertEquals("α₂", expected[COLUMN_α2], c.getEndingAzimuth(), azimuthTolerance / cosφ2); - assertEquals("∆s", expected[COLUMN_Δs], c.getGeodesicDistance(), linearTolerance); + assertEquals("∆s", expected[COLUMN_Δs], c.getGeodesicDistance(), linearTolerance * + (isLongArcOnEquator && iterationReachedPrecisionLimit() ? 2 : 1)); + clear(); } catch (GeodeticException | AssertionError e) { if (!isTestingInverse || e instanceof AssertionError || isFailure(expected)) { out.printf("Test failure at line %d: %s%n" @@ -544,6 +558,21 @@ public strictfp class GeodeticCalculatorTest extends TestCase { } /** + * Returns {@code true} if iteration stopped before to reach the desired accuracy because the correction + * to apply in iterative steps is smaller than what can be applied using IEEE 754 double arithmetic. + * It never happen in the spherical case because there is no iteration in that case. + * + * <p>This method is invoked only in situations empirically identified as susceptible to have this problem:</p> + * <uL> + * <li>When start point and end point are both within about 3 meters from equator + * and their distance is more than 10,000 km.</li> + * </ul> + */ + boolean iterationReachedPrecisionLimit() { + return false; + } + + /** * Tells whether failure to compute geodesic for the given data should cause the test case to fail. * The default implementation always return {@code true}. Subclass can override if some points are * known to fail. @@ -557,7 +586,7 @@ public strictfp class GeodeticCalculatorTest extends TestCase { } /** - * Tests {@link GeodesicsOnEllipsoid#getRhumblineLength()} using points given by Bennett (1996). + * Tests {@link GeodeticCalculator#getRhumblineLength()} using points given by Bennett (1996). * This is an anti-regression test since the result was computed by SIS for the spherical case. */ @Test @@ -570,7 +599,7 @@ public strictfp class GeodeticCalculatorTest extends TestCase { } /** - * Tests {@link GeodesicsOnEllipsoid#getRhumblineLength()} using points given by Bennett (1996). + * Tests {@link GeodeticCalculator#getRhumblineLength()} using points given by Bennett (1996). * This is an anti-regression test since the result was computed by SIS for the spherical case. */ @Test @@ -583,7 +612,7 @@ public strictfp class GeodeticCalculatorTest extends TestCase { } /** - * Tests {@link GeodesicsOnEllipsoid#getRhumblineLength()} using points given by Bennett (1996). + * Tests {@link GeodeticCalculator#getRhumblineLength()} using points given by Bennett (1996). * This is an anti-regression test since the result was computed by SIS for the spherical case. */ @Test @@ -602,7 +631,7 @@ public strictfp class GeodeticCalculatorTest extends TestCase { * @throws TransformException if an error occurred while projection the test point. */ @Test - public void test() throws TransformException { + public void testProjectionAroundStart() throws TransformException { final GeodeticCalculator c = create(false); final double distance = 600000; // In metres. final double azimuth = 37; // Geographic angle (degrees relative to North).
