Scientific comparison reports

The comparison oracle in hebog.validation.comparison is independent of the PyBDSF product readers and of Hebog's algorithms, so catalogue, RMS and mask equivalence are testable before any frozen external product counts as evidence. The frozen PyBDSF products under config/baselines/ and every Hebog comparison use the same typed reports; stratification by source class is layered on top and never changes these calculations.

Catalogue matching

CatalogueSource stores ICRS positions in degrees, peak flux density in Jy/beam and integrated flux density in Jy, and can carry candidate-reported one-sigma errors, fitted and deconvolved CatalogueEllipse records, an explicit resolved, major-axis-only, unresolved or unavailable state, parent-island identity, component count and canonical quality flags; from_units accepts degrees or arcseconds and Jy or mJy. The caller decides whether it is comparing sources, Gaussian components or another row type first; the concepts are not interchangeable.

compare_catalogues forms all pairs inside the caller's beam-normalized separation gate, boundary included, and solves one global one-to-one assignment with lexicographic objectives: most pairs inside the gate, then the smallest sum of great-circle separations, then the smallest sum of the symmetric flux difference |candidate - reference| / (candidate + reference). Position decides before flux: two objects sharing a gate are paired by separation even when the crossed pairing agrees better in flux, and the flux errors then show the disagreement. Flux enters each pair's cost with a weight of 10⁻⁹ beam, so it decides only between pairings whose separations differ by less than that, which covers coincident rows and rounding. Ties are resolved deterministically for a given input order, and right ascension wraps at 360°.

The report gives separations in beam FWHM and flux differences as signed fractions of the reference, with median and 95th-percentile summaries. Completeness is matched/reference and reliability matched/candidate; an empty denominator counts as 1.0, so two empty catalogues agree and a candidate-only catalogue has zero reliability. Match-only metrics are None without pairs. For matched rows it adds signed fitted and deconvolved axis differences (the worse of major and minor per row), shortest position-angle differences modulo 180°, classification accuracy, component-count agreement and exact and Jaccard quality-flag agreement, keeping fitted and deconvolved angles separate. The caller supplies the minimum reference axis ratio for position-angle evidence (the compact contract uses 1.1). Reference or injected truth alone selects every population, so a missing candidate shape, classification or parent identity counts as unavailable rather than removing a difficult row; each gate reports eligible, available and availability-fraction fields.

Candidate uncertainties are normalized against the reference value, with eligible and available counts, one-sigma coverage, mean normalized residual and sample standard deviation per metric; a missing uncertainty is an availability failure, not a silent omission. Calibration claims therefore need governed truth as the reference; a PyBDSF comparison measures compatibility only.

Parent-island association is compared after assignment through linear-size contingency summaries: pairwise true-positive, false-positive and false-negative co-association counts with precision, recall, intersection over union and agreement, so the far larger set of unrelated pairs cannot hide a split or merge. A missing candidate parent identity stays in the population as a unique unavailable identity. CatalogueOutlierThresholds has no defaults: callers supply their frozen position, flux and shape limits to obtain a catastrophic-outlier fraction, and the report retains them.

RMS maps

compare_rms_maps takes equal-shaped non-negative RMS arrays in Jy/beam and an optional boolean valid mask. Non-finite values and pixels outside the mask are excluded, and a negative finite RMS inside it is invalid input. Absolute differences use every comparable pixel; fractional differences use only pixels with positive reference RMS, with zero-reference pixels counted explicitly. Empty comparisons return counts and None metrics.

Masks and island labels

compare_masks takes boolean arrays of equal shape and reports true-positive, true-negative, false-positive and false-negative counts with agreement, precision, recall and intersection over union. Two empty masks score 1.0 on all four; a missing candidate-positive class has precision 0.0 when the reference has positive pixels and 1.0 otherwise.

compare_island_labels compares non-negative integer label planes independently of their label values, with zero as background. It builds the sparse positive-overlap graph, separates its components and assigns objects by the number of overlapping pairs and then their intersecting pixels, retaining per-match intersection over union, completeness, reliability, unmatched labels and every split or merge visible in the graph. An optional valid mask excludes pixels first. This object report stops high background agreement from concealing a topologically wrong mask.

Released and pinned-master PyBDSF do not produce identical catalogues on the same image; that divergence is resolved against injected truth, not by treating either as authoritative.

hebog.validation.comparison

Independent scientific comparison reports for governed validation.

Catalogue matching uses canonical degrees and janskys. It maximizes the number of pairs inside the separation gate, then minimizes their total angular separation, and uses integrated-flux agreement only to break ties. Array reports never broadcast inputs and state how many pixels were excluded.

AssociationComparisonReport dataclass

Pairwise parent-association confusion counts and rates.

Source code in src/hebog/validation/comparison.py

@dataclass(frozen=True, slots=True)
class AssociationComparisonReport:
    """Pairwise parent-association confusion counts and rates."""

    matched_source_count: int
    identity_eligible_count: int
    identity_available_count: int
    identity_availability_fraction: float | None
    compared_pair_count: int
    true_positive_pair_count: int
    false_positive_pair_count: int
    false_negative_pair_count: int
    disagreement_pair_count: int
    agreement_fraction: float | None
    precision: float
    recall: float
    intersection_over_union: float

BinomialConfidenceInterval dataclass

Two-sided Wilson score interval for one observed success rate.

Source code in src/hebog/validation/comparison.py

@dataclass(frozen=True, slots=True)
class BinomialConfidenceInterval:
    """Two-sided Wilson score interval for one observed success rate."""

    confidence_level: float
    lower: float
    upper: float

CatalogueComparisonReport dataclass

Assignment, completeness, reliability, position, and flux metrics.

Source code in src/hebog/validation/comparison.py

