Plotting and diagnostic reports

Plotting functions import Matplotlib lazily and return (figure, axes). Install the plot extra before using them. Diagnostic I/O accepts either a recorder dictionary or a JSON path where documented.

plot

Matplotlib rendering helpers for reconstruction results and diagnostics.

plot_reconstruction(truth: np.ndarray, result: ReconstructionResult, *, layout: Sequence[Sequence[str]] | str | None = None, figsize: tuple[float, float] = (13.0, 7.0)) -> PlotResult

Plot ground truth, reconstruction, and loss history.

Source code in target/docs-python-api/fpm_rs/plot.py

def plot_reconstruction(
    truth: np.ndarray,
    result: ReconstructionResult,
    *,
    layout: Sequence[Sequence[str]] | str | None = None,
    figsize: tuple[float, float] = (13.0, 7.0),
) -> PlotResult:
    """Plot ground truth, reconstruction, and loss history."""
    plt = _import_pyplot()

    truth_array = np.asarray(truth)
    amplitude = np.asarray(result.amplitude)
    phase = np.asarray(result.phase)
    if truth_array.ndim != 2:
        raise ValueError(
            f"truth must be two-dimensional, got shape {truth_array.shape}"
        )
    if truth_array.size == 0:
        raise ValueError("truth must not be empty")
    if amplitude.shape != truth_array.shape or phase.shape != truth_array.shape:
        raise ValueError(
            "truth, reconstructed amplitude, and reconstructed phase must have "
            f"the same shape; got {truth_array.shape}, {amplitude.shape}, and "
            f"{phase.shape}"
        )

    figure, axes = plt.subplot_mosaic(
        _RECONSTRUCTION_LAYOUT if layout is None else layout,
        figsize=figsize,
        constrained_layout=True,
    )
    missing = _RECONSTRUCTION_REQUIRED_AXES.difference(axes)
    if missing:
        plt.close(figure)
        raise ValueError(
            f"layout is missing required axes: {', '.join(sorted(missing))}"
        )

    amplitude_max = max(float(np.abs(truth_array).max()), float(amplitude.max()))
    axes["true_amplitude"].imshow(
        np.abs(truth_array), cmap="gray", vmin=0.0, vmax=amplitude_max
    )
    axes["true_amplitude"].set_title("Original amplitude")
    axes["reconstructed_amplitude"].imshow(
        amplitude, cmap="gray", vmin=0.0, vmax=amplitude_max
    )
    axes["reconstructed_amplitude"].set_title("Reconstructed amplitude")

    axes["true_phase"].imshow(
        np.angle(truth_array), cmap="twilight", vmin=-np.pi, vmax=np.pi
    )
    axes["true_phase"].set_title("Original phase")
    axes["reconstructed_phase"].imshow(phase, cmap="twilight", vmin=-np.pi, vmax=np.pi)
    axes["reconstructed_phase"].set_title("Reconstructed phase")

    iterations = [record[0] for record in result.history]
    losses = [record[1] for record in result.history]
    axes["loss"].semilogy(iterations, losses, marker="o")
    axes["loss"].set(
        title="Reconstruction error (loss)", xlabel="Iteration", ylabel="Loss"
    )
    axes["loss"].grid(True, which="both", alpha=0.3)

    for name in _RECONSTRUCTION_IMAGE_AXES:
        axes[name].set_axis_off()

    return figure, axes

plot_convergence(diagnostics: DiagnosticsData, *, figsize: tuple[float, float] = (13.0, 9.0)) -> PlotResult | None

Plot convergence diagnostics and return the figure and named axes.

Source code in target/docs-python-api/fpm_rs/plot.py

