Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Outlier detection and robust estimation #22

Open
wants to merge 102 commits into
base: master
Choose a base branch
from
Open
Changes from 1 commit
Commits
Show all changes
102 commits
Select commit Hold shift + click to select a range
e8f1628
fix edge case
MerlinDumeur Jun 27, 2023
b2d1182
Merge branch 'robust' of github.com:MerlinDumeur/pymultifracs into ro…
MerlinDumeur Jun 27, 2023
b6ab233
Prepared figures
MerlinDumeur Jun 27, 2023
f0cf3dc
Improved robust algorithm
MerlinDumeur Aug 8, 2023
839514b
Merge branch 'robust' of github.com:MerlinDumeur/pymultifracs into ro…
MerlinDumeur Sep 1, 2023
248b35f
fix edge case
MerlinDumeur Jun 27, 2023
9a23e56
Prepared figures
MerlinDumeur Jun 27, 2023
4e9a8a6
Improved robust algorithm
MerlinDumeur Aug 8, 2023
d3cf506
Merge branch 'robust' of github.com:MerlinDumeur/pymultifracs into ro…
MerlinDumeur Sep 1, 2023
9fe0cca
EURASIP modifs
MerlinDumeur Sep 9, 2023
b3aae5f
Tiny fix
MerlinDumeur Sep 11, 2023
29330cc
Added custom width of leader blocks + on the fly rejection of leaders
MerlinDumeur Oct 4, 2023
0ff3844
Fixed beta estimation
MerlinDumeur Oct 5, 2023
f707405
Fixing functions
MerlinDumeur Oct 17, 2023
48c30bf
fix edge case
MerlinDumeur Jun 27, 2023
f4007bd
Prepared figures
MerlinDumeur Jun 27, 2023
805b739
Improved robust algorithm
MerlinDumeur Aug 8, 2023
48723c9
fix edge case
MerlinDumeur Jun 27, 2023
7a3399c
Prepared figures
MerlinDumeur Jun 27, 2023
1c4b122
Improved robust algorithm
MerlinDumeur Aug 8, 2023
2eba38a
EURASIP modifs
MerlinDumeur Sep 9, 2023
f41cdb6
Tiny fix
MerlinDumeur Sep 11, 2023
da8c1af
Added custom width of leader blocks + on the fly rejection of leaders
MerlinDumeur Oct 4, 2023
de84045
Fixed beta estimation
MerlinDumeur Oct 5, 2023
8200a03
Fixing functions
MerlinDumeur Oct 17, 2023
b7c3444
Fixing tests
MerlinDumeur Oct 26, 2023
a11ab88
Merge
MerlinDumeur Oct 26, 2023
a785166
Prepared last method
MerlinDumeur Dec 4, 2023
b1d3f33
Fix cumulant plotting bug
MerlinDumeur Dec 19, 2023
bf2e05c
Fix PSD plotting function
MerlinDumeur Dec 21, 2023
d6f63c9
Fixed robust cumulant estimation
MerlinDumeur Jan 25, 2024
95b244e
Added multifractal spectra with rejected coefficients; individual gam…
MerlinDumeur Feb 23, 2024
6b39772
Plotting function allow shift in gamint
MerlinDumeur Feb 23, 2024
999fd63
plot improvement
MerlinDumeur Feb 29, 2024
51d094a
Pep8 improvement, doc, prepare robust
MerlinDumeur Mar 1, 2024
2ddc566
Fixed debugging
MerlinDumeur Mar 1, 2024
cf5a209
Fixing wavelet function
MerlinDumeur Mar 4, 2024
2b80d83
Faster cumulants
MerlinDumeur Mar 5, 2024
a74f943
Update doc
MerlinDumeur Mar 5, 2024
f0d6524
Fixed type linking
MerlinDumeur Mar 5, 2024
3430416
Modifying API
MerlinDumeur Mar 7, 2024
038475d
Progress on new API
MerlinDumeur Mar 18, 2024
2a5b168
API rework
MerlinDumeur Mar 21, 2024
1225ec1
API rework 2 -- not tested
MerlinDumeur Mar 22, 2024
e1efddb
Fix imports
MerlinDumeur Mar 22, 2024
3f58573
Fixing WT transform and plotting
MerlinDumeur Mar 25, 2024
8bd85d8
Fixed plotting and bootstrapping
MerlinDumeur Mar 26, 2024
8e869f0
Fixing tests
MerlinDumeur Mar 27, 2024
f2cea67
Fixed regression
MerlinDumeur Mar 28, 2024
6a04f8b
Started bivariate analysis
MerlinDumeur Apr 2, 2024
97e9f15
Bivariate sf WIP
MerlinDumeur Apr 2, 2024
df3bf43
Bivariate structure function done
MerlinDumeur Apr 3, 2024
5b73fe5
Restructure dataclasses, add bivariate structure plotting
MerlinDumeur Apr 3, 2024
24dd51a
Added bivariate cumulants
MerlinDumeur Apr 4, 2024
87bc901
Fixing cumulants
MerlinDumeur Apr 5, 2024
517ad58
Fixing cumulants
MerlinDumeur Apr 5, 2024
59d0227
Fixed bivariate -- passing tests
MerlinDumeur Apr 8, 2024
b866240
Bootstrap WIP
MerlinDumeur Apr 8, 2024
f41cdfb
Bootstrap WIP
MerlinDumeur Apr 9, 2024
db5b47e
Bootstrap WIP
MerlinDumeur Apr 9, 2024
4c9a422
Bootstrap test successful
MerlinDumeur Apr 9, 2024
91db93f
Last test fix
MerlinDumeur Apr 10, 2024
0fca2bb
Refresh doc
MerlinDumeur Apr 10, 2024
b54e531
Integrated robust estimation
MerlinDumeur Apr 11, 2024
2d142c4
Doc update, WT analysis fixed, added example notebook
MerlinDumeur Apr 12, 2024
6be815e
Trying to fix CI
MerlinDumeur Apr 12, 2024
9e96f8f
Fixed 3.10 compatibility
MerlinDumeur Apr 13, 2024
72d2319
Fixed CI error, conda -> mamba
MerlinDumeur Apr 14, 2024
6ed5150
Micromamba CI fix
MerlinDumeur Apr 14, 2024
ef96994
CI fix 2
MerlinDumeur Apr 14, 2024
761fddb
Fixing bootstrapping bug
MerlinDumeur Apr 19, 2024
1198bf2
Bootstrap demo notebook
MerlinDumeur May 13, 2024
ffb6513
Benchmarking code improvement
MerlinDumeur May 25, 2024
ecd4f5c
Fixing WT transform for outlier detection
MerlinDumeur May 27, 2024
5a6f9aa
Fixing plot function
MerlinDumeur May 28, 2024
8b768dd
Fix WSE
MerlinDumeur Nov 27, 2024
501b3c8
Update installation process
MerlinDumeur Jan 21, 2025
7bb3aa4
update meta.yml
MerlinDumeur Jan 21, 2025
612acf2
update ci,meta.yml
MerlinDumeur Jan 21, 2025
499d7be
Update CI
MerlinDumeur Jan 21, 2025
c73dcef
doc update
MerlinDumeur Jan 22, 2025
851cad5
updated codecov and linting actions
MerlinDumeur Jan 22, 2025
0463e8b
CI fix
MerlinDumeur Jan 22, 2025
18a36ff
Mamba CI fix
MerlinDumeur Jan 22, 2025
249d72a
CI fix
MerlinDumeur Jan 22, 2025
414e102
CI fix
MerlinDumeur Jan 22, 2025
310f207
Hopefully fixes CI
MerlinDumeur Jan 22, 2025
421e959
Hopefully fixes CI-2
MerlinDumeur Jan 22, 2025
774bb7a
Hopefully fixes CI-3
MerlinDumeur Jan 22, 2025
38955a9
Hopefully fixes CI-4
MerlinDumeur Jan 22, 2025
304616b
Add tests
MerlinDumeur Jan 22, 2025
6cc81bd
Doc progress
MerlinDumeur Jan 23, 2025
6e6e110
Test fix and doc fix
MerlinDumeur Jan 23, 2025
3536979
Bivar test and benchmark coverage
MerlinDumeur Jan 23, 2025
7acbf14
Fixed missing tqdm dependency
MerlinDumeur Jan 23, 2025
70e275f
Doc update
MerlinDumeur Jan 27, 2025
b0a4e3a
Final doc changes
MerlinDumeur Jan 27, 2025
f7656be
Formatting fix
MerlinDumeur Jan 27, 2025
69f5581
Merge triton local
MerlinDumeur Jan 29, 2025
444663e
Added python version and doc badges
MerlinDumeur Jan 29, 2025
3b89efe
Fixed merge
MerlinDumeur Jan 29, 2025
8a3e38d
Fix badge
MerlinDumeur Jan 29, 2025
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
Final doc changes
MerlinDumeur committed Jan 27, 2025
commit b0a4e3aafb16d533ace44fa0163ab1826aa65a33
26 changes: 13 additions & 13 deletions examples/01-simul.py
Original file line number Diff line number Diff line change
@@ -61,7 +61,7 @@
WT.plot(j1=6, j2=11)