@dataclass(frozen=True, slots=True)
class CatalogueComparisonReport:
    """Assignment, completeness, reliability, position, and flux metrics."""

    reference_count: int
    candidate_count: int
    matches: tuple[CatalogueMatch, ...]
    unmatched_reference_identifiers: tuple[str, ...]
    unmatched_candidate_identifiers: tuple[str, ...]
    completeness: float
    reliability: float
    median_separation_beam_fwhm: float | None
    percentile_95_separation_beam_fwhm: float | None
    median_absolute_peak_flux_fractional_difference: float | None
    percentile_95_absolute_peak_flux_fractional_difference: float | None
    median_absolute_integrated_flux_fractional_difference: float | None
    percentile_95_absolute_integrated_flux_fractional_difference: float | None
    median_absolute_fitted_axis_fractional_difference: float | None
    percentile_95_absolute_fitted_axis_fractional_difference: float | None
    median_absolute_deconvolved_axis_fractional_difference: float | None
    percentile_95_absolute_deconvolved_axis_fractional_difference: float | None
    median_absolute_fitted_position_angle_difference_degrees: float | None
    percentile_95_absolute_fitted_position_angle_difference_degrees: (
        float | None
    )
    median_absolute_deconvolved_position_angle_difference_degrees: float | None
    percentile_95_absolute_deconvolved_position_angle_difference_degrees: (
        float | None
    )
    unresolved_classification_count: int
    unresolved_classification_accuracy: float | None
    field_availability: tuple[CatalogueFieldAvailabilityReport, ...]
    association: AssociationComparisonReport
    component_count_comparison_count: int
    component_count_agreement_fraction: float | None
    quality_flag_exact_agreement_fraction: float | None
    median_quality_flag_jaccard: float | None
    uncertainty_calibration: tuple[UncertaintyCalibrationReport, ...]
    catastrophic_outlier_thresholds: CatalogueOutlierThresholds | None
    catastrophic_outlier_reference_identifiers: tuple[str, ...]
    catastrophic_outlier_fraction: float | None

CatalogueFieldAvailabilityReport dataclass

Reference-selected availability for one gated candidate field.

Source code in src/hebog/validation/comparison.py

@dataclass(frozen=True, slots=True)
class CatalogueFieldAvailabilityReport:
    """Reference-selected availability for one gated candidate field."""

    metric: _AvailabilityMetric
    eligible_count: int
    available_count: int
    availability_fraction: float | None

CatalogueMatch dataclass

One assigned reference/candidate source pair.

Source code in src/hebog/validation/comparison.py

@dataclass(frozen=True, slots=True)
class CatalogueMatch:
    """One assigned reference/candidate source pair."""

    reference_identifier: str
    candidate_identifier: str
    separation_beam_fwhm: float
    peak_flux_fractional_difference: float
    integrated_flux_fractional_difference: float
    fitted_major_axis_fractional_difference: float | None
    fitted_minor_axis_fractional_difference: float | None
    fitted_position_angle_difference_degrees: float | None
    deconvolved_major_axis_fractional_difference: float | None
    deconvolved_minor_axis_fractional_difference: float | None
    deconvolved_position_angle_difference_degrees: float | None
    unresolved_classification_agrees: bool | None
    component_count_agrees: bool | None
    quality_flag_jaccard: float
    quality_flags_agree: bool

CatalogueOutlierThresholds dataclass

Caller-frozen thresholds defining catastrophic matched-row errors.

Source code in src/hebog/validation/comparison.py

@dataclass(frozen=True, slots=True)
class CatalogueOutlierThresholds:
    """Caller-frozen thresholds defining catastrophic matched-row errors."""

    position_beams: float
    peak_flux_fractional_difference: float
    integrated_flux_fractional_difference: float
    fitted_axis_fractional_difference: float
    deconvolved_axis_fractional_difference: float

    def __post_init__(self) -> None:
        """Require finite positive thresholds without hidden defaults."""
        if any(
            not np.isfinite(value) or value <= 0
            for value in (
                self.position_beams,
                self.peak_flux_fractional_difference,
                self.integrated_flux_fractional_difference,
                self.fitted_axis_fractional_difference,
                self.deconvolved_axis_fractional_difference,
            )
        ):
            raise ValueError("catalogue outlier thresholds must be positive")

__post_init__()

Require finite positive thresholds without hidden defaults.

Source code in src/hebog/validation/comparison.py

def __post_init__(self) -> None:
    """Require finite positive thresholds without hidden defaults."""
    if any(
        not np.isfinite(value) or value <= 0
        for value in (
            self.position_beams,
            self.peak_flux_fractional_difference,
            self.integrated_flux_fractional_difference,
            self.fitted_axis_fractional_difference,
            self.deconvolved_axis_fractional_difference,
        )
    ):
        raise ValueError("catalogue outlier thresholds must be positive")

IslandComparisonReport dataclass

Object matches, splits, merges, and overlap for labelled masks.

Source code in src/hebog/validation/comparison.py

@dataclass(frozen=True, slots=True)
class IslandComparisonReport:
    """Object matches, splits, merges, and overlap for labelled masks."""

    compared_pixel_count: int
    excluded_pixel_count: int
    reference_count: int
    candidate_count: int
    matches: tuple[IslandLabelMatch, ...]
    unmatched_reference_labels: tuple[int, ...]
    unmatched_candidate_labels: tuple[int, ...]
    split_reference_labels: tuple[int, ...]
    merged_candidate_labels: tuple[int, ...]
    completeness: float
    reliability: float
    median_matched_intersection_over_union: float | None
    minimum_matched_intersection_over_union: float | None

IslandLabelMatch dataclass

One overlap-based assignment between two labelled regions.

Source code in src/hebog/validation/comparison.py

