Repository: JuliaDynamics/ComplexityMeasures.jl (julia)

Issue:
Implement `outcome_space` and adjust what is returned by `outcomes` for `SpatialSymbolicPermutation`
#172 does the implementation of `outcome_space`, as well as actually returning the permutations themselves (not integers), for the 1D permutation based probability estimators. This is a reminder to do the same for the spatial estimators once #126 is merged.

EDIT: also test `outcome_space` and `total_outcomes`.

Pull request:
Fix #173

Patch:
diff --git a/src/probabilities.jl b/src/probabilities.jl
index bd6034d48..37cd93744 100644
--- a/src/probabilities.jl
+++ b/src/probabilities.jl
@@ -78,6 +78,7 @@ All currently implemented probability estimators are:
 - [`ValueHistogram`](@ref).
 - [`TransferOperator`](@ref).
 - [`Dispersion`](@ref).
+- [`SpatialDispersion`](@ref).
 - [`WaveletOverlap`](@ref).
 - [`PowerSpectrum`](@ref).
 - [`SymbolicPermutation`](@ref).
diff --git a/src/probabilities_estimators/dispersion/dispersion.jl b/src/probabilities_estimators/dispersion/dispersion.jl
index a901e1a52..64dd80d04 100644
--- a/src/probabilities_estimators/dispersion/dispersion.jl
+++ b/src/probabilities_estimators/dispersion/dispersion.jl
@@ -49,10 +49,10 @@ always get non-zero probabilities using the `Dispersion` probability estimator.
 
 ## Outcome space
 
-The outcome space for `Dispersion` is the unique delay vectors with elements the
+The outcome space for `Dispersion` is the unique delay vectors whose elements are the
 the symbols (integers) encoded by the Gaussian CDF.
 Hence, the outcome space is all `m`-dimensional delay vectors whose elements
-are all possible values in `1:c`.
+are all possible values in `1:c`. There are `c ^ m` such vectors.
 
 ## Data requirements and parameters
 
diff --git a/src/probabilities_estimators/spatial/spatial_dispersion/SpatialDispersion.jl b/src/probabilities_estimators/spatial/spatial_dispersion/SpatialDispersion.jl
index 300f70028..8f0f3ee56 100644
--- a/src/probabilities_estimators/spatial/spatial_dispersion/SpatialDispersion.jl
+++ b/src/probabilities_estimators/spatial/spatial_dispersion/SpatialDispersion.jl
@@ -11,16 +11,26 @@ import Base.maximum
         L = nothing,
     )
 
-A dispersion-based probabilities/entropy estimator for `N`-dimensional spatiotemporal
-systems, based on Azami et al. (2019)'s 2D square dispersion entropy estimator,
-but here generalized for `N`-dimensional input data `x`.
+A dispersion-based probabilities estimator that generalises [`Dispersion`](@ref) for
+input data that are high-dimensional arrays.
+
+`SpatialDispersion` is based on Azami et al. (2019)[^Azami2019]'s 2D square dispersion
+(Shannon) entropy estimator, but is here implemented as a pure probabilities
+probabilities estimator that is generalized for `N`-dimensional input data `x`,
+with arbitrary neighborhood regions (stencils) and (optionally) periodic boundary
+conditions.
+
+In combination with [`entropy`](@ref) and [`entropy_normalized`](@ref), this probabilities
+estimator can be used to compute (normalized) generalized spatiotemporal dispersion
+[`Entropy`](@ref) of any type.
 
 ## Arguments
 
 - `stencil`. Defines what local area (hyperrectangle), or which points within this area,
     to include around each hypervoxel (i.e. pixel in 2D). The examples below demonstrate
     different ways of specifying stencils. For details, see
-    [`SpatialSymbolicPermutation`](@ref).
+    [`SpatialSymbolicPermutation`](@ref). See [`SpatialSymbolicPermutation`](@ref) for
+    more information about stencils.
 -  `x::AbstractArray`. The input data. Must be provided because we need to know its size
     for optimization and bound checking.
 
@@ -29,8 +39,7 @@ but here generalized for `N`-dimensional input data `x`.
 - `periodic::Bool`. If `periodic == true`, then the stencil should wrap around at the
     end of the array. If `periodic = false`, then pixels whose stencil exceeds the array
     bounds are skipped.
