API Reference

Notes

FIF_1D and FIF_ND no longer accept the old combined kernel_construction_method argument. Use kernel_construction_method_flux and kernel_construction_method_observable instead.

Naive FIF kernels remain available for comparison only and emit a warning because their outputs are not remotely accurate.

Simulation

Simulation module for multifractal fields

Contains methods for generating synthetic multifractal time series and fields.

scaleinvariance.fBm_1D_circulant(size, H, periodic=True)[source]

Generate fractional Brownian motion by spectral synthesis.

Creates fBm with power spectrum |k|^(-β) where β = 2H + 1. Uses random phases and inverse FFT.

Parameters:

  • size (int) – Number of points (must be power of 2)

  • H (float) – Hurst parameter. Typical range is (0, 1) for fBm, but negative values are also supported (H < 0 gives blue noise with β < 1).

  • periodic (bool, optional) – If True (default), returns full periodic simulation suitable for periodic boundary conditions. If False, doubles simulation size internally then returns only the first half to eliminate periodicity artifacts.

Returns:

Generated 1D fBm, zero mean, unit std

Return type:

numpy.ndarray

scaleinvariance.fBm_ND_circulant(size, H, periodic=True)[source]

Generate N-dimensional fractional Brownian motion by spectral synthesis.

Creates N-D fBm with rotationally symmetric power spectrum |k|^(-β) where β = 2H + D. Uses random phases and inverse FFT to generate N-D random field.

Parameters:

  • size (tuple of ints) – Size of N-D array. Each dimension must be power of 2.

  • H (float) – Hurst parameter. Typical range is (0, 1) for fBm, but negative values are also supported (H < 0 gives blue noise with β < D).

  • periodic (bool or tuple of bool, optional) – If True (default), returns full periodic simulation suitable for periodic boundary conditions. If False, doubles simulation size internally along each dimension then returns only the first octant to eliminate periodicity artifacts. Can also be a tuple of bool for per-axis control.

Returns:

Generated N-D fBm field, zero mean, unit std

Return type:

numpy.ndarray

scaleinvariance.fBm_1D(size, H, causal=False, outer_scale=None, outer_scale_width_factor=2.0, periodic=True, gaussian_noise=None, kernel_construction_method='LS2010')[source]

Generate 1D fractional Brownian motion by fractional integration of Gaussian noise.

Supports extended Hurst range H in (-0.5, 1.5).

  • causal=False (default): spectral fractional integration of order H + 1/2. White noise has flat PSD, so multiplying its spectrum by |f|^(-(H+1/2)) gives fBm with target power spectrum |f|^(-(2H+1)). Handles every H in one pass — no branching on H < 0.5 / H == 0.5 / H > 0.5.

  • causal=True: real-space convolution with a corrected causal power-law kernel (LS2010 or naive). Branches on H: * H in (-0.5, 0.5): convolve noise with |x|^(H - 1/2) * H == 0.5: cumsum of noise (standard Brownian motion) * H in (0.5, 1.5): cumsum then convolve with |x|^(H - 3/2)

Parameters:

  • size (int) – Length of simulation (must be even, power of 2 recommended)

  • H (float) – Hurst exponent in (-0.5, 1.5). Controls correlation structure.

  • causal (bool, optional) – Use causal kernels (future doesn’t affect past). Default False. When False, the fast spectral fractional-integral path is used and kernel_construction_method is ignored.

  • outer_scale (int or None, optional) – Large-scale cutoff. If None (default), no outer-scale cutoff is applied: the causal real-space kernels skip the Hanning window and the non-causal spectral path imposes no low-frequency cutoff beyond unavoidable DC handling. If set, a Hanning-window cutoff centered at this scale is applied to the causal kernels, and it sets the low-frequency regularization (min_freq = 1/outer_scale) on the spectral path.

  • outer_scale_width_factor (float, optional) – Controls transition width for outer scale (causal path only). Default is 2.0. Has no effect when outer_scale is None.

  • periodic (bool, optional) – If True (default), returns full periodic simulation. If False, doubles the simulation size internally then returns the first half to eliminate periodicity artifacts.

  • gaussian_noise (ndarray, optional) – Pre-generated Gaussian noise for reproducibility. Must have same size as the requested output.

  • kernel_construction_method ({'LS2010', 'naive'}, optional) – Kernel construction method for the causal path. Ignored when causal=False.

Returns:

1D array of simulated fBm values (zero mean, unit std).

Return type:

numpy.ndarray

Raises:

ValueError – If size is not even, H is outside (-0.5, 1.5), or an unknown causal kernel_construction_method is given.

Examples

>>> # Standard acausal fBm with H=0.7
>>> fbm = fBm_1D(2048, H=0.7)
>>>
>>> # Anti-correlated process with H=-0.3
>>> fbm_blue = fBm_1D(2048, H=-0.3)
>>>
>>> # Causal LS2010 path
>>> fbm_causal = fBm_1D(2048, H=0.7, causal=True)

Notes

  • Computational complexity is O(N log N).

  • For the acausal path, the fractional-integration kernel is |f|^(-(H+1/2)); for causal it is the real-space Riesz kernel.

  • fBm_1D_circulant is an alternative (phases + power spectrum) that is only slightly cheaper but produces statistically equivalent output.

scaleinvariance.FIF_1D(size, alpha, C1, H, levy_noise=None, causal=False, outer_scale=None, outer_scale_width_factor=2.0, kernel_construction_method_flux='LS2010', kernel_construction_method_observable='spectral', periodic=True, observable_kernel_odd_axes=None)[source]

Generate a 1D Fractionally Integrated Flux (FIF) multifractal simulation.

FIF is a multiplicative cascade model that generates multifractal fields with specified Hurst exponent H, intermittency parameter C1, and Levy index alpha. The method follows Lovejoy & Schertzer 2010, including finite-size corrections.

Returns field normalized by mean (or zero mean, unit std for odd-axis output).

Algorithm:

  1. Generate extremal Lévy noise with stability parameter alpha

  2. Convolve with corrected kernel (finite-size corrected |k|^(-1/alpha)) to create log-flux

  3. Scale by (C1)^(1/alpha) and exponentiate to get conserved flux

  4. For H ≠ 0: Convolve flux with kernel |k|^(-1+H) to get observable

  5. For H < 0: Apply differencing to handle negative Hurst exponents