@dataclass(frozen=True, slots=True)
class IslandLabelMatch:
    """One overlap-based assignment between two labelled regions."""

    reference_label: int
    candidate_label: int
    intersection_pixel_count: int
    union_pixel_count: int
    intersection_over_union: float

IslandPopulationReport dataclass

Aggregate object non-inferiority metrics across a dataset matrix.

Source code in src/hebog/validation/comparison.py

@dataclass(frozen=True, slots=True)
class IslandPopulationReport:
    """Aggregate object non-inferiority metrics across a dataset matrix."""

    case_count: int
    reference_count: int
    candidate_count: int
    matched_count: int
    completeness: float
    reliability: float
    completeness_confidence_interval: BinomialConfidenceInterval | None
    reliability_confidence_interval: BinomialConfidenceInterval | None
    split_count: int
    merge_count: int
    median_matched_intersection_over_union: float | None
    minimum_matched_intersection_over_union: float | None

MaskComparisonReport dataclass

Pixel confusion matrix and derived mask agreement metrics.

Source code in src/hebog/validation/comparison.py

@dataclass(frozen=True, slots=True)
class MaskComparisonReport:
    """Pixel confusion matrix and derived mask agreement metrics."""

    compared_pixel_count: int
    excluded_pixel_count: int
    true_positive_count: int
    true_negative_count: int
    false_positive_count: int
    false_negative_count: int
    agreement_fraction: float
    precision: float
    recall: float
    intersection_over_union: float

NumericConfidenceInterval dataclass

Two-sided confidence interval for one real-valued statistic.

Source code in src/hebog/validation/comparison.py

@dataclass(frozen=True, slots=True)
class NumericConfidenceInterval:
    """Two-sided confidence interval for one real-valued statistic."""

    confidence_level: float
    lower: float
    upper: float

RmsComparisonReport dataclass

Absolute and relative RMS-map differences in canonical units.

Source code in src/hebog/validation/comparison.py

@dataclass(frozen=True, slots=True)
class RmsComparisonReport:
    """Absolute and relative RMS-map differences in canonical units."""

    compared_pixel_count: int
    relative_pixel_count: int
    excluded_pixel_count: int
    zero_reference_pixel_count: int
    median_absolute_difference_jy_per_beam: float | None
    percentile_95_absolute_difference_jy_per_beam: float | None
    median_absolute_fractional_difference: float | None
    percentile_95_absolute_fractional_difference: float | None

UncertaintyCalibrationDecision dataclass

Reviewed entire-interval decision for one uncertainty stratum.

Source code in src/hebog/validation/comparison.py

@dataclass(frozen=True, slots=True)
class UncertaintyCalibrationDecision:
    """Reviewed entire-interval decision for one uncertainty stratum."""

    status: Literal["pass", "fail", "report-only"]
    failed_metrics: tuple[_UncertaintyDecisionFailure, ...]

UncertaintyCalibrationReport dataclass

Bias, dispersion, and one-sigma coverage for one reported error.

Source code in src/hebog/validation/comparison.py

@dataclass(frozen=True, slots=True)
class UncertaintyCalibrationReport:
    """Bias, dispersion, and one-sigma coverage for one reported error."""

    metric: UncertaintyMetric
    eligible_count: int
    sample_count: int
    availability_fraction: float
    within_one_sigma_count: int
    coverage_fraction: float | None
    mean_normalized_residual: float | None
    sample_standard_deviation: float | None
    coverage_confidence_interval: (
        BinomialConfidenceInterval | NumericConfidenceInterval | None
    )
    mean_confidence_interval: NumericConfidenceInterval | None
    dispersion_confidence_interval: NumericConfidenceInterval | None

aggregate_island_comparisons(reports, *, confidence_level=0.95)

Aggregate a generated or governed island-comparison matrix.

Source code in src/hebog/validation/comparison.py

def aggregate_island_comparisons(
    reports: Sequence[IslandComparisonReport],
    *,
    confidence_level: float = 0.95,
) -> IslandPopulationReport:
    """Aggregate a generated or governed island-comparison matrix."""
    reference_count = sum(report.reference_count for report in reports)
    candidate_count = sum(report.candidate_count for report in reports)
    matched_count = sum(len(report.matches) for report in reports)
    overlaps = np.asarray(
        [
            match.intersection_over_union
            for report in reports
            for match in report.matches
        ],
        dtype=np.float64,
    )
    return IslandPopulationReport(
        case_count=len(reports),
        reference_count=reference_count,
        candidate_count=candidate_count,
        matched_count=matched_count,
        completeness=_safe_classification_fraction(
            matched_count,
            reference_count,
            empty_value=1.0,
        ),
        reliability=_safe_classification_fraction(
            matched_count,
            candidate_count,
            empty_value=1.0 if not reference_count else 0.0,
        ),
        completeness_confidence_interval=wilson_score_interval(
            matched_count,
            reference_count,
            confidence_level=confidence_level,
        ),
        reliability_confidence_interval=wilson_score_interval(
            matched_count,
            candidate_count,
            confidence_level=confidence_level,
        ),
        split_count=sum(
            len(report.split_reference_labels) for report in reports
        ),
        merge_count=sum(
            len(report.merged_candidate_labels) for report in reports
        ),
        median_matched_intersection_over_union=_array_median(overlaps),
        minimum_matched_intersection_over_union=(
            float(np.min(overlaps)) if overlaps.size else None
        ),
    )

compare_catalogues(reference, candidate, *, beam_fwhm_degrees, maximum_separation_beams, outlier_thresholds=None, position_angle_minimum_axis_ratio=None)

Match two catalogues and calculate scientific comparison metrics.

Source code in src/hebog/validation/comparison.py