# %% [markdown]
# Multi-resolution quantities derived from the wavelet transform can be obtained using the associated methods ``get_leaders`` and ``get_wse``:
# Multi-resolution quantities derived from the wavelet transform can be obtained using the associated methods :meth:`.WaveletDec.get_leaders` and :meth:`.WaveletDec.get_wse`:

# %%
WTL = WT.get_leaders(p_exp=np.inf)
@@ -85,7 +85,7 @@

# %% [markdown]
# In order for the analysis to be meaningful under the chosen multifractal formalism (wavelet coefficient, wavelet (p-)leader, etc.) it may be necessary to verify a minimum regularity condition.
# The method :func:`.check_regularity` is available with all multi-resolution quantities, and takes ``scaling_ranges`` as an argument:
# The method :meth:`.WaveletDec.check_regularity` is available with all multi-resolution quantities, and takes ``scaling_ranges`` as an argument:

# %%
WT.check_regularity(scaling_ranges)
@@ -96,13 +96,13 @@
# %% [markdown]
# In case the minimal regularity is too low, it may be necessary to fractionally integrate the time series.
#
# A simple approach is provided in the `.auto_integrate` method, which will try to find a fractional integration coefficient large enough that all signals may be analyzed, and return the properly integrated multi-resolution quantity.
# A simple approach is provided in the :meth:`.WaveletDec.auto_integrate` method, which will try to find a fractional integration coefficient large enough that all signals may be analyzed, and return the properly integrated multi-resolution quantity.

