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]

Calculate the Hurst exponent using Haar fluctuation analysis.

For a self-similar process with Hurst exponent H, the Haar fluctuation scales as: F_H(r) ~ r^H

This function estimates H by performing a linear regression on log(F_H(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.

  • max_sep (int, optional) – Maximum separation/lag to include in the fit. Defaults to array_size - 1 (a warning is emitted for small arrays).

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

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

  • periodic (bool, optional) – If True, use periodic (toroidal) circular convolution along axis (see haar_fluctuation). Unlike the structure function, no L // 2 cap applies; the full lag range is retained.

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.

    haar_valuesnp.ndarray

    The Haar fluctuation values.

    fit_linenp.ndarray

    The fitted line used to estimate H.

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.analysis.haar_fluctuation.haar_fluctuation(data, order=1, max_sep=None, axis=0, lags='powers of 1.2', nan_behavior='raise', periodic=False)[source]

Calculate the mean absolute Haar fluctuation for regularly spaced scalar data along a given axis.

Uses FFT-based convolution via the backend for efficient computation, with a direct fallback when NaNs are present.

Parameters:

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

  • order (float or 1-D array-like of float, optional) – Order(s) of the fluctuation (default: 1). Each order must be positive. When an array is given, the absolute Haar convolution |conv(data, kernel)| is computed once per lag and reused across all orders (avoids redundant convolutions).

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

  • axis (int, optional) – Axis along which to compute the Haar fluctuation (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’).

  • nan_behavior (str, optional) –

    Behavior when encountering NaN values in the data. Options are:

    • ’raise’: Raise an error if NaN values are found (default).

    • ’ignore’: Ignore NaN values in the calculation.

  • periodic (bool, optional) – If True, treat the domain as periodic (toroidal) along axis: the Haar fluctuation is computed by circular convolution, giving all L wraparound windows per lag rather than dropping the L - r + 1 edge windows. Unlike the structure function, the Haar fluctuation has no F(r) = F(L - r) reflection symmetry (its half-window width r/2 differs from (L - r)/2), so the full lag range up to the array size is kept; the only constraint is r <= L (the kernel must fit the domain). NaN handling under nan_behavior='ignore' is respected (a wraparound window containing any NaN yields NaN at that position).

Returns:

lags : 1-D array of lag values. haar : mean absolute haar fluctuation 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.haar_fluctuation.haar_fluctuation_analysis(*args, **kwargs)[source]

Deprecated: use haar_fluctuation() instead.

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

Calculate the Hurst exponent using Haar fluctuation analysis.

For a self-similar process with Hurst exponent H, the Haar fluctuation scales as: F_H(r) ~ r^H

This function estimates H by performing a linear regression on log(F_H(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.

  • max_sep (int, optional) – Maximum separation/lag to include in the fit. Defaults to array_size - 1 (a warning is emitted for small arrays).

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

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

  • periodic (bool, optional) – If True, use periodic (toroidal) circular convolution along axis (see haar_fluctuation). Unlike the structure function, no L // 2 cap applies; the full lag range is retained.

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.

    haar_valuesnp.ndarray

    The Haar fluctuation values.

    fit_linenp.ndarray

    The fitted line used to estimate H.

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