def compare_catalogues(  # noqa: PLR0913
    reference: Sequence[CatalogueSource],
    candidate: Sequence[CatalogueSource],
    *,
    beam_fwhm_degrees: float,
    maximum_separation_beams: float,
    outlier_thresholds: CatalogueOutlierThresholds | None = None,
    position_angle_minimum_axis_ratio: float | None = None,
) -> CatalogueComparisonReport:
    """Match two catalogues and calculate scientific comparison metrics."""
    if not np.isfinite(beam_fwhm_degrees) or beam_fwhm_degrees <= 0:
        raise ValueError("beam_fwhm_degrees must be positive and finite")
    if (
        not np.isfinite(maximum_separation_beams)
        or maximum_separation_beams <= 0
    ):
        raise ValueError(
            "maximum_separation_beams must be positive and finite"
        )
    shape_is_present = any(
        source.fitted_shape is not None or source.deconvolved_shape is not None
        for source in (*reference, *candidate)
    )
    if position_angle_minimum_axis_ratio is None:
        if shape_is_present:
            raise ValueError(
                "position_angle_minimum_axis_ratio is required when "
                "comparing shapes"
            )
        position_angle_axis_ratio = np.inf
    else:
        if (
            not np.isfinite(position_angle_minimum_axis_ratio)
            or position_angle_minimum_axis_ratio <= 1.0
        ):
            raise ValueError(
                "position_angle_minimum_axis_ratio must be finite and "
                "greater than one"
            )
        position_angle_axis_ratio = position_angle_minimum_axis_ratio
    _validate_catalogue_identifiers(reference, catalogue_name="reference")
    _validate_catalogue_identifiers(candidate, catalogue_name="candidate")

    assignments = _match_indices(
        reference,
        candidate,
        beam_fwhm_degrees=beam_fwhm_degrees,
        maximum_separation_beams=maximum_separation_beams,
    )
    matches: list[CatalogueMatch] = []
    uncertainty_samples: dict[str, list[float]] = {}
    uncertainty_eligible_counts: Counter[str] = Counter()
    field_eligible_counts: Counter[str] = Counter()
    field_available_counts: Counter[str] = Counter()
    matched_reference_indices: set[int] = set()
    matched_candidate_indices: set[int] = set()
    for reference_index, candidate_index, separation_beams in assignments:
        reference_source = reference[reference_index]
        candidate_source = candidate[candidate_index]
        matched_reference_indices.add(reference_index)
        matched_candidate_indices.add(candidate_index)
        fitted_shape_differences = _shape_differences(
            reference_source.fitted_shape,
            candidate_source.fitted_shape,
            position_angle_minimum_axis_ratio=position_angle_axis_ratio,
        )
        deconvolved_shape_differences = _deconvolved_shape_differences(
            reference_source,
            candidate_source,
            position_angle_minimum_axis_ratio=position_angle_axis_ratio,
        )
        _record_field_availability(
            reference_source,
            candidate_source,
            field_eligible_counts,
            field_available_counts,
        )
        uncertainty_eligible_counts.update(
            _uncertainty_eligible_metrics(
                reference_source,
                position_angle_minimum_axis_ratio=position_angle_axis_ratio,
            )
        )
        for metric, sample in normalized_uncertainty_samples(
            reference_source,
            candidate_source,
            position_angle_minimum_axis_ratio=position_angle_axis_ratio,
        ):
            uncertainty_samples.setdefault(metric, []).append(sample)
        matches.append(
            CatalogueMatch(
                reference_identifier=reference_source.identifier,
                candidate_identifier=candidate_source.identifier,
                separation_beam_fwhm=separation_beams,
                peak_flux_fractional_difference=(
                    candidate_source.peak_flux_jy_per_beam
                    - reference_source.peak_flux_jy_per_beam
                )
                / reference_source.peak_flux_jy_per_beam,
                integrated_flux_fractional_difference=(
                    candidate_source.integrated_flux_jy
                    - reference_source.integrated_flux_jy
                )
                / reference_source.integrated_flux_jy,
                fitted_major_axis_fractional_difference=(
                    fitted_shape_differences[0]
                ),
                fitted_minor_axis_fractional_difference=(
                    fitted_shape_differences[1]
                ),
                fitted_position_angle_difference_degrees=(
                    fitted_shape_differences[2]
                ),
                deconvolved_major_axis_fractional_difference=(
                    deconvolved_shape_differences[0]
                ),
                deconvolved_minor_axis_fractional_difference=(
                    deconvolved_shape_differences[1]
                ),
                deconvolved_position_angle_difference_degrees=(
                    deconvolved_shape_differences[2]
                ),
                unresolved_classification_agrees=(
                    _deconvolution_classification_agrees(
                        reference_source,
                        candidate_source,
                    )
                ),
                component_count_agrees=(
                    None
                    if reference_source.component_count is None
                    or candidate_source.component_count is None
                    else reference_source.component_count
                    == candidate_source.component_count
                ),
                quality_flag_jaccard=_quality_flag_jaccard(
                    reference_source,
                    candidate_source,
                ),
                quality_flags_agree=(
                    reference_source.quality_flags
                    == candidate_source.quality_flags
                ),
            )
        )

    separations = [match.separation_beam_fwhm for match in matches]
    peak_flux_differences = [
        abs(match.peak_flux_fractional_difference) for match in matches
    ]
    integrated_flux_differences = [
        abs(match.integrated_flux_fractional_difference) for match in matches
    ]
    fitted_axis_differences = [
        value
        for match in matches
        if (
            value := _maximum_absolute_available(
                (
                    match.fitted_major_axis_fractional_difference,
                    match.fitted_minor_axis_fractional_difference,
                )
            )
        )
        is not None
    ]
    deconvolved_axis_differences = [
        value
        for match in matches
        if (
            value := _maximum_absolute_available(
                (
                    match.deconvolved_major_axis_fractional_difference,
                    match.deconvolved_minor_axis_fractional_difference,
                )
            )
        )
        is not None
    ]
    fitted_position_angle_differences = [
        abs(match.fitted_position_angle_difference_degrees)
        for match in matches
        if match.fitted_position_angle_difference_degrees is not None
    ]
    deconvolved_position_angle_differences = [
        abs(match.deconvolved_position_angle_difference_degrees)
        for match in matches
        if match.deconvolved_position_angle_difference_degrees is not None
    ]
    unresolved_classifications = [
        match.unresolved_classification_agrees
        for match in matches
        if match.unresolved_classification_agrees is not None
    ]
    component_comparisons = [
        match.component_count_agrees
        for match in matches
        if match.component_count_agrees is not None
    ]
    association_report = _association_report(reference, candidate, assignments)
    catastrophic_outliers = (
        tuple(
            match.reference_identifier
            for match in matches
            if _is_catastrophic_outlier(match, outlier_thresholds)
        )
        if outlier_thresholds is not None
        else ()
    )
    match_count = len(matches)
    reference_count = len(reference)
    candidate_count = len(candidate)
    return CatalogueComparisonReport(
        reference_count=reference_count,
        candidate_count=candidate_count,
        matches=tuple(matches),
        unmatched_reference_identifiers=tuple(
            source.identifier
            for index, source in enumerate(reference)
            if index not in matched_reference_indices
        ),
        unmatched_candidate_identifiers=tuple(
            source.identifier
            for index, source in enumerate(candidate)
            if index not in matched_candidate_indices
        ),
        completeness=(
            match_count / reference_count if reference_count else 1.0
        ),
        reliability=(
            match_count / candidate_count if candidate_count else 1.0
        ),
        median_separation_beam_fwhm=_median(separations),
        percentile_95_separation_beam_fwhm=_percentile_95(separations),
        median_absolute_peak_flux_fractional_difference=_median(
            peak_flux_differences
        ),
        percentile_95_absolute_peak_flux_fractional_difference=(
            _percentile_95(peak_flux_differences)
        ),
        median_absolute_integrated_flux_fractional_difference=_median(
            integrated_flux_differences
        ),
        percentile_95_absolute_integrated_flux_fractional_difference=(
            _percentile_95(integrated_flux_differences)
        ),
        median_absolute_fitted_axis_fractional_difference=_median(
            fitted_axis_differences
        ),
        percentile_95_absolute_fitted_axis_fractional_difference=(
            _percentile_95(fitted_axis_differences)
        ),
        median_absolute_deconvolved_axis_fractional_difference=_median(
            deconvolved_axis_differences
        ),
        percentile_95_absolute_deconvolved_axis_fractional_difference=(
            _percentile_95(deconvolved_axis_differences)
        ),
        median_absolute_fitted_position_angle_difference_degrees=_median(
            fitted_position_angle_differences
        ),
        percentile_95_absolute_fitted_position_angle_difference_degrees=(
            _percentile_95(fitted_position_angle_differences)
        ),
        median_absolute_deconvolved_position_angle_difference_degrees=(
            _median(deconvolved_position_angle_differences)
        ),
        percentile_95_absolute_deconvolved_position_angle_difference_degrees=(
            _percentile_95(deconvolved_position_angle_differences)
        ),
        unresolved_classification_count=len(unresolved_classifications),
        unresolved_classification_accuracy=(
            sum(unresolved_classifications) / len(unresolved_classifications)
            if unresolved_classifications
            else None
        ),
        field_availability=_field_availability_reports(
            field_eligible_counts,
            field_available_counts,
        ),
        association=association_report,
        component_count_comparison_count=len(component_comparisons),
        component_count_agreement_fraction=(
            sum(component_comparisons) / len(component_comparisons)
            if component_comparisons
            else None
        ),
        quality_flag_exact_agreement_fraction=(
            sum(match.quality_flags_agree for match in matches) / match_count
            if match_count
            else None
        ),
        median_quality_flag_jaccard=_median(
            [match.quality_flag_jaccard for match in matches]
        ),
        uncertainty_calibration=_uncertainty_calibration_reports(
            uncertainty_samples,
            uncertainty_eligible_counts,
        ),
        catastrophic_outlier_thresholds=outlier_thresholds,
        catastrophic_outlier_reference_identifiers=catastrophic_outliers,
        catastrophic_outlier_fraction=(
            len(catastrophic_outliers) / match_count
            if outlier_thresholds is not None and match_count
            else None
        ),
    )