# %%
WTpL = WTpL.auto_integrate(scaling_ranges)

# %% [markdown]
# Otherwise, and for instance in the case where multiple sets of data need to be compared using the same integration coefficient, the fractional integration can be set using the :func:`.integrate` method on a MRQ object by passing the fractional integration coefficient :math:`\gamma`:
# Otherwise, and for instance in the case where multiple sets of data need to be compared using the same integration coefficient, the fractional integration can be set using the :meth:`.WaveletDec.integrate` method on a MRQ object by passing the fractional integration coefficient :math:`\gamma`:

# %%
WT_int = WT.integrate(.5)
@@ -114,15 +114,15 @@
# %% [markdown]
# Multifractal analysis is carried out using the :func:`mfa` function.
#
# Parameters:
# Basic parameters:
#
# - ``mrq``: Multi-resolution quantity (:class:`.WaveletDec`, :class:`.WaveletLeader`, :class:`WaveletWSE`) on which to perform the analysis.
# - ``mrq``: Multi-resolution quantity (:class:`.WaveletDec`, :class:`.WaveletLeader`, :class:`.WaveletWSE`) on which to perform the analysis.
#
# - ``weighted``: whether the linear regressions should be weighted. Defaults to None, which means no weighting is performed. ``"Nj"`` indicates that the weights are determined from the number of coefficients at each scale.
#
# - ``q``: list of moments.
#
# .. note:: by default, `mfa` checks the regularity of the time series. It is possible to disable this by using the `check_regularity` argument.
# .. note:: by default, :func:`mfa` checks the regularity of the time series. It is possible to disable this by passing ``check_regularity=False``.

# %%
from pymultifracs import mfa
@@ -145,30 +145,30 @@
# **Structure functions**

# %% [markdown]
# The structure functions :math:`S_q(j)` and their associated exponents may be visualized using the `.plot()` method
# The structure functions :math:`S_q(j)` and their associated exponents may be visualized using the :meth:`.StructureFunction.plot()` method

# %%
pwt.structure.plot(figsize=(10, 4), nrow=2)

# %% [markdown]
# We can plot :math:`\zeta(q)` using the :func:`.plot_scaling` method
# We can plot :math:`\zeta(q)` using the :meth:`StructureFunction.plot_scaling` method

# %%
pwt.structure.plot_scaling()

# %% [markdown]
# **Cumulants**
# **Cumulants**

# %% [markdown]
# The cumulant scaling functions may be visualized using
# The cumulant scaling functions may be visualized using :meth:`.Cumulants.plot`

# %%
pwt.cumulants.plot()

# %% [markdown]
# **Multifractal spectrum**
# Visualizing the multifractal spectrum requires more densely sampled values of $q$:
# Visualizing the multifractal spectrum requires more densely sampled values of :math:`q`:

# %%
pwt = mfa(WTpL, scaling_ranges, weighted='Nj', q=build_q_log(.1, 5, 20))
pwt.spectrum.plot()
pwt.spectrum.plot()
42 changes: 42 additions & 0 deletions pymultifracs/multiresquantity.py
Original file line number Diff line number Diff line change
@@ -325,6 +325,12 @@ def plot(self, j1, j2, ax=None, vmin=None, vmax=None, cbar=True,
figsize=figsize, gamma=gamma, nan_idx=nan_idx,
signal_idx=signal_idx, cbar_kw=cbar_kw, cmap=cmap)

def get_variable_name(self):
return "d"

def get_suffix(self):
return "", ""

def get_formalism(self):
"""
Obtains the fomalism of the multi-resolution quantity
@@ -423,6 +429,25 @@ def get_wse(self, theta=.5, gamint=None):
return integ.get_wse(theta, gamint)

def auto_integrate(self, scaling_ranges, weighted=None, idx_reject=None):
"""
Automatically integrates the signal to match the requirements of
the multifractal formalism associated to the multi-resolution
quantity.

Parameters
----------

scaling_ranges : list[tuple[int, int]]
List of pairs of :math:`(j_1, j_2)` ranges of scales for the analysis.
weighted : str | None
Weighting mode for the linear regressions. Defaults to None, which is
no weighting. Possible values are 'Nj' which weighs by number of
coefficients, and 'bootstrap' which weights by bootstrap-derived
estimates of variance.
idx_reject : dict[int, ndarray of bool]
Dictionary associating each scale to a boolean array indicating whether
certain coefficients should be removed.
"""

hmin, _ = estimation.estimate_hmin(
self, scaling_ranges, weighted, idx_reject)
@@ -574,6 +599,17 @@ def bootstrap(self, R, min_scale=1, idx_reject=None):

# return self.bootstrapped_obj

def get_var_name(self):

return r"\ell"

def get_suffix(self):

if self.p_exp == np.inf:
return "", ""

return f"^{{({self.p_exp})}}", f"^{{({self.p_exp})}}"

def get_formalism(self):
"""
Obtains the fomalism of the multi-resolution quantity
@@ -756,6 +792,12 @@ class WaveletWSE(WaveletDec):
"""
theta: float

def get_variable_name(self):
return r"\ell"

def get_suffix(self):
return rf"^{{({self.theta}, 1)}}", "^{ws}"

def get_formalism(self):
return 'weak scaling exponent'

21 changes: 13 additions & 8 deletions pymultifracs/scalingfunction.py
Original file line number Diff line number Diff line change
@@ -12,7 +12,7 @@
linear_regression, compute_R2
from .autorange import compute_Lambda, compute_R, find_max_lambda
from .utils import fast_power, mask_reject, isclose, fixednansum, \
AbstractDataclass
AbstractDataclass, Formalism
from . import multiresquantity, viz


@@ -23,7 +23,7 @@ class AbstractScalingFunction(AbstractDataclass):
weighted: str | None = None
n_sig: int = field(init=False)
j: np.ndarray = field(init=False)
formalism: str = field(init=False)
formalism: Formalism = field(init=False)
nj: dict[int, np.ndarray] = field(init=False, repr=False)
values: np.ndarray = field(init=False, repr=False)
slope: np.ndarray = field(init=False, repr=False)
@@ -83,13 +83,16 @@ def __getattr__(self, name):
@dataclass(kw_only=True)
class ScalingFunction(AbstractScalingFunction):
mrq: InitVar[multiresquantity.WaveletDec]
variable_suffix: str = field(init=False)
regularity_suffix: str = field(init=False)
gamint: float = field(init=False)

def __post_init__(self, idx_reject, mrq):

self.gamint = mrq.gamint
self.n_sig = mrq.n_sig
self.formalism = mrq.get_formalism()
self.variable_suffix, self.regularity_suffix = mrq.get_suffix()
self.j = np.array(list(mrq.values))

self.nj = mrq.get_nj_interv()
@@ -412,7 +415,7 @@ def plot(self, nrow=4, filename=None, ignore_q0=True, figsize=None,
ax = axes[counter % nrow][counter // nrow]
ax.errorbar(x, y, CI, fmt='r--.', zorder=4)
ax.set_xlabel('Temporal scale $j$')
ax.set_ylabel(f'$S_{{{q:.1g}}}(j)$')
ax.set_ylabel(f'$S_{{{q:.1g}}}{self.variable_suffix}(j)$')
ax.tick_params(bottom=False, top=False, which='minor')

counter += 1
@@ -433,7 +436,7 @@ def plot(self, nrow=4, filename=None, ignore_q0=True, figsize=None,
else:
CI_legend = ""

legend = rf'$s_{{{q:.1g}}}$ = {slope:.2f}' + CI_legend
legend = rf'$s_{{{q:.1g}}}{self.variable_suffix}$ = {slope:.2f}' + CI_legend

ax.plot([x0, x1], [y0, y1], color='k',
linestyle='-', linewidth=2, label=legend, zorder=5)
@@ -473,11 +476,12 @@ def plot_scaling(self, filename=None, ax=None, signal_idx=0, range_idx=0,
" value is used")

if ax is None:
_, ax = plt.subplots()
_, ax = plt.subplots(figsize=(4, 2.5), layout='tight')

ax.plot(self.q, self.slope[:, range_idx, signal_idx], **plot_kw)

ax.set(
xlabel = '$q$', ylabel=r'$\zeta(q)$',
xlabel = 'Moment $q$', ylabel=rf'Scaling function $\zeta{self.variable_suffix}(q)$',
# title=self.formalism + ' - scaling function'
)

@@ -588,7 +592,7 @@ def _compute_robust(self, mrq, idx_reject, **robust_kwargs):
T_X_j = np.abs(mrq.values[j])
T_X_j = T_X_j[:, None, :]

if self.formalism == 'wavelet p-leader':
if self.formalism == Formalism.wavelet_pleader:
T_X_j = T_X_j * mrq.ZPJCorr[None, :, :, ind_j]

log_T_X_j = np.log(T_X_j)
@@ -929,7 +933,8 @@ def plot(self, filename=None, ax=None, fmt='ko-', range_idx=0,
CI_Dq, CI_hq, fmt,
**plot_kwargs)

ax.set(xlabel='Regularity $h$', ylabel='Fractal dimension $D(h)$',
ax.set(xlabel=f'Regularity $h{self.regularity_suffix}$',
ylabel=rf'Fractal dimension $\mathcal{{L}}{self.variable_suffix}(h)$',
ylim=(0, 1.1), xlim=(0, 1.5),
title=self.formalism + ' - multifractal spectrum')

8 changes: 8 additions & 0 deletions pymultifracs/utils.py
Original file line number Diff line number Diff line change
@@ -3,6 +3,7 @@
Merlin Dumeur <merlin@dumeur.net>
"""

