From 810142b5a9fdb98d2165bf507bf8f953e13de263 Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Sun, 4 Sep 2022 08:43:36 +0100
Subject: [PATCH 01/11] Add docstring of spatiotemporal permutation

---
 src/core.jl                       | 14 +++++-----
 src/symbolic/symbolic_stencils.jl | 46 +++++++++++++++++++++++++++++++
 2 files changed, 53 insertions(+), 7 deletions(-)
 create mode 100644 src/symbolic/symbolic_stencils.jl

diff --git a/src/core.jl b/src/core.jl
index 6c9850ae3..9fe79a638 100644
--- a/src/core.jl
+++ b/src/core.jl
@@ -6,7 +6,7 @@ export probabilities, probabilities!
 export genentropy, genentropy!
 export Dataset, dimension
 
-const Vector_or_Dataset = Union{AbstractVector, AbstractDataset}
+const Array_or_Dataset = Union{AbstractArray, AbstractDataset}
 
 """
     Probabilities(x) → p
@@ -49,7 +49,7 @@ abstract type ProbabilitiesEstimator end
 const ProbEst = ProbabilitiesEstimator # shorthand
 
 """
-    probabilities(x::Vector_or_Dataset, est::ProbabilitiesEstimator) → p::Probabilities
+    probabilities(x::Array_or_Dataset, est::ProbabilitiesEstimator) → p::Probabilities
 
 Calculate probabilities representing `x` based on the provided
 estimator and return them as a [`Probabilities`](@ref) container (`Vector`-like).
@@ -58,7 +58,7 @@ documentation of the individual estimators for more.
 
 The configuration options are always given as arguments to the chosen estimator.
 
-    probabilities(x::Vector_or_Dataset, ε::AbstractFloat) → p::Probabilities
+    probabilities(x::Array_or_Dataset, ε::AbstractFloat) → p::Probabilities
 
 Convenience syntax which provides probabilities for `x` based on rectangular binning
 (i.e. performing a histogram). In short, the state space is divided into boxes of length
@@ -71,11 +71,11 @@ This allows computation of probabilities (histograms) of high-dimensional
 datasets and with small box sizes `ε` without memory overflow and with maximum performance.
 To obtain the bin information along with `p`, use [`binhist`](@ref).
 
-    probabilities(x::Vector_or_Dataset, n::Integer) → p::Probabilities
+    probabilities(x::Array_or_Dataset, n::Integer) → p::Probabilities
 Same as the above method, but now each dimension of the data is binned into `n::Int` equal
 sized bins instead of bins of length `ε::AbstractFloat`.
 
