Small improvements to Haversine (#226)

Co-authored-by: Milan Bouchet-Valat <[email protected]>

Author: torfjelde

src/haversine.jl

@@ -11,7 +11,8 @@ The computed distance has the unit of the radius.
 struct Haversine{T<:Number} <: Metric
     radius::T
 end
-Haversine() = Haversine(Float32(6_371_000))
+Haversine{T}() where {T<:Number} = Haversine(T(6_371_000))
+Haversine() = Haversine{Int}()
 
 function (dist::Haversine)(x, y)
     length(x) == length(y) == 2 || haversine_error(dist)
@@ -26,13 +27,16 @@ function (dist::Haversine)(x, y)
     a = sind(Δφ/2)^2 + cosd(φ₁)*cosd(φ₂)*sind(Δλ/2)^2
 
     # distance on the sphere
-    2 * dist.radius * asin( min(√a, one(a)) ) # take care of floating point errors
+    2 * (dist.radius * asin( min(√a, one(a)) )) # take care of floating point errors
 end
 
-haversine(x, y, radius::Number=Float32(6_371_000)) = Haversine(radius)(x, y)
+haversine(x, y, radius::Number=6_371_000) = Haversine(radius)(x, y)
 
 @noinline haversine_error(dist) = throw(ArgumentError("expected both inputs to have length 2 in $dist distance"))
 
+result_type(::Haversine{T1}, ::Type{T2}, ::Type{T3}) where {T1<:Number,T2<:Number,T3<:Number} =
+    float(promote_type(T1, T2, T3))
+
 """
     SphericalAngle()
 
@@ -61,4 +65,4 @@ end
 
 const spherical_angle = SphericalAngle()
 
-result_type(::Union{Haversine, SphericalAngle}, ::Type, ::Type) = Float64
+result_type(::SphericalAngle, ::Type{T1}, ::Type{T2}) where {T1<:Number,T2<:Number} = float(promote_type(T1, T2))

test/test_dists.jl

@@ -60,6 +60,36 @@ function test_metricity(dist, x, y, z)
     end
 end
 
+function test_batch(dist, _x, _y, args=(); Tins=(Float32, Float64), Touts=Tins)
+    _xs = hcat(_x, _x)
+    _ys = hcat(_y, _y)
+
+    for (Tin, Tout) in zip(Tins, Touts)
+        x, y = Tin.(_x), Tin.(_y)
+        xs, ys = Tin.(_xs), Tin.(_ys)
+
+        # Instantiate using `args` if not already instantiated.
+        d = if dist isa PreMetric
+            dist
+        else
+            dist(map(a -> Tin.(a), args)...)
+        end
+
+        # Ensure that `result_type` preserves the type.
+        @test result_type(d, xs, ys) === Tout
+
+        # Ensure that `pairwise` preserves the type,
+        # e.g. if someone decides to implement their own.
+        colpairs = pairwise(d, xs, ys; dims=2)
+        @test eltype(colpairs) === Tout
+        @test all(colpairs .== d(x, y))
+
+        cols = colwise(d, xs, ys)
+        @test eltype(cols) === Tout
+        @test all(cols .== d(x, y))
+    end
+end
+
 @testset "PreMetric, SemiMetric, Metric on $T" for T in (Float64, F64)
     Random.seed!(123)
     n = 100
@@ -424,6 +454,17 @@ end
 end #testset
 
 @testset "haversine" begin
+    let x = [1.,-15.625], y = [-179.,15.625]
+        # Type used in default constructor should work for both `Float32` and `Float64`.
+        test_batch(Haversine(), x, y)
+
+        # `Float32` should be promoted to `Float64`.
+        test_batch(Haversine{Float32}(), x, y)
+
+        # `Haversine(T(radius))` should result in `T` when inputs also has eltype `T`.
+        test_batch(Haversine, x, y, (6_371_000, ))
+    end
+
     for T in (Float64, F64)
         @test haversine([-180.,0.], [180.,0.], 1.) ≈ 0 atol=1e-10
         @test haversine([0.,-90.],  [0.,90.],  1.) ≈ π atol=1e-10