Add distances, divergences, and deviations from StatsBase (#58)

Author: ararslan

README.md

@@ -23,12 +23,17 @@ This package also provides optimized functions to compute column-wise and pairwi
 * Correlation distance
 * Chi-square distance
 * Kullback-Leibler divergence
+* Generalized Kullback-Leibler divergence
 * Rényi divergence
 * Jensen-Shannon divergence
 * Mahalanobis distance
 * Squared Mahalanobis distance
 * Bhattacharyya distance
 * Hellinger distance
+* Mean absolute deviation
+* Mean squared deviation
+* Root mean squared deviation
+* Normalized root mean squared deviation
 
 For ``Euclidean distance``, ``Squared Euclidean distance``, ``Cityblock distance``, ``Minkowski distance``, and ``Hamming distance``, a weighted version is also provided.
 
@@ -138,13 +143,18 @@ Each distance corresponds to a distance type. The type name and the correspondin
 |  CorrDist            |  `corr_dist(x, y)`         | `cosine_dist(x - mean(x), y - mean(y))` |
 |  ChiSqDist           |  `chisq_dist(x, y)`        | `sum((x - y).^2 / (x + y))` |
 |  KLDivergence        |  `kl_divergence(p, q)`     | `sum(p .* log(p ./ q))` |
+|  GenKLDivergence     |  `gkl_divergence(x, y)`    | `sum(p .* log(p ./ q) - p + q)` |
 |  RenyiDivergence     | `renyi_divergence(p, q, k)`| `log(sum( p .* (p ./ q) .^ (k - 1))) / (k - 1)` |
 |  JSDivergence        |  `js_divergence(p, q)`     | `KL(p, m) / 2 + KL(p, m) / 2 with m = (p + q) / 2` |
 |  SpanNormDist        |  `spannorm_dist(x, y)`     | `max(x - y) - min(x - y )` |
 |  BhattacharyyaDist   |  `bhattacharyya(x, y)`     | `-log(sum(sqrt(x .* y) / sqrt(sum(x) * sum(y)))` |
 |  HellingerDist       |  `hellinger(x, y) `        | `sqrt(1 - sum(sqrt(x .* y) / sqrt(sum(x) * sum(y))))` |
 |  Mahalanobis         |  `mahalanobis(x, y, Q)`    | `sqrt((x - y)' * Q * (x - y))` |
 |  SqMahalanobis       |  `sqmahalanobis(x, y, Q)`  | ` (x - y)' * Q * (x - y)`  |
+|  MeanAbsDeviation    |  `meanad(x, y)`            | `mean(abs.(x - y))` |
+|  MeanSqDeviation     |  `msd(x, y)`               | `mean(abs2.(x - y))` |
+|  RMSDeviation        |  `rmsd(x, y)`              | `sqrt(msd(x, y))` |
+|  NormRMSDeviation    |  `nrmsd(x, y)`             | `rmsd(x, y) / (maximum(x) - minimum(x))` |
 |  WeightedEuclidean   |  `weuclidean(x, y, w)`     | `sqrt(sum((x - y).^2 .* w))`  |
 |  WeightedSqEuclidean |  `wsqeuclidean(x, y, w)`   | `sum((x - y).^2 .* w)`  |
 |  WeightedCityblock   |  `wcityblock(x, y, w)`     | `sum(abs(x - y) .* w)`  |

src/Distances.jl

@@ -30,6 +30,7 @@ export
     CorrDist,
     ChiSqDist,
     KLDivergence,
+    GenKLDivergence,
     JSDivergence,
     RenyiDivergence,
     SpanNormDist,
@@ -44,6 +45,11 @@ export
     BhattacharyyaDist,
     HellingerDist,
 
+    MeanAbsDeviation,
+    MeanSqDeviation,
+    RMSDeviation,
+    NormRMSDeviation,
+
     # convenient functions
     euclidean,
     sqeuclidean,
@@ -59,6 +65,7 @@ export
     corr_dist,
     chisq_dist,
     kl_divergence,
+    gkl_divergence,
     js_divergence,
     renyi_divergence,
     spannorm_dist,
@@ -71,7 +78,12 @@ export
     sqmahalanobis,
     mahalanobis,
     bhattacharyya,
-    hellinger
+    hellinger,
+
+    meanad,
+    msd,
+    rmsd,
+    nrmsd
 
 include("common.jl")
 include("generic.jl")

src/metrics.jl

@@ -28,6 +28,8 @@ struct CorrDist <: SemiMetric end
 
 struct ChiSqDist <: SemiMetric end
 struct KLDivergence <: PreMetric end
+struct KLDivergence <: PreMetric end
+struct GenKLDivergence <: PreMetric end
 
 """
     RenyiDivergence(α::Real)
@@ -89,8 +91,15 @@ struct JSDivergence <: SemiMetric end
 
 struct SpanNormDist <: SemiMetric end
 
+# Deviations are handled separately from the other distances/divergences and
+# are excluded from `UnionMetrics`
+struct MeanAbsDeviation <: Metric end
+struct MeanSqDeviation <: SemiMetric end
+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}
+const UnionMetrics = Union{Euclidean, SqEuclidean, Chebyshev, Cityblock, Minkowski, Hamming, Jaccard, RogersTanimoto, CosineDist, CorrDist, ChiSqDist, KLDivergence, RenyiDivergence, JSDivergence, SpanNormDist, GenKLDivergence}
 
 """
     Euclidean([thresh])
@@ -240,6 +249,11 @@ chisq_dist(a::AbstractArray, b::AbstractArray) = evaluate(ChiSqDist(), a, b)
 @inline eval_reduce(::KLDivergence, s1, s2) = s1 + s2
 kl_divergence(a::AbstractArray, b::AbstractArray) = evaluate(KLDivergence(), a, b)
 
+# GenKLDivergence
+@inline eval_op(::GenKLDivergence, ai, bi) = ai > 0 ? ai * log(ai / bi) - ai + bi : bi
+@inline eval_reduce(::GenKLDivergence, s1, s2) = s1 + s2
+gkl_divergence(a::AbstractArray, b::AbstractArray) = evaluate(GenKLDivergence(), a, b)
+
 # RenyiDivergence
 function eval_start(::RenyiDivergence, a::AbstractArray{T}, b::AbstractArray{T}) where {T <: AbstractFloat}
     zero(T), zero(T), sum(a), sum(b)
@@ -369,6 +383,23 @@ end
 end
 rogerstanimoto(a::AbstractArray{T}, b::AbstractArray{T}) where {T <: Bool} = evaluate(RogersTanimoto(), a, b)
 
+# Deviations
+
+evaluate(::MeanAbsDeviation, a, b) = cityblock(a, b) / length(a)
+meanad(a, b) = evaluate(MeanAbsDeviation(), a, b)
+
+evaluate(::MeanSqDeviation, a, b) = sqeuclidean(a, b) / length(a)
+msd(a, b) = evaluate(MeanSqDeviation(), a, b)
+
+evaluate(::RMSDeviation, a, b) = sqrt(evaluate(MeanSqDeviation(), a, b))
+rmsd(a, b) = evaluate(RMSDeviation(), a, b)
+
+function evaluate(::NormRMSDeviation, a, b)
+    amin, amax = extrema(a)
+    return evaluate(RMSDeviation(), a, b) / (amax - amin)
+end
+nrmsd(a, b) = evaluate(NormRMSDeviation(), a, b)
+
 
 ###########################################################
 #

test/test_dists.jl

@@ -134,6 +134,13 @@ end
         @test chisq_dist(x, y) == sum((x - vec(y)).^2 ./ (x + vec(y)))
         @test spannorm_dist(x, y) == maximum(x - vec(y)) - minimum(x - vec(y))
 
+        @test gkl_divergence(x, y) ≈ sum(i -> x[i] * log(x[i] / y[i]) - x[i] + y[i], 1:length(x))
+
+        @test meanad(x, y) ≈ mean(Float64[abs(x[i] - y[i]) for i in 1:length(x)])
+        @test msd(x, y) ≈ mean(Float64[abs2(x[i] - y[i]) for i in 1:length(x)])
+        @test rmsd(x, y) ≈ sqrt(msd(x, y))
+        @test nrmsd(x, y) ≈ sqrt(msd(x, y)) / (maximum(x) - minimum(x))
+
         w = ones(4)
         @test sqeuclidean(x, y) ≈ wsqeuclidean(x, y, w)