1.1 module-scoped API¶

def sparse_fpca_score_uncertainty(
    fit: SparseFPCAResult,
    trajectories: IrregularTrajectorySet,
    *,
    condition_limit: float | None = None,
    failure_action: str = "error",
) -> SparseFPCAScoreUncertaintyResult:
    """Compute conditional covariance of native sparse PACE scores.

    For retained score vector ``xi_i`` and native observation times ``T_i``,
    this evaluates

    ``Lambda - Lambda Phi_i.T Sigma_i^{-1} Phi_i Lambda``

    using the same full fitted covariance-plus-noise system as native PACE
    scoring. The calculation is conditional on the fitted population objects;
    population-estimation uncertainty is deliberately excluded.
    """

    sparse = _native_sparse_provenance(fit)
    _validate_alignment(fit, trajectories)
    if failure_action not in {"error", "retain_nan"}:
        raise ValueError("failure_action must be 'error' or 'retain_nan'")

    grid, fitted_covariance, eigenvalues, eigenfunctions, noise_variance = (
        _required_population_arrays(fit)
    )
    n_components = int(fit.n_components)
    retained_values = eigenvalues[:n_components]
    retained_functions = eigenfunctions[:n_components]

    score_ridge = float(sparse.get("score_ridge", 0.0))
    if not np.isfinite(score_ridge) or score_ridge < 0:
        raise ValueError("fit score_ridge must be finite and non-negative")

    fitted_condition_limit = float(sparse.get("score_condition_limit", 1e12))
    if condition_limit is None:
        selected_condition_limit = fitted_condition_limit
        condition_limit_source = "fitted_score_condition_limit"
    else:
        selected_condition_limit = float(condition_limit)
        condition_limit_source = "explicit_override"
    if not np.isfinite(selected_condition_limit) or selected_condition_limit <= 1:
        raise ValueError("condition_limit must be finite and greater than 1")

    min_score_samples = int(sparse.get("min_score_samples", 2))
    if min_score_samples < 1:
        raise ValueError("fit min_score_samples must be positive")

    times = _effective_times(fit, trajectories, sparse)
    fitted_sample_counts = sparse.get("sample_counts")
    if fitted_sample_counts is not None:
        expected_counts = tuple(int(value) for value in fitted_sample_counts)
        actual_counts = tuple(int(time.size) for time in times)
        if expected_counts != actual_counts:
            raise ValueError(
                "trajectories do not reproduce the effective sample counts used by the fit"
            )

    covariance_out = np.full(
        (len(fit.curve_ids), n_components, n_components),
        np.nan,
        dtype=float,
    )
    standard_errors = np.full(
        (len(fit.curve_ids), n_components),
        np.nan,
        dtype=float,
    )
    prior_covariance = np.diag(retained_values)
    posterior_psd_tolerance = 1e-10 * max(
        1.0,
        float(np.max(np.abs(retained_values))),
    )

    rows: list[dict[str, Any]] = []
    first_failure: tuple[str, str, dict[str, Any]] | None = None

    for index, (curve_id, time) in enumerate(
        zip(fit.curve_ids, times, strict=True)
    ):
        status_code = "ok"
        solve_status = "not_attempted"
        condition_number = np.nan
        minimum_eigenvalue = np.nan
        maximum_eigenvalue = np.nan
        posterior_minimum_eigenvalue = np.nan
        solve_relative_residual = np.nan
        posterior_psd_repair_applied = False

        if time.size < min_score_samples:
            status_code = "curve_too_sparse_for_score_system"
        else:
            covariance_i = evaluate_fitted_covariance(
                grid,
                fitted_covariance,
                time,
            )
            sigma_i = covariance_i + (
                noise_variance + score_ridge
            ) * np.eye(time.size)
            sigma_i = 0.5 * (sigma_i + sigma_i.T)
            sigma_eigenvalues = np.linalg.eigvalsh(sigma_i)
            minimum_eigenvalue = float(np.min(sigma_eigenvalues))
            maximum_eigenvalue = float(np.max(sigma_eigenvalues))
            condition_number = (
                np.inf
                if minimum_eigenvalue <= 0
                else maximum_eigenvalue / minimum_eigenvalue
            )

            if minimum_eigenvalue <= 0:
                status_code = "score_covariance_not_positive_definite"
            elif condition_number > selected_condition_limit:
                status_code = "score_covariance_ill_conditioned"
            else:
                phi_i = np.column_stack(
                    [
                        evaluate_fitted_function(
                            grid,
                            retained_functions[component],
                            time,
                        )
                        for component in range(n_components)
                    ]
                )
                rhs = phi_i @ prior_covariance
                solved = np.linalg.solve(sigma_i, rhs)
                denominator = max(1.0, float(np.linalg.norm(rhs)))
                solve_relative_residual = float(
                    np.linalg.norm(sigma_i @ solved - rhs) / denominator
                )
                posterior = prior_covariance - rhs.T @ solved
                posterior = 0.5 * (posterior + posterior.T)
                posterior_eigenvalues = np.linalg.eigvalsh(posterior)
                posterior_minimum_eigenvalue = float(
                    np.min(posterior_eigenvalues)
                )
                if posterior_minimum_eigenvalue < -posterior_psd_tolerance:
                    status_code = "conditional_score_covariance_not_psd"
                else:
                    if posterior_minimum_eigenvalue < 0:
                        values, vectors = np.linalg.eigh(posterior)
                        values = np.maximum(values, 0.0)
                        posterior = (vectors * values) @ vectors.T
                        posterior = 0.5 * (posterior + posterior.T)
                        posterior_psd_repair_applied = True
                    covariance_out[index] = posterior
                    standard_errors[index] = np.sqrt(
                        np.maximum(np.diag(posterior), 0.0)
                    )
                    solve_status = "solved"

        row = {
            "curve_id": str(curve_id),
            "n_samples": int(time.size),
            "status_code": status_code,
            "solve_status": solve_status,
            "condition_number": float(condition_number),
            "minimum_eigenvalue": float(minimum_eigenvalue),
            "maximum_eigenvalue": float(maximum_eigenvalue),
            "posterior_minimum_eigenvalue": float(
                posterior_minimum_eigenvalue
            ),
            "posterior_psd_repair_applied": bool(
                posterior_psd_repair_applied
            ),
            "solve_relative_residual": float(solve_relative_residual),
            "score_ridge": score_ridge,
            "condition_limit": selected_condition_limit,
        }
        rows.append(row)
        if status_code != "ok" and first_failure is None:
            first_failure = (status_code, str(curve_id), dict(row))

    diagnostics = pd.DataFrame(rows)
    if first_failure is not None and failure_action == "error":
        code, curve_id, details = first_failure
        raise SparseNativeError(
            code,
            f"conditional PACE uncertainty system failed for curve {curve_id!r}",
            details={
                **details,
                "all_curve_diagnostics": diagnostics.to_dict(orient="records"),
            },
        )

    provenance = {
        "method": _METHOD,
        "covariance_source": "full_fitted_covariance_plus_noise_and_declared_score_ridge",
        "fitted_covariance_source": "SparseFPCAResult.covariance",
        "fitted_noise_source": "SparseFPCAResult.noise_variance",
        "eigensystem_source": "SparseFPCAResult.eigenvalues_and_eigenfunctions",
        "native_observation_times": True,
        "raw_sparse_trajectory_interpolation_performed": False,
        "score_ridge": score_ridge,
        "condition_limit": selected_condition_limit,
        "condition_limit_source": condition_limit_source,
        "failure_action": failure_action,
        "posterior_psd_tolerance": posterior_psd_tolerance,
        "population_estimation_uncertainty_included": False,
        "mean_estimation_uncertainty_included": False,
        "covariance_estimation_uncertainty_included": False,
        "eigensystem_estimation_uncertainty_included": False,
        "noise_variance_estimation_uncertainty_included": False,
        "bandwidth_uncertainty_included": False,
        "interpretation": (
            "conditional on fitted mean/covariance/eigensystem and measurement-error "
            "variance; excludes population-estimation uncertainty"
        ),
    }
    return SparseFPCAScoreUncertaintyResult(
        reference=fit,
        curve_ids=tuple(fit.curve_ids),
        covariance=covariance_out,
        standard_errors=standard_errors,
        diagnostics=diagnostics,
        n_components=n_components,
        provenance=provenance,
    )

