Allow generators and iterators

src/Distances.jl

@@ -97,6 +97,15 @@ export
     rmsd,
     nrmsd
 
+if VERSION < v"1.2-"
+    import Base: has_offset_axes
+    require_one_based_indexing(A...) = 
+        !has_offset_axes(A...) || 
+            throw(ArgumentError("offset arrays are not supported but got an array with index other than 1"))
+else
+    import Base: require_one_based_indexing
+end
+
 include("common.jl")
 include("generic.jl")
 include("metrics.jl")

src/bhattacharyya.jl

@@ -9,39 +9,54 @@ struct HellingerDist <: Metric end
 
 # Bhattacharyya coefficient
 
-function bhattacharyya_coeff(a::AbstractVector{T}, b::AbstractVector{T}) where {T <: Number}
-    if length(a) != length(b)
-        throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b))."))
+function bhattacharyya_coeff(a, b)
+    n = length(a)
+    if n != length(b)
+        throw(DimensionMismatch("first argument has length $n which does not match the length of the second, $(length(b))."))
     end
+    sqab, asum, bsum = _bhattacharyya_coeff(a, b)
+    # We must normalize since we cannot assume that the vectors are normalized to probability vectors.
+    return sqab / sqrt(asum * bsum)
+end
 
-    n = length(a)
+@inline function _bhattacharyya_coeff(a, b)
+    Ta = _eltype(a)
+    Tb = _eltype(b)
+    T = typeof(sqrt(zero(promote_type(Ta, Tb))))
     sqab = zero(T)
-    # We must normalize since we cannot assume that the vectors are normalized to probability vectors.
-    asum = zero(T)
-    bsum = zero(T)
+    asum = zero(Ta)
+    bsum = zero(Tb)
+
+    for (ai, bi) in zip(a, b)
+        sqab += sqrt(ai * bi)
+        asum += ai
+        bsum += bi
+    end
+    return sqab, asum, bsum
+end
+@inline function _bhattacharyya_coeff(a::AbstractVector{Ta}, b::AbstractVector{Tb}) where {Ta<:Number,Tb<:Number}
+    T = typeof(sqrt(oneunit(Ta)*oneunit(Tb)))
+    sqab = zero(T)
+    asum = zero(Ta)
+    bsum = zero(Tb)
 
-    @simd for i = 1:n
+    @simd for i = 1:length(a)
         @inbounds ai = a[i]
         @inbounds bi = b[i]
         sqab += sqrt(ai * bi)
         asum += ai
         bsum += bi
     end
-
-    sqab / sqrt(asum * bsum)
+    return sqab, asum, bsum
 end
 
-bhattacharyya_coeff(a::T, b::T) where {T <: Number} = throw("Bhattacharyya coefficient cannot be calculated for scalars")
-
 # Faster pair- and column-wise versions TBD...
 
 
 # Bhattacharyya distance
-(::BhattacharyyaDist)(a::AbstractVector{T}, b::AbstractVector{T}) where {T <: Number} = -log(bhattacharyya_coeff(a, b))
-(::BhattacharyyaDist)(a::T, b::T) where {T <: Number} = throw("Bhattacharyya distance cannot be calculated for scalars")
+(::BhattacharyyaDist)(a, b) = -log(bhattacharyya_coeff(a, b))
 bhattacharyya(a, b) = BhattacharyyaDist()(a, b)
 
 # Hellinger distance
-(::HellingerDist)(a::AbstractVector{T}, b::AbstractVector{T}) where {T <: Number} = sqrt(1 - bhattacharyya_coeff(a, b))
-(::HellingerDist)(a::T, b::T) where {T <: Number} = throw("Hellinger distance cannot be calculated for scalars")
+(::HellingerDist)(a, b) = sqrt(1 - bhattacharyya_coeff(a, b))
 hellinger(a, b) = HellingerDist()(a, b)

src/bregman.jl

@@ -22,26 +22,26 @@ end
 Bregman(F, ∇) =  Bregman(F, ∇, LinearAlgebra.dot)
 
 # Evaluation fuction 
-function (dist::Bregman)(p::AbstractVector, q::AbstractVector)
+function (dist::Bregman)(p, q)
     # Create cache vals.
-    FP_val = dist.F(p);
-    FQ_val = dist.F(q); 
-    DQ_val = dist.∇(q);
-    p_size = size(p);
+    FP_val = dist.F(p)
+    FQ_val = dist.F(q) 
+    DQ_val = dist.∇(q)
+    p_size = length(p)
     # Check F codomain. 
     if !(isa(FP_val, Real) && isa(FQ_val, Real))
         throw(ArgumentError("F Codomain Error: F doesn't map the vectors to real numbers"))
     end 
     # Check vector size. 
-    if !(p_size == size(q))
+    if p_size != length(q)
         throw(DimensionMismatch("The vector p ($(size(p))) and q ($(size(q))) are different sizes."))
     end
     # Check gradient size. 
-    if !(size(DQ_val) == p_size)
+    if length(DQ_val) != p_size
         throw(DimensionMismatch("The gradient result is not the same size as p and q"))
     end 
     # Return the Bregman divergence. 
