fitting ¤

def circuit_from_coefficients(
    coefficients: ModelCoefficients | str | Path, *, name: str = "FittedModel", holomorphic: bool = True
) -> Circuit:
    """Build a simulation-ready :class:`circulax.Circuit` from coefficients.

    The return type is the same for rational cores, delayed models, and ideal
    through lines. Fitting and passivity enforcement never occur here.
    """
    from circulax import Circuit, compile_circuit

    realized = _component_from_coefficients(coefficients, name=name, holomorphic=holomorphic)
    if isinstance(realized, Circuit):
        return realized
    ports = {port: f"model,{port}" for port in realized.ports}
    net = {"instances": {"model": {"component": "model"}}, "connections": {}, "ports": ports}
    return compile_circuit(net, {"model": realized}, g_leak=0)

def fit_model(
    S: np.ndarray,
    freqs: np.ndarray,
    *,
    options: ModelFitOptions | None = None,
    z0: float = 50.0,
    initial_poles: np.ndarray | None = None,
) -> ModelCoefficients:
    """Fit S data and return coefficients only, without making a component.

    Automatic order selection is the default. Optional enforcement is sampled,
    not globally certified; independent rational testing remains necessary.
    """
    options = options or ModelFitOptions()
    S, freqs = np.asarray(S, complex), np.asarray(freqs, float)
    if freqs.ndim != 1 or not len(freqs) or not np.all(np.isfinite(freqs)) or freqs[0] < 0 or np.any(np.diff(freqs) <= 0):
        raise ValueError("frequencies must be finite, nonnegative and increasing")
    if S.ndim != 3 or S.shape[0] != len(freqs) or S.shape[1] != S.shape[2] or not S.shape[1] or not np.all(np.isfinite(S)):
        raise ValueError("S must be finite with shape (frequencies, ports, ports)")
    # Validate impedance before doing any potentially expensive fitting.
    ModelCoefficients(np.array([], complex), np.zeros((*S.shape[1:], 0)), np.zeros(S.shape[1:]), z0)
    if initial_poles is not None and options.vector_fit_order is not None:
        raise ValueError("Choose initial_poles or vector_fit_order, not both")
    if options.delay_mode == "infer":
        if initial_poles is not None:
            raise ValueError("delay inference cannot be combined with supplied poles")
        from .delay_selection import fit_auto_delay

        return fit_auto_delay(S, freqs, z0, options)
    if options.reciprocal and not np.allclose(S, S.swapaxes(1, 2), atol=1e-8, rtol=0):
        raise ValueError("reciprocal fitting requires symmetric data")
    fit_started = time.perf_counter()
    original = S
    delays = np.zeros(S.shape[1]) if options.port_delays is None else np.asarray(options.port_delays)
    if delays.shape != (S.shape[1],):
        raise ValueError("port_delays must match the number of ports")
    phase = np.exp(2j * np.pi * freqs[:, None] * delays)
    S = S * phase[:, :, None] * phase[:, None, :]
    static = VFModel(np.empty(0, complex), np.empty((*S.shape[1:], 0), complex), S.real.mean(axis=0), np.zeros(S.shape[1:]))
    if S.shape[1] == 2 and np.max(np.abs(S - np.array([[0, 1], [1, 0]]))) <= min(
        1e-14, options.normalized_rmse / 2, options.max_absolute_error / 2
    ):
        static = VFModel(static.poles, static.residues, np.array([[0.0, 1.0], [1.0, 0.0]]), static.E)
    static_error, static_max = errors_numpy(static, S, freqs)
    if (
        options.delay_mode == "supplied"
        and initial_poles is None
        and static_error <= min(options.normalized_rmse, 1e-12)
        and static_max <= min(options.max_absolute_error, 1e-12)
    ):
        model, diagnostics = static, {"fitter": "static", "pole_count": 0}
    elif initial_poles is None and options.method == "vector_fitting":
        model, diagnostics = _fit_vector_fitting(S, freqs, z0, options)
    elif initial_poles is None:
        model, diagnostics = fit_s_numpy(
            S,
            freqs,
            tol=options.tol,
            mmax=options.mmax,
            iterations=options.iterations,
            reciprocal=options.reciprocal,
            screening=options.screening,
            reduction_stage=options.reduction_stage,
            aaa_backend=options.aaa_backend,
            normalized_rmse=options.normalized_rmse,
            max_absolute_error=options.max_absolute_error,
        )
    else:
        poles = np.asarray(initial_poles, complex)
        if poles.ndim != 1 or not np.all(np.isfinite(poles)) or np.any(poles.real >= 0):
            raise ValueError("initial poles must be finite and strictly stable")
        pole_groups(poles)
        model = refine_numpy(S, freqs, poles, iterations=options.iterations, reciprocal=options.reciprocal)
        diagnostics = {
            "fitter": "numpy-supplied-poles",
            "order_selection": "supplied",
            "initial_poles": [[float(p.real), float(p.imag)] for p in poles],
        }
    if options.delay_mode == "supplied" and initial_poles is None and options.vector_fit_order is None and len(model.poles):
        diagnostics["unreduced_core_poles"] = len(model.poles)
        for subset in candidate_subsets(model, freqs):
            candidate = refit_numpy(S, 2j * np.pi * freqs, model.poles[subset])
            candidate_error, candidate_max = errors_numpy(candidate, S, freqs)
            if candidate_error <= options.normalized_rmse and candidate_max <= options.max_absolute_error:
                model = candidate
                break
    fitting_seconds = time.perf_counter() - fit_started
    enforcement_started = time.perf_counter()
    if options.enforce_passivity:
        model, report = enforce_s_passivity_numpy(
            model,
            freqs,
            enforcement_freqs=options.enforcement_freqs,
            limit=options.passivity_limit,
            max_iterations=options.enforcement_iterations,
        )
        diagnostics["enforcement"] = report
        if not report["converged"]:
            raise ValueError(f"Passivity enforcement failed: {report['optimizer_message']}")
    enforcement_seconds = time.perf_counter() - enforcement_started if options.enforce_passivity else 0.0
    error, maximum = errors_numpy(model, S, freqs)
    if error > options.normalized_rmse or maximum > options.max_absolute_error:
        raise ValueError(f"Final fit exceeds accuracy limits: NRMSE={error:.6g}, max error={maximum:.6g}")
    diagnostics.update(
        training_nrmse=error,
        training_max_error=maximum,
        options=json.loads(json.dumps(asdict(options))),
        global_passivity_certified=False,
    )
    result = ModelCoefficients(model.poles, model.residues, model.D, z0, (float(freqs[0]), float(freqs[-1])), diagnostics, delays)
    difference = result.evaluate(freqs) - original
    error = float(np.linalg.norm(difference) / max(np.linalg.norm(original), 1e-30))
    maximum = float(np.max(np.abs(difference)))
    if error > options.normalized_rmse or maximum > options.max_absolute_error:
        raise ValueError("Final fit exceeds original-domain accuracy limits")
    diagnostics.update(
        training_nrmse=error,
        training_max_error=maximum,
        core_pole_count=len(model.poles),
        delay_provenance=options.delay_mode,
        port_delays_seconds=delays.tolist(),
    )
    diagnostics.update(fitting_seconds=fitting_seconds, enforcement_seconds=enforcement_seconds, pole_count=len(model.poles))
    return result