def plot_convergence(
    diagnostics: DiagnosticsData,
    *,
    figsize: tuple[float, float] = (13.0, 9.0),
) -> PlotResult | None:
    """Plot convergence diagnostics and return the figure and named axes."""
    history = diagnostics.get("iteration_history")
    if not isinstance(history, list) or not history:
        return None

    iterations = _series(history, "iteration")
    if np.isnan(iterations).all():
        iterations = np.arange(len(history), dtype=float)
    total_loss = _series(history, "total_loss")
    data_loss = _series(history, "data_loss")
    regularization_loss = _series(history, "regularization_loss")
    object_change = _series(history, "object_relative_change")
    pupil_change = _series(history, "pupil_relative_change")
    median_frame_loss = _series(history, "median_frame_loss")
    worst_frame_loss = _series(history, "worst_frame_loss")
    elapsed_ms = _series(history, "elapsed_ms")

    plt = _import_pyplot()
    figure, raw_axes = plt.subplots(2, 2, figsize=figsize, constrained_layout=True)
    raw_axes = raw_axes.ravel()
    axes = {
        "loss": raw_axes[0],
        "relative_change": raw_axes[1],
        "frame_loss": raw_axes[2],
        "elapsed_time": raw_axes[3],
    }

    _plot_lines(
        axes["loss"],
        iterations,
        [
            ("total loss", total_loss),
            ("data loss", data_loss),
            ("regularization loss", regularization_loss),
        ],
        title="Loss history",
        xlabel="Iteration",
        ylabel="Loss",
        log_scale=True,
    )
    _plot_lines(
        axes["relative_change"],
        iterations,
        [("object change", object_change), ("pupil change", pupil_change)],
        title="Relative change",
        xlabel="Iteration",
        ylabel="Relative change",
        log_scale=True,
    )
    _plot_lines(
        axes["frame_loss"],
        iterations,
        [
            ("median frame loss", median_frame_loss),
            ("worst frame loss", worst_frame_loss),
        ],
        title="Per-frame loss",
        xlabel="Iteration",
        ylabel="Loss",
        log_scale=True,
    )
    _plot_lines(
        axes["elapsed_time"],
        iterations,
        [("elapsed ms", elapsed_ms)],
        title="Elapsed time",
        xlabel="Iteration",
        ylabel="Milliseconds",
        log_scale=False,
    )
    return figure, axes

plot_frame_residuals(diagnostics: DiagnosticsData, *, figsize: tuple[float, float] = (13.0, 9.0)) -> PlotResult | None

Plot latest per-frame residual summaries and return named axes.

Source code in target/docs-python-api/fpm_rs/plot.py

def plot_frame_residuals(
    diagnostics: DiagnosticsData,
    *,
    figsize: tuple[float, float] = (13.0, 9.0),
) -> PlotResult | None:
    """Plot latest per-frame residual summaries and return named axes."""
    frame_diagnostics = latest_frame_diagnostics(diagnostics.get("frame_diagnostics"))
    if not frame_diagnostics:
        return None

    illumination_index = _series(frame_diagnostics, "illumination_index")
    if np.isnan(illumination_index).all():
        illumination_index = _series(frame_diagnostics, "frame_index")
    residual_l2 = _series(frame_diagnostics, "residual_l2")
    normalized_l2 = _series(frame_diagnostics, "normalized_l2")
    residual_mean = _series(frame_diagnostics, "residual_mean")
    residual_std = _series(frame_diagnostics, "residual_std")
    residual_max_abs = _series(frame_diagnostics, "residual_max_abs")
    measured_sum = _series(frame_diagnostics, "measured_sum")
    predicted_sum = _series(frame_diagnostics, "predicted_sum")

    plt = _import_pyplot()
    figure, raw_axes = plt.subplots(2, 2, figsize=figsize, constrained_layout=True)
    raw_axes = raw_axes.ravel()
    axes = {
        "residual_magnitude": raw_axes[0],
        "residual_moments": raw_axes[1],
        "frame_sums": raw_axes[2],
        "measured_vs_predicted": raw_axes[3],
    }

    _plot_lines(
        axes["residual_magnitude"],
        illumination_index,
        [("residual L2", residual_l2), ("normalized L2", normalized_l2)],
        title="Residual magnitude",
        xlabel="Illumination index",
        ylabel="Residual",
        log_scale=True,
    )
    _plot_lines(
        axes["residual_moments"],
        illumination_index,
        [
            ("residual mean", residual_mean),
            ("residual std", residual_std),
            ("residual max abs", residual_max_abs),
        ],
        title="Residual moments",
        xlabel="Illumination index",
        ylabel="Value",
        log_scale=False,
    )
    _plot_lines(
        axes["frame_sums"],
        illumination_index,
        [("measured sum", measured_sum), ("predicted sum", predicted_sum)],
        title="Frame sums",
        xlabel="Illumination index",
        ylabel="Sum",
        log_scale=False,
    )
    finite_mask = np.isfinite(measured_sum) & np.isfinite(predicted_sum)
    if np.any(finite_mask):
        scatter = axes["measured_vs_predicted"].scatter(
            measured_sum[finite_mask],
            predicted_sum[finite_mask],
            c=normalized_l2[finite_mask],
            cmap="viridis",
            edgecolor="none",
        )
        minimum = float(
            np.nanmin(
                [measured_sum[finite_mask].min(), predicted_sum[finite_mask].min()]
            )
        )
        maximum = float(
            np.nanmax(
                [measured_sum[finite_mask].max(), predicted_sum[finite_mask].max()]
            )
        )
        axes["measured_vs_predicted"].plot(
            [minimum, maximum],
            [minimum, maximum],
            color="0.4",
            linestyle="--",
            linewidth=1.0,
        )
        figure.colorbar(
            scatter,
            ax=axes["measured_vs_predicted"],
            label="Normalized L2",
        )
    axes["measured_vs_predicted"].set(
        title="Measured vs predicted sums",
        xlabel="Measured sum",
        ylabel="Predicted sum",
    )
    axes["measured_vs_predicted"].grid(True, alpha=0.25)
    return figure, axes

