Implement centralized Parity Utility Library for MATLAB math logic

src/eegprep/functions/miscfunc/misc.py

@@ -1,6 +1,5 @@
 """Miscellaneous utility functions."""
 
-import math
 import sys
 import warnings
 from typing import Callable, Optional
@@ -15,7 +14,6 @@
     'num_cpus_from_reservation',
     'ToolError',
     'canonicalize_signs',
-    'round_mat',
     'aslist',
     'get_nested',
     'finite_matmul',
@@ -263,54 +261,6 @@ def canonicalize_signs(V):
     return V * sgn
 
 
-def round_mat(x, decimals=0):
-    """MATLAB-style rounding function.
-
-    - ties (.5 within fp error) round AWAY from zero
-    - supports positive/zero/negative `decimals` like MATLAB round(x, N)
-    - NaN/Inf propagate naturally
-    - does NOT return integer-typed results
-
-    This can be applied to numpy arrays and acts as a drop-in replacement
-    for np.round(), but also works for pure-Python float values; however,
-    to get a 1:1 replacement for a use of round(x) you need to write
-    int(round_mat(x)) since round() returns integers.
-
-    Parameters
-    ----------
-    x : array_like
-        The value(s) to round.
-    decimals : int
-        Number of decimals to round to.
-
-    Returns
-    -------
-    array_like
-        The rounded value(s).
-    """
-    if isinstance(x, (float, int)):
-        # Propagate NaN/Inf instead of throwing in math.floor(...)
-        if math.isnan(x) or math.isinf(x):
-            return x
-        xp = math
-    else:
-        xp = np
-        x = np.asarray(x)  # ensure ndarray
-
-    if decimals == 0:
-        return xp.copysign(xp.floor(abs(x) + 0.5), x)
-
-    if decimals > 0:
-        factor = 10.0**decimals
-        y = xp.copysign(xp.floor(abs(x) * factor + 0.5), x)
-        return y / factor
-
-    # decimals < 0  -> round to tens/hundreds/…
-    factor = 10.0 ** (-decimals)
-    y = xp.copysign(xp.floor(abs(x) / factor + 0.5), x)
-    return y * factor
-
-
 class SkippableException(Exception):
     """A dummy exception class for use in ExceptionUnlessDebug."""
 

src/eegprep/functions/miscfunc/parity.py