-    probabilities(x::Vector_or_Dataset) → p::Probabilities
+    probabilities(x::Array_or_Dataset) → p::Probabilities
 Directly count probabilities from the elements of `x` without any discretization,
 binning, or other processing (mostly useful when `x` contains categorical or integer data).
 """
@@ -97,7 +97,7 @@ Compute the generalized order-`q` entropy of some probabilities
 returned by the [`probabilities`](@ref) function. Alternatively, compute entropy
 from pre-computed `Probabilities`.
 
-    genentropy(x::Vector_or_Dataset, est; q = 1.0, base)
+    genentropy(x::Array_or_Dataset, est; q = 1.0, base)
 
 A convenience syntax, which calls first `probabilities(x, est)`
 and then calculates the entropy of the result (and thus `est` can be a
@@ -152,7 +152,7 @@ end
 genentropy(x::AbstractArray{<:Real}) =
     error("For single-argument input, do `genentropy(Probabilities(x))` instead.")
 
-function genentropy(x::Vector_or_Dataset, est; q = 1.0, α = nothing, base = MathConstants.e)
+function genentropy(x::Array_or_Dataset, est; q = 1.0, α = nothing, base = MathConstants.e)
     if α ≠ nothing
         @warn "Keyword `α` is deprecated in favor of `q`."
         q = α
diff --git a/src/symbolic/symbolic_stencils.jl b/src/symbolic/symbolic_stencils.jl
new file mode 100644
index 000000000..a43ac2254
--- /dev/null
+++ b/src/symbolic/symbolic_stencils.jl
@@ -0,0 +1,46 @@
+"""
+    SpatiotemporalPermutation(stencil, periodic = true)
+A symbolic, permutation-based probabilities/entropy estimator for higher-dimensional
+data, typically representing spatiotemporal systems. The data are higher-dimensional
+arrays, such as 2D [^Ribeiro2012] or 3D [^Schlemmer2018] `Matrix`.
+This approach is also known as _Spatiotemporal Permutation Entropy_.
+
+A _stencil_ defines what local area around each pixel to
+consider, and compute the ordinal pattern within the stencil. Stencils are given as
+vectors of `CartesianIndex` which encode the _offsets_ of the pixes to include in the
+stencil, with respect to the current pixel. For example
+```julia
+stencil = CartesianIndex.([(0,1), (0,1), (1,1)])
+est = StencilSymbolicPermutation(stencil)
+```
+Here the stencil creates a 2x2 square extending to the bottom and right of the pixel.
+Notice that offset `0` (meaning the pixel itself) is alwaus inclided automatically.
+The length of the stencil decides the order of the permutation entropy, and the ordering
+within the stencil lifts any ambiguity regarding the order that pixels are compared to.
+The pixel without any offset is always first in the order.
+
+After having defined `est`, one calculates the permutation entropy of ordinal patterns
+by calling [`genentropy`](@ref) with `est`, and with the data that will be a matrix.
+To apply this to timeseries of spatial patterns, simply loop over the call, e.g.:
+```julia
+data = [rand(50,50) for _ in 1:50]
+entropies = [genentropy(x, est) for x in data]
+```
+
+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.
+
+[^Ribeiro2012]:
+    Ribeiro et al. (2012). Complexity-entropy causality plane as a complexity measure
+    for two-dimensional patterns. https://doi.org/10.1371/journal.pone.0040689
+
+[^Schlemmer2018]:
+    Schlemmer et al. (2012). Spatiotemporal Permutation Entropy as a Measure for
+    Complexity of Cardiac Arrhythmia. https://doi.org/10.3389/fphy.2018.00039
+"""
+struct SpatiotemporalPermutation{D,P} <: ProbabilityEstimator
+    stencil::Vector{CartesianIndex{D}}
+end
+function SpatiotemporalPermutation(stencil::Vector{CartesianIndex{D}}, p::Bool) where {D}
+    SpatiotemporalPermutation{D, p}(stencil)
+end

From ef8926c3f883c9e4c28fdf68e31c5437fb873c43 Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Sun, 4 Sep 2022 09:37:24 +0100
Subject: [PATCH 02/11] Working implementation for periodic

---
 src/symbolic/symbolic_stencils.jl | 68 ++++++++++++++++++++++++++-----
 1 file changed, 58 insertions(+), 10 deletions(-)

diff --git a/src/symbolic/symbolic_stencils.jl b/src/symbolic/symbolic_stencils.jl
index a43ac2254..7ec8e5767 100644
--- a/src/symbolic/symbolic_stencils.jl
+++ b/src/symbolic/symbolic_stencils.jl
@@ -2,7 +2,7 @@
     SpatiotemporalPermutation(stencil, periodic = true)
 A symbolic, permutation-based probabilities/entropy estimator for higher-dimensional
 data, typically representing spatiotemporal systems. The data are higher-dimensional
-arrays, such as 2D [^Ribeiro2012] or 3D [^Schlemmer2018] `Matrix`.
+arrays, such as 2D [^Ribeiro2012] or 3D [^Schlemmer2018].
 This approach is also known as _Spatiotemporal Permutation Entropy_.
 
 A _stencil_ defines what local area around each pixel to
@@ -11,19 +11,19 @@ vectors of `CartesianIndex` which encode the _offsets_ of the pixes to include i
 stencil, with respect to the current pixel. For example
 ```julia
 stencil = CartesianIndex.([(0,1), (0,1), (1,1)])
-est = StencilSymbolicPermutation(stencil)
+est = SpatiotemporalPermutation(stencil)
 ```
 Here the stencil creates a 2x2 square extending to the bottom and right of the pixel.
-Notice that offset `0` (meaning the pixel itself) is alwaus inclided automatically.
+Notice that no offset (meaning the pixel itself) is always included automatically.
 The length of the stencil decides the order of the permutation entropy, and the ordering
-within the stencil lifts any ambiguity regarding the order that pixels are compared to.
+within the stencil dictates the order that pixels are compared with.
 The pixel without any offset is always first in the order.
 
 After having defined `est`, one calculates the permutation entropy of ordinal patterns
-by calling [`genentropy`](@ref) with `est`, and with the data that will be a matrix.
-To apply this to timeseries of spatial patterns, simply loop over the call, e.g.:
+by calling [`genentropy`](@ref) with `est`, and with the array data.
+To apply this to timeseries of spatial data, simply loop over the call, e.g.:
 ```julia
-data = [rand(50,50) for _ in 1:50]
+data = [rand(50, 50) for _ in 1:50]
 entropies = [genentropy(x, est) for x in data]
 ```
 
@@ -38,9 +38,57 @@ array. If `periodic = false`, pixels whose stencil exceeds the array bounds are
     Schlemmer et al. (2012). Spatiotemporal Permutation Entropy as a Measure for
     Complexity of Cardiac Arrhythmia. https://doi.org/10.3389/fphy.2018.00039
 """
-struct SpatiotemporalPermutation{D,P} <: ProbabilityEstimator
+struct SpatiotemporalPermutation{D,P} <: ProbabilitiesEstimator
     stencil::Vector{CartesianIndex{D}}
+    viewer::Vector{CartesianIndex{D}}
 end
-function SpatiotemporalPermutation(stencil::Vector{CartesianIndex{D}}, p::Bool) where {D}
-    SpatiotemporalPermutation{D, p}(stencil)
+function SpatiotemporalPermutation(
+        stencil::Vector{CartesianIndex{D}}, p::Bool = true
+    ) where {D}
+    # Ensure that no offset is part of the stencil
+    stencil = pushfirst!(copy(stencil), CartesianIndex{D}(zeros(Int, D)...))
+    SpatiotemporalPermutation{D, p}(stencil, copy(stencil))
 end