Parameters:

  • size (int) – Length of simulation (must be even, power of 2 recommended)

  • alpha (float) – Lévy stability parameter in [0.5, 2] and != 1. Controls noise distribution. Values < 0.5 are currently not implemented for FIF simulation.

  • C1 (float) – Codimension of the mean, controls intermittency strength. Must be > 0 for multifractal behavior. Special case: C1 = 0 routes to fBm (no intermittency), but requires causal=False since fBm cannot be causal.

  • H (float) – Hurst exponent. For kernel_construction_method_observable='spectral' (default), any finite real is accepted (modulo an overflow guard for very large positive H at large sizes; see fractional_integral_spectral()). For the 'LS2010' and 'naive' paths, H must lie in [-1, 1]; values in [-1, 0) are handled by integrating with H + 1 and then taking a first difference.

  • levy_noise (torch.Tensor, optional) – Pre-generated Lévy noise for reproducibility. Must have same size as simulation.

  • causal (bool, optional) – Use causal kernels (future doesn’t affect past). Default False. False gives symmetric (non-causal) kernels.

  • outer_scale (int or None, optional) – Large-scale cutoff. If None (default), no outer-scale cutoff is applied: the real-space kernels skip the Hanning window and the spectral path imposes no low-frequency cutoff beyond unavoidable DC handling. The LS2010 kernel’s intrinsic finite-size correction still acts. If set, a Hanning-window cutoff centered at this scale is applied to the real-space kernels, and it sets the low-frequency regularization (min_freq = 1/outer_scale) on the spectral path.

  • outer_scale_width_factor (float, optional) – Controls transition width for outer scale. Transition width = outer_scale * width_factor. Default is 2.0. Has no effect when outer_scale is None.

  • kernel_construction_method_flux (str, optional) – Method for constructing flux (cascade) kernel. Options: - ‘LS2010’: Lovejoy & Schertzer 2010 finite-size corrections (default) - ‘naive’: Simple power-law kernels without corrections

  • kernel_construction_method_observable (str, optional) –

    Method for constructing observable (H) kernel. Options: - ‘spectral’: Fourier-space power-law transfer function (default, non-causal only) - ‘LS2010’: Lovejoy & Schertzer 2010 finite-size corrections - ‘naive’: Simple power-law kernels without corrections - ‘spectral_odd’: Deprecated — equivalent to 'spectral' with

    observable_kernel_odd_axes=True. Will be removed in a future release; use the explicit form.

  • periodic (bool, optional) – If True (default), returns full periodic simulation suitable for periodic boundary conditions. If False, doubles simulation size internally then returns only the first half to eliminate periodicity artifacts.

  • observable_kernel_odd_axes (bool or tuple of bool, optional) – Make the observable kernel antisymmetric along the specified axes. Implemented as a Fourier-space phase factor i^k · Π sign(f_j), which leaves the magnitude spectrum unchanged (so Hurst recovery works identically to the even kernel) and yields a zero-mean output by construction. The output is normalized to zero mean and unit standard deviation rather than unit mean. Requires ``kernel_construction_method_observable=’spectral’`` — real-space sign-flipping of the LS2010 or naive kernel does not preserve the magnitude spectrum because the smooth correction term contributes a non-trivial low-frequency FT. The same axis cannot be both causal and odd; flux kernels cannot be made odd.

Returns:

1D array of simulated multifractal field values

Return type:

numpy.ndarray

Raises:

  • ValueError – If size is not even or H is outside valid range (-1, 1)

  • ValueError – If C1 < 0 or if C1 = 0 with causal=True (fBm cannot be causal)

  • ValueError – If provided levy_noise doesn’t match specified size

  • ValueError – If alpha < 0.5 or alpha == 1 (not implemented)

  • ValueError – If causal and observable_kernel_odd_axes are both set on the same axis.

Examples

>>> # Basic multifractal with strong intermittency
>>> fif = FIF_1D(1024, alpha=1.8, C1=0.1, H=0.3)
>>>
>>> # Monofractal case (C1=0 routes to fBm, requires causal=False)
>>> fbm = FIF_1D(1024, alpha=1.8, C1=0, H=0.3, causal=False)
>>>
>>> # Odd (antisymmetric) observable kernel — zero-mean output
>>> fif_odd = FIF_1D(1024, alpha=1.8, C1=0.1, H=0.3,
...                  observable_kernel_odd_axes=True)

Notes

  • Computational complexity is O(N log N) due to FFT-based convolutions

  • Large C1 values (> 0.5) can produce extreme values requiring careful handling

  • alpha < 0.5 and alpha = 1 are not currently implemented

  • C1 = 0: Routes internally to fBm_1D_circulant() for monofractal case (requires causal=False since fBm cannot be causal)

scaleinvariance.FIF_ND(size, alpha, C1, H, levy_noise=None, outer_scale=None, outer_scale_width_factor=2.0, kernel_construction_method_flux='LS2010', kernel_construction_method_observable='spectral', periodic=False, scale_metric=None, elliptical_dim=None, causal=False, observable_kernel_odd_axes=None, scale_metric_dim=None)[source]

Generate an N-D Fractionally Integrated Flux (FIF) multifractal simulation.

FIF is a multiplicative cascade model that generates multifractal fields with specified Hurst exponent H, intermittency parameter C1, and Levy index alpha. The method follows Lovejoy & Schertzer 2010, including finite-size corrections. Returns field normalized by mean.

Algorithm:

  1. Generate N-D extremal Lévy noise with stability parameter alpha

  2. Convolve with corrected kernel (finite-size corrected |r|^(-d/alpha)) to create log-flux

  3. Scale by (C1)^(1/alpha) and exponentiate to get conserved flux

  4. For H ≠ 0: Convolve flux with kernel |r|^(-d+H) to get observable