-- `encoding::Encoding`. Determines how input data is mapped to discrete categories. Must be
-    a valid [`Encoding`](@ref).
+- `c::Int`. Determines how many discrete categories to use for the Gaussian encoding.
 - `skip_encoding`. If `skip_encoding == true`, `encoding` is ignored, and dispersion
     patterns are computed directly from `x`, under the assumption that `L` is the alphabet
     length for `x` (useful for categorical or integer data). Thus, if
@@ -40,6 +49,13 @@ but here generalized for `N`-dimensional input data `x`.
     `stencil` and `encoding`. If `L` is set to an integer, then the data is considered
     pre-encoded and the number of total outcomes is set to `L`.
 
+## Outcome space
+
+The outcome space for `SpatialDispersion` is the unique delay vectors whose elements are the
+the symbols (integers) encoded by the Gaussian CDF. Hence, the outcome space is all
+`m`-dimensional delay vectors whose elements are all possible values in `1:c`.
+There are `c ^ m` such vectors.
+
 ## Description
 
 Estimating probabilities/entropies from higher-dimensional data is conceptually simple.
@@ -47,8 +63,8 @@ Estimating probabilities/entropies from higher-dimensional data is conceptually
 1. Discretize each value (hypervoxel) in `x` relative to all other values `xᵢ ∈ x` using the
     provided `encoding` scheme.
 2. Use `stencil` to extract relevant (discretized) points around each hypervoxel.
-3. Construct a symbol string from these points.
-4. Take the sum-normalized histogram of the symbol strings as a probability distribution.
+3. Construct a symbol these points.
+4. Take the sum-normalized histogram of the symbol as a probability distribution.
 5. Optionally, compute [`entropy`](@ref) or [`entropy_normalized`](@ref) from this
     probability distribution.
 
