Rework numerical tolerances in test suite

[ettersi]

trafficstars

This PR reworks the numerical tolerances in the test suite as suggested in https://github.com/JuliaEarth/CoordRefSystems.jl/pull/240#discussion_r1925135197 and discussed further in https://github.com/JuliaEarth/CoordRefSystems.jl/issues/242.

There's two main changes here.

1) Redefine allapprox()

The test suite relies mainly on allapprox(c1, c2) to check for numerical accuracy. Currently, this function checks that c1.X ≈ c2.X for every component X of the coordinate type in question.

https://github.com/JuliaEarth/CoordRefSystems.jl/blob/8380fe8d5ed0e92c66352bb080ad9c2a6d888bec/test/testutils.jl#L1-L18

Such a component-wise accuracy test isn't particularly meaningful. For example, under the above definition we have

julia> allapprox(LatLon(0.0, 1.0), LatLon(0.0, 1.0 + eps()))
true

but

julia> allapprox(LatLon(0.0, 0.0), LatLon(0.0, eps()))
false

even though the only difference between the two comparisons is an arbitrary choice of where the prime meridian should be.

This PR redefines allapprox() as follows.

• For Cartesian coordinates, do a vector-wise, relative comparison.

  svec(coords::Cartesian) = SVector(getfield(coords, :coords))
  allapprox(coords₁::Cartesian{Datum}, coords₂::Cartesian{Datum}) where {Datum} = svec(coords₁) ≈ svec(coords₂)
  

• For radial coordinates like Polar, Cylindrical or LatLon, do a relative comparison for the length components and an absolute comparison for the angular components. For example:

  function is_approx_angle(α, β)
      # Check that abs(α - β) is small relative to and modulo 2π / 180°.
      # Details are a bit messy, so I'm omitting them here. 
  end
  allapprox(coords₁::Polar{Datum}, coords₂::Polar{Datum}) where {Datum} = coords₁.ρ ≈ coords₂.ρ && is_approx_angle(coords₁.ϕ, coords₂.ϕ)
  

• For Projected CRS, do a relative comparison lower-bounded by an absolute comparison with the radius of the earth as reference.

  svec(coords::CoordRefSystems.Projected) = SVector(coords.x, coords.y)
  function allapprox(coords₁::C, coords₂::C) where {C <: CoordRefSystems.Projected}
    T = promote_type(Unitful.numtype.((coords₁.x, coords₂.x))...)
    return isapprox(
      svec(coords₁),
      svec(coords₂);
      rtol = sqrt(eps(T)),
      atol = sqrt(eps(T(ustrip(m, majoraxis(ellipsoid(datum(C))))))) * m,
    )
  end
  

Honestly, I can't really explain why exactly this is the right thing to do - as far as tangible explanations go, I think the best I can offer is to point out that by doing this, I managed to simplify and streamline several test cases. More abstractly speaking, I guess the point is that these new definitions reflect the underlying symmetries (scale and rotation) better than the old ones did, but this is fairly hand-wavy explanation so make of it what you will.

2) Rework the LatLon round-trip tests

The above redefinition of allapprox() broke some of the tests in test/crs/domains.jl, so I rewrote that file quite a bit. As you will have noticed, this work uncovered several accuracy bugs that I've filed in the GitHub issues mentioned in the source code comments.

One thing that I should probably explain in a bit more detail here is the story behind the new sqrt_tol() function in domains.jl.