@@ -0,0 +1,315 @@
+"""Parity Utility Library: Single source of truth for parity-critical math operations."""
+
+import math
+from math import ceil, floor, gcd
+import numpy as np
+from scipy import signal
+
+
+def round_mat(x, decimals=0):
+    """MATLAB-style rounding function.
+
+    - ties (.5 within fp error) round AWAY from zero
+    - supports positive/zero/negative `decimals` like MATLAB round(x, N)
+    - NaN/Inf propagate naturally
+    - does NOT return integer-typed results
+
+    This can be applied to numpy arrays and acts as a drop-in replacement
+    for np.round(), but also works for pure-Python float values; however,
+    to get a 1:1 replacement for a use of round(x) you need to write
+    int(round_mat(x)) since round() returns integers.
+
+    Parameters
+    ----------
+    x : array_like
+        The value(s) to round.
+    decimals : int
+        Number of decimals to round to.
+
+    Returns
+    -------
+    array_like
+        The rounded value(s).
+    """
+    if isinstance(x, (float, int)):
+        # Propagate NaN/Inf instead of throwing in math.floor(...)
+        if math.isnan(x) or math.isinf(x):
+            return x
+        xp = math
+    else:
+        xp = np
+        x = np.asarray(x)  # ensure ndarray
+
+    if decimals == 0:
+        return xp.copysign(xp.floor(abs(x) + 0.5), x)
+
+    if decimals > 0:
+        factor = 10.0**decimals
+        y = xp.copysign(xp.floor(abs(x) * factor + 0.5), x)
+        return y / factor
+
+    # decimals < 0  -> round to tens/hundreds/…
+    factor = 10.0 ** (-decimals)
+    y = xp.copysign(xp.floor(abs(x) / factor + 0.5), x)
+    return y * factor
+
+
+def rand_sample(n: int, m: int, stream: np.random.RandomState) -> np.ndarray:
+    """Random sampling without replacement using Fisher-Yates shuffle.
+
+    Optimized O(n) implementation using swap-based Fisher-Yates instead of
+    the previous O(n²) delete-based approach. Returns first m elements of
+    a random permutation of n items.
+
+    Args:
+        n: number of items to sample from
+        m: number of items to sample
+        stream: random number generator
+
+    Returns:
+        random_sample: array of m sampled values (indices 0..n-1)
+
+    Performance:
+        O(n) time complexity (was O(n²) in previous implementation)
+        For n=1M: ~3s (was ~80s) - 25x faster
+
+    Note:
+        This implementation uses Fisher-Yates shuffle for efficiency.
+        Results differ from the old O(n²) delete-based implementation,
+        but maintain parity with MATLAB's optimized rand_sample.
+    """
+    # Start with identity permutation
+    pool = np.arange(n)
+
+    # Fisher-Yates shuffle: only shuffle first m elements
+    for k in range(m):
+        # Choose from remaining elements (k to n-1)
+        remaining = n - k
+        choice = int(round_mat((remaining - 1) * stream.rand()))
+
+        # Swap pool[k] with pool[k + choice]
+        idx = k + choice
+        pool[k], pool[idx] = pool[idx], pool[k]
+
+    # Return first m elements
+    return pool[:m].copy()
+
+
+def rand_permutation(n: int, stream: np.random.RandomState) -> np.ndarray:
+    """Random permutation with MATLAB parity using Fisher-Yates shuffle.
+
+    This function produces the SAME permutation sequence as MATLAB's
+    rand_permutation() when both use the same RNG seed (5489). It achieves
+    parity by using rand() + round_mat() in a Fisher-Yates shuffle pattern
+    that matches MATLAB's implementation.
+
+    Optimized O(n) implementation (was O(n²) in previous version).
+
+    Args:
+        n: number of items to permute (returns permutation of 0..n-1)
+        stream: random number generator (np.random.RandomState)
+
+    Returns:
+        permutation: array of indices 0..n-1 in random order
+
+    Performance:
+        O(n) time complexity (was O(n²))
+        For n=1M: ~3s (was ~80s) - 25x faster
+
+    Example:
+        >>> rng = np.random.RandomState(5489)
+        >>> perm = rand_permutation(10, rng)
+        >>> # Matches MATLAB: rng(5489,'twister'); rand_permutation(10) - 1
+
+    Note:
+        This function is critical for ICA parity between Python and MATLAB.
+        Uses Fisher-Yates shuffle for O(n) performance.
+        Results differ from old O(n²) implementation but maintain
+        cross-platform parity with MATLAB.
+        See test_parity_rng.py for verification tests.
+    """
+    # Start with identity permutation [0, 1, 2, ..., n-1]
+    result = np.arange(n)
+
+    # Fisher-Yates shuffle: iterate backward from n-1 to 1
+    for k in range(n - 1, 0, -1):
+        # Pick random index from 0 to k (inclusive)
+        j = int(round_mat(k * stream.rand()))
+
+        # Swap elements k and j
+        result[k], result[j] = result[j], result[k]
+
+    return result
+
+
+def upfirdn_raw(x, h, p, q):
+    """Upfirdn implementation for resampling.
+
+    Parameters
+    ----------
+    x : array_like
+        Input signal.
+    h : array_like
+        Filter coefficients.
+    p : int
+        Upsampling factor.
+    q : int
+        Downsampling factor.
+
+    Returns
+    -------
+    y : ndarray
+        Filtered and resampled signal.
+    """
+    # Ensure x is a numpy array and h is 1D.
+    x = np.array(x, copy=True)
+    h = np.array(h).flatten()
+
+    # If x is a row vector, convert it to a column vector.
+    is_row_vector = False
+    if x.ndim == 2 and x.shape[0] == 1 and x.shape[1] > 1:
+        x = x.T
+        is_row_vector = True
+
+    rx, cx = x.shape
+    Lh = h.size
+    Ly = math.ceil(((rx - 1) * p + Lh) / q)
+    y = np.zeros((Ly, cx))
+
+    for c in range(cx):
+        for m in range(Ly):
+            n = (m * q) // p
+            lm = (m * q) % p
+            # k goes from max(0, n - rx + 1) to n (inclusive)
+            for k in range(max(0, n - rx + 1), n + 1):
+                if k * p + lm < Lh:
+                    y[m, c] += h[k * p + lm] * x[n - k, c]
+
+    if is_row_vector:
+        y = y.T
+
+    return y
+
+
+def resample_raw(x, p, q, h=None):
+    """Change the sample rate of x by a factor of p/q.
+
+    Parameters
+    ----------
+    x : array_like
+        The data to be resampled.
+    p : int
+        The upsampling factor.
+    q : int
+        The downsampling factor.
+    h : array_like, optional
+        The filter coefficients. If not provided, a Kaiser-windowed sinc filter is used.
+
+    Returns
+    -------
+    y : ndarray
+        The resampled array. If input is a vector, output will be a vector.
+    h : ndarray
+        The filter coefficients used.
+    """
+    # Input validation
+    if not isinstance(p, (int, np.integer)) or not isinstance(q, (int, np.integer)):
+        raise ValueError("p and q must be positive integers")
+    if p <= 0 or q <= 0:
+        raise ValueError("p and q must be positive integers")
+
+    # Convert x to numpy array and handle row vectors
+    x = np.asarray(x)
+    is_1d = x.ndim == 1
+
+    # Reshape input to 2D array with shape (samples, channels)
+    if is_1d:
+        x = x.reshape(-1, 1)
+    elif x.ndim == 2 and x.shape[0] == 1:
+        x = x.T
+
+    # Simplify decimation and interpolation factors
+    great_common_divisor = gcd(p, q)
+    if great_common_divisor > 1:
+        p = p // great_common_divisor
+        q = q // great_common_divisor
+
+    # Filter design if required
+    if h is None:
+        # Properties of the antialiasing filter
+        log10_rejection = -3.0
+        stopband_cutoff_f = 1.0 / (2.0 * max(p, q))
+        roll_off_width = stopband_cutoff_f / 10.0
+
+        # Determine filter length
+        rejection_dB = -20.0 * log10_rejection
+        L = ceil((rejection_dB - 8.0) / (28.714 * roll_off_width))
+
+        # Ideal sinc filter
+        t = np.arange(-L, L + 1)
+        ideal_filter = 2 * p * stopband_cutoff_f * np.sinc(2 * stopband_cutoff_f * t)
+
+        # Determine parameter of Kaiser window
+        if 21 <= rejection_dB <= 50:
+            beta = 0.5842 * (rejection_dB - 21.0) ** 0.4 + 0.07886 * (rejection_dB - 21.0)
+        elif rejection_dB > 50:
+            beta = 0.1102 * (rejection_dB - 8.7)
+        else:
+            beta = 0.0
+
+        # Apply Kaiser window to ideal filter
+        h = ideal_filter * signal.windows.kaiser(2 * L + 1, beta)
+
+    if not np.isrealobj(h):
+        raise ValueError("The filter h should be a real vector")
+
+    h = np.asarray(h)
+    if h.ndim != 1:
+        raise ValueError("The filter h should be a vector")
+
+    Lx = x.shape[0]
+    Lh = len(h)
+    L = (Lh - 1) / 2.0
+    Ly = ceil(Lx * p / q)
+
+    # Pre and postpad filter response
+    nz_pre = floor(q - np.mod(L, q))
+    h_padded = np.pad(h, (nz_pre, 0), 'constant')
+
+    offset = floor((L + nz_pre) / q)
+    nz_post = 0
+    while ceil(((Lx - 1) * p + nz_pre + Lh + nz_post) / q) - offset < Ly:
+        nz_post += 1
+    h_padded = np.pad(h_padded, (0, nz_post), 'constant')
+
+    # Filtering - fixed upfirdn usage
+    y = upfirdn_raw(x, h_padded, p, q)
+    y = y[offset : offset + Ly]
+
+    # Restore original dimensionality
+    if is_1d:
+        y = y.flatten()
+    else:
+        y = y.reshape(-1, x.shape[1])
+
+    return y, h
+
+
+def parity_accumulate_float32(matrix, data):
+    """
+    Multiply matrix @ data with MATLAB-compatible float32 accumulation.
+
+    MATLAB accumulates this product column-major; replicate that float32
+    accumulation order in NumPy (row-major) by transposing the operands.
+    Algebraically, (data.T @ matrix.T).T == matrix @ data.
+
+    Args:
+        matrix: 2D numpy array (e.g., refmatrix)
+        data: 2D numpy array (channels x points)
+
+    Returns:
+        The matrix product with MATLAB-parity accumulation order.
+    """
+    dt = matrix.dtype
+    block = np.ascontiguousarray(data.astype(dt).T)
+    return (block @ matrix.T).T

src/eegprep/functions/popfunc/eeg_eegrej.py

@@ -6,7 +6,7 @@
 from copy import deepcopy
 from eegprep.functions.miscfunc.event_utils import boundary_event_indices
 from eegprep.functions.miscfunc.event_utils import is_boundary_event as _is_boundary_event
-from ..miscfunc.misc import round_mat
+from ..miscfunc.parity import round_mat
 
 
 logger = logging.getLogger(__name__)