compare_island_labels(reference, candidate, *, valid_mask=None)

Compare non-negative integer island labels by pixel overlap.

Source code in src/hebog/validation/comparison.py

def compare_island_labels(
    reference: npt.ArrayLike,
    candidate: npt.ArrayLike,
    *,
    valid_mask: npt.ArrayLike | None = None,
) -> IslandComparisonReport:
    """Compare non-negative integer island labels by pixel overlap."""
    reference_array, candidate_array = _as_same_shape_arrays(
        reference,
        candidate,
    )
    reference_labels = _as_island_labels(reference_array)
    candidate_labels = _as_island_labels(candidate_array)
    valid = _valid_region(reference_array.shape, valid_mask)
    compared_pixel_count = int(np.count_nonzero(valid))
    reference_counts = _positive_label_counts(reference_labels, valid)
    candidate_counts = _positive_label_counts(candidate_labels, valid)
    intersections = _label_intersections(
        reference_labels,
        candidate_labels,
        valid,
    )
    assigned = _match_label_overlaps(
        intersections,
        compared_pixel_count=compared_pixel_count,
    )
    matches = tuple(
        IslandLabelMatch(
            reference_label=reference_label,
            candidate_label=candidate_label,
            intersection_pixel_count=intersection,
            union_pixel_count=(
                reference_counts[reference_label]
                + candidate_counts[candidate_label]
                - intersection
            ),
            intersection_over_union=(
                intersection
                / (
                    reference_counts[reference_label]
                    + candidate_counts[candidate_label]
                    - intersection
                )
            ),
        )
        for reference_label, candidate_label, intersection in assigned
    )
    matched_reference = {match.reference_label for match in matches}
    matched_candidate = {match.candidate_label for match in matches}
    reference_degrees: dict[int, int] = {}
    candidate_degrees: dict[int, int] = {}
    for reference_label, candidate_label in intersections:
        reference_degrees[reference_label] = (
            reference_degrees.get(reference_label, 0) + 1
        )
        candidate_degrees[candidate_label] = (
            candidate_degrees.get(candidate_label, 0) + 1
        )
    overlap_values = np.asarray(
        [match.intersection_over_union for match in matches],
        dtype=np.float64,
    )
    return IslandComparisonReport(
        compared_pixel_count=compared_pixel_count,
        excluded_pixel_count=reference_labels.size - compared_pixel_count,
        reference_count=len(reference_counts),
        candidate_count=len(candidate_counts),
        matches=matches,
        unmatched_reference_labels=tuple(
            sorted(set(reference_counts) - matched_reference)
        ),
        unmatched_candidate_labels=tuple(
            sorted(set(candidate_counts) - matched_candidate)
        ),
        split_reference_labels=tuple(
            sorted(
                label
                for label, degree in reference_degrees.items()
                if degree > 1
            )
        ),
        merged_candidate_labels=tuple(
            sorted(
                label
                for label, degree in candidate_degrees.items()
                if degree > 1
            )
        ),
        completeness=_safe_classification_fraction(
            len(matches),
            len(reference_counts),
            empty_value=1.0,
        ),
        reliability=_safe_classification_fraction(
            len(matches),
            len(candidate_counts),
            empty_value=1.0 if not reference_counts else 0.0,
        ),
        median_matched_intersection_over_union=_array_median(overlap_values),
        minimum_matched_intersection_over_union=(
            float(np.min(overlap_values)) if overlap_values.size else None
        ),
    )

