"""Fractal functions"""
import numpy as np
from numba import jit
from math import log, floor
from .entropy import num_zerocross
from .utils import _linear_regression, _log_n
all = ['petrosian_fd', 'katz_fd', 'higuchi_fd', 'detrended_fluctuation']
[docs]def petrosian_fd(x, axis=-1):
"""Petrosian fractal dimension.
Parameters
----------
x : list or np.array
1D or N-D data.
axis : int
The axis along which the FD is calculated. Default is -1 (last).
Returns
-------
pfd : float
Petrosian fractal dimension.
Notes
-----
The Petrosian fractal dimension of a time-series :math:`x` is defined by:
.. math:: P = \\frac{\\log_{10}(N)}{\\log_{10}(N) +
\\log_{10}(\\frac{N}{N+0.4N_{\\delta}})}
where :math:`N` is the length of the time series, and
:math:`N_{\\delta}` is the number of sign changes in the signal derivative.
Original code from the `pyrem <https://github.com/gilestrolab/pyrem>`_
package by Quentin Geissmann.
References
----------
* A. Petrosian, Kolmogorov complexity of finite sequences and
recognition of different preictal EEG patterns, in , Proceedings of the
Eighth IEEE Symposium on Computer-Based Medical Systems, 1995,
pp. 212-217.
* Goh, Cindy, et al. "Comparison of fractal dimension algorithms for
the computation of EEG biomarkers for dementia." 2nd International
Conference on Computational Intelligence in Medicine and Healthcare
(CIMED2005). 2005.
Examples
--------
>>> import numpy as np
>>> import entropy as ent
>>> import stochastic.processes.noise as sn
>>> rng = np.random.default_rng(seed=42)
>>> x = sn.FractionalGaussianNoise(hurst=0.5, rng=rng).sample(10000)
>>> print(f"{ent.petrosian_fd(x):.4f}")
1.0264
Fractional Gaussian noise with H = 0.9
>>> rng = np.random.default_rng(seed=42)
>>> x = sn.FractionalGaussianNoise(hurst=0.9, rng=rng).sample(10000)
>>> print(f"{ent.petrosian_fd(x):.4f}")
1.0235
Fractional Gaussian noise with H = 0.1
>>> rng = np.random.default_rng(seed=42)
>>> x = sn.FractionalGaussianNoise(hurst=0.1, rng=rng).sample(10000)
>>> print(f"{ent.petrosian_fd(x):.4f}")
1.0283
Random
>>> rng = np.random.default_rng(seed=42)
>>> print(f"{ent.petrosian_fd(rng.random(1000)):.4f}")
1.0350
Pure sine wave
>>> x = np.sin(2 * np.pi * 1 * np.arange(3000) / 100)
>>> print(f"{ent.petrosian_fd(x):.4f}")
1.0010
Linearly-increasing time-series (should be 1)
>>> x = np.arange(1000)
>>> print(f"{ent.petrosian_fd(x):.4f}")
1.0000
"""
x = np.asarray(x)
N = x.shape[axis]
# Number of sign changes in the first derivative of the signal
nzc_deriv = num_zerocross(np.diff(x, axis=axis), axis=axis)
pfd = np.log10(N) / (np.log10(N) + np.log10(N / (N + 0.4 * nzc_deriv)))
return pfd
[docs]def katz_fd(x, axis=-1):
"""Katz Fractal Dimension.
Parameters
----------
x : list or np.array
1D or N-D data.
axis : int
The axis along which the FD is calculated. Default is -1 (last).
Returns
-------
kfd : float
Katz fractal dimension.
Notes
-----
Katz’s method calculates the fractal dimension of a sample as follows:
the sum and average of the Euclidean distances between the successive
points of the sample (:math:`L` and :math:`a` , resp.) are calculated as
well as the maximum distance between the first point and any other point
of the sample (:math:`d`). The fractal dimension of the sample (:math:`D`)
then becomes:
.. math::
D = \\frac{\\log_{10}(L/a)}{\\log_{10}(d/a)} =
\\frac{\\log_{10}(n)}{\\log_{10}(d/L)+\\log_{10}(n)}
where :math:`n` is :math:`L` divided by :math:`a`.
Original code from the `mne-features <https://mne.tools/mne-features/>`_
package by Jean-Baptiste Schiratti and Alexandre Gramfort.
References
----------
* https://ieeexplore.ieee.org/abstract/document/904882
* https://hal.inria.fr/inria-00442374/
* https://www.hindawi.com/journals/ddns/2011/724697/
Examples
--------
>>> import numpy as np
>>> import entropy as ent
>>> import stochastic.processes.noise as sn
>>> rng = np.random.default_rng(seed=42)
>>> x = sn.FractionalGaussianNoise(hurst=0.5, rng=rng).sample(10000)
>>> print(f"{ent.katz_fd(x):.4f}")
6.4713
Fractional Gaussian noise with H = 0.9
>>> rng = np.random.default_rng(seed=42)
>>> x = sn.FractionalGaussianNoise(hurst=0.9, rng=rng).sample(10000)
>>> print(f"{ent.katz_fd(x):.4f}")
4.5720
Fractional Gaussian noise with H = 0.1
>>> rng = np.random.default_rng(seed=42)
>>> x = sn.FractionalGaussianNoise(hurst=0.1, rng=rng).sample(10000)
>>> print(f"{ent.katz_fd(x):.4f}")
7.6540
Random
>>> rng = np.random.default_rng(seed=42)
>>> print(f"{ent.katz_fd(rng.random(1000)):.4f}")
8.1531
Pure sine wave
>>> x = np.sin(2 * np.pi * 1 * np.arange(3000) / 100)
>>> print(f"{ent.katz_fd(x):.4f}")
2.4871
Linearly-increasing time-series (should be 1)
>>> x = np.arange(1000)
>>> print(f"{ent.katz_fd(x):.4f}")
1.0000
"""
x = np.asarray(x)
dists = np.abs(np.diff(x, axis=axis))
ll = dists.sum(axis=axis)
ln = np.log10(ll / dists.mean(axis=axis))
aux_d = x - np.take(x, indices=[0], axis=axis)
d = np.max(np.abs(aux_d), axis=axis)
kfd = np.squeeze(ln / (ln + np.log10(d / ll)))
if not kfd.ndim:
kfd = kfd.item()
return kfd
@jit('float64(float64[:], int32)')
def _higuchi_fd(x, kmax):
"""Utility function for `higuchi_fd`.
"""
n_times = x.size
lk = np.empty(kmax)
x_reg = np.empty(kmax)
y_reg = np.empty(kmax)
for k in range(1, kmax + 1):
lm = np.empty((k,))
for m in range(k):
ll = 0
n_max = floor((n_times - m - 1) / k)
n_max = int(n_max)
for j in range(1, n_max):
ll += abs(x[m + j * k] - x[m + (j - 1) * k])
ll /= k
ll *= (n_times - 1) / (k * n_max)
lm[m] = ll
# Mean of lm
m_lm = 0
for m in range(k):
m_lm += lm[m]
m_lm /= k
lk[k - 1] = m_lm
x_reg[k - 1] = log(1. / k)
y_reg[k - 1] = log(m_lm)
higuchi, _ = _linear_regression(x_reg, y_reg)
return higuchi
[docs]def higuchi_fd(x, kmax=10):
"""Higuchi Fractal Dimension.
Parameters
----------
x : list or np.array
One dimensional time series.
kmax : int
Maximum delay/offset (in number of samples).
Returns
-------
hfd : float
Higuchi fractal dimension.
Notes
-----
Original code from the `mne-features <https://mne.tools/mne-features/>`_
package by Jean-Baptiste Schiratti and Alexandre Gramfort.
This function uses Numba to speed up the computation.
References
----------
Higuchi, Tomoyuki. "Approach to an irregular time series on the
basis of the fractal theory." Physica D: Nonlinear Phenomena 31.2
(1988): 277-283.
Examples
--------
>>> import numpy as np
>>> import entropy as ent
>>> import stochastic.processes.noise as sn
>>> rng = np.random.default_rng(seed=42)
>>> x = sn.FractionalGaussianNoise(hurst=0.5, rng=rng).sample(10000)
>>> print(f"{ent.higuchi_fd(x):.4f}")
1.9983
Fractional Gaussian noise with H = 0.9
>>> rng = np.random.default_rng(seed=42)
>>> x = sn.FractionalGaussianNoise(hurst=0.9, rng=rng).sample(10000)
>>> print(f"{ent.higuchi_fd(x):.4f}")
1.8517
Fractional Gaussian noise with H = 0.1
>>> rng = np.random.default_rng(seed=42)
>>> x = sn.FractionalGaussianNoise(hurst=0.1, rng=rng).sample(10000)
>>> print(f"{ent.higuchi_fd(x):.4f}")
2.0581
Random
>>> rng = np.random.default_rng(seed=42)
>>> print(f"{ent.higuchi_fd(rng.random(1000)):.4f}")
2.0013
Pure sine wave
>>> x = np.sin(2 * np.pi * 1 * np.arange(3000) / 100)
>>> print(f"{ent.higuchi_fd(x):.4f}")
1.0091
Linearly-increasing time-series
>>> x = np.arange(1000)
>>> print(f"{ent.higuchi_fd(x):.4f}")
1.0040
"""
x = np.asarray(x, dtype=np.float64)
kmax = int(kmax)
return _higuchi_fd(x, kmax)
@jit('f8(f8[:])', nopython=True)
def _dfa(x):
"""
Utility function for detrended fluctuation analysis
"""
N = len(x)
nvals = _log_n(4, 0.1 * N, 1.2)
walk = np.cumsum(x - x.mean())
fluctuations = np.zeros(len(nvals))
for i_n, n in enumerate(nvals):
d = np.reshape(walk[:N - (N % n)], (N // n, n))
ran_n = np.array([float(na) for na in range(n)])
d_len = len(d)
trend = np.empty((d_len, ran_n.size))
for i in range(d_len):
slope, intercept = _linear_regression(ran_n, d[i])
trend[i, :] = intercept + slope * ran_n
# Calculate root mean squares of walks in d around trend
# Note that np.mean on specific axis is not supported by Numba
flucs = np.sum((d - trend) ** 2, axis=1) / n
# https://github.com/neuropsychology/NeuroKit/issues/206
fluctuations[i_n] = np.sqrt(np.mean(flucs))
# Filter zero
nonzero = np.nonzero(fluctuations)[0]
fluctuations = fluctuations[nonzero]
nvals = nvals[nonzero]
if len(fluctuations) == 0:
# all fluctuations are zero => we cannot fit a line
dfa = np.nan
else:
dfa, _ = _linear_regression(np.log(nvals), np.log(fluctuations))
return dfa
[docs]def detrended_fluctuation(x):
"""
Detrended fluctuation analysis (DFA).
Parameters
----------
x : list or np.array
One-dimensional time-series.
Returns
-------
alpha : float
the estimate alpha (:math:`\\alpha`) for the Hurst parameter.
:math:`\\alpha < 1`` indicates a
stationary process similar to fractional Gaussian noise with
:math:`H = \\alpha`.
:math:`\\alpha > 1`` indicates a non-stationary process similar to
fractional Brownian motion with :math:`H = \\alpha - 1`
Notes
-----
`Detrended fluctuation analysis
<https://en.wikipedia.org/wiki/Detrended_fluctuation_analysis>`_
is used to find long-term statistical dependencies in time series.
The idea behind DFA originates from the definition of self-affine
processes. A process :math:`X` is said to be self-affine if the standard
deviation of the values within a window of length n changes with the window
length factor :math:`L` in a power law:
.. math:: \\text{std}(X, L * n) = L^H * \\text{std}(X, n)
where :math:`\\text{std}(X, k)` is the standard deviation of the process
:math:`X` calculated over windows of size :math:`k`. In this equation,
:math:`H` is called the Hurst parameter, which behaves indeed very similar
to the Hurst exponent.
For more details, please refer to the excellent documentation of the
`nolds <https://cschoel.github.io/nolds/>`_
Python package by Christopher Scholzel, from which this function is taken:
https://cschoel.github.io/nolds/nolds.html#detrended-fluctuation-analysis
Note that the default subseries size is set to
entropy.utils._log_n(4, 0.1 * len(x), 1.2)). The current implementation
does not allow to manually specify the subseries size or use overlapping
windows.
The code is a faster (Numba) adaptation of the original code by Christopher
Scholzel.
References
----------
* C.-K. Peng, S. V. Buldyrev, S. Havlin, M. Simons,
H. E. Stanley, and A. L. Goldberger, “Mosaic organization of
DNA nucleotides,” Physical Review E, vol. 49, no. 2, 1994.
* R. Hardstone, S.-S. Poil, G. Schiavone, R. Jansen,
V. V. Nikulin, H. D. Mansvelder, and K. Linkenkaer-Hansen,
“Detrended fluctuation analysis: A scale-free view on neuronal
oscillations,” Frontiers in Physiology, vol. 30, 2012.
Examples
--------
Fractional Gaussian noise with H = 0.5
>>> import numpy as np
>>> import entropy as ent
>>> import stochastic.processes.noise as sn
>>> rng = np.random.default_rng(seed=42)
>>> x = sn.FractionalGaussianNoise(hurst=0.5, rng=rng).sample(10000)
>>> print(f"{ent.detrended_fluctuation(x):.4f}")
0.5216
Fractional Gaussian noise with H = 0.9
>>> rng = np.random.default_rng(seed=42)
>>> x = sn.FractionalGaussianNoise(hurst=0.9, rng=rng).sample(10000)
>>> print(f"{ent.detrended_fluctuation(x):.4f}")
0.8833
Fractional Gaussian noise with H = 0.1
>>> rng = np.random.default_rng(seed=42)
>>> x = sn.FractionalGaussianNoise(hurst=0.1, rng=rng).sample(10000)
>>> print(f"{ent.detrended_fluctuation(x):.4f}")
0.1262
Random
>>> rng = np.random.default_rng(seed=42)
>>> print(f"{ent.detrended_fluctuation(rng.random(1000)):.4f}")
0.5276
Pure sine wave
>>> x = np.sin(2 * np.pi * 1 * np.arange(3000) / 100)
>>> print(f"{ent.detrended_fluctuation(x):.4f}")
1.5848
Linearly-increasing time-series
>>> x = np.arange(1000)
>>> print(f"{ent.detrended_fluctuation(x):.4f}")
2.0390
"""
x = np.asarray(x, dtype=np.float64)
return _dfa(x)
