1+ export SpatialSymbolicPermutation
2+
3+ """
4+ SpatialSymbolicPermutation(stencil, x, periodic = true)
5+ A symbolic, permutation-based probabilities/entropy estimator for spatiotemporal systems.
6+ The data are a high-dimensional array `x`, such as 2D [^Ribeiro2012] or 3D [^Schlemmer2018].
7+ This approach is also known as _Spatiotemporal Permutation Entropy_.
8+ `x` is given because we need to know its size for optimization and bound checking.
9+
10+ A _stencil_ defines what local area around each pixel to
11+ consider, and compute the ordinal pattern within the stencil. Stencils are given as
12+ vectors of `CartesianIndex` which encode the _offsets_ of the pixes to include in the
13+ stencil, with respect to the current pixel. For example
14+ ```julia
15+ data = [rand(50, 50) for _ in 1:50]
16+ x = data[1] # first "time slice" of a spatial system evolution
17+ stencil = CartesianIndex.([(0,1), (1,1), (1,0)])
18+ est = SpatialSymbolicPermutation(stencil, x)
19+ ```
20+ Here the stencil creates a 2x2 square extending to the bottom and right of the pixel
21+ (directions here correspond to the way Julia prints matrices by default).
22+ Notice that no offset (meaning the pixel itself) is always included automatically.
23+ The length of the stencil decides the order of the permutation entropy, and the ordering
24+ within the stencil dictates the order that pixels are compared with.
25+ The pixel without any offset is always first in the order.
26+
27+ After having defined `est`, one calculates the permutation entropy of ordinal patterns
28+ by calling [`genentropy`](@ref) with `est`, and with the array data.
29+ To apply this to timeseries of spatial data, simply loop over the call, e.g.:
30+ ```julia
31+ entropy = genentropy(x, est)
32+ entropy_vs_time = genentropy.(data, est) # broadcasting with `.`
33+ ```
34+
35+ The argument `periodic` decides whether the stencil should wrap around at the end of the
36+ array. If `periodic = false`, pixels whose stencil exceeds the array bounds are skipped.
37+
38+ [^Ribeiro2012]:
39+ Ribeiro et al. (2012). Complexity-entropy causality plane as a complexity measure
40+ for two-dimensional patterns. https://doi.org/10.1371/journal.pone.0040689
41+
42+ [^Schlemmer2018]:
43+ Schlemmer et al. (2018). Spatiotemporal Permutation Entropy as a Measure for
44+ Complexity of Cardiac Arrhythmia. https://doi.org/10.3389/fphy.2018.00039
45+ """
46+ struct SpatialSymbolicPermutation{D,P,V} <: ProbabilitiesEstimator
47+ stencil:: Vector{CartesianIndex{D}}
48+ viewer:: Vector{CartesianIndex{D}}
49+ arraysize:: Dims{D}
50+ valid:: V
51+ end
52+ function SpatialSymbolicPermutation (
53+ stencil:: Vector{CartesianIndex{D}} , x:: AbstractArray , p:: Bool = true
54+ ) where {D}
55+ # Ensure that no offset is part of the stencil
56+ stencil = pushfirst! (copy (stencil), CartesianIndex {D} (zeros (Int, D)... ))
57+ arraysize = size (x)
58+ @assert length (arraysize) == D " Indices and input array must match dimensionality!"
59+ # Store valid indices for later iteration
60+ if p
61+ valid = CartesianIndices (x)
62+ else
63+ # collect maximum offsets in each dimension for limiting ranges
64+ maxoffsets = [maximum (s[i] for s in stencil) for i in 1 : D]
65+ # Safety check
66+ minoffsets = [min (0 , minimum (s[i] for s in stencil)) for i in 1 : D]
67+ ranges = Iterators. product (
68+ [(1 - minoffsets[i]): (arraysize[i]- maxoffsets[i]) for i in 1 : D]. ..
69+ )
70+ valid = Base. Generator (idxs -> CartesianIndex {D} (idxs), ranges)
71+ end
72+ SpatialSymbolicPermutation {D, p, typeof(valid)} (stencil, copy (stencil), arraysize, valid)
73+ end
74+
75+ # This source code is a modification of the code of Agents.jl that finds neighbors
76+ # in grid-like spaces. It's the code of `nearby_positions` in `grid_general.jl`.
77+ function pixels_in_stencil (pixel, spatperm:: SpatialSymbolicPermutation{D,false} ) where {D}
78+ @inbounds for i in eachindex (spatperm. stencil)
79+ spatperm. viewer[i] = spatperm. stencil[i] + pixel
80+ end
81+ return spatperm. viewer
82+ end
83+
84+ function pixels_in_stencil (pixel, spatperm:: SpatialSymbolicPermutation{D,true} ) where {D}
85+ @inbounds for i in eachindex (spatperm. stencil)
86+ # It's annoying that we have to change to tuple and then to CartesianIndex
87+ # because iteration over cartesian indices is not allowed. But oh well.
88+ spatperm. viewer[i] = CartesianIndex {D} (
89+ mod1 .(Tuple (spatperm. stencil[i] + pixel), spatperm. arraysize)
90+ )
91+ end
92+ return spatperm. viewer
93+ end
94+
95+ function Entropies. probabilities (x, est:: SpatialSymbolicPermutation )
96+ # TODO : This can be literally a call to `symbolize` and then
97+ # calling probabilities on it. Should do once the `symbolize` refactoring is done.
98+ s = zeros (Int, length (est. valid))
99+ probabilities! (s, x, est)
100+ end
101+
102+ function Entropies. probabilities! (s, x, est:: SpatialSymbolicPermutation )
103+ m = length (est. stencil)
104+ for (i, pixel) in enumerate (est. valid)
105+ pixels = pixels_in_stencil (pixel, est)
106+ s[i] = Entropies. encode_motif (view (x, pixels), m)
107+ end
108+ @show s
109+ return probabilities (s)
110+ end
111+
112+ # Pretty printing
113+ function Base. show (io:: IO , est:: SpatialSymbolicPermutation{D} ) where {D}
114+ print (io, " Spatial permutation estimator for $D -dimensional data. Stencil:"
115+ print (io, " \n "
116+ print (io, est. stencil)
117+ end
0 commit comments