def sparse_fpca_score_uncertainty_frame(
    result: SparseFPCAScoreUncertaintyResult,
) -> pd.DataFrame:
    """Return one tidy row per curve and retained score component."""

    if not isinstance(result, SparseFPCAScoreUncertaintyResult):
        raise TypeError("result must be a SparseFPCAScoreUncertaintyResult")
    diagnostic = result.diagnostics.set_index("curve_id", drop=False)
    rows: list[dict[str, Any]] = []
    for curve_index, curve_id in enumerate(result.curve_ids):
        curve_diagnostic = diagnostic.loc[str(curve_id)]
        for component in range(result.n_components):
            rows.append(
                {
                    "curve_id": str(curve_id),
                    "component": component + 1,
                    "conditional_variance": float(
                        result.covariance[curve_index, component, component]
                    ),
                    "conditional_standard_error": float(
                        result.standard_errors[curve_index, component]
                    ),
                    "status_code": str(curve_diagnostic["status_code"]),
                    "condition_number": float(
                        curve_diagnostic["condition_number"]
                    ),
                    "method": result.method,
                }
            )
    return pd.DataFrame(rows)

def sparse_fpca_score_uncertainty_reporting_text(
    result: SparseFPCAScoreUncertaintyResult,
) -> str:
    """Return manuscript-ready scope wording for conditional score uncertainty."""

    if not isinstance(result, SparseFPCAScoreUncertaintyResult):
        raise TypeError("result must be a SparseFPCAScoreUncertaintyResult")

    diagnostics = result.diagnostics
    if "status_code" not in diagnostics.columns:
        raise ValueError("result diagnostics must contain status_code")
    failed = int(np.count_nonzero(diagnostics["status_code"].to_numpy() != "ok"))
    n_curves = len(result.curve_ids)
    score_ridge = float(result.provenance.get("score_ridge", 0.0))

    text = (
        f"Conditional PACE score uncertainty was computed for {n_curves} curves "
        f"and {result.n_components} retained component(s) using the full fitted "
        "covariance evaluated at native observation times. The covariance is "
        "conditional on the fitted mean, covariance surface, eigensystem, and "
        "measurement-error variance and therefore excludes population-estimation "
        "uncertainty, including uncertainty from estimating smoothing bandwidths "
        "and the noise variance."
    )
    if score_ridge > 0:
        text += (
            f" The fitted PACE score ridge ({score_ridge:g}) was inherited as "
            "numerical regularization of the score system and was not interpreted "
            "as additional measurement-error variance."
        )
    else:
        text += " No score ridge was added to the conditional system."

    if failed:
        text += (
            f" {failed} of {n_curves} curve-level uncertainty systems failed the "
            "declared numerical diagnostics and were retained explicitly rather "
            "than silently dropped."
        )
    else:
        text += " All curve-level conditional uncertainty systems passed the declared numerical diagnostics."
    return text