compare_masks(reference, candidate, *, valid_mask=None)

Compare two boolean masks over an optional valid pixel region.

Source code in src/hebog/validation/comparison.py

def compare_masks(
    reference: npt.ArrayLike,
    candidate: npt.ArrayLike,
    *,
    valid_mask: npt.ArrayLike | None = None,
) -> MaskComparisonReport:
    """Compare two boolean masks over an optional valid pixel region."""
    reference_array, candidate_array = _as_same_shape_arrays(
        reference,
        candidate,
    )
    if not np.issubdtype(reference_array.dtype, np.bool_) or not np.issubdtype(
        candidate_array.dtype,
        np.bool_,
    ):
        raise TypeError("mask products must be boolean arrays")
    valid = _valid_region(reference_array.shape, valid_mask)
    reference_mask = np.asarray(reference_array, dtype=np.bool_)
    candidate_mask = np.asarray(candidate_array, dtype=np.bool_)

    true_positive = int(
        np.count_nonzero(valid & reference_mask & candidate_mask)
    )
    true_negative = int(
        np.count_nonzero(valid & ~reference_mask & ~candidate_mask)
    )
    false_positive = int(
        np.count_nonzero(valid & ~reference_mask & candidate_mask)
    )
    false_negative = int(
        np.count_nonzero(valid & reference_mask & ~candidate_mask)
    )
    compared = int(np.count_nonzero(valid))
    reference_positive = true_positive + false_negative
    candidate_positive = true_positive + false_positive
    return MaskComparisonReport(
        compared_pixel_count=compared,
        excluded_pixel_count=reference_mask.size - compared,
        true_positive_count=true_positive,
        true_negative_count=true_negative,
        false_positive_count=false_positive,
        false_negative_count=false_negative,
        agreement_fraction=_safe_classification_fraction(
            true_positive + true_negative,
            compared,
            empty_value=1.0,
        ),
        precision=_safe_classification_fraction(
            true_positive,
            candidate_positive,
            empty_value=1.0 if reference_positive == 0 else 0.0,
        ),
        recall=_safe_classification_fraction(
            true_positive,
            reference_positive,
            empty_value=1.0,
        ),
        intersection_over_union=_safe_classification_fraction(
            true_positive,
            true_positive + false_positive + false_negative,
            empty_value=1.0,
        ),
    )

compare_rms_maps(reference_rms_jy_per_beam, candidate_rms_jy_per_beam, *, valid_mask=None)

Compare non-negative RMS maps over finite, optionally masked pixels.

Source code in src/hebog/validation/comparison.py

def compare_rms_maps(
    reference_rms_jy_per_beam: npt.ArrayLike,
    candidate_rms_jy_per_beam: npt.ArrayLike,
    *,
    valid_mask: npt.ArrayLike | None = None,
) -> RmsComparisonReport:
    """Compare non-negative RMS maps over finite, optionally masked pixels."""
    reference_input, candidate_input = _as_same_shape_arrays(
        reference_rms_jy_per_beam,
        candidate_rms_jy_per_beam,
    )
    reference = np.asarray(reference_input, dtype=np.float64)
    candidate = np.asarray(candidate_input, dtype=np.float64)
    selected = _valid_region(reference.shape, valid_mask)
    if np.any(selected & np.isfinite(reference) & (reference < 0)) or np.any(
        selected & np.isfinite(candidate) & (candidate < 0)
    ):
        raise ValueError("finite RMS values must be non-negative")

    comparable = selected & np.isfinite(reference) & np.isfinite(candidate)
    absolute_differences = np.asarray(
        np.abs(candidate[comparable] - reference[comparable]),
        dtype=np.float64,
    )
    positive_reference = comparable & (reference > 0)
    relative_differences = np.asarray(
        np.abs(
            (candidate[positive_reference] - reference[positive_reference])
            / reference[positive_reference]
        ),
        dtype=np.float64,
    )
    compared_count = int(np.count_nonzero(comparable))
    relative_count = int(np.count_nonzero(positive_reference))
    zero_reference_count = int(np.count_nonzero(comparable & (reference == 0)))
    return RmsComparisonReport(
        compared_pixel_count=compared_count,
        relative_pixel_count=relative_count,
        excluded_pixel_count=reference.size - compared_count,
        zero_reference_pixel_count=zero_reference_count,
        median_absolute_difference_jy_per_beam=_array_median(
            absolute_differences
        ),
        percentile_95_absolute_difference_jy_per_beam=(
            _array_percentile_95(absolute_differences)
        ),
        median_absolute_fractional_difference=_array_median(
            relative_differences
        ),
        percentile_95_absolute_fractional_difference=(
            _array_percentile_95(relative_differences)
        ),
    )