-    return FP_val - FQ_val - dist.inner(DQ_val, p-q);
+    return FP_val - FQ_val - dist.inner(DQ_val, p .- q)
 end 
 
 # Convenience function. 

src/common.jl

@@ -12,6 +12,12 @@ function get_common_ncols(a::AbstractMatrix, b::AbstractMatrix)
     return na
 end
 
+function get_common_length(a, b)
+    n = length(a)
+    length(b) == n || throw(DimensionMismatch("The lengths of a and b must match."))
+    return n
+end
+
 function get_pairwise_dims(r::AbstractMatrix, a::AbstractMatrix, b::AbstractMatrix)
     ma, na = size(a)
     mb, nb = size(b)

src/generic.jl

@@ -26,40 +26,64 @@ evaluate(dist::PreMetric, a, b) = dist(a, b)
 # Generic functions
 
 """
-    result_type(dist::PreMetric, Ta::Type, Tb::Type) -> T
-    result_type(dist::PreMetric, a::AbstractArray, b::AbstractArray) -> T
+    result_type(dist, Ta::Type, Tb::Type) -> T
+    result_type(dist, a, b) -> T
 
 Infer the result type of metric `dist` with input type `Ta` and `Tb`, or input
 data `a` and `b`.
 """
-result_type(::PreMetric, ::Type, ::Type) = Float64 # fallback
-result_type(dist::PreMetric, a::AbstractArray, b::AbstractArray) = result_type(dist, eltype(a), eltype(b))
+# result_type(::PreMetric, ::Type, ::Type) = Float64 # fallback in Distances
+result_type(f, a::Type, b::Type) = typeof(f(oneunit(a), oneunit(b))) # don't require `PreMetric` subtyping
 
+result_type(dist, a, b) = result_type(dist, _eltype(a), _eltype(b))
+# description of approach:
+# (a) for generic iterators, rely on Base.IteratorEltype(a)
+_eltype(a) = __eltype(Base.IteratorEltype(a), a)
+__eltype(::Base.HasEltype, a) = _eltype(eltype(a))
+__eltype(::Base.EltypeUnknown, a) = _eltype(typeof(first(a)))
+# (b) for arrays (of arrays), go down the eltype route, until you hit a non-array eltype and return that
+_eltype(a::AbstractArray) = _eltype(eltype(a))
+_eltype(::Type{T}) where {T<:AbstractArray} = _eltype(eltype(T))
+_eltype(T::Type) = T
 
 # Generic column-wise evaluation
 
+function colwise!(r::AbstractArray, metric::PreMetric, a, b)
+    require_one_based_indexing(r)
+    n = length(a)
+    length(b) == n || throw(DimensionMismatch("iterators have different lengths"))
+    length(r) == n || throw(DimensionMismatch("Incorrect size of r."))
+    @inbounds for (j, ab) in enumerate(zip(a, b))
+        r[j] = metric(ab...)
+    end
+    r
+end
+
 function colwise!(r::AbstractArray, metric::PreMetric, a::AbstractVector, b::AbstractMatrix)
+    require_one_based_indexing(r)
     n = size(b, 2)
     length(r) == n || throw(DimensionMismatch("Incorrect size of r."))
-    @inbounds for j = 1:n
-        r[j] = metric(a, view(b, :, j))
+    @inbounds for (rj, bj) in enumerate(axes(b, 2))
+        r[rj] = metric(a, view(b, :, bj))
     end
     r
 end
 
 function colwise!(r::AbstractArray, metric::PreMetric, a::AbstractMatrix, b::AbstractVector)
+    require_one_based_indexing(r)
     n = size(a, 2)
     length(r) == n || throw(DimensionMismatch("Incorrect size of r."))
-    @inbounds for j = 1:n
-        r[j] = metric(view(a, :, j), b)
+    @inbounds for (rj, aj) in enumerate(axes(a, 2))
+        r[rj] = metric(view(a, :, aj), b)
     end
     r
 end
 
 function colwise!(r::AbstractArray, metric::PreMetric, a::AbstractMatrix, b::AbstractMatrix)
+    require_one_based_indexing(r, a, b)
     n = get_common_ncols(a, b)
     length(r) == n || throw(DimensionMismatch("Incorrect size of r."))
-    @inbounds for j = 1:n
+    @inbounds for j in 1:n
         r[j] = metric(view(a, :, j), view(b, :, j))
     end
     r
