Add BrayCurtis (#76)

add BrayCurtis dissimilarity (add tests, benchmarks, update README)

Author: kescobo

README.md

@@ -36,6 +36,7 @@ This package also provides optimized functions to compute column-wise and pairwi
 * Mean squared deviation
 * Root mean squared deviation
 * Normalized root mean squared deviation
+* Bray-Curtis dissimilarity
 
 For ``Euclidean distance``, ``Squared Euclidean distance``, ``Cityblock distance``, ``Minkowski distance``, and ``Hamming distance``, a weighted version is also provided.
 
@@ -139,6 +140,7 @@ Each distance corresponds to a distance type. The type name and the correspondin
 |  Hamming             |  `hamming(k, l)`           | `sum(k .!= l)` |
 |  RogersTanimoto      |  `rogerstanimoto(a, b)`    | `2(sum(a&!b) + sum(!a&b)) / (2(sum(a&!b) + sum(!a&b)) + sum(a&b) + sum(!a&!b))` |
 |  Jaccard             |  `jaccard(x, y)`           | `1 - sum(min(x, y)) / sum(max(x, y))` |
+|  BrayCurtis          |  `braycurtis(x, y)`        | `sum(abs(x - y)) / sum(abs(x + y))`  |
 |  CosineDist          |  `cosine_dist(x, y)`       | `1 - dot(x, y) / (norm(x) * norm(y))` |
 |  CorrDist            |  `corr_dist(x, y)`         | `cosine_dist(x - mean(x), y - mean(y))` |
 |  ChiSqDist           |  `chisq_dist(x, y)`        | `sum((x - y).^2 / (x + y))` |

benchmark/benchmarks.jl

@@ -20,6 +20,8 @@ function create_distances(w, Q)
         HellingerDist(),
 
         Haversine(6371.),
+  
+        BrayCurtis(),
 
         WeightedSqEuclidean(w),
         WeightedEuclidean(w),

benchmark/print_table.jl

@@ -32,6 +32,7 @@ order = [
     :SqMahalanobis,
     :Mahalanobis,
     :Haversine,
+    :BrayCurtis,
 ]
 
 BenchmarkTools.DEFAULT_PARAMETERS.seconds = 2.0 # Long enough

src/Distances.jl

@@ -23,6 +23,7 @@ export
     Chebyshev,
     Minkowski,
     Jaccard,
+    BrayCurtis,
     RogersTanimoto,
 
     Hamming,
@@ -57,6 +58,7 @@ export
     sqeuclidean,
     cityblock,
     jaccard,
+    braycurtis,
     rogerstanimoto,
     chebyshev,
     minkowski,

src/metrics.jl

@@ -25,6 +25,7 @@ struct Hamming <: Metric end
 
 struct CosineDist <: SemiMetric end
 struct CorrDist <: SemiMetric end
+struct BrayCurtis <: SemiMetric end
 
 struct ChiSqDist <: SemiMetric end
 struct KLDivergence <: PreMetric end
@@ -98,7 +99,7 @@ struct RMSDeviation <: Metric end
 struct NormRMSDeviation <: Metric end
 
 
-const UnionMetrics = Union{Euclidean,SqEuclidean,Chebyshev,Cityblock,Minkowski,Hamming,Jaccard,RogersTanimoto,CosineDist,CorrDist,ChiSqDist,KLDivergence,RenyiDivergence,JSDivergence,SpanNormDist,GenKLDivergence}
+const UnionMetrics = Union{Euclidean,SqEuclidean,Chebyshev,Cityblock,Minkowski,Hamming,Jaccard,RogersTanimoto,CosineDist,CorrDist,ChiSqDist,KLDivergence,RenyiDivergence,BrayCurtis,JSDivergence,SpanNormDist,GenKLDivergence}
 
 """
     Euclidean([thresh])
@@ -356,6 +357,27 @@ end
 end
 jaccard(a::AbstractArray, b::AbstractArray) = evaluate(Jaccard(), a, b)
 
+# BrayCurtis
+
+@inline eval_start(::BrayCurtis, a::AbstractArray{Bool}, b::AbstractArray{Bool}) = 0, 0
+@inline eval_start(::BrayCurtis, a::AbstractArray{T}, b::AbstractArray{T}) where {T} = zero(T), zero(T)
+@inline function eval_op(::BrayCurtis, s1, s2)
+    abs_m = abs(s1 - s2)
+    abs_p = abs(s1 + s2)
+    abs_m, abs_p
+end
+@inline function eval_reduce(::BrayCurtis, s1, s2)
+    @inbounds a = s1[1] + s2[1]
+    @inbounds b = s1[2] + s2[2]
+    a, b
+end
+@inline function eval_end(::BrayCurtis, a)
+    @inbounds v = a[1] / a[2]
+    return v
+end
+braycurtis(a::AbstractArray, b::AbstractArray) = evaluate(BrayCurtis(), a, b)
+
+
 # Tanimoto
 
 @inline eval_start(::RogersTanimoto, a::AbstractArray, b::AbstractArray) = 0, 0, 0, 0

test/test_dists.jl

@@ -125,6 +125,8 @@ end
     bt = [true, false, true]
     bf = [false, true, true]
     @test rogerstanimoto(bt, bf) == 4.0 / 5.0
+    @test braycurtis(bt, bf) == 0.5
+
     for T in (Float64, F64)
 
         for (_x, _y) in (([4.0, 5.0, 6.0, 7.0], [3.0, 9.0, 8.0, 1.0]),
@@ -135,6 +137,7 @@ end
             @test jaccard(x, y) == 13.0 / 28
             @test cityblock(x, y) == 13.0
             @test chebyshev(x, y) == 6.0
+            @test braycurtis(x, y) == 1. - (30./43.)
             @test minkowski(x, y, 2) == sqrt(57.0)
             @test_throws DimensionMismatch cosine_dist(1.0:2, 1.0:3)
             @test cosine_dist(x, y) ≈ (1.0 - 112. / sqrt(19530.0))
@@ -177,6 +180,11 @@ end
         @inferred evaluate(Jaccard(), [1, 2, 3], [1, 2, 3])
         @inferred evaluate(Jaccard(), [true, false, true], [false, true, true])
 
+        # Test Bray-Curtis. Should be 1 if no elements are shared, 0 if all are the same
+        @test braycurtis([1,0,3],[0,1,0]) == 1.0
+        @test braycurtis(rand(10), zeros(10)) == 1.0
+        @test braycurtis([1,0],[1,0]) == 0.0
+
         # Test KL, Renyi and JS divergences
         r = rand(T, 12)
         p = copy(r)
@@ -232,12 +240,16 @@ end #testset
         @test isa(cityblock(a, b), T)
         @test chebyshev(a, b) == 0.0
         @test isa(chebyshev(a, b), T)
+        @test braycurtis(a, b) == 0.0
+        @test isa(braycurtis(a, b), T)
         @test minkowski(a, b, 2) == 0.0
         @test isa(minkowski(a, b, 2), T)
         @test hamming(a, b) == 0.0
         @test isa(hamming(a, b), Int)
         @test renyi_divergence(a, b, 1.0) == 0.0
         @test isa(renyi_divergence(a, b, 2.0), T)
+        @test braycurtis(a, b) == 0.0
+        @test isa(braycurtis(a, b), T)
 
         w = T[]
         @test isa(whamming(a, b, w), T)
@@ -384,6 +396,7 @@ end
 
         test_colwise(BhattacharyyaDist(), X, Y, T)
         test_colwise(HellingerDist(), X, Y, T)
+        test_colwise(BrayCurtis(), X, Y, T)
 
         w = rand(T, m)
 
@@ -456,6 +469,7 @@ end
 
         test_pairwise(BhattacharyyaDist(), X, Y, T)
         test_pairwise(HellingerDist(), X, Y, T)
+        test_pairwise(BrayCurtis(), X, Y, T)
 
         w = rand(m)