plot_raw_stack_stats(diagnostics: DiagnosticsData, *, figsize: tuple[float, float] = (13.0, 9.0)) -> PlotResult | None

Plot raw measurement-stack statistics and return named axes.

Source code in target/docs-python-api/fpm_rs/plot.py

def plot_raw_stack_stats(
    diagnostics: DiagnosticsData,
    *,
    figsize: tuple[float, float] = (13.0, 9.0),
) -> PlotResult | None:
    """Plot raw measurement-stack statistics and return named axes."""
    raw_stats = diagnostics.get("raw_frame_stats")
    if not isinstance(raw_stats, list) or not raw_stats:
        return None

    frame_index = _series(raw_stats, "frame_index")
    if np.isnan(frame_index).all():
        frame_index = np.arange(len(raw_stats), dtype=float)
    mean = _series(raw_stats, "mean")
    std = _series(raw_stats, "std")
    minimum = _series(raw_stats, "min")
    maximum = _series(raw_stats, "max")
    saturated = _series(raw_stats, "saturated_pixels")
    zero_pixels = _series(raw_stats, "zero_pixels")

    plt = _import_pyplot()
    figure, raw_axes = plt.subplots(2, 2, figsize=figsize, constrained_layout=True)
    raw_axes = raw_axes.ravel()
    axes = {
        "mean_std": raw_axes[0],
        "min_max": raw_axes[1],
        "saturated_pixels": raw_axes[2],
        "zero_pixels": raw_axes[3],
    }

    _plot_lines(
        axes["mean_std"],
        frame_index,
        [("mean", mean), ("std", std)],
        title="Raw-frame mean and std",
        xlabel="Frame index",
        ylabel="Value",
        log_scale=False,
    )
    _plot_lines(
        axes["min_max"],
        frame_index,
        [("min", minimum), ("max", maximum)],
        title="Raw-frame min and max",
        xlabel="Frame index",
        ylabel="Value",
        log_scale=False,
    )
    _plot_lines(
        axes["saturated_pixels"],
        frame_index,
        [("saturated pixels", saturated)],
        title="Saturated pixels",
        xlabel="Frame index",
        ylabel="Pixels",
        log_scale=False,
    )
    _plot_lines(
        axes["zero_pixels"],
        frame_index,
        [("zero pixels", zero_pixels)],
        title="Zero pixels",
        xlabel="Frame index",
        ylabel="Pixels",
        log_scale=False,
    )
    return figure, axes

plot_fourier_coverage(diagnostics: DiagnosticsData, *, figsize: tuple[float, float] = (8.5, 8.5)) -> PlotResult | None

Plot shifted pupils in Fourier space and return the figure and axis.

Source code in target/docs-python-api/fpm_rs/plot.py