+
+# This source code is a modification of the code of Agents.jl that finds neighbors
+# in grid-like spaces. It's the code of `nearby_positions` in `grid_general.jl`.
+function pixels_in_stencil(array, pixel, spatperm::SpatiotemporalPermutation{D,false}) where {D}
+    # TODO: Not clear yet how out of bounds is handled in this case.
+    @inbounds for i in eachindex(spatperm.stencil)
+        statperm.viewer[i] = spatperm.stencil[i] + pixel
+    end
+    return statperm.viewer
+    # space_size = size(array)
+    # positions_iterator = (n + pixel for n in spatperm.stencil)
+    # return Base.Iterators.filter(
+    #     pos -> checkbounds(Bool, array, pos...), positions_iterator
+    # )
+end
+
+function pixels_in_stencil(array, pixel, spatperm::SpatiotemporalPermutation{D,true}) where {D}
+    space_size = size(array)
+    @inbounds for i in eachindex(spatperm.stencil)
+        spatperm.viewer[i] = CartesianIndex{D}(mod1.(Tuple(spatperm.stencil[i] + pixel), space_size))
+    end
+    return spatperm.viewer
+end
+
+function Entropies.probabilities(x, est::SpatiotemporalPermutation) where {m, T}
+    s = zeros(Int, length(x))
+    probabilities!(s, x, est)
+end
+
+function Entropies.probabilities!(s::AbstractVector{Int}, x, est::SpatiotemporalPermutation)
+    m = length(est.stencil)
+    for (i, pixel) in enumerate(CartesianIndices(x))
+        pixels = pixels_in_stencil(x, pixel, est)
+        s[i] = Entropies.encode_motif(view(x, pixels), m)
+    end
+    return probabilities(s)
+end
+
+# %%
+stencil = CartesianIndex.([(0,1), (0,1), (1,1)])
+est = SpatiotemporalPermutation(stencil)
+x = rand(50, 50)
+p = probabilities(x, est)

From b531961e0c49b9b65f172d1871867e49d0ea863a Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Sun, 4 Sep 2022 10:25:28 +0100
Subject: [PATCH 03/11] Working code both periodic and not

---
 src/symbolic/symbolic_stencils.jl | 74 +++++++++++++++++++------------
 1 file changed, 45 insertions(+), 29 deletions(-)

diff --git a/src/symbolic/symbolic_stencils.jl b/src/symbolic/symbolic_stencils.jl
index 7ec8e5767..7a1ce336a 100644
--- a/src/symbolic/symbolic_stencils.jl
+++ b/src/symbolic/symbolic_stencils.jl
@@ -1,19 +1,23 @@
 """
-    SpatiotemporalPermutation(stencil, periodic = true)
-A symbolic, permutation-based probabilities/entropy estimator for higher-dimensional
-data, typically representing spatiotemporal systems. The data are higher-dimensional
-arrays, such as 2D [^Ribeiro2012] or 3D [^Schlemmer2018].
+    SpatiotemporalPermutation(stencil, x, periodic = true)
+A symbolic, permutation-based probabilities/entropy estimator for spatiotemporal systems.
+The data are a high-dimensional array `x`, such as 2D [^Ribeiro2012] or 3D [^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.
 
 A _stencil_ defines what local area around each pixel to
 consider, and compute the ordinal pattern within the stencil. Stencils are given as
 vectors of `CartesianIndex` which encode the _offsets_ of the pixes to include in the
-stencil, with respect to the current pixel. For example
+stencil, with respect to the current pixel. For simplicity, only zero or positive offsets
+are allowed. For example
 ```julia
+data = [rand(50, 50) for _ in 1:50]
+x = data[1] # first "time slice" of a spatial system evolution
 stencil = CartesianIndex.([(0,1), (0,1), (1,1)])
-est = SpatiotemporalPermutation(stencil)
+est = SpatiotemporalPermutation(stencil, x)
 ```
-Here the stencil creates a 2x2 square extending to the bottom and right of the pixel.
+Here the stencil creates a 2x2 square extending to the bottom and right of the pixel
+(directions here correspond to the way Julia prints matrices by default).
 Notice that no offset (meaning the pixel itself) is always included automatically.
 The length of the stencil decides the order of the permutation entropy, and the ordering
 within the stencil dictates the order that pixels are compared with.
@@ -23,8 +27,8 @@ After having defined `est`, one calculates the permutation entropy of ordinal pa
 by calling [`genentropy`](@ref) with `est`, and with the array data.
 To apply this to timeseries of spatial data, simply loop over the call, e.g.:
 ```julia
-data = [rand(50, 50) for _ in 1:50]
-entropies = [genentropy(x, est) for x in data]
+entropy = genentropy(x, est)
+entropy_vs_time = genentropy.(data, est) # broadcasting with `.`
 ```
 
 The argument `periodic` decides whether the stencil should wrap around at the end of the
@@ -38,50 +42,58 @@ array. If `periodic = false`, pixels whose stencil exceeds the array bounds are
     Schlemmer et al. (2012). Spatiotemporal Permutation Entropy as a Measure for
     Complexity of Cardiac Arrhythmia. https://doi.org/10.3389/fphy.2018.00039
 """
-struct SpatiotemporalPermutation{D,P} <: ProbabilitiesEstimator
+struct SpatiotemporalPermutation{D,P,V} <: ProbabilitiesEstimator
     stencil::Vector{CartesianIndex{D}}
     viewer::Vector{CartesianIndex{D}}
+    arraysize::Dims{D}
+    valid::V
 end
 function SpatiotemporalPermutation(
-        stencil::Vector{CartesianIndex{D}}, p::Bool = true
+        stencil::Vector{CartesianIndex{D}}, x::AbstractArray, p::Bool = true
     ) where {D}
     # Ensure that no offset is part of the stencil
     stencil = pushfirst!(copy(stencil), CartesianIndex{D}(zeros(Int, D)...))