@@ -85,6 +101,16 @@ est = SpatialDispersion(stencil, first(imgs))
 h_vs_t = entropy_normalized.(Ref(est), imgs)
 ```
 
+
+Computing generalized spatiotemporal dispersion entropy is trivial, e.g. with
+[`Renyi`](@ref):
+
+```julia
+x = reshape(repeat(1:5, 500) .+ 0.1*rand(500*5), 50, 50)
+est = SpatialDispersion(stencil, x)
+entropy(Renyi(q = 2), est, x)
+```
+
 See also: [`SpatialSymbolicPermutation`](@ref), [`GaussianCDFEncoding`](@ref),
 [`symbolize`](@ref).
 
@@ -183,7 +209,4 @@ function total_outcomes(est::SpatialDispersion)::Int
     end
 end
 
-# TODO: how to represent the outcomes? We have to think about this...
-function outcome_space(est::SpatialDispersion)
-
-end
+outcome_space(est::SpatialDispersion) = outcome_space(Dispersion(; c = est.c, m = est.m))
diff --git a/src/probabilities_estimators/spatial/spatial_permutation/SpatialSymbolicPermutation.jl b/src/probabilities_estimators/spatial/spatial_permutation/SpatialSymbolicPermutation.jl
index c4531ab63..fad309527 100644
--- a/src/probabilities_estimators/spatial/spatial_permutation/SpatialSymbolicPermutation.jl
+++ b/src/probabilities_estimators/spatial/spatial_permutation/SpatialSymbolicPermutation.jl
@@ -6,36 +6,44 @@ export SpatialSymbolicPermutation
 ###########################################################################################
 
 """
-    SpatialSymbolicPermutation(stencil, x, periodic = true)
+    SpatialSymbolicPermutation <: ProbabilitiesEstimator
+    SpatialSymbolicPermutation(stencil, x; periodic = true)
 
-A symbolic, permutation-based probabilities estimator for spatiotemporal systems.
+A symbolic, permutation-based probabilities estimator for spatiotemporal systems that
+generalises [`SymbolicPermutation`](@ref) to high-dimensional arrays.
 
-The input data `x` are high-dimensional arrays, for example 2D arrays [^Ribeiro2012] or 3D arrays
-[^Schlemmer2018]. This approach is also known as _spatiotemporal permutation entropy_.
-`x` is given because we need to know its size for optimization and bound checking.
+`SpatialSymbolicPermutation` is based on the 2D and 3D *spatiotemporal permutation entropy*
+estimators by by Ribeiro et al. (2012)[^Ribeiro2012] and Schlemmer et al.
+(2018)[^Schlemmer2018]), respectively, but is here implemented as a pure probabilities
+probabilities estimator that is generalized for `N`-dimensional input data `x`,
+with arbitrary neighborhood regions (stencils) and periodic boundary conditions.
 
-A _stencil_ defines what local area (which points) around each pixel to
-consider, and compute ordinal patterns from.
-The number of included points in a stencil (`m`) determines the length of the vectors
-to be discretized, i.e. there are `m!` possible ordinal patterns around each pixel.
+In combination with [`entropy`](@ref) and [`entropy_normalized`](@ref), this
+probabilities estimator can be used to compute (normalized) generalized spatiotemporal
+permutation [`Entropy`](@ref) of any type.
 
-Example usage:
-```julia
-data = [rand(50, 50) for _ in 1:50]
-x = data[1] # first "time slice" of a spatial system evolution
-stencil = ...
-est = SpatialSymbolicPermutation(stencil, x)
-```
+## Arguments
+
+- `stencil`. Defines what local area (hyperrectangle), or which points within this area,
+    to include around each hypervoxel (i.e. pixel in 2D). The examples below demonstrate
+    different ways of specifying stencils. For details, see
+    [`SpatialSymbolicPermutation`](@ref).
+-  `x::AbstractArray`. The input data. Must be provided because we need to know its size
+    for optimization and bound checking.
+
+## Keyword arguments
+
+- `periodic::Bool`. If `periodic == true`, then the stencil should wrap around at the
+    end of the array. If `periodic = false`, then pixels whose stencil exceeds the array
+    bounds are skipped.
+
+## Stencils
 
 Stencils are passed in one of the following three ways:
 
 1. As vectors of `CartesianIndex` which encode the pixels to include in the
     stencil, with respect to the current pixel, or integer arrays of the same dimensionality
-    as the data. For example
-
-    ```julia
-    stencil = CartesianIndex.([(0,0), (0,1), (1,1), (1,0)])
-    ```
+    as the data. For example `stencil = CartesianIndex.([(0,0), (0,1), (1,1), (1,0)])`.
     Don't forget to include the zero offset index if you want to include the point itself,
     which is almost always the case.
     Here the stencil creates a 2x2 square extending to the bottom and right of the pixel
@@ -45,11 +53,7 @@ Stencils are passed in one of the following three ways:
 2. As a `D`-dimensional array (where `D` matches the dimensionality of the input data)
     containing `0`s and `1`s, where if `stencil[index] == 1`, the corresponding pixel is
     included, and if `stencil[index] == 0`, it is not included.
-    To generate the same estimator as in 1., use
-
-    ```julia
-    stencil = [1 1; 1 1]
-    ```
+    To generate the same estimator as in 1., use `stencil = [1 1; 1 1]`.
     When passing a stencil as a `D`-dimensional array, `m = sum(stencil)`
 
 3. As a `Tuple` containing two `Tuple`s, both of length `D`, for `D`-dimensional data.
@@ -57,24 +61,46 @@ Stencils are passed in one of the following three ways:
     dictates the number of pixels to be included along the `i`th axis and `lag[i]`
     the separation of pixels along the same axis.
     This method can only generate (hyper)rectangular stencils. To create the same estimator as
-    in the previous examples, use here
-
-    ```julia
-    stencil = ((2, 2), (1, 1))
-    ```
+    in the previous examples, use here `stencil = ((2, 2), (1, 1))`.
     When passing a stencil using `extent` and `lag`, `m = prod(extent)!`.
 
-After having defined `est`, one calculates the spatial permutation entropy
-by calling [`entropy`](@ref) with `est`, and with the array data.
+## Example: spatiotemporal entropy for time series
+
+Usage is simple. First, define a `SpatialSymbolicPermutation` estimator by specifying
+a stencil and giving some input data (a matrix with the same dimensions as the data
+as you're going to analyse). Then simply call [`entropy`](@ref) with the estimator.
+
+```julia
+using Entropies
+x = rand(50, 50) # first "time slice" of a spatial system evolution
+stencil = [1 1; 0 1] # or one of the other ways of specifying stencils
+est = SpatialSymbolicPermutation(stencil, x)
+h = entropy(est, x)
+```
+
 To apply this to timeseries of spatial data, simply loop over the call, e.g.:
 
 ```julia