evaluate_uncertainty_calibration(report, *, minimum_samples, nominal_coverage, maximum_absolute_coverage_difference, maximum_absolute_mean, minimum_standard_deviation, maximum_standard_deviation)

Apply the reviewed entire-confidence-interval decision rule.

Source code in src/hebog/validation/comparison.py

def evaluate_uncertainty_calibration(  # noqa: PLR0913
    report: UncertaintyCalibrationReport,
    *,
    minimum_samples: int,
    nominal_coverage: float,
    maximum_absolute_coverage_difference: float,
    maximum_absolute_mean: float,
    minimum_standard_deviation: float,
    maximum_standard_deviation: float,
) -> UncertaintyCalibrationDecision:
    """Apply the reviewed entire-confidence-interval decision rule."""
    if minimum_samples < _MINIMUM_MEAN_INTERVAL_SAMPLES:
        raise ValueError("minimum_samples must be at least two")
    if report.sample_count < minimum_samples:
        return UncertaintyCalibrationDecision(
            status="report-only",
            failed_metrics=("insufficient-samples",),
        )
    coverage = report.coverage_confidence_interval
    mean = report.mean_confidence_interval
    dispersion = report.dispersion_confidence_interval
    if coverage is None or mean is None or dispersion is None:
        raise ValueError(
            "eligible calibration report lacks confidence intervals"
        )
    failures: list[_UncertaintyDecisionFailure] = []
    if (
        coverage.lower
        < nominal_coverage - maximum_absolute_coverage_difference
        or coverage.upper
        > nominal_coverage + maximum_absolute_coverage_difference
    ):
        failures.append("coverage")
    if (
        mean.lower < -maximum_absolute_mean
        or mean.upper > maximum_absolute_mean
    ):
        failures.append("mean")
    if (
        dispersion.lower < minimum_standard_deviation
        or dispersion.upper > maximum_standard_deviation
    ):
        failures.append("dispersion")
    return UncertaintyCalibrationDecision(
        status="fail" if failures else "pass",
        failed_metrics=tuple(failures),
    )

normalized_uncertainty_samples(reference, candidate, *, position_angle_minimum_axis_ratio)

Return normalized candidate-minus-reference residuals.

A residual divides an offset by an uncertainty, so both must be the same kind of angle. E_RA is a great-circle error, while differencing two right ascensions gives coordinate degrees, so the RA offset is converted before the division; without that the residual would be 1/cos(dec) too large.

Source code in src/hebog/validation/comparison.py

def normalized_uncertainty_samples(
    reference: CatalogueSource,
    candidate: CatalogueSource,
    *,
    position_angle_minimum_axis_ratio: float,
) -> tuple[tuple[UncertaintyMetric, float], ...]:
    """Return normalized candidate-minus-reference residuals.

    A residual divides an offset by an uncertainty, so both must be the same
    kind of angle. `E_RA` is a great-circle error, while differencing two right
    ascensions gives coordinate degrees, so the RA offset is converted before
    the division; without that the residual would be 1/cos(dec) too large.
    """
    samples: list[tuple[UncertaintyMetric, float]] = []

    def add(
        metric: UncertaintyMetric,
        difference: float,
        uncertainty: float | None,
    ) -> None:
        if uncertainty is not None:
            samples.append((metric, difference / uncertainty))

    add(
        "right-ascension",
        _signed_periodic_difference(
            candidate.right_ascension_degrees,
            reference.right_ascension_degrees,
            period=_FULL_CIRCLE_DEGREES,
        )
        * cos(radians(candidate.declination_degrees)),
        candidate.right_ascension_error_degrees,
    )
    add(
        "declination",
        candidate.declination_degrees - reference.declination_degrees,
        candidate.declination_error_degrees,
    )
    add(
        "peak-flux",
        candidate.peak_flux_jy_per_beam - reference.peak_flux_jy_per_beam,
        candidate.peak_flux_error_jy_per_beam,
    )
    add(
        "integrated-flux",
        candidate.integrated_flux_jy - reference.integrated_flux_jy,
        candidate.integrated_flux_error_jy,
    )
    _add_shape_uncertainty_samples(
        samples,
        prefix="fitted",
        reference=reference.fitted_shape,
        candidate=candidate.fitted_shape,
        position_angle_minimum_axis_ratio=position_angle_minimum_axis_ratio,
    )
    _add_shape_uncertainty_samples(
        samples,
        prefix="deconvolved",
        reference=reference.deconvolved_shape,
        candidate=candidate.deconvolved_shape,
        position_angle_minimum_axis_ratio=position_angle_minimum_axis_ratio,
    )
    return tuple(samples)