-    SpatiotemporalPermutation{D, p}(stencil, copy(stencil))
+    arraysize = size(x)
+    # Store valid indices for later iteration
+    if p
+        valid = CartesianIndices(x)
+    else
+        # collect maximum offsets in each dimension for limiting ranges
+        maxoffsets = [maximum(s[i] for s in stencil) for i in 1:D]
+        ranges = Iterators.product([1:(arraysize[i]-maxoffsets[i]) for i in 1:D]...)
+        valid = Base.Generator(idxs -> CartesianIndex{D}(idxs), ranges)
+    end
+    SpatiotemporalPermutation{D, p, typeof(valid)}(stencil, copy(stencil), arraysize, valid)
 end
 
 # This source code is a modification of the code of Agents.jl that finds neighbors
 # in grid-like spaces. It's the code of `nearby_positions` in `grid_general.jl`.
-function pixels_in_stencil(array, pixel, spatperm::SpatiotemporalPermutation{D,false}) where {D}
-    # TODO: Not clear yet how out of bounds is handled in this case.
+function pixels_in_stencil(pixel, spatperm::SpatiotemporalPermutation{D,false}) where {D}
     @inbounds for i in eachindex(spatperm.stencil)
-        statperm.viewer[i] = spatperm.stencil[i] + pixel
+        spatperm.viewer[i] = spatperm.stencil[i] + pixel
     end
-    return statperm.viewer
-    # space_size = size(array)
-    # positions_iterator = (n + pixel for n in spatperm.stencil)
-    # return Base.Iterators.filter(
-    #     pos -> checkbounds(Bool, array, pos...), positions_iterator
-    # )
+    return spatperm.viewer
 end
 
-function pixels_in_stencil(array, pixel, spatperm::SpatiotemporalPermutation{D,true}) where {D}
-    space_size = size(array)
+function pixels_in_stencil(pixel, spatperm::SpatiotemporalPermutation{D,true}) where {D}
+    arraysize = spatperm.arraysize
     @inbounds for i in eachindex(spatperm.stencil)
-        spatperm.viewer[i] = CartesianIndex{D}(mod1.(Tuple(spatperm.stencil[i] + pixel), space_size))
+        # It's annoying that we have to change to tuple and then to CartesianIndex
+        # because iteration over cartesian indices is not allowed. But oh well.
+        spatperm.viewer[i] = CartesianIndex{D}(mod1.(Tuple(spatperm.stencil[i] + pixel), arraysize))
     end
     return spatperm.viewer
 end
 
-function Entropies.probabilities(x, est::SpatiotemporalPermutation) where {m, T}
-    s = zeros(Int, length(x))
+function Entropies.probabilities(x, est::SpatiotemporalPermutation)
+    s = zeros(Int, length(est.valid))
     probabilities!(s, x, est)
 end
 
 function Entropies.probabilities!(s::AbstractVector{Int}, x, est::SpatiotemporalPermutation)
     m = length(est.stencil)
-    for (i, pixel) in enumerate(CartesianIndices(x))
-        pixels = pixels_in_stencil(x, pixel, est)
+    for (i, pixel) in enumerate(est.valid)
+        pixels = pixels_in_stencil(pixel, est)
         s[i] = Entropies.encode_motif(view(x, pixels), m)
     end
     return probabilities(s)
@@ -89,6 +101,10 @@ end
 
 # %%
 stencil = CartesianIndex.([(0,1), (0,1), (1,1)])
-est = SpatiotemporalPermutation(stencil)
 x = rand(50, 50)
+est = SpatiotemporalPermutation(stencil, x, false)
 p = probabilities(x, est)
+
+using BenchmarkTools
+s = zeros(Int, length(est.valid))
+@btime probabilities!($s, $x, $est)
\ No newline at end of file

From 676c1440c6300aa1a488ceb9f069ad8acf7cfcd9 Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Sun, 4 Sep 2022 17:07:40 +0100
Subject: [PATCH 04/11] add sanity check for negative ofsets

---
 src/symbolic/symbolic_stencils.jl | 20 +++++++-------------
 1 file changed, 7 insertions(+), 13 deletions(-)

