Bivariate doc example and fix (#23)

* Bivariate doc example

README.rst

@@ -30,7 +30,7 @@ Implemented features:
 * Outlier detection.
 
 
-The initial implementation of the code in this package was based on the Wavelet p-Leader and Bootstrap based MultiFractal analysis (PLBMF) `Matlab toolbox <http://www.ens-lyon.fr/PHYSIQUE/Equipe3/MultiFracs/software.html>`_ written by Patrice Abry, Herwig Wendt and colleagues. For a thorough introduction to multifractal analysis, you may access H. Wendt's PhD thesis available in `his webiste <https://www.irit.fr/~Herwig.Wendt/data/ThesisWendt.pdf)>`_.
+The initial implementation of the code in this package was based on the Wavelet p-Leader and Bootstrap based MultiFractal analysis (PLBMF) `Matlab toolbox <http://www.ens-lyon.fr/PHYSIQUE/Equipe3/MultiFracs/software.html>`_ written by Patrice Abry, Herwig Wendt and colleagues. For a thorough introduction to multifractal analysis, you may access H. Wendt's PhD thesis available in `his website <https://www.irit.fr/~Herwig.Wendt/data/ThesisWendt.pdf)>`_.
 
 
 .. For a brief introduction to multifractal analysis, see the file THEORY.ipynb

examples/01-simul.py

@@ -6,15 +6,14 @@
 ==========================
 """
 
-# pylint: disable=C0413
-# flake8: disabel=E402
-
 # %% [markdown]
 # The ``pymultifracs`` package estimates the multifractal properties of time series. This is an example using simulated Multifractal Random Walks which covers all steps of the multifractal analysis procedure and gives an overview of the toolbox's features.
 #
 # Let us first generate a few Multifractal Random Walks with parameters :math:`H = 0.8` and :math:`\lambda=\sqrt{0.05}`
 
 # %% [python]
+# pylint: disable=C0413
+# flake8: disabel=E402
 
 # sphinx_gallery_start_ignore
 import seaborn as sns

examples/10-bivariate.py

@@ -0,0 +1,86 @@
+"""
+.. _analysis_bivariate:
+
+==========================
+Analysis of bivariate data
+==========================
+"""
+
+# #############################################################################
+
+# The ``pymultifracs`` package can perform bivariate multifractal analysis,
+# and produces bivariate cumulants and structure functions. The purpose of
+# bivariate MFA is to determine the relationship between the regularity
+# fluctuations of a pair of time series.
+#
+# We will be using pre-generated data to serve as an example of bivariate
+# multi-fractal time series. Let us first load the data; those example time
+# series are also found in the ``tests/data/`` folder of the repository on
+# github.
+
+import pooch
+from scipy.io import loadmat
+
+url = ('https://github.com/neurospin/pymultifracs/raw/refs/heads/master/tests/'
+       'data/DataSet_ssMF.mat')
+
+fname = pooch.retrieve(url=url, known_hash=None, path=pooch.os_cache("bivar"))
+
+X = loadmat(fname)['data'].T
+
+# %%
+# The data is a two-dimensional array, with both self-similarity and
+# multifractal correlation.
+
+# sphinx_gallery_start_ignore
+# pylint: disable=C0413
+# flake8: disable=E402
+import seaborn as sns
+import matplotlib as mpl
+mpl.rcParams['figure.dpi'] = 600
+# sns.set_theme(style="whitegrid")
+# sns.set_context('paper')
+# sphinx_gallery_end_ignore
+
+import matplotlib.pyplot as plt
+
+fig, ax = plt.subplots(2, 1, sharex=True)
+ax[0].plot(X[:, 0])
+ax[1].plot(X[:, 1], c='C1')
+
+# %%
+# First we need :func:`wavelet_analysis` as ususal, to obtain p-leaders in this case
+
+from pymultifracs import wavelet_analysis
+
+WTpL = wavelet_analysis(X).integrate(.75).get_leaders(2)
+
+# %%
+# Then we use :func:`.bivariate.bimfa` on the wavelet p-leaders. We provide
+# the same MRQ twice, in order to compute the estimates for all signal pairs.
+
+import numpy as np
+from pymultifracs.bivariate import bimfa
+
+q1 = np.array([0, 1, 2])
+q2 = np.array([0, 1, 2])
+
+pwt = bimfa(WTpL, WTpL, scaling_ranges=[(3, 9)], q1=q1, q2=q2)
+
+# %%
+
+# We can obtain the multifractal correlation matrix :math:`\rho_{mf}`:
+
+print(pwt.cumulants.rho_mf.squeeze())
+
+# %%
+# We can plot the structure function:
+pwt.structure.plot()
+
+# %%
+# We can plot the bivariate cumulants:
+pwt.cumulants.plot()
+
+# %%
+# We can plot the Legendre spectrum reconstituted from the log-cumulants:
+pwt.cumulants.plot_legendre(h_support=(.1, .95))

pymultifracs/bivariate/biscaling_function.py

@@ -7,6 +7,7 @@
 
 import numpy as np
 import matplotlib.pyplot as plt
+import matplotlib as mpl
 from scipy import special
 
 
@@ -312,7 +313,7 @@ def plot(self, figsize=None, scaling_range=0, signal_idx1=0,
         idx = np.s_[j_min:j_max]
 
         _, axes = plt.subplots(len(self.q1), len(self.q2), sharex=True,
-                               figsize=figsize)
+                               figsize=figsize, squeeze=False)
 
         x = self.j[idx]
 
@@ -343,7 +344,7 @@ def plot(self, figsize=None, scaling_range=0, signal_idx1=0,
                 ax.errorbar(x, y, CI)
                 # ax.tick_params(bottom=False, top=False, which='minor')
                 ax.set(xlabel='Temporal scale $j$',
-                       ylabel=f'$q_1={q1:.1f}$, $q_2={q2:.1f}$')
+                       ylabel=f'$q_1={q1:.1g}$, $q_2={q2:.1f}$')
 
                 x0, x1 = self.scaling_ranges[scaling_range]
 
@@ -407,7 +408,6 @@ def __post_init__(self, idx_reject, mrq1, mrq2, min_j):
 
         self._compute(mrq1, mrq2, idx_reject)
 
-
         self.slope = np.zeros(
             (self.n_cumul+1, self.n_cumul+1, len(self.scaling_ranges),
              *self.values.shape[4:]))
@@ -422,23 +422,19 @@ def __post_init__(self, idx_reject, mrq1, mrq2, min_j):
         self.slope[idx_margin1] = slope1[:, 0, :, :, None]
         self.intercept[idx_margin1] = intercept1[:, 0, :, :, None]
 
-
         slope2, intercept2, _ = self._compute_fit(
             mrq1, mrq2, margin=1, value_name='margin2_values')
         self.slope[idx_margin2] = slope2[:, 0, :, None]
         self.intercept[idx_margin2] = intercept2[:, 0, :, None]
 
-
         idx = np.s_[1:, 1:]
         self.slope[idx], self.intercept[idx], _ = self._compute_fit(mrq1, mrq2)
 
         # self.margin1_log_cumulants = self.margin1_slope * np.log2(np.e)
         # self.margin2_log_cumulants = self.margin2_slope * np.log2(np.e)
 
-
         self.log_cumulants = self.slope * np.log2(np.e)
 
-
         self._compute_rho()
 
     def _compute(self, mrq1, mrq2, idx_reject):
@@ -624,7 +620,7 @@ def plot(self, j1=None, j2=None, figsize=None, scaling_range=0,
         for m1 in range(self.n_cumul + 1):
             for m2 in range(self.n_cumul + 1):
 
-                if m1 != 0 and m2 !=0 and (m1, m2) not in self.m:
+                if m1 != 0 and m2 != 0 and (m1, m2) not in self.m:
                     fig.delaxes(axes[m1][m2])
                     continue
 
@@ -637,37 +633,51 @@ def plot(self, j1=None, j2=None, figsize=None, scaling_range=0,
         if filename is not None:
             plt.savefig(filename)
 
-    def _compute_legendre(self, h_support=(0, 1.5), resolution=100):
+    def _compute_legendre(self, h_support=(0, 1.5), resolution=100,
+                          signal_idx1=0, signal_idx2=1, idx_range=0):
         """
         Compute the bivariate Legendre spectrum.
         """
 
+        if signal_idx1 == signal_idx2:
+            raise ValueError('signal_idx1 should be different from signal_idx2')
+
         h_support = np.linspace(*h_support, resolution)
 
-        b = (self.c20 * self.c02) - (self.c11 ** 2)
+        sl_ = np.s_[:, signal_idx1, signal_idx2, idx_range]
+
+        c11 = self.c11[sl_]
+        c10 = self.c10[sl_]
+        c01 = self.c01[sl_]
+        c20 = self.c20[sl_]
+        c02 = self.c02[sl_]
+
+        b = (c20 * c02) - (c11 ** 2)
 
         L = np.ones((resolution, resolution))
 
         for i, h in enumerate(h_support):
-            L[i, :] += self.c02 * b / 2 * (((h - self.c10) / b) ** 2)
-            L[:, i] += self.c20 * b / 2 * (((h - self.c01) / b) ** 2)
+            L[i, :] += c02 * b / 2 * (((h - c10) / b) ** 2)
+            L[:, i] += c20 * b / 2 * (((h - c01) / b) ** 2)
 
         for i, h1 in enumerate(h_support):
             for j, h2 in enumerate(h_support):
-                L[i, j] -= (self.c11 * b
-                            * ((h1 - self.c10) / b)
-                            * ((h2 - self.c01) / b))
+                L[i, j] -= (c11 * b
+                            * ((h1 - c10) / b)
+                            * ((h2 - c01) / b))
 
         return h_support, L
 
-    def plot_legendre(self, h_support=(0, 1.5), resolution=30,
-                      figsize=(10, 10), cmap=None):
+    def plot_legendre(self, h_support=(0, 1.5), resolution=100,
+                      figsize=(10, 10), cmap=None, signal_idx1=0,
+                      signal_idx2=1, idx_range=0):
         """
         Plot the bivariate Legendre spectrum
         """
 
-        h, L = self._compute_legendre(h_support=h_support,
-                                      resolution=200)
+        h, L = self._compute_legendre(
+            h_support=h_support, resolution=200, signal_idx1=signal_idx1,
+            signal_idx2=signal_idx2, idx_range=idx_range)
 
         h_x = h[L.max(axis=0) >= 0]
         h_y = h[L.max(axis=1) >= 0]
@@ -677,34 +687,58 @@ def plot_legendre(self, h_support=(0, 1.5), resolution=30,
 
         h, L = self._compute_legendre((hmin, hmax), resolution)
 
-        cmap = cmap or plt.cm.coolwarm  # pylint: disable=no-member
+        cmap = cmap or plt.cm.viridis  # pylint: disable=no-member
+
+        fig, ax = plt.subplots(
+            figsize=figsize, subplot_kw={'projection': '3d'})
+        # ax = fig.add_subplot(1, 1, 1, projection='3d')
 
-        fig = plt.figure(figsize=figsize)
-        ax = fig.gca(projection='3d')
+        ax.set_proj_type('persp', focal_length=.15)
 
         h_x = h[L.max(axis=0) >= 0]
         h_y = h[L.max(axis=1) >= 0]
 
         L = L[:, L.max(axis=0) >= 0]
         L = L[L.max(axis=1) >= 0, :]
 
+        L[L < 0] = None
+
         colors = cmap(L)
         colors[L < 0] = 0
 
         X, Y = np.meshgrid(h_x, h_y)
 
+        light = mpl.colors.LightSource(azdeg=60, altdeg=60)
+
         # Plot the surface.
-        surf = ax.plot_surface(X, Y, L, cmap=cmap,
-                               linewidth=1, antialiased=False,
-                               vmin=0, vmax=1,
-                               rstride=1, cstride=1,
-                               facecolors=colors, shade=False, linestyle='-',
-                               edgecolors='black', zorder=1)
+        surf = ax.plot_surface(X, Y, L, alpha=.95, cmap=cmap,
+                            #   facecolors=colors,
+                            # lightsource=light,
+                            #    linewidth=1, vmin=0, vmax=1,
+                            #    rstride=1, cstride=1,
+                            # linestyle='-',
+                            #    zorder=1)
+        )
+
+        # argmax = np.argmax(L, axis=-1)
+        # ax.contour(X[..., argmax], Y[..., argmax], L[..., argmax], zdir='x')
+
+        # ax.contour(X, Y, L > 1, zdir='x', offset=h_x.min(), levels=0)
+
+        ax.set(xlabel='$h_1$', ylabel='$h_2$', zlabel='$\mathcal{L}(h_1, h_2)$')
+
+        # ax.contour(X, Y, L, zdir='y', offset=h_y.max(), levels=[0])
+
+        # ax.contour(X, Y, L, zdir='x', levels=10)
+        # ax.contour(X, Y, L, zdir='y', levels=10)
+        # ax.contour(X, Y, L, zdir='z', levels=10)
+
+        # surf.set_edgecolors((0.1, 0.2, 0.5, 1))
         ax.set_zlim(0, 1)
         ax.view_init(elev=45)
 
         # TODO manage to plot the contours or switch to 3D plotting libs
-        fig.colorbar(surf, shrink=0.6, aspect=10)
+        # fig.colorbar(surf, shrink=0.6, aspect=10)
 
     def plot_legendre_pv(self, resolution=30, figsize=(10, 10), cmap=None,
                          use_ipyvtk=False):

setup.py

@@ -20,7 +20,8 @@
     'test': [
         'pytest', 'pytest-xdist', 'pytest-cov', 'statsmodels', 'recombinator',
         'joblib'],
-    'doc': ['sphinx', 'numpydoc', 'pydata-sphinx-theme', 'sphinx-gallery', 'joblib'],
+    'doc': ['sphinx', 'numpydoc', 'pydata-sphinx-theme', 'sphinx-gallery', 'joblib',
+            'pooch'],
 }
 
 setup(name='pymultifracs',

tests/bivar_test.py

@@ -48,11 +48,11 @@ def test_key(key):
         lwt = bimfa(
             WTpL, WTpL, scaling_ranges, weighted=None, n_cumul=2,
             q1=np.array([0, 1, 2]), q2=np.array([0, 1, 2]), R=1)
-        
+
         lwt.structure.get_jrange()
         lwt.structure.plot()
         lwt.cumulants.plot()
-        # lwt.cumulants.plot_legendre()
+        lwt.cumulants.plot_legendre()
 
         p = param[key]