conditioning_numpy ¤

def condition_sparameters(
    S: np.ndarray,
    freqs: np.ndarray,
    *,
    passivity: bool = True,
    reciprocity: bool = True,
    causality: bool = True,
    preserve_dc: bool = False,
    passivity_limit: float = 1.0,
    tolerance: float = 1e-8,
    max_iterations: int = 100,
) -> tuple[np.ndarray, ConditioningReport]:
    """Alternate sample-domain projections, checking every final constraint.

    Causality requires an explicitly supplied DC-to-fmax uniform grid. No
    extrapolation or resampling happens implicitly. The finite FFT assumes a
    periodic impulse window; its causality check is a grid-dependent proxy.
    DC preservation retains the original real DC matrix and rejects infeasible
    endpoints. Disabled constraints are reported but do not gate convergence.
    An iteration limit returns ``converged=False``, never a claimed success.
    """
    original = _matrices(S)
    _check_grid(original, freqs, causality)
    if not np.isfinite(tolerance) or tolerance <= 0 or max_iterations < 1:
        raise ValueError("positive tolerance and iteration limit are required")
    if not np.isfinite(passivity_limit) or not 0 < passivity_limit <= 1:
        raise ValueError("passivity limit must be in (0, 1]")
    dc = None
    if preserve_dc:
        if not causality:
            raise ValueError("DC preservation requires the causal projection")
        dc = original[..., 0, :, :]
        if np.max(np.abs(dc.imag)) > tolerance:
            raise ValueError("preserved DC matrix must be real")
        dc = dc.real
        if passivity and np.max(np.linalg.svd(dc, compute_uv=False)) > passivity_limit + tolerance:
            raise ValueError("preserved DC matrix violates the passivity limit")
        if reciprocity and np.max(np.abs(dc - dc.swapaxes(-1, -2))) > tolerance:
            raise ValueError("preserved DC matrix violates reciprocity")
    result = original.copy()
    converged = False
    for iteration in range(1, max_iterations + 1):
        if reciprocity:
            result = project_s_reciprocal(result)
        if passivity:
            result = project_s_passive(result, passivity_limit)
        if causality:
            result = _causal_projection(result, dc)
        sigma = float(np.max(np.linalg.svd(result, compute_uv=False)))
        symmetry = float(np.max(np.abs(result - result.swapaxes(-1, -2))))
        negative = None
        if causality:
            impulse = _impulse(result)
            negative = float(np.linalg.norm(impulse[..., result.shape[-3] :, :, :]) / max(np.linalg.norm(impulse), 1e-30))
        dc_error = None if dc is None else float(np.max(np.abs(result[..., 0, :, :] - dc)))
        converged = (
            (not passivity or sigma <= passivity_limit + tolerance)
            and (not reciprocity or symmetry <= tolerance)
            and (not causality or negative <= tolerance)
            and (dc_error is None or dc_error <= tolerance)
        )
        if converged:
            break
    correction = result - original
    report = ConditioningReport(
        converged,
        iteration,
        sigma,
        symmetry,
        negative,
        dc_error,
        float(np.linalg.norm(correction) / max(np.linalg.norm(original), 1e-30)),
        float(np.max(np.abs(correction))),
    )
    return result, report

def project_s_passive(S: np.ndarray, limit: float = 1.0) -> np.ndarray:
    """Nearest matrix in Frobenius norm with singular values at most limit.

    Applies to power-normalized S data (e.g. real positive reference
    impedances). All frequency and batch axes share one stacked SVD call.
    """
    S = _matrices(S)
    if not np.isfinite(limit) or not 0 < limit <= 1:
        raise ValueError("passivity limit must be in (0, 1]")
    u, singular, vh = np.linalg.svd(S, full_matrices=False)
    return (u * np.minimum(singular, limit)[..., None, :]) @ vh

def project_s_reciprocal(S: np.ndarray) -> np.ndarray:
    """Project onto transpose symmetry without complex conjugation."""
    S = _matrices(S)
    return 0.5 * (S + S.swapaxes(-1, -2))