-h = entropy(x, est)
-h_vs_t = [entropy(d, est) for d in data]
+data = [rand(50, 50) for i in 1:50]
+est = SpatialSymbolicPermutation(stencil, first(data))
+h_vs_t = [entropy(est, d) for d in data]
+```
+
+Computing generalized spatiotemporal permutation entropy is trivial, e.g. with
+[`Renyi`](@ref):
+
+```julia
+x = reshape(repeat(1:5, 500) .+ 0.1*rand(500*5), 50, 50)
+est = SpatialSymbolicPermutation(stencil, x)
+entropy(Renyi(q = 2), est, x)
 ```
 
-The argument `periodic` decides whether the stencil should wrap around at the end of the
-array. If `periodic = false`, pixels whose stencil exceeds the array bounds are skipped.
+## Outcome space
+
+The outcome space `Ω` for `SpatialSymbolicPermutation` is the set of length-`m` ordinal
+patterns (i.e. permutations) that can be formed by the integers `1, 2, …, m`,
+ordered lexicographically. There are `factorial(m)` such patterns.
+Here, `m` refers to the number of points included by `stencil`.
 
 [^Ribeiro2012]:
     Ribeiro et al. (2012). Complexity-entropy causality plane as a complexity measure
@@ -103,7 +129,6 @@ function SpatialSymbolicPermutation(stencil, x::AbstractArray{T, D};
     )
 end
 
-
 function probabilities(est::SpatialSymbolicPermutation, x)
     # TODO: This can be literally a call to `symbolize` and then
     # calling probabilities on it. Should do once the `symbolize` refactoring is done.
@@ -139,6 +164,11 @@ function Base.show(io::IO, est::SpatialSymbolicPermutation{D}) where {D}
     show(io, MIME"text/plain"(), est.stencil)
 end
 
+function outcome_space(est::SpatialSymbolicPermutation)
+    encoding = OrdinalPatternEncoding(; est.m, est.lt)
+    decode.(Ref(encoding), 1:factorial(est.m))
+end
+
 function total_outcomes(est::SpatialSymbolicPermutation)
     return factorial(est.m)
 end


Test patch:
diff --git a/test/probabilities/estimators/spatial/spatial_dispersion.jl b/test/probabilities/estimators/spatial/spatial_dispersion.jl
index f575372ac..b3a8212b1 100644
--- a/test/probabilities/estimators/spatial/spatial_dispersion.jl
+++ b/test/probabilities/estimators/spatial/spatial_dispersion.jl
@@ -53,3 +53,5 @@ est_y_presymb = SpatialDispersion(stencil, y; skip_encoding =  true, L = 2)
 
 @test 0.99 <= round(entropy_normalized(est_y, y), digits = 2) <= 1.0
 @test 0.99 <= round(entropy_normalized(est_y_presymb, y), digits = 2) <= 1.0
+
+@test outcome_space(est_y) == outcome_space(Dispersion(c = est_y.c, m = est_y.m))
diff --git a/test/probabilities/estimators/spatial/spatial_permutation.jl b/test/probabilities/estimators/spatial/spatial_permutation.jl
index ca4c0e7e2..fe731a4f6 100644
--- a/test/probabilities/estimators/spatial/spatial_permutation.jl
+++ b/test/probabilities/estimators/spatial/spatial_permutation.jl
@@ -80,3 +80,5 @@ stencil = [1 1; 1 1];
 est = SpatialSymbolicPermutation(stencil, x)
 hsp = entropy_normalized(Renyi(), est, x)
 @test round(hsp, digits = 2) == 1.00
+
+@test outcome_space(est) == outcome_space(SymbolicPermutation(m = est.m))