def validate_model(
    coefficients: ModelCoefficients | str | Path,
    *,
    measured_S: np.ndarray | None = None,
    freqs: np.ndarray | None = None,
    validation_S: np.ndarray | None = None,
    validation_freqs: np.ndarray | None = None,
    simulation_frequency_range: tuple[float, float] | None = None,
    expected_passive: bool = True,
    expected_reciprocal: bool = True,
    normalized_rmse: float = 0.02,
    max_absolute_error: float = 0.05,
    passivity_tolerance: float = 1e-9,
) -> ModelValidationReport:
    """Validate fitted coefficients without compiling a circuit.

    Accuracy checks are performed when measured data are supplied. Independent
    validation requires both ``validation_S`` and ``validation_freqs``.
    Passivity is required by default; active models must opt out explicitly.
    """
    source = ModelCoefficients.load(coefficients) if isinstance(coefficients, (str, Path)) else coefficients
    source = ModelCoefficients(
        source.poles, source.residues, source.D, source.z0, source.frequency_range, source.metadata, source.port_delays
    )
    if (measured_S is None) != (freqs is None):
        raise ValueError("measured_S and freqs must be supplied together")
    if (validation_S is None) != (validation_freqs is None):
        raise ValueError("validation_S and validation_freqs must be supplied together")
    limits = (normalized_rmse, max_absolute_error, passivity_tolerance)
    if not np.all(np.isfinite(limits)) or normalized_rmse <= 0 or max_absolute_error <= 0 or passivity_tolerance < 0:
        raise ValueError("validation limits must be finite; error limits must be positive and passivity tolerance nonnegative")

    ports = len(source.D)
    training_data = None
    validation_data = None
    if measured_S is not None:
        training_data = _validation_dataset(measured_S, freqs, ports, "training")
    if validation_S is not None:
        validation_data = _validation_dataset(validation_S, validation_freqs, ports, "held-out")

    findings: list[str] = []
    physical_grids = [data[1] for data in (training_data, validation_data) if data is not None]
    evaluation_freqs = np.unique(np.concatenate(physical_grids)) if physical_grids else np.asarray(source.frequency_range, float)
    response = source.evaluate(evaluation_freqs)
    training_nrmse = training_max = validation_nrmse = validation_max = None
    if training_data is not None:
        training_S, training_freqs = training_data
        training_nrmse, training_max = _validation_errors(source.evaluate(training_freqs), training_S)
        if training_nrmse > normalized_rmse or training_max > max_absolute_error:
            findings.append("training accuracy exceeds the requested limits")
    else:
        findings.append("no measured training data supplied")
    if validation_data is not None:
        heldout_S, heldout_freqs = validation_data
        validation_nrmse, validation_max = _validation_errors(source.evaluate(heldout_freqs), heldout_S)
        if validation_nrmse > normalized_rmse or validation_max > max_absolute_error:
            findings.append("held-out accuracy exceeds the requested limits")
    else:
        findings.append("no independent validation data supplied")

    if expected_reciprocal and np.max(np.abs(response - response.swapaxes(1, 2))) > 1e-8:
        findings.append("model is not reciprocal")
    maximum_singular = float(np.max(np.linalg.svd(response, compute_uv=False)))
    if expected_passive and maximum_singular > 1 + passivity_tolerance:
        findings.append("model fails sampled passivity")
    maximum_s_pole = float(np.max(source.poles.real)) if len(source.poles) else -np.inf
    if maximum_s_pole >= 0:
        findings.append("model has unstable S poles")

    n = len(source.D)
    through = not len(source.poles) and n == 2 and np.array_equal(source.D, [[0, 1], [1, 0]])
    maximum_y_pole = -np.inf
    if not through:
        if np.linalg.matrix_rank(np.eye(n) + source.D) < n:
            findings.append("model has an unsupported singular I+D core")
        elif len(source.poles):
            model = VFModel(source.poles, source.residues, source.D.real, np.zeros((n, n)))
            try:
                ss, _ = scattering_state_space_to_admittance(vfmodel_to_ss(model, n), z0=source.z0)
                maximum_y_pole = float(np.max(np.asarray(ss.A).real)) if len(ss.A) else -np.inf
                if maximum_y_pole >= 0:
                    findings.append("model has unstable Y poles")
            except np.linalg.LinAlgError:
                findings.append("S-to-Y realization failed")

    if simulation_frequency_range is not None:
        if len(simulation_frequency_range) != 2:
            raise ValueError("simulation_frequency_range must contain (minimum, maximum) Hz")
        low, high = simulation_frequency_range
        if not np.all(np.isfinite((low, high))) or not 0 <= low <= high:
            raise ValueError("simulation_frequency_range must be finite, nonnegative, and ordered")
        if low < source.frequency_range[0] or high > source.frequency_range[1]:
            findings.append("requested simulation range exceeds the fitted frequency range")

    failures = tuple(item for item in findings if not item.startswith("no measured") and not item.startswith("no independent"))
    status: Literal["pass", "warn", "fail"] = "fail" if failures else ("warn" if findings else "pass")
    return ModelValidationReport(
        status,
        training_nrmse,
        training_max,
        validation_nrmse,
        validation_max,
        maximum_singular,
        maximum_s_pole,
        maximum_y_pole,
        tuple(findings),
    )