def select_sparse_fpca_bandwidths(
    trajectories: IrregularTrajectorySet,
    *,
    dimension: str,
    evaluation_grid: np.ndarray,
    mean_bandwidths: Sequence[float],
    covariance_bandwidths: Sequence[float],
    noise_bandwidth: float | None = None,
    noise_support: tuple[float, float] | None = None,
    analysis_support_action: str = "error",
    noise_variance_method: str = "diagonal_difference",
    measurement_error_variance: float | None = None,
    n_splits: int = 5,
    resampling_unit: str = "curve",
    group_column: str | None = None,
    shuffle: bool = True,
    random_state: int | None = 0,
    min_valid_folds: int | None = None,
    failure_action: str = "retain",
    predictive_condition_limit: float = 1e12,
    psd_action: str = "error",
    psd_tolerance: float = 1e-8,
    positive_eigen_tolerance: float = 1e-10,
    mean_smoother: str = "local_linear",
    covariance_smoother: str = "local_linear",
    kernel: str = "epanechnikov",
    mean_min_local_points: int = 3,
    covariance_min_local_pairs: int = 6,
    noise_min_local_points: int = 3,
) -> SparseFPCABandwidthSelectionResult:
    """Select sparse-FPCA mean/covariance bandwidths by held-out likelihood.

    Each candidate is fitted from scratch inside every training fold. The loss
    for a held-out curve is its per-observation Gaussian negative log predictive
    density under the training-only fitted mean, full latent covariance, and
    measurement-error variance. Fold loss is the mean of the held-out curve
    losses, so curves rather than individual observations are the primary
    validation unit.

    ``noise_bandwidth`` is not tuned in this first qualified tranche. When
    diagonal-difference noise estimation is requested it must be supplied once
    and is held fixed across all mean/covariance candidates and folds.

    This selector never changes :func:`fit_sparse_fpca` defaults and never
    invokes automatic tuning from inside the fitter. The selected numeric
    bandwidths must be passed explicitly to a subsequent fit.
    """

    mean_grid = _positive_grid(mean_bandwidths, name="mean_bandwidths")
    covariance_grid = _positive_grid(
        covariance_bandwidths,
        name="covariance_bandwidths",
    )
    grid, required_folds, resolved_noise_bandwidth = _validate_selector_inputs(
        trajectories,
        dimension=dimension,
        evaluation_grid=evaluation_grid,
        noise_variance_method=noise_variance_method,
        measurement_error_variance=measurement_error_variance,
        noise_bandwidth=noise_bandwidth,
        noise_support=noise_support,
        analysis_support_action=analysis_support_action,
        n_splits=n_splits,
        resampling_unit=resampling_unit,
        group_column=group_column,
        min_valid_folds=min_valid_folds,
        predictive_condition_limit=predictive_condition_limit,
        failure_action=failure_action,
    )
    dimension_index = trajectories.dimension_names.index(dimension)
    splits, assignments, groups = _make_folds(
        trajectories,
        n_splits=n_splits,
        resampling_unit=resampling_unit,
        group_column=group_column,
        shuffle=shuffle,
        random_state=random_state,
    )
    candidates = _candidate_table(
        mean_grid,
        covariance_grid,
        resolved_noise_bandwidth,
    )

    fold_rows: list[dict[str, Any]] = []
    curve_rows: list[dict[str, Any]] = []

    for candidate in candidates.to_dict(orient="records"):
        candidate_id = str(candidate["candidate_id"])
        mean_bandwidth = float(candidate["mean_bandwidth"])
        covariance_bandwidth = float(candidate["covariance_bandwidth"])

        for fold, (train_idx, validation_idx) in enumerate(splits):
            train = trajectories.subset(train_idx)
            validation = trajectories.subset(validation_idx)
            raw_train_observations = int(np.sum(train.sample_counts))
            raw_validation_observations = int(np.sum(validation.sample_counts))
            effective_train_observations = _effective_observation_count(
                train,
                grid,
                action=analysis_support_action,
            )
            effective_validation_observations = _effective_observation_count(
                validation,
                grid,
                action=analysis_support_action,
            )
            base_row: dict[str, Any] = {
                "candidate_id": candidate_id,
                "fold": int(fold),
                "mean_bandwidth": mean_bandwidth,
                "covariance_bandwidth": covariance_bandwidth,
                "noise_bandwidth": resolved_noise_bandwidth,
                "n_train_curves": int(train.n_curves),
                "n_validation_curves": int(validation.n_curves),
                "n_train_observations_raw": raw_train_observations,
                "n_validation_observations_raw": raw_validation_observations,
                "n_train_observations_effective": effective_train_observations,
                "n_validation_observations_effective": effective_validation_observations,
            }
            try:
                fit = fit_sparse_fpca(
                    train,
                    dimension=dimension,
                    n_components=1,
                    evaluation_grid=grid,
                    mean_bandwidth=mean_bandwidth,
                    covariance_bandwidth=covariance_bandwidth,
                    noise_bandwidth=resolved_noise_bandwidth,
                    noise_support=noise_support,
                    analysis_support_action=analysis_support_action,
                    mean_smoother=mean_smoother,
                    covariance_smoother=covariance_smoother,
                    kernel=kernel,
                    noise_variance_method=noise_variance_method,
                    measurement_error_variance=measurement_error_variance,
                    psd_action=psd_action,
                    psd_tolerance=psd_tolerance,
                    positive_eigen_tolerance=positive_eigen_tolerance,
                    score_failure_action="retain_nan",
                    mean_min_local_points=mean_min_local_points,
                    covariance_min_local_pairs=covariance_min_local_pairs,
                    noise_min_local_points=noise_min_local_points,
                )
            except (SparseNativeError, ValueError, np.linalg.LinAlgError) as exc:
                code, message = _failure_details(exc)
                fold_rows.append(
                    {
                        **base_row,
                        "status_code": "fit_failure",
                        "failure_code": code,
                        "failure_message": message,
                        "mean_curve_gaussian_nll": np.nan,
                        "median_curve_gaussian_nll": np.nan,
                        "n_evaluated_curves": 0,
                        "n_evaluated_observations": 0,
                        "fitted_noise_variance": np.nan,
                    }
                )
                if failure_action == "error":
                    raise
                continue

            validation_losses: list[float] = []
            validation_observations = 0
            validation_failed = False
            first_validation_failure: tuple[str, str] | None = None
            for curve_id, time, values in zip(
                validation.curve_ids,
                validation.time,
                validation.values,
                strict=True,
            ):
                observed = np.asarray(values[:, dimension_index], dtype=float)
                raw_curve_observations = int(len(time))
                try:
                    effective_time, effective_observed, outside = (
                        _effective_validation_curve(
                            time,
                            observed,
                            grid,
                            action=analysis_support_action,
                        )
                    )
                    loss, condition, minimum, maximum = _curve_gaussian_nll(
                        fit,
                        effective_time,
                        effective_observed,
                        predictive_condition_limit=predictive_condition_limit,
                    )
                    validation_losses.append(loss)
                    validation_observations += int(effective_time.size)
                    curve_rows.append(
                        {
                            "candidate_id": candidate_id,
                            "fold": int(fold),
                            "curve_id": str(curve_id),
                            "status_code": "ok",
                            "failure_code": None,
                            "failure_message": None,
                            "gaussian_nll_per_observation": float(loss),
                            "n_observations_raw": raw_curve_observations,
                            "n_evaluated_observations": int(effective_time.size),
                            "outside_support_observations": int(outside),
                            "predictive_condition_number": float(condition),
                            "predictive_minimum_eigenvalue": float(minimum),
                            "predictive_maximum_eigenvalue": float(maximum),
                        }
                    )
                except (SparseNativeError, ValueError, np.linalg.LinAlgError) as exc:
                    validation_failed = True
                    code, message = _failure_details(exc)
                    if first_validation_failure is None:
                        first_validation_failure = (code, message)
                    curve_rows.append(
                        {
                            "candidate_id": candidate_id,
                            "fold": int(fold),
                            "curve_id": str(curve_id),
                            "status_code": "validation_failure",
                            "failure_code": code,
                            "failure_message": message,
                            "gaussian_nll_per_observation": np.nan,
                            "n_observations_raw": raw_curve_observations,
                            "n_evaluated_observations": 0,
                            "outside_support_observations": np.nan,
                            "predictive_condition_number": np.nan,
                            "predictive_minimum_eigenvalue": np.nan,
                            "predictive_maximum_eigenvalue": np.nan,
                        }
                    )
                    if failure_action == "error":
                        raise

            if validation_failed:
                code, message = first_validation_failure or (
                    "validation_failure",
                    "held-out validation failed",
                )
                fold_rows.append(
                    {
                        **base_row,
                        "status_code": "validation_failure",
                        "failure_code": code,
                        "failure_message": message,
                        "mean_curve_gaussian_nll": np.nan,
                        "median_curve_gaussian_nll": np.nan,
                        "n_evaluated_curves": int(len(validation_losses)),
                        "n_evaluated_observations": int(validation_observations),
                        "fitted_noise_variance": float(fit.noise_variance),
                    }
                )
                continue

            losses = np.asarray(validation_losses, dtype=float)
            fold_rows.append(
                {
                    **base_row,
                    "status_code": "ok",
                    "failure_code": None,
                    "failure_message": None,
                    "mean_curve_gaussian_nll": float(np.mean(losses)),
                    "median_curve_gaussian_nll": float(np.median(losses)),
                    "n_evaluated_curves": int(losses.size),
                    "n_evaluated_observations": int(validation_observations),
                    "fitted_noise_variance": float(fit.noise_variance),
                }
            )

    fold_results = pd.DataFrame(fold_rows)
    curve_losses = pd.DataFrame(curve_rows)
    summary = _summarise_candidates(
        candidates,
        fold_results,
        min_valid_folds=required_folds,
    )
    selected = _select_minimum_candidate(summary)
    status_code = "ok" if selected is not None else "no_eligible_candidate"

    provenance = {
        "method": _METHOD,
        "criterion": _CRITERION,
        "criterion_formula": (
            "curve mean of 0.5*(logdet(Sigma_i)+r_i^T Sigma_i^-1 r_i+"
            "m_i*log(2*pi))/m_i; fold mean weights held-out curves equally"
        ),
        "validation_distribution": (
            "Gaussian observation model from training-only fitted mean, full latent "
            "covariance, and measurement-error variance"
        ),
        "fit_inside_fold": True,
        "validation_observations_used_for_training": False,
        "validation_observations_used_for_pace_scoring": False,
        "full_fitted_covariance_used_for_validation": True,
        "score_ridge_included_in_validation_covariance": False,
        "resampling_unit": resampling_unit,
        "group_column": group_column if resampling_unit == "group" else None,
        "group_leakage_prevented": resampling_unit == "group",
        "n_splits": int(n_splits),
        "shuffle": bool(shuffle) if resampling_unit == "curve" else False,
        "random_state": (
            random_state if resampling_unit == "curve" and shuffle else None
        ),
        "min_valid_folds": int(required_folds),
        "failure_action": failure_action,
        "predictive_condition_limit": float(predictive_condition_limit),
        "candidate_order": "sorted_unique_mean_x_covariance_cartesian_product",
        "selection_rule": "minimum_mean_fold_loss",
        "tie_break_rule": (
            "minimum loss, then larger mean/covariance bandwidths "
            "lexicographically, then candidate_id"
        ),
        "one_se_rule_implemented": False,
        "one_se_reason": (
            "no predeclared scalar simplicity ordering for a multidimensional "
            "bandwidth tuple"
        ),
        "noise_bandwidth_tuned": False,
        "noise_bandwidth_policy": (
            "fixed_declared_across_candidates_and_folds"
            if resolved_noise_bandwidth is not None
            else "not_applicable_fixed_measurement_error_variance"
        ),
        "automatic_fit_bandwidth_selection_performed": False,
        "selected_values_must_be_passed_explicitly_to_fit_sparse_fpca": True,
        "evaluation_grid": grid.tolist(),
        "analysis_support_action": analysis_support_action,
        "noise_variance_method": noise_variance_method,
        "noise_bandwidth": resolved_noise_bandwidth,
        "noise_support": None if noise_support is None else list(noise_support),
        "measurement_error_variance": (
            None
            if measurement_error_variance is None
            else float(measurement_error_variance)
        ),
        "psd_action": psd_action,
        "psd_tolerance": float(psd_tolerance),
        "positive_eigen_tolerance": float(positive_eigen_tolerance),
        "mean_smoother": mean_smoother,
        "covariance_smoother": covariance_smoother,
        "kernel": kernel,
        "mean_min_local_points": int(mean_min_local_points),
        "covariance_min_local_pairs": int(covariance_min_local_pairs),
        "noise_min_local_points": int(noise_min_local_points),
        "candidate_count": int(len(candidates)),
        "fold_assignment_count": int(len(assignments)),
        "training_population_fit_n_components": 1,
        "training_score_failure_action": "retain_nan",
    }
    if groups is not None:
        provenance["group_count"] = int(len(pd.unique(groups)))

    return SparseFPCABandwidthSelectionResult(
        assignments=assignments,
        fold_results=fold_results,
        curve_losses=curve_losses,
        candidate_summary=summary,
        candidates=candidates,
        selected_bandwidths=selected,
        criterion=_CRITERION,
        resampling_unit=resampling_unit,
        n_splits=n_splits,
        group_column=group_column if resampling_unit == "group" else None,
        min_valid_folds=required_folds,
        status_code=status_code,
        provenance=provenance,
    )