@@ -69,6 +93,12 @@ function colwise!(r::AbstractArray, metric::SemiMetric, a::AbstractMatrix, b::Ab
     colwise!(r, metric, b, a)
 end
 
+function colwise(metric::PreMetric, a, b)
+    n = get_common_length(a, b)
+    r = Vector{result_type(metric, a, b)}(undef, n)
+    colwise!(r, metric, a, b)
+end
+
 function colwise(metric::PreMetric, a::AbstractMatrix, b::AbstractMatrix)
     n = get_common_ncols(a, b)
     r = Vector{result_type(metric, a, b)}(undef, n)
@@ -90,8 +120,20 @@ end
 
 # Generic pairwise evaluation
 
+function _pairwise!(r::AbstractMatrix, metric::PreMetric, a, b=a)
+    require_one_based_indexing(r)
+    na = length(a)
+    nb = length(b)
+    size(r) == (na, nb) || throw(DimensionMismatch("Incorrect size of r."))
+    @inbounds for (j, bj) in enumerate(b), (i, ai) in enumerate(a)
+        r[i, j] = metric(ai, bj)
+    end
+    r
+end
+
 function _pairwise!(r::AbstractMatrix, metric::PreMetric,
                     a::AbstractMatrix, b::AbstractMatrix=a)
+    require_one_based_indexing(r, a, b)
     na = size(a, 2)
     nb = size(b, 2)
     size(r) == (na, nb) || throw(DimensionMismatch("Incorrect size of r."))
@@ -104,18 +146,34 @@ function _pairwise!(r::AbstractMatrix, metric::PreMetric,
     r
 end
 
+function _pairwise!(r::AbstractMatrix, metric::SemiMetric, a)
+    require_one_based_indexing(r)
+    n = length(a)
+    size(r) == (n, n) || throw(DimensionMismatch("Incorrect size of r."))
+    itr = Iterators.product(enumerate(a), enumerate(a))
+    @inbounds for ((i, ai), (j, aj)) in itr
+        r[i, j] = i > j ? 
+            metric(ai, aj) :
+            i == j ?
+                0 :
+                r[j,i]
+    end
+    r
+end
+
 function _pairwise!(r::AbstractMatrix, metric::SemiMetric, a::AbstractMatrix)
+    require_one_based_indexing(r)
     n = size(a, 2)
     size(r) == (n, n) || throw(DimensionMismatch("Incorrect size of r."))
     @inbounds for j = 1:n
+        for i = 1:(j - 1)
+            r[i, j] = r[j, i]   # leveraging the symmetry of SemiMetric
+        end
+        r[j, j] = 0
         aj = view(a, :, j)
         for i = (j + 1):n
             r[i, j] = metric(view(a, :, i), aj)
         end
-        r[j, j] = 0
-        for i = 1:(j - 1)
-            r[i, j] = r[j, i]   # leveraging the symmetry of SemiMetric
-        end
     end
     r
 end
@@ -204,6 +262,13 @@ function pairwise(metric::PreMetric, a::AbstractMatrix, b::AbstractMatrix;
     pairwise!(r, metric, a, b, dims=dims)
 end
 
+function pairwise(metric::PreMetric, a, b)
+    m = length(a)
+    n = length(b)
+    r = Matrix{result_type(metric, a, b)}(undef, m, n)
+    _pairwise!(r, metric, a, b)
+end
+
 function pairwise(metric::PreMetric, a::AbstractMatrix;
                   dims::Union{Nothing,Integer}=nothing)
     dims = deprecated_dims(dims)
@@ -212,3 +277,9 @@ function pairwise(metric::PreMetric, a::AbstractMatrix;
     r = Matrix{result_type(metric, a, a)}(undef, n, n)
     pairwise!(r, metric, a, dims=dims)
 end
+
+function pairwise(metric::PreMetric, a)
+    n = length(a)
+    r = Matrix{result_type(metric, a, a)}(undef, n, n)
+    _pairwise!(r, metric, a)
+end

src/haversine.jl

@@ -10,20 +10,17 @@ struct Haversine{T<:Real} <: Metric
     radius::T
 end
 
-const VecOrLengthTwoTuple{T} = Union{AbstractVector{T}, NTuple{2, T}}
-
-function (dist::Haversine)(x::VecOrLengthTwoTuple, y::VecOrLengthTwoTuple)
+function (dist::Haversine)(x, y)
     length(x) == length(y) == 2 || haversine_error()
 
-    @inbounds begin
-        # longitudes
-        Δλ = deg2rad(y[1] - x[1])
-
-        # latitudes
-        φ₁ = deg2rad(x[2])
-        φ₂ = deg2rad(y[2])
-    end
+    @inbounds x1, x2 = x
+    @inbounds y1, y2 = y
+    # longitudes
+    Δλ = deg2rad(y1 - x1)
 
+    # latitudes
+    φ₁ = deg2rad(x2)
+    φ₂ = deg2rad(y2)
     Δφ = φ₂ - φ₁
 
     # haversine formula
@@ -33,6 +30,6 @@ function (dist::Haversine)(x::VecOrLengthTwoTuple, y::VecOrLengthTwoTuple)
     2 * dist.radius * asin( min(√a, one(a)) ) # take care of floating point errors
 end
 
-haversine(x::VecOrLengthTwoTuple, y::VecOrLengthTwoTuple, radius::Real) = Haversine(radius)(x, y)
+haversine(x, y, radius::Real) = Haversine(radius)(x, y)
 
 @noinline haversine_error() = throw(ArgumentError("expected both inputs to have length 2 in Haversine distance"))