Parameters:

  • size (tuple of ints) – Size of simulation as tuple specifying dimensions (e.g., (height, width) for 2D, (depth, height, width) for 3D). All dimensions must be even numbers.

  • alpha (float) – Lévy stability parameter in [0.5, 2] and != 1. Controls noise distribution. Values < 0.5 are currently not implemented for FIF simulation.

  • C1 (float) – Codimension of the mean, controls intermittency strength. Must be > 0 for multifractal behavior. Special case: C1 = 0 routes to fBm_ND_circulant() (no intermittency).

  • H (float) – Hurst exponent. For kernel_construction_method_observable='spectral' (default), any finite real is accepted (modulo an overflow guard for very large positive H at large sizes). For the 'LS2010' path, H must lie in [0, 1] — negative H is not supported on that path.

  • levy_noise (ndarray, optional) – Pre-generated N-D Lévy noise for reproducibility. Must have same shape as simulation.

  • outer_scale (int or None, optional) – Large-scale cutoff. If None (default), no outer-scale cutoff is applied: the real-space kernels skip the Hanning window and the spectral path imposes no low-frequency cutoff beyond unavoidable DC handling. The LS2010 kernel’s intrinsic finite-size correction still acts. If set, a Hanning-window cutoff centered at this scale is applied to the real-space kernels, and it sets the low-frequency regularization (min_freq = 1/outer_scale) on the spectral path.

  • outer_scale_width_factor (float, optional) – Controls transition width for outer scale. Transition width = outer_scale * width_factor. Default is 2.0. Has no effect when outer_scale is None.

  • kernel_construction_method_flux (str, optional) – Method for constructing flux (cascade) kernel. Options: - ‘LS2010’: Lovejoy & Schertzer 2010 finite-size corrections (default)

  • kernel_construction_method_observable (str, optional) – Method for constructing observable (H) kernel. Options: - ‘spectral’: Fourier-space power-law transfer function (default, requires no scale_metric) - ‘LS2010’: Lovejoy & Schertzer 2010 finite-size corrections

  • periodic (bool or tuple of bool, optional) – Controls periodicity behavior for each axis. - If bool: applies same periodicity to all axes (default False) - If tuple: specifies per-axis periodicity, must have length matching size tuple When an axis is non-periodic (False), the domain is doubled along that axis internally and one hyper-octant is returned to suppress periodic artifacts.

  • scale_metric (ndarray, optional) – Custom N-D array defining the scale metric (generalized distance) for GSI norms. If None (default), standard Euclidean distance is computed. Shape requirement: Must match simulation domain shape after accounting for periodicity (doubled along non-periodic axes). For example, if size=(256, 256) and periodic=False, scale_metric.shape must be (512, 512). Note: For LS2010 method, the metric should use dx=2 spacing to match the kernel construction grid.

  • elliptical_dim (float, optional) – Elliptical dimension D_el of the GSI scale metric — the effective scaling dimension that may differ from the embedding spatial dimension. For canonical anisotropy with stratification exponent Hz: D_el = 1 + Hz (2D), 2 + Hz (3D), d - 1 + Hz in d-D. If None (default), uses the spatial dimension (no anisotropy correction).

  • scale_metric_dim (float, optional) – Deprecated — renamed to elliptical_dim. If both are passed, a ValueError is raised. scale_metric_dim will be removed in a future release.

  • causal (bool or tuple of bool, optional) – Make the FIF kernels causal along the specified axes. False (default) gives the fully acausal field. True makes every axis causal; a tuple of bool specifies per-axis. Causality zeros the half of the kernel where the axis coordinate is negative. Variance compensation scales as 2^k for k causal axes (the log-flux is multiplied by (2^k · C1)^(1/α)). Causality is applied to both flux and observable kernels for consistency. Not supported with kernel_construction_method_observable='spectral'.

  • observable_kernel_odd_axes (bool or tuple of bool, optional) – Make the observable kernel antisymmetric along the specified axes. Implemented as a Fourier-space phase factor i^k · Π sign(f_j), which preserves the magnitude spectrum and yields zero-mean output by construction. Output is normalized to zero mean and unit standard deviation rather than unit mean. Requires ``kernel_construction_method_observable=’spectral’`` — the real-space LS2010 kernel mixes a power law with a smooth correction, and sign-flipping the corrected kernel does not preserve the magnitude spectrum. As a consequence, odd_axes is also not available alongside a custom scale_metric (which itself requires the LS2010 path). The same axis cannot be both causal and odd; flux kernels cannot be made odd.

Returns:

N-D array of simulated multifractal field values

Return type:

numpy.ndarray

Raises:

  • ValueError – If size dimensions are not even or H is outside valid range (0, 1)

  • ValueError – If C1 < 0 (must be non-negative; C1 = 0 routes to fBm)

  • ValueError – If provided levy_noise doesn’t match specified size

  • ValueError – If alpha < 0.5 or alpha == 1 (not implemented)

  • ValueError – If causal and observable_kernel_odd_axes are both set on the same axis, or if either is given a tuple whose length does not match size.

Examples

>>> # 2D multifractal with strong intermittency
>>> fif = FIF_ND((512, 512), alpha=1.8, C1=0.1, H=0.3)
>>>
>>> # 3D multifractal
>>> fif = FIF_ND((128, 128, 128), alpha=1.5, C1=0.05, H=0.2)
>>>
>>> # 2D with per-axis periodicity (periodic in x, non-periodic in y)
>>> fif = FIF_ND((256, 256), alpha=1.6, C1=0.08, H=0.4, periodic=(True, False))
>>>
>>> # Monofractal case (C1=0 routes to fBm_ND_circulant)
>>> fbm = FIF_ND((512, 512), alpha=1.8, C1=0, H=0.3)
>>>
>>> # Using custom GSI norm (example with anisotropic scaling)
>>> import numpy as np
>>> size = (256, 256)
>>> sim_size = (512, 512)  # doubled for non-periodic
>>> # Create coordinate arrays with dx=2 spacing (to match LS2010)
>>> x = np.arange(-(sim_size[0]-1), sim_size[0], 2)
>>> y = np.arange(-(sim_size[1]-1), sim_size[1], 2)
>>> X, Y = np.meshgrid(x, y, indexing='ij')
>>> # Define anisotropic GSI metric with different scaling in x and y
>>> scale_metric = (X**2 + (Y/2)**2)**0.5  # y-direction scales faster
>>> # Use non-integer dimension for GSI
>>> fif = FIF_ND(size, alpha=1.7, C1=0.1, H=0.3, periodic=False,
...              scale_metric=scale_metric, elliptical_dim=2.3,
...              kernel_construction_method_observable='LS2010')

Notes

  • Computational complexity is O(N log N) due to FFT-based convolutions

  • Large C1 values (> 0.5) can produce extreme values requiring careful handling

  • alpha < 0.5 and alpha = 1 are not currently implemented

  • Causality is opt-in (causal parameter); the default is fully acausal. Causality is applied to both flux and observable kernels.

  • Spectral N-D kernels currently assume the default isotropic Fourier metric and do not support custom scale_metric. They also do not support causal.

  • The 'LS2010' observable path does not support negative H values; for negative H, use kernel_construction_method_observable='spectral' (the default) which handles any real H.

  • C1 = 0: Routes internally to fBm_ND_circulant() for monofractal case (causal is not supported with C1=0 because fBm cannot be causal).

Generalized Scale Invariance (GSI)

FIF_ND supports anisotropic GSI scaling via a user-supplied scale_metric (an N-D array of generalized distances) together with the elliptical_dim parameter (the GSI elliptical dimension D_el). For the canonical self-affine anisotropy with stratification exponent Hz, D_el = d - 1 + Hz in d dimensions. The canonical_scale_metric helper builds the standard stratified metric (last axis anisotropic, all others isotropic) on the dx=2 grid expected by the LS2010 kernel construction.

GSI in FIF_ND currently requires kernel_construction_method_observable='LS2010' — the spectral observable path does not accept a custom scale_metric.

scaleinvariance.canonical_scale_metric(size, ls, Hz=0.5555555555555556)[source]

Compute the canonical anisotropic scale metric for GSI multifractal models.

This metric introduces anisotropy in the vertical (z) direction, allowing different scaling behavior along different axes. The metric is defined as:

For 3D:

r = ls * sqrt((dx/ls)^2 + (dy/ls)^2 + (dz/ls)^(2/Hz))

For 2D (x-z plane, where axis 1 is vertical):

r = ls * sqrt((dx/ls)^2 + (dz/ls)^(2/Hz))

where ls is the spheroscale (characteristic horizontal scale) and Hz controls the degree of anisotropy in the vertical direction.

IMPORTANT: Uses dx=2 grid spacing to match LS2010 kernel construction. Coordinate arrays are created as arange(-(size[i]-1), size[i], 2) for each dimension.

