From 3fabb54155b6fcb12b855fb18519f22d1f74ffb8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Tue, 28 Oct 2025 18:38:32 +0100 Subject: [PATCH] Remove undesired promotions with `Float64` literals --- src/common.jl | 2 + src/deviation.jl | 88 +++++----- src/moments.jl | 324 +++++++++++++++++++++-------------- src/weights.jl | 5 +- test/deviation.jl | 125 +++++++++----- test/moments.jl | 427 ++++++++++++++++++++++++++++++---------------- 6 files changed, 621 insertions(+), 350 deletions(-) diff --git a/src/common.jl b/src/common.jl index e4d4f712b..4dc268c6c 100644 --- a/src/common.jl +++ b/src/common.jl @@ -9,3 +9,5 @@ function depcheck(fname::Symbol, varname::Symbol, b::Union{Bool, Nothing}) b end end + +_add((x1, x2)::Tuple{<:Real,<:Real}, (y1, y2)::Tuple{<:Real,<:Real}) = (x1 + y1, x2 + y2) diff --git a/src/deviation.jl b/src/deviation.jl index 82d68e839..cfbf96c49 100644 --- a/src/deviation.jl +++ b/src/deviation.jl @@ -10,7 +10,7 @@ Count the number of indices at which the elements of the arrays """ function counteq(a::AbstractArray, b::AbstractArray) n = length(a) - length(b) == n || throw(DimensionMismatch("Inconsistent lengths.")) + length(b) == n || throw(DimensionMismatch("Inconsistent array lengths.")) c = 0 for i in eachindex(a, b) if a[i] == b[i] @@ -29,7 +29,7 @@ Count the number of indices at which the elements of the arrays """ function countne(a::AbstractArray, b::AbstractArray) n = length(a) - length(b) == n || throw(DimensionMismatch("Inconsistent lengths.")) + length(b) == n || throw(DimensionMismatch("Inconsistent array lengths.")) c = 0 for i in eachindex(a, b) if a[i] != b[i] @@ -46,12 +46,14 @@ end Compute the squared L2 distance between two arrays: ``\\sum_{i=1}^n |a_i - b_i|^2``. Efficient equivalent of `sum(abs2, a - b)`. """ -function sqL2dist(a::AbstractArray{T}, b::AbstractArray{T}) where T<:Number +function sqL2dist(a::AbstractArray{<:Number}, b::AbstractArray{<:Number}) n = length(a) - length(b) == n || throw(DimensionMismatch("Input dimension mismatch")) - r = 0.0 - for i in eachindex(a, b) - r += abs2(a[i] - b[i]) + length(b) == n || throw(DimensionMismatch("Inconsistent array lengths.")) + if iszero(n) + r = zero(abs2(zero(eltype(a)) - zero(eltype(b)))) + else + broadcasted = Broadcast.broadcasted((ai, bi) -> abs2(ai - bi), vec(a), vec(b)) + r = sum(Broadcast.instantiate(broadcasted)) end return r end @@ -64,7 +66,7 @@ end Compute the L2 distance between two arrays: ``\\sqrt{\\sum_{i=1}^n |a_i - b_i|^2}``. Efficient equivalent of `sqrt(sum(abs2, a - b))`. """ -L2dist(a::AbstractArray{T}, b::AbstractArray{T}) where {T<:Number} = sqrt(sqL2dist(a, b)) +L2dist(a::AbstractArray{<:Number}, b::AbstractArray{<:Number}) = sqrt(sqL2dist(a, b)) # L1 distance @@ -74,12 +76,14 @@ L2dist(a::AbstractArray{T}, b::AbstractArray{T}) where {T<:Number} = sqrt(sqL2di Compute the L1 distance between two arrays: ``\\sum_{i=1}^n |a_i - b_i|``. Efficient equivalent of `sum(abs, a - b)`. """ -function L1dist(a::AbstractArray{T}, b::AbstractArray{T}) where T<:Number +function L1dist(a::AbstractArray{<:Number}, b::AbstractArray{<:Number}) n = length(a) - length(b) == n || throw(DimensionMismatch("Input dimension mismatch")) - r = 0.0 - for i in eachindex(a, b) - r += abs(a[i] - b[i]) + length(b) == n || throw(DimensionMismatch("Inconsistent array lengths.")) + if iszero(n) + r = zero(abs(zero(eltype(a)) - zero(eltype(b)))) + else + broadcasted = Broadcast.broadcasted((ai, bi) -> abs(ai - bi), vec(a), vec(b)) + r = sum(Broadcast.instantiate(broadcasted)) end return r end @@ -93,15 +97,14 @@ Compute the L∞ distance, also called the Chebyshev distance, between two arrays: ``\\max_{1≤i≤n} |a_i - b_i|``. Efficient equivalent of `maxabs(a - b)`. """ -function Linfdist(a::AbstractArray{T}, b::AbstractArray{T}) where T<:Number +function Linfdist(a::AbstractArray{<:Number}, b::AbstractArray{<:Number}) n = length(a) - length(b) == n || throw(DimensionMismatch("Input dimension mismatch")) - r = 0.0 - for i in eachindex(a, b) - v = abs(a[i] - b[i]) - if r < v - r = v - end + length(b) == n || throw(DimensionMismatch("Inconsistent array lengths.")) + if iszero(n) + r = zero(abs(zero(eltype(a)) - zero(eltype(b)))) + else + broadcasted = Broadcast.broadcasted((ai, bi) -> abs(ai - bi), vec(a), vec(b)) + r = maximum(Broadcast.instantiate(broadcasted)) end return r end @@ -115,19 +118,20 @@ Compute the generalized Kullback-Leibler divergence between two arrays: ``\\sum_{i=1}^n (a_i \\log(a_i/b_i) - a_i + b_i)``. Efficient equivalent of `sum(a*log(a/b)-a+b)`. """ -function gkldiv(a::AbstractArray{T}, b::AbstractArray{T}) where T<:AbstractFloat +function gkldiv(a::AbstractArray{<:Real}, b::AbstractArray{<:Real}) n = length(a) - r = 0.0 - for i in eachindex(a, b) - ai = a[i] - bi = b[i] - if ai > 0 - r += (ai * log(ai / bi) - ai + bi) - else - r += bi + length(b) == n || throw(DimensionMismatch("Inconsistent array lengths.")) + if iszero(n) + za = zero(eltype(a)) + zb = zero(eltype(b)) + r = zero(xlogy(za, za / zb) + (zb - za)) + else + broadcasted = Broadcast.broadcasted(vec(a), vec(b)) do ai, bi + return xlogy(ai, ai / bi) + (bi - ai) end + return sum(Broadcast.instantiate(broadcasted)) end - return r::Float64 + return r end @@ -137,8 +141,7 @@ end Return the mean absolute deviation between two arrays: `mean(abs, a - b)`. """ -meanad(a::AbstractArray{T}, b::AbstractArray{T}) where {T<:Number} = - L1dist(a, b) / length(a) +meanad(a::AbstractArray{<:Number}, b::AbstractArray{<:Number}) = L1dist(a, b) / length(a) # MaxAD: maximum absolute deviation @@ -147,7 +150,7 @@ meanad(a::AbstractArray{T}, b::AbstractArray{T}) where {T<:Number} = Return the maximum absolute deviation between two arrays: `maxabs(a - b)`. """ -maxad(a::AbstractArray{T}, b::AbstractArray{T}) where {T<:Number} = Linfdist(a, b) +maxad(a::AbstractArray{<:Number}, b::AbstractArray{<:Number}) = Linfdist(a, b) # MSD: mean squared deviation @@ -156,8 +159,7 @@ maxad(a::AbstractArray{T}, b::AbstractArray{T}) where {T<:Number} = Linfdist(a, Return the mean squared deviation between two arrays: `mean(abs2, a - b)`. """ -msd(a::AbstractArray{T}, b::AbstractArray{T}) where {T<:Number} = - sqL2dist(a, b) / length(a) +msd(a::AbstractArray{<:Number}, b::AbstractArray{<:Number}) = sqL2dist(a, b) / length(a) # RMSD: root mean squared deviation @@ -168,13 +170,14 @@ Return the root mean squared deviation between two optionally normalized arrays. The root mean squared deviation is computed as `sqrt(msd(a, b))`. """ -function rmsd(a::AbstractArray{T}, b::AbstractArray{T}; normalize::Bool=false) where T<:Number +function rmsd(a::AbstractArray{<:Number}, b::AbstractArray{<:Number}; normalize::Bool=false) v = sqrt(msd(a, b)) if normalize - amin, amax = extrema(a) - v /= (amax - amin) + amin, amax = isempty(a) ? (zero(eltype(a)), zero(eltype(a))) : extrema(a) + return v / (amax - amin) + else + return v end - return v end @@ -186,6 +189,7 @@ Compute the peak signal-to-noise ratio between two arrays `a` and `b`. `maxv` is the maximum possible value either array can take. The PSNR is computed as `10 * log10(maxv^2 / msd(a, b))`. """ -function psnr(a::AbstractArray{T}, b::AbstractArray{T}, maxv::Real) where T<:Real - 20. * log10(maxv) - 10. * log10(msd(a, b)) +function psnr(a::AbstractArray{<:Real}, b::AbstractArray{<:Real}, maxv::Real) + msd_a_b, _maxv = promote(msd(a, b), maxv) + return 20 * log10(_maxv) - 10 * log10(msd_a_b) end diff --git a/src/moments.jl b/src/moments.jl index 9e8f670e8..04de97a83 100644 --- a/src/moments.jl +++ b/src/moments.jl @@ -20,13 +20,18 @@ replacing ``\\frac{1}{\\sum{w}}`` with a factor dependent on the type of weights """ function var(v::AbstractArray{<:Real}, w::AbstractWeights; mean=nothing, corrected::Union{Bool, Nothing}=nothing) + length(w) == length(v) || throw(DimensionMismatch("Inconsistent array lengths.")) corrected = depcheck(:var, :corrected, corrected) - - if mean == nothing - _moment2(v, w, Statistics.mean(v, w); corrected=corrected) - else - _moment2(v, w, mean; corrected=corrected) + if mean === nothing + mean = Statistics.mean(v, w) end + return _moment2(v, w, mean; corrected) +end +function var(v::AbstractArray{<:Real}, w::UnitWeights; mean=nothing, + corrected::Union{Bool, Nothing}=nothing) + length(w) == length(v) || throw(DimensionMismatch("Inconsistent array lengths.")) + corrected = depcheck(:var, :corrected, corrected) + return var(v; mean, corrected) end ## var along dim @@ -58,9 +63,20 @@ end function var(A::AbstractArray{<:Real}, w::AbstractWeights, dim::Int; mean=nothing, corrected::Union{Bool, Nothing}=nothing) corrected = depcheck(:var, :corrected, corrected) - var!(similar(A, Float64, Base.reduced_indices(axes(A), dim)), A, w, dim; + if mean === nothing + z = (zero(eltype(w)) * zero(eltype(A))^2) / zero(eltype(w)) + else + z = (zero(eltype(w)) * zero(zero(eltype(A)) - zero(eltype(mean)))^2) / zero(eltype(w)) + end + var!(similar(A, typeof(z), Base.reduced_indices(axes(A), dim)), A, w, dim; mean=mean, corrected=corrected) end +function var(v::AbstractArray{<:Real}, w::UnitWeights, dim::Int; mean=nothing, + corrected::Union{Bool, Nothing}=nothing) + length(w) == length(v) || throw(DimensionMismatch("Inconsistent array lengths.")) + corrected = depcheck(:var, :corrected, corrected) + return var(v; mean, corrected, dims=dim) +end ## std """ @@ -160,83 +176,130 @@ end ##### General central moment -function _moment2(v::AbstractArray{<:Real}, m::Real; corrected=false) +function _moment2(v::AbstractArray{<:Real}, m::Real; corrected::Bool) n = length(v) - s = 0.0 - for i = 1:n - z = v[i] - m - s += z * z + if iszero(n) + z = zero(zero(eltype(v)) - m) + s = z^2 + else + s = let m = m + sum(v) do vi + zi = vi - m + return zi^2 + end + end end - varcorrection(n, corrected) * s + return oftype(float(s), varcorrection(n, corrected)) * s end -function _moment2(v::AbstractArray{<:Real}, wv::AbstractWeights, m::Real; corrected=false) - n = length(v) - s = 0.0 - for i = 1:n - z = v[i] - m - s += (z * z) * wv[i] +function _moment2(v::AbstractArray{<:Real}, wv::AbstractWeights, m::Real; corrected::Bool) + if isempty(v) + z = zero(zero(eltype(v)) - m) + s = z^2 * zero(eltype(wv)) + else + broadcasted = let m = m + Broadcast.broadcasted(vec(v), wv) do vi, wvi + zi = vi - m + return zi^2 * wvi + end + end + s = sum(Broadcast.instantiate(broadcasted)) end - - varcorrection(wv, corrected) * s + return oftype(float(s), varcorrection(wv, corrected)) * s end function _moment3(v::AbstractArray{<:Real}, m::Real) n = length(v) - s = 0.0 - for i = 1:n - z = v[i] - m - s += z * z * z + if iszero(n) + z = zero(zero(eltype(v)) - m) + s = z^3 + else + s = let m = m + sum(v) do vi + zi = vi - m + return zi^3 + end + end end s / n end function _moment3(v::AbstractArray{<:Real}, wv::AbstractWeights, m::Real) - n = length(v) - s = 0.0 - for i = 1:n - z = v[i] - m - s += (z * z * z) * wv[i] + if isempty(v) + z = zero(zero(eltype(v)) - m) + s = zero(z^3 * zero(eltype(wv))) + else + broadcasted = let m = m + Broadcast.broadcasted(vec(v), wv) do vi, wvi + zi = vi - m + return zi^3 * wvi + end + end + s = sum(Broadcast.instantiate(broadcasted)) end s / sum(wv) end function _moment4(v::AbstractArray{<:Real}, m::Real) n = length(v) - s = 0.0 - for i = 1:n - z = v[i] - m - s += abs2(z * z) + if iszero(n) + z = zero(zero(eltype(v)) - m) + s = zero(z^4) + else + s = let m = m + sum(v) do vi + zi = vi - m + return zi^4 + end + end end s / n end function _moment4(v::AbstractArray{<:Real}, wv::AbstractWeights, m::Real) - n = length(v) - s = 0.0 - for i = 1:n - z = v[i] - m - s += abs2(z * z) * wv[i] + if isempty(v) + z = zero(zero(eltype(v)) - m) + s = zero(z^4 * zero(eltype(wv))) + else + broadcasted = let m = m + Broadcast.broadcasted(vec(v), wv) do vi, wvi + zi = vi - m + return zi^4 * wvi + end + end + s = sum(Broadcast.instantiate(broadcasted)) end s / sum(wv) end function _momentk(v::AbstractArray{<:Real}, k::Int, m::Real) n = length(v) - s = 0.0 - for i = 1:n - z = v[i] - m - s += (z ^ k) + if iszero(n) + z = zero(zero(eltype(v)) - m) + s = zero(z^k) + else + s = let m = m, k = k + sum(v) do vi + zi = vi - m + return zi^k + end + end end s / n end function _momentk(v::AbstractArray{<:Real}, k::Int, wv::AbstractWeights, m::Real) - n = length(v) - s = 0.0 - for i = 1:n - z = v[i] - m - s += (z ^ k) * wv[i] + if isempty(v) + z = zero(zero(eltype(v)) - m) + s = zero(z^k * zero(eltype(wv))) + else + broadcasted = let m = m, k = k + Broadcast.broadcasted(vec(v), wv) do vi, wvi + zi = vi - m + return zi^k * wvi + end + end + s = sum(Broadcast.instantiate(broadcasted)) end s / sum(wv) end @@ -248,26 +311,25 @@ end Return the `k`th order central moment of a real-valued array `v`, optionally specifying a weighting vector `wv` and a center `m`. """ -function moment(v::AbstractArray{<:Real}, k::Int, m::Real) - k == 2 ? _moment2(v, m) : +function moment(v::AbstractArray{<:Real}, k::Int, m::Real=mean(v)) + k == 2 ? _moment2(v, m; corrected = false) : k == 3 ? _moment3(v, m) : k == 4 ? _moment4(v, m) : _momentk(v, k, m) end -function moment(v::AbstractArray{<:Real}, k::Int, wv::AbstractWeights, m::Real) - k == 2 ? _moment2(v, wv, m) : +function moment(v::AbstractArray{<:Real}, k::Int, wv::AbstractWeights, m::Real=mean(v, wv)) + length(wv) == length(v) || throw(DimensionMismatch("Inconsistent array lengths.")) + k == 2 ? _moment2(v, wv, m; corrected = false) : k == 3 ? _moment3(v, wv, m) : k == 4 ? _moment4(v, wv, m) : _momentk(v, k, wv, m) end - -moment(v::AbstractArray{<:Real}, k::Int) = moment(v, k, mean(v)) -function moment(v::AbstractArray{<:Real}, k::Int, wv::AbstractWeights) - moment(v, k, wv, mean(v, wv)) +function moment(v::AbstractArray{<:Real}, k::Int, wv::UnitWeights, m::Real) + length(wv) == length(v) || throw(DimensionMismatch("Inconsistent array lengths.")) + return moment(v, k, m) end - ##### Skewness and Kurtosis # Skewness @@ -278,44 +340,50 @@ end Compute the standardized skewness of a real-valued array `v`, optionally specifying a weighting vector `wv` and a center `m`. """ -function skewness(v::AbstractArray{<:Real}, m::Real) +function skewness(v::AbstractArray{<:Real}, m::Real=mean(v)) n = length(v) - cm2 = 0.0 # empirical 2nd centered moment (variance) - cm3 = 0.0 # empirical 3rd centered moment - for i = 1:n - z = v[i] - m - z2 = z * z - - cm2 += z2 - cm3 += z2 * z + if iszero(n) + z = zero(zero(eltype(v)) - m) + cm2 = z^2 # empirical 2nd centered moment (variance) + cm3 = cm2 * z # empirical 3rd centered moment + else + cm2, cm3 = let m = m + mapreduce(_add, v) do vi + zi = vi - m + z2i = zi^2 + z3i = z2i * zi + return z2i, z3i + end + end end - cm3 /= n - cm2 /= n - return cm3 / sqrt(cm2 * cm2 * cm2) # this is much faster than cm2^1.5 + return (cm3/n) / sqrt((cm2/n)^3) # this is much faster than cm2^1.5 end -function skewness(v::AbstractArray{<:Real}, wv::AbstractWeights, m::Real) +function skewness(v::AbstractArray{<:Real}, wv::AbstractWeights, m::Real=mean(v, wv)) n = length(v) length(wv) == n || throw(DimensionMismatch("Inconsistent array lengths.")) - cm2 = 0.0 # empirical 2nd centered moment (variance) - cm3 = 0.0 # empirical 3rd centered moment - - for i = 1:n - x_i = v[i] - w_i = wv[i] - z = x_i - m - z2w = z * z * w_i - cm2 += z2w - cm3 += z2w * z + if iszero(n) + z = zero(zero(eltype(v)) - m) + cm2 = z^2 * zero(eltype(wv)) # empirical 2nd centered moment (variance) + cm3 = cm2 * z # empirical 3rd centered moment + else + broadcasted = let m = m + Broadcast.broadcasted(vec(v), wv) do vi, wvi + zi = vi - m + z2wi = zi^2 * wvi + z3wi = z2wi * zi + return z2wi, z3wi + end + end + cm2, cm3 = reduce(_add, Broadcast.instantiate(broadcasted)) end sw = sum(wv) - cm3 /= sw - cm2 /= sw - return cm3 / sqrt(cm2 * cm2 * cm2) # this is much faster than cm2^1.5 + return (cm3/sw) / sqrt((cm2/sw)^3) # this is much faster than cm2^1.5 +end +function skewness(v::AbstractArray{<:Real}, wv::UnitWeights, m::Real) + length(wv) == length(v) || throw(DimensionMismatch("Inconsistent array lengths.")) + return skewness(v, m) end - -skewness(v::AbstractArray{<:Real}) = skewness(v, mean(v)) -skewness(v::AbstractArray{<:Real}, wv::AbstractWeights) = skewness(v, wv, mean(v, wv)) # (excessive) Kurtosis # This is Type 1 definition according to Joanes and Gill (1998) @@ -325,44 +393,52 @@ skewness(v::AbstractArray{<:Real}, wv::AbstractWeights) = skewness(v, wv, mean(v Compute the excess kurtosis of a real-valued array `v`, optionally specifying a weighting vector `wv` and a center `m`. """ -function kurtosis(v::AbstractArray{<:Real}, m::Real) +function kurtosis(v::AbstractArray{<:Real}, m::Real=mean(v)) n = length(v) - cm2 = 0.0 # empirical 2nd centered moment (variance) - cm4 = 0.0 # empirical 4th centered moment - for i = 1:n - z = v[i] - m - z2 = z * z - cm2 += z2 - cm4 += z2 * z2 + if iszero(n) + z = zero(zero(eltype(v)) - m) + cm2 = z^2 # empirical 2nd centered moment (variance) + cm4 = cm2^2 # empirical 4th centered moment + else + cm2, cm4 = let m = m + mapreduce(_add, v) do vi + zi = vi - m + z2i = zi^2 + z4i = z2i^2 + return z2i, z4i + end + end end - cm4 /= n - cm2 /= n - return (cm4 / (cm2 * cm2)) - 3.0 + return (cm4/n) / (cm2/n)^2 - 3 end -function kurtosis(v::AbstractArray{<:Real}, wv::AbstractWeights, m::Real) +function kurtosis(v::AbstractArray{<:Real}, wv::AbstractWeights, m::Real=mean(v, wv)) n = length(v) length(wv) == n || throw(DimensionMismatch("Inconsistent array lengths.")) - cm2 = 0.0 # empirical 2nd centered moment (variance) - cm4 = 0.0 # empirical 4th centered moment - - for i = 1 : n - x_i = v[i] - w_i = wv[i] - z = x_i - m - z2 = z * z - z2w = z2 * w_i - cm2 += z2w - cm4 += z2w * z2 + if iszero(n) + z = zero(zero(eltype(v)) - m) + z2 = z^2 + cm2 = z2 * zero(eltype(wv)) # empirical 2nd centered moment (variance) + cm4 = cm2 * z2 # empirical 4th centered moment + else + broadcasted = let m = m + Broadcast.broadcasted(vec(v), wv) do vi, wvi + zi = vi - m + z2i = zi^2 + z2wi = z2i * wvi + z4wi = z2wi * z2i + return z2wi, z4wi + end + end + cm2, cm4 = reduce(_add, Broadcast.instantiate(broadcasted)) end sw = sum(wv) - cm4 /= sw - cm2 /= sw - return (cm4 / (cm2 * cm2)) - 3.0 + return (cm4/sw) / (cm2/sw)^2 - 3 +end +function kurtosis(v::AbstractArray{<:Real}, wv::UnitWeights, m::Real) + length(wv) == length(v) || throw(DimensionMismatch("Inconsistent array lengths.")) + return kurtosis(v, m) end - -kurtosis(v::AbstractArray{<:Real}) = kurtosis(v, mean(v)) -kurtosis(v::AbstractArray{<:Real}, wv::AbstractWeights) = kurtosis(v, wv, mean(v, wv)) """ cumulant(v, k, [wv::AbstractWeights], m=mean(v)) @@ -380,19 +456,19 @@ https://doi.org/10.2307/2684642 """ function cumulant(v::AbstractArray{<:Real}, krange::Union{Integer, AbstractRange{<:Integer}}, wv::AbstractWeights, m::Real=mean(v, wv)) - if minimum(krange) <= 0 + n = length(v) + length(wv) == n || throw(DimensionMismatch("Inconsistent array lengths.")) + kmin, kmax = extrema(krange) + if kmin <= 0 throw(ArgumentError("Cumulant orders must be strictly positive.")) end - k = maximum(krange) - cmoms = zeros(typeof(m), k) - cumls = zeros(typeof(m), k) - cmoms[1] = 0 + cmoms = [moment(v, i, wv, m) for i in 2:kmax] + cumls = Vector{eltype(cmoms)}(undef, kmax) cumls[1] = m - for i = 2:k - kn = wv isa UnitWeights ? moment(v, i, m) : moment(v, i, wv, m) - cmoms[i] = kn - for j = 2:i-2 - kn -= binomial(i-1, j)*cmoms[j]*cumls[i-j] + for i = 2:kmax + kn = cmoms[i-1] + for j = 2:(i-2) + kn -= binomial(i-1, j)*cmoms[j-1]*cumls[i-j] end cumls[i] = kn end diff --git a/src/weights.jl b/src/weights.jl index dabe7d60d..ec374d776 100644 --- a/src/weights.jl +++ b/src/weights.jl @@ -409,7 +409,10 @@ wsum(v::AbstractArray, w::AbstractVector, dims::Colon=:) = transpose(w) * vec(v) # Optimized methods (to ensure we use BLAS when possible) for W in (AnalyticWeights, FrequencyWeights, ProbabilityWeights, Weights) @eval begin - wsum(v::AbstractArray, w::$W, dims::Colon) = transpose(w.values) * vec(v) + function wsum(v::AbstractArray, w::$W, dims::Colon) + length(w) == length(v) || throw(DimensionMismatch("Inconsistent array lengths.")) + return transpose(w.values) * vec(v) + end end end diff --git a/test/deviation.jl b/test/deviation.jl index c01fbb64e..4ace5846d 100644 --- a/test/deviation.jl +++ b/test/deviation.jl @@ -1,48 +1,93 @@ using StatsBase, OffsetArrays using Test -a = [1, 2, 3, 4, 5, 6, 7] -b = [1, 3, 3, 4, 6, 7, 8] +@testset "counting (arrays with element types $T1 and $T2)" for T1 in (Int, Float32, Float64), T2 in (Int, Float32, Float64) + a = T1[1, 2, 3, 4, 5, 6, 7] + b = T2[1, 3, 3, 4, 6, 7, 8] + a_offset = OffsetArray(a, -5:1) + b_offset = OffsetArray(b, -5:1) + for (a, b) in ((a, b), (a_offset, b_offset)) + @test @inferred(counteq(a, b))::Int == 3 + @test @inferred(countne(a, b))::Int == 4 + end -@test counteq(a, b) == 3 -@test countne(a, b) == 4 + # Empty arrays + empty_a = T1[] + empty_b = T2[] + empty_a_offset = OffsetArray(empty_a, -5) + empty_b_offset = OffsetArray(empty_b, -5) + for (a, b) in ((empty_a, empty_b), (empty_a_offset, empty_b_offset)) + @test @inferred(counteq(a, b))::Int == 0 + @test @inferred(countne(a, b))::Int == 0 + end -a = OffsetArray(a, -5:1) -b = OffsetArray(b, -5:1) + # Inconsistent lengths + err = DimensionMismatch("Inconsistent array lengths.") + for (a, b) in ((a, empty_b), (empty_a, b), (a_offset, empty_b_offset), (empty_a_offset, b_offset)) + @test_throws err counteq(a, b) + @test_throws err countne(a, b) + end +end -@test counteq(a, b) == 3 -@test countne(a, b) == 4 +@testset "deviation (arrays with element types $T1 and $T2)" for T1 in (Float32, Float64), T2 in (Float32, Float64) + T = promote_type(T1, T2) + a = rand(T1, 5, 6) + b = rand(T2, 5, 6) + a_offset = OffsetArray(a, 5, -10) + b_offset = OffsetArray(b, 5, -10) + for (a, b) in ((a, b), (a_offset, b_offset)) + @test @inferred(sqL2dist(a, b))::T ≈ sum(abs2.(a - b)) + @test @inferred(L2dist(a, b))::T ≈ sqrt(sqL2dist(a, b)) + @test @inferred(L1dist(a, b))::T ≈ sum(abs.(a - b)) + @test @inferred(Linfdist(a, b))::T ≈ maximum(abs.(a - b)) + @test @inferred(gkldiv(a, b))::T ≈ sum(a .* log.(a ./ b) - a + b) + @test @inferred(meanad(a, b))::T ≈ mean(abs.(a - b)) + @test @inferred(maxad(a, b))::T ≈ maximum(abs.(a - b)) + @test @inferred(msd(a, b))::T ≈ mean(abs2.(a - b)) + @test @inferred(rmsd(a, b))::T ≈ sqrt(msd(a, b)) + @test @inferred(rmsd(a, b; normalize=true))::T ≈ rmsd(a, b) / (maximum(a) - minimum(a)) + for T2 in (Int, Float32, Float64) + S = promote_type(T, T2) + @test @inferred(psnr(a, b, T2(2)))::S ≈ 10 * log10(4 / msd(a, b)) + end + end -a = rand(5, 6) -b = rand(5, 6) + # Empty arrays + empty_a = T1[] + empty_b = T2[] + empty_a_offset = OffsetArray(empty_a, 5) + empty_b_offset = OffsetArray(empty_b, 5) + for (a, b) in ((empty_a, empty_b), (empty_a_offset, empty_b_offset)) + @test iszero(@inferred(sqL2dist(a, b))::T) + @test iszero(@inferred(L2dist(a, b))::T) + @test iszero(@inferred(L1dist(a, b))::T) + @test iszero(@inferred(Linfdist(a, b))::T) + @test iszero(@inferred(gkldiv(a, b))::T) + @test isnan(@inferred(meanad(a, b))::T) + @test iszero(@inferred(maxad(a, b))::T) + @test isnan(@inferred(msd(a, b))::T) + @test isnan(@inferred(rmsd(a, b))::T) + @test isnan(@inferred(rmsd(a, b; normalize=true))::T) + for T2 in (Int, Float32, Float64) + S = promote_type(T, T2) + @test isnan(@inferred(psnr(a, b, T2(2)))::S) + end + end -@test sqL2dist(a, b) ≈ sum(abs2.(a - b)) -@test L2dist(a, b) ≈ sqrt(sqL2dist(a, b)) -@test L1dist(a, b) ≈ sum(abs.(a - b)) -@test Linfdist(a, b) ≈ maximum(abs.(a - b)) - -@test gkldiv(a, b) ≈ sum(a .* log.(a ./ b) - a + b) - -@test meanad(a, b) ≈ mean(abs.(a - b)) -@test maxad(a, b) ≈ maximum(abs.(a - b)) -@test msd(a, b) ≈ mean(abs2.(a - b)) -@test rmsd(a, b) ≈ sqrt(msd(a, b)) -@test rmsd(a, b; normalize=true) ≈ rmsd(a, b) / (maximum(a) - minimum(a)) -@test psnr(a, b, 2) ≈ 10 * log10(4 / msd(a, b)) - -a = OffsetArray(a, 5, -10) -b = OffsetArray(b, 5, -10) - -@test sqL2dist(a, b) ≈ sum(abs2.(a - b)) -@test L2dist(a, b) ≈ sqrt(sqL2dist(a, b)) -@test L1dist(a, b) ≈ sum(abs.(a - b)) -@test Linfdist(a, b) ≈ maximum(abs.(a - b)) - -@test gkldiv(a, b) ≈ sum(a .* log.(a ./ b) - a + b) - -@test meanad(a, b) ≈ mean(abs.(a - b)) -@test maxad(a, b) ≈ maximum(abs.(a - b)) -@test msd(a, b) ≈ mean(abs2.(a - b)) -@test rmsd(a, b) ≈ sqrt(msd(a, b)) -@test rmsd(a, b; normalize=true) ≈ rmsd(a, b) / (maximum(a) - minimum(a)) -@test psnr(a, b, 2) ≈ 10 * log10(4 / msd(a, b)) \ No newline at end of file + err = DimensionMismatch("Inconsistent array lengths.") + for (a, b) in ((a, empty_b), (empty_a, b), (a_offset, empty_b_offset), (empty_a_offset, b_offset)) + @test_throws err sqL2dist(a, b) + @test_throws err L2dist(a, b) + @test_throws err L1dist(a, b) + @test_throws err Linfdist(a, b) + @test_throws err gkldiv(a, b) + @test_throws err meanad(a, b) + @test_throws err maxad(a, b) + @test_throws err msd(a, b) + @test_throws err rmsd(a, b) + @test_throws err rmsd(a, b; normalize=true) + for T2 in (Int, Float32, Float64) + @test_throws err psnr(a, b, T2(2)) + end + end +end diff --git a/test/moments.jl b/test/moments.jl index 55f1c8552..a3402b965 100644 --- a/test/moments.jl +++ b/test/moments.jl @@ -1,4 +1,5 @@ using StatsBase +using Statistics using Test @testset "StatsBase.Moments" begin @@ -7,44 +8,48 @@ weight_funcs = (weights, aweights, fweights, pweights) ##### weighted var & std x = [0.57, 0.10, 0.91, 0.72, 0.46, 0.0] +xf0 = Float32.(x) w = [3.84, 2.70, 8.29, 8.91, 9.71, 0.0] +wf0 = Float32.(w) -@testset "Uncorrected with $f" for f in weight_funcs +@testset "Uncorrected with $f (values of type $(eltype(x)), weights of type $(eltype(w)))" for f in weight_funcs, x in (x, xf0), w in (w, wf0) + TX = eltype(x) + T = promote_type(TX, eltype(w)) wv = f(w) - m = mean(x, wv) + m = @inferred(mean(x, wv))::T # expected uncorrected output expected_var = sum(abs2.(x .- m), wv) / sum(wv) expected_std = sqrt.(expected_var) @testset "Variance" begin - @test var(x, wv; corrected=false) ≈ expected_var - @test var(x, wv; mean=m, corrected=false) ≈ expected_var - @test varm(x, wv, m; corrected=false) ≈ expected_var + @test @inferred(var(x, wv; corrected=false))::T ≈ expected_var + @test @inferred(var(x, wv; mean=m, corrected=false))::T ≈ expected_var + @test @inferred(varm(x, wv, m; corrected=false))::T ≈ expected_var end @testset "Standard Deviation" begin - @test std(x, wv; corrected=false) ≈ expected_std - @test std(x, wv; mean=m, corrected=false) ≈ expected_std - @test stdm(x, wv, m; corrected=false) ≈ expected_std + @test @inferred(std(x, wv; corrected=false))::T ≈ expected_std + @test @inferred(std(x, wv; mean=m, corrected=false))::T ≈ expected_std + @test @inferred(stdm(x, wv, m; corrected=false))::T ≈ expected_std end @testset "Mean and Variance" begin - (m, v) = mean_and_var(x; corrected=false) + (m, v) = @inferred(mean_and_var(x; corrected=false))::Tuple{TX,TX} @test m == mean(x) @test v == var(x; corrected=corrected=false) - (m, v) = mean_and_var(x, wv; corrected=false) + (m, v) = @inferred(mean_and_var(x, wv; corrected=false))::Tuple{T,T} @test m == mean(x, wv) @test v == var(x, wv; corrected=false) end @testset "Mean and Standard Deviation" begin - (m, s) = mean_and_std(x; corrected=false) + (m, s) = @inferred(mean_and_std(x; corrected=false))::Tuple{TX,TX} @test m == mean(x) @test s == std(x; corrected=false) - (m, s) = mean_and_std(x, wv; corrected=false) + (m, s) = @inferred(mean_and_std(x, wv; corrected=false))::Tuple{T,T} @test m == mean(x, wv) @test s == std(x, wv; corrected=false) end @@ -54,17 +59,21 @@ end expected_var = [NaN, 0.0694434191182236, 0.05466601256158146, 0.06628969012045285] expected_std = sqrt.(expected_var) -@testset "Corrected with $(weight_funcs[i])" for i in eachindex(weight_funcs) +@testset "Corrected with $(weight_funcs[i]) (values of type $(eltype(x)), weights of type $(eltype(w)))" for i in eachindex(weight_funcs), x in (x, xf0), w in (w, wf0) + TX = eltype(x) + TW = eltype(w) + T = promote_type(TX, TW) + TR = TX === Float32 || TW === Float32 ? Float32 : Float64 wv = weight_funcs[i](w) - m = mean(x, wv) + m = @inferred(mean(x, wv))::T @testset "Variance" begin if isa(wv, Weights) @test_throws ArgumentError var(x, wv; corrected=true) else - @test var(x, wv; corrected=true) ≈ expected_var[i] - @test var(x, wv; mean=m, corrected=true) ≈ expected_var[i] - @test varm(x, wv, m; corrected=true) ≈ expected_var[i] + @test @inferred(var(x, wv; corrected=true))::T ≈ TR(expected_var[i]) + @test @inferred(var(x, wv; mean=m, corrected=true))::T ≈ TR(expected_var[i]) + @test @inferred(varm(x, wv, m; corrected=true))::T ≈ TR(expected_var[i]) end end @@ -72,35 +81,35 @@ expected_std = sqrt.(expected_var) if isa(wv, Weights) @test_throws ArgumentError std(x, wv; corrected=true) else - @test std(x, wv; corrected=true) ≈ expected_std[i] - @test std(x, wv; mean=m, corrected=true) ≈ expected_std[i] - @test stdm(x, wv, m; corrected=true) ≈ expected_std[i] + @test @inferred(std(x, wv; corrected=true))::T ≈ TR(expected_std[i]) + @test @inferred(std(x, wv; mean=m, corrected=true))::T ≈ TR(expected_std[i]) + @test @inferred(stdm(x, wv, m; corrected=true))::T ≈ TR(expected_std[i]) end end @testset "Mean and Variance" begin - (m, v) = mean_and_var(x; corrected=true) + (m, v) = @inferred(mean_and_var(x; corrected=true))::Tuple{TX,TX} @test m == mean(x) @test v == var(x; corrected=true) if isa(wv, Weights) @test_throws ArgumentError mean_and_var(x, wv; corrected=true) else - (m, v) = mean_and_var(x, wv; corrected=true) + (m, v) = @inferred(mean_and_var(x, wv; corrected=true))::Tuple{T,T} @test m == mean(x, wv) @test v == var(x, wv; corrected=true) end end @testset "Mean and Standard Deviation" begin - (m, s) = mean_and_std(x; corrected=true) + (m, s) = @inferred(mean_and_std(x; corrected=true))::Tuple{TX,TX} @test m == mean(x) @test s == std(x; corrected=true) if isa(wv, Weights) @test_throws ArgumentError mean_and_std(x, wv; corrected=true) else - (m, s) = mean_and_std(x, wv; corrected=true) + (m, s) = @inferred(mean_and_std(x, wv; corrected=true))::Tuple{T,T} @test m == mean(x, wv) @test s == std(x, wv; corrected=true) end @@ -108,14 +117,22 @@ expected_std = sqrt.(expected_var) end x = rand(5, 6) +xf0 = Float32.(x) w1 = [0.57, 5.10, 0.91, 1.72, 0.0] +w1f0 = Float32.(w1) w2 = [3.84, 2.70, 8.29, 8.91, 9.71, 0.0] - -@testset "Uncorrected with $f" for f in weight_funcs +w2f0 = Float32.(w2) + +@testset "Uncorrected with $f (values of type $(eltype(x)), 1st weights of type $(eltype(w1)), 2nd weights of type $(eltype(w2)))" for f in weight_funcs, x in (x, xf0), w1 in (w1, w1f0), w2 in (w2, w2f0) + TX = eltype(x) + TW1 = eltype(w1) + TW2 = eltype(w2) + T1 = promote_type(TX, TW1) + T2 = promote_type(TX, TW2) wv1 = f(w1) wv2 = f(w2) - m1 = mean(x, wv1, dims=1) - m2 = mean(x, wv2, dims=2) + m1 = @inferred(mean(x, wv1, dims=1))::Matrix{T1} + m2 = @inferred(mean(x, wv2, dims=2))::Matrix{T2} expected_var1 = sum(abs2.(x .- m1) .* w1, dims = 1) ./ sum(wv1) expected_var2 = sum(abs2.(x .- m2) .* w2', dims = 2) ./ sum(wv2) @@ -123,61 +140,66 @@ w2 = [3.84, 2.70, 8.29, 8.91, 9.71, 0.0] expected_std2 = sqrt.(expected_var2) @testset "Variance" begin - @test var(x, wv1, 1; corrected=false) ≈ expected_var1 - @test var(x, wv2, 2; corrected=false) ≈ expected_var2 - @test var(x, wv1, 1; mean=m1, corrected=false) ≈ expected_var1 - @test var(x, wv2, 2; mean=m2, corrected=false) ≈ expected_var2 - @test varm(x, wv1, m1, 1; corrected=false) ≈ expected_var1 - @test varm(x, wv2, m2, 2; corrected=false) ≈ expected_var2 + @test @inferred(var(x, wv1, 1; corrected=false))::Matrix{T1} ≈ expected_var1 + @test @inferred(var(x, wv2, 2; corrected=false))::Matrix{T2} ≈ expected_var2 + @test @inferred(var(x, wv1, 1; mean=m1, corrected=false))::Matrix{T1} ≈ expected_var1 + @test @inferred(var(x, wv2, 2; mean=m2, corrected=false))::Matrix{T2} ≈ expected_var2 + @test @inferred(varm(x, wv1, m1, 1; corrected=false))::Matrix{T1} ≈ expected_var1 + @test @inferred(varm(x, wv2, m2, 2; corrected=false))::Matrix{T2} ≈ expected_var2 end @testset "Standard Deviation" begin - @test std(x, wv1, 1; corrected=false) ≈ expected_std1 - @test std(x, wv2, 2; corrected=false) ≈ expected_std2 - @test std(x, wv1, 1; mean=m1, corrected=false) ≈ expected_std1 - @test std(x, wv2, 2; mean=m2, corrected=false) ≈ expected_std2 - @test stdm(x, wv1, m1, 1; corrected=false) ≈ expected_std1 - @test stdm(x, wv2, m2, 2; corrected=false) ≈ expected_std2 + @test @inferred(std(x, wv1, 1; corrected=false))::Matrix{T1} ≈ expected_std1 + @test @inferred(std(x, wv2, 2; corrected=false))::Matrix{T2} ≈ expected_std2 + @test @inferred(std(x, wv1, 1; mean=m1, corrected=false))::Matrix{T1} ≈ expected_std1 + @test @inferred(std(x, wv2, 2; mean=m2, corrected=false))::Matrix{T2} ≈ expected_std2 + @test @inferred(stdm(x, wv1, m1, 1; corrected=false))::Matrix{T1} ≈ expected_std1 + @test @inferred(stdm(x, wv2, m2, 2; corrected=false))::Matrix{T2} ≈ expected_std2 end @testset "Mean and Variance" begin for d in 1:2 - (m, v) = mean_and_var(x, d; corrected=false) + (m, v) = @inferred(mean_and_var(x, d; corrected=false))::Tuple{Matrix{TX},Matrix{TX}} @test m == mean(x, dims=d) @test v == var(x, dims=d, corrected=false) end - (m, v) = mean_and_var(x, wv1, 1; corrected=false) + (m, v) = @inferred(mean_and_var(x, wv1, 1; corrected=false))::Tuple{Matrix{T1},Matrix{T1}} @test m == mean(x, wv1, dims=1) @test v == var(x, wv1, 1; corrected=false) - (m, v) = mean_and_var(x, wv2, 2; corrected=false) + (m, v) = @inferred(mean_and_var(x, wv2, 2; corrected=false))::Tuple{Matrix{T2},Matrix{T2}} @test m == mean(x, wv2, dims=2) @test v == var(x, wv2, 2; corrected=false) end @testset "Mean and Standard Deviation" begin for d in 1:2 - (m, s) = mean_and_std(x, d; corrected=false) + (m, s) = @inferred(mean_and_std(x, d; corrected=false))::Tuple{Matrix{TX},Matrix{TX}} @test m == mean(x, dims=d) @test s == std(x, dims=d; corrected=false) end - (m, s) = mean_and_std(x, wv1, 1; corrected=false) + (m, s) = @inferred(mean_and_std(x, wv1, 1; corrected=false))::Tuple{Matrix{T1},Matrix{T1}} @test m == mean(x, wv1, dims=1) @test s == std(x, wv1, 1; corrected=false) - (m, s) = mean_and_std(x, wv2, 2; corrected=false) + (m, s) = @inferred(mean_and_std(x, wv2, 2; corrected=false))::Tuple{Matrix{T2},Matrix{T2}} @test m == mean(x, wv2, dims=2) @test s == std(x, wv2, 2; corrected=false) end end -@testset "Corrected with $f" for f in weight_funcs +@testset "Corrected with $f (values of type $(eltype(x)), weights of type $(eltype(w1)))" for f in weight_funcs, x in (Float32.(x), Float64.(x)), (w1, w2) in ((Float32.(w1), Float32.(w2)), (Float64.(w1), Float64.(w2))) + TX = eltype(x) + TW1 = eltype(w1) + TW2 = eltype(w2) + T1 = promote_type(TX, TW1) + T2 = promote_type(TX, TW2) wv1 = f(w1) wv2 = f(w2) - m1 = mean(x, wv1, dims=1) - m2 = mean(x, wv2, dims=2) + m1 = @inferred(mean(x, wv1, dims=1))::Matrix{T1} + m2 = @inferred(mean(x, wv2, dims=2))::Matrix{T2} if !isa(wv1, Weights) expected_var1 = sum(abs2.(x .- m1) .* w1, dims = 1) .* StatsBase.varcorrection(wv1, true) @@ -190,12 +212,12 @@ end if isa(wv1, Weights) @test_throws ArgumentError var(x, wv1, 1; corrected=true) else - @test var(x, wv1, 1; corrected=true) ≈ expected_var1 - @test var(x, wv2, 2; corrected=true) ≈ expected_var2 - @test var(x, wv1, 1; mean=m1, corrected=true) ≈ expected_var1 - @test var(x, wv2, 2; mean=m2, corrected=true) ≈ expected_var2 - @test varm(x, wv1, m1, 1; corrected=true) ≈ expected_var1 - @test varm(x, wv2, m2, 2; corrected=true) ≈ expected_var2 + @test @inferred(var(x, wv1, 1; corrected=true))::Matrix{T1} ≈ expected_var1 + @test @inferred(var(x, wv2, 2; corrected=true))::Matrix{T2} ≈ expected_var2 + @test @inferred(var(x, wv1, 1; mean=m1, corrected=true))::Matrix{T1} ≈ expected_var1 + @test @inferred(var(x, wv2, 2; mean=m2, corrected=true))::Matrix{T2} ≈ expected_var2 + @test @inferred(varm(x, wv1, m1, 1; corrected=true))::Matrix{T1} ≈ expected_var1 + @test @inferred(varm(x, wv2, m2, 2; corrected=true))::Matrix{T2} ≈ expected_var2 end end @@ -203,18 +225,18 @@ end if isa(wv1, Weights) @test_throws ArgumentError std(x, wv1, 1; corrected=true) else - @test std(x, wv1, 1; corrected=true) ≈ expected_std1 - @test std(x, wv2, 2; corrected=true) ≈ expected_std2 - @test std(x, wv1, 1; mean=m1, corrected=true) ≈ expected_std1 - @test std(x, wv2, 2; mean=m2, corrected=true) ≈ expected_std2 - @test stdm(x, wv1, m1, 1; corrected=true) ≈ expected_std1 - @test stdm(x, wv2, m2, 2; corrected=true) ≈ expected_std2 + @test @inferred(std(x, wv1, 1; corrected=true))::Matrix{T1} ≈ expected_std1 + @test @inferred(std(x, wv2, 2; corrected=true))::Matrix{T2} ≈ expected_std2 + @test @inferred(std(x, wv1, 1; mean=m1, corrected=true))::Matrix{T1} ≈ expected_std1 + @test @inferred(std(x, wv2, 2; mean=m2, corrected=true))::Matrix{T2} ≈ expected_std2 + @test @inferred(stdm(x, wv1, m1, 1; corrected=true))::Matrix{T1} ≈ expected_std1 + @test @inferred(stdm(x, wv2, m2, 2; corrected=true))::Matrix{T2} ≈ expected_std2 end end @testset "Mean and Variance" begin for d in 1:2 - (m, v) = mean_and_var(x, d; corrected=true) + (m, v) = @inferred(mean_and_var(x, d; corrected=true))::Tuple{Matrix{TX},Matrix{TX}} @test m == mean(x, dims=d) @test v == var(x, dims=d, corrected=true) end @@ -222,11 +244,11 @@ end if isa(wv1, Weights) @test_throws ArgumentError mean_and_var(x, wv1, 1; corrected=true) else - (m, v) = mean_and_var(x, wv1, 1; corrected=true) + (m, v) = @inferred(mean_and_var(x, wv1, 1; corrected=true))::Tuple{Matrix{T1},Matrix{T1}} @test m == mean(x, wv1, dims=1) @test v == var(x, wv1, 1; corrected=true) - (m, v) = mean_and_var(x, wv2, 2; corrected=true) + (m, v) = @inferred(mean_and_var(x, wv2, 2; corrected=true))::Tuple{Matrix{T2},Matrix{T2}} @test m == mean(x, wv2, dims=2) @test v == var(x, wv2, 2; corrected=true) end @@ -234,7 +256,7 @@ end @testset "Mean and Standard Deviation" begin for d in 1:2 - (m, s) = mean_and_std(x, d; corrected=true) + (m, s) = @inferred(mean_and_std(x, d; corrected=true))::Tuple{Matrix{TX},Matrix{TX}} @test m == mean(x, dims=d) @test s == std(x, dims=d, corrected=true) end @@ -242,11 +264,11 @@ end if isa(wv1, Weights) @test_throws ArgumentError mean_and_std(x, wv1, 1; corrected=true) else - (m, s) = mean_and_std(x, wv1, 1; corrected=true) + (m, s) = @inferred(mean_and_std(x, wv1, 1; corrected=true))::Tuple{Matrix{T1},Matrix{T1}} @test m == mean(x, wv1, dims=1) @test s == std(x, wv1, 1; corrected=true) - (m, s) = mean_and_std(x, wv2, 2; corrected=true) + (m, s) = @inferred(mean_and_std(x, wv2, 2; corrected=true))::Tuple{Matrix{T2},Matrix{T2}} @test m == mean(x, wv2, dims=2) @test s == std(x, wv2, 2; corrected=true) end @@ -254,90 +276,209 @@ end end @testset "Skewness and Kurtosis with $f" for f in weight_funcs - wv = f(ones(5) * 2.0) - - @test skewness(1:5) ≈ 0.0 - @test skewness([1, 2, 3, 4, 5]) ≈ 0.0 - @test skewness([1, 2, 2, 2, 5]) ≈ 1.1731251294063556 - @test skewness([1, 4, 4, 4, 5]) ≈ -1.1731251294063556 - - @test skewness([1, 2, 2, 2, 5], wv) ≈ 1.1731251294063556 + for T in (Int, Float32, Float64) + for v in (T(1):T(5), collect(T, 1:5)) + s = @inferred(skewness(v)) + @test s isa float(T) + @test iszero(s) + + k = @inferred(kurtosis(v)) + @test k isa float(T) + @test k ≈ oftype(k, -1.3) + end - @test kurtosis(1:5) ≈ -1.3 - @test kurtosis([1, 2, 3, 4, 5]) ≈ -1.3 - @test kurtosis([1, 2, 3, 3, 2]) ≈ -1.1530612244897953 + v = T[1, 2, 2, 2, 5] + s = @inferred(skewness(v)) + @test s isa float(T) + @test s ≈ oftype(s, 1.1731251294063556) + + v = T[1, 4, 4, 4, 5] + s = @inferred(skewness(v)) + @test s isa float(T) + @test s ≈ oftype(s, -1.1731251294063556) + + v = T[1, 2, 3, 3, 2] + k = @inferred(kurtosis(v)) + @test k isa float(T) + @test k ≈ oftype(k, -1.1530612244897953) + + # Empty arrays + s = @inferred(skewness(T[])) + @test s isa float(T) + @test isnan(s) + k = @inferred(kurtosis(T[])) + @test k isa float(T) + @test isnan(k) + + for T2 in (Int, Float32, Float64) + wv = f(fill(T2(2), 5)) + v = T[1, 2, 2, 2, 5] + s = @inferred(skewness(v, wv)) + @test s isa float(promote_type(T, T2)) + @test s ≈ oftype(s, 1.1731251294063556) + + v = collect(T, 1:5) + k = @inferred(kurtosis(v, wv)) + @test k isa float(promote_type(T, T2)) + @test k ≈ oftype(k, -1.3) + + # Empty arrays + wv = f(T2[]) + s = @inferred(skewness(T[], wv)) + @test s isa float(promote_type(T, T2)) + @test isnan(s) + k = @inferred(kurtosis(T[], wv)) + @test k isa float(promote_type(T, T2)) + @test isnan(k) + end - @test kurtosis([1, 2, 3, 4, 5], wv) ≈ -1.3 + # Invalid arguments + v = collect(T, 1:5) + for n in (length(x) - 1, length(x) + 1) + @test_throws DimensionMismatch("Inconsistent array lengths.") kurtosis(v, f(ones(T, n))) + @test_throws DimensionMismatch("Inconsistent array lengths.") skewness(v, f(ones(T, n))) + end + end end @testset "General Moments with $f" for f in weight_funcs - x = collect(2.0:8.0) - @test moment(x, 2) ≈ sum((x .- 5).^2) / length(x) - @test moment(x, 3) ≈ sum((x .- 5).^3) / length(x) - @test moment(x, 4) ≈ sum((x .- 5).^4) / length(x) - @test moment(x, 5) ≈ sum((x .- 5).^5) / length(x) - - @test moment(x, 2, 4.0) ≈ sum((x .- 4).^2) / length(x) - @test moment(x, 3, 4.0) ≈ sum((x .- 4).^3) / length(x) - @test moment(x, 4, 4.0) ≈ sum((x .- 4).^4) / length(x) - @test moment(x, 5, 4.0) ≈ sum((x .- 4).^5) / length(x) - - w = f([1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0]) - x2 = collect(2.0:6.0) - @test moment(x, 2, w) ≈ sum((x2 .- 4).^2) / 5 - @test moment(x, 3, w) ≈ sum((x2 .- 4).^3) / 5 - @test moment(x, 4, w) ≈ sum((x2 .- 4).^4) / 5 - @test moment(x, 5, w) ≈ sum((x2 .- 4).^5) / 5 + for T in (Int, Float32, Float64) + x = collect(T, 2:8) + for k in 2:5 + momk = @inferred(moment(x, k)) + @test momk isa float(T) + @test momk ≈ sum((x .- 5).^k) / length(x) + + # Empty array + momk = @inferred(moment(T[], k)) + @test momk isa float(T) + @test isnan(momk) + + for TM in (Int, Float32, Float64) + m = TM(4) + momk = @inferred(moment(x, k, m)) + @test momk isa float(promote_type(T, TM)) + @test momk ≈ sum((x .- 4).^k) / length(x) + + # Empty array + momk = @inferred(moment(T[], k, zero(TM))) + @test momk isa float(promote_type(T, TM)) + @test isnan(momk) + end + end + + for T2 in (Int, Float32, Float64) + wv = f(T2[1, 1, 1, 1, 1, 0, 0]) + x2 = collect(T, 2:6) + for k in 2:5 + momk = @inferred(moment(x, k, wv)) + @test momk isa float(promote_type(T, T2)) + @test momk ≈ sum((x2 .- 4).^k) / 5 + + # Empty array + momk = @inferred(moment(T[], k, f(T2[]))) + @test momk isa float(promote_type(T, T2)) + @test isnan(momk) + + for TM in (Int, Float32, Float64) + m = TM(3) + momk = @inferred(moment(x, k, wv, m)) + @test momk isa float(promote_type(T, T2, TM)) + @test momk ≈ sum((x2 .- 3).^k) / 5 + + # Empty array + momk = @inferred(moment(T[], k, f(T2[]), zero(TM))) + @test momk isa float(promote_type(T, T2, TM)) + @test isnan(momk) + end + end + end + end end @testset "Cumulants with $f" for f in weight_funcs - x = collect(2.0:8.0) - @test cumulant(x, 2) ≈ moment(x, 2) - @test cumulant(x, 3) ≈ moment(x, 3) - @test cumulant(x, 4) ≈ moment(x, 4) - 3*moment(x, 2)^2 - @test cumulant(x, 5) ≈ moment(x, 5) - 10*moment(x, 3)*moment(x, 2) - @test cumulant(x, 6) ≈ - moment(x, 6) - 15*moment(x, 4)*moment(x, 2) - 10*moment(x, 3)^2 + 30*moment(x, 2)^3 - - @test cumulant(x, 1:6) == [cumulant(x, i) for i in 1:6] - - @test cumulant(x, 2, 4.0) ≈ moment(x, 2, 4.0) - @test cumulant(x, 3, 4.0) ≈ moment(x, 3, 4.0) - @test cumulant(x, 4, 4.0) ≈ moment(x, 4, 4.0) - 3*moment(x, 2, 4.0)^2 - @test cumulant(x, 5, 4.0) ≈ moment(x, 5, 4.0) - 10*moment(x, 3, 4.0)*moment(x,2, 4.0) - @test cumulant(x, 6, 4.0) ≈ - moment(x, 6, 4.0) - 15*moment(x, 4, 4.0)*moment(x,2, 4.0) - - 10*moment(x, 3, 4.0)^2 + 30*moment(x, 2, 4.0)^3 - - @test cumulant(x, 1:6, 4.0) == [cumulant(x, i, 4.0) for i in 1:6] - - w1 = f([1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0]) - x2 = collect(2.0:6.0) - @test cumulant(x, 2, w) ≈ moment(x2, 2) - @test cumulant(x, 3, w) ≈ moment(x2, 3) - @test cumulant(x, 4, w) ≈ moment(x2, 4) - 3*moment(x2, 2)^2 - @test cumulant(x, 5, w) ≈ moment(x2, 5) - 10*moment(x2,3)*moment(x2,2) - @test cumulant(x, 6, w) ≈ - moment(x2, 6) - 15*moment(x2,4)*moment(x2,2) + 10*moment(x2,3)^2 + 30*moment(x2,2)^3 - - @test cumulant(x, 1:6, w) == [cumulant(x2, i) for i in 1:6] - - x3 = collect(2:8) - @test cumulant(x3, 6) ≈ - moment(x3, 6) - 15*moment(x3, 4)*moment(x3, 2) - 10*moment(x3, 3)^2 + - 30*moment(x3, 2)^3 - - w2 = f([1, 1, 1, 1, 1, 0, 0]) - x4 = collect(2:6) - @test cumulant(x3, 6, w) ≈ - moment(x4, 6) - 15*moment(x4, 4)*moment(x4, 2) + - 10*moment(x4, 3)^2 + 30*moment(x4, 2)^3 - - @test_throws ArgumentError cumulant(x, -1) - @test_throws ArgumentError cumulant(x, 0) - @test_throws ArgumentError cumulant(x, 0:3) - @test_throws ArgumentError cumulant(x, -1:3) - @test_throws ArgumentError cumulant(x, 1:0) + for T in (Int, Float32, Float64) + x = collect(T, 2:8) + for k in 1:6 + cumk = @inferred(cumulant(x, k)) + @test cumk isa float(T) + if k == 1 + @test cumk ≈ mean(x) + elseif k == 2 || k == 3 + @test cumk ≈ moment(x, k) + elseif k == 4 + @test cumk ≈ moment(x, 4) - 3*moment(x, 2)^2 + elseif k == 5 + @test cumk ≈ moment(x, 5) - 10*moment(x, 3)*moment(x, 2) + else + @assert k == 6 + @test cumk ≈ moment(x, 6) - 15*moment(x, 4)*moment(x, 2) - 10*moment(x, 3)^2 + 30*moment(x, 2)^3 + end + end + cumks = @inferred(cumulant(x, 1:6)) + @test cumks isa Vector{float(T)} + @test cumks == [cumulant(x, i) for i in 1:6] + + for TM in (Int, Float32, Float64) + m = TM(4) + for k in 1:6 + cumk = @inferred(cumulant(x, k, m)) + @test cumk isa float(promote_type(T, TM)) + if k == 1 + @test cumk ≈ m + elseif k == 2 || k == 3 + @test cumk ≈ moment(x, k, m) + elseif k == 4 + @test cumk ≈ moment(x, 4, m) - 3*moment(x, 2, m)^2 + elseif k == 5 + @test cumk ≈ moment(x, 5, m) - 10*moment(x, 3, m)*moment(x, 2, m) + else + @assert k == 6 + @test cumk ≈ moment(x, 6, m) - 15*moment(x, 4, m)*moment(x, 2, m) - 10*moment(x, 3, m)^2 + 30*moment(x, 2, m)^3 + end + end + cumks = @inferred(cumulant(x, 1:6, m)) + @test cumks isa Vector{float(promote_type(T, TM))} + @test cumks == [cumulant(x, i, m) for i in 1:6] + end + + for T2 in (Int, Float32, Float64) + wv = f(T2[1, 1, 1, 1, 1, 0, 0]) + x2 = collect(T, 2:6) + for k in 1:6 + cumk = @inferred(cumulant(x, k, wv)) + @test cumk isa float(promote_type(T, T2)) + @test cumk ≈ cumulant(x2, k) rtol = cbrt(eps(typeof(cumk))) + end + cumks = @inferred(cumulant(x, 1:6, wv)) + @test cumks isa Vector{float(promote_type(T, T2))} + @test cumks == [cumulant(x, i, wv) for i in 1:6] + + for TM in (Int, Float32, Float64) + m = TM(3) + for k in 1:6 + cumk = @inferred(cumulant(x, k, wv, m)) + @test cumk isa float(promote_type(T, T2, TM)) + @test cumk ≈ cumulant(x2, k, m) rtol = cbrt(eps(typeof(cumk))) + end + cumks = @inferred(cumulant(x, 1:6, wv, m)) + @test cumks isa Vector{float(promote_type(T, T2, TM))} + @test cumks == [cumulant(x, i, wv, m) for i in 1:6] + end + end + + # Invalid arguments + @test_throws ArgumentError cumulant(x, -1) + @test_throws ArgumentError cumulant(x, 0) + @test_throws ArgumentError cumulant(x, 0:3) + @test_throws ArgumentError cumulant(x, -1:3) + @test_throws ArgumentError cumulant(x, 1:0) + + for n in (length(x) - 1, length(x) + 1), krange in (1, 1:3) + @test_throws DimensionMismatch("Inconsistent array lengths.") cumulant(x, krange, f(ones(n))) + @test_throws DimensionMismatch("Inconsistent array lengths.") cumulant(x, krange, f(ones(n)), 0.0) + end + end end end # @testset "StatsBase.Moments"