def plot_fourier_coverage(
    diagnostics: DiagnosticsData,
    *,
    figsize: tuple[float, float] = (8.5, 8.5),
) -> PlotResult | None:
    """Plot shifted pupils in Fourier space and return the figure and axis."""
    coverage = diagnostics.get("coverage")
    if not isinstance(coverage, dict):
        return None
    centers = _array2d(coverage.get("pupil_centers_px"))
    radius = _float_or_nan(coverage.get("pupil_radius_px"))
    if centers.size == 0 or not np.isfinite(radius) or radius <= 0.0:
        return None

    illumination_na = _series_from_sequence(coverage.get("illumination_na"))
    synthetic_na = _float_or_nan(coverage.get("synthetic_na"))

    plt = _import_pyplot()
    from matplotlib.patches import Circle

    figure, ax = plt.subplots(figsize=figsize, constrained_layout=True)

    centers_x = centers[:, 0]
    centers_y = centers[:, 1]
    if illumination_na.size == centers.shape[0] and np.isfinite(illumination_na).any():
        scatter = ax.scatter(
            centers_x,
            centers_y,
            c=illumination_na,
            cmap="viridis",
            s=30,
            zorder=3,
        )
        figure.colorbar(scatter, ax=ax, label="Illumination NA")
    else:
        ax.scatter(centers_x, centers_y, color="C0", s=30, zorder=3)

    for center_x, center_y in centers:
        ax.add_patch(
            Circle((center_x, center_y), radius, fill=False, edgecolor="C1", alpha=0.35)
        )

    if np.isfinite(synthetic_na):
        ax.set_title(f"Fourier coverage (synthetic NA: {synthetic_na:.4g})")
    else:
        ax.set_title("Fourier coverage")
    ax.set_xlabel("Fourier x (px)")
    ax.set_ylabel("Fourier y (px)")
    ax.set_aspect("equal", adjustable="datalim")
    ax.grid(True, alpha=0.2)
    _set_limits_from_centers(ax, centers, radius)
    return figure, {"coverage": ax}

plot_crop_indices(diagnostics: DiagnosticsData, *, figsize: tuple[float, float] = (8.5, 8.5)) -> PlotResult | None

Plot Fourier crop boxes and return the figure and axis.

Source code in target/docs-python-api/fpm_rs/plot.py

def plot_crop_indices(
    diagnostics: DiagnosticsData,
    *,
    figsize: tuple[float, float] = (8.5, 8.5),
) -> PlotResult | None:
    """Plot Fourier crop boxes and return the figure and axis."""
    coverage = diagnostics.get("coverage")
    if not isinstance(coverage, dict):
        return None
    crop_indices = coverage.get("crop_indices")
    if not isinstance(crop_indices, list) or not crop_indices:
        return None

    plt = _import_pyplot()
    from matplotlib.patches import Rectangle

    figure, ax = plt.subplots(figsize=figsize, constrained_layout=True)
    centers_x: list[float] = []
    centers_y: list[float] = []
    x_bounds: list[float] = []
    y_bounds: list[float] = []
    for crop in crop_indices:
        if not isinstance(crop, dict):
            continue
        x_start = _float_or_nan(crop.get("x_start"))
        y_start = _float_or_nan(crop.get("y_start"))
        x_end = _float_or_nan(crop.get("x_end"))
        y_end = _float_or_nan(crop.get("y_end"))
        center_x = _float_or_nan(crop.get("center_x"))
        center_y = _float_or_nan(crop.get("center_y"))
        if not np.isfinite([x_start, y_start, x_end, y_end, center_x, center_y]).all():
            continue
        ax.add_patch(
            Rectangle(
                (x_start, y_start),
                x_end - x_start,
                y_end - y_start,
                fill=False,
                edgecolor="C0",
                alpha=0.5,
            )
        )
        centers_x.append(center_x)
        centers_y.append(center_y)
        x_bounds.extend([x_start, x_end])
        y_bounds.extend([y_start, y_end])
    if not centers_x:
        plt.close(figure)
        return None
    ax.scatter(centers_x, centers_y, color="C1", s=25, zorder=3)
    ax.set_title("Crop indices")
    ax.set_xlabel("Fourier x (px)")
    ax.set_ylabel("Fourier y (px)")
    ax.set_aspect("equal", adjustable="datalim")
    ax.grid(True, alpha=0.2)
    _set_limits_from_bounds(ax, np.asarray(x_bounds), np.asarray(y_bounds))
    return figure, {"crop_indices": ax}

plot_residuals_on_fourier_centers(diagnostics: DiagnosticsData, *, figsize: tuple[float, float] = (8.5, 8.5)) -> PlotResult | None