uncertainty_calibration_report(metric, normalized_residuals, *, eligible_count, confidence_level, bootstrap_resamples, bootstrap_seed, cluster_ids=None, coverage_interval='wilson-score', mean_interval='student-t', dispersion_interval='scipy-bca-bootstrap-fixed-seed')

Summarize one uncertainty stratum with reviewed intervals.

Source code in src/hebog/validation/comparison.py

def uncertainty_calibration_report(  # noqa: PLR0913
    metric: UncertaintyMetric,
    normalized_residuals: npt.ArrayLike,
    *,
    eligible_count: int,
    confidence_level: float,
    bootstrap_resamples: int,
    bootstrap_seed: int,
    cluster_ids: npt.ArrayLike | None = None,
    coverage_interval: CoverageIntervalMethod = "wilson-score",
    mean_interval: MeanIntervalMethod = "student-t",
    dispersion_interval: DispersionIntervalMethod = (
        "scipy-bca-bootstrap-fixed-seed"
    ),
) -> UncertaintyCalibrationReport:
    """Summarize one uncertainty stratum with reviewed intervals."""
    if metric not in _UNCERTAINTY_METRICS:
        raise ValueError("unsupported uncertainty metric")
    values = np.asarray(normalized_residuals, dtype=np.float64)
    if values.ndim != 1:
        raise ValueError("normalized residuals must be one-dimensional")
    if not np.all(np.isfinite(values)):
        raise ValueError("normalized residuals must be finite")
    if (
        isinstance(eligible_count, bool)
        or eligible_count < 0
        or values.size > eligible_count
    ):
        raise ValueError("eligible count must include every residual sample")
    if not np.isfinite(confidence_level) or not 0 < confidence_level < 1:
        raise ValueError("confidence_level must be finite and in (0, 1)")
    if isinstance(bootstrap_resamples, bool) or bootstrap_resamples < 1:
        raise ValueError("bootstrap_resamples must be positive")
    if isinstance(bootstrap_seed, bool) or bootstrap_seed < 0:
        raise ValueError("bootstrap_seed must be non-negative")

    clustered = (
        coverage_interval == "cluster-robust-student-t"
        or mean_interval == "cluster-robust-student-t"
        or dispersion_interval == "cluster-percentile-bootstrap-fixed-seed"
    )
    cluster_codes = _cluster_codes(values, cluster_ids) if clustered else None

    sample_count = int(values.size)
    within_one_sigma_count = int(np.count_nonzero(np.abs(values) <= 1.0))
    coverage_values = np.asarray(np.abs(values) <= 1.0, dtype=np.float64)
    return UncertaintyCalibrationReport(
        metric=metric,
        eligible_count=eligible_count,
        sample_count=sample_count,
        availability_fraction=(
            sample_count / eligible_count if eligible_count else 0.0
        ),
        within_one_sigma_count=within_one_sigma_count,
        coverage_fraction=(
            within_one_sigma_count / sample_count if sample_count else None
        ),
        mean_normalized_residual=(
            float(np.mean(values)) if sample_count else None
        ),
        sample_standard_deviation=(
            _sample_standard_deviation(values) if sample_count > 1 else None
        ),
        coverage_confidence_interval=(
            _cluster_robust_mean_interval(
                coverage_values,
                cast(npt.NDArray[np.int64], cluster_codes),
                confidence_level=confidence_level,
                bounds=(0.0, 1.0),
            )
            if coverage_interval == "cluster-robust-student-t"
            else wilson_score_interval(
                within_one_sigma_count,
                sample_count,
                confidence_level=confidence_level,
            )
        ),
        mean_confidence_interval=(
            _cluster_robust_mean_interval(
                values,
                cast(npt.NDArray[np.int64], cluster_codes),
                confidence_level=confidence_level,
            )
            if mean_interval == "cluster-robust-student-t"
            else _student_mean_interval(
                values,
                confidence_level=confidence_level,
            )
        ),
        dispersion_confidence_interval=(
            _cluster_dispersion_interval(
                values,
                cast(npt.NDArray[np.int64], cluster_codes),
                confidence_level=confidence_level,
                bootstrap_resamples=bootstrap_resamples,
                bootstrap_seed=bootstrap_seed,
            )
            if dispersion_interval == "cluster-percentile-bootstrap-fixed-seed"
            else _dispersion_interval(
                values,
                confidence_level=confidence_level,
                bootstrap_resamples=bootstrap_resamples,
                bootstrap_seed=bootstrap_seed,
            )
        ),
    )

wilson_score_interval(successes, total, *, confidence_level=0.95)

Return a Wilson interval, or None for an empty population.

Source code in src/hebog/validation/comparison.py

def wilson_score_interval(
    successes: int,
    total: int,
    *,
    confidence_level: float = 0.95,
) -> BinomialConfidenceInterval | None:
    """Return a Wilson interval, or ``None`` for an empty population."""
    if (
        isinstance(successes, bool)
        or isinstance(total, bool)
        or successes < 0
        or total < 0
        or successes > total
    ):
        raise ValueError(
            "binomial counts must satisfy 0 <= successes <= total"
        )
    if not np.isfinite(confidence_level) or not 0 < confidence_level < 1:
        raise ValueError("confidence_level must be finite and in (0, 1)")
    if total == 0:
        return None
    probability = successes / total
    z_score = NormalDist().inv_cdf(0.5 + confidence_level / 2.0)
    z_squared = z_score * z_score
    denominator = 1.0 + z_squared / total
    centre = probability + z_squared / (2.0 * total)
    radius = z_score * np.sqrt(
        probability * (1.0 - probability) / total
        + z_squared / (4.0 * total * total)
    )
    return BinomialConfidenceInterval(
        confidence_level=confidence_level,
        lower=float(max(0.0, (centre - radius) / denominator)),
        upper=float(min(1.0, (centre + radius) / denominator)),
    )