fix: expose effective probability cutoff for PPCR operating points

Expose producer-owned probability_threshold column in prepare_performance_data()
derived via Type-7 linear quantile mapping for PPCR stratification, and populate
operating_points.cutoff from probability_threshold in _prepare_probs_distribution_data().

Co-authored-by: uriahf <11351434+uriahf@users.noreply.github.com>

Author: google-labs-jules[bot]

src/rtichoke/performance_data/performance_data.py

@@ -23,6 +23,7 @@
     "reference_group",
     "stratified_by",
     "chosen_cutoff",
+    "probability_threshold",
     "true_positives",
     "true_negatives",
     "false_positives",
@@ -244,9 +245,11 @@ def prepare_performance_data(
     Returns
     -------
     pl.DataFrame
-        A Polars DataFrame where each row corresponds to a probability cutoff
-        for a given model/dataset. Columns include the cutoff value and a rich
-        set of performance metrics (e.g., `tpr`, `fpr`, `ppv`, `net_benefit`).
+        A Polars DataFrame where each row corresponds to an operating point for a
+        given model/dataset. Columns include `chosen_cutoff` (requested grid value),
+        `probability_threshold` (effective score boundary corresponding to the grid
+        point), `ppcr`, and a rich set of performance metrics (e.g., `sensitivity`,
+        `specificity`, `ppv`, `net_benefit`).
 
     Examples
     --------
@@ -274,6 +277,44 @@ def prepare_performance_data(
 
     performance_data = _turn_cumulative_aj_to_performance_data(cumulative_aj_data)
 
+    prob_threshold_rows = []
+    for group_name, probs_array in probs.items():
+        probs_g = np.asarray(probs_array, dtype=float)
+        unique_cutoffs = performance_data.filter(
+            pl.col("reference_group") == group_name
+        )["chosen_cutoff"].unique()
+        for c in unique_cutoffs:
+            eff_cutoff = float(np.quantile(probs_g, 1.0 - float(c), method="linear"))
+            prob_threshold_rows.append(
+                {
+                    "reference_group": group_name,
+                    "stratified_by": "ppcr",
+                    "chosen_cutoff": float(c),
+                    "probability_threshold": eff_cutoff,
+                }
+            )
+
+    ppcr_threshold_df = pl.DataFrame(
+        prob_threshold_rows,
+        schema={
+            "reference_group": performance_data.schema["reference_group"],
+            "stratified_by": performance_data.schema["stratified_by"],
+            "chosen_cutoff": pl.Float64,
+            "probability_threshold": pl.Float64,
+        },
+    )
+
+    performance_data = performance_data.join(
+        ppcr_threshold_df,
+        on=["reference_group", "stratified_by", "chosen_cutoff"],
+        how="left",
+    ).with_columns(
+        pl.when(pl.col("stratified_by") == "probability_threshold")
+        .then(pl.col("chosen_cutoff"))
+        .otherwise(pl.col("probability_threshold"))
+        .alias("probability_threshold")
+    )
+
     return performance_data.select(_PERFORMANCE_DATA_COLUMNS).sort(
         ["reference_group", "stratified_by", "chosen_cutoff"]
     )

src/rtichoke/performance_data/probs_distribution.py

@@ -202,7 +202,11 @@ def _prepare_probs_distribution_data(
     Returns
     -------
     _PredictionDistributionData
-        TypedDict containing ``bins`` and ``operating_points`` Polars DataFrames.
+        TypedDict containing ``bins``, ``operating_points``, and ``rank_bins``
+        Polars DataFrames. In ``operating_points``, ``value`` is the requested grid
+        point (requested PPCR or probability threshold), ``cutoff`` is the effective
+        predicted-probability boundary (from ``probability_threshold``), and
+        ``realized_ppcr`` is the empirical proportion classified positive.
     """
     if not isinstance(stratified_by, (list, tuple)) or len(stratified_by) != 1:
         raise ValueError(
@@ -270,7 +274,7 @@ def _prepare_probs_distribution_data(
             requested_value = float(
                 row["ppcr"] if stratification_type == "ppcr" else row["chosen_cutoff"]
             )
-            effective_cutoff = float(row["chosen_cutoff"])
+            effective_cutoff = float(row["probability_threshold"])
             n_observations = int(row["n"])
             predicted_positives = int(row["predicted_positives"])
             realized_ppcr = (
@@ -292,7 +296,11 @@ def _prepare_probs_distribution_data(
             )
 
         # Build interval boundaries from effective cutoffs
-        cutoffs = evaluation_performance_data["chosen_cutoff"].to_numpy().astype(float)
+        cutoffs = (
+            evaluation_performance_data["probability_threshold"]
+            .to_numpy()
+            .astype(float)
+        )
         interval_boundaries = np.unique(np.concatenate(([0.0, 1.0], cutoffs)))
         interval_boundaries.sort()
 

tests/test_performance_schema_order.py

@@ -13,6 +13,7 @@
     "reference_group",
     "stratified_by",
     "chosen_cutoff",
+    "probability_threshold",
     "true_positives",
     "true_negatives",
     "false_positives",

tests/test_probs_distribution.py

@@ -512,3 +512,229 @@ def test_primary_golden_fixture_order_invariance():
     res_perm = _prepare_probs_distribution_data({"m1": probs_perm}, reals_perm, by=0.20)
 
     assert res_orig["rank_bins"].equals(res_perm["rank_bins"])
+
+
+def test_frozen_r_oracle_ppcr_tied_fixture():
+    probs = {"m1": np.array([0.00, 0.15, 0.30, 0.50, 0.50, 0.50, 0.65, 0.80, 1.00])}
+    reals = np.array([0, 1, 0, 1, 0, 1, 1, 0, 1])
+
+    perf_df = prepare_performance_data(probs, reals, stratified_by=("ppcr",), by=0.20)
+    dist = _prepare_probs_distribution_data(
+        probs, reals, stratified_by=("ppcr",), by=0.20
+    )
+    ops = dist["operating_points"]
+
+    # Independently hardcoded frozen R oracle expected values
+    expected = [
+        # req_ppcr, eff_cutoff, realized_ppcr, tp, fp, tn, fn, sens, spec, ppv, npv, lift
+        (0.00, 1.00, 0.00 / 9.0, 0, 0, 4, 5, 0.0, 1.0, np.nan, 4 / 9, np.nan),
+        (0.20, 0.71, 2.00 / 9.0, 1, 1, 3, 4, 0.2, 0.75, 0.5, 3 / 7, 0.9),
+        (0.40, 0.50, 3.00 / 9.0, 2, 1, 3, 3, 0.4, 0.75, 2 / 3, 0.5, 1.2),
+        (0.60, 0.50, 3.00 / 9.0, 2, 1, 3, 3, 0.4, 0.75, 2 / 3, 0.5, 1.2),
+        (0.80, 0.24, 9.00 / 9.0, 5, 4, 0, 0, 1.0, 0.0, 5 / 9, np.nan, 1.0),
+        (1.00, 0.00, 9.00 / 9.0, 5, 4, 0, 0, 1.0, 0.0, 5 / 9, np.nan, 1.0),
+    ]
+
+    for (
+        req_ppcr,
+        eff_cutoff,
+        real_ppcr,
+        tp,
+        fp,
+        tn,
+        fn,
+        sens,
+        spec,
+        ppv,
+        npv,
+        lift,
+    ) in expected:
+        p_row = perf_df.filter(pl.col("chosen_cutoff") == req_ppcr).row(0, named=True)
+        op_row = ops.filter(pl.col("value") == req_ppcr).row(0, named=True)
+
+        assert p_row["ppcr"] == pytest.approx(req_ppcr)
+        assert p_row["chosen_cutoff"] == pytest.approx(req_ppcr)
+        assert p_row["probability_threshold"] == pytest.approx(eff_cutoff, abs=1e-4)
+
+        assert op_row["type"] == "ppcr"
+        assert op_row["value"] == pytest.approx(req_ppcr)
+        assert op_row["cutoff"] == pytest.approx(eff_cutoff, abs=1e-4)
+        assert op_row["realized_ppcr"] == pytest.approx(real_ppcr)
+
+        assert p_row["true_positives"] == tp
+        assert p_row["false_positives"] == fp
+        assert p_row["true_negatives"] == tn
+        assert p_row["false_negatives"] == fn
+
+        assert p_row["sensitivity"] == pytest.approx(sens)
+        assert p_row["specificity"] == pytest.approx(spec)
+
+        if np.isnan(ppv):
+            assert np.isnan(p_row["ppv"])
+        else:
+            assert p_row["ppv"] == pytest.approx(ppv)
+
+        if np.isnan(npv):
+            assert np.isnan(p_row["npv"])
+        else:
+            assert p_row["npv"] == pytest.approx(npv)
+
+        if np.isnan(lift):
+            assert np.isnan(p_row["lift"])
+        else:
+            assert p_row["lift"] == pytest.approx(lift)
+
+    # Concrete assertion that operating_point.cutoff != operating_point.value for tied PPCR point
+    op_02 = ops.filter(pl.col("value") == 0.20).row(0, named=True)
+    assert op_02["cutoff"] != op_02["value"]
+    assert op_02["cutoff"] == pytest.approx(0.71)
+    assert op_02["value"] == pytest.approx(0.20)
+
+
+def test_ppcr_repeated_effective_cutoffs_under_ties():
+    probs = {"m1": np.array([0.00, 0.15, 0.30, 0.50, 0.50, 0.50, 0.65, 0.80, 1.00])}
+    reals = np.array([0, 1, 0, 1, 0, 1, 1, 0, 1])
+
+    dist = _prepare_probs_distribution_data(
+        probs, reals, stratified_by=("ppcr",), by=0.20
+    )
+    ops = dist["operating_points"]
+
+    op_04 = ops.filter(pl.col("value") == 0.40).row(0, named=True)
+    op_06 = ops.filter(pl.col("value") == 0.60).row(0, named=True)
+
+    # Different requested value
+    assert op_04["value"] == pytest.approx(0.40)
+    assert op_06["value"] == pytest.approx(0.60)
+    assert op_04["value"] != op_06["value"]
+
+    # Same effective cutoff, realized PPCR, and confusion matrix
+    assert op_04["cutoff"] == pytest.approx(op_06["cutoff"]) == pytest.approx(0.50)
+    assert (
+        op_04["realized_ppcr"]
+        == pytest.approx(op_06["realized_ppcr"])
+        == pytest.approx(3 / 9)
+    )
+
+
+def test_ppcr_n_less_than_q():
+    probs = {"m1": np.array([0.10, 0.50, 0.90])}
+    reals = np.array([0, 1, 1])
+
+    res = _prepare_probs_distribution_data(
+        probs, reals, stratified_by=("ppcr",), by=0.20
+    )
+    rb = res["rank_bins"]
+    ops = res["operating_points"]
+
+    # Explicit empty rank strata preserved (q = 5)
+    assert len(rb) == 5
+    assert rb["n_positive"].to_list() == [0, 0, 1, 0, 1]
+    assert rb["n_negative"].to_list() == [1, 0, 0, 0, 0]
+
+    # Verify operating point cutoffs for N < q
+    cutoffs = ops["cutoff"].to_list()
+    assert len(cutoffs) == 6
+    assert np.allclose(cutoffs, [0.9, 0.74, 0.58, 0.42, 0.26, 0.1])
+
+
+def test_ppcr_distinct_predictions_type7_boundaries():
+    probs = {"m1": np.array([0.10, 0.30, 0.70, 0.90])}
+    reals = np.array([0, 1, 0, 1])
+
+    res = _prepare_probs_distribution_data(
+        probs, reals, stratified_by=("ppcr",), by=0.25
+    )
+    ops = res["operating_points"]
+
+    # Expected Type-7 linear quantiles at 1.0, 0.75, 0.50, 0.25, 0.0
+    expected_cutoffs = [0.90, 0.75, 0.50, 0.25, 0.10]
+    for row, exp_c in zip(ops.iter_rows(named=True), expected_cutoffs):
+        assert row["cutoff"] == pytest.approx(exp_c)
+
+
+def test_ppcr_exact_scores_zero_and_one_endpoints():
+    probs = {"m1": np.array([0.0, 0.0, 0.5, 1.0, 1.0])}
+    reals = np.array([0, 1, 0, 1, 0])
+
+    res = _prepare_probs_distribution_data(
+        probs, reals, stratified_by=("ppcr",), by=0.20
+    )
+    ops = res["operating_points"]
+
+    op_0 = ops.filter(pl.col("value") == 0.0).row(0, named=True)
+    assert op_0["cutoff"] == pytest.approx(1.0)
+    assert op_0["realized_ppcr"] == pytest.approx(0.0)
+
+    op_1 = ops.filter(pl.col("value") == 1.0).row(0, named=True)
+    assert op_1["cutoff"] == pytest.approx(0.0)
+    assert op_1["realized_ppcr"] == pytest.approx(1.0)
+
+
+def test_ppcr_multiple_models_shared_outcomes():
+    probs = {
+        "m1": np.array([0.1, 0.4, 0.6, 0.9]),
+        "m2": np.array([0.2, 0.5, 0.8, 1.0]),
+    }
+    reals = np.array([0, 1, 0, 1])
+
+    res = _prepare_probs_distribution_data(
+        probs, reals, stratified_by=("ppcr",), by=0.5
+    )
+    ops = res["operating_points"]
+
+    ops_m1 = ops.filter(pl.col("model") == "m1")
+    ops_m2 = ops.filter(pl.col("model") == "m2")
+
+    assert len(ops_m1) == 3
+    assert len(ops_m2) == 3
+
+    # Per-model effective cutoffs computed independently
+    m1_c05 = ops_m1.filter(pl.col("value") == 0.5).row(0, named=True)["cutoff"]
+    m2_c05 = ops_m2.filter(pl.col("value") == 0.5).row(0, named=True)["cutoff"]
+
+    assert m1_c05 == pytest.approx(0.5)
+    assert m2_c05 == pytest.approx(0.65)
+
+
+def test_ppcr_multiple_populations_unequal_sample_sizes():
+    probs = {
+        "pop1": np.array([0.1, 0.5, 0.9]),
+        "pop2": np.array([0.2, 0.4, 0.6, 0.8, 1.0]),
+    }
+    reals = {
+        "pop1": np.array([0, 1, 1]),
+        "pop2": np.array([0, 1, 0, 1, 0]),
+    }
+
+    res = _prepare_probs_distribution_data(
+        probs, reals, stratified_by=("ppcr",), by=0.5
+    )
+    ops = res["operating_points"]
+
+    p1_ops = ops.filter(pl.col("population") == "pop1")
+    p2_ops = ops.filter(pl.col("population") == "pop2")
+
+    assert len(p1_ops) == 3
+    assert len(p2_ops) == 3
+
+
+def test_ppcr_pre_post_non_regression():
+    probs = {"m1": np.array([0.00, 0.15, 0.30, 0.50, 0.50, 0.50, 0.65, 0.80, 1.00])}
+    reals = np.array([0, 1, 0, 1, 0, 1, 1, 0, 1])
+
+    perf_df = prepare_performance_data(probs, reals, stratified_by=("ppcr",), by=0.20)
+    dist = _prepare_probs_distribution_data(
+        probs, reals, stratified_by=("ppcr",), by=0.20
+    )
+
+    # Rank bins schema and masses frozen
+    rb = dist["rank_bins"]
+    assert rb["n_positive"].to_list() == [1, 2, 0, 1, 1]
+    assert rb["n_negative"].to_list() == [1, 2, 0, 0, 1]
+
+    # Metrics frozen
+    assert perf_df["true_positives"].to_list() == [0, 1, 2, 2, 5, 5]
+    assert perf_df["false_positives"].to_list() == [0, 1, 1, 1, 4, 4]
+    assert perf_df["true_negatives"].to_list() == [4, 3, 3, 3, 0, 0]
+    assert perf_df["false_negatives"].to_list() == [5, 4, 3, 3, 0, 0]