diff --git a/Project.toml b/Project.toml index 548ed44..26369f4 100644 --- a/Project.toml +++ b/Project.toml @@ -1,20 +1,19 @@ name = "MeanSquaredDisplacement" uuid = "13c93d70-909c-440c-af92-39d48ffa2ba2" authors = ["Riccardo Foffi "] -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] diff --git a/src/MeanSquaredDisplacement.jl b/src/MeanSquaredDisplacement.jl index bcd7ab9..5b65e82 100644 --- a/src/MeanSquaredDisplacement.jl +++ b/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 diff --git a/src/acf.jl b/src/acf.jl index e2517d5..e55a85e 100644 --- a/src/acf.jl +++ b/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)) - - - diff --git a/src/autodot.jl b/src/autodot.jl new file mode 100644 index 0000000..c8e3e94 --- /dev/null +++ b/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)) diff --git a/test/runtests.jl b/test/runtests.jl index f6b2bfa..162eede 100644 --- a/test/runtests.jl +++ b/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)]