Parameters:

  • size (tuple of ints) – Size of the simulation domain (e.g., (512, 512) for 2D, (256, 256, 128) for 3D). Must be a 2D or 3D tuple.

  • ls (float) – Spheroscale - characteristic horizontal scale length. Must be > 0. This sets the fundamental length scale for the metric.

  • Hz (float, optional) – Anisotropy exponent for vertical (z) direction. Default is 5/9 ≈ 0.556 (empirical value for atmospheric stratification). - Hz = 1.0: Isotropic (standard Euclidean metric) - Hz < 1.0: Vertical direction scales slower (stratified/layered structures) - Hz > 1.0: Vertical direction scales faster Applied to the last axis in 2D and 3D cases.

Returns:

N-D array of scale metric values with shape matching size. Always returns a NumPy array regardless of backend.

Return type:

numpy.ndarray

Raises:

ValueError – If size is not a 2D or 3D tuple If ls <= 0 If Hz <= 0

Examples

>>> # 2D x-z plane with stratification (Hz < 1, vertical scales slower)
>>> metric_2d = canonical_scale_metric((512, 256), ls=100.0, Hz=0.556)
>>>
>>> # 2D isotropic (set Hz=1.0)
>>> metric_2d_iso = canonical_scale_metric((512, 512), ls=100.0, Hz=1.0)
>>>
>>> # 3D with stratification (Hz < 1, vertical scales slower)
>>> metric_3d = canonical_scale_metric((256, 256, 128), ls=50.0, Hz=0.5)
>>>
>>> # Use with FIF_ND for GSI multifractal simulation
>>> from scaleinvariance import FIF_ND
>>> size = (256, 256, 128)
>>> ls = 50.0
>>> Hz = 0.6
>>> metric = canonical_scale_metric(size, ls, Hz)
>>> # For non-periodic simulation, need to double the domain
>>> sim_size = tuple(s * 2 for s in size)
>>> metric_doubled = canonical_scale_metric(sim_size, ls, Hz)
>>> fif = FIF_ND(size, alpha=1.7, C1=0.1, H=0.3, periodic=False,
...              scale_metric=metric_doubled, elliptical_dim=2.5,
...              kernel_construction_method_observable='LS2010')

Notes

  • The metric reduces to standard Euclidean distance when Hz=1.0 and ls=1.0

  • For atmospheric applications, Hz < 1 is typical due to stratification

  • The ls parameter can be tuned to match physical scales in the system

  • Grid spacing is dx=2 to match LS2010 kernel construction in FIF_ND

References

Lovejoy, S., & Schertzer, D. (2007). Scaling and multifractal fields in the solid earth and topography. Nonlinear Processes in Geophysics, 14(4), 465-502.

Analysis

Analysis module for multifractal fields

Contains methods for estimating scaling exponents and analyzing multifractal properties.

Hurst Estimation

scaleinvariance.haar_fluctuation_hurst(data, min_sep=None, max_sep=None, axis=0, return_fit=False, periodic=False)[source]

Hurst exponent from Haar-fluctuation scaling.

Thin wrapper for wavelet_fluctuation_hurst(..., wavelet='haar').

scaleinvariance.wavelet_fluctuation_hurst(data, wavelet='haar', min_sep=None, max_sep=None, axis=0, return_fit=False, periodic=False)[source]

Estimate the Hurst exponent from first-order wavelet-fluctuation scaling.

For a self-similar process the order-1 fluctuation scales as F(r) ~ r^H; H is the slope of log F(r) vs log r. The accessible H range depends on the wavelet’s vanishing moments (Haar: roughly 0 < H < 1; the Mexican hat and Morlet reach higher H).

Parameters:

  • data (array-like) – N-dimensional scalar data array.

  • wavelet (str or Wavelet, optional) – Wavelet to use (default: ‘haar’).

  • min_sep (int, optional) – Lag window for the linear fit. Defaults match haar_fluctuation_hurst.

  • max_sep (int, optional) – Lag window for the linear fit. Defaults match haar_fluctuation_hurst.

  • axis (int, optional) – Axis along which to compute (default: 0).

  • return_fit (bool, optional) – If True, also return (lags, fluctuation values, fit_line).

  • periodic (bool, optional) – Periodic (toroidal) circular convolution along axis.

Return type:

H, uncertainty[, lags, values, fit_line]

scaleinvariance.structure_function_hurst(data, min_sep=None, max_sep=None, axis=0, return_fit=False, periodic=False)[source]

Calculate the Hurst exponent using structure functions.

For a self-similar process with Hurst exponent H, the first order structure function scales as: S_1(r) ~ r^H

This function estimates H by performing a linear regression on log(S_1(r)) vs log(r).

Parameters:

  • data (array-like) – N-dimensional scalar data array.

  • min_sep (int, optional) – Minimum separation/lag to include in the fit. If None, defaults to 1.

  • max_sep (int, optional) – Maximum separation/lag to include in the fit. If None, defaults to array_size/8 or array_size if array_size < 128 (with warning).

  • axis (int, optional) – Axis along which to compute the structure function (default: 0).

  • return_fit (bool, optional) – If True, return the lags, structure function values, and fitted line.

  • periodic (bool, optional) – If True, use periodic (toroidal) increments along axis (see structure_function). The default max_sep becomes array_size // 2.

Returns:

  • If return_fit is False –

    Hfloat

    Estimated Hurst exponent.

    uncertaintyfloat

    Standard error of the Hurst exponent estimate.

  • If return_fit is True –

    Hfloat

    Estimated Hurst exponent.

    uncertaintyfloat

    Standard error of the Hurst exponent estimate.

    lagsnp.ndarray

    The lags used in the calculation.

    sf_valuesnp.ndarray

    The structure function values.

    fit_linenp.ndarray

    The fitted line used to estimate H.

scaleinvariance.spectral_hurst(data, max_wavelength=None, min_wavelength=None, nbins=50, axis=0, return_fit=False)[source]

Estimate the Hurst exponent from the power spectral density along specified axis.

For fractional Brownian motion, PSD ~ f^(-β) with β = 2H + 1. This function estimates H by performing linear regression on log(PSD) vs log(frequency).

Note: this estimator is FFT-based and therefore inherently treats the domain as periodic, so there is no periodic flag (unlike structure_function_hurst / haar_fluctuation_hurst).

Parameters:

  • data (array-like) – N-dimensional input data array

  • max_wavelength (float, optional) – Maximum wavelength to include in analysis. If None, defaults to domain_size/8 or domain_size if domain_size < 512 (with warning).

  • min_wavelength (float, optional) – Minimum wavelength to include in analysis. If None, defaults to 4*Nyquist or Nyquist if domain_size < 512.

  • nbins (int, optional) – Number of logarithmic bins for PSD (default: 50)

  • axis (int, optional) – Axis along which to compute FFT (default: 0)

  • return_fit (bool, optional) – If True, return frequencies, PSD values, and fitted line

Returns:

  • If return_fit is False –

    Hfloat

    Estimated Hurst exponent

    uncertaintyfloat

    Standard error of the Hurst exponent estimate

  • If return_fit is True –

    Hfloat

    Estimated Hurst exponent

    uncertaintyfloat

    Standard error of the Hurst exponent estimate

    freqsnp.ndarray

    The frequencies used in the calculation

    psd_valuesnp.ndarray

    The power spectral density values

    fit_linenp.ndarray

    The fitted line used to estimate H