sqrt_tol(x, xlim) = sqrt_tol(abs(x - xlim) / xlim) * xlim
function sqrt_tol(e)
  T = typeof(e)
  if e >= eps(T)^(1//4)
    return eps(T)^(1//2)
  elseif e >= eps(T)^(1//2)
    return eps(T)^(3//4) / e
  else
    return eps(T)^(1//4)
  end
end

The fundamental problem addressed by this function is that when doing arithmetic on a circle / sphere, you sometimes encounter operations that fundamentally can be accurate only up to O(sqrt(eps())), not O(eps()). The prototypical example of such an operation is the function f(x) = sqrt(1 - x^2) which maps an x-coordinate to the y-coordinate f(x) such that [x, f(x)] lies on the unit circle. The problem with this function is that it maps all the points in the interval x in 0 .. sqrt(eps())/2 to the image f(x) in 1-eps()/4 .. 1 which contains only one floating-point number (namely 1.0). Therefore, all these inputs are mapped to the same output,

julia> x = LinRange(0, sqrt(eps())/2, 100_000)
       unique(f.(x))
1-element Vector{Float64}:
 1.0

and so when you apply the inverse function (which in this case happens to be the same function f(x)), you can't recover the information that was lost in the forward pass, resulting in an error that can be up to sqrt(eps())/2:

julia> x = sqrt(eps())/2
       abs(f(f(x)) - x) == sqrt(eps())/2
true

Situations analogous to this arise in the round-trip conversions between LatLon and several Projected coordinates, and in these cases I'm using sqrt_tol() to compute relaxed tolerances as we approach the pathological limits.

---

Beyond these to main changes, I've removed a couple of if T == Float64 checks in test/crs/conversions.jl since with the new definition of allapprox(), these tests now pass also for T == Float32. It's the fact that I could do this which in my opinion serves as a strong justification that the new definitions are "better" than the old ones.

[codecov-commenter]

Codecov Report

All modified and coverable lines are covered by tests :white_check_mark:

Project coverage is 93.03%. Comparing base (0dc546e) to head (a45d558).

Additional details and impacted files

@@           Coverage Diff           @@
##             main     #269   +/-   ##
=======================================
  Coverage   93.03%   93.03%           
=======================================
  Files          35       35           
  Lines        1652     1652           
=======================================
  Hits         1537     1537           
  Misses        115      115           

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

:rocket: New features to boost your workflow:

• :snowflake: Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[ettersi]

Not sure what's up with the failed test. It only fails on Windows and I believe doesn't involve code I touched, so I'd like to think this is unrelated..

[juliohm]

Amazing work as usual ❤️ I will try to review it by the end of the day.

Not sure what's up with the failed test. It only fails on Windows and I believe doesn't involve code I touched, so I'd like to think this is unrelated..

All tests are passing after a second try 👍🏽

[ettersi]

Tests pass locally 🤷

[juliohm]

Any idea of what might be causing the test failure on Julia LTS? I opened #270, which when solved will help split the tests further.

[juliohm]

I fixed that issue. We can now rebase this PR to refine the conditions where the changes fail.

[juliohm]

So the tests are only failing with Float32 on Julia LTS according to CI.

[ettersi]

Found the issue. In Julia 1.11.5:

julia> promote(π, 1f0°)[1]
3.1415927f0

In Julia v1.10.9 (LTS):

julia> promote(π, 1f0°)[1]
3.141592653589793

[juliohm]

Perhaps we should file a bug in Julia LTS? Or is this coming from an old Unitful.jl version on LTS?

[juliohm]

The tests sometimes fail due to some randomness in isapproxtest3D, which is unrelated to the changes here. I restarted the tests and all of them are green again.

I will review the changes again tonight (BRT time zone) to see if there is something else we can improve before merging.

[juliohm]

I've started looking more carefully into the changes and realized that we already define isapprox as part of the public API. I believe we just need to specialize the methods for different CRS types instead of defining allapprox (renamed to isapproxcoords) in testutils.jl. That will give us full access to the internals of each CRS type in the implementations.

Can we experiment with this approach? It should also improve the situation in downstream tests in CoordGridTransforms.jl, which also define allapprox (code duplication).

[ettersi]

Perhaps we should file a bug in Julia LTS? Or is this coming from an old Unitful.jl version on LTS?

~~Will do as soon as I find the time. Also don't mind if someone else beats me to it~~ https://github.com/PainterQubits/Unitful.jl/issues/778

The tests sometimes fail due to some randomness in isapproxtest3D, which is unrelated to the changes here. I restarted the tests and all of them are green again.

You can avoid such non-reproducible failures by putting Random.seed!(42) at the beginning of your runtests.jl. Also, I believe the failures happen because some of the tests apply tolerances intended to be used relative to the radius of the earth (CoordRefSystems.atol()) to quantities that have O(1) scale (e.g. rand(Cartesian3D) which produces a point in the unit cube [0,1]^3). Haven't looked into this in detail, so might be wrong.

I've started looking more carefully into the changes and realized that we already define isapprox as part of the public API. I believe we just need to specialize the methods for different CRS types instead of defining allapprox (renamed to isapproxcoords) in testutils.jl. That will give us full access to the internals of each CRS type in the implementations.

Personally, I would say the (conceptual) definition of isapprox() should be

isapprox(c1::CRS, c2::CRS) = isapprox(
    svec(convert(Cartesian3D, c1)), 
    svec(convert(Cartesian3D, c2)),
)

and it seems you agree since you've added the docstring

  isapprox(coords₁, coords₂; kwargs...)

  Converts coords₁ and coords₂ to Cartesian coordinates and compare the coordinate values with the isapprox method of vectors. The conversion to Cartesian coordinates takes care of possibly different Datum.

The isapproxcoords() methods defined in this PR do not satisfy this definition (e.g. isapproxcoords(LatLon(90,0), LatLon(90, 180)) -> false). Therefore, I don't think they should be moved to isapprox().

[juliohm]

You can avoid such non-reproducible failures by putting Random.seed!(42) at the beginning of your runtests.jl

It is in the plans. Just trying to sort out this isapprox definition first.

The isapproxcoords() methods defined in this PR do not satisfy this definition (e.g. isapproxcoords(LatLon(90,0), LatLon(90, 180)) -> false). Therefore, I don't think they should be moved to isapprox().

Apparently we have a couple of isapprox definitions that we need to review carefully:

• CoordRefSystems.isapprox for CRS
• Meshes.isapprox for Point
• isapproxcoords in testutils.jl

I believe we just need two definitions. The one in Meshes.jl seems to be correct. It compares two points in the physical world via a distance in Cartesian coordinates using the vector from the origin to the point.

The one in CoordRefSystems.isapprox may be wrong. We just want to compare the coordinate values as a tuple of numbers. In the case of Cartesian CRS, the definition coincides with the one for SVector. The other cases are the ones that are problematic.

The isapproxcoords definition in testutils.jl seems to be doing what we need in CoordRefSystems.isapprox, so we could probably rewrite CoordRefSystems.isapprox to have this new meaning. It is a breaking change.

Please let me know if you see some flaw with this design. I will look into the definitions closely again.

[ettersi]

Please let me know if you see some flaw with this design. I will look into the definitions closely again.

Makes sense. I don't think I have any further inputs on this.

[juliohm]

@ettersi can we port some of these changes to the latest main branch?

[ettersi]

Sure! Please feel free to scavenge from here whatever you find useful!

[juliohm]

Thank you @ettersi ! Closing this in favor of #294