from enum import Enum
from dataclasses import dataclass
from typing import Any
import inspect
@@ -33,6 +34,13 @@
"""


class Formalism(Enum):
wavelet_coef = 1
wavelet_leader = 2
wavelet_pleader = 3
weak_scaling_exponent = 4


@dataclass
class AbstractDataclass:
bootstrapped_obj: Any | None = None
9 changes: 6 additions & 3 deletions pymultifracs/viz.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
import time
import os
from collections import namedtuple
from dataclasses import dataclass

import matplotlib as mpl
import matplotlib.pyplot as plt
@@ -303,7 +304,7 @@ def plot_cm(cm, ind_m, j1, j2, range_idx, ax, C_color='grey',

ax.errorbar(x, y, CI, fmt=C_fmt, color=C_color, lw=lw_C, **errobar_params)

ax.set(xlabel='Temporal scale $j$', ylabel=f'$C_{m}(j)$')
ax.set(xlabel='Temporal scale $j$', ylabel=f'$C_{m}{cm.variable_suffix}(j)$')

if len(cm.log_cumulants) > 0 and plot_fit:

@@ -330,7 +331,7 @@ def plot_cm(cm, ind_m, j1, j2, range_idx, ax, C_color='grey',
else:
CI_legend = ""

legend = (rf'$c_{m}$ = {cp_string_format(slope_log2_e)}'
legend = (rf'$c_{m}{cm.variable_suffix}$ = {cp_string_format(slope_log2_e)}'
+ CI_legend)

ax.plot([x0, x1], [y0, y1], color=fit_color,
@@ -505,13 +506,15 @@ def plot_coef(mrq, j1, j2, ax=None, vmin=None, vmax=None, cbar=True,
if cbar:
# formatter = mpl.ticker.LogFormatterSciNotation(labelOnlyBase=False, minor_thresholds=(np.inf, np.inf))
# formatter = mpl.ticker.LogFormatterSciNotation(labelOnlyBase=False, minor_thresholds=(np.inf, np.inf))
#

if cbar_kw is None:
cbar_kw = dict(
fraction=.1, aspect=8,
ticks=mpl.ticker.MaxNLocator(4, symmetric=False))

if 'label' not in cbar_kw:
cbar_kw['label'] = f"${mrq.get_variable_name()}{mrq.get_suffix()[0]}(j, k)$"

cb = plt.colorbar(qm, ax=ax, **cbar_kw)
# plt.colorbar(qm, ax=ax es[0], ticks=locator, aspect=1)
cb.ax.tick_params(which='major', size=3)