Structure Functions

scaleinvariance.structure_function(data, order=1, max_sep=None, axis=0, lags='powers of 1.2', periodic=False)[source]

Calculate the mean structure function for regularly spaced scalar data along a given axis.

The structure function S_n(r) of order n is defined as:

S_n(r) = <|u(x + r) - u(x)|^n>

where <> denotes the average over all valid pairs separated by lag r.

Parameters:

  • data (array-like) – N-dimensional scalar data array.

  • order (float or 1-D array-like of float, optional) – Order(s) of the structure function (default: 1). Each order must be positive. When an array is given, the absolute increments |u(x+r) - u(x)| are computed once per lag and reused across all orders (avoids redundant work).

  • max_sep (int, optional) – Maximum lag to compute. If None, defaults to (size along axis) - 1, or (size along axis) // 2 when periodic=True.

  • axis (int, optional) – Axis along which to compute the structure function (default: 0).

  • lags (str, list, or np.ndarray, optional) –

    Option for selecting lag values. Can be:

    • A list or array of lags.

    • ’all’: all integer lags from 1 to max_sep.

    • ’powers of X’: lags as integers of powers of X, where X is any number (e.g., ‘powers of 2’, ‘powers of 1.2’).

  • periodic (bool, optional) – If True, treat the domain as periodic (toroidal) along axis: the increment at lag r uses all L pairs u((x + r) mod L) - u(x) rather than dropping the L - r edge pairs. On a periodic domain S(r) = S(L - r) exactly, so lags above L // 2 are redundant; max_sep is therefore capped at L // 2 and an explicit max_sep > L // 2 raises ValueError.

Returns:

lags : 1-D array of lag values. sf : structure function values.

Shape is (n_lags,) if order is a scalar, or (n_orders, n_lags) if order is a 1-D array. The lag axis is always last.

Return type:

tuple (np.ndarray, np.ndarray)

scaleinvariance.costructure_function(data1, data2, order1=1, order2=1, max_sep=None, axis=0, lags='powers of 1.2', periodic=False)[source]

Calculate the mixed (two-field, two-exponent) structure function.

Defined as:

S_{p,q}(r) = <|u_1(x + r) - u_1(x)|^p * |u_2(x + r) - u_2(x)|^q>

where <> denotes the average over positions where neither increment is NaN. (Matches the NaN-only masking convention of structure_function and haar_fluctuation; inf values are not specially handled.)

The diagonal p = q case recovers a single-exponent co-structure function; the p -> 0 or q -> 0 slices recover the single-field structure functions of the other field (up to the caveat that 0 is not a permitted order here).

Parameters:

  • data1 (array-like) – N-dimensional scalar data arrays. Must have identical shape.

  • data2 (array-like) – N-dimensional scalar data arrays. Must have identical shape.

  • order1 (float or 1-D array-like of float, optional) – Exponent(s) applied to |delta f_1| and |delta f_2| respectively (default: 1). Each order must be positive. Arrays are allowed on either or both; the absolute increments of each field are computed once per lag and reused across all (order1, order2) combinations.

  • order2 (float or 1-D array-like of float, optional) – Exponent(s) applied to |delta f_1| and |delta f_2| respectively (default: 1). Each order must be positive. Arrays are allowed on either or both; the absolute increments of each field are computed once per lag and reused across all (order1, order2) combinations.

  • max_sep (int, optional) – Maximum lag to compute. If None, defaults to (size along axis) - 1.

  • axis (int, optional) – Axis along which to compute increments (default: 0).

  • lags (str, list, or np.ndarray, optional) – Lag selection (see structure_function for accepted forms).

  • periodic (bool, optional) – If True, treat the domain as periodic (toroidal) along axis: both increments wrap around the array end (all L pairs per lag). max_sep is capped at L // 2 and an explicit max_sep > L // 2 raises ValueError (see structure_function).

Returns:

lags : 1-D array of lag values. cosf : co-structure function values.

Shape depends on order1/order2 being scalar vs array: - (n_lags,) both scalar - (n_order1, n_lags) order1 array, order2 scalar - (n_order2, n_lags) order1 scalar, order2 array - (n_order1, n_order2, n_lags) both arrays The lag axis is always last.

Return type:

tuple (np.ndarray, np.ndarray)

Notes

The average is over positions where neither increment is NaN. This is asymmetric with structure_function, which averages over positions where the single field’s increment is non-NaN. If one field has NaNs and the other does not, the effective sample count here is the intersection of valid positions.

inf values are not specially handled and will propagate through to the result (consistent with the NaN-only masking of the single-field structure_function).

Analysis Functions

scaleinvariance.wavelet_fluctuation(data, wavelet='haar', order=1, max_sep=None, axis=0, lags='powers of 1.2', nan_behavior='raise', periodic=False)[source]

Mean absolute wavelet fluctuation of regularly spaced data along an axis.

At each scale (lag) r the data is convolved with the chosen wavelet’s kernel g_r and the q-th moment of the modulus is averaged over all positions (and over every other axis, treated as independent realizations):

F_q(r) = < |g_r * data|^q >.

For a self-affine field with Hurst exponent H this scales as F_q(r) ~ r^{xi(q)} with xi(q) = qH - K(q). All kernels share the L1 (“fluctuation”) normalization (sum|g_r| fixed across scales, mean removed), so the slope is the same physical exponent for every wavelet, while the amplitude (intercept) differs by wavelet. See wavelets for the kernel definitions and the lag->scale convention.

Parameters:

  • data (array-like) – N-dimensional scalar data array.

  • wavelet (str or Wavelet, optional) – Which wavelet to use (default: ‘haar’). One of the names in scaleinvariance.analysis.wavelets.WAVELETS: ‘haar’, ‘structure_function’, ‘mexican_hat’, ‘morlet’. An unknown name raises. (‘haar’ reproduces the previous haar_fluctuation exactly; ‘structure_function’ reproduces structure_function to numerical precision via the general path – the dedicated functions remain the faster route.)

  • order (float or 1-D array-like of float, optional) – Moment order(s) (default: 1). Each must be positive. With an array the modulus convolution |g_r * data| is computed once per lag and reused across all orders.

  • max_sep (int, optional) – Maximum lag to compute. If None, defaults to (size along axis) - 1.

  • axis (int, optional) – Axis along which to compute (default: 0).

  • lags (str, list, or np.ndarray, optional) – ‘all’, ‘powers of X’, or an explicit array of lags. Even-width wavelets (Haar) keep only even lags.

  • nan_behavior (str, optional) –

    ‘raise’ (default) errors on NaN input; ‘ignore’ propagates NaN through a NaN-safe convolution (a window containing any NaN yields NaN there), then the moment is a nanmean over the surviving positions. The rule is identical for every wavelet (and bit-identical to the old haar_fluctuation), but two consequences follow from the kernel widths (see wavelets):

    • A single NaN poisons a span as wide as the kernel: ~r for ‘haar’ / ‘structure_function’, but ~5.8 r (Mexican hat) and ~15 r (Morlet). Gappy data therefore loses far more output to the localized wavelets at the same lag, and large lags can go entirely NaN before r approaches the domain length.

    • wavelet='structure_function' would be stricter under NaN than the dedicated structure_function: the validity mask counts the full width-r+1 window (including the zero taps between the two +/-1 endpoints), so an interior NaN would poison the output even though it is multiplied by zero. To avoid silently discarding those pairs, when NaNs are present this case delegates to structure_function (which needs only the two endpoints finite) and emits a UserWarning saying so. The delegate reuses the same resolved lag array, so the only difference from the wavelet path is the (more correct) NaN handling.

  • periodic (bool, optional) – If True, treat the domain as periodic (toroidal) along axis and use circular convolution, giving all L wraparound windows per lag. Unlike the structure function there is no F(r) = F(L - r) symmetry, so the full lag range is kept; the constraint is simply that the kernel fit the domain.

