From 0f3924467b6ef7ac6623e8b1a0da166057f0f6c0 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Thu, 13 Oct 2022 16:19:34 +0200 Subject: [PATCH 01/27] Multiscale analysis --- docs/src/multiscale.md | 15 ++ docs/toc.jl | 3 +- src/Entropies.jl | 1 + src/multiscale.jl | 213 ++++++++++++++++++++++++++ test/multiscale/multiscale_entropy.jl | 7 + test/runtests.jl | 6 +- 6 files changed, 243 insertions(+), 2 deletions(-) create mode 100644 docs/src/multiscale.md create mode 100644 src/multiscale.jl create mode 100644 test/multiscale/multiscale_entropy.jl diff --git a/docs/src/multiscale.md b/docs/src/multiscale.md new file mode 100644 index 000000000..8c0db7f01 --- /dev/null +++ b/docs/src/multiscale.md @@ -0,0 +1,15 @@ +# Multiscale + +## API + +```@docs +multiscale +downsample +``` + +## Algorithms + +```@docs +Regular +Composite +``` diff --git a/docs/toc.jl b/docs/toc.jl index 89adeb08f..55e9b10ee 100644 --- a/docs/toc.jl +++ b/docs/toc.jl @@ -3,6 +3,7 @@ ENTROPIES_PAGES = [ "Probabilities" => "probabilities.md", "Entropies" => "entropies.md", "Complexity measures" => "complexity_measures.md", + "Multiscale" => "multiscale.md", "Entropies.jl examples" => "examples.md", "Utility methods" => "utils.md", -] \ No newline at end of file +] diff --git a/src/Entropies.jl b/src/Entropies.jl index d4f286776..e7bdb3d89 100644 --- a/src/Entropies.jl +++ b/src/Entropies.jl @@ -22,6 +22,7 @@ include("symbolization/symbolize.jl") include("probabilities_estimators/probabilities_estimators.jl") include("entropies/entropies.jl") include("complexity_measures/complexity_measures.jl") +include("multiscale.jl") include("deprecations.jl") diff --git a/src/multiscale.jl b/src/multiscale.jl new file mode 100644 index 000000000..61ef6dfad --- /dev/null +++ b/src/multiscale.jl @@ -0,0 +1,213 @@ +# This file contains an API for multiscale (coarse-grained/downsampled) computations of +# entropy and various complexity measures on time series. +# Can be used both to compute actual entropies (i.e. diversity entropy), or +# complexity measures (sample entropy, approximate entropy). + +using Statistics +export multiscale +export downsample +export Regular +export Composite +export MultiScaleAlgorithm + +""" + MultiScaleAlgorithm + +An abstract type for multiscale methods. +""" +abstract type MultiScaleAlgorithm end + +""" + Regular <: MultiScaleAlgorithm + Regular(; f::Function = Statistics.mean) + +The original multi-scale algorithm for multiscale entropy analysis (Costa et al., +2022)[^Costa2002], which yields a single downsampled time series per scale `s`. + +## Description + +Given a scalar-valued input time series `x`, the `Regular` multiscale algorithm downsamples +and coarse-grains `x` by splitting it into non-overlapping windows of length `s`, and +then constructing a new downsampled time series `D_t(s, f)` by applying the function `f` +to each of the resulting length-`s` windows. + +The downsampled time series `D_t(s)` with `t ∈ [1, 2, …, L]`, where `L = floor(N / s)`, +is given by: + +```math +\\{ D_t(s, f) \\}_{t = 1}^{L} = \\left\\{ f \\left( \\bf x_t \\right) \\right\\}_{t = 1}^{L} = +\\left\\{ + {f\\left( (x_i)_{i = (t - 1)s + 1}^{ts} \\right)} +\\right\\}_{t = 1}^{L} +``` + +where `f` is some summary statistic applied to the length-`ts-((t - 1)s + 1)` tuples `xₖ`. +Different choices of `f` have yield different multiscale methods appearing in the +literature. For example: + +- `f == Statistics.mean` yields the original first-moment multiscale entropy (Costa + et al., 2002)[^Costa2002]. +- `f == Statistics.var` yields the "generalized multiscale entropy" (Costa & Goldberger, + 2015)[^Costa2015], which uses the second-moment (variance) instead of the mean. + +[^Costa2002]: Costa, M., Goldberger, A. L., & Peng, C. K. (2002). Multiscale entropy + analysis of complex physiologic time series. Physical review letters, 89(6), 068102. +[^Costa2015]: Costa, M. D., & Goldberger, A. L. (2015). Generalized multiscale entropy + analysis: Application to quantifying the complex volatility of human heartbeat time + series. Entropy, 17(3), 1197-1203. +""" +Base.@kwdef struct Regular <: MultiScaleAlgorithm + f::Function = Statistics.mean +end + +""" + Composite <: MultiScaleAlgorithm + Composite(; f::Function = Statistics.mean) + +Composite multi-scale algorithm for multiscale entropy analysis (Wu et al., +2013)[^Wu2013], used, with [`multiscale`](@ref) to compute, for example, composite +multiscale entropy (CMSE). + +## Description + +Given a scalar-valued input time series `x`, the composite multiscale algorithm, +like [`Regular`](@ref), downsamples and coarse-grains `x` by splitting it into +non-overlapping windows of length `s`, and then constructing downsampled time series by +applying the function `f` to each of the resulting length-`s` windows. + +However, Wu et al. (2013)[^Wu2013] realized that for each scale `s`, there are actually `s` +different ways of selecting windows, depending on where indexing starts/ends. +These `s` different downsampled time series `D_t(s, f)` at each scale `s` are +constructed as follows: + +```math +\\{ D_{k}(s) \\} = \\{ D_{t, k}(s) \\}_{t = 1}^{L}, = \\{ f \\left( \\bf x_{t, k} \\right) \\} = +\\left\\{ + {f\\left( (x_i)_{i = (t - 1)s + k}^{ts + k - 1} \\right)} +\\right\\}_{t = 1}^{L}, +``` + +where `L = floor((N - s + 1) / s)` and `1 ≤ k ≤ s`, such that ``D_{i, k}(s)`` is the `i`-th +element of the `k`-th downsampled time series at scale `s`. + +Finally, compute ``\\dfrac{1}{s} \\sum_{k = 1}^s g(D_{k}(s))``, where `g` is some summary +function, for example [`entropy`](@ref). + +!!! note "Relation to [`Regular`](@ref)" + The downsampled time series ``D_{t, 1}(s)`` constructed using the composite + multiscale method is equivalent to the downsampled time series `D_{t}(s)` constructed + using the [`Regular`](@ref) method, for which `k == 1` is fixed, such that only + a single time series is returned. + +[^Wu2013]: Wu, S. D., Wu, C. W., Lin, S. G., Wang, C. C., & Lee, K. Y. (2013). Time series + analysis using composite multiscale entropy. Entropy, 15(3), 1069-1084. +""" +Base.@kwdef struct Composite <: MultiScaleAlgorithm + f::Function = Statistics.mean +end + +""" + downsample(algorithm::MultiScaleAlgorithm, x::AbstractVector{T}, s::Int, args...; + kwargs...) + +Downsample and coarse-grain `x` to scale `s` according to the given multiscale `algorithm`. + +Positional arguments `args` and keyword arguments `kwargs` are propagated to relevant +functions in `algorithm`, if applicable. + +`The return type depends on `algorithm`. For example: + +- [`Regular`](@ref) yields a single `Vector` per scale. +- [`Composite`](@ref) yields a `Vector{Vector}` per scale. +""" +function downsample(method::MultiScaleAlgorithm, x::AbstractVector{T}, s::Int, args...; kwargs...) where T + f = method.f + + if s == 1 + return x + else + N = length(x) + L = floor(Int, N / s) + ys = zeros(T, L) + + for t = 1:L + inds = ((t - 1)*s + 1):(t * s) + ys[t] = @views f(x[inds], args...; kwargs...) + end + return ys + end +end + +function downsample(method::Composite, x::AbstractVector{T}, s::Int, args...; kwargs...) where T + f = method.f + + if s == 1 + return [x] + else + N = length(x) + # note: there must be a typo or error in Wu et al. (2013), because if we use + # floor(N / s), as indicated in their paper, we'll get out of bounds errors. + # The number of samples needs to take into consideratio how many unique ways + # there are of selecting non-overlapping windows of length `s`. Hence, + # we use floor((N - s + 1) / s) instead. + L = floor(Int, (N - s + 1) / s) + ys = [zeros(T, L) for i = 1:s] + for k = 1:s + for t = 1:L + inds = ((t - 1)*s + k):(t * s + k - 1) + ys[k][t] = @views f(x[inds], args...; kwargs...) + end + end + return ys + end +end + +""" + multiscale(e::Entropy, alg::MultiScaleAlgorithm, x::AbstractVector, est::ProbabilitiesEstimator; + scalemax::Int = 8, normalize = false) + +Compute the multi-scale entropy (Costa et al., 2002)[^Costa2002] of type `e` of +timeseries `x` using the given probabilities estimator `est` + +Utilizes [`downsample`](@ref) to compute the entropy of coarse-grained, downsampled +versions of `x` for scale factors `1:maxscale`. The length of the most severely +downsampled version of `x` is `N ÷ maxscale`, while for scale factor `1`, the original +time series is considered. + +If `normalize == true`, then compute normalized entropy (if that is possible for this +particular combination of entropy type and probability estimator). + +[^Costa2002]: Costa, M., Goldberger, A. L., & Peng, C. K. (2002). Multiscale entropy + analysis of complex physiologic time series. Physical review letters, 89(6), 068102. +""" +function multiscale end + +function multiscale(e::Entropy, alg::Regular, x::AbstractVector, est::ProbabilitiesEstimator; + maxscale::Int = 8, normalize = false) + + downscaled_timeseries = [downsample(alg, x, s) for s in 1:maxscale] + hs = zeros(Float64, maxscale) + if normalize + hs = entropy_normalized.(Ref(e), downscaled_timeseries, Ref(est)) + else + hs .= entropy.(Ref(e), downscaled_timeseries, Ref(est)) + end + + return hs +end + +function multiscale(e::Entropy, alg::Composite, x::AbstractVector, est::ProbabilitiesEstimator; + maxscale::Int = 10, normalize = false) + + downscaled_timeseries = [downsample(alg, x, s) for s in 1:maxscale] + hs = zeros(Float64, maxscale) + for s in 1:maxscale + if normalize + hs[s] = mean(entropy_normalized.(Ref(e), downscaled_timeseries[s], Ref(est))) + else + hs[s] = mean(entropy.(Ref(e), downscaled_timeseries[s], Ref(est))) + end + end + + return hs +end diff --git a/test/multiscale/multiscale_entropy.jl b/test/multiscale/multiscale_entropy.jl new file mode 100644 index 000000000..4693d44e2 --- /dev/null +++ b/test/multiscale/multiscale_entropy.jl @@ -0,0 +1,7 @@ +est = Dispersion() +x = rand(100) +e = Shannon() + +# TODO: test actual outcomes too. Make some analytical examples. +@test multiscale(e, Regular(), x, est) isa Vector{T} where T <: Real +@test multiscale(e, Composite(), x, est) isa Vector{T} where T <: Real diff --git a/test/runtests.jl b/test/runtests.jl index 9aac82312..9ff89b589 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -24,8 +24,12 @@ testfile(file, testname=defaultname(file)) = @testset "$testname" begin; include testfile("entropies/tsallis.jl") testfile("entropies/curado.jl") testfile("entropies/stretched_exponential.jl") - + testfile("entropies/nearest_neighbors_direct.jl") + + # Multiscale analysis + testfile("multiscale/multiscale_entropy.jl") + # Various testfile("complexity_measures/complexity_measures.jl") testfile("utils/utils.jl") From 944727c9e09fe8cbd07d7fd94a71ab51ab63c865 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Tue, 18 Oct 2022 14:12:47 +0200 Subject: [PATCH 02/27] Complexity API, and rework reverse dispersion --- docs/src/complexity_measures.md | 9 +- src/Entropies.jl | 2 +- src/complexity.jl | 36 +++++++ .../reverse_dispersion_entropy.jl | 98 +++++++++++-------- src/symbolization/GaussianSymbolization.jl | 2 + src/symbolization/symbolize.jl | 7 ++ .../complexity_measures/reverse_dispersion.jl | 8 +- 7 files changed, 112 insertions(+), 50 deletions(-) create mode 100644 src/complexity.jl diff --git a/docs/src/complexity_measures.md b/docs/src/complexity_measures.md index c063d28c1..4aacfc7e4 100644 --- a/docs/src/complexity_measures.md +++ b/docs/src/complexity_measures.md @@ -1,8 +1,13 @@ -# [Complexity measures](@id complexity_measures) +# [Complexity API](@id complexity_measures) + +```@docs +complexity +complexity_normalized! +``` ## Reverse dispersion entropy ```@docs -reverse_dispersion +ReverseDispersion distance_to_whitenoise ``` diff --git a/src/Entropies.jl b/src/Entropies.jl index d4f286776..b4d0a8b4b 100644 --- a/src/Entropies.jl +++ b/src/Entropies.jl @@ -21,7 +21,7 @@ include("entropy.jl") include("symbolization/symbolize.jl") include("probabilities_estimators/probabilities_estimators.jl") include("entropies/entropies.jl") -include("complexity_measures/complexity_measures.jl") +include("complexity.jl") include("deprecations.jl") diff --git a/src/complexity.jl b/src/complexity.jl new file mode 100644 index 000000000..33d687622 --- /dev/null +++ b/src/complexity.jl @@ -0,0 +1,36 @@ +export ComplexityMeasure +export complexity +export complexity_normalized + +""" + ComplexityMeasure + +Abstract type for (entropy-like) complexity measures. +""" +abstract type ComplexityMeasure end + +""" + complexity(c::ComplexityMeasure, x) + +Estimate the complexity measure `c` for input data `x`, where `c` can be any of the +following measures: + +- [`ReverseDispersion`](@ref). + +""" +function complexity end + +""" + complexity_normalized(c::ComplexityMeasure, x) → m ∈ [a, b] + +Estimate the normalized complexity measure `c` for input data `x`, where `c` can +can be any of the following measures: + +- [`ReverseDispersion`](@ref). + +The potential range `[a, b]` of the output value depends on `c`. See the documentation +strings for the individual measures to get the normalized ranges. +""" +function complexity_normalized end + +include("complexity_measures/complexity_measures.jl") diff --git a/src/complexity_measures/reverse_dispersion_entropy.jl b/src/complexity_measures/reverse_dispersion_entropy.jl index 7b66b39ae..9768adcc9 100644 --- a/src/complexity_measures/reverse_dispersion_entropy.jl +++ b/src/complexity_measures/reverse_dispersion_entropy.jl @@ -1,42 +1,13 @@ -export reverse_dispersion +export ReverseDispersion export distance_to_whitenoise -# Note: this is not an entropy estimator, so we don't use the entropy_xxx_norm interface -# for normalization, even though we rely on `alphabet_length`. """ - distance_to_whitenoise(p::Probabilities, estimator::Dispersion; normalize = false) + ReverseDispersion <: ComplexityMeasure + ReverseDispersion(; m = 2, τ = 1, check_unique = true, + symbolization::SymbolizationScheme = GaussianSymbolization(c = 5) + ) -Compute the distance of the probability distribution `p` from a uniform distribution, -given the parameters of `estimator` (which must be known beforehand). - -If `normalize == true`, then normalize the value to the interval `[0, 1]` by using the -parameters of `estimator`. - -Used to compute reverse dispersion entropy([`reverse_dispersion`](@ref); -Li et al., 2019[^Li2019]). - -[^Li2019]: Li, Y., Gao, X., & Wang, L. (2019). Reverse dispersion entropy: a new - complexity measure for sensor signal. Sensors, 19(23), 5203. -""" -function distance_to_whitenoise(p::Probabilities, est::Dispersion; normalize = false) - # We can safely skip non-occurring symbols, because they don't contribute - # to the sum in eq. 3 in Li et al. (2019) - Hrde = sum(abs2, p) - (1 / alphabet_length(est)) - - if normalize - return Hrde / (1 - (1 / alphabet_length(est))) - else - return Hrde - end -end - -# Note again: this is a *complexity measure*, not an entropy estimator, so we don't use -# the entropy_xxx_norm interface for normalization, even though we rely on `alphabet_length`. -""" - reverse_dispersion(x::AbstractVector{T}, est::Dispersion = Dispersion(); - normalize = true) where T <: Real - -Compute the reverse dispersion entropy complexity measure (Li et al., 2019)[^Li2019]. +Estimator for the reverse dispersion entropy complexity measure (Li et al., 2019)[^Li2019]. ## Description @@ -54,23 +25,64 @@ embedding dimension `m` and embedding delay `τ`. Recommended parameter values[^Li2018] are `m ∈ [2, 3]`, `τ = 1` for the embedding, and `c ∈ [3, 4, …, 8]` categories for the Gaussian mapping. -If `normalize == true`, then the reverse dispersion entropy is normalized to `[0, 1]`. +If normalizing, then the reverse dispersion entropy is normalized to `[0, 1]`. The minimum value of ``H_{rde}`` is zero and occurs precisely when the dispersion pattern distribution is flat, which occurs when all ``p_i``s are equal to ``1/c^m``. Because ``H_{rde} \\geq 0``, ``H_{rde}`` can therefore be said to be a measure of how far the dispersion pattern probability distribution is from white noise. +## Data requirements + +Like for [`Dispersion`](@ref), the input must have more than one unique element for the +default Gaussian mapping symbolization to be well-defined. Li et al. (2018) recommends +that `x` has at least 1000 data points. + +If `check_unique == true` (default), then it is checked that the input has +more than one unique value. If `check_unique == false` and the input only has one +unique element, then a `InexactError` is thrown when trying to compute probabilities. + [^Li2019]: Li, Y., Gao, X., & Wang, L. (2019). Reverse dispersion entropy: a new complexity measure for sensor signal. Sensors, 19(23), 5203. """ -function reverse_dispersion(x::AbstractVector{T}, est::Dispersion = Dispersion(); - normalize = true) where T <: Real +struct ReverseDispersion{S <: SymbolizationScheme} <: ComplexityMeasure + symbolization::S = GaussianSymbolization(c = 5) + m::Int = 2 + τ::Int = 1 + check_unique::Bool = false +end + +alphabet_length(est::ReverseDispersion)::Int = alphabet_length(est.symbolization) ^ est.m + +""" + distance_to_whitenoise(p::Probabilities, estimator::ReverseDispersion; normalize = false) - p = probabilities(x, est) +Compute the distance of the probability distribution `p` from a uniform distribution, +given the parameters of `estimator` (which must be known beforehand). + +If `normalize == true`, then normalize the value to the interval `[0, 1]` by using the +parameters of `estimator`. + +Used to compute reverse dispersion entropy([`reverse_dispersion`](@ref); +Li et al., 2019[^Li2019]). + +[^Li2019]: Li, Y., Gao, X., & Wang, L. (2019). Reverse dispersion entropy: a new + complexity measure for sensor signal. Sensors, 19(23), 5203. +""" +function distance_to_whitenoise(p::Probabilities, est::ReverseDispersion; normalize = false) + # We can safely skip non-occurring symbols, because they don't contribute + # to the sum in eq. 3 in Li et al. (2019) + return sum(abs2, p) - (1 / alphabet_length(est)) +end + +function complexity(c::ReverseDispersion, x) + p = probabilities(x, c) + # The following step combines distance information with the probabilities, yielding + # something which is no longer a set of probabilities. + Hrde = distance_to_whitenoise(p, c, normalize = c.normalize) + return Hrde +end - # The following step combines distance information with the probabilities, so - # from here on, it is not possible to use `renyi_entropy` or similar methods, because - # we're not dealing with probabilities anymore. - Hrde = distance_to_whitenoise(p, est, normalize = normalize) +function complexity_normalized(c::ReverseDispersion, x) + return complexity(c, x) / (1 - (1 / alphabet_length(c))) end diff --git a/src/symbolization/GaussianSymbolization.jl b/src/symbolization/GaussianSymbolization.jl index 6ceb4d11b..48063517c 100644 --- a/src/symbolization/GaussianSymbolization.jl +++ b/src/symbolization/GaussianSymbolization.jl @@ -53,6 +53,8 @@ Base.@kwdef struct GaussianSymbolization{I <: Integer} <: SymbolizationScheme c::I = 3 end +alphabet_length(symbolization::GaussianSymbolization) = s.c + g(xᵢ, μ, σ) = exp((-(xᵢ - μ)^2)/(2σ^2)) """ diff --git a/src/symbolization/symbolize.jl b/src/symbolization/symbolize.jl index 52bf93bd1..fe7af4612 100644 --- a/src/symbolization/symbolize.jl +++ b/src/symbolization/symbolize.jl @@ -5,6 +5,13 @@ An abstract type for symbolization schemes. """ abstract type SymbolizationScheme end +# The internal structure of different symbolization schemes may be different, so use +# `alphabet_length` to have a consistent way of getting the total number of possible states. +# Thus, the default behaviour is to throw an ArgumentError when computing some normalized +# quantity depending on `alphabet_length` of the symbolization scheme. +alphabet_length(s::S) where S <: SymbolizationScheme = + throw(ArgumentError("`alphabet_length` not defined for $S.")) + """ symbolize(x, scheme::SymbolizationScheme) → Vector{Int} symbolize!(s, x, scheme::SymbolizationScheme) → Vector{Int} diff --git a/test/complexity_measures/reverse_dispersion.jl b/test/complexity_measures/reverse_dispersion.jl index d9423c7c4..7112962fe 100644 --- a/test/complexity_measures/reverse_dispersion.jl +++ b/test/complexity_measures/reverse_dispersion.jl @@ -1,10 +1,10 @@ x = rand(100) -@test reverse_dispersion(x) isa Real -@test 0.0 <= reverse_dispersion(x, normalize = true) <= 1.0 +@test complexity(ReverseDispersion(), x) isa Real +@test 0.0 <= complexity_normalized(ReverseDispersion(), x) <= 1.0 @testset "Distance to whitenoise" begin m, n_classes = 2, 2 - est = Dispersion(m = m, symbolization = GaussianSymbolization(c = n_classes)) + est = ReverseDispersion(m = m, symbolization = GaussianSymbolization(c = n_classes)) # Reverse dispersion entropy is 0 when all probabilities are identical and equal # to 1/(n_classes^m). @@ -12,7 +12,7 @@ x = rand(100) Hrde_minimal = distance_to_whitenoise(flat_dist, est, normalize = false) @test round(Hrde_minimal, digits = 7) ≈ 0.0 - # Reverse dispersion entropy is maximal when there is only one non-zero dispersal + # Reverse dispersion entropy is maximal when there is only one non-zero dispersal # pattern. Then reverse dispersion entropy is # 1 - 1/(n_classes^m). When normalizing to this value, the RDE should be 1.0. m, n_classes = 2, 2 From c3a70080c185c581581274a0476a6097f7a2850e Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Tue, 18 Oct 2022 14:15:18 +0200 Subject: [PATCH 03/27] Remember to use Base.@kwdef --- src/complexity_measures/reverse_dispersion_entropy.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/complexity_measures/reverse_dispersion_entropy.jl b/src/complexity_measures/reverse_dispersion_entropy.jl index 9768adcc9..fd328102b 100644 --- a/src/complexity_measures/reverse_dispersion_entropy.jl +++ b/src/complexity_measures/reverse_dispersion_entropy.jl @@ -45,7 +45,7 @@ unique element, then a `InexactError` is thrown when trying to compute probabili [^Li2019]: Li, Y., Gao, X., & Wang, L. (2019). Reverse dispersion entropy: a new complexity measure for sensor signal. Sensors, 19(23), 5203. """ -struct ReverseDispersion{S <: SymbolizationScheme} <: ComplexityMeasure +Base.@kwdef struct ReverseDispersion{S <: SymbolizationScheme} <: ComplexityMeasure symbolization::S = GaussianSymbolization(c = 5) m::Int = 2 τ::Int = 1 From 949ad52bc35a89fd8201a7bfe24f3eacee99fd2a Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Tue, 18 Oct 2022 14:24:59 +0200 Subject: [PATCH 04/27] Fix normalization --- .../reverse_dispersion_entropy.jl | 20 ++++++++++++------- src/symbolization/GaussianSymbolization.jl | 2 +- 2 files changed, 14 insertions(+), 8 deletions(-) diff --git a/src/complexity_measures/reverse_dispersion_entropy.jl b/src/complexity_measures/reverse_dispersion_entropy.jl index fd328102b..09ca1c03a 100644 --- a/src/complexity_measures/reverse_dispersion_entropy.jl +++ b/src/complexity_measures/reverse_dispersion_entropy.jl @@ -72,17 +72,23 @@ Li et al., 2019[^Li2019]). function distance_to_whitenoise(p::Probabilities, est::ReverseDispersion; normalize = false) # We can safely skip non-occurring symbols, because they don't contribute # to the sum in eq. 3 in Li et al. (2019) - return sum(abs2, p) - (1 / alphabet_length(est)) + Hrde = sum(abs2, p) - (1 / alphabet_length(est)) + + if normalize + return Hrde / (1 - (1 / alphabet_length(est))) + else + return Hrde + end end function complexity(c::ReverseDispersion, x) - p = probabilities(x, c) - # The following step combines distance information with the probabilities, yielding - # something which is no longer a set of probabilities. - Hrde = distance_to_whitenoise(p, c, normalize = c.normalize) - return Hrde + (; symbolization, m, τ, check_unique) = c + p = probabilities(x, Dispersion(; symbolization, m, τ, check_unique)) + return distance_to_whitenoise(p, c, normalize = false) end function complexity_normalized(c::ReverseDispersion, x) - return complexity(c, x) / (1 - (1 / alphabet_length(c))) + (; symbolization, m, τ, check_unique) = c + p = probabilities(x, Dispersion(; symbolization, m, τ, check_unique)) + return distance_to_whitenoise(p, c, normalize = true) end diff --git a/src/symbolization/GaussianSymbolization.jl b/src/symbolization/GaussianSymbolization.jl index 48063517c..afe45e732 100644 --- a/src/symbolization/GaussianSymbolization.jl +++ b/src/symbolization/GaussianSymbolization.jl @@ -53,7 +53,7 @@ Base.@kwdef struct GaussianSymbolization{I <: Integer} <: SymbolizationScheme c::I = 3 end -alphabet_length(symbolization::GaussianSymbolization) = s.c +alphabet_length(symbolization::GaussianSymbolization) = symbolization.c g(xᵢ, μ, σ) = exp((-(xᵢ - μ)^2)/(2σ^2)) From 0c7958d7d6a5958d51291a5e77036328051f095a Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Tue, 18 Oct 2022 14:33:15 +0200 Subject: [PATCH 05/27] Fix docs --- src/complexity_measures/reverse_dispersion_entropy.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/complexity_measures/reverse_dispersion_entropy.jl b/src/complexity_measures/reverse_dispersion_entropy.jl index 09ca1c03a..d98eb9900 100644 --- a/src/complexity_measures/reverse_dispersion_entropy.jl +++ b/src/complexity_measures/reverse_dispersion_entropy.jl @@ -63,7 +63,7 @@ given the parameters of `estimator` (which must be known beforehand). If `normalize == true`, then normalize the value to the interval `[0, 1]` by using the parameters of `estimator`. -Used to compute reverse dispersion entropy([`reverse_dispersion`](@ref); +Used to compute reverse dispersion entropy([`ReverseDispersion`](@ref); Li et al., 2019[^Li2019]). [^Li2019]: Li, Y., Gao, X., & Wang, L. (2019). Reverse dispersion entropy: a new From 9b9b073838e880acce3fdcb2ef3dc7edc82ed4ba Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Tue, 18 Oct 2022 14:50:31 +0200 Subject: [PATCH 06/27] Typo --- docs/src/complexity_measures.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/complexity_measures.md b/docs/src/complexity_measures.md index 4aacfc7e4..e983b58a4 100644 --- a/docs/src/complexity_measures.md +++ b/docs/src/complexity_measures.md @@ -2,7 +2,7 @@ ```@docs complexity -complexity_normalized! +complexity_normalized ``` ## Reverse dispersion entropy From 39affdf72436ef36bd0544799724bce8ec4f0f9a Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Tue, 18 Oct 2022 15:35:30 +0200 Subject: [PATCH 07/27] Multiscale complexity --- src/Entropies.jl | 1 + src/multiscale.jl | 33 +++++++++++++++++++++++++++++++++ 2 files changed, 34 insertions(+) diff --git a/src/Entropies.jl b/src/Entropies.jl index b4d0a8b4b..6b613cbdc 100644 --- a/src/Entropies.jl +++ b/src/Entropies.jl @@ -22,6 +22,7 @@ include("symbolization/symbolize.jl") include("probabilities_estimators/probabilities_estimators.jl") include("entropies/entropies.jl") include("complexity.jl") +include("multiscale.jl") include("deprecations.jl") diff --git a/src/multiscale.jl b/src/multiscale.jl index 61ef6dfad..7e8c8a3ef 100644 --- a/src/multiscale.jl +++ b/src/multiscale.jl @@ -211,3 +211,36 @@ function multiscale(e::Entropy, alg::Composite, x::AbstractVector, est::Probabil return hs end + + + +function multiscale(e::ComplexityMeasure, alg::Regular, x::AbstractVector; + maxscale::Int = 8, normalize = false) + + downscaled_timeseries = [downsample(alg, x, s) for s in 1:maxscale] + complexities = zeros(Float64, maxscale) + if normalize + complexities = complexity_normalized.(Ref(e), downscaled_timeseries) + else + complexities .= complexity.(Ref(e), downscaled_timeseries) + end + + return complexities +end + +function multiscale(e::ComplexityMeasure, alg::Composite, x::AbstractVector; + maxscale::Int = 10, normalize = false) + + downscaled_timeseries = [downsample(alg, x, s) for s in 1:maxscale] + complexities = zeros(Float64, maxscale) + for s in 1:maxscale + if normalize + complexities[s] = + mean(complexity_normalized.(Ref(e), downscaled_timeseries[s], Ref(est))) + else + complexities[s] = mean(entropy.(Ref(e), downscaled_timeseries[s], Ref(est))) + end + end + + return complexities +end From b73383eedc2fb7a9f0bf3c3056d187780a5b5fe2 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Tue, 18 Oct 2022 22:57:30 +0200 Subject: [PATCH 08/27] Update docs --- src/multiscale.jl | 32 ++++++++++++++++++++++---------- 1 file changed, 22 insertions(+), 10 deletions(-) diff --git a/src/multiscale.jl b/src/multiscale.jl index 7e8c8a3ef..7c7dfbe62 100644 --- a/src/multiscale.jl +++ b/src/multiscale.jl @@ -163,19 +163,31 @@ function downsample(method::Composite, x::AbstractVector{T}, s::Int, args...; kw end """ - multiscale(e::Entropy, alg::MultiScaleAlgorithm, x::AbstractVector, est::ProbabilitiesEstimator; - scalemax::Int = 8, normalize = false) + multiscale(e::Entropy, alg::MultiScaleAlgorithm, x::AbstractVector, + est::ProbabilitiesEstimator; + scalemax::Int = 8, + normalize = false) -Compute the multi-scale entropy (Costa et al., 2002)[^Costa2002] of type `e` of -timeseries `x` using the given probabilities estimator `est` +Compute the multi-scale entropy (Costa et al., 2002)[^Costa2002] of type `e` for +timeseries `x`, using downsampling algorithm `alg`, and probabilities estimator `est`. -Utilizes [`downsample`](@ref) to compute the entropy of coarse-grained, downsampled -versions of `x` for scale factors `1:maxscale`. The length of the most severely -downsampled version of `x` is `N ÷ maxscale`, while for scale factor `1`, the original -time series is considered. +If `normalize == true`, then compute normalized entropy. + + multiscale(c::ComplexityMeasure, alg::MultiScaleAlgorithm, x::AbstractVector; + scalemax::Int = 8, + normalize = false) + +Compute the multi-scale complexity measure `c` for the timeseries `x`, using +using downsampling algorithm `alg`. + +If `normalize == true`, then compute the normalized version of the complexity measure. -If `normalize == true`, then compute normalized entropy (if that is possible for this -particular combination of entropy type and probability estimator). +## Description + +Utilizes [`downsample`](@ref) to compute the entropy of coarse-grained, downsampled +versions of `x` for scale factors `1:maxscale`. If `N = length(x)`, then the length of the +most severely downsampled version of `x` is `N ÷ scalemax`, while for scale factor `1`, +the original time series is considered. [^Costa2002]: Costa, M., Goldberger, A. L., & Peng, C. K. (2002). Multiscale entropy analysis of complex physiologic time series. Physical review letters, 89(6), 068102. From bd4e4dcdeb5267b55d458578bc107173991d8e13 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Sat, 22 Oct 2022 22:30:55 +0200 Subject: [PATCH 09/27] remove unnecessary spaces --- src/multiscale.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/multiscale.jl b/src/multiscale.jl index 7c7dfbe62..647cf0ceb 100644 --- a/src/multiscale.jl +++ b/src/multiscale.jl @@ -224,8 +224,6 @@ function multiscale(e::Entropy, alg::Composite, x::AbstractVector, est::Probabil return hs end - - function multiscale(e::ComplexityMeasure, alg::Regular, x::AbstractVector; maxscale::Int = 8, normalize = false) From 21483ecb298582007d6e8cd45f38b95f9b5a74fe Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Sun, 23 Oct 2022 19:06:41 +0200 Subject: [PATCH 10/27] Docs --- docs/src/complexity_measures.md | 2 +- docs/src/entropies.md | 2 +- docs/src/multiscale.md | 40 ++++- src/multiscale.jl | 249 ++++---------------------- src/multiscale/composite.jl | 107 +++++++++++ src/multiscale/regular.jl | 93 ++++++++++ test/multiscale/multiscale_entropy.jl | 34 +++- 7 files changed, 304 insertions(+), 223 deletions(-) create mode 100644 src/multiscale/composite.jl create mode 100644 src/multiscale/regular.jl diff --git a/docs/src/complexity_measures.md b/docs/src/complexity_measures.md index 3e6b341d4..5262e75d8 100644 --- a/docs/src/complexity_measures.md +++ b/docs/src/complexity_measures.md @@ -1,4 +1,4 @@ -# [Complexity measures API](@id complexity_measures) +## [Complexity measures API](@id complexity_measures) ```@docs complexity diff --git a/docs/src/entropies.md b/docs/src/entropies.md index c85b3bc49..b74e3f449 100644 --- a/docs/src/entropies.md +++ b/docs/src/entropies.md @@ -1,4 +1,4 @@ -## Entropy API +## [Entropy API](@id entropies_list) ```@docs entropy diff --git a/docs/src/multiscale.md b/docs/src/multiscale.md index 8c0db7f01..cab5f78d1 100644 --- a/docs/src/multiscale.md +++ b/docs/src/multiscale.md @@ -1,15 +1,47 @@ -# Multiscale - -## API +## Multiscale API ```@docs multiscale downsample ``` -## Algorithms +## Available literature methods + +A non-exhaustive list of literature methods, and the syntax to compute them, are listed +below. Please open an issue or make a pull-request to +[Entropies.jl](https://github.com/JuliaDynamics/Entropies.jl) if you find a literature +method missing from this list, or if you publish a paper based on some new multiscale +combination. + +### Multiscale complexity measures + +| Method | Syntax | Reference | +| ------------- | ------------- | ------------- | +| Multiscale sample entropy (first moment) | `multiscale(SampleEntropy(), Regular(f = mean), x)` |Costa et al. (2002)[^Costa2002] | +| Generalized multiscale sample entropy (second moment)| `multiscale(SampleEntropy(), Regular(f = std), x)` | Costa et al. (2015)[^Costa2015] | + +### Multiscale entropy + +| Method | Syntax | Reference | +| ------------- | ------------- | ------------- | +| Refined composite multiscale dispersion entropy | `multiscale(Dispersion(), Composite(), x, normalized = true)` | Azami et al. (2017)[^Azami2017] | + +## [Algorithms](@id multiscale_algorithms) ```@docs +MultiScaleAlgorithm Regular Composite ``` + +[^Costa2002]: Costa, M., Goldberger, A. L., & Peng, C. K. (2002). Multiscale entropy + analysis of complex physiologic time series. Physical review letters, 89(6), 068102. +[^Costa2015]: Costa, M. D., & Goldberger, A. L. (2015). Generalized multiscale entropy + analysis: Application to quantifying the complex volatility of human heartbeat time + series. Entropy, 17(3), 1197-1203. +[^Azami2017]: Azami, H., Rostaghi, M., Abásolo, D., & Escudero, J. (2017). Refined + composite multiscale dispersion entropy and its application to biomedical signals. + IEEE Transactions on Biomedical Engineering, 64(12), 2872-2879. +[^Morabito2012]: Morabito, F. C., Labate, D., Foresta, F. L., Bramanti, A., Morabito, G., + & Palamara, I. (2012). Multivariate multi-scale permutation entropy for complexity + analysis of Alzheimer’s disease EEG. Entropy, 14(7), 1186-1202. diff --git a/src/multiscale.jl b/src/multiscale.jl index 647cf0ceb..12cb9ae29 100644 --- a/src/multiscale.jl +++ b/src/multiscale.jl @@ -6,251 +6,72 @@ using Statistics export multiscale export downsample -export Regular -export Composite export MultiScaleAlgorithm """ MultiScaleAlgorithm -An abstract type for multiscale methods. +The supertype for all multiscale algorithms. """ abstract type MultiScaleAlgorithm end """ - Regular <: MultiScaleAlgorithm - Regular(; f::Function = Statistics.mean) - -The original multi-scale algorithm for multiscale entropy analysis (Costa et al., -2022)[^Costa2002], which yields a single downsampled time series per scale `s`. - -## Description - -Given a scalar-valued input time series `x`, the `Regular` multiscale algorithm downsamples -and coarse-grains `x` by splitting it into non-overlapping windows of length `s`, and -then constructing a new downsampled time series `D_t(s, f)` by applying the function `f` -to each of the resulting length-`s` windows. - -The downsampled time series `D_t(s)` with `t ∈ [1, 2, …, L]`, where `L = floor(N / s)`, -is given by: - -```math -\\{ D_t(s, f) \\}_{t = 1}^{L} = \\left\\{ f \\left( \\bf x_t \\right) \\right\\}_{t = 1}^{L} = -\\left\\{ - {f\\left( (x_i)_{i = (t - 1)s + 1}^{ts} \\right)} -\\right\\}_{t = 1}^{L} -``` - -where `f` is some summary statistic applied to the length-`ts-((t - 1)s + 1)` tuples `xₖ`. -Different choices of `f` have yield different multiscale methods appearing in the -literature. For example: - -- `f == Statistics.mean` yields the original first-moment multiscale entropy (Costa - et al., 2002)[^Costa2002]. -- `f == Statistics.var` yields the "generalized multiscale entropy" (Costa & Goldberger, - 2015)[^Costa2015], which uses the second-moment (variance) instead of the mean. - -[^Costa2002]: Costa, M., Goldberger, A. L., & Peng, C. K. (2002). Multiscale entropy - analysis of complex physiologic time series. Physical review letters, 89(6), 068102. -[^Costa2015]: Costa, M. D., & Goldberger, A. L. (2015). Generalized multiscale entropy - analysis: Application to quantifying the complex volatility of human heartbeat time - series. Entropy, 17(3), 1197-1203. -""" -Base.@kwdef struct Regular <: MultiScaleAlgorithm - f::Function = Statistics.mean -end - -""" - Composite <: MultiScaleAlgorithm - Composite(; f::Function = Statistics.mean) - -Composite multi-scale algorithm for multiscale entropy analysis (Wu et al., -2013)[^Wu2013], used, with [`multiscale`](@ref) to compute, for example, composite -multiscale entropy (CMSE). - -## Description - -Given a scalar-valued input time series `x`, the composite multiscale algorithm, -like [`Regular`](@ref), downsamples and coarse-grains `x` by splitting it into -non-overlapping windows of length `s`, and then constructing downsampled time series by -applying the function `f` to each of the resulting length-`s` windows. - -However, Wu et al. (2013)[^Wu2013] realized that for each scale `s`, there are actually `s` -different ways of selecting windows, depending on where indexing starts/ends. -These `s` different downsampled time series `D_t(s, f)` at each scale `s` are -constructed as follows: - -```math -\\{ D_{k}(s) \\} = \\{ D_{t, k}(s) \\}_{t = 1}^{L}, = \\{ f \\left( \\bf x_{t, k} \\right) \\} = -\\left\\{ - {f\\left( (x_i)_{i = (t - 1)s + k}^{ts + k - 1} \\right)} -\\right\\}_{t = 1}^{L}, -``` - -where `L = floor((N - s + 1) / s)` and `1 ≤ k ≤ s`, such that ``D_{i, k}(s)`` is the `i`-th -element of the `k`-th downsampled time series at scale `s`. - -Finally, compute ``\\dfrac{1}{s} \\sum_{k = 1}^s g(D_{k}(s))``, where `g` is some summary -function, for example [`entropy`](@ref). - -!!! note "Relation to [`Regular`](@ref)" - The downsampled time series ``D_{t, 1}(s)`` constructed using the composite - multiscale method is equivalent to the downsampled time series `D_{t}(s)` constructed - using the [`Regular`](@ref) method, for which `k == 1` is fixed, such that only - a single time series is returned. - -[^Wu2013]: Wu, S. D., Wu, C. W., Lin, S. G., Wang, C. C., & Lee, K. Y. (2013). Time series - analysis using composite multiscale entropy. Entropy, 15(3), 1069-1084. -""" -Base.@kwdef struct Composite <: MultiScaleAlgorithm - f::Function = Statistics.mean -end - -""" - downsample(algorithm::MultiScaleAlgorithm, x::AbstractVector{T}, s::Int, args...; - kwargs...) + downsample(algorithm::MultiScaleAlgorithm, x, s::Int) Downsample and coarse-grain `x` to scale `s` according to the given multiscale `algorithm`. Positional arguments `args` and keyword arguments `kwargs` are propagated to relevant functions in `algorithm`, if applicable. -`The return type depends on `algorithm`. For example: +The return type depends on `algorithm`. For example: - [`Regular`](@ref) yields a single `Vector` per scale. - [`Composite`](@ref) yields a `Vector{Vector}` per scale. """ -function downsample(method::MultiScaleAlgorithm, x::AbstractVector{T}, s::Int, args...; kwargs...) where T - f = method.f - - if s == 1 - return x - else - N = length(x) - L = floor(Int, N / s) - ys = zeros(T, L) - - for t = 1:L - inds = ((t - 1)*s + 1):(t * s) - ys[t] = @views f(x[inds], args...; kwargs...) - end - return ys - end -end - -function downsample(method::Composite, x::AbstractVector{T}, s::Int, args...; kwargs...) where T - f = method.f - - if s == 1 - return [x] - else - N = length(x) - # note: there must be a typo or error in Wu et al. (2013), because if we use - # floor(N / s), as indicated in their paper, we'll get out of bounds errors. - # The number of samples needs to take into consideratio how many unique ways - # there are of selecting non-overlapping windows of length `s`. Hence, - # we use floor((N - s + 1) / s) instead. - L = floor(Int, (N - s + 1) / s) - ys = [zeros(T, L) for i = 1:s] - for k = 1:s - for t = 1:L - inds = ((t - 1)*s + k):(t * s + k - 1) - ys[k][t] = @views f(x[inds], args...; kwargs...) - end - end - return ys - end -end +downsample(method::MultiScaleAlgorithm, x, s::Int) """ - multiscale(e::Entropy, alg::MultiScaleAlgorithm, x::AbstractVector, - est::ProbabilitiesEstimator; - scalemax::Int = 8, - normalize = false) + multiscale(e::Entropy, alg, x, est; maxscale = 8, normalize = false) + multiscale(c::ComplexityMeasure, alg, x; maxscale = 8, normalize = false) -Compute the multi-scale entropy (Costa et al., 2002)[^Costa2002] of type `e` for -timeseries `x`, using downsampling algorithm `alg`, and probabilities estimator `est`. +Compute the multi-scale entropy `e` with probabilities estimator `est`, or the multi-scale +complexity measure `c`, for timeseries `x`. -If `normalize == true`, then compute normalized entropy. +## Description - multiscale(c::ComplexityMeasure, alg::MultiScaleAlgorithm, x::AbstractVector; - scalemax::Int = 8, - normalize = false) +Utilizes [`downsample`](@ref) to compute the entropy/complexity of coarse-grained, +downsampled versions of `x` for scale factors `1:maxscale`. If `N = length(x)`, then the +length of the most severely downsampled version of `x` is `N ÷ scalemax`, while for scale +factor `1`, the original time series is considered. -Compute the multi-scale complexity measure `c` for the timeseries `x`, using -using downsampling algorithm `alg`. +If relevant [`MultiscaleAlgorithm`](@ref)s and corresponding [`downsample`](@ref) +methods are defined, +- the first method generalizes all multi-scale entropy estimators where actual entropies + (i.e. functionals of explicitly estimated probability distributions) are computed, and +- the second method generalizes all multi-scale complexity ("entropy-like") measures. -If `normalize == true`, then compute the normalized version of the complexity measure. +## Arguments -## Description +- `e::Entropy`. A valid [entropy type](@ref entropies_list), i.e. `Shannon()` or `Renyi()`. +- `e::Complexity`. A valid [complexity measure](@ref complexity_measures), i.e. + `ReverseDispersion()` +- `alg::MultiScaleAlgorithm`. A valid [multiscale algorithm](@ref multiscale_algorithms), + i.e. `Regular()` or `Composite()`, which determines how down-sampling/coarse-graining + is performed. +- `x`. The input data. Usually a timeseries. +- `est::ProbabilitiesEstimator`. A [probabilities estimator](@ref probabilities_estimators), + which determines how probabilities are estimated for each downsampled time series. + +## Keyword Arguments -Utilizes [`downsample`](@ref) to compute the entropy of coarse-grained, downsampled -versions of `x` for scale factors `1:maxscale`. If `N = length(x)`, then the length of the -most severely downsampled version of `x` is `N ÷ scalemax`, while for scale factor `1`, -the original time series is considered. +- `scalemax::Int`. The maximum number of scales (i.e. levels of downsampling). +- `normalize::Bool`. If `normalize == true`, then (if possible) compute normalized + entropy/complexity. If `normalize == false`, then compute the non-normalized measure. [^Costa2002]: Costa, M., Goldberger, A. L., & Peng, C. K. (2002). Multiscale entropy analysis of complex physiologic time series. Physical review letters, 89(6), 068102. """ function multiscale end -function multiscale(e::Entropy, alg::Regular, x::AbstractVector, est::ProbabilitiesEstimator; - maxscale::Int = 8, normalize = false) - - downscaled_timeseries = [downsample(alg, x, s) for s in 1:maxscale] - hs = zeros(Float64, maxscale) - if normalize - hs = entropy_normalized.(Ref(e), downscaled_timeseries, Ref(est)) - else - hs .= entropy.(Ref(e), downscaled_timeseries, Ref(est)) - end - - return hs -end - -function multiscale(e::Entropy, alg::Composite, x::AbstractVector, est::ProbabilitiesEstimator; - maxscale::Int = 10, normalize = false) - - downscaled_timeseries = [downsample(alg, x, s) for s in 1:maxscale] - hs = zeros(Float64, maxscale) - for s in 1:maxscale - if normalize - hs[s] = mean(entropy_normalized.(Ref(e), downscaled_timeseries[s], Ref(est))) - else - hs[s] = mean(entropy.(Ref(e), downscaled_timeseries[s], Ref(est))) - end - end - - return hs -end - -function multiscale(e::ComplexityMeasure, alg::Regular, x::AbstractVector; - maxscale::Int = 8, normalize = false) - - downscaled_timeseries = [downsample(alg, x, s) for s in 1:maxscale] - complexities = zeros(Float64, maxscale) - if normalize - complexities = complexity_normalized.(Ref(e), downscaled_timeseries) - else - complexities .= complexity.(Ref(e), downscaled_timeseries) - end - - return complexities -end - -function multiscale(e::ComplexityMeasure, alg::Composite, x::AbstractVector; - maxscale::Int = 10, normalize = false) - - downscaled_timeseries = [downsample(alg, x, s) for s in 1:maxscale] - complexities = zeros(Float64, maxscale) - for s in 1:maxscale - if normalize - complexities[s] = - mean(complexity_normalized.(Ref(e), downscaled_timeseries[s], Ref(est))) - else - complexities[s] = mean(entropy.(Ref(e), downscaled_timeseries[s], Ref(est))) - end - end - - return complexities -end +include("multiscale/regular.jl") +include("multiscale/composite.jl") diff --git a/src/multiscale/composite.jl b/src/multiscale/composite.jl new file mode 100644 index 000000000..04e1bb79c --- /dev/null +++ b/src/multiscale/composite.jl @@ -0,0 +1,107 @@ +export Composite + +""" + Composite <: MultiScaleAlgorithm + Composite(; f::Function = Statistics.mean) + +Composite multi-scale algorithm for multiscale entropy analysis (Wu et al., +2013)[^Wu2013], used, with [`multiscale`](@ref) to compute, for example, composite +multiscale entropy (CMSE). + +## Description + +Given a scalar-valued input time series `x`, the composite multiscale algorithm, +like [`Regular`](@ref), downsamples and coarse-grains `x` by splitting it into +non-overlapping windows of length `s`, and then constructing downsampled time series by +applying the function `f` to each of the resulting length-`s` windows. + +However, Wu et al. (2013)[^Wu2013] realized that for each scale `s`, there are actually `s` +different ways of selecting windows, depending on where indexing starts/ends. +These `s` different downsampled time series `D_t(s, f)` at each scale `s` are +constructed as follows: + +```math +\\{ D_{k}(s) \\} = \\{ D_{t, k}(s) \\}_{t = 1}^{L}, = \\{ f \\left( \\bf x_{t, k} \\right) \\} = +\\left\\{ + {f\\left( (x_i)_{i = (t - 1)s + k}^{ts + k - 1} \\right)} +\\right\\}_{t = 1}^{L}, +``` + +where `L = floor((N - s + 1) / s)` and `1 ≤ k ≤ s`, such that ``D_{i, k}(s)`` is the `i`-th +element of the `k`-th downsampled time series at scale `s`. + +Finally, compute ``\\dfrac{1}{s} \\sum_{k = 1}^s g(D_{k}(s))``, where `g` is some summary +function, for example [`entropy`](@ref) or [`complexity`](@ref). + +!!! note "Relation to Regular" + The downsampled time series ``D_{t, 1}(s)`` constructed using the composite + multiscale method is equivalent to the downsampled time series ``D_{t}(s)`` constructed + using the [`Regular`](@ref) method, for which `k == 1` is fixed, such that only + a single time series is returned. + +See also: [`Regular`](@ref). + +[^Wu2013]: Wu, S. D., Wu, C. W., Lin, S. G., Wang, C. C., & Lee, K. Y. (2013). Time series + analysis using composite multiscale entropy. Entropy, 15(3), 1069-1084. +""" +Base.@kwdef struct Composite <: MultiScaleAlgorithm + f::Function = Statistics.mean +end + +function downsample(method::Composite, x::AbstractVector{T}, s::Int, args...; kwargs...) where T + f = method.f + + if s == 1 + return [x] + else + N = length(x) + # note: there must be a typo or error in Wu et al. (2013), because if we use + # floor(N / s), as indicated in their paper, we'll get out of bounds errors. + # The number of samples needs to take into consideratio how many unique ways + # there are of selecting non-overlapping windows of length `s`. Hence, + # we use floor((N - s + 1) / s) instead. + L = floor(Int, (N - s + 1) / s) + ys = [zeros(T, L) for i = 1:s] + for k = 1:s + for t = 1:L + inds = ((t - 1)*s + k):(t * s + k - 1) + ys[k][t] = @views f(x[inds], args...; kwargs...) + end + end + return ys + end +end + +function multiscale(e::Entropy, alg::Composite, x::AbstractVector, est::ProbabilitiesEstimator; + maxscale::Int = 10, normalize = false) + + downscaled_timeseries = [downsample(alg, x, s) for s in 1:maxscale] + hs = zeros(Float64, maxscale) + for s in 1:maxscale + if normalize + hs[s] = mean(entropy_normalized.(Ref(e), downscaled_timeseries[s], Ref(est))) + else + hs[s] = mean(entropy.(Ref(e), downscaled_timeseries[s], Ref(est))) + end + end + + return hs +end + + +function multiscale(e::ComplexityMeasure, alg::Composite, x::AbstractVector; + maxscale::Int = 10, normalize = false) + + downscaled_timeseries = [downsample(alg, x, s) for s in 1:maxscale] + complexities = zeros(Float64, maxscale) + for s in 1:maxscale + if normalize + complexities[s] = + mean(complexity_normalized.(Ref(e), downscaled_timeseries[s], Ref(est))) + else + complexities[s] = mean(entropy.(Ref(e), downscaled_timeseries[s], Ref(est))) + end + end + + return complexities +end diff --git a/src/multiscale/regular.jl b/src/multiscale/regular.jl new file mode 100644 index 000000000..506f338be --- /dev/null +++ b/src/multiscale/regular.jl @@ -0,0 +1,93 @@ +export Regular + +""" + Regular <: MultiScaleAlgorithm + Regular(; f::Function = Statistics.mean) + +The original multi-scale algorithm for multiscale entropy analysis (Costa et al., +2022)[^Costa2002], which yields a single downsampled time series per scale `s`. + +## Description + +Given a scalar-valued input time series `x`, the `Regular` multiscale algorithm downsamples +and coarse-grains `x` by splitting it into non-overlapping windows of length `s`, and +then constructing a new downsampled time series ``D_t(s, f)`` by applying the function `f` +to each of the resulting length-`s` windows. + +The downsampled time series `D_t(s)` with `t ∈ [1, 2, …, L]`, where `L = floor(N / s)`, +is given by: + +```math +\\{ D_t(s, f) \\}_{t = 1}^{L} = \\left\\{ f \\left( \\bf x_t \\right) \\right\\}_{t = 1}^{L} = +\\left\\{ + {f\\left( (x_i)_{i = (t - 1)s + 1}^{ts} \\right)} +\\right\\}_{t = 1}^{L} +``` + +where `f` is some summary statistic applied to the length-`ts-((t - 1)s + 1)` tuples `xₖ`. +Different choices of `f` have yield different multiscale methods appearing in the +literature. For example: + +- `f == Statistics.mean` yields the original first-moment multiscale sample entropy (Costa + et al., 2002)[^Costa2002]. +- `f == Statistics.var` yields the generalized multiscale sample entropy (Costa & + Goldberger, 2015)[^Costa2015], which uses the second-moment (variance) instead of the + mean. + +See also: [`Composite`](@ref). + +[^Costa2002]: Costa, M., Goldberger, A. L., & Peng, C. K. (2002). Multiscale entropy + analysis of complex physiologic time series. Physical review letters, 89(6), 068102. +[^Costa2015]: Costa, M. D., & Goldberger, A. L. (2015). Generalized multiscale entropy + analysis: Application to quantifying the complex volatility of human heartbeat time + series. Entropy, 17(3), 1197-1203. +""" +Base.@kwdef struct Regular <: MultiScaleAlgorithm + f::Function = Statistics.mean +end + +function downsample(method::Regular, x::AbstractVector{T}, s::Int, args...; kwargs...) where T + f = method.f + + if s == 1 + return x + else + N = length(x) + L = floor(Int, N / s) + ys = zeros(T, L) + + for t = 1:L + inds = ((t - 1)*s + 1):(t * s) + ys[t] = @views f(x[inds], args...; kwargs...) + end + return ys + end +end + +function multiscale(e::Entropy, alg::Regular, x::AbstractVector, est::ProbabilitiesEstimator; + maxscale::Int = 8, normalize = false) + + downscaled_timeseries = [downsample(alg, x, s) for s in 1:maxscale] + hs = zeros(Float64, maxscale) + if normalize + hs = entropy_normalized.(Ref(e), downscaled_timeseries, Ref(est)) + else + hs .= entropy.(Ref(e), downscaled_timeseries, Ref(est)) + end + + return hs +end + +function multiscale(e::ComplexityMeasure, alg::Regular, x::AbstractVector; + maxscale::Int = 8, normalize = false) + + downscaled_timeseries = [downsample(alg, x, s) for s in 1:maxscale] + complexities = zeros(Float64, maxscale) + if normalize + complexities = complexity_normalized.(Ref(e), downscaled_timeseries) + else + complexities .= complexity.(Ref(e), downscaled_timeseries) + end + + return complexities +end diff --git a/test/multiscale/multiscale_entropy.jl b/test/multiscale/multiscale_entropy.jl index 4693d44e2..264ea597a 100644 --- a/test/multiscale/multiscale_entropy.jl +++ b/test/multiscale/multiscale_entropy.jl @@ -1,7 +1,35 @@ est = Dispersion() x = rand(100) +maxscale = 5 + e = Shannon() +mr = multiscale(e, Regular(), x, est; maxscale, normalize = false) +mc = multiscale(e, Composite(), x, est; maxscale, normalize = false) +mrn = multiscale(e, Regular(), x, est; maxscale, normalize = true) +mcn = multiscale(e, Composite(), x, est; maxscale, normalize = true) +@test mr isa Vector{T} where T <: Real +@test mc isa Vector{T} where T <: Real +@test mrn isa Vector{T} where T <: Real +@test mcn isa Vector{T} where T <: Real +@test length(mr) == 5 +@test length(mc) == 5 +@test length(mrn) == 5 +@test length(mcn) == 5 + +c = ReverseDispersion() +mcr = multiscale(c, Regular(), x; maxscale, normalize = false) +mcc = multiscale(c, Composite(), x, est; maxscale, normalize = false) +mcrn = multiscale(c, Regular(), x; maxscale, normalize = true) +mccn = multiscale(c, Composite(), x, est; maxscale, normalize = true) +@test mcr isa Vector{T} where T <: Real +@test mcc isa Vector{T} where T <: Real +@test mcrn isa Vector{T} where T <: Real +@test mccn isa Vector{T} where T <: Real +@test length(mr) == 5 +@test length(mc) == 5 +@test length(mrn) == 5 +@test length(mcn) == 5 -# TODO: test actual outcomes too. Make some analytical examples. -@test multiscale(e, Regular(), x, est) isa Vector{T} where T <: Real -@test multiscale(e, Composite(), x, est) isa Vector{T} where T <: Real +# TODO: To verify that downsampling algorithms work as expected, we need to make a few +# simple analytical examples too. There are no such examples in the original papers, +# so we'll have to make some ourselves. From 109864eb2f26eb0705550b002c2e0096dbaab4e1 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Sun, 23 Oct 2022 20:19:37 +0200 Subject: [PATCH 11/27] Fix tests --- src/multiscale/composite.jl | 4 ++-- test/multiscale/multiscale_entropy.jl | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/multiscale/composite.jl b/src/multiscale/composite.jl index 04e1bb79c..25615e825 100644 --- a/src/multiscale/composite.jl +++ b/src/multiscale/composite.jl @@ -97,9 +97,9 @@ function multiscale(e::ComplexityMeasure, alg::Composite, x::AbstractVector; for s in 1:maxscale if normalize complexities[s] = - mean(complexity_normalized.(Ref(e), downscaled_timeseries[s], Ref(est))) + mean(complexity_normalized.(Ref(e), downscaled_timeseries[s])) else - complexities[s] = mean(entropy.(Ref(e), downscaled_timeseries[s], Ref(est))) + complexities[s] = mean(complexity.(Ref(e), downscaled_timeseries[s])) end end diff --git a/test/multiscale/multiscale_entropy.jl b/test/multiscale/multiscale_entropy.jl index 264ea597a..f048d45cc 100644 --- a/test/multiscale/multiscale_entropy.jl +++ b/test/multiscale/multiscale_entropy.jl @@ -17,10 +17,10 @@ mcn = multiscale(e, Composite(), x, est; maxscale, normalize = true) @test length(mcn) == 5 c = ReverseDispersion() -mcr = multiscale(c, Regular(), x; maxscale, normalize = false) -mcc = multiscale(c, Composite(), x, est; maxscale, normalize = false) -mcrn = multiscale(c, Regular(), x; maxscale, normalize = true) -mccn = multiscale(c, Composite(), x, est; maxscale, normalize = true) +mcr = multiscale(c, Regular(), x; maxscale, normalize = false) +mcc = multiscale(c, Composite(), x; maxscale, normalize = false) +mcrn = multiscale(c, Regular(), x; maxscale, normalize = true) +mccn = multiscale(c, Composite(), x; maxscale, normalize = true) @test mcr isa Vector{T} where T <: Real @test mcc isa Vector{T} where T <: Real @test mcrn isa Vector{T} where T <: Real From 1916c2e90d3949e65ccff5072dddfbe4c93f85ba Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Sun, 23 Oct 2022 20:19:52 +0200 Subject: [PATCH 12/27] Add generic error messages --- src/multiscale.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/multiscale.jl b/src/multiscale.jl index 12cb9ae29..4b7721eec 100644 --- a/src/multiscale.jl +++ b/src/multiscale.jl @@ -73,5 +73,17 @@ methods are defined, """ function multiscale end +function multiscale(e::Entropy, x, est::ProbabilitiesEstimator) + msg = "`multiscale` entropy not implemented for $e $est on data type $(typeof(x))" + throw(ArgumentError(msg)) +end + +function multiscale(c::ComplexityMeasure, x) + msg = "`multiscale` complexity not implemented for $c and data type $(typeof(x))" + throw(ArgumentError(msg)) +end + + + include("multiscale/regular.jl") include("multiscale/composite.jl") From 89b7bbc18bee6234130bd70cda453ec83fd044c7 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Sun, 23 Oct 2022 20:20:16 +0200 Subject: [PATCH 13/27] Multivariate downsampling --- src/multiscale.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/multiscale.jl b/src/multiscale.jl index 4b7721eec..a24d0f45c 100644 --- a/src/multiscale.jl +++ b/src/multiscale.jl @@ -30,6 +30,9 @@ The return type depends on `algorithm`. For example: """ downsample(method::MultiScaleAlgorithm, x, s::Int) +downsample(alg::MultiScaleAlgorithm, x::AbstractDataset, args...; kwargs...) = + Dataset(map(t -> downsample(alg, t, args...; kwargs...), columns(x))...) + """ multiscale(e::Entropy, alg, x, est; maxscale = 8, normalize = false) multiscale(c::ComplexityMeasure, alg, x; maxscale = 8, normalize = false) From 9c8b5591754271e56b517343563b737628a8f7d9 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Mon, 24 Oct 2022 00:02:09 +0200 Subject: [PATCH 14/27] Multiscale normalized --- src/multiscale.jl | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/multiscale.jl b/src/multiscale.jl index a24d0f45c..27388ecd3 100644 --- a/src/multiscale.jl +++ b/src/multiscale.jl @@ -76,6 +76,15 @@ methods are defined, """ function multiscale end +""" + multiscale(e::Entropy, alg, x, est; maxscale = 8) + multiscale(c::ComplexityMeasure, alg, x; maxscale = 8) + +The same as [`multiscale`](@ref), but normalizes the entropy or complexity measure. +""" +function multiscale_normalized end + + function multiscale(e::Entropy, x, est::ProbabilitiesEstimator) msg = "`multiscale` entropy not implemented for $e $est on data type $(typeof(x))" throw(ArgumentError(msg)) From a757a42318e756a00ede2cfb893b0bea5ad73a2b Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Mon, 24 Oct 2022 00:02:20 +0200 Subject: [PATCH 15/27] Test downsampling explicitly --- test/multiscale/downsampling.jl | 36 +++++++++++++++++++ .../{multiscale_entropy.jl => multiscale.jl} | 0 test/runtests.jl | 3 +- 3 files changed, 38 insertions(+), 1 deletion(-) create mode 100644 test/multiscale/downsampling.jl rename test/multiscale/{multiscale_entropy.jl => multiscale.jl} (100%) diff --git a/test/multiscale/downsampling.jl b/test/multiscale/downsampling.jl new file mode 100644 index 000000000..82d242d4e --- /dev/null +++ b/test/multiscale/downsampling.jl @@ -0,0 +1,36 @@ +using Statistics +x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] + +############################################################## +# Regular downsampling +############################################################## +r = Regular(f = mean) + +# just returns original time series +@test downsample(r, x, 1) == x + +# for scale 2, downsampling indices for `x` are 1:2, 3:4, 5:6, 7:8, 9:10, so we should get +# the means of consecutive running pairs +@test downsample(r, x, 2) == [1.5, 3.5, 5.5, 7.5, 9.5] + +# for scale 3, downsampling indices for `x` are 1:3, 4:6, 7:9. Take the means: +@test downsample(r, x, 3) == [2.0, 5.0, 8.0] + +# for scale 4, downsampling indices for `x` are 1:4, 5:8. Take the means: +@test downsample(r, x, 4) == [2.5, 6.5] + +# for scale 5, downsampling indices for `x` are 1:5, 6:10. Take the means: +@test downsample(r, x, 5) == [3.0, 8.0] + +# for scale 6, only a single set of sampling indices is possible - 1:6 +@test downsample(r, x, 6) == [3.5] +@test downsample(r, x, 6) == [3.5] + +# For scales larger than the number of elements, there should be no values in the +# downsampled time series. +@test isempty(downsample(r, x, length(x) + 1)) + +############################################################## +# Regular downsampling +############################################################## +c = Composite(f = mean) diff --git a/test/multiscale/multiscale_entropy.jl b/test/multiscale/multiscale.jl similarity index 100% rename from test/multiscale/multiscale_entropy.jl rename to test/multiscale/multiscale.jl diff --git a/test/runtests.jl b/test/runtests.jl index 9ff89b589..1334a545f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -28,7 +28,8 @@ testfile(file, testname=defaultname(file)) = @testset "$testname" begin; include testfile("entropies/nearest_neighbors_direct.jl") # Multiscale analysis - testfile("multiscale/multiscale_entropy.jl") + testfile("multiscale/downsampling.jl") + testfile("multiscale/multiscale.jl") # Various testfile("complexity_measures/complexity_measures.jl") From 00b12458177d6863105548015620badf26d0319f Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Mon, 24 Oct 2022 00:02:42 +0200 Subject: [PATCH 16/27] Always return float-vectors --- src/multiscale/composite.jl | 6 ++++-- src/multiscale/regular.jl | 5 ++++- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/src/multiscale/composite.jl b/src/multiscale/composite.jl index 25615e825..42282f9a9 100644 --- a/src/multiscale/composite.jl +++ b/src/multiscale/composite.jl @@ -50,9 +50,10 @@ end function downsample(method::Composite, x::AbstractVector{T}, s::Int, args...; kwargs...) where T f = method.f + ET = eltype(one(1.0)) # consistently return floats, even if input is e.g. integer-valued if s == 1 - return [x] + return ET.(x) else N = length(x) # note: there must be a typo or error in Wu et al. (2013), because if we use @@ -61,7 +62,7 @@ function downsample(method::Composite, x::AbstractVector{T}, s::Int, args...; kw # there are of selecting non-overlapping windows of length `s`. Hence, # we use floor((N - s + 1) / s) instead. L = floor(Int, (N - s + 1) / s) - ys = [zeros(T, L) for i = 1:s] + ys = [zeros(ET, L) for i = 1:s] for k = 1:s for t = 1:L inds = ((t - 1)*s + k):(t * s + k - 1) @@ -72,6 +73,7 @@ function downsample(method::Composite, x::AbstractVector{T}, s::Int, args...; kw end end +# TODO: make a separate multiscale_normalized? function multiscale(e::Entropy, alg::Composite, x::AbstractVector, est::ProbabilitiesEstimator; maxscale::Int = 10, normalize = false) diff --git a/src/multiscale/regular.jl b/src/multiscale/regular.jl index 506f338be..199e207b3 100644 --- a/src/multiscale/regular.jl +++ b/src/multiscale/regular.jl @@ -49,15 +49,17 @@ end function downsample(method::Regular, x::AbstractVector{T}, s::Int, args...; kwargs...) where T f = method.f + ET = eltype(one(1.0)) # consistently return floats, even if input is e.g. integer-valued if s == 1 return x else N = length(x) L = floor(Int, N / s) - ys = zeros(T, L) + ys = zeros(ET, L) for t = 1:L inds = ((t - 1)*s + 1):(t * s) + @show inds ys[t] = @views f(x[inds], args...; kwargs...) end return ys @@ -77,6 +79,7 @@ function multiscale(e::Entropy, alg::Regular, x::AbstractVector, est::Probabilit return hs end +# TODO: make a separate multiscale_normalized? function multiscale(e::ComplexityMeasure, alg::Regular, x::AbstractVector; maxscale::Int = 8, normalize = false) From f4f98986cf707ba971b1b652996a75bfdeb4c5e3 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Thu, 3 Nov 2022 16:06:25 +0100 Subject: [PATCH 17/27] Update API Multiscale API --- docs/make.jl | 1 + docs/src/entropies.md | 2 +- docs/src/multiscale.md | 22 ++++++++++++++---- docs/src/probabilities.md | 1 + src/multiscale.jl | 42 +++++++++++++++-------------------- src/multiscale/composite.jl | 30 +++++++++---------------- src/multiscale/regular.jl | 30 +++++++------------------ test/multiscale/multiscale.jl | 28 +++++++---------------- 8 files changed, 66 insertions(+), 90 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 73a219ccb..6178452fd 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -37,6 +37,7 @@ ENTROPIES_PAGES = [ "probabilities.md", "entropies.md", "complexity_measures.md", + "multiscale.md", "examples.md", "utils.md", "devdocs.md", diff --git a/docs/src/entropies.md b/docs/src/entropies.md index 1dbc6bd92..b445e829f 100644 --- a/docs/src/entropies.md +++ b/docs/src/entropies.md @@ -1,4 +1,4 @@ -# Entropies +# [Entropies](@id entropies) ```@docs entropy diff --git a/docs/src/multiscale.md b/docs/src/multiscale.md index cab5f78d1..9e238aceb 100644 --- a/docs/src/multiscale.md +++ b/docs/src/multiscale.md @@ -1,7 +1,17 @@ +# Multiscale + ## Multiscale API +The multiscale API is defined by the [multiscale algorithms](@ref multiscale_algorithms) +defined below, and the following functions. + +- [`multiscale`](@ref) +- [`multiscale_normalized`](@ref) +- [`downsample`](@ref) + ```@docs multiscale +multiscale_normalized downsample ``` @@ -34,14 +44,18 @@ Regular Composite ``` -[^Costa2002]: Costa, M., Goldberger, A. L., & Peng, C. K. (2002). Multiscale entropy +[^Costa2002]: + Costa, M., Goldberger, A. L., & Peng, C. K. (2002). Multiscale entropy analysis of complex physiologic time series. Physical review letters, 89(6), 068102. -[^Costa2015]: Costa, M. D., & Goldberger, A. L. (2015). Generalized multiscale entropy +[^Costa2015]: + Costa, M. D., & Goldberger, A. L. (2015). Generalized multiscale entropy analysis: Application to quantifying the complex volatility of human heartbeat time series. Entropy, 17(3), 1197-1203. -[^Azami2017]: Azami, H., Rostaghi, M., Abásolo, D., & Escudero, J. (2017). Refined +[^Azami2017]: + Azami, H., Rostaghi, M., Abásolo, D., & Escudero, J. (2017). Refined composite multiscale dispersion entropy and its application to biomedical signals. IEEE Transactions on Biomedical Engineering, 64(12), 2872-2879. -[^Morabito2012]: Morabito, F. C., Labate, D., Foresta, F. L., Bramanti, A., Morabito, G., +[^Morabito2012]: + Morabito, F. C., Labate, D., Foresta, F. L., Bramanti, A., Morabito, G., & Palamara, I. (2012). Multivariate multi-scale permutation entropy for complexity analysis of Alzheimer’s disease EEG. Entropy, 14(7), 1186-1202. diff --git a/docs/src/probabilities.md b/docs/src/probabilities.md index f6a739eaf..51dd6af8d 100644 --- a/docs/src/probabilities.md +++ b/docs/src/probabilities.md @@ -3,6 +3,7 @@ ## Probabilities API The probabilities API is defined by + - [ProbabilitiesEstimator](@ref) - [probabilities](@ref) - [probabilities_and_outcomes](@ref) diff --git a/src/multiscale.jl b/src/multiscale.jl index 27388ecd3..5d7349a78 100644 --- a/src/multiscale.jl +++ b/src/multiscale.jl @@ -1,10 +1,8 @@ -# This file contains an API for multiscale (coarse-grained/downsampled) computations of -# entropy and various complexity measures on time series. -# Can be used both to compute actual entropies (i.e. diversity entropy), or -# complexity measures (sample entropy, approximate entropy). +# This file contains an API for multiscale (coarse-grained/downsampled) computations. using Statistics export multiscale +export multiscale_normalized export downsample export MultiScaleAlgorithm @@ -33,12 +31,14 @@ downsample(method::MultiScaleAlgorithm, x, s::Int) downsample(alg::MultiScaleAlgorithm, x::AbstractDataset, args...; kwargs...) = Dataset(map(t -> downsample(alg, t, args...; kwargs...), columns(x))...) + +function multiscale end +function multiscale_normalized end + """ multiscale(e::Entropy, alg, x, est; maxscale = 8, normalize = false) - multiscale(c::ComplexityMeasure, alg, x; maxscale = 8, normalize = false) -Compute the multi-scale entropy `e` with probabilities estimator `est`, or the multi-scale -complexity measure `c`, for timeseries `x`. +Compute the multi-scale entropy `e` with probabilities estimator `est` for timeseries `x`. ## Description @@ -47,7 +47,7 @@ downsampled versions of `x` for scale factors `1:maxscale`. If `N = length(x)`, length of the most severely downsampled version of `x` is `N ÷ scalemax`, while for scale factor `1`, the original time series is considered. -If relevant [`MultiscaleAlgorithm`](@ref)s and corresponding [`downsample`](@ref) +If relevant [`MultiScaleAlgorithm`](@ref)s and corresponding [`downsample`](@ref) methods are defined, - the first method generalizes all multi-scale entropy estimators where actual entropies (i.e. functionals of explicitly estimated probability distributions) are computed, and @@ -55,9 +55,7 @@ methods are defined, ## Arguments -- `e::Entropy`. A valid [entropy type](@ref entropies_list), i.e. `Shannon()` or `Renyi()`. -- `e::Complexity`. A valid [complexity measure](@ref complexity_measures), i.e. - `ReverseDispersion()` +- `e::Entropy`. A valid [entropy type](@ref entropies), i.e. `Shannon()` or `Renyi()`. - `alg::MultiScaleAlgorithm`. A valid [multiscale algorithm](@ref multiscale_algorithms), i.e. `Regular()` or `Composite()`, which determines how down-sampling/coarse-graining is performed. @@ -74,26 +72,22 @@ methods are defined, [^Costa2002]: Costa, M., Goldberger, A. L., & Peng, C. K. (2002). Multiscale entropy analysis of complex physiologic time series. Physical review letters, 89(6), 068102. """ -function multiscale end +function multiscale(e::Entropy, alg::MultiScaleAlgorithm, x, est::ProbabilitiesEstimator) + msg = "`multiscale` entropy not implemented for $e $est on data type $(typeof(x))" + throw(ArgumentError(msg)) +end """ - multiscale(e::Entropy, alg, x, est; maxscale = 8) - multiscale(c::ComplexityMeasure, alg, x; maxscale = 8) + multiscale_normalized(e::Entropy, alg, x, est; maxscale = 8) -The same as [`multiscale`](@ref), but normalizes the entropy or complexity measure. +The same as [`multiscale`](@ref), but normalizes the estimated quantity according """ -function multiscale_normalized end - - -function multiscale(e::Entropy, x, est::ProbabilitiesEstimator) - msg = "`multiscale` entropy not implemented for $e $est on data type $(typeof(x))" +function multiscale_normalized(e::Entropy, alg::MultiScaleAlgorithm, x, + est::ProbabilitiesEstimator) + msg = "`multiscale_normalized` not implemented for $e $est on data type $(typeof(x))" throw(ArgumentError(msg)) end -function multiscale(c::ComplexityMeasure, x) - msg = "`multiscale` complexity not implemented for $c and data type $(typeof(x))" - throw(ArgumentError(msg)) -end diff --git a/src/multiscale/composite.jl b/src/multiscale/composite.jl index 42282f9a9..6ebf95cd7 100644 --- a/src/multiscale/composite.jl +++ b/src/multiscale/composite.jl @@ -53,7 +53,7 @@ function downsample(method::Composite, x::AbstractVector{T}, s::Int, args...; kw ET = eltype(one(1.0)) # consistently return floats, even if input is e.g. integer-valued if s == 1 - return ET.(x) + return [ET.(x)] else N = length(x) # note: there must be a typo or error in Wu et al. (2013), because if we use @@ -74,36 +74,28 @@ function downsample(method::Composite, x::AbstractVector{T}, s::Int, args...; kw end # TODO: make a separate multiscale_normalized? -function multiscale(e::Entropy, alg::Composite, x::AbstractVector, est::ProbabilitiesEstimator; - maxscale::Int = 10, normalize = false) +function multiscale(e::Entropy, alg::Composite, x::AbstractVector, + est::ProbabilitiesEstimator; + maxscale::Int = 10) downscaled_timeseries = [downsample(alg, x, s) for s in 1:maxscale] hs = zeros(Float64, maxscale) for s in 1:maxscale - if normalize - hs[s] = mean(entropy_normalized.(Ref(e), downscaled_timeseries[s], Ref(est))) - else - hs[s] = mean(entropy.(Ref(e), downscaled_timeseries[s], Ref(est))) - end + hs[s] = mean(entropy.(Ref(e), downscaled_timeseries[s], Ref(est))) end return hs end - -function multiscale(e::ComplexityMeasure, alg::Composite, x::AbstractVector; - maxscale::Int = 10, normalize = false) +function multiscale_normalized(e::Entropy, alg::Composite, x::AbstractVector, + est::ProbabilitiesEstimator; + maxscale::Int = 10) downscaled_timeseries = [downsample(alg, x, s) for s in 1:maxscale] - complexities = zeros(Float64, maxscale) + hs = zeros(Float64, maxscale) for s in 1:maxscale - if normalize - complexities[s] = - mean(complexity_normalized.(Ref(e), downscaled_timeseries[s])) - else - complexities[s] = mean(complexity.(Ref(e), downscaled_timeseries[s])) - end + #hs[s] = mean(entropy_normalized.(Ref(e), downscaled_timeseries[s], Ref(est))) end - return complexities + return hs end diff --git a/src/multiscale/regular.jl b/src/multiscale/regular.jl index 199e207b3..fbb03fa9d 100644 --- a/src/multiscale/regular.jl +++ b/src/multiscale/regular.jl @@ -59,38 +59,24 @@ function downsample(method::Regular, x::AbstractVector{T}, s::Int, args...; kwar for t = 1:L inds = ((t - 1)*s + 1):(t * s) - @show inds ys[t] = @views f(x[inds], args...; kwargs...) end return ys end end -function multiscale(e::Entropy, alg::Regular, x::AbstractVector, est::ProbabilitiesEstimator; - maxscale::Int = 8, normalize = false) +function multiscale(e::Entropy, alg::Regular, x::AbstractVector, + est::ProbabilitiesEstimator; + maxscale::Int = 8) downscaled_timeseries = [downsample(alg, x, s) for s in 1:maxscale] - hs = zeros(Float64, maxscale) - if normalize - hs = entropy_normalized.(Ref(e), downscaled_timeseries, Ref(est)) - else - hs .= entropy.(Ref(e), downscaled_timeseries, Ref(est)) - end - - return hs + return entropy.(Ref(e), downscaled_timeseries, Ref(est)) end -# TODO: make a separate multiscale_normalized? -function multiscale(e::ComplexityMeasure, alg::Regular, x::AbstractVector; - maxscale::Int = 8, normalize = false) +function multiscale_normalized(e::Entropy, alg::Regular, x::AbstractVector, + est::ProbabilitiesEstimator; + maxscale::Int = 8) downscaled_timeseries = [downsample(alg, x, s) for s in 1:maxscale] - complexities = zeros(Float64, maxscale) - if normalize - complexities = complexity_normalized.(Ref(e), downscaled_timeseries) - else - complexities .= complexity.(Ref(e), downscaled_timeseries) - end - - return complexities + return entropy_normalized.(Ref(e), downscaled_timeseries, Ref(est)) end diff --git a/test/multiscale/multiscale.jl b/test/multiscale/multiscale.jl index f048d45cc..ae9a02446 100644 --- a/test/multiscale/multiscale.jl +++ b/test/multiscale/multiscale.jl @@ -1,33 +1,21 @@ est = Dispersion() x = rand(100) maxscale = 5 - e = Shannon() -mr = multiscale(e, Regular(), x, est; maxscale, normalize = false) -mc = multiscale(e, Composite(), x, est; maxscale, normalize = false) -mrn = multiscale(e, Regular(), x, est; maxscale, normalize = true) -mcn = multiscale(e, Composite(), x, est; maxscale, normalize = true) + +# Generic tests. +mr = multiscale(e, Regular(), x, est; maxscale) +mrn = multiscale_normalized(e, Regular(), x, est; maxscale) @test mr isa Vector{T} where T <: Real -@test mc isa Vector{T} where T <: Real @test mrn isa Vector{T} where T <: Real -@test mcn isa Vector{T} where T <: Real @test length(mr) == 5 -@test length(mc) == 5 @test length(mrn) == 5 -@test length(mcn) == 5 -c = ReverseDispersion() -mcr = multiscale(c, Regular(), x; maxscale, normalize = false) -mcc = multiscale(c, Composite(), x; maxscale, normalize = false) -mcrn = multiscale(c, Regular(), x; maxscale, normalize = true) -mccn = multiscale(c, Composite(), x; maxscale, normalize = true) -@test mcr isa Vector{T} where T <: Real -@test mcc isa Vector{T} where T <: Real -@test mcrn isa Vector{T} where T <: Real -@test mccn isa Vector{T} where T <: Real -@test length(mr) == 5 +mc = multiscale(e, Composite(), x, est; maxscale) +mcn = multiscale_normalized(e, Composite(), x, est; maxscale) +@test mc isa Vector{T} where T <: Real +@test mcn isa Vector{T} where T <: Real @test length(mc) == 5 -@test length(mrn) == 5 @test length(mcn) == 5 # TODO: To verify that downsampling algorithms work as expected, we need to make a few From 683e304551f21139831b0c246cfe4bc3ef4a5cf6 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Thu, 3 Nov 2022 17:22:22 +0100 Subject: [PATCH 18/27] Correct analytical tests --- src/multiscale.jl | 18 +++++++++++---- src/multiscale/composite.jl | 39 +++++++++++++++++++++++---------- src/multiscale/regular.jl | 1 + test/multiscale/downsampling.jl | 36 ++++++++++++++++++++++++------ test/multiscale/multiscale.jl | 7 ++---- 5 files changed, 74 insertions(+), 27 deletions(-) diff --git a/src/multiscale.jl b/src/multiscale.jl index 5d7349a78..e7c50ac09 100644 --- a/src/multiscale.jl +++ b/src/multiscale.jl @@ -36,7 +36,7 @@ function multiscale end function multiscale_normalized end """ - multiscale(e::Entropy, alg, x, est; maxscale = 8, normalize = false) + multiscale(e::Entropy, alg, x, est; maxscale = length(x) ÷ 5) Compute the multi-scale entropy `e` with probabilities estimator `est` for timeseries `x`. @@ -65,7 +65,9 @@ methods are defined, ## Keyword Arguments -- `scalemax::Int`. The maximum number of scales (i.e. levels of downsampling). +- `scalemax::Int`. The maximum number of scales (i.e. levels of downsampling). The actual + maximum scale level is `length(x) ÷ 2`, but the default is `length(x) ÷ 5`, to avoid + computing the entropies for time series that are extremely short. - `normalize::Bool`. If `normalize == true`, then (if possible) compute normalized entropy/complexity. If `normalize == false`, then compute the non-normalized measure. @@ -80,7 +82,8 @@ end """ multiscale_normalized(e::Entropy, alg, x, est; maxscale = 8) -The same as [`multiscale`](@ref), but normalizes the estimated quantity according +The same as [`multiscale`](@ref), but normalizes the entropy if [`entropy_maximum`](@ref) +is implemented for `e`. """ function multiscale_normalized(e::Entropy, alg::MultiScaleAlgorithm, x, est::ProbabilitiesEstimator) @@ -88,7 +91,14 @@ function multiscale_normalized(e::Entropy, alg::MultiScaleAlgorithm, x, throw(ArgumentError(msg)) end - +max_scale_level(method::MultiScaleAlgorithm, x) = length(x) ÷ 2 +function verify_scale_level(method, x, s::Int) + err = DomainError( + "Maximum scale for length-$(length(x)) timeseries is "* + "`s = $(max_scale_level(method, x))`. Got s = $s" + ) + length(x) ÷ s >= 2 || throw(err) +end include("multiscale/regular.jl") diff --git a/src/multiscale/composite.jl b/src/multiscale/composite.jl index 6ebf95cd7..8c6fe4d27 100644 --- a/src/multiscale/composite.jl +++ b/src/multiscale/composite.jl @@ -48,23 +48,41 @@ Base.@kwdef struct Composite <: MultiScaleAlgorithm f::Function = Statistics.mean end -function downsample(method::Composite, x::AbstractVector{T}, s::Int, args...; kwargs...) where T +function downsample(method::Composite, x::AbstractVector{T}, s::Int, args...; + kwargs...) where T + verify_scale_level(method, x, s) + f = method.f - ET = eltype(one(1.0)) # consistently return floats, even if input is e.g. integer-valued + ET = eltype(one(1.0)) # always return floats, even if input is e.g. integer-valued if s == 1 return [ET.(x)] else N = length(x) - # note: there must be a typo or error in Wu et al. (2013), because if we use - # floor(N / s), as indicated in their paper, we'll get out of bounds errors. - # The number of samples needs to take into consideratio how many unique ways - # there are of selecting non-overlapping windows of length `s`. Hence, - # we use floor((N - s + 1) / s) instead. - L = floor(Int, (N - s + 1) / s) - ys = [zeros(ET, L) for i = 1:s] + # Because input time series are finite, there is always a minimum number of windows + # that we can construct at a given scale. We restrict the number of windows + # considered at each scale to this minimum to ensure windows are well-defined, + # i.e. we're not trying to summarize data at indices outside the input data, + # which would give out-of-bounds errors. + # + # For example, if the input is [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], then we sample + # the following subvectors at different scales: + # Scale 3: + # start index 1: [1, 2, 3], [4, 5, 6], [7, 8, 9] + # start index 2: [2, 3, 4], [5, 6, 7], [8, 9, 10] + # start index 3: [3, 4, 5], [6, 7, 8] + # Only two windows are possible for start index 3, so `min_possible_windows = 2` + # Scale 4: + # start index 1: [1, 2, 3, 4], [5, 6, 7, 8] + # start index 2: [2, 3, 4, 5], [6, 7, 8, 9] + # start index 3: [3, 4, 5, 6], [7, 8, 9, 10] + # start index 4: [4, 5, 6, 7] + # Only one window is possible for start index 4, so `min_possible_windows = 1` + min_possible_windows = floor(Int, (N - s + 1) / s) + + ys = [zeros(ET, min_possible_windows) for i = 1:s] for k = 1:s - for t = 1:L + for t = 1:min_possible_windows inds = ((t - 1)*s + k):(t * s + k - 1) ys[k][t] = @views f(x[inds], args...; kwargs...) end @@ -73,7 +91,6 @@ function downsample(method::Composite, x::AbstractVector{T}, s::Int, args...; kw end end -# TODO: make a separate multiscale_normalized? function multiscale(e::Entropy, alg::Composite, x::AbstractVector, est::ProbabilitiesEstimator; maxscale::Int = 10) diff --git a/src/multiscale/regular.jl b/src/multiscale/regular.jl index fbb03fa9d..fcff2ab6b 100644 --- a/src/multiscale/regular.jl +++ b/src/multiscale/regular.jl @@ -48,6 +48,7 @@ end function downsample(method::Regular, x::AbstractVector{T}, s::Int, args...; kwargs...) where T f = method.f + verify_scale_level(method, x, s) ET = eltype(one(1.0)) # consistently return floats, even if input is e.g. integer-valued if s == 1 diff --git a/test/multiscale/downsampling.jl b/test/multiscale/downsampling.jl index 82d242d4e..bbfb3af27 100644 --- a/test/multiscale/downsampling.jl +++ b/test/multiscale/downsampling.jl @@ -22,15 +22,37 @@ r = Regular(f = mean) # for scale 5, downsampling indices for `x` are 1:5, 6:10. Take the means: @test downsample(r, x, 5) == [3.0, 8.0] -# for scale 6, only a single set of sampling indices is possible - 1:6 -@test downsample(r, x, 6) == [3.5] -@test downsample(r, x, 6) == [3.5] - -# For scales larger than the number of elements, there should be no values in the -# downsampled time series. -@test isempty(downsample(r, x, length(x) + 1)) +@test_throws DomainError downsample(r, x, 6) ############################################################## # Regular downsampling ############################################################## c = Composite(f = mean) +# for scale 2, downsampling indices for `x` are 1:2, 3:4, 5:6, 7:8, 9:10, so we should get +# the means of consecutive running pairs +@test downsample(c, x, 2) == [ + [1.5, 3.5, 5.5, 7.5], # indices 1:2, 3:4, 5:6, 7:8 + [2.5, 4.5, 6.5, 8.5], # indices 2:3, 4:5, 6:7, 8:9 +] +@test downsample(c, x, 3) == [ + [2.0, 5.0], # indices 1:3, 4:6 + [3.0, 6.0], # indices 2:4, 5:7 + [4.0, 7.0], # indices 3:5, 6:8 +] + +@test downsample(c, x, 4) == [ + [2.5], # indices 1:4 + [3.5], # indices 2:5 + [4.5], # indices 3:6 + [5.5], # indices 4:7 +] + +@test downsample(c, x, 5) == [ + [3.0], # indices 1:5 + [4.0], # indices 2:6 + [5.0], # indices 3:7 + [6.0], # indices 4:8 + [7.0], # indices 5:9 +] + +@test_throws DomainError downsample(c, x, 6) diff --git a/test/multiscale/multiscale.jl b/test/multiscale/multiscale.jl index ae9a02446..342cb15a1 100644 --- a/test/multiscale/multiscale.jl +++ b/test/multiscale/multiscale.jl @@ -3,7 +3,8 @@ x = rand(100) maxscale = 5 e = Shannon() -# Generic tests. +# Generic tests is all we need here. The tests that make sure that the entropies are +# computed for the correctly sampled timeseries are in `/test/multiscale/downsampling.jl` mr = multiscale(e, Regular(), x, est; maxscale) mrn = multiscale_normalized(e, Regular(), x, est; maxscale) @test mr isa Vector{T} where T <: Real @@ -17,7 +18,3 @@ mcn = multiscale_normalized(e, Composite(), x, est; maxscale) @test mcn isa Vector{T} where T <: Real @test length(mc) == 5 @test length(mcn) == 5 - -# TODO: To verify that downsampling algorithms work as expected, we need to make a few -# simple analytical examples too. There are no such examples in the original papers, -# so we'll have to make some ourselves. From 2b98410b2681bb9e82646984dfb8cada1177c235 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Thu, 3 Nov 2022 17:32:34 +0100 Subject: [PATCH 19/27] Cleaner docs --- docs/src/multiscale.md | 38 ++++++++++---------------------------- src/multiscale.jl | 2 +- 2 files changed, 11 insertions(+), 29 deletions(-) diff --git a/docs/src/multiscale.md b/docs/src/multiscale.md index 9e238aceb..4ebff0ab1 100644 --- a/docs/src/multiscale.md +++ b/docs/src/multiscale.md @@ -9,12 +9,22 @@ defined below, and the following functions. - [`multiscale_normalized`](@ref) - [`downsample`](@ref) +## Multiscale entropy + ```@docs multiscale multiscale_normalized downsample ``` +### [Algorithms](@id multiscale_algorithms) + +```@docs +MultiScaleAlgorithm +Regular +Composite +``` + ## Available literature methods A non-exhaustive list of literature methods, and the syntax to compute them, are listed @@ -23,39 +33,11 @@ below. Please open an issue or make a pull-request to method missing from this list, or if you publish a paper based on some new multiscale combination. -### Multiscale complexity measures - -| Method | Syntax | Reference | -| ------------- | ------------- | ------------- | -| Multiscale sample entropy (first moment) | `multiscale(SampleEntropy(), Regular(f = mean), x)` |Costa et al. (2002)[^Costa2002] | -| Generalized multiscale sample entropy (second moment)| `multiscale(SampleEntropy(), Regular(f = std), x)` | Costa et al. (2015)[^Costa2015] | - -### Multiscale entropy - | Method | Syntax | Reference | | ------------- | ------------- | ------------- | | Refined composite multiscale dispersion entropy | `multiscale(Dispersion(), Composite(), x, normalized = true)` | Azami et al. (2017)[^Azami2017] | -## [Algorithms](@id multiscale_algorithms) - -```@docs -MultiScaleAlgorithm -Regular -Composite -``` - -[^Costa2002]: - Costa, M., Goldberger, A. L., & Peng, C. K. (2002). Multiscale entropy - analysis of complex physiologic time series. Physical review letters, 89(6), 068102. -[^Costa2015]: - Costa, M. D., & Goldberger, A. L. (2015). Generalized multiscale entropy - analysis: Application to quantifying the complex volatility of human heartbeat time - series. Entropy, 17(3), 1197-1203. [^Azami2017]: Azami, H., Rostaghi, M., Abásolo, D., & Escudero, J. (2017). Refined composite multiscale dispersion entropy and its application to biomedical signals. IEEE Transactions on Biomedical Engineering, 64(12), 2872-2879. -[^Morabito2012]: - Morabito, F. C., Labate, D., Foresta, F. L., Bramanti, A., Morabito, G., - & Palamara, I. (2012). Multivariate multi-scale permutation entropy for complexity - analysis of Alzheimer’s disease EEG. Entropy, 14(7), 1186-1202. diff --git a/src/multiscale.jl b/src/multiscale.jl index e7c50ac09..4395f4d07 100644 --- a/src/multiscale.jl +++ b/src/multiscale.jl @@ -36,7 +36,7 @@ function multiscale end function multiscale_normalized end """ - multiscale(e::Entropy, alg, x, est; maxscale = length(x) ÷ 5) + multiscale(e::Entropy, alg, x, est; kwargs...) Compute the multi-scale entropy `e` with probabilities estimator `est` for timeseries `x`. From 073b185b53c53fde2a581d23d3effb1791e679f4 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Thu, 3 Nov 2022 18:06:48 +0100 Subject: [PATCH 20/27] Switch order --- docs/src/multiscale.md | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/docs/src/multiscale.md b/docs/src/multiscale.md index 4ebff0ab1..fe09d9cb8 100644 --- a/docs/src/multiscale.md +++ b/docs/src/multiscale.md @@ -9,20 +9,20 @@ defined below, and the following functions. - [`multiscale_normalized`](@ref) - [`downsample`](@ref) -## Multiscale entropy +## [Algorithms](@id multiscale_algorithms) ```@docs -multiscale -multiscale_normalized -downsample +MultiScaleAlgorithm +Regular +Composite ``` -### [Algorithms](@id multiscale_algorithms) +## Multiscale entropy ```@docs -MultiScaleAlgorithm -Regular -Composite +multiscale +multiscale_normalized +downsample ``` ## Available literature methods From a53da06d005aedb4c6b0158a4ecb6b31dcde7256 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Fri, 4 Nov 2022 08:59:35 +0100 Subject: [PATCH 21/27] Shorten docs. --- docs/src/multiscale.md | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/docs/src/multiscale.md b/docs/src/multiscale.md index fe09d9cb8..1fb0708aa 100644 --- a/docs/src/multiscale.md +++ b/docs/src/multiscale.md @@ -2,14 +2,13 @@ ## Multiscale API -The multiscale API is defined by the [multiscale algorithms](@ref multiscale_algorithms) -defined below, and the following functions. +The multiscale API is defined by the functions - [`multiscale`](@ref) - [`multiscale_normalized`](@ref) - [`downsample`](@ref) -## [Algorithms](@id multiscale_algorithms) +which dispatch any of the [`MultiScaleAlgorithm`](@ref)s listed below. ```@docs MultiScaleAlgorithm From cafdcefd764b8dee2b1de1a8d8bce86beff77034 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Fri, 4 Nov 2022 11:19:59 +0100 Subject: [PATCH 22/27] Fix cross-references --- docs/src/probabilities.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/probabilities.md b/docs/src/probabilities.md index 51dd6af8d..41b3af086 100644 --- a/docs/src/probabilities.md +++ b/docs/src/probabilities.md @@ -4,9 +4,9 @@ The probabilities API is defined by -- [ProbabilitiesEstimator](@ref) -- [probabilities](@ref) -- [probabilities_and_outcomes](@ref) +- [`ProbabilitiesEstimator`](@ref) +- [`probabilities`](@ref) +- [`probabilities_and_outcomes`](@ref) ```@docs ProbabilitiesEstimator From e744a64ffbeae0c2f70340a6e3b201dc7f2979f2 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Fri, 16 Dec 2022 20:43:01 +0100 Subject: [PATCH 23/27] Auto stash before merge of "multiscale_entropy" and "main" --- src/multiscale.jl | 111 +++++++++++++++++++++++++++++++++ src/multiscale/composite.jl | 118 ++++++++++++++++++++++++++++++++++++ 2 files changed, 229 insertions(+) create mode 100644 src/multiscale.jl create mode 100644 src/multiscale/composite.jl diff --git a/src/multiscale.jl b/src/multiscale.jl new file mode 100644 index 000000000..8adbcf67a --- /dev/null +++ b/src/multiscale.jl @@ -0,0 +1,111 @@ +# This file contains an API for multiscale (coarse-grained/downsampled) computations. + +using Statistics +export multiscale +export multiscale_normalized +export downsample +export MultiScaleAlgorithm + +""" + MultiScaleAlgorithm + +The supertype for all multiscale algorithms. +""" +abstract type MultiScaleAlgorithm end + +""" + downsample(algorithm::MultiScaleAlgorithm, x, s::Int) + +Downsample and coarse-grain `x` to scale `s` according to the given multiscale `algorithm`. + +Positional arguments `args` and keyword arguments `kwargs` are propagated to relevant +functions in `algorithm`, if applicable. + +The return type depends on `algorithm`. For example: + +- [`Regular`](@ref) yields a single `Vector` per scale. +- [`Composite`](@ref) yields a `Vector{Vector}` per scale. +""" +downsample(method::MultiScaleAlgorithm, s::Int, x) + +downsample(alg::MultiScaleAlgorithm, s::Int, x::AbstractDataset) = + Dataset(map(t -> downsample(alg, s, t)), columns(x)...) + + +function multiscale end +function multiscale_normalized end + +""" + multiscale(e::Entropy, alg, est, x; maxscale::Int = 10) + +Compute the multi-scale entropy `e` with probabilities estimator `est` for timeseries `x`. + +## Description + +Utilizes [`downsample`](@ref) to compute the entropy/complexity of coarse-grained, +downsampled versions of `x` for scale factors `1:maxscale`. If `N = length(x)`, then the +length of the most severely downsampled version of `x` is `N ÷ scalemax`, while for scale +factor `1`, the original time series is considered. + +Provided the [`MultiScaleAlgorithm`](@ref), [`downsample`](@ref) +method and [`EntropyEstimator`](@ref) or [`ProbabilitiesEstimator`](@ref) is defined, +this function generalizes all multi-scale entropy estimators (discrete or continuous). + +Multi-scale complexity ("entropy-like") measures are found in the Complexity.jl package. + +## Arguments + +- `e::Entropy`. A valid [entropy type](@ref entropies), i.e. `Shannon()` or `Renyi()`. +- `alg::MultiScaleAlgorithm`. A valid [multiscale algorithm](@ref multiscale_algorithms), + i.e. `Regular()` or `Composite()`, which determines how down-sampling/coarse-graining + is performed. +- `x`. The input data. Usually a timeseries. +- `est`. For discrete entropy, `est` is a [`ProbabilitiesEstimator`](@ref), which determines + how probabilities are estimated for each sampled time series. Alternatively,for + continuous/differential entropy, `est` can be a [`EntropyEstimator`](@ref), + which dictates the entropy estimation method for each downsampled time series. + +## Keyword Arguments + +- `maxscale::Int`. The maximum number of scales (i.e. levels of downsampling). The actual + maximum scale level is `length(x) ÷ 2`, but the default is `length(x) ÷ 5`, to avoid + computing the entropies for time series that are extremely short. + +[^Costa2002]: Costa, M., Goldberger, A. L., & Peng, C. K. (2002). Multiscale entropy + analysis of complex physiologic time series. Physical review letters, 89(6), 068102. +""" +function multiscale(e::Entropy, alg::MultiScaleAlgorithm, + est::Union{ProbabilitiesEstimator, EntropyEstimator}, x) + msg = "`multiscale` entropy not implemented for $e $est on data type $(typeof(x))" + throw(ArgumentError(msg)) +end + +""" + multiscale_normalized(e::Entropy, alg, est, x) + +The same as [`multiscale`](@ref), but normalizes the entropy if [`entropy_maximum`](@ref) +is implemented for `e`. + +Only works if `est` is an [`EntropyEstimator`](@ref). +""" +function multiscale_normalized(e::Entropy, alg::MultiScaleAlgorithm, + est::ProbabilitiesEstimator, x) + msg = "`multiscale_normalized` not implemented for $e $est on data type $(typeof(x))" + throw(ArgumentError(msg)) +end + +multiscale_normalized(e::Entropy, alg::MultiScaleAlgorithm, est::EntropyEstimator, x) = + error("multiscale_normalized not defined for $(typeof(est))") + +max_scale_level(method::MultiScaleAlgorithm, x) = length(x) ÷ 2 +function verify_scale_level(method, s::Int, x) + err = DomainError( + "Maximum scale for length-$(length(x)) timeseries is "* + "`s = $(max_scale_level(method, x))`. Got s = $s" + ) + length(x) ÷ s >= 2 || throw(err) +end + + +include("multiscale/regular.jl") +include("multiscale/composite.jl") diff --git a/src/multiscale/composite.jl b/src/multiscale/composite.jl new file mode 100644 index 000000000..9eaaffa9f --- /dev/null +++ b/src/multiscale/composite.jl @@ -0,0 +1,118 @@ +export Composite + +""" + Composite <: MultiScaleAlgorithm + Composite(; f::Function = Statistics.mean) + +Composite multi-scale algorithm for multiscale entropy analysis (Wu et al., +2013)[^Wu2013], used, with [`multiscale`](@ref) to compute, for example, composite +multiscale entropy (CMSE). + +## Description + +Given a scalar-valued input time series `x`, the composite multiscale algorithm, +like [`Regular`](@ref), downsamples and coarse-grains `x` by splitting it into +non-overlapping windows of length `s`, and then constructing downsampled time series by +applying the function `f` to each of the resulting length-`s` windows. + +However, Wu et al. (2013)[^Wu2013] realized that for each scale `s`, there are actually `s` +different ways of selecting windows, depending on where indexing starts/ends. +These `s` different downsampled time series `D_t(s, f)` at each scale `s` are +constructed as follows: + +```math +\\{ D_{k}(s) \\} = \\{ D_{t, k}(s) \\}_{t = 1}^{L}, = \\{ f \\left( \\bf x_{t, k} \\right) \\} = +\\left\\{ + {f\\left( (x_i)_{i = (t - 1)s + k}^{ts + k - 1} \\right)} +\\right\\}_{t = 1}^{L}, +``` + +where `L = floor((N - s + 1) / s)` and `1 ≤ k ≤ s`, such that ``D_{i, k}(s)`` is the `i`-th +element of the `k`-th downsampled time series at scale `s`. + +Finally, compute ``\\dfrac{1}{s} \\sum_{k = 1}^s g(D_{k}(s))``, where `g` is some summary +function, for example [`entropy`](@ref) or [`complexity`](@ref). + +!!! note "Relation to Regular" + The downsampled time series ``D_{t, 1}(s)`` constructed using the composite + multiscale method is equivalent to the downsampled time series ``D_{t}(s)`` constructed + using the [`Regular`](@ref) method, for which `k == 1` is fixed, such that only + a single time series is returned. + +See also: [`Regular`](@ref). + +[^Wu2013]: Wu, S. D., Wu, C. W., Lin, S. G., Wang, C. C., & Lee, K. Y. (2013). Time series + analysis using composite multiscale entropy. Entropy, 15(3), 1069-1084. +""" +Base.@kwdef struct Composite <: MultiScaleAlgorithm + f::Function = Statistics.mean +end + +function downsample(method::Composite, s::Int, x::AbstractVector{T}, args...; + kwargs...) where T + verify_scale_level(method, x, s) + + f = method.f + ET = eltype(one(1.0)) # always return floats, even if input is e.g. integer-valued + + if s == 1 + return [ET.(x)] + else + N = length(x) + # Because input time series are finite, there is always a minimum number of windows + # that we can construct at a given scale. We restrict the number of windows + # considered at each scale to this minimum to ensure windows are well-defined, + # i.e. we're not trying to summarize data at indices outside the input data, + # which would give out-of-bounds errors. + # + # For example, if the input is [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], then we sample + # the following subvectors at different scales: + # Scale 3: + # start index 1: [1, 2, 3], [4, 5, 6], [7, 8, 9] + # start index 2: [2, 3, 4], [5, 6, 7], [8, 9, 10] + # start index 3: [3, 4, 5], [6, 7, 8] + # Only two windows are possible for start index 3, so `min_possible_windows = 2` + # Scale 4: + # start index 1: [1, 2, 3, 4], [5, 6, 7, 8] + # start index 2: [2, 3, 4, 5], [6, 7, 8, 9] + # start index 3: [3, 4, 5, 6], [7, 8, 9, 10] + # start index 4: [4, 5, 6, 7] + # Only one window is possible for start index 4, so `min_possible_windows = 1` + min_possible_windows = floor(Int, (N - s + 1) / s) + + ys = [zeros(ET, min_possible_windows) for i = 1:s] + for k = 1:s + for t = 1:min_possible_windows + inds = ((t - 1)*s + k):(t * s + k - 1) + ys[k][t] = @views f(x[inds], args...; kwargs...) + end + end + return ys + end +end + +function multiscale(e::Entropy, alg::Composite, est::ProbabilitiesEstimator, + x::AbstractVector; + maxscale::Int = 10) + + downscaled_timeseries = [downsample(alg, x, s) for s in 1:maxscale] + hs = zeros(Float64, maxscale) + for s in 1:maxscale + hs[s] = mean(entropy.(Ref(e), Ref(est), downscaled_timeseries[s])) + end + + return hs +end + +function multiscale_normalized(e::Entropy, alg::Composite, est::ProbabilitiesEstimator, + x::AbstractVector; + maxscale::Int = 10) + + downscaled_timeseries = [downsample(alg, x, s) for s in 1:maxscale] + hs = zeros(Float64, maxscale) + for s in 1:maxscale + hs[s] = mean(entropy_normalized.(Ref(e), Ref(est), downscaled_timeseries[s])) + end + + return hs +end From 6003f5e896c434315eec44ed6326fa1873339fda Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Fri, 16 Dec 2022 21:02:42 +0100 Subject: [PATCH 24/27] Re-introduce regular.jl file that was lost during merge --- src/multiscale/regular.jl | 83 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 83 insertions(+) create mode 100644 src/multiscale/regular.jl diff --git a/src/multiscale/regular.jl b/src/multiscale/regular.jl new file mode 100644 index 000000000..f05f47251 --- /dev/null +++ b/src/multiscale/regular.jl @@ -0,0 +1,83 @@ +export Regular + +""" + Regular <: MultiScaleAlgorithm + Regular(; f::Function = Statistics.mean) + +The original multi-scale algorithm for multiscale entropy analysis (Costa et al., +2022)[^Costa2002], which yields a single downsampled time series per scale `s`. + +## Description + +Given a scalar-valued input time series `x`, the `Regular` multiscale algorithm downsamples +and coarse-grains `x` by splitting it into non-overlapping windows of length `s`, and +then constructing a new downsampled time series ``D_t(s, f)`` by applying the function `f` +to each of the resulting length-`s` windows. + +The downsampled time series `D_t(s)` with `t ∈ [1, 2, …, L]`, where `L = floor(N / s)`, +is given by: + +```math +\\{ D_t(s, f) \\}_{t = 1}^{L} = \\left\\{ f \\left( \\bf x_t \\right) \\right\\}_{t = 1}^{L} = +\\left\\{ + {f\\left( (x_i)_{i = (t - 1)s + 1}^{ts} \\right)} +\\right\\}_{t = 1}^{L} +``` + +where `f` is some summary statistic applied to the length-`ts-((t - 1)s + 1)` tuples `xₖ`. +Different choices of `f` have yield different multiscale methods appearing in the +literature. For example: + +- `f == Statistics.mean` yields the original first-moment multiscale sample entropy (Costa + et al., 2002)[^Costa2002]. +- `f == Statistics.var` yields the generalized multiscale sample entropy (Costa & + Goldberger, 2015)[^Costa2015], which uses the second-moment (variance) instead of the + mean. + +See also: [`Composite`](@ref). + +[^Costa2002]: Costa, M., Goldberger, A. L., & Peng, C. K. (2002). Multiscale entropy + analysis of complex physiologic time series. Physical review letters, 89(6), 068102. +[^Costa2015]: Costa, M. D., & Goldberger, A. L. (2015). Generalized multiscale entropy + analysis: Application to quantifying the complex volatility of human heartbeat time + series. Entropy, 17(3), 1197-1203. +""" +Base.@kwdef struct Regular <: MultiScaleAlgorithm + f::Function = Statistics.mean +end + +function downsample(method::Regular, x::AbstractVector{T}, s::Int, args...; kwargs...) where T + f = method.f + verify_scale_level(method, s, x) + + ET = eltype(one(1.0)) # consistently return floats, even if input is e.g. integer-valued + if s == 1 + return x + else + N = length(x) + L = floor(Int, N / s) + ys = zeros(ET, L) + + for t = 1:L + inds = ((t - 1)*s + 1):(t * s) + ys[t] = @views f(x[inds], args...; kwargs...) + end + return ys + end +end + +function multiscale(e::Entropy, alg::Regular, x::AbstractVector, + est::ProbabilitiesEstimator; + maxscale::Int = 8) + + downscaled_timeseries = [downsample(alg, x, s) for s in 1:maxscale] + return entropy.(Ref(e), Ref(est), downscaled_timeseries) +end + +function multiscale_normalized(e::Entropy, alg::Regular, x::AbstractVector, + est::ProbabilitiesEstimator; + maxscale::Int = 8) + + downscaled_timeseries = [downsample(alg, x, s) for s in 1:maxscale] + return entropy_normalized.(Ref(e), Ref(est), downscaled_timeseries) +end From 8c717f0074850ba15e97c49387f4e16c4e8be410 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Fri, 16 Dec 2022 21:02:51 +0100 Subject: [PATCH 25/27] Correct argument order --- src/multiscale/composite.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/multiscale/composite.jl b/src/multiscale/composite.jl index 9eaaffa9f..bdccafd66 100644 --- a/src/multiscale/composite.jl +++ b/src/multiscale/composite.jl @@ -50,7 +50,7 @@ end function downsample(method::Composite, s::Int, x::AbstractVector{T}, args...; kwargs...) where T - verify_scale_level(method, x, s) + verify_scale_level(method, s, x) f = method.f ET = eltype(one(1.0)) # always return floats, even if input is e.g. integer-valued From 9f0be0e114b7e1e575601ece12125c8031bcdcb9 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Fri, 16 Dec 2022 21:30:48 +0100 Subject: [PATCH 26/27] Update multiscale and tests --- docs/make.jl | 1 + docs/src/multiscale.md | 42 ++++++++++++++++++++++++ src/Entropies.jl | 2 ++ src/multiscale.jl | 37 +++++++++++---------- src/multiscale/composite.jl | 13 ++++---- src/multiscale/regular.jl | 15 +++++---- test/multiscale/Composite.jl | 19 +++++++++++ test/multiscale/Regular.jl | 19 +++++++++++ test/multiscale/downsampling.jl | 58 +++++++++++++++++++++++++++++++++ test/multiscale/multiscale.jl | 3 ++ test/runtests.jl | 1 + 11 files changed, 180 insertions(+), 30 deletions(-) create mode 100644 docs/src/multiscale.md create mode 100644 test/multiscale/Composite.jl create mode 100644 test/multiscale/Regular.jl create mode 100644 test/multiscale/downsampling.jl create mode 100644 test/multiscale/multiscale.jl diff --git a/docs/make.jl b/docs/make.jl index bb7120f7a..ed0f4a495 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -36,6 +36,7 @@ ENTROPIES_PAGES = [ "index.md", "probabilities.md", "entropies.md", + "multiscale.md", "examples.md", "devdocs.md", ] diff --git a/docs/src/multiscale.md b/docs/src/multiscale.md new file mode 100644 index 000000000..8d3a69ea6 --- /dev/null +++ b/docs/src/multiscale.md @@ -0,0 +1,42 @@ +# Multiscale + +## Multiscale API + +The multiscale API is defined by the functions + +- [`multiscale`](@ref) +- [`multiscale_normalized`](@ref) +- [`downsample`](@ref) + +which dispatch any of the [`MultiScaleAlgorithm`](@ref)s listed below. + +```@docs +MultiScaleAlgorithm +Regular +Composite +``` + +## Multiscale entropy + +```@docs +multiscale +multiscale_normalized +downsample +``` + +## Available literature methods + +A non-exhaustive list of literature methods, and the syntax to compute them, are listed +below. Please open an issue or make a pull-request to +[Entropies.jl](https://github.com/JuliaDynamics/Entropies.jl) if you find a literature +method missing from this list, or if you publish a paper based on some new multiscale +combination. + +| Method | Syntax | Reference | +| ----------------------------------------------- | ------------------------------------------------------------------ | ------------------------------- | +| Refined composite multiscale dispersion entropy | `multiscale(Composite(), Dispersion(), est, x, normalized = true)` | Azami et al. (2017)[^Azami2017] | + +[^Azami2017]: + Azami, H., Rostaghi, M., Abásolo, D., & Escudero, J. (2017). Refined + composite multiscale dispersion entropy and its application to biomedical signals. + IEEE Transactions on Biomedical Engineering, 64(12), 2872-2879. \ No newline at end of file diff --git a/src/Entropies.jl b/src/Entropies.jl index 9d12ab4e3..a022330ae 100644 --- a/src/Entropies.jl +++ b/src/Entropies.jl @@ -20,6 +20,8 @@ const Vector_or_Dataset = Union{<:AbstractVector{<:Real}, <:AbstractDataset} include("probabilities.jl") include("entropy.jl") include("encodings.jl") +include("multiscale.jl") + # Library implementations (files include other files) include("probabilities_estimators/probabilities_estimators.jl") include("entropies/entropies.jl") diff --git a/src/multiscale.jl b/src/multiscale.jl index 8adbcf67a..076a16346 100644 --- a/src/multiscale.jl +++ b/src/multiscale.jl @@ -14,7 +14,7 @@ The supertype for all multiscale algorithms. abstract type MultiScaleAlgorithm end """ - downsample(algorithm::MultiScaleAlgorithm, x, s::Int) + downsample(algorithm::MultiScaleAlgorithm, s::Int, x) Downsample and coarse-grain `x` to scale `s` according to the given multiscale `algorithm`. @@ -36,23 +36,26 @@ function multiscale end function multiscale_normalized end """ - multiscale(e::Entropy, alg, est, x; maxscale::Int = 10) + multiscale(alg::MultiScaleAlgorithm, e::Entropy, est::EntropyEstimator, x; kwargs...) + multiscale(alg::MultiScaleAlgorithm, e::Entropy, est::ProbabilitiesEstimator, x; kwargs...) -Compute the multi-scale entropy `e` with probabilities estimator `est` for timeseries `x`. +Compute the multi-scale entropy `e` with estimator `est` for timeseries `x`. + +The first signature estimates differential/continuous multiscale entropy. The second +signature estimates discrete multiscale entropy. + +This function generalizes *all* multi-scale entropy estimators, as long as a relevant +[`MultiScaleAlgorithm`](@ref), [`downsample`](@ref) method and estimator is defined. +Multi-scale complexity ("entropy-like") measures, such as "sample entropy", are found in +the Complexity.jl package. ## Description Utilizes [`downsample`](@ref) to compute the entropy/complexity of coarse-grained, downsampled versions of `x` for scale factors `1:maxscale`. If `N = length(x)`, then the -length of the most severely downsampled version of `x` is `N ÷ scalemax`, while for scale +length of the most severely downsampled version of `x` is `N ÷ maxscale`, while for scale factor `1`, the original time series is considered. -Provided the [`MultiScaleAlgorithm`](@ref), [`downsample`](@ref) -method and [`EntropyEstimator`](@ref) or [`ProbabilitiesEstimator`](@ref) is defined, -this function generalizes all multi-scale entropy estimators (discrete or continuous). - -Multi-scale complexity ("entropy-like") measures are found in the Complexity.jl package. - ## Arguments - `e::Entropy`. A valid [entropy type](@ref entropies), i.e. `Shannon()` or `Renyi()`. @@ -74,27 +77,27 @@ Multi-scale complexity ("entropy-like") measures are found in the Complexity.jl [^Costa2002]: Costa, M., Goldberger, A. L., & Peng, C. K. (2002). Multiscale entropy analysis of complex physiologic time series. Physical review letters, 89(6), 068102. """ -function multiscale(e::Entropy, alg::MultiScaleAlgorithm, +function multiscale(alg::MultiScaleAlgorithm, e::Entropy, est::Union{ProbabilitiesEstimator, EntropyEstimator}, x) msg = "`multiscale` entropy not implemented for $e $est on data type $(typeof(x))" throw(ArgumentError(msg)) end """ - multiscale_normalized(e::Entropy, alg, est, x) + multiscale_normalized(alg::MultiScaleAlgorithm, e::Entropy, + est::ProbabilitiesEstimator, x) The same as [`multiscale`](@ref), but normalizes the entropy if [`entropy_maximum`](@ref) is implemented for `e`. -Only works if `est` is an [`EntropyEstimator`](@ref). +Note: this doesn't work if `est` is an [`EntropyEstimator`](@ref). """ -function multiscale_normalized(e::Entropy, alg::MultiScaleAlgorithm, - est::ProbabilitiesEstimator, x) - msg = "`multiscale_normalized` not implemented for $e $est on data type $(typeof(x))" +function multiscale_normalized(alg::MultiScaleAlgorithm, e::Entropy, est, x) + msg = "`multiscale_normalized` not implemented for $e $(typeof(est)) on data type $(typeof(x))" throw(ArgumentError(msg)) end -multiscale_normalized(e::Entropy, alg::MultiScaleAlgorithm, est::EntropyEstimator, x) = +multiscale_normalized(alg::MultiScaleAlgorithm, e::Entropy, est::EntropyEstimator, x) = error("multiscale_normalized not defined for $(typeof(est))") max_scale_level(method::MultiScaleAlgorithm, x) = length(x) ÷ 2 diff --git a/src/multiscale/composite.jl b/src/multiscale/composite.jl index bdccafd66..1907037a7 100644 --- a/src/multiscale/composite.jl +++ b/src/multiscale/composite.jl @@ -91,11 +91,12 @@ function downsample(method::Composite, s::Int, x::AbstractVector{T}, args...; end end -function multiscale(e::Entropy, alg::Composite, est::ProbabilitiesEstimator, +function multiscale(alg::Composite, e::Entropy, + est::Union{ProbabilitiesEstimator, EntropyEstimator}, x::AbstractVector; - maxscale::Int = 10) + maxscale::Int = 8) - downscaled_timeseries = [downsample(alg, x, s) for s in 1:maxscale] + downscaled_timeseries = [downsample(alg, s, x) for s in 1:maxscale] hs = zeros(Float64, maxscale) for s in 1:maxscale hs[s] = mean(entropy.(Ref(e), Ref(est), downscaled_timeseries[s])) @@ -104,11 +105,11 @@ function multiscale(e::Entropy, alg::Composite, est::ProbabilitiesEstimator, return hs end -function multiscale_normalized(e::Entropy, alg::Composite, est::ProbabilitiesEstimator, +function multiscale_normalized(alg::Composite, e::Entropy, est::ProbabilitiesEstimator, x::AbstractVector; - maxscale::Int = 10) + maxscale::Int = 8) - downscaled_timeseries = [downsample(alg, x, s) for s in 1:maxscale] + downscaled_timeseries = [downsample(alg, s, x) for s in 1:maxscale] hs = zeros(Float64, maxscale) for s in 1:maxscale hs[s] = mean(entropy_normalized.(Ref(e), Ref(est), downscaled_timeseries[s])) diff --git a/src/multiscale/regular.jl b/src/multiscale/regular.jl index f05f47251..3cbcb2c89 100644 --- a/src/multiscale/regular.jl +++ b/src/multiscale/regular.jl @@ -46,7 +46,7 @@ Base.@kwdef struct Regular <: MultiScaleAlgorithm f::Function = Statistics.mean end -function downsample(method::Regular, x::AbstractVector{T}, s::Int, args...; kwargs...) where T +function downsample(method::Regular, s::Int, x::AbstractVector{T}, args...; kwargs...) where T f = method.f verify_scale_level(method, s, x) @@ -66,18 +66,19 @@ function downsample(method::Regular, x::AbstractVector{T}, s::Int, args...; kwar end end -function multiscale(e::Entropy, alg::Regular, x::AbstractVector, - est::ProbabilitiesEstimator; +function multiscale(alg::Regular, e::Entropy, + est::Union{ProbabilitiesEstimator, EntropyEstimator}, + x::AbstractVector; maxscale::Int = 8) - downscaled_timeseries = [downsample(alg, x, s) for s in 1:maxscale] + downscaled_timeseries = [downsample(alg, s, x) for s in 1:maxscale] return entropy.(Ref(e), Ref(est), downscaled_timeseries) end -function multiscale_normalized(e::Entropy, alg::Regular, x::AbstractVector, - est::ProbabilitiesEstimator; +function multiscale_normalized(alg::Regular, e::Entropy, + est::ProbabilitiesEstimator, x::AbstractVector,; maxscale::Int = 8) - downscaled_timeseries = [downsample(alg, x, s) for s in 1:maxscale] + downscaled_timeseries = [downsample(alg, s, x) for s in 1:maxscale] return entropy_normalized.(Ref(e), Ref(est), downscaled_timeseries) end diff --git a/test/multiscale/Composite.jl b/test/multiscale/Composite.jl new file mode 100644 index 000000000..b4ce894d3 --- /dev/null +++ b/test/multiscale/Composite.jl @@ -0,0 +1,19 @@ +using Entropies + +est = Dispersion() +x = rand(100) +maxscale = 5 +e = Shannon() + +# Generic tests is all we need here. The tests that make sure that the entropies are +# computed for the correctly sampled timeseries are in `/test/multiscale/downsampling.jl` +mc = multiscale(Composite(), e, est, x; maxscale) +mcn = multiscale_normalized(Composite(), e, est, x; maxscale) +@test mc isa Vector{T} where T <: Real +@test mcn isa Vector{T} where T <: Real +@test length(mc) == 5 +@test length(mcn) == 5 + +# `EntropyEstimator`s` should work for `multiscale`, but not `multiscale_normalized` +@test multiscale(Composite(), e, Kraskov(), x) isa Vector{T} where T <: Real +@test_throws ErrorException multiscale_normalized(Composite(), e, Kraskov(), x) diff --git a/test/multiscale/Regular.jl b/test/multiscale/Regular.jl new file mode 100644 index 000000000..684d47213 --- /dev/null +++ b/test/multiscale/Regular.jl @@ -0,0 +1,19 @@ +using Entropies + +est = Dispersion() +x = rand(100) +maxscale = 5 +e = Shannon() + +# Generic tests is all we need here. The tests that make sure that the entropies are +# computed for the correctly sampled timeseries are in `/test/multiscale/downsampling.jl` +mr = multiscale(Regular(), e, est, x; maxscale) +mrn = multiscale_normalized(Regular(), e, est, x; maxscale) +@test mr isa Vector{T} where T <: Real +@test mrn isa Vector{T} where T <: Real +@test length(mr) == 5 +@test length(mrn) == 5 + +# `EntropyEstimator`s` should work for `multiscale`, but not `multiscale_normalized` +@test multiscale(Regular(), e, Kraskov(), x) isa Vector{T} where T <: Real +@test_throws ErrorException multiscale_normalized(Regular(), e, Kraskov(), x) diff --git a/test/multiscale/downsampling.jl b/test/multiscale/downsampling.jl new file mode 100644 index 000000000..b1817efd3 --- /dev/null +++ b/test/multiscale/downsampling.jl @@ -0,0 +1,58 @@ +using Statistics +x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] + +############################################################## +# Regular downsampling +############################################################## +r = Regular(f = mean) + +# just returns original time series +@test downsample(r, 1, x) == x + +# for scale 2, downsampling indices for `x` are 1:2, 3:4, 5:6, 7:8, 9:10, so we should get +# the means of consecutive running pairs +@test downsample(r, 2, x) == [1.5, 3.5, 5.5, 7.5, 9.5] + +# for scale 3, downsampling indices for `x` are 1:3, 4:6, 7:9. Take the means: +@test downsample(r, 3, x) == [2.0, 5.0, 8.0] + +# for scale 4, downsampling indices for `x` are 1:4, 5:8. Take the means: +@test downsample(r, 4, x) == [2.5, 6.5] + +# for scale 5, downsampling indices for `x` are 1:5, 6:10. Take the means: +@test downsample(r, 5, x) == [3.0, 8.0] + +@test_throws DomainError downsample(r, 6, x) + +############################################################## +# Regular downsampling +############################################################## +c = Composite(f = mean) +# for scale 2, downsampling indices for `x` are 1:2, 3:4, 5:6, 7:8, 9:10, so we should get +# the means of consecutive running pairs +@test downsample(c, 2, x) == [ + [1.5, 3.5, 5.5, 7.5], # indices 1:2, 3:4, 5:6, 7:8 + [2.5, 4.5, 6.5, 8.5], # indices 2:3, 4:5, 6:7, 8:9 +] +@test downsample(c, 3, x) == [ + [2.0, 5.0], # indices 1:3, 4:6 + [3.0, 6.0], # indices 2:4, 5:7 + [4.0, 7.0], # indices 3:5, 6:8 +] + +@test downsample(c, 4, x) == [ + [2.5], # indices 1:4 + [3.5], # indices 2:5 + [4.5], # indices 3:6 + [5.5], # indices 4:7 +] + +@test downsample(c, 5, x) == [ + [3.0], # indices 1:5 + [4.0], # indices 2:6 + [5.0], # indices 3:7 + [6.0], # indices 4:8 + [7.0], # indices 5:9 +] + +@test_throws DomainError downsample(c, 6, x) diff --git a/test/multiscale/multiscale.jl b/test/multiscale/multiscale.jl new file mode 100644 index 000000000..f83a41ac4 --- /dev/null +++ b/test/multiscale/multiscale.jl @@ -0,0 +1,3 @@ +testfile("downsampling.jl") +testfile("Regular.jl") +testfile("Composite.jl") diff --git a/test/runtests.jl b/test/runtests.jl index 153abe8e6..3332a9967 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,6 +6,7 @@ testfile(file, testname=defaultname(file)) = @testset "$testname" begin; include @testset "Entropies.jl" begin include("probabilities/probabilities.jl") include("entropies/entropies.jl") + testfile("multiscale/multiscale.jl") # Various testfile("utils/fasthist.jl") From 54c1c3eefa47489e5ed89047d0c86e6f85a5d0f2 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Fri, 16 Dec 2022 22:19:56 +0100 Subject: [PATCH 27/27] Update gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index cc8db9b16..fcd1cceec 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,4 @@ docs/build/* Manifest.toml *.scss *.css +.DS_Store \ No newline at end of file