Map latest frame residuals onto Fourier centers and return the axis.

Source code in target/docs-python-api/fpm_rs/plot.py

def plot_residuals_on_fourier_centers(
    diagnostics: DiagnosticsData,
    *,
    figsize: tuple[float, float] = (8.5, 8.5),
) -> PlotResult | None:
    """Map latest frame residuals onto Fourier centers and return the axis."""
    coverage = diagnostics.get("coverage")
    frame_diagnostics = latest_frame_diagnostics(diagnostics.get("frame_diagnostics"))
    if not isinstance(coverage, dict) or not frame_diagnostics:
        return None
    centers = _array2d(coverage.get("pupil_centers_px"))
    if centers.size == 0:
        return None
    residual_by_illumination: dict[int, float] = {}
    for entry in frame_diagnostics:
        if not isinstance(entry, dict):
            continue
        illumination_index = entry.get("illumination_index")
        residual = _float_or_nan(entry.get("normalized_l2"))
        if isinstance(illumination_index, int) and np.isfinite(residual):
            residual_by_illumination[illumination_index] = residual
    if not residual_by_illumination:
        return None

    xs = []
    ys = []
    colors = []
    for illumination_index, residual in residual_by_illumination.items():
        if 0 <= illumination_index < centers.shape[0]:
            xs.append(float(centers[illumination_index, 0]))
            ys.append(float(centers[illumination_index, 1]))
            colors.append(float(residual))
    if not xs:
        return None

    plt = _import_pyplot()
    figure, ax = plt.subplots(figsize=figsize, constrained_layout=True)
    scatter = ax.scatter(xs, ys, c=colors, cmap="magma", s=35)
    figure.colorbar(scatter, ax=ax, label="Normalized L2 residual")
    ax.set_title("Residuals on Fourier centers")
    ax.set_xlabel("Fourier x (px)")
    ax.set_ylabel("Fourier y (px)")
    ax.set_aspect("equal", adjustable="datalim")
    ax.grid(True, alpha=0.2)
    _set_limits_from_points(ax, np.asarray(xs), np.asarray(ys))
    return figure, {"residuals_on_fourier_centers": ax}

diagnostics

Diagnostics loading, reporting, and lightweight text summaries.

__all__ module-attribute

coerce_diagnostics(value: Mapping[str, Any] | str | Path) -> dict[str, Any]

Source code in target/docs-python-api/fpm_rs/diagnostics/io.py

def coerce_diagnostics(value: Mapping[str, Any] | str | Path) -> dict[str, Any]:
    if isinstance(value, (str, Path)):
        return load_diagnostics(value)
    if isinstance(value, Mapping):
        return dict(value)
    raise TypeError(
        "diagnostics must be a mapping or a path to a diagnostics JSON file"
    )

ensure_output_dir(path: str | Path) -> Path

Source code in target/docs-python-api/fpm_rs/diagnostics/io.py

def ensure_output_dir(path: str | Path) -> Path:
    output_dir = Path(path)
    output_dir.mkdir(parents=True, exist_ok=True)
    return output_dir

load_diagnostics(path: str | Path) -> dict[str, Any]

Source code in target/docs-python-api/fpm_rs/diagnostics/io.py

def load_diagnostics(path: str | Path) -> dict[str, Any]:
    diagnostics_path = Path(path)
    if not diagnostics_path.exists():
        raise FileNotFoundError(f"diagnostics file not found: {diagnostics_path}")
    try:
        with diagnostics_path.open("r", encoding="utf-8") as handle:
            payload = json.load(handle)
    except json.JSONDecodeError as error:
        raise ValueError(
            f"invalid diagnostics JSON in {diagnostics_path}: {error}"
        ) from error
    if not isinstance(payload, dict):
        raise ValueError(
            f"diagnostics JSON must contain an object at the top level, got {type(payload).__name__}"
        )
    return payload

make_diagnostic_report(diagnostics: Mapping[str, Any] | str | Path, output_dir: str | Path = 'diagnostic_report') -> None

Source code in target/docs-python-api/fpm_rs/diagnostics/report.py

