diff --git a/.travis.yml b/.travis.yml index 030de14..476a426 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,10 +1,17 @@ language: python python: -# - "3.4" -# - "3.5" + - "3.4" + - "3.5" - "3.6" - "3.7" -# - "nightly" + - "3.8" + - "3.9" + - "nightly" + +arch: + - amd64 + - arm64 + # command to install dependencies before_install: @@ -14,9 +21,25 @@ install: - pip install coverage - pip install -r requirements.txt -# command to run tests +jobs: + include: + - stage: "Extra packages" + name: "Testing extra packages" + python: 3.6 + install: + - pip install coverage + - pip install -r requirements.txt + - pip install EMD-signal + script: + - coverage run -m pytest test/test_EMD.py + - coverage run -m pytest test/test_MFDFA_extras.py + allow_failures: + - python: "nightly" + + script: - coverage run -m pytest test/ + - coverage run -m pytest test/test_fgn.py + - coverage run -m pytest test/test_MFDFA.py after_success: bash <(curl -s https://codecov.io/bash) diff --git a/LICENSE b/LICENSE index 2bfa978..490b4e9 100644 --- a/LICENSE +++ b/LICENSE @@ -1,6 +1,6 @@ MIT License -Copyright (c) 2019 Rydin +Copyright (c) 2019-2021 Rydin Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal diff --git a/MFDFA/MFDFA.py b/MFDFA/MFDFA.py index c122809..ac4c223 100644 --- a/MFDFA/MFDFA.py +++ b/MFDFA/MFDFA.py @@ -11,9 +11,9 @@ def MFDFA(timeseries: np.ndarray, lag: np.ndarray=None, order: int=1, q: np.ndarray=2, stat: bool=False, modified: bool=False, - extensions: dict={"EMD":False, "eDFA":False}) -> np.ndarray: + extensions: dict={'EMD':False, 'eDFA':False, 'window':False}) -> np.ndarray: """ - Multi-Fractal Detrended Fluctuation Analysis of timeseries. MFDFA generates + Multifractal Detrended Fluctuation Analysis of timeseries. MFDFA generates a fluctuation function F²(q,s), with s the segment size and q the q-powers, Take a timeseries Xₜ, find the integral Yₜ = cumsum(Xₜ), and segment the timeseries into Nₛ segments of size s. @@ -86,6 +86,11 @@ def MFDFA(timeseries: np.ndarray, lag: np.ndarray=None, order: int=1, - ``eDFA``: bool (default ``False``) A method to evaluate the strength of multifractality. Calls function `eDFA()`. + - ``window``: bool (default ``False``) + A moving window for smaller timeseries. Set ``window`` as int > 0 with + the number of steps the window shoud move over the data. ``window = 1`` + will move window by ``1`` step. Since the timeseries is segmented at + each lag lenght, any window choise > lag is only segmented once. Returns ------- @@ -117,9 +122,19 @@ def MFDFA(timeseries: np.ndarray, lag: np.ndarray=None, order: int=1, assert timeseries.shape[1] == 1, "Timeseries needs to be 1 dimensional" timeseries = timeseries.reshape(-1,1) + # Size of array N = timeseries.shape[0] + # Assert if window is given, that it is int and > 0 + window = False + if 'window' in extensions: + if extensions['window'] != False: + window = extensions['window'] + assert isinstance(window, int), "'window' is not integer" + assert window > 0, "'window' is not > 0" + + # Fractal powers as floats q = np.asarray_chkfinite(q, dtype = float) @@ -129,7 +144,7 @@ def MFDFA(timeseries: np.ndarray, lag: np.ndarray=None, order: int=1, # Reshape q to perform np.float_power q = q.reshape(-1, 1) - # x-axis + # x-axis needed for polyfit X = np.linspace(1, lag.max(), lag.max()) # "Profile" of the series @@ -145,7 +160,7 @@ def MFDFA(timeseries: np.ndarray, lag: np.ndarray=None, order: int=1, if stat == True: f_std = np.empty((0, q.size)) - + # Check which extensions are requested if ('eDFA', True) in extensions.items(): f_eDFA = np.empty((0, q.size)) @@ -155,7 +170,7 @@ def MFDFA(timeseries: np.ndarray, lag: np.ndarray=None, order: int=1, assert isinstance(extensions['EMD'], list), 'list IMFs to detrend' # Detrending of the timeseries using EMD with given IMFs in a list - Y = detrendedtimeseries(Y, extensions["EMD"]) + Y = detrendedtimeseries(Y, extensions['EMD']) # Force order = 0 since the data is detrended with EMD, i.e., no # need to do polynomial fittings anymore @@ -164,32 +179,63 @@ def MFDFA(timeseries: np.ndarray, lag: np.ndarray=None, order: int=1, # Loop over elements in lag # Notice that given one has to split the timeseries into different segments # of length lag(), some elements at the end of the array might be - # missing. The same procedure is run in reverse, where elements at the - # begining of the series are discarded instead. + # missing. The same procedure is run in reverse—if not using an moving + # window—where elements at the begining of the series are discarded instead. for i in lag: - # Reshape into (N/lag, lag) - Y_ = Y[:N - N % i].reshape((N - N % i) // i, i) - Y_r = Y[N % i:].reshape((N - N % i) // i, i) - - # If order = 0 one gets simply a Fluctuation Analysis (FA), or one is - # using the EMD setting, and the data is detrended. - if order == 0: - # Skip detrending - F = np.append( np.var(Y_, axis=1), np.var(Y_r, axis=1) ) - - else: - # Perform a polynomial fit to each segments - p = polyfit(X[:i], Y_.T, order) - p_r = polyfit(X[:i], Y_r.T, order) - - # Subtract the trend from the fit and calculate the variance - F = np.append( - np.var(Y_ - polyval(X[:i], p), axis = 1), - np.var(Y_r - polyval(X[:i], p_r), axis = 1) - ) - - # Caculate the Multi-Fractal (Non)-Detrended Fluctuation Analysis + # Standard option + if window == False: + # Reshape into (N/lag, lag) + Y_ = Y[:N - N % i].reshape((N - N % i) // i, i) + Y_r = Y[N % i:].reshape((N - N % i) // i, i) + + # If order = 0 one gets simply Fluctuation Analysis (FA), or if one is + # using the EMD setting the data is detrended and no polynomial fitting + # is needed. + if order == 0: + # Skip detrending + F = np.append(np.var(Y_, axis=1), np.var(Y_r, axis=1)) + + else: + # Perform a polynomial fit to each segments + p = polyfit(X[:i], Y_.T, order) + p_r = polyfit(X[:i], Y_r.T, order) + + # Subtract the trend from the fit and calculate the variance + F = np.append( + np.var(Y_ - polyval(X[:i], p), axis = 1), + np.var(Y_r - polyval(X[:i], p_r), axis = 1) + ) + + # For short timeseries, using a moving window instead of segmenting the + # timeseries. Notice the number of operations is considerably larger + # depending on the moving window displacement. + if window != False: + + F = np.empty(0) + for j in range(0, i-1, window): + + # subtract j points as the moving window shortens the data + N_0 = N - j + + # Reshape into (N_0/lag, lag) + Y_ = Y[j:N - N_0 % i].reshape((N - N_0 % i) // i, i) + + # If order = 0 one gets simply Fluctuation Analysis (FA), or if one + # is using the EMD setting the data is detrended and no polynomial + # fitting is needed. + if order == 0: + # Skip detrending + F = np.append(F, np.var(Y_, axis=1)) + + else: + # Perform a polynomial fit to each segments + p = polyfit(X[:i], Y_.T, order) + + # Subtract the trend from the fit and calculate the variance + F = np.append(F, np.var(Y_ - polyval(X[:i], p), axis = 1)) + + # Caculate the Multifractal (Non)-Detrended Fluctuation Analysis f = np.append(f, np.float_power( np.mean(np.float_power(F, q / 2), axis = 1), diff --git a/MFDFA/emddetrender.py b/MFDFA/emddetrender.py index 35e979d..ea8ff1b 100644 --- a/MFDFA/emddetrender.py +++ b/MFDFA/emddetrender.py @@ -2,9 +2,11 @@ # Empirical Mode Decomposition algorithm', https://github.com/laszukdawid/PyEMD, # licenced under the Apache 2.0 Licencing. -from PyEMD import EMD +# from PyEMD import EMD import numpy as np +# Import of PyEMD is called inside the function + def detrendedtimeseries(timeseries: np.ndarray, modes: list): """ The function calculates the Intrinsic Mode Functions (IMFs) of a given @@ -50,6 +52,7 @@ def detrendedtimeseries(timeseries: np.ndarray, modes: list): return detrendedTimeseries + def IMFs(timeseries: np.ndarray): """ Extract the Intrinsic Mode Functions (IMFs) of a given timeseries. @@ -62,7 +65,7 @@ def IMFs(timeseries: np.ndarray): Notes ----- .. versionadded:: 0.3 - + Returns ------- IMFs: np.ndarray @@ -70,6 +73,11 @@ def IMFs(timeseries: np.ndarray): These are of shape ``(..., timeseries.size)``, with the first dimension varying depending on the data. Last entry is the residuals. """ + + # Check if EMD-signal is installed + missing_library() + from PyEMD import EMD + # Initiate pyEMD's EMD function emd = EMD() @@ -78,3 +86,13 @@ def IMFs(timeseries: np.ndarray): # Returns the IMFs as a (..., timeseries.size) numpy array. return IMFs + + +def missing_library(): + try: + import PyEMD.EMD as _EMD + except ImportError: + raise ImportError(("PyEMD is required to do Empirical Mode " + "decomposition. Please install PyEMD with 'pip " + "install EMD-signal'.") + ) diff --git a/README.md b/README.md index 389336d..34e7ac1 100644 --- a/README.md +++ b/README.md @@ -7,17 +7,12 @@ [![Documentation Status](https://readthedocs.org/projects/mfdfa/badge/?version=latest)](https://mfdfa.readthedocs.io/en/latest/?badge=latest) -### Preparation for version 0.4 -I have been deliberating over the presence of the Empirical Mode Decomposition [EMD](https://en.wikipedia.org/wiki/Hilbert%E2%80%93Huang_transform) in v0.3, i.e., the [PyEMD](https://github.com/laszukdawid/PyEMD) by [Dawid Laszuk](https://github.com/laszukdawid). -The EMD package is phenomenal, but comprises an overhead of ~700 kb (and usually some additional packages). -My plan now is to create a separate repository (likely to be named MFDFA-EMD) which adds the EMD functionality included in v0.3. - -The cause is simple: the original MFDFA is meant to be lightweight, i.e., depend only of `numpy`, so that one day maybe someone can make a parallel-isable version of it. - # MFDFA Multifractal Detrended Fluctuation Analysis `MFDFA` is a model-independent method to uncover the self-similarity of a stochastic process or auto-regressive model. `DFA` was first developed by Peng *et al.*1 and later extended to study multifractality `MFDFA` by Kandelhardt *et al.*2. +In the latest release there is as well added a moving window system, especially useful for short timeseries, a recent extension to DFA called *extended DFA*, and the extra feature of Empirical Mode Decomposition as detrending method. + # Installation To install MFDFA you can simply use @@ -36,6 +31,8 @@ There is an added library `fgn` to generate fractional Gaussian noise. # Employing the `MFDFA` library ## An exemplary one-dimensional fractional Ornstein–Uhlenbeck process +The rationale here is simple: Numerically integrate a stochastic process in which we know exactly the fractal properties, characterised by the Hurst coefficient, and recover this with MFDFA. +We will use a fractional Ornstein–Uhlenbeck, a commonly employ stochastic process with mean-reverting properties. For a more detailed explanation on how to integrate an Ornstein–Uhlenbeck process, see the [kramersmoyal's package](https://github.com/LRydin/KramersMoyal#a-one-dimensional-stochastic-process). You can also follow the [fOU.ipynb](/examples/fOU.ipynb) @@ -77,7 +74,7 @@ To now utilise the `MFDFA`, we take this exemplary process and run the (multifra ```python # Select a band of lags, which usually ranges from # very small segments of data, to very long ones, as -lag = np.logspace(0.7, 4, 30).astype(int) +lag = np.unique(np.logspace(0.5, 3, 100).astype(int)) # Notice these must be ints, since these will segment # the data into chucks of lag size @@ -105,15 +102,14 @@ np.polyfit(np.log(lag[:15]), np.log(dfa[:15]),1)[0] # Now what you should obtain is: slope = H + 1 ``` - + ## Uncovering multifractality in stochastic processes -`MFDFA`, as an extension to `DFA`, was developed to uncover if a given process is mono or multi fractal. -Let `Xₜ` be a multi fractal stochastic process. This mean `Xₜ` scales with some function alpha(t) as `Xcₜ = |c|alpha(t) Xₜ`. -With the help of taking different powers variations of the `DFA`, one we can distinguish monofractal and multifractal processes. +You can find more about multifractality in the [documentation](https://mfdfa.readthedocs.io/en/latest/1dLevy.html). # Changelog +- Version 0.4 - EMD is now optional. Restored back compatibility: py3.3 to py3.9. For EMD py3.6 or larger is needed. - Version 0.3 - Adding EMD detrending. First release. PyPI code. - Version 0.2 - Removed experimental features. Added documentation - Version 0.1 - Uploaded initial working code @@ -128,11 +124,11 @@ This package abides to a [Conduct of Fairness](contributions.md). This library has been submitted for publication at [The Journal of Open Source Software](https://joss.theoj.org/) in December 2019. It was rejected. The review process can be found [here on GitHub](https://github.com/openjournals/joss-reviews/issues/1966). The plan is to extend the library and find another publisher. ### History -This project was started in 2019 at the [Faculty of Mathematics, University of Oslo](https://www.mn.uio.no/math/english/research/groups/risk-stochastics/) in the Risk and Stochastics section by Leonardo Rydin Gorjão and is supported by Dirk Witthaut and the [Institute of Energy and Climate Research Systems Analysis and Technology Evaluation](https://www.fz-juelich.de/iek/iek-ste/EN/Home/home_node.html). I'm very thankful to all the folk in Section 3 in the Faculty of Mathematics, University of Oslo, for helping me getting around the world of stochastic processes: Dennis, Anton, Michele, Fabian, Marc, Prof. Benth and Prof. di Nunno. In April 2020 Galib Hassan joined in extending `MFDFA`, particularly the implementation of EMD. +This project was started in 2019 at the [Faculty of Mathematics, University of Oslo](https://www.mn.uio.no/math/english/research/groups/risk-stochastics/) in the Risk and Stochastics section by Leonardo Rydin Gorjão and is supported by Dirk Witthaut and the [Institute of Energy and Climate Research Systems Analysis and Technology Evaluation](https://www.fz-juelich.de/iek/iek-ste/EN/Home/home_node.html). I'm very thankful to all the folk in Section 3 in the Faculty of Mathematics, University of Oslo, for helping me getting around the world of stochastic processes: Dennis, Anton, Michele, Fabian, Marc, Prof. Benth and Prof. di Nunno. In April 2020 Galib Hassan joined in extending `MFDFA`, particularly the implementation of `EMD`. ### Funding -Helmholtz Association Initiative *Energy System 2050 - A Contribution of the Research Field Energy* and the grant No. VH-NG-1025, *STORM - Stochastics for Time-Space Risk Models* project of the Research Council of Norway (RCN) No. 274410, and the *E-ON Stipendienfonds*. +Helmholtz Association Initiative *Energy System 2050 - A Contribution of the Research Field Energy* and the grant No. VH-NG-1025; *STORM - Stochastics for Time-Space Risk Models* project of the Research Council of Norway (RCN) No. 274410, and the *E-ON Stipendienfonds*. ### References 1Peng, C.-K., Buldyrev, S. V., Havlin, S., Simons, M., Stanley, H. E., & Goldberger, A. L. (1994). *Mosaic organization of DNA nucleotides*. [Physical Review E, 49(2), 1685–1689](https://doi.org/10.1103/PhysRevE.49.1685)\ diff --git a/contributions.md b/contributions.md index 2765e1b..5e72fbf 100644 --- a/contributions.md +++ b/contributions.md @@ -1,5 +1,5 @@ # Contributions -`MDFA` has lived so far out of collaborations. It welcomes reviews and ideas from everyone. If you want to share your ideas or report a bug, open an [issue](https://github.com/LRydin/KramersMoyal/issues) here on GitHub, or contact me directly: [leonardo.rydin"at"gmail.com](mailto:leonardo.rydin@gmail.com) +`MFDFA` has lived so far out of collaborations. It welcomes reviews and ideas from everyone. If you want to share your ideas or report a bug, open an [issue](https://github.com/LRydin/KramersMoyal/issues) here on GitHub, or contact me directly: [leonardo.rydin"at"gmail.com](mailto:leonardo.rydin@gmail.com) I'm are happy to have you help me in furthering this project, helping me improve the code and documentation, make it more user-friendly, broadening its applicability, or devising other methodologies of implementation. diff --git a/docs/1dLevy.rst b/docs/1dLevy.rst new file mode 100644 index 0000000..6b1f901 --- /dev/null +++ b/docs/1dLevy.rst @@ -0,0 +1,47 @@ +Multifractality in one dimensional distributions +================================================ + +To show how multifractality can be studied, let us take a sample of random numbers of a symmetric `Lévy distribution `_. + + +Univariate random numbers from a Lévy stable distribution +--------------------------------------------------------- + +To obtain a sample of random numbers of Lévy stable distributions, use `scipy`'s `levy_stable`. In particular, take an :math:`\alpha`-stable distribution, with :math:`\alpha=1.5` + +.. code:: python + + # Imports + from MFDFA import MFDFA + from scipy.stats import levy_stable + + # Generate 100000 points + alpha = 1.5 + X = levy_stable.rvs(alpha=alpha, beta = 0, size=10000) + + +For :code:`MFDFA` to detect the multifractal spectrum of the data, we need to vary the parameter :math:`q\in[-10,10]` and exclude :math:`0`. Let us also use a quadratic polynomial fitting by setting :code:`order=2` + +.. code:: python + + # Select a band of lags, which are ints + lag = np.unique(np.logspace(0.5, 3, 100).astype(int)) + + # Select a list of powers q + q_list = np.linspace(-10,10,41) + q_list = q_list[q_list!=0.0] + + # The order of the polynomial fitting + order = 2 + + # Obtain the (MF)DFA as + lag, dfa = MFDFA(y, lag = lag, q = q_list, order = order) + + +Again, we plot this in a double logarithmic scale, but now we include 6 curves, from 6 selected :math:`q={-10,-5-2,2,5,10}`. Include as well are the theoretical curves for :math:`q=-10`, with a slope of :math:`1/\alpha=1/1.5` and :math:`q=10`, with a slope of :math:`1/q=1/10` + + +.. image:: /_static/fig2.png + :height: 450 + :align: center + :alt: Plot of the MFDFA of a Lévy stable distribution for a few q values. diff --git a/docs/1dprocess.rst b/docs/1dprocess.rst index fa02843..1ea3271 100644 --- a/docs/1dprocess.rst +++ b/docs/1dprocess.rst @@ -49,7 +49,7 @@ To now utilise the :code:`MFDFA`, we take this exemplary process and run the (mu # Select a band of lags, which usually ranges from # very small segments of data, to very long ones, as - lag = np.logspace(0.7, 4, 30).astype(int) + lag = np.unique(np.logspace(0.5, 3, 100).astype(int)) # Notice these must be ints, since these will segment # the data into chucks of lag size diff --git a/docs/_static/fig1.png b/docs/_static/fig1.png index 5e0e01a..d3c6975 100644 Binary files a/docs/_static/fig1.png and b/docs/_static/fig1.png differ diff --git a/docs/_static/fig2.png b/docs/_static/fig2.png new file mode 100644 index 0000000..fbc39eb Binary files /dev/null and b/docs/_static/fig2.png differ diff --git a/docs/conf.py b/docs/conf.py index 94e21e6..e84e6a7 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -22,12 +22,12 @@ # -- Project information ----------------------------------------------------- project = 'MFDFA' -copyright = '2019-2020 Leonardo Rydin' +copyright = '2019-2021 Leonardo Rydin' author = 'Leonardo Rydin' # The full version, including alpha/beta/rc tags -release = '0.3' -version = '0.3' +release = '0.4' +version = '0.4' # -- General configuration --------------------------------------------------- @@ -89,7 +89,7 @@ 'style_external_links': False, 'style_nav_header_background': 'black', # Toc options - 'collapse_navigation': True, + 'collapse_navigation': False, 'sticky_navigation': True, 'navigation_depth': 4, 'titles_only': False diff --git a/docs/emd_detrending.rst b/docs/emd_detrending.rst index 305b2bb..f78e077 100644 --- a/docs/emd_detrending.rst +++ b/docs/emd_detrending.rst @@ -4,6 +4,14 @@ Employing Empirical Mode Decompositions for detrending `Empirical Mode Decomposition `_ (EMD), or maybe more correctly described, the Hilbert─Huang transform is a transformation analogous to a Fourier or Hilbert transform that decomposes a one-dimensional timeseries or signal into its Intrinsic Mode Functions (IMFs). For our purposes, we simply want to employ EMD to detrend a timeseries. +.. warning:: + + To use this feature, you need to first install `PyEMD `_ (EMD-signal) with + + :: + + pip install EMD-signal + Understanding :code:`MFDFA`'s :code:`EMD` detrender ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/docs/extendedDFA.rst b/docs/extendedDFA.rst new file mode 100644 index 0000000..b859e35 --- /dev/null +++ b/docs/extendedDFA.rst @@ -0,0 +1,27 @@ +Extended Detrended Fluctuation Analysis +--------------------------------------- + +In the publication `Detrended fluctuation analysis of cerebrovascular responses to abrupt changes in peripheral arterial pressure in rats. `_ the authors introduce a new metric similar to the conventional Detrended Fluctuation Analysis (DFA) which they denote *Extended* Detrended Fluctuation Analysis (eDFA), which relies on extracting the difference of the minima and maxima for each segmentation of the data, granting a new power-law exponent to study, i.e., as in eq. (5) in the paper + +.. math:: + \mathrm{d}F (n) = \mathrm{max}[F(n)] - \mathrm{min}[F(n)], + +which in turn results in + +.. math:: + \mathrm{d}F(n) \sim n^\beta. + + +Using :code:`MFDFA`'s :code:`eDFA` extension +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +To obtain the :code:`eDFA`, simply set the extension to :code:`True` and add a new output function, here denoted :code:`edfa` + +.. code:: python + + # Select a band of lags, which usually ranges from + # very small segments of data, to very long ones, as + lag = np.logspace(0.7, 4, 30).astype(int) + + # Obtain the (MF)DFA by declaring the IMFs to subtract + # in a list in the dictionary of the extensions + lag, dfa, edfa = MFDFA(y, lag = lag, extensions = {'eDFA': True}) diff --git a/docs/extensions.rst b/docs/extensions.rst index ce12552..3f778fe 100644 --- a/docs/extensions.rst +++ b/docs/extensions.rst @@ -4,3 +4,7 @@ Extensions MFDFA as seen since its development a set of enhancements. In particular the usage of Empirical Mode Decomposition as a source of detrending, instead of polynomial fittings, which allows for a more precise removal of known trends in the timeseries. .. include:: emd_detrending.rst + +.. include:: extendedDFA.rst + +.. include:: moving_window.rst diff --git a/docs/index.rst b/docs/index.rst index 1bafb25..e80e3bb 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -10,6 +10,8 @@ MFDFA's documentation .. include:: 1dprocess.rst +.. include:: 1dLevy.rst + .. include:: extensions.rst Literature @@ -30,9 +32,11 @@ Table of Content .. toctree:: :maxdepth: 3 + :hidden: installation 1dprocess + 1dLevy extensions functions/index license diff --git a/docs/installation.rst b/docs/installation.rst index 2964dec..9fc6da9 100644 --- a/docs/installation.rst +++ b/docs/installation.rst @@ -1,14 +1,21 @@ Installation ============ -For the moment the library is available from TestPyPI, so you can use +:code:`MFDFA` is available from PyPI, so you can use :: - pip install -i https://test.pypi.org/simple/ MFDFA + pip install MFDFA Then on your favourite editor just use .. code:: python from MFDFA import MFDFA + +.. warning:: + To use the extension to include Empirical Mode Decomposition detrending you will also need + + :: + + pip install EMD-signal diff --git a/docs/license.rst b/docs/license.rst index d107f45..2ea3045 100644 --- a/docs/license.rst +++ b/docs/license.rst @@ -3,7 +3,7 @@ License MIT License -Copyright (c) 2019-2020 Rydin +Copyright (c) 2019-2021 Rydin Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal diff --git a/docs/moving_window.rst b/docs/moving_window.rst new file mode 100644 index 0000000..776f13a --- /dev/null +++ b/docs/moving_window.rst @@ -0,0 +1,18 @@ +Moving window for segmentation +------------------------------ + +For short timeseries the segmentation of the data—especially for large lags—results in bad statistics, e.g. if a timeseries has `2048` datapoints and one wishes to study the flucutation analysis up to a lag of `512`, only 4 segmentations of the data are possible for the lag `512`. Instead one can use an moving window over the timeseries to obtain better statistics at large lags. + +Using :code:`MFDFA`'s :code:`window` extension +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +To utilise a moving window one has to declare the moving windows step-size, i.e., the number of data points the window will move over the data. Say we wish to increase the statistics of the aforementioned example to include a moving window moving `32` steps (so one has 64 segments at a lag of `512`) + +.. code:: python + + # Select a band of lags, which usually ranges from + # very small segments of data, to very long ones, as + lag = np.logspace(0.7, 4, 30).astype(int) + + # Obtain the (MF)DFA by declaring the IMFs to subtract + # in a list in the dictionary of the extensions + lag, dfa, edfa = MFDFA(y, lag = lag, extensions = {'window': 32}) diff --git a/other/fig1.png b/other/fig1.png deleted file mode 100644 index 5e0e01a..0000000 Binary files a/other/fig1.png and /dev/null differ diff --git a/requirements.txt b/requirements.txt index 64a39a7..24ce15a 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,2 +1 @@ numpy -EMD-signal diff --git a/setup.py b/setup.py index 80dedea..5bf7cad 100644 --- a/setup.py +++ b/setup.py @@ -5,7 +5,7 @@ setuptools.setup( name="MFDFA", - version="0.3", + version="0.4", author="Leonardo Rydin Gorjao", author_email="leonardo.rydin@gmail.com", description="Multifractal Detrended Fluctuation Analysis in Python", @@ -13,10 +13,12 @@ long_description_content_type="text/markdown", url="https://github.com/LRydin/MFDFA", packages=setuptools.find_packages(), + install_requires = ["numpy"], + extras_require = {"EMD-signal": ["EMD-signal"]}, classifiers=[ "Programming Language :: Python :: 3", "License :: OSI Approved :: MIT License", "Operating System :: OS Independent", ], - python_requires='>=3.6', + python_requires='>=3.3', ) diff --git a/test/test_MFDFA.py b/test/test_MFDFA.py index 2bf7622..87d6c85 100644 --- a/test/test_MFDFA.py +++ b/test/test_MFDFA.py @@ -22,28 +22,38 @@ def test_MFDFA(): assert dfa.ndim == 2, "Output is not 2 dimensional" assert dfa.shape[1] <= q.shape[0], "Output shape mismatch" + # Testing conventional MFDFA with stats + lag, dfa, dfa_std = MFDFA(X, lag = lag, q = q, order = 2, + stat = True) + + assert dfa.ndim == 2, "Output is not 2 dimensional" + assert dfa.shape[1] <= q.shape[0], "Output shape mismatch" + + # Testing modified = True option lag, dfa, dfa_std = MFDFA(X, lag = lag, q = q, order = 2, modified = True, stat = True) assert dfa.ndim == 2, "Output is not 2 dimensional" assert dfa.shape[1] <= q.shape[0], "Output shape mismatch" + # Testing eDFA extension lag, dfa, edfa = MFDFA(X, lag = lag, q = q, order = 3, - modified = True, stat = False, extensions = {'eDFA':True}) + stat = False, extensions = {'eDFA':True}) assert dfa.ndim == 2, "Output is not 2 dimensional" assert edfa.ndim == 2, "Output is not 2 dimensional" assert dfa.shape[1] <= q.shape[0], "Output shape mismatch" - lag, dfa = MFDFA(X, lag = lag, q = q, order = 0, - modified = True, stat = False, extensions = {'EMD': [0]}) + # Testing moving window + lag, dfa = MFDFA(X, lag = lag, q = q, order = 1, + stat = False, extensions = {'window': 5}) assert dfa.ndim == 2, "Output is not 2 dimensional" assert dfa.shape[1] <= q.shape[0], "Output shape mismatch" - lag, dfa, edfa = MFDFA(X, lag = lag, q = q, order = 1, - modified = True, stat = False, - extensions = {'eDFA':True, 'EMD': [0]}) + # Testing moving window and eDFA with stats = True + lag, dfa, dfa_std, edfa = MFDFA(X, lag = lag, q = q, order = 1, + stat = True, extensions = {'eDFA':True, 'window': 5}) assert dfa.ndim == 2, "Output is not 2 dimensional" assert edfa.ndim == 2, "Output is not 2 dimensional" diff --git a/test/test_MFDFA_extras.py b/test/test_MFDFA_extras.py new file mode 100644 index 0000000..bedf2de --- /dev/null +++ b/test/test_MFDFA_extras.py @@ -0,0 +1,32 @@ +import numpy as np + +import sys +sys.path.append("../") +from MFDFA import MFDFA + +def test_MFDFA(): + for N in [1000, 10000]: + for q_list in [1, 2, 6, 21]: + + X = np.random.normal(size = N, loc = 0) + q = np.linspace(-10, 10, q_list) + + lag = np.unique( + np.logspace( + 0, np.log10(X.size // 4), 25 + ).astype(int) + 1 + ) + + lag, dfa = MFDFA(X, lag = lag, q = q, order = 0, + modified = True, stat = False, extensions = {'EMD': [0]}) + + assert dfa.ndim == 2, "Output is not 2 dimensional" + assert dfa.shape[1] <= q.shape[0], "Output shape mismatch" + + lag, dfa, edfa = MFDFA(X, lag = lag, q = q, order = 1, + modified = True, stat = False, + extensions = {'eDFA':True, 'EMD': [0]}) + + assert dfa.ndim == 2, "Output is not 2 dimensional" + assert edfa.ndim == 2, "Output is not 2 dimensional" + assert dfa.shape[1] <= q.shape[0], "Output shape mismatch" diff --git a/tester.py b/tester.py index b0d7e48..c205479 100644 --- a/tester.py +++ b/tester.py @@ -3,7 +3,8 @@ from numpy.polynomial.polynomial import polyfit import matplotlib.pyplot as plt -import MFDFA as MFDFA +from MFDFA import MFDFA +from MFDFA import fgn # %% Lets take a fractional Ornstein–Uhlenbeck process Xₜ with a time-dependent # diffusion or volatility σₜ, some drift of mean reverting term θ(Xₜ), and @@ -31,7 +32,7 @@ H = 0.5 # Generate the fractional Gaussian noise -dw = (t_final ** H) * MFDFA.fgn(time.size, H = H) +dw = (t_final ** H) * fgn(time.size, H = H) # Integrate the process for i in range(1,time.size): @@ -49,7 +50,7 @@ q = 2 # dfa records the fluctuation function using the EMD as a detrending mechanims. -lag, dfa, dfa_std, e_dfa = MFDFA.MFDFA(X, lag, q = q, order = 1, stat = True, extensions = {"eDFA": True}) +lag, dfa, dfa_std, e_dfa = MFDFA(X, lag, q = q, order = 1, stat = True, extensions = {"eDFA": True}) # %% Plots # Visualise the results in a log-log plot