Returns:

  • lags (np.ndarray) – 1-D int64 array of lag values.

  • fluctuations (np.ndarray) – Shape (n_lags,) for scalar order, else (n_orders, n_lags); the lag axis is last. Lags whose kernel exceeds the domain (the localized wavelets have kernels several times wider than their nominal scale) are returned as NaN.

scaleinvariance.analysis.wavelet_fluctuation.haar_fluctuation(data, order=1, max_sep=None, axis=0, lags='powers of 1.2', nan_behavior='raise', periodic=False)[source]

Mean absolute Haar fluctuation.

Thin wrapper for wavelet_fluctuation(..., wavelet='haar') (the Haar fluctuation is the difference of the means of the two halves of a width-r window). Kept for convenience and backward compatibility; see wavelet_fluctuation for the full parameter documentation.

scaleinvariance.analysis.wavelet_fluctuation.haar_fluctuation_analysis(*args, **kwargs)[source]

Deprecated: use haar_fluctuation() instead.

scaleinvariance.analysis.wavelet_fluctuation.haar_fluctuation_hurst(data, min_sep=None, max_sep=None, axis=0, return_fit=False, periodic=False)[source]

Hurst exponent from Haar-fluctuation scaling.

Thin wrapper for wavelet_fluctuation_hurst(..., wavelet='haar').

scaleinvariance.analysis.wavelet_fluctuation.wavelet_fluctuation(data, wavelet='haar', order=1, max_sep=None, axis=0, lags='powers of 1.2', nan_behavior='raise', periodic=False)[source]

Mean absolute wavelet fluctuation of regularly spaced data along an axis.

At each scale (lag) r the data is convolved with the chosen wavelet’s kernel g_r and the q-th moment of the modulus is averaged over all positions (and over every other axis, treated as independent realizations):

F_q(r) = < |g_r * data|^q >.

For a self-affine field with Hurst exponent H this scales as F_q(r) ~ r^{xi(q)} with xi(q) = qH - K(q). All kernels share the L1 (“fluctuation”) normalization (sum|g_r| fixed across scales, mean removed), so the slope is the same physical exponent for every wavelet, while the amplitude (intercept) differs by wavelet. See wavelets for the kernel definitions and the lag->scale convention.

Parameters:

  • data (array-like) – N-dimensional scalar data array.

  • wavelet (str or Wavelet, optional) – Which wavelet to use (default: ‘haar’). One of the names in scaleinvariance.analysis.wavelets.WAVELETS: ‘haar’, ‘structure_function’, ‘mexican_hat’, ‘morlet’. An unknown name raises. (‘haar’ reproduces the previous haar_fluctuation exactly; ‘structure_function’ reproduces structure_function to numerical precision via the general path – the dedicated functions remain the faster route.)

  • order (float or 1-D array-like of float, optional) – Moment order(s) (default: 1). Each must be positive. With an array the modulus convolution |g_r * data| is computed once per lag and reused across all orders.

  • max_sep (int, optional) – Maximum lag to compute. If None, defaults to (size along axis) - 1.

  • axis (int, optional) – Axis along which to compute (default: 0).

  • lags (str, list, or np.ndarray, optional) – ‘all’, ‘powers of X’, or an explicit array of lags. Even-width wavelets (Haar) keep only even lags.

  • nan_behavior (str, optional) –

    ‘raise’ (default) errors on NaN input; ‘ignore’ propagates NaN through a NaN-safe convolution (a window containing any NaN yields NaN there), then the moment is a nanmean over the surviving positions. The rule is identical for every wavelet (and bit-identical to the old haar_fluctuation), but two consequences follow from the kernel widths (see wavelets):

    • A single NaN poisons a span as wide as the kernel: ~r for ‘haar’ / ‘structure_function’, but ~5.8 r (Mexican hat) and ~15 r (Morlet). Gappy data therefore loses far more output to the localized wavelets at the same lag, and large lags can go entirely NaN before r approaches the domain length.

    • wavelet='structure_function' would be stricter under NaN than the dedicated structure_function: the validity mask counts the full width-r+1 window (including the zero taps between the two +/-1 endpoints), so an interior NaN would poison the output even though it is multiplied by zero. To avoid silently discarding those pairs, when NaNs are present this case delegates to structure_function (which needs only the two endpoints finite) and emits a UserWarning saying so. The delegate reuses the same resolved lag array, so the only difference from the wavelet path is the (more correct) NaN handling.

  • periodic (bool, optional) – If True, treat the domain as periodic (toroidal) along axis and use circular convolution, giving all L wraparound windows per lag. Unlike the structure function there is no F(r) = F(L - r) symmetry, so the full lag range is kept; the constraint is simply that the kernel fit the domain.

Returns:

  • lags (np.ndarray) – 1-D int64 array of lag values.

  • fluctuations (np.ndarray) – Shape (n_lags,) for scalar order, else (n_orders, n_lags); the lag axis is last. Lags whose kernel exceeds the domain (the localized wavelets have kernels several times wider than their nominal scale) are returned as NaN.

scaleinvariance.analysis.wavelet_fluctuation.wavelet_fluctuation_hurst(data, wavelet='haar', min_sep=None, max_sep=None, axis=0, return_fit=False, periodic=False)[source]

Estimate the Hurst exponent from first-order wavelet-fluctuation scaling.

For a self-similar process the order-1 fluctuation scales as F(r) ~ r^H; H is the slope of log F(r) vs log r. The accessible H range depends on the wavelet’s vanishing moments (Haar: roughly 0 < H < 1; the Mexican hat and Morlet reach higher H).

Parameters:

  • data (array-like) – N-dimensional scalar data array.

  • wavelet (str or Wavelet, optional) – Wavelet to use (default: ‘haar’).

  • min_sep (int, optional) – Lag window for the linear fit. Defaults match haar_fluctuation_hurst.

  • max_sep (int, optional) – Lag window for the linear fit. Defaults match haar_fluctuation_hurst.

  • axis (int, optional) – Axis along which to compute (default: 0).

  • return_fit (bool, optional) – If True, also return (lags, fluctuation values, fit_line).

  • periodic (bool, optional) – Periodic (toroidal) circular convolution along axis.

Return type:

H, uncertainty[, lags, values, fit_line]