def make_diagnostic_report(
    diagnostics: Mapping[str, Any] | str | Path,
    output_dir: str | Path = "diagnostic_report",
) -> None:
    diagnostics = coerce_diagnostics(diagnostics)
    output_path = ensure_output_dir(output_dir)
    for filename, plotter in (
        ("convergence.png", plot_convergence),
        ("frame_residuals.png", plot_frame_residuals),
        ("raw_stack_stats.png", plot_raw_stack_stats),
        ("fourier_coverage.png", plot_fourier_coverage),
        ("crop_indices.png", plot_crop_indices),
        ("residuals_on_fourier_centers.png", plot_residuals_on_fourier_centers),
    ):
        result = plotter(diagnostics)
        if result is not None:
            result[0].savefig(output_path / filename, dpi=180, bbox_inches="tight")
            _close_figure(result[0])
    write_ground_truth_metrics(diagnostics, output_dir)
    write_summary(diagnostics, output_dir)

savefig(path: str | Path, dpi: int = 180) -> None

Source code in target/docs-python-api/fpm_rs/diagnostics/io.py

def savefig(path: str | Path, dpi: int = 180) -> None:
    import matplotlib.pyplot as plt

    figure = plt.gcf()
    figure.savefig(Path(path), dpi=dpi, bbox_inches="tight")

write_ground_truth_metrics(diag: dict[str, Any], output_dir: str | Path) -> None

Source code in target/docs-python-api/fpm_rs/diagnostics/report.py

def write_ground_truth_metrics(diag: dict[str, Any], output_dir: str | Path) -> None:
    metrics = diag.get("ground_truth_metrics")
    if not isinstance(metrics, dict) or not metrics:
        return
    output_path = ensure_output_dir(output_dir) / "ground_truth_metrics.txt"
    lines = []
    for key in (
        "amplitude_rmse",
        "amplitude_nrmse",
        "complex_rmse",
        "complex_nrmse",
        "phase_rmse",
        "phase_mae",
        "fourier_nrmse",
    ):
        value = metrics.get(key)
        if value is None:
            continue
        lines.append(f"{key}: {value}")
    output_path.write_text("\n".join(lines) + ("\n" if lines else ""), encoding="utf-8")

write_summary(diag: dict[str, Any], output_dir: str | Path) -> None

Source code in target/docs-python-api/fpm_rs/diagnostics/report.py

def write_summary(diag: dict[str, Any], output_dir: str | Path) -> None:
    output_path = ensure_output_dir(output_dir) / "summary.txt"
    iteration_history = diag.get("iteration_history")
    all_frame_diagnostics = diag.get("frame_diagnostics")
    frame_diagnostics = latest_frame_diagnostics(all_frame_diagnostics)
    raw_frame_stats = diag.get("raw_frame_stats")
    coverage = diag.get("coverage")
    ground_truth_metrics = diag.get("ground_truth_metrics")

    lines = [
        f"recorded_iterations: {len(iteration_history) if isinstance(iteration_history, list) else 0}",
        f"frame_diagnostics: {len(all_frame_diagnostics) if isinstance(all_frame_diagnostics, list) else 0}",
        f"raw_frame_stats: {len(raw_frame_stats) if isinstance(raw_frame_stats, list) else 0}",
        f"coverage_available: {bool(isinstance(coverage, dict) and coverage)}",
        f"ground_truth_metrics_available: {bool(isinstance(ground_truth_metrics, dict) and ground_truth_metrics)}",
    ]

    if isinstance(iteration_history, list) and iteration_history:
        last = iteration_history[-1]
        if isinstance(last, dict):
            lines.append("last_iteration:")
            for key in (
                "iteration",
                "total_loss",
                "data_loss",
                "regularization_loss",
                "object_relative_change",
                "pupil_relative_change",
                "median_frame_loss",
                "worst_frame_loss",
                "elapsed_ms",
            ):
                value = last.get(key)
                if value is not None:
                    lines.append(f"  {key}: {value}")

    worst_frame = _worst_frame(frame_diagnostics)
    if worst_frame is not None:
        frame_index, illumination_index, value = worst_frame
        lines.append("worst_frame:")
        lines.append(f"  frame_index: {frame_index}")
        lines.append(f"  illumination_index: {illumination_index}")
        lines.append(f"  normalized_l2: {value}")

    output_path.write_text("\n".join(lines) + "\n", encoding="utf-8")