Autodot (#5)

* rewrite _autodot and check_lags

* swap StatsBase dependency for Statistics

* update compat

* bump patch version

* remove Statistics compat for julia 1.0 compatibility

Project.toml

@@ -1,20 +1,19 @@
 name = "MeanSquaredDisplacement"
 uuid = "13c93d70-909c-440c-af92-39d48ffa2ba2"
 authors = ["Riccardo Foffi <riccardo.foffi@protonmail.com>"]
-version = "0.1.0"
+version = "0.1.1"
 
 [deps]
 AxisArrays = "39de3d68-74b9-583c-8d2d-e117c070f3a9"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
 AxisArrays = "0.4"
 DSP = "0.7"
 OffsetArrays = "1"
-StatsBase = "0.33, 0.34"
 julia = "1"
 
 [extras]

src/MeanSquaredDisplacement.jl

@@ -4,7 +4,7 @@ using AxisArrays
 using LinearAlgebra: dot
 using DSP
 using OffsetArrays
-using StatsBase
+using Statistics: mean
 
 export emsd, imsd, unfold!
 
@@ -38,7 +38,7 @@ Return the ensemble average of the mean squared displacement of each column
 of `x` at lag times `lags`.
 If not specified `lags` defaults to `0:size(x,1)-1`.
 """
-function emsd(x::AbstractMatrix, lags::AbstractVector{<:Integer}=0:size(x,1)-1)
+function emsd(x::AbstractMatrix, lags=0:size(x,1)-1)
     vec(mean(imsd(x, lags), dims=2))
 end
 
@@ -49,16 +49,16 @@ end
 Return the mean squared displacement of each column of `x` at lag times `lags`.
 If not specified `lags` defaults to `0:size(x,1)-1`.
 """
-function imsd(x::AbstractMatrix, lags::AbstractVector{<:Integer}=0:size(x,1)-1)
+function imsd(x::AbstractMatrix, lags=0:size(x,1)-1)
     mapslices(y -> imsd(y, lags), x, dims=1)
 end
 
 """
     imsd(x::AbstractVector, [lags])
 Return the mean squared displacement of `x` at lag times `lags`.
-If not specified `lags` defaults to `0:length(x)-1`.
+If not specified `lags` defaults to `0:size(x,1)-1`.
 """
-function imsd(x::AbstractVector, lags::AbstractVector{<:Integer}=0:size(x,1)-1)
+function imsd(x::AbstractVector, lags=0:size(x,1)-1)
     l = length(x)
     S₂ = acf(x, lags)
     D = OffsetArray([0.0; dot.(x,x); 0.0], -1:l)
@@ -75,6 +75,7 @@ end
 
 
 include("acf.jl")
+include("autodot.jl")
 include("unfold.jl")
 
 end # module MeanSquaredDisplacement

src/acf.jl

@@ -16,7 +16,7 @@ function fftacf!(r::AbstractVector, x::AbstractVector, lags)
     lx = length(x)
     m = length(lags)
     @assert length(r) == m
-    StatsBase.check_lags(lx, lags)
+    check_lags(lx, lags)
     A = conv(x, reverse(x))
     for k in 1:m
         δ = lags[k]
@@ -34,7 +34,7 @@ function fftacf!(r::AbstractVector, x::AbstractVector{T}, lags) where {T<:Union{
     lx = length(x)
     m = length(lags)
     @assert length(r) == m
-    StatsBase.check_lags(lx, lags)
+    check_lags(lx, lags)
     y = [getindex.(x,i) for i in eachindex(first(x))]
     A = sum([conv(s, reverse(s)) for s in y])
     for k in 1:m
@@ -44,6 +44,11 @@ function fftacf!(r::AbstractVector, x::AbstractVector{T}, lags) where {T<:Union{
     return r
 end
 
+check_lags(lx::Int, lags::AbstractVector) = (
+    maximum(lags) < lx ||
+    error("lags must be less than the sample length.")
+)
+
 #== Autocovariance function with autodot ==#
 # We follow the structure of `StatsBase.autocov` but
 # 1. allow for non-scalar timeseries
@@ -63,17 +68,10 @@ function rawacf!(r::AbstractVector, x::AbstractVector, lags)
     lx = length(x)
     m = length(lags)
     @assert length(r) == m
-    StatsBase.check_lags(lx, lags)
+    check_lags(lx, lags)
     for k in 1:m
         δ = lags[k]
-        r[k] = StatsBase._autodot(x, lx, δ) / (lx-δ)
+        r[k] = _autodot(x, lx, δ) / (lx-δ)
     end
     return r
 end
-
-# Extend to accept non-scalar timeseries
-StatsBase._autodot(x::AbstractVector{T}, lx::Int, l::Int) where {T<:Union{AbstractVector,NTuple}} =
-    dot(view(x, 1:(lx-l)), view(x, (1+l):lx))
-
-
-

src/autodot.jl

@@ -0,0 +1,4 @@
+_autodot(x::AbstractVector{<:Union{Float32, Float64}}, lx::Int, l::Int) = dot(x, 1:(lx-l), x, (1+l):lx)
+_autodot(x::AbstractVector{<:Real}, lx::Int, l::Int) = dot(view(x, 1:(lx-l)), view(x, (1+l):lx))
+_autodot(x::AbstractVector{<:AbstractVector}, lx::Int, l::Int) = dot(view(x, 1:(lx-l)), view(x, (1+l):lx))
+_autodot(x::AbstractVector{<:NTuple{N,<:Real}}, lx::Int, l::Int) where N = dot(view(x, 1:(lx-l)), view(x, (1+l):lx))

test/runtests.jl

@@ -30,7 +30,7 @@ using Test
     @test m ≈ m2
 
     # emsd is the average of imsds
-    using StatsBase: mean
+    using Statistics: mean
     @test emsd(x) ≈ vec(mean(m2, dims=2))
 
     x = [collect(1:10), collect(1:10)]