def sparse_fpca_bandwidth_selection_reporting_text(
    result: SparseFPCABandwidthSelectionResult,
) -> str:
    """Return compact manuscript-ready wording for an audited selector result."""

    if not isinstance(result, SparseFPCABandwidthSelectionResult):
        raise TypeError("result must be a SparseFPCABandwidthSelectionResult")
    n_candidates = len(result.candidates)
    failures = int(np.count_nonzero(result.fold_results["status_code"] != "ok"))
    unit = (
        f"grouped by {result.group_column!r}"
        if result.resampling_unit == "group"
        else "at the curve level"
    )
    prefix = (
        f"Sparse-FPCA mean/covariance bandwidths were evaluated using "
        f"{result.n_splits}-fold cross-validation {unit}. Population mean, covariance, "
        "and measurement-noise objects were refitted using training curves only. The "
        "prespecified criterion was mean held-out-curve Gaussian negative log predictive "
        "density, using the full fitted latent covariance plus measurement-error "
        "variance; held-out observations were not used for population fitting or PACE "
        "score estimation. "
        f"The audit evaluated {n_candidates} bandwidth candidate(s) and retained "
        f"{failures} failed candidate-fold evaluation(s)."
    )
    if result.provenance.get("noise_bandwidth") is not None:
        prefix += (
            f" The diagonal-difference noise bandwidth "
            f"{float(result.provenance['noise_bandwidth']):g} was declared in advance "
            "and held fixed rather than tuned by this criterion."
        )
    if result.selected_bandwidths is None:
        return prefix + " No candidate satisfied the declared valid-fold requirement."
    selected = result.selected_bandwidths
    return (
        prefix
        + " The minimum mean fold loss selected mean bandwidth "
        + f"{selected['mean_bandwidth']:g} and covariance bandwidth "
        + f"{selected['covariance_bandwidth']:g}. The selected numeric values were "
        "not applied automatically to the fitter."
    )