Wavelet definitions for wavelet-fluctuation analysis.

A wavelet here is a discrete kernel g_r built at an integer scale (lag) r; convolving a field with g_r and taking |.| gives the fluctuation at scale r. For a self-affine field with Hurst exponent H, <|g_r * f|^q> ~ r^{xi(q)} with xi(q) = qH - K(q).

Normalization (the L1 / “fluctuation” convention)

Every kernel is built with the same L1 norm across scales, sum |g_r| = L1_NORM (= 2, matching the historical Haar and structure- function kernels), and with the discrete mean removed so sum g_r = 0 (one exact vanishing moment). Holding the L1 norm fixed is what makes the fluctuation scale as r^H rather than r^{H+1/2} (the L2 / energy convention used by the usual continuous wavelet transform). The absolute constant only shifts the intercept of the scaling, never the slope, so it is a convention, not a tuned parameter; different wavelets therefore share the same slope but not the same amplitude.

Scale convention (lag -> sigma)

The localized wavelets are anchored so that lag r equals the distance from the central peak to the first trough; at the smallest lag r = 1 that trough sits on the immediate grid neighbour. For the Mexican hat the trough is at sqrt(3) sigma so sigma = r / sqrt(3); for the Morlet (omega0 = 6) the first trough of the real part is at sigma * pi / omega0 so sigma = omega0 r / pi. Kernels are truncated at +/- 5 sigma for the Mexican hat (its polynomial-times-Gaussian tail is slightly heavier) and +/- 4 sigma for the Morlet.

Because a localized wavelet’s support is several sigma wide, its kernel is wider than its nominal scale (~5.8 r for the Mexican hat, ~15 r for the Morlet). A lag whose kernel would exceed the domain cannot be evaluated and is returned as NaN by wavelet_fluctuation – the largest reachable lag is correspondingly smaller for the wider wavelets.

The same kernel width governs NaN contamination under nan_behavior='ignore': an output position is NaN whenever its full kernel window contains any NaN, so a single NaN poisons a span as wide as the kernel (~r for Haar / structure function, ~5.8 r / ~15 r for the localized wavelets). The structure-function wavelet would therefore be stricter under NaN than the dedicated structure_function: the mask counts the whole width-r+1 window, including the zero taps between the two +/-1 endpoints, whereas the dedicated increment only needs the two endpoints finite. Because of this, wavelet_fluctuation auto-delegates to structure_function (with a warning) when the data contain NaNs and wavelet='structure_function'.

class scaleinvariance.analysis.wavelets.Wavelet(name: str, build_kernel: Callable[[int], object], even_only: bool, is_complex: bool, valid_H_range: Tuple[float, float], description: str)[source]

Bases: object

A wavelet usable for fluctuation analysis.

name
Type:

str

build_kernel

build_kernel(lag) -> backend 1-D array at integer scale lag.

Type:

callable

even_only

If True, only even lags are valid (the kernel needs an even width).

Type:

bool

is_complex

Whether the kernel is complex (the fluctuation is then the modulus of the complex convolution). Informational; the convolution backend detects complexity itself.

Type:

bool

valid_H_range

Approximate range of Hurst exponents the wavelet can resolve, set by its number of vanishing moments. Documentation only – not enforced.

Type:

tuple of float

description
Type:

str

build_kernel: Callable[[int], object]
description: str
even_only: bool
is_complex: bool
name: str
valid_H_range: Tuple[float, float]
scaleinvariance.analysis.wavelets.get_wavelet(name)[source]

Return the Wavelet registered under name (raises on unknown name).

scaleinvariance.analysis.structure_function.costructure_function(data1, data2, order1=1, order2=1, max_sep=None, axis=0, lags='powers of 1.2', periodic=False)[source]

Calculate the mixed (two-field, two-exponent) structure function.

Defined as:

S_{p,q}(r) = <|u_1(x + r) - u_1(x)|^p * |u_2(x + r) - u_2(x)|^q>

where <> denotes the average over positions where neither increment is NaN. (Matches the NaN-only masking convention of structure_function and haar_fluctuation; inf values are not specially handled.)

The diagonal p = q case recovers a single-exponent co-structure function; the p -> 0 or q -> 0 slices recover the single-field structure functions of the other field (up to the caveat that 0 is not a permitted order here).

Parameters:

  • data1 (array-like) – N-dimensional scalar data arrays. Must have identical shape.

  • data2 (array-like) – N-dimensional scalar data arrays. Must have identical shape.

  • order1 (float or 1-D array-like of float, optional) – Exponent(s) applied to |delta f_1| and |delta f_2| respectively (default: 1). Each order must be positive. Arrays are allowed on either or both; the absolute increments of each field are computed once per lag and reused across all (order1, order2) combinations.

  • order2 (float or 1-D array-like of float, optional) – Exponent(s) applied to |delta f_1| and |delta f_2| respectively (default: 1). Each order must be positive. Arrays are allowed on either or both; the absolute increments of each field are computed once per lag and reused across all (order1, order2) combinations.

  • max_sep (int, optional) – Maximum lag to compute. If None, defaults to (size along axis) - 1.

  • axis (int, optional) – Axis along which to compute increments (default: 0).

  • lags (str, list, or np.ndarray, optional) – Lag selection (see structure_function for accepted forms).

  • periodic (bool, optional) – If True, treat the domain as periodic (toroidal) along axis: both increments wrap around the array end (all L pairs per lag). max_sep is capped at L // 2 and an explicit max_sep > L // 2 raises ValueError (see structure_function).

Returns:

lags : 1-D array of lag values. cosf : co-structure function values.

Shape depends on order1/order2 being scalar vs array: - (n_lags,) both scalar - (n_order1, n_lags) order1 array, order2 scalar - (n_order2, n_lags) order1 scalar, order2 array - (n_order1, n_order2, n_lags) both arrays The lag axis is always last.

Return type:

tuple (np.ndarray, np.ndarray)

Notes

The average is over positions where neither increment is NaN. This is asymmetric with structure_function, which averages over positions where the single field’s increment is non-NaN. If one field has NaNs and the other does not, the effective sample count here is the intersection of valid positions.

inf values are not specially handled and will propagate through to the result (consistent with the NaN-only masking of the single-field structure_function).

scaleinvariance.analysis.structure_function.structure_function(data, order=1, max_sep=None, axis=0, lags='powers of 1.2', periodic=False)[source]

Calculate the mean structure function for regularly spaced scalar data along a given axis.

The structure function S_n(r) of order n is defined as:

S_n(r) = <|u(x + r) - u(x)|^n>

where <> denotes the average over all valid pairs separated by lag r.

