Coordinate Transformations

[mileslucas]

This PR addresses extra comments from #48

Right now I am doing timing/accuracy tests for coords2cart and cart2coords to assess sqrt vs hypot vs CoordinateTransformations.jl

[codecov]

Codecov Report

Patch and project coverage have no change.

Comparison is base (f69923e) 90.47% compared to head (56910f2) 90.47%.

Additional details and impacted files

@@           Coverage Diff           @@
##           master      #51   +/-   ##
=======================================
  Coverage   90.47%   90.47%           
=======================================
  Files           3        3           
  Lines         126      126           
=======================================
  Hits          114      114           
  Misses         12       12           

Impacted Files	Coverage Δ
src/SkyCoords.jl	98.68% <ø> (ø)

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

☔ View full report at Codecov.
📢 Do you have feedback about the report comment? Let us know in this issue.

[mileslucas]

Benchmarks for coordinate transformations

Setup:

using BenchmarkTools
lon = randn() * pi + pi
lat = randn() * pi
function maxerr(func, invfunc)
	test_lon = range(0, pi, length=100)
	test_lat = range(-pi, pi, length=100)
	test_big_lon = range(big(0), big(pi), length=100)
	test_big_lat = range(-big(pi), big(pi), length=100)
	
	res = func.(test_lon, test_lat)
	res_big = func.(test_big_lon, test_big_lat)
	delt = mapreduce((a,b) -> sum(a .- b), max, res, res_big)

	invres = invfunc.(res)
	invres_big = invfunc.(res_big)
	invdelt = mapreduce((a,b) -> sum(a .- b), max, invres, invres_big)
	return delt, invdelt
end

sqrt

This is the current implementation

function coords2cart(lon, lat)
    sinlon, coslon = sincos(lon)
    sinlat, coslat = sincos(lat)
    SVector{3}(coslat * coslon, coslat * sinlon, sinlat)
end
function cart2coords(v)
    lon = atan(v[begin + 1], v[begin])
    xy_norm = sqrt(v[begin]^2 + v[begin + 1]^2)
    lat = atan(v[begin + 2], xy_norm)
    return lon, lat
end

Results

julia> @benchmark cart2coords(coords2cart($lon, $lat))
BenchmarkTools.Trial: 10000 samples with 991 evaluations.
 Range (min … max):  40.406 ns … 65.464 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     40.573 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   40.780 ns ±  1.135 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▅█▇▅▂                             ▁▁                        ▂
  █████▇▅▇▅▇▆▅▇▆▅████▇▅▆▄▆▆▁▃▅▅▅▅▆▅▅██▇▇▄▄▅▅▅▅▃▄▃▃▃▅▄▄▁▃▄▁▁▁▃ █
  40.4 ns      Histogram: log(frequency) by time      45.4 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> maxerr(coords2cart, cart2coords)
(2.471323154795485647167061374096440014070139517448629832542271562123740152545622e-16, 3.551145583454658824002860858772591006876497547505050342092755672005132716512651e-16)

hypot

function cart2coords_hypot(v)
    lon = atan(v[begin + 1], v[begin])
    xy_norm = hypot(v[begin], v[begin + 1])
    lat = atan(v[begin + 2], xy_norm)
    return lon, lat
end

Results

julia> @benchmark cart2coords_hypot(coords2cart($lon, $lat))
BenchmarkTools.Trial: 10000 samples with 986 evaluations.
 Range (min … max):  52.146 ns … 78.939 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     52.274 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   52.425 ns ±  0.892 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▅▇█▄▄▃▂                                          ▁ ▁        ▂
  ████████▅▅▄▃▄▄▄▃▁▁▁▅▃▄▄▄▃▄▃▅▅▆▄▃▄▄▄▄▃▃▁▁▃▃▃▁▄▁▃▅████▆▅▃▄▄▄▅ █
  52.1 ns      Histogram: log(frequency) by time      55.7 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> maxerr(coords2cart, cart2coords_hypot)
(2.471323154795485647167061374096440014070139517448629832542271562123740152545622e-16, 3.551145583454658824002860858772591006876497547505050342092755672005132716512651e-16)

CoordinateTransformations.jl

using CoordinateTransformations
function coords2cart_ct(lon, lat)
	sph = Spherical(1, lon, lat)
    return CartesianFromSpherical()(sph)
end
function cart2coords_ct(v)
	sph = SphericalFromCartesian()(v)
    return sph.θ, sph.ϕ
end

Results

julia> @benchmark cart2coords_ct(coords2cart_ct($lon, $lat))
BenchmarkTools.Trial: 10000 samples with 982 evaluations.
 Range (min … max):  57.663 ns … 81.635 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     57.832 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   58.198 ns ±  1.126 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▂██▄▄▂                ▃▃▂▂▁       ▁▂▁                       ▂
  ███████▇▆▅▅▅▆▆▆▄▃▁▅▄▅▆██████▇▆▆▆▅▇███▇▆▇▅▅▅▅▆▆▆▆▆▅▆▆▆▆▆▄▃▄▆ █
  57.7 ns      Histogram: log(frequency) by time      62.6 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> maxerr(coords2cart_ct, cart2coords_ct)
(2.471323154795485647167061374096440014070139517448629832542271562123740152545622e-16, 3.551145583454658824002860858772591006876497547505050342092755672005132716512651e-16)

[mileslucas]

cc @giordano any comments on these results? If you have suggestions for testing differently I'm happy to try them, too.

[giordano]

The problem of sqrt+squares is that it under/overflows:

julia> sqrt(1e300 ^ 2 + 1e300 ^ 2)
Inf

julia> hypot(1e300, 1e300)
1.4142135623730952e300

julia> sqrt(1e-300 ^ 2 + 1e-300 ^ 2)
0.0

julia> hypot(1e-300, 1e-300)
1.414213562373095e-300

However, as you can see, hypot is slower because it does some more stuff to avoid that. So there's always a trade-off between speed and accuracy here, someone would need to make a call about what to favour 🙂

[mileslucas]

I think a 20 nanosecond hit is worth the accuracy, personally.