When a Fit Fails¶

def synthesize_data() -> tuple[np.ndarray, np.ndarray]:
    """One clean, well-sampled Gaussian peak with light noise.

    Deliberately easy data: nothing below fails because the measurement is
    hard. Every failure on this page is caused by where the fit was told to
    start, not by the data.
    """
    rng = np.random.default_rng(7)
    x = np.linspace(-4, 4, 120)
    y_true = AMPLITUDE_TRUE * np.exp(-0.5 * ((x - CENTER_TRUE) / SIGMA_TRUE) ** 2)
    return x, y_true + rng.normal(0, 0.05, len(x))


x, y = synthesize_data()

def build_graph(amplitude: float, center: float, sigma: float) -> FitGraph:
    """One Gaussian, with the *initial guess* supplied by the caller."""
    return FitGraph(
        nodes=[
            ModelNodeSpec(
                id="peak",
                model_type=ModelType.GAUSSIAN,
                parameters={
                    "amplitude": Parameter(value=amplitude),
                    "center": Parameter(value=center),
                    "sigma": Parameter(value=sigma, min=1e-3),
                },
            ),
        ],
    )

def three_fits(x: np.ndarray, y: np.ndarray) -> dict[str, FitResult]:
    """Same data, same model, three outcomes.

    ``reported`` and ``silent`` differ only in their starting guess; the
    solver is local ``"lm"`` in both. ``recovered`` reuses ``silent``'s bad
    guess and changes only the solver.
    """
    data = MeasurementData(x=x.tolist(), y=y.tolist())
    return {
        "reported": fit(build_graph(50.0, -3.9, 0.05), data, FitOptions(solver="lm")),
        "silent": fit(build_graph(3.0, -3.9, 0.8), data, FitOptions(solver="lm")),
        "recovered": fit(build_graph(3.0, -3.9, 0.8), data, FitOptions(solver="global")),
    }


results = three_fits(x, y)

def inspect(results: dict[str, FitResult], y: np.ndarray) -> None:
    """Print what each fit reported, and check the three outcomes hold.

    The ``|A|/max|y|`` column is the quantity the degeneracy guard actually
    tests. Watch it straddle the 1% threshold between the first two rows
    while ``r^2`` stays essentially identical -- that is the guard's
    boundary, made visible.
    """
    y_max_abs = float(np.abs(y).max())
    header = (
        f"{'case':11s} {'success':>8s} {'r^2':>9s} {'amplitude':>10s} {'|A|/max|y|':>11s}  message"
    )
    print(header)
    for name, result in results.items():
        amp = result.parameters["peak.amplitude"].value
        print(
            f"{name:11s} {result.success!s:>8s} {result.r_squared:9.4f} "
            f"{amp:10.4f} {abs(amp) / y_max_abs:11.4f}  {result.message}",
        )

    reported, silent, recovered = results["reported"], results["silent"], results["recovered"]

    # The guard fires: amplitude collapsed under 1% of max|y| with r^2 < 0.
    assert not reported.success
    assert "degenerate_fit" in reported.message
    assert reported.r_squared < 0

    # The guard stays silent on a fit that is just as wrong, because its
    # amplitude landed above the 1% threshold. This is the case the page exists
    # for: `success` alone would have told a caller everything was fine.
    assert silent.success
    assert silent.r_squared < 0
    assert abs(silent.parameters["peak.amplitude"].value) / y_max_abs > 1e-2

    # Same bad guess, global solver, truth recovered.
    assert recovered.success
    assert recovered.r_squared > 0.99
    assert abs(recovered.parameters["peak.amplitude"].value - AMPLITUDE_TRUE) < 0.05
    assert abs(recovered.parameters["peak.center"].value - CENTER_TRUE) < 0.05


inspect(results, y)