Merge pull request #839 from LilithHafner/patch-2

Combine and optimize `mad!` and `mad`

Author: mschauer

Project.toml

@@ -26,11 +26,12 @@ StatsAPI = "1.2"
 julia = "1"
 
 [extras]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Dates", "DelimitedFiles", "OffsetArrays", "StableRNGs", "Test"]
+test = ["BenchmarkTools", "Dates", "DelimitedFiles", "OffsetArrays", "StableRNGs", "Test"]

src/scalarstats.jl

@@ -532,41 +532,14 @@ If `normalize` is set to `true`, the MAD is multiplied by
 of the standard deviation under the assumption that the data is normally distributed.
 """
 function mad(x; center=nothing, normalize::Union{Bool, Nothing}=nothing, constant=nothing)
-    isempty(x) && throw(ArgumentError("mad is not defined for empty arrays"))
-    T = eltype(x)
-    # Knowing the eltype allows allocating a single array able to hold both original values
-    # and differences from the center, instead of two arrays
-    S = isconcretetype(T) ? promote_type(T, typeof(middle(zero(T)))) : T
-    x2 = x isa AbstractArray ? copyto!(similar(x, S), x) : collect(S, x)
-    c = center === nothing ? median!(x2) : center
-    if isconcretetype(T)
-        x2 .= abs.(x2 .- c)
-    else
-        x2 = abs.(x2 .- c)
-    end
-    m = median!(x2)
-    if normalize isa Nothing
-        Base.depwarn("the `normalize` keyword argument will be false by default in future releases: set it explicitly to silence this deprecation", :mad)
-        normalize = true
-    end
-    if !isa(constant, Nothing)
-        Base.depwarn("keyword argument `constant` is deprecated, use `normalize` instead or apply the multiplication directly", :mad)
-        m * constant
-    elseif normalize
-        m * mad_constant
-    else
-        m
-    end
+    mad!(Base.copymutable(x); center=center, normalize=normalize, constant=constant)
 end
 
 """
     StatsBase.mad!(x; center=median!(x), normalize=true)
 
 Compute the median absolute deviation (MAD) of array `x` around `center`
 (by default, around the median), overwriting `x` in the process.
-`x` must be able to hold values of generated by calling `middle` on its elements
-(for example an integer vector is not appropriate since `middle` can produce
-non-integer values).
 
 If `normalize` is set to `true`, the MAD is multiplied by
 `1 / quantile(Normal(), 3/4) ≈ 1.4826`, in order to obtain a consistent estimator
@@ -577,8 +550,12 @@ function mad!(x::AbstractArray;
               normalize::Union{Bool,Nothing}=true,
               constant=nothing)
     isempty(x) && throw(ArgumentError("mad is not defined for empty arrays"))
-    x .= abs.(x .- center)
-    m = median!(x)
+    c = center === nothing ? median!(x) : center
+    T = promote_type(typeof(c), eltype(x))
+    U = eltype(x)
+    x2 = U == T ? x : isconcretetype(U) && isconcretetype(T) && sizeof(U) == sizeof(T) ? reinterpret(T, x) : similar(x, T)
+    x2 .= abs.(x .- c)
+    m = median!(x2)
     if normalize isa Nothing
         Base.depwarn("the `normalize` keyword argument will be false by default in future releases: set it explicitly to silence this deprecation", :mad)
         normalize = true

test/scalarstats.jl

@@ -2,6 +2,7 @@ using StatsBase
 using Test
 using DelimitedFiles
 using Statistics
+using BenchmarkTools: @benchmark
 
 ##### Location
 
@@ -212,6 +213,9 @@ x = sort!(vcat([5:-1:i for i in 1:5]...))
 @test mad(Any[1, 2.1], normalize=false) ≈ 0.55
 @test mad(Union{Int,Missing}[1, 2], normalize=false) ≈ 0.5
 @test_throws ArgumentError mad(Int[], normalize = true)
+@test mad(Iterators.repeated(4, 10)) == 0
+@test mad(Integer[1,2,3,4]) === mad(1:4)
+@test (@benchmark mad((i for i in 1:10000))).allocs < 200
 
 # Issue 197
 @test mad(1:2, normalize=true) ≈ 0.7413011092528009