diff --git a/src/symbolic/symbolic_stencils.jl b/src/symbolic/symbolic_stencils.jl
index 7a1ce336a..705914c55 100644
--- a/src/symbolic/symbolic_stencils.jl
+++ b/src/symbolic/symbolic_stencils.jl
@@ -60,6 +60,9 @@ function SpatiotemporalPermutation(
     else
         # collect maximum offsets in each dimension for limiting ranges
         maxoffsets = [maximum(s[i] for s in stencil) for i in 1:D]
+        # Safety check
+        minoffsets = [minimum(s[i] for s in stencil) for i in 1:D]
+        @assert all(≥(0), minoffsets) "Offsets in stencil must be ≥ 0!"
         ranges = Iterators.product([1:(arraysize[i]-maxoffsets[i]) for i in 1:D]...)
         valid = Base.Generator(idxs -> CartesianIndex{D}(idxs), ranges)
     end
@@ -76,11 +79,12 @@ function pixels_in_stencil(pixel, spatperm::SpatiotemporalPermutation{D,false})
 end
 
 function pixels_in_stencil(pixel, spatperm::SpatiotemporalPermutation{D,true}) where {D}
-    arraysize = spatperm.arraysize
     @inbounds for i in eachindex(spatperm.stencil)
         # It's annoying that we have to change to tuple and then to CartesianIndex
         # because iteration over cartesian indices is not allowed. But oh well.
-        spatperm.viewer[i] = CartesianIndex{D}(mod1.(Tuple(spatperm.stencil[i] + pixel), arraysize))
+        spatperm.viewer[i] = CartesianIndex{D}(
+            mod1.(Tuple(spatperm.stencil[i] + pixel), spatperm.arraysize)
+        )
     end
     return spatperm.viewer
 end
@@ -97,14 +101,4 @@ function Entropies.probabilities!(s::AbstractVector{Int}, x, est::Spatiotemporal
         s[i] = Entropies.encode_motif(view(x, pixels), m)
     end
     return probabilities(s)
-end
-
-# %%
-stencil = CartesianIndex.([(0,1), (0,1), (1,1)])
-x = rand(50, 50)
-est = SpatiotemporalPermutation(stencil, x, false)
-p = probabilities(x, est)
-
-using BenchmarkTools
-s = zeros(Int, length(est.valid))
-@btime probabilities!($s, $x, $est)
\ No newline at end of file
+end
\ No newline at end of file

From 1575ddd8dd5f1c719dd404d2e2c0f366b8425e61 Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Sun, 4 Sep 2022 17:34:55 +0100
Subject: [PATCH 05/11] add tests

---
 src/symbolic/symbolic_stencils.jl | 22 +++++++++++++---------
 test/runtests.jl                  |  2 ++
 test/spatial_permutation_tests.jl | 16 ++++++++++++++++
 3 files changed, 31 insertions(+), 9 deletions(-)
 create mode 100644 test/spatial_permutation_tests.jl

diff --git a/src/symbolic/symbolic_stencils.jl b/src/symbolic/symbolic_stencils.jl
index 705914c55..3658224d2 100644
--- a/src/symbolic/symbolic_stencils.jl
+++ b/src/symbolic/symbolic_stencils.jl
@@ -1,5 +1,7 @@
+export SpatialSymbolicPermutation
+
 """
-    SpatiotemporalPermutation(stencil, x, periodic = true)
+    SpatialSymbolicPermutation(stencil, x, periodic = true)
 A symbolic, permutation-based probabilities/entropy estimator for spatiotemporal systems.
 The data are a high-dimensional array `x`, such as 2D [^Ribeiro2012] or 3D [^Schlemmer2018].
 This approach is also known as _Spatiotemporal Permutation Entropy_.
@@ -14,7 +16,7 @@ are allowed. For example
 data = [rand(50, 50) for _ in 1:50]
 x = data[1] # first "time slice" of a spatial system evolution
 stencil = CartesianIndex.([(0,1), (0,1), (1,1)])
-est = SpatiotemporalPermutation(stencil, x)
+est = SpatialSymbolicPermutation(stencil, x)
 ```
 Here the stencil creates a 2x2 square extending to the bottom and right of the pixel
 (directions here correspond to the way Julia prints matrices by default).
@@ -42,13 +44,13 @@ array. If `periodic = false`, pixels whose stencil exceeds the array bounds are
     Schlemmer et al. (2012). Spatiotemporal Permutation Entropy as a Measure for
     Complexity of Cardiac Arrhythmia. https://doi.org/10.3389/fphy.2018.00039
 """
-struct SpatiotemporalPermutation{D,P,V} <: ProbabilitiesEstimator
+struct SpatialSymbolicPermutation{D,P,V} <: ProbabilitiesEstimator
     stencil::Vector{CartesianIndex{D}}
     viewer::Vector{CartesianIndex{D}}
     arraysize::Dims{D}
     valid::V
 end
-function SpatiotemporalPermutation(
+function SpatialSymbolicPermutation(
         stencil::Vector{CartesianIndex{D}}, x::AbstractArray, p::Bool = true
     ) where {D}
     # Ensure that no offset is part of the stencil
@@ -66,19 +68,19 @@ function SpatiotemporalPermutation(
         ranges = Iterators.product([1:(arraysize[i]-maxoffsets[i]) for i in 1:D]...)
         valid = Base.Generator(idxs -> CartesianIndex{D}(idxs), ranges)
     end
-    SpatiotemporalPermutation{D, p, typeof(valid)}(stencil, copy(stencil), arraysize, valid)
+    SpatialSymbolicPermutation{D, p, typeof(valid)}(stencil, copy(stencil), arraysize, valid)
 end
 
 # This source code is a modification of the code of Agents.jl that finds neighbors
 # in grid-like spaces. It's the code of `nearby_positions` in `grid_general.jl`.
-function pixels_in_stencil(pixel, spatperm::SpatiotemporalPermutation{D,false}) where {D}
+function pixels_in_stencil(pixel, spatperm::SpatialSymbolicPermutation{D,false}) where {D}
     @inbounds for i in eachindex(spatperm.stencil)
         spatperm.viewer[i] = spatperm.stencil[i] + pixel
     end
     return spatperm.viewer
 end
 
-function pixels_in_stencil(pixel, spatperm::SpatiotemporalPermutation{D,true}) where {D}
+function pixels_in_stencil(pixel, spatperm::SpatialSymbolicPermutation{D,true}) where {D}
     @inbounds for i in eachindex(spatperm.stencil)
         # It's annoying that we have to change to tuple and then to CartesianIndex
         # because iteration over cartesian indices is not allowed. But oh well.
@@ -89,12 +91,14 @@ function pixels_in_stencil(pixel, spatperm::SpatiotemporalPermutation{D,true}) w
     return spatperm.viewer
 end
 
-function Entropies.probabilities(x, est::SpatiotemporalPermutation)
+function Entropies.probabilities(x, est::SpatialSymbolicPermutation)
+    # TODO: This can be literally a call to symbolize and then
+    # calling probabilities on it. Should do once the symbolize refactoring is done.
     s = zeros(Int, length(est.valid))
     probabilities!(s, x, est)
 end
 
-function Entropies.probabilities!(s::AbstractVector{Int}, x, est::SpatiotemporalPermutation)
+function Entropies.probabilities!(s::AbstractVector{Int}, x, est::SpatialSymbolicPermutation)
     m = length(est.stencil)
     for (i, pixel) in enumerate(est.valid)
         pixels = pixels_in_stencil(pixel, est)
diff --git a/test/runtests.jl b/test/runtests.jl
index a9392b602..a028ea0b4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -377,3 +377,5 @@ end
         @test de >= 0.0
     end
 end
+
+include("spatial_permutation_tests.jl")
\ No newline at end of file
diff --git a/test/spatial_permutation_tests.jl b/test/spatial_permutation_tests.jl
new file mode 100644
index 000000000..5906efaad
--- /dev/null
+++ b/test/spatial_permutation_tests.jl
@@ -0,0 +1,16 @@
+using Entropies, Test
+
+@testset "Spatiotemporal Permutation Entropy" begin
+    x = [1 2 1; 8 3 4; 6 7 5]
+    # Re-create the Ribeiro et al, 2012 using stencil instead of specifying m
+    stencil = CartesianIndex.([(0,1), (0,1),(1,1)])
+    est = SpatiotemporalPermutation(stencil, x, false)
+
+    # Generic tests
+    ps = probabilities(x, est)
+    @test ps isa Probabilities
+
+    # test case from Ribeiro et al, 2012
+    h = genentropy(x, est, base = MathConstants.e)
+    @test round(h, digits = 2) ≈ 0.33
+end
\ No newline at end of file

From e56fc07eafa6b4b1a2b73ba67e31d3b9a7af4456 Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Mon, 5 Sep 2022 07:35:59 +0100
Subject: [PATCH 06/11] allow also negative offsets

---
 src/symbolic/symbolic_stencils.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/symbolic/symbolic_stencils.jl b/src/symbolic/symbolic_stencils.jl
index 3658224d2..56c24303a 100644
--- a/src/symbolic/symbolic_stencils.jl
+++ b/src/symbolic/symbolic_stencils.jl
@@ -10,8 +10,7 @@ This approach is also known as _Spatiotemporal Permutation Entropy_.
 A _stencil_ defines what local area around each pixel to
 consider, and compute the ordinal pattern within the stencil. Stencils are given as
 vectors of `CartesianIndex` which encode the _offsets_ of the pixes to include in the
-stencil, with respect to the current pixel. For simplicity, only zero or positive offsets
-are allowed. For example
+stencil, with respect to the current pixel. For example
 ```julia
 data = [rand(50, 50) for _ in 1:50]
 x = data[1] # first "time slice" of a spatial system evolution
@@ -63,9 +62,10 @@ function SpatialSymbolicPermutation(
         # collect maximum offsets in each dimension for limiting ranges
         maxoffsets = [maximum(s[i] for s in stencil) for i in 1:D]
         # Safety check
-        minoffsets = [minimum(s[i] for s in stencil) for i in 1:D]
-        @assert all(≥(0), minoffsets) "Offsets in stencil must be ≥ 0!"
-        ranges = Iterators.product([1:(arraysize[i]-maxoffsets[i]) for i in 1:D]...)
+        minoffsets = [min(0, minimum(s[i] for s in stencil)) for i in 1:D]
+        ranges = Iterators.product(
+            [(1-minoffsets[i]):(arraysize[i]-maxoffsets[i]) for i in 1:D]...
+        )
         valid = Base.Generator(idxs -> CartesianIndex{D}(idxs), ranges)
     end
     SpatialSymbolicPermutation{D, p, typeof(valid)}(stencil, copy(stencil), arraysize, valid)

From d25f7f94ee86c24cd76fb303de222304db83ac92 Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Mon, 5 Sep 2022 07:36:43 +0100
Subject: [PATCH 07/11] correct age for schlemmer reference

---
 src/symbolic/symbolic_stencils.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/symbolic/symbolic_stencils.jl b/src/symbolic/symbolic_stencils.jl
index 56c24303a..0cca0f050 100644
--- a/src/symbolic/symbolic_stencils.jl
+++ b/src/symbolic/symbolic_stencils.jl
@@ -40,7 +40,7 @@ array. If `periodic = false`, pixels whose stencil exceeds the array bounds are
     for two-dimensional patterns. https://doi.org/10.1371/journal.pone.0040689
 
 [^Schlemmer2018]:
-    Schlemmer et al. (2012). Spatiotemporal Permutation Entropy as a Measure for
+    Schlemmer et al. (2018). Spatiotemporal Permutation Entropy as a Measure for
     Complexity of Cardiac Arrhythmia. https://doi.org/10.3389/fphy.2018.00039
 """
 struct SpatialSymbolicPermutation{D,P,V} <: ProbabilitiesEstimator

From 9902107a81dc03729d90d003622042c1e0c02e7a Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Mon, 5 Sep 2022 08:02:20 +0100
Subject: [PATCH 08/11] refer to symbolize in code

---
 src/symbolic/symbolic_stencils.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/symbolic/symbolic_stencils.jl b/src/symbolic/symbolic_stencils.jl
index 0cca0f050..d42922e21 100644
--- a/src/symbolic/symbolic_stencils.jl
+++ b/src/symbolic/symbolic_stencils.jl
@@ -92,8 +92,8 @@ function pixels_in_stencil(pixel, spatperm::SpatialSymbolicPermutation{D,true})
 end
 
 function Entropies.probabilities(x, est::SpatialSymbolicPermutation)
-    # TODO: This can be literally a call to symbolize and then
-    # calling probabilities on it. Should do once the symbolize refactoring is done.
+    # TODO: This can be literally a call to `symbolize` and then
+    # calling probabilities on it. Should do once the `symbolize` refactoring is done.
     s = zeros(Int, length(est.valid))
     probabilities!(s, x, est)
 end

From 982ec33e2f4f88dc5dadd2e2fda7730352550ab4 Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Tue, 6 Sep 2022 18:43:51 +0100
Subject: [PATCH 09/11] fix everything and the tests

---
 src/core.jl                                     |  4 ++--
 ...bolic_stencils.jl => spatial_permutation.jl} |  5 +++--
 src/symbolic/symbolic.jl                        |  1 +
 test/spatial_permutation_tests.jl               | 17 ++++++++++-------
 4 files changed, 16 insertions(+), 11 deletions(-)
 rename src/symbolic/{symbolic_stencils.jl => spatial_permutation.jl} (97%)

diff --git a/src/core.jl b/src/core.jl
index 9fe79a638..ddf66da7b 100644
--- a/src/core.jl
+++ b/src/core.jl
@@ -6,7 +6,7 @@ export probabilities, probabilities!
 export genentropy, genentropy!
 export Dataset, dimension
 
-const Array_or_Dataset = Union{AbstractArray, AbstractDataset}
+const Array_or_Dataset = Union{<:AbstractArray, <:AbstractDataset}
 
 """
     Probabilities(x) → p
@@ -149,7 +149,7 @@ function genentropy(prob::Probabilities; q = 1.0, α = nothing, base = MathConst
     end
 end
 
-genentropy(x::AbstractArray{<:Real}) =
+genentropy(::AbstractArray{<:Real}) =
     error("For single-argument input, do `genentropy(Probabilities(x))` instead.")
 
 function genentropy(x::Array_or_Dataset, est; q = 1.0, α = nothing, base = MathConstants.e)
diff --git a/src/symbolic/symbolic_stencils.jl b/src/symbolic/spatial_permutation.jl
similarity index 97%
rename from src/symbolic/symbolic_stencils.jl
rename to src/symbolic/spatial_permutation.jl
index d42922e21..301b56e4e 100644
--- a/src/symbolic/symbolic_stencils.jl
+++ b/src/symbolic/spatial_permutation.jl
@@ -14,7 +14,7 @@ stencil, with respect to the current pixel. For example
 ```julia
 data = [rand(50, 50) for _ in 1:50]
 x = data[1] # first "time slice" of a spatial system evolution
-stencil = CartesianIndex.([(0,1), (0,1), (1,1)])
+stencil = CartesianIndex.([(0,1), (1,1), (1,0)])
 est = SpatialSymbolicPermutation(stencil, x)
 ```
 Here the stencil creates a 2x2 square extending to the bottom and right of the pixel
@@ -98,11 +98,12 @@ function Entropies.probabilities(x, est::SpatialSymbolicPermutation)
     probabilities!(s, x, est)
 end
 
-function Entropies.probabilities!(s::AbstractVector{Int}, x, est::SpatialSymbolicPermutation)
+function Entropies.probabilities!(s, x, est::SpatialSymbolicPermutation)
     m = length(est.stencil)
     for (i, pixel) in enumerate(est.valid)
         pixels = pixels_in_stencil(pixel, est)
         s[i] = Entropies.encode_motif(view(x, pixels), m)
     end
+    @show s
     return probabilities(s)
 end
\ No newline at end of file
diff --git a/src/symbolic/symbolic.jl b/src/symbolic/symbolic.jl
index 99be6dbf1..dec62c7aa 100644
--- a/src/symbolic/symbolic.jl
+++ b/src/symbolic/symbolic.jl
@@ -10,6 +10,7 @@ include("utils.jl")
 include("SymbolicPermutation.jl")
 include("SymbolicWeightedPermutation.jl")
 include("SymbolicAmplitudeAware.jl")
+include("spatial_permutation.jl")
 
 """
     permentropy(x; τ = 1, m = 3, base = MathConstants.e)
diff --git a/test/spatial_permutation_tests.jl b/test/spatial_permutation_tests.jl
index 5906efaad..07454fa98 100644
--- a/test/spatial_permutation_tests.jl
+++ b/test/spatial_permutation_tests.jl
@@ -1,16 +1,19 @@
 using Entropies, Test
 
-@testset "Spatiotemporal Permutation Entropy" begin
+@testset "Spatial Permutation Entropy" begin
     x = [1 2 1; 8 3 4; 6 7 5]
-    # Re-create the Ribeiro et al, 2012 using stencil instead of specifying m
-    stencil = CartesianIndex.([(0,1), (0,1),(1,1)])
-    est = SpatiotemporalPermutation(stencil, x, false)
+    # Re-create the Ribeiro et al, 2012 using stencil
+    # (you get 4 symbols in a 3x3 matrix. For this matrix, the upper left
+    # and bottom right are the same symbol. So three probabilities in the end).
+    stencil = CartesianIndex.([(1,0), (0,1), (1,1)])
+    est = SpatialSymbolicPermutation(stencil, x, false)
 
     # Generic tests
     ps = probabilities(x, est)
     @test ps isa Probabilities
+    @test length(ps) == 3
+    @test sort(ps) == [0.25, 0.25, 0.5]
 
-    # test case from Ribeiro et al, 2012
-    h = genentropy(x, est, base = MathConstants.e)
-    @test round(h, digits = 2) ≈ 0.33
+    h = genentropy(x, est, base = 2)
+    @test h == 1.5
 end
\ No newline at end of file

From 9cfc7dcb0c878c8a57a146df9269c875e9722a2b Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Tue, 6 Sep 2022 18:45:15 +0100
Subject: [PATCH 10/11] add changelog entry

---
 CHANGELOG.md | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index be15e51a4..05d4cc4b1 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,6 +1,9 @@
 # CHANGELOG
 Changelog is kept with respect to version 0.11 of Entropies.jl.
 
+## main
+* New probability estimator `SpatialSymbolicPermutation` suitable for computing spatial permutation entropies
+
 ## 1.1
 * Introduce convenience function `permentropy`.
 * Several type instabilities fixed.

From b7e2391126dd2b7dc45b244ee5207093d293c4f2 Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Tue, 6 Sep 2022 19:02:08 +0100
Subject: [PATCH 11/11] add one more test and pretty printing

---
 src/symbolic/spatial_permutation.jl |  8 ++++++++
 test/spatial_permutation_tests.jl   | 30 +++++++++++++++++++++++++++++
 2 files changed, 38 insertions(+)

diff --git a/src/symbolic/spatial_permutation.jl b/src/symbolic/spatial_permutation.jl
index 301b56e4e..fa6ca1e96 100644
--- a/src/symbolic/spatial_permutation.jl
+++ b/src/symbolic/spatial_permutation.jl
@@ -55,6 +55,7 @@ function SpatialSymbolicPermutation(
     # Ensure that no offset is part of the stencil
     stencil = pushfirst!(copy(stencil), CartesianIndex{D}(zeros(Int, D)...))
     arraysize = size(x)
+    @assert length(arraysize) == D "Indices and input array must match dimensionality!"
     # Store valid indices for later iteration
     if p
         valid = CartesianIndices(x)
@@ -106,4 +107,11 @@ function Entropies.probabilities!(s, x, est::SpatialSymbolicPermutation)
     end
     @show s
     return probabilities(s)
+end
+
+# Pretty printing
+function Base.show(io::IO, est::SpatialSymbolicPermutation{D}) where {D}
+    print(io, "Spatial permutation estimator for $D-dimensional data. Stencil:")
+    print(io, "\n")
+    print(io, est.stencil)
 end
\ No newline at end of file
diff --git a/test/spatial_permutation_tests.jl b/test/spatial_permutation_tests.jl
index 07454fa98..ee49f08ea 100644
--- a/test/spatial_permutation_tests.jl
+++ b/test/spatial_permutation_tests.jl
@@ -16,4 +16,34 @@ using Entropies, Test
 
     h = genentropy(x, est, base = 2)
     @test h == 1.5
+
+    # In fact, doesn't matter how we order the stencil,
+    # the symbols will always be equal in top-left and bottom-right
+    stencil = CartesianIndex.([(1,0), (1,1), (0,1)])
+    est = SpatialSymbolicPermutation(stencil, x, false)
+    @test genentropy(x, est, base = 2) == 1.5
+
+    # But for sanity, let's ensure we get a different number
+    # for a different stencil
+    stencil = CartesianIndex.([(1,0), (2,0)])
+    est = SpatialSymbolicPermutation(stencil, x, false)
+    ps = sort(probabilities(x, est))
+    @test ps[1] == 1/3
+    @test ps[2] == 2/3
+
+    # TODO: Symbolize tests once its part of the API.
+
+    # Also test that it works for arbitrarily high-dimensional arrays
+    stencil = CartesianIndex.([(0,1,0), (0,0,1), (1,0,0)])
+    z = reshape(1:125, (5,5,5))
+    est = SpatialSymbolicPermutation(stencil, z, false)
+    # Analytically the total stencils are of length 4*4*4 = 64
+    # but all of them given the same probabilities because of the layout
+    ps = probabilities(z, est)
+    @test ps == [1]
+    # if we shuffle, we get random stuff
+    using Random
+    w = shuffle!(Random.MersenneTwister(42), collect(z))
+    ps = probabilities(w, est)
+    @test length(ps) > 1
 end
\ No newline at end of file