Multiscale analysis

.gitignore

@@ -2,3 +2,4 @@ docs/build/*
 Manifest.toml
 *.scss
 *.css
+.DS_Store

docs/make.jl

@@ -36,6 +36,7 @@ ENTROPIES_PAGES = [
     "index.md",
     "probabilities.md",
     "entropies.md",
+    "multiscale.md",
     "examples.md",
     "devdocs.md",
 ]

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.

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")

src/multiscale.jl

@@ -0,0 +1,114 @@
+# 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, s::Int, x)
+
+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(alg::MultiScaleAlgorithm, e::Entropy, est::EntropyEstimator, x; kwargs...)
+    multiscale(alg::MultiScaleAlgorithm, e::Entropy, est::ProbabilitiesEstimator, x; kwargs...)
+
+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 ÷ maxscale`, while for scale
+factor `1`, the original time series is considered.
+
+## 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(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(alg::MultiScaleAlgorithm, e::Entropy,
+        est::ProbabilitiesEstimator, x)
+
+The same as [`multiscale`](@ref), but normalizes the entropy if [`entropy_maximum`](@ref)
+is implemented for `e`.
+
+Note: this doesn't work if `est` is an [`EntropyEstimator`](@ref).
+"""
+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(alg::MultiScaleAlgorithm, e::Entropy, 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")

src/multiscale/composite.jl

@@ -0,0 +1,119 @@
+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, s, x)
+
+    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(alg::Composite, e::Entropy,
+        est::Union{ProbabilitiesEstimator, EntropyEstimator},
+        x::AbstractVector;
+        maxscale::Int = 8)
+
+    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]))
+    end
+
+    return hs
+end
+
+function multiscale_normalized(alg::Composite, e::Entropy, est::ProbabilitiesEstimator,
+        x::AbstractVector;
+        maxscale::Int = 8)
+
+    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]))
+    end
+
+    return hs
+end

src/multiscale/regular.jl

@@ -0,0 +1,84 @@
+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, s::Int, x::AbstractVector{T}, 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(alg::Regular, e::Entropy,
+        est::Union{ProbabilitiesEstimator, EntropyEstimator},
+        x::AbstractVector;
+        maxscale::Int = 8)
+
+    downscaled_timeseries = [downsample(alg, s, x) for s in 1:maxscale]
+    return entropy.(Ref(e), Ref(est), downscaled_timeseries)
+end
+
+function multiscale_normalized(alg::Regular, e::Entropy,
+        est::ProbabilitiesEstimator, x::AbstractVector,;
+        maxscale::Int = 8)
+
+    downscaled_timeseries = [downsample(alg, s, x) for s in 1:maxscale]
+    return entropy_normalized.(Ref(e), Ref(est), downscaled_timeseries)
+end