Parameters:

  • data (array-like) – N-dimensional scalar data array.

  • order (float or 1-D array-like of float, optional) – Order(s) of the structure function (default: 1). Each order must be positive. When an array is given, the absolute increments |u(x+r) - u(x)| are computed once per lag and reused across all orders (avoids redundant work).

  • max_sep (int, optional) – Maximum lag to compute. If None, defaults to (size along axis) - 1, or (size along axis) // 2 when periodic=True.

  • axis (int, optional) – Axis along which to compute the structure function (default: 0).

  • lags (str, list, or np.ndarray, optional) –

    Option for selecting lag values. Can be:

    • A list or array of lags.

    • ’all’: all integer lags from 1 to max_sep.

    • ’powers of X’: lags as integers of powers of X, where X is any number (e.g., ‘powers of 2’, ‘powers of 1.2’).

  • periodic (bool, optional) – If True, treat the domain as periodic (toroidal) along axis: the increment at lag r uses all L pairs u((x + r) mod L) - u(x) rather than dropping the L - r edge pairs. On a periodic domain S(r) = S(L - r) exactly, so lags above L // 2 are redundant; max_sep is therefore capped at L // 2 and an explicit max_sep > L // 2 raises ValueError.

Returns:

lags : 1-D array of lag values. sf : structure function values.

Shape is (n_lags,) if order is a scalar, or (n_orders, n_lags) if order is a 1-D array. The lag axis is always last.

Return type:

tuple (np.ndarray, np.ndarray)

scaleinvariance.analysis.structure_function.structure_function_analysis(*args, **kwargs)[source]

Deprecated: use structure_function() instead.

scaleinvariance.analysis.structure_function.structure_function_hurst(data, min_sep=None, max_sep=None, axis=0, return_fit=False, periodic=False)[source]

Calculate the Hurst exponent using structure functions.

For a self-similar process with Hurst exponent H, the first order structure function scales as: S_1(r) ~ r^H

This function estimates H by performing a linear regression on log(S_1(r)) vs log(r).

Parameters:

  • data (array-like) – N-dimensional scalar data array.

  • min_sep (int, optional) – Minimum separation/lag to include in the fit. If None, defaults to 1.

  • max_sep (int, optional) – Maximum separation/lag to include in the fit. If None, defaults to array_size/8 or array_size if array_size < 128 (with warning).

  • axis (int, optional) – Axis along which to compute the structure function (default: 0).

  • return_fit (bool, optional) – If True, return the lags, structure function values, and fitted line.

  • periodic (bool, optional) – If True, use periodic (toroidal) increments along axis (see structure_function). The default max_sep becomes array_size // 2.

Returns:

  • If return_fit is False –

    Hfloat

    Estimated Hurst exponent.

    uncertaintyfloat

    Standard error of the Hurst exponent estimate.

  • If return_fit is True –

    Hfloat

    Estimated Hurst exponent.

    uncertaintyfloat

    Standard error of the Hurst exponent estimate.

    lagsnp.ndarray

    The lags used in the calculation.

    sf_valuesnp.ndarray

    The structure function values.

    fit_linenp.ndarray

    The fitted line used to estimate H.

scaleinvariance.analysis.spectral.power_spectrum_binned(data, max_wavelength=None, min_wavelength=None, nbins=50, axis=0)[source]

Compute binned power spectral density for data along specified axis.

Note: this estimator is FFT-based and therefore inherently treats the domain as periodic, so there is no periodic flag (unlike structure_function / haar_fluctuation).

Parameters:

  • data (array-like) – N-dimensional input data array

  • max_wavelength (float, optional) – Maximum wavelength to include in analysis

  • min_wavelength (float, optional) – Minimum wavelength to include in analysis

  • nbins (int, optional) – Number of logarithmic bins for PSD (default: 50)

  • axis (int, optional) – Axis along which to compute FFT (default: 0)

Returns:

binned_freq : 1-D array of binned frequencies binned_psd : 1-D array of binned power spectral density values

Return type:

tuple (np.ndarray, np.ndarray)

scaleinvariance.analysis.spectral.spectral_analysis(*args, **kwargs)[source]

Deprecated: use power_spectrum_binned() instead.

scaleinvariance.analysis.spectral.spectral_hurst(data, max_wavelength=None, min_wavelength=None, nbins=50, axis=0, return_fit=False)[source]

Estimate the Hurst exponent from the power spectral density along specified axis.

For fractional Brownian motion, PSD ~ f^(-β) with β = 2H + 1. This function estimates H by performing linear regression on log(PSD) vs log(frequency).

Note: this estimator is FFT-based and therefore inherently treats the domain as periodic, so there is no periodic flag (unlike structure_function_hurst / haar_fluctuation_hurst).

Parameters:

  • data (array-like) – N-dimensional input data array

  • max_wavelength (float, optional) – Maximum wavelength to include in analysis. If None, defaults to domain_size/8 or domain_size if domain_size < 512 (with warning).

  • min_wavelength (float, optional) – Minimum wavelength to include in analysis. If None, defaults to 4*Nyquist or Nyquist if domain_size < 512.

  • nbins (int, optional) – Number of logarithmic bins for PSD (default: 50)

  • axis (int, optional) – Axis along which to compute FFT (default: 0)

  • return_fit (bool, optional) – If True, return frequencies, PSD values, and fitted line

Returns:

  • If return_fit is False –

    Hfloat

    Estimated Hurst exponent

    uncertaintyfloat

    Standard error of the Hurst exponent estimate

  • If return_fit is True –

    Hfloat

    Estimated Hurst exponent

    uncertaintyfloat

    Standard error of the Hurst exponent estimate

    freqsnp.ndarray

    The frequencies used in the calculation

    psd_valuesnp.ndarray

    The power spectral density values

    fit_linenp.ndarray

    The fitted line used to estimate H

Utilities

Utility functions for scaleinvariance package

Common utilities used across analysis and simulation modules.

scaleinvariance.utils.estimate_hurst_from_scaling(lags, values, min_sep=None, max_sep=None, return_fit=False)[source]

Estimate Hurst exponent from scaling relationship.

Performs linear regression on log(values) vs log(lags) to estimate scaling exponent.

Parameters:

  • lags (array-like) – Lag values.

  • values (array-like) – Corresponding values (structure function or fluctuation values).

  • min_sep (int, optional) – Minimum separation/lag to include in the fit.

  • max_sep (int, optional) – Maximum separation/lag to include in the fit.

  • return_fit (bool, optional) – If True, return additional fit information.

Returns:

  • H (float) – Estimated Hurst exponent.

  • uncertainty (float) – Standard error of the Hurst exponent estimate.

  • lags_fit (np.ndarray (if return_fit=True)) – The lags used in the fit.

  • values_fit (np.ndarray (if return_fit=True)) – The values used in the fit.

  • fit_line (np.ndarray (if return_fit=True)) – The fitted line.

scaleinvariance.utils.process_lags(lags, max_sep, even_only=False)[source]

Process lag specification into array of lag values.

Parameters:

  • lags (str, list, or np.ndarray) –

    Option for selecting lag values. Can be:

    • A list or array of lags.

    • ’all’: all integer lags from 1 to max_sep.

    • ’powers of XXXX’: lags as closest intergers to powers of XXXX, e.g. 2 (~ 3 per decade) or 1.2 (~ 8 per decade).

  • max_sep (int) – Maximum separation/lag value.

  • even_only (bool, optional) – If True, filter to only even lags (default: False).

Returns:

Array of processed lag values.

Return type:

np.ndarray