feat: add parallel uncertainty fraction mapping

src/pynamicalsys/common/basin_analysis.py

@@ -19,6 +19,7 @@
 from numba import njit, prange
 from typing import Tuple
 from numpy.typing import NDArray
+from joblib import delayed, Parallel
 
 
 def uncertainty_fraction(x, y, basin, epsilon_max, n_eps=100, epsilon_min=0):
@@ -168,3 +169,107 @@ def basin_entropy(
     Sbb = np.mean(S[boundary_mask]) if boundary_mask.any() else 0.0
 
     return float(Sb), float(Sbb)
+
+def uncertainty_fraction_mapping(
+    X: NDArray[np.float64],
+    Y: NDArray[np.float64],
+    basin: NDArray[np.float64],
+    mapping: object,
+    parameters: object,
+    exits: object,
+    escape: str = "exiting",
+    n_samples: int = 120_000,
+    p_samples: int = 7,
+    threshold: float = 0.1,
+    n_eps: int = 100,
+    max_time: int = 1000,
+    seed: int = 13,
+    n_jobs: int = -1,
+) -> Tuple[NDArray[np.float64], NDArray[np.float64]]:
+    """
+    Wrapper to compute the uncertainty fraction f(epsilon) for a given basin of attraction with mapping.
+    """
+
+    rng = np.random.default_rng(seed)
+
+    flat_x = X.flatten()
+    flat_y = Y.flatten()
+    flat_basin = basin.flatten()
+
+    epsilons = np.logspace(np.log10(0.1), np.log10(1 / X.size), n_eps)
+
+    idx = rng.choice(len(flat_x), size=n_samples, replace=False)
+    xs = flat_x[idx]
+    ys = flat_y[idx]
+    labels = flat_basin[idx]
+
+    angles = np.linspace(0, 2 * np.pi, p_samples, endpoint=False)
+    dcos = np.cos(angles)
+    dsin = np.sin(angles)
+
+    needed = max(1, int(np.ceil(threshold * p_samples)))
+
+    ss = np.random.SeedSequence(seed)
+    child_seeds = ss.spawn(n_eps)
+
+    f_eps = Parallel(n_jobs=n_jobs)(
+        delayed(_process_single_epsilon)(
+            eps, xs, ys, labels, p_samples, dcos, dsin, needed,
+            mapping, max_time, parameters, exits, escape, n_samples,
+            child_seeds[i],
+        )
+        for i, eps in enumerate(epsilons)
+    )
+
+    return epsilons, f_eps
+
+
+def _process_single_epsilon(
+    eps: float,
+    xs: NDArray[np.float64],
+    ys: NDArray[np.float64],
+    labels: NDArray[np.float64],
+    p_samples: int,
+    dcos: NDArray[np.float64],
+    dsin: NDArray[np.float64],
+    needed: int,
+    mapping: object,
+    max_time: int,
+    parameters: object,
+    exits: object,
+    escape: str,
+    n_samples: int,
+    child_seed: np.random.SeedSequence,
+) -> float:
+    """Worker function that computes the uncertainty fraction for a single
+    epsilon value.
+    """
+
+    rng = np.random.default_rng(child_seed)
+    uncertain = 0
+
+    for x, y, side_center in zip(xs, ys, labels):
+        u_rand = rng.random(p_samples)
+        r = eps * np.sqrt(u_rand)
+        x_eps = x + r * dcos
+        y_eps = y + r * dsin
+
+        different = 0
+        for xp, yp in zip(x_eps, y_eps):
+            side_p, _ = mapping.escape_analysis(
+                u=(xp, yp),
+                max_time=max_time,
+                parameters=parameters,
+                exits=exits,
+                escape=escape,
+            )
+            if side_p != side_center:
+                different += 1
+                if different >= needed:
+                    break
+
+        if different >= needed:
+            uncertain += 1
+
+    return uncertain / n_samples
+

src/pynamicalsys/core/basin_metrics.py

@@ -19,7 +19,7 @@
 from numbers import Integral, Real
 from typing import Optional, Tuple
 from numpy.typing import NDArray
-from pynamicalsys.common.basin_analysis import basin_entropy, uncertainty_fraction
+from pynamicalsys.common.basin_analysis import basin_entropy, uncertainty_fraction, uncertainty_fraction_mapping
 
 
 class BasinMetrics:
@@ -204,3 +204,207 @@ def uncertainty_fraction(
             n_eps=n_eps,
             epsilon_min=epsilon_min,
         )
+
+    def uncertainty_fraction_mapping(
+            self,
+            X,
+            Y,
+            mapping,
+            parameters,
+            exits,
+            escape='exiting',
+            n_samples=120_000,
+            p_samples=7,
+            threshold=0.1,
+            n_eps = 100,
+            max_time=1000,
+            seed=13,
+            n_jobs=-1):
+        """
+        Estimate the uncertainty fraction of a dynamical mapping using Monte Carlo sampling.
+
+        The computation is performed through random sampling of initial conditions and
+        perturbations of size ε, allowing the estimation of the scaling law:
+
+        f(ε) ~ ε^{⍺},
+
+        where ⍺ is the uncertainty exponent. For a D-dimensional phase space, the
+        dimension d of the basin boundary is related to ⍺ by:
+
+        d = D - ⍺.
+
+        Parameters
+        ----------
+        X : NDArray[np.float64]
+            2D array containing the x-coordinates of grid.
+        Y : NDArray[np.float64]
+            2D array containing the y-coordinates of grid.
+        mapping : callable
+            Dynamical mapping function describing the discrete-time system evolution.
+            The function must accept the current state and the system parameters as input.
+        parameters : list
+            Parameters passed to the mapping function.
+        exits : list
+            List containing the exit regions or exit conditions of the system.
+        escape : str, optional
+            Escape criterion type (default: ``"exiting"``).
+
+            Currently, only:
+            - ``"exiting"`` : trajectories are classified according to the exit reached.
+            is implemented.
+        n_samples : int, optional
+            Number of random initial conditions sampled for the Monte Carlo estimation
+            (default: 120000).
+        p_samples : int, optional
+            Number of perturbed neighbors generated for each sampled point
+            (default: 7).
+        threshold : float, optional
+            Fraction threshold used to classify a point as uncertain
+            (default: 0.1).
+
+            Must satisfy:
+
+            0 <= threshold <= 1
+        n_eps : int, optional
+            Number of epsilon values used in the uncertainty scaling analysis
+            (default: 100).
+        max_time : int, optional
+            Maximum number of iterations allowed for trajectory evolution
+            (default: 1000).
+        seed : int, optional
+            Seed for the random number generator used during sampling
+            (default: 13).
+        n_jobs : int, optional
+            Number of parallel jobs used during computation
+            (default: -1).
+
+            Common values:
+            - ``-1`` : use all available CPU cores
+
+        Returns
+        -------
+        Tuple[NDArray[np.float64], NDArray[np.float64]]
+            A tuple containing:
+
+            - epsilons : Array of epsilon values.
+            - uncertainty_fraction : Array containing the estimated uncertainty
+            fraction corresponding to each epsilon value.
+
+        Raises
+        ------
+        ValueError
+            If:
+
+            - ``X`` or ``Y`` are not 2-dimensional arrays.
+            - ``X``, ``Y``, and ``basin`` do not have the same shape.
+            - ``parameters`` is ``None``.
+            - ``escape`` is not ``"exiting"``.
+            - ``n_samples`` is not a positive integer.
+            - ``p_samples`` is not a positive integer.
+            - ``n_eps`` is not a positive integer.
+            - ``max_time`` is not a positive integer.
+            - ``seed`` is negative.
+            - ``threshold`` is not between 0 and 1.
+            - ``n_jobs`` is zero.
+
+        Notes
+        -----
+        - This implementation uses Monte Carlo sampling to reduce computational cost
+        when dealing with large phase spaces.
+        - Parallel execution is supported through the ``n_jobs`` parameter.
+
+        Examples
+        --------
+        >>> import numpy as np
+        >>> from numba import njit
+        >>> from joblib import Parallel, delayed
+        >>> from pynamicalsys import BasinMetrics, DiscreteDynamicalSystem as dds
+        >>>
+        >>> @njit
+        ... def henon_map(state, parameters):
+        ...     x, y = state
+        ...     a, b = parameters
+        ...     return np.array([1 - a*x**2 + y, b*x])
+        >>>
+        >>> henon = dds(
+        ...     mapping=henon_map,
+        ...     number_of_parameters=2,
+        ...     system_dimension=2,
+        ... )
+        >>>
+        >>> grid_size = 1000
+        >>> x = np.linspace(-2, 2, grid_size)
+        >>> y = np.linspace(-2, 2, grid_size)
+        >>> X, Y = np.meshgrid(x, y, indexing='ij')
+        >>> grid_points = np.column_stack((X.ravel(), Y.ravel()))
+        >>>
+        >>> exits = np.array([[-10, 10], [-10, 10]], dtype=np.float64)
+        >>> N = 2000
+        >>>
+        >>> escape = np.array(Parallel(n_jobs=-1)(
+        ...     delayed(henon.escape_analysis)(
+        ...         u, N, exits, parameters=[1.45, 0.3], escape="exiting"
+        ...     )
+        ...     for u in grid_points
+        ... ))
+        >>>
+        >>> basin = escape[:, 0].reshape(grid_size, grid_size)
+        >>>
+        >>> metrics = BasinMetrics(basin)
+        >>>
+        >>> epsilons, f = metrics.uncertainty_fraction_mapping(
+        ...     X,
+        ...     Y,
+        ...     mapping=henon,
+        ...     parameters=[1.40, 0.3],
+        ...     exits=exits,
+        ...     n_samples=10_000,
+        ...     p_samples=7,
+        ...     max_time=1000,
+        ... )
+        """
+
+        X = np.asarray(X, dtype=np.float64)
+        Y = np.asarray(Y, dtype=np.float64)
+
+        if X.ndim != 2 or Y.ndim != 2:
+            raise ValueError("X, Y, and basin must be 2-dimensional arrays")
+        if X.shape != Y.shape or X.shape != self.basin.shape:
+            raise ValueError("X, Y, and basin must have the same shape")
+        if parameters is None:
+            raise ValueError("parameters cannot be None")
+        if escape != "exiting":
+            raise ValueError("escape must be 'exiting'. Option 'entering' is not implemented yet")
+        if not isinstance(n_samples, int) or n_samples <= 0:
+            raise ValueError("n_samples must be a positive integer")
+        if not isinstance(p_samples, int) or p_samples <= 0:
+            raise ValueError("p_samples must be a positive integer")
+        if not isinstance(n_eps, int) or n_eps <= 0:
+            raise ValueError("n_eps must be a positive integer")
+        if not isinstance(max_time, int) or max_time <= 0:
+            raise ValueError("max_time must be a positive integer")
+        if not isinstance(seed, int) or seed < 0:
+            raise ValueError("seed must be a non-negative integer")
+        if not isinstance(threshold, (int, float)) or not (0 <= threshold <= 1):
+            raise ValueError("threshold must be a float between 0 and 1")
+        if n_jobs == 0:
+            raise ValueError("n_jobs cannot be zero")
+
+
+        return uncertainty_fraction_mapping(
+            X,
+            Y,
+            self.basin,
+            mapping,
+            parameters,
+            exits,
+            escape,
+            n_samples,
+            p_samples,
+            threshold,
+            n_eps,
+            max_time,
+            seed,
+            n_jobs
+        )
+