def observation_process_data(
    frame: pd.DataFrame,
    *,
    curve_column: str,
    time_column: str,
    observed_column: str,
    group_column: str | None = None,
    candidate_predictors: Sequence[str] = (),
    predictor_sources: Mapping[str, str] | None = None,
    coordinate_columns: Sequence[str] = (),
    time_unit: str = "unknown",
    coordinate_system: str = "unknown",
    provenance: Mapping[str, Any] | None = None,
) -> ObservationProcessData:
    """Validate one explicit row per candidate/scheduled sample.

    Candidate predictors must be available on every candidate row. Raw gaze
    coordinates, when supplied, must be finite on retained rows and genuinely
    missing (NaN) on unretained rows. Complete external state belongs in
    ``candidate_predictors`` instead of being presented as recovered gaze.
    """

    if not isinstance(frame, pd.DataFrame):
        raise TypeError("frame must be a pandas DataFrame")
    if frame.empty:
        raise ValueError("frame must contain at least one candidate sample")

    predictors = tuple(dict.fromkeys(map(str, candidate_predictors)))
    coordinates = tuple(map(str, coordinate_columns))
    if len(coordinates) not in {0, 2}:
        raise ValueError("coordinate_columns must be empty or contain exactly two columns")

    required = [curve_column, time_column, observed_column, *predictors, *coordinates]
    if group_column is not None:
        required.append(group_column)
    _require_columns(frame, required)

    reserved = {curve_column, time_column, observed_column}
    if len(reserved) != 3:
        raise ValueError("curve_column, time_column, and observed_column must be distinct")
    if group_column is not None and group_column in reserved:
        raise ValueError("group_column must be distinct from core denominator columns")

    data = frame.copy(deep=True)
    if data[curve_column].isna().any():
        raise ValueError("curve identifiers must not be missing")
    data[curve_column] = data[curve_column].astype(str)
    if group_column is not None:
        if data[group_column].isna().any():
            raise ValueError("group identifiers must not be missing")
        data[group_column] = data[group_column].astype(str)

    try:
        candidate_time = pd.to_numeric(data[time_column], errors="raise").to_numpy(dtype=float)
    except (TypeError, ValueError) as exc:
        raise ValueError("candidate times must be numeric") from exc
    if not np.all(np.isfinite(candidate_time)):
        raise ValueError("candidate times must be finite")
    data[time_column] = candidate_time
    data[observed_column] = _validate_binary(data[observed_column], name=observed_column)

    duplicates = data.duplicated([curve_column, time_column], keep=False)
    if duplicates.any():
        examples = data.loc[duplicates, [curve_column, time_column]].head(5).to_dict("records")
        raise ValueError(
            "candidate rows must be unique within curve/time; "
            f"duplicate examples: {examples}"
        )
    for curve_id, group in data.groupby(curve_column, sort=False):
        times = group[time_column].to_numpy(dtype=float)
        if times.size > 1 and not np.all(np.diff(times) > 0):
            raise ValueError(
                "candidate times must be strictly increasing in supplied row order "
                f"within curve {curve_id!r}"
            )

    source_map = dict(predictor_sources or {})
    unknown_sources = sorted(set(source_map) - set(predictors))
    if unknown_sources:
        raise ValueError(
            "predictor_sources contains names absent from candidate_predictors: "
            f"{unknown_sources}"
        )
    resolved_sources: dict[str, str] = {}
    for predictor in predictors:
        source = source_map.get(predictor, "candidate_time")
        if source not in _ALLOWED_SOURCES:
            raise ValueError(
                f"predictor source for {predictor!r} must be one of "
                f"{sorted(_ALLOWED_SOURCES)}"
            )
        if data[predictor].isna().any():
            raise ValueError(
                f"candidate-time predictor {predictor!r} is unavailable on some "
                "candidate rows; do not fill contemporaneous missing gaze implicitly"
            )
        resolved_sources[predictor] = source

    observed = data[observed_column].to_numpy(dtype=int).astype(bool)
    for coordinate in coordinates:
        numeric = pd.to_numeric(data[coordinate], errors="coerce").to_numpy(dtype=float)
        if np.any(~np.isfinite(numeric[observed])):
            raise ValueError(
                f"coordinate column {coordinate!r} must be finite on observed rows"
            )
        if np.any(~np.isnan(numeric[~observed])):
            raise ValueError(
                f"coordinate column {coordinate!r} contains values on unobserved rows; "
                "label complete external state as a candidate predictor instead"
            )
        data[coordinate] = numeric

    data = data.sort_values([curve_column, time_column], kind="stable").reset_index(drop=True)
    merged_provenance = dict(provenance or {})
    merged_provenance.update(
        {
            "denominator_source": "explicit_candidate_rows",
            "denominator_reconstructed_from_retained_timestamps": False,
            "nominal_sampling_rate_used_to_infer_missing_rows": False,
            "current_missing_gaze_imputed": False,
            "candidate_predictor_sources": dict(resolved_sources),
        }
    )
    return ObservationProcessData(
        frame=data,
        curve_column=str(curve_column),
        time_column=str(time_column),
        observed_column=str(observed_column),
        group_column=None if group_column is None else str(group_column),
        candidate_predictors=predictors,
        predictor_sources=resolved_sources,
        coordinate_columns=coordinates,
        time_unit=str(time_unit),
        coordinate_system=str(coordinate_system),
        provenance=merged_provenance,
    )

def diagnose_observation_process(
    process: ObservationProcessData,
    *,
    predictors: Sequence[str] | None = None,
    history_predictors: Sequence[str] = (),
    eccentricity_reference: tuple[float, float] | None = None,
    risk_set: str = "all_candidates",
    time_basis: str | None = "linear",
    time_bins: int | Sequence[float] | None = None,
    association_bins: int | Sequence[float] | None = 4,
    low_support_threshold: float = 0.10,
) -> ObservationProcessDiagnosticResult:
    """Describe dependence and support in an explicit observation process.

    No iid sample-level standard errors or p-values are reported. Association
    summaries are descriptive and do not identify MAR/MNAR mechanisms.
    """

    if not isinstance(process, ObservationProcessData):
        raise TypeError("process must be an ObservationProcessData")
    if time_basis not in {None, "linear"}:
        raise ValueError("time_basis must be None or 'linear' in the A3 tranche")
    low_support_threshold = float(low_support_threshold)
    if not np.isfinite(low_support_threshold) or not 0.0 <= low_support_threshold <= 1.0:
        raise ValueError("low_support_threshold must be finite and within [0, 1]")

    requested_predictors = (
        process.candidate_predictors
        if predictors is None
        else tuple(dict.fromkeys(map(str, predictors)))
    )
    unknown = sorted(set(requested_predictors) - set(process.candidate_predictors))
    if unknown:
        raise ValueError(
            "predictors must have been declared as complete candidate-time predictors: "
            f"{unknown}"
        )
    requested_history = tuple(dict.fromkeys(map(str, history_predictors)))
    unknown_history = sorted(set(requested_history) - _HISTORY_NAMES)
    if unknown_history:
        raise ValueError(f"unknown history predictor(s): {unknown_history}")

    if eccentricity_reference is not None:
        reference = np.asarray(eccentricity_reference, dtype=float)
        if reference.shape != (2,) or not np.all(np.isfinite(reference)):
            raise ValueError("eccentricity_reference must contain two finite coordinates")
        resolved_reference = (float(reference[0]), float(reference[1]))
    else:
        resolved_reference = None

    candidate_frame = _derive_history(process, eccentricity_reference=resolved_reference)
    risk_mask = _risk_mask(process, candidate_frame, risk_set)
    evaluated = candidate_frame.loc[risk_mask].copy()
    if evaluated.empty:
        raise ValueError("selected risk_set contains no candidate rows")
    observed = evaluated[process.observed_column].to_numpy(dtype=int)

    curve_summary = _summarize_curves(process)
    group_summary = _summarize_groups(process)
    global_summary = pd.DataFrame(
        [
            {
                "candidate_count": process.n_candidates,
                "observed_count": process.n_observed,
                "missing_count": process.n_missing,
                "observed_fraction": process.observed_fraction,
                "missing_fraction": 1.0 - process.observed_fraction,
                "curve_count": len(process.curve_ids),
                "group_count": (
                    int(process.frame[process.group_column].nunique())
                    if process.group_column is not None
                    else np.nan
                ),
                "risk_set": risk_set,
                "risk_candidate_count": int(len(evaluated)),
                "risk_observed_count": int(observed.sum()),
                "risk_missing_count": int(len(evaluated) - observed.sum()),
            }
        ]
    )

    resolved_time_edges = _resolve_edges(
        process.frame[process.time_column].to_numpy(dtype=float),
        time_bins,
        name="time_bins",
    )
    support_summary = _regional_support(
        process,
        resolved_time_edges,
        low_support_threshold=low_support_threshold,
    )
    if resolved_time_edges is None:
        time_summary = pd.DataFrame(
            columns=[
                "predictor",
                "source",
                "profile_type",
                "level",
                "lower",
                "upper",
                "candidate_count",
                "observed_count",
                "observed_fraction",
            ]
        )
    else:
        time_summary = _binned_summary(
            evaluated[process.time_column].to_numpy(dtype=float),
            observed,
            resolved_time_edges,
            predictor=process.time_column,
            source="candidate_time",
        )

    failures: list[dict[str, Any]] = []
    association_rows: list[dict[str, Any]] = []
    profile_frames: list[pd.DataFrame] = []
    if time_basis == "linear":
        association_rows.append(
            _numeric_association(
                evaluated[process.time_column].to_numpy(dtype=float),
                observed,
                predictor=process.time_column,
                source="candidate_time",
            )
        )

    for predictor in requested_predictors:
        series = evaluated[predictor]
        source = process.predictor_sources[predictor]
        if is_numeric_dtype(series):
            values = pd.to_numeric(series, errors="coerce").to_numpy(dtype=float)
            association_rows.append(
                _numeric_association(values, observed, predictor=predictor, source=source)
            )
            try:
                edges = _resolve_edges(values, association_bins, name="association_bins")
            except ValueError as exc:
                failures.append(
                    {
                        "item": predictor,
                        "source": source,
                        "status_code": "profile_unavailable",
                        "message": str(exc),
                    }
                )
                edges = None
            if edges is not None:
                profile_frames.append(
                    _binned_summary(
                        values,
                        observed,
                        edges,
                        predictor=predictor,
                        source=source,
                    )
                )
        else:
            association_rows.append(
                _categorical_association(series, observed, predictor=predictor, source=source)
            )
            profile_frames.append(
                _categorical_profile(series, observed, predictor=predictor, source=source)
            )

    history_requires_coordinates = {
        "previous_observed_x",
        "previous_observed_y",
        "previous_observed_speed",
        "previous_observed_eccentricity",
    }
    for predictor in requested_history:
        if predictor in history_requires_coordinates and not process.coordinate_columns:
            failures.append(
                {
                    "item": predictor,
                    "source": "history",
                    "status_code": "history_coordinates_unavailable",
                    "message": "raw observed coordinate columns were not supplied",
                }
            )
            continue
        if predictor == "previous_observed_eccentricity" and resolved_reference is None:
            failures.append(
                {
                    "item": predictor,
                    "source": "history",
                    "status_code": "reference_definition_required",
                    "message": "eccentricity_reference is required for this history predictor",
                }
            )
            continue
        values = evaluated[predictor].to_numpy(dtype=float)
        finite = np.isfinite(values)
        if not np.any(finite):
            failures.append(
                {
                    "item": predictor,
                    "source": "history",
                    "status_code": "history_unavailable",
                    "message": "no candidate rows have sufficient prior observed history",
                }
            )
            continue
        history_values = values[finite]
        history_observed = observed[finite]
        row = _numeric_association(
            history_values,
            history_observed,
            predictor=predictor,
            source="history",
        )
        row["n_excluded_insufficient_history"] = int((~finite).sum())
        association_rows.append(row)
        try:
            edges = _resolve_edges(
                history_values,
                association_bins,
                name="association_bins",
            )
        except ValueError as exc:
            failures.append(
                {
                    "item": predictor,
                    "source": "history",
                    "status_code": "profile_unavailable",
                    "message": str(exc),
                }
            )
            edges = None
        if edges is not None:
            profile_frames.append(
                _binned_summary(
                    history_values,
                    history_observed,
                    edges,
                    predictor=predictor,
                    source="history",
                )
            )

    associations = pd.DataFrame(association_rows)
    profiles = (
        pd.concat(profile_frames, ignore_index=True)
        if profile_frames
        else pd.DataFrame(
            columns=[
                "predictor",
                "source",
                "profile_type",
                "level",
                "lower",
                "upper",
                "candidate_count",
                "observed_count",
                "observed_fraction",
            ]
        )
    )
    failure_frame = pd.DataFrame(failures) if failures else _empty_failure_frame()

    provenance = {
        "scope": _SCOPE,
        "denominator_source": "explicit_candidate_rows",
        "time_basis": time_basis,
        "time_bins": None if resolved_time_edges is None else list(resolved_time_edges),
        "association_bins": (
            association_bins
            if isinstance(association_bins, int) or association_bins is None
            else list(map(float, association_bins))
        ),
        "risk_set": risk_set,
        "low_support_threshold": low_support_threshold,
        "time_summary_uses_selected_risk_set": True,
        "support_summary_uses_full_candidate_denominator": True,
        "predictor_sources": dict(process.predictor_sources),
        "history_predictors": list(requested_history),
        "history_uses_past_information_only": True,
        "eccentricity_reference": (
            None if resolved_reference is None else list(resolved_reference)
        ),
        "iid_sample_level_inference_performed": False,
        "p_values_reported": False,
        "missing_mechanism_classification_performed": False,
        "inverse_probability_weighting_performed": False,
        "inverse_intensity_correction_performed": False,
        "current_missing_gaze_imputed": False,
        "sparse_estimator_modified": False,
    }
    return ObservationProcessDiagnosticResult(
        process=process,
        candidate_frame=candidate_frame,
        global_summary=global_summary,
        curve_summary=curve_summary,
        group_summary=group_summary,
        time_summary=time_summary,
        support_summary=support_summary,
        associations=associations,
        profiles=profiles,
        failures=failure_frame,
        risk_set=risk_set,
        time_basis=time_basis,
        association_bins=(
            association_bins
            if isinstance(association_bins, int) or association_bins is None
            else tuple(map(float, association_bins))
        ),
        time_bins=(
            time_bins
            if isinstance(time_bins, int) or time_bins is None
            else tuple(map(float, time_bins))
        ),
        low_support_threshold=low_support_threshold,
        provenance=provenance,
    )