Buckets:

AsianPlayer
/

TriDrive

Files

xet

AsianPlayer/TriDrive / code /eval /power_gold.py

AsianPlayer

16 days ago

download

raw

8.89 kB

"""cwm3/eval/power_gold.py — how many clean-test POSITIVES the warning-head-vs-DMS claim needs (roadmap §B3 item 4; CPU, seconds).
Hanley–McNeil (1982) variance of an AUROC A with n_pos positives and n_neg negatives:
    Var(A) = [A(1−A) + (n_pos−1)(Q1 − A²) + (n_neg−1)(Q2 − A²)] / (n_pos·n_neg),   Q1 = A/(2−A),   Q2 = 2A²/(1+A)
and for the PAIRED difference of two AUROCs on the same rows Var(Δ) = Var(A1) + Var(A2) − 2·r·SD(A1)·SD(A2) with r the between-score correlation.
The observed device-cluster paired CI of the head − DMS Δ on the 509-row / 63-positive gold has half-width 0.088 (warning_head_oof.md, pooled; manual
0.142, engaged 0.115); the empirical resamples include the device clustering that the formula does not, so the table below is CALIBRATED to that anchor:
    half_width(n_pos) = 0.088 · sqrt(63 / n_pos)     (the Hanley–McNeil positive-limited scaling; n_neg grows in the same proportion and matters little)
and reported next to the uncalibrated formula (r fitted so that the formula reproduces 0.088 at 63 positives, A1 = 0.712 / A2 = 0.637, n_neg = 446).
A paired Δ is declared at the 95 % level when |Δ| > half_width.  Extending the gold with round-2 labels: n_pos = 63 + n_new_pos, where round-2 clips
are drawn from the active strata (results/mining/round2_draw_v1.csv), so the warranted share is expected ABOVE the gold's 12.4 % (63 / 509); the table
gives the number of extra labelled clips needed at 12.4 % (gold rate) and at 25 % (active-strata expectation) for each n_pos.
  python eval/power_gold.py  -> results/tables/power_gold.md (append-only section) + results/power_gold/power_gold_<stamp>.csv
No openpilot DMS field is read.
"""
from __future__ import annotations

import argparse
import math
import sys
import time
from pathlib import Path

import numpy as np
import pandas as pd

CWM3 = Path(__file__).resolve().parents[1]
sys.path.insert(0, str(CWM3 / "common"))
from runmeta import write_run_json                          # noqa: E402

TABLE = CWM3 / "results" / "tables" / "power_gold.md"
OUT = CWM3 / "results" / "power_gold"


def hm_var(A, n_pos, n_neg):
    Q1, Q2 = A / (2 - A), 2 * A * A / (1 + A)
    return (A * (1 - A) + (n_pos - 1) * (Q1 - A * A) + (n_neg - 1) * (Q2 - A * A)) / (n_pos * n_neg)


def paired_half_width(A1, A2, n_pos, n_neg, r):
    v = hm_var(A1, n_pos, n_neg) + hm_var(A2, n_pos, n_neg) - 2 * r * math.sqrt(hm_var(A1, n_pos, n_neg) * hm_var(A2, n_pos, n_neg))
    return 1.96 * math.sqrt(max(v, 0.0))


def fit_r(A1, A2, n_pos, n_neg, target):
    """between-score correlation that makes the Hanley–McNeil paired half-width equal the observed one."""
    lo, hi = -0.99, 0.99
    for _ in range(60):
        mid = 0.5 * (lo + hi)
        if paired_half_width(A1, A2, n_pos, n_neg, mid) > target:
            lo = mid
        else:
            hi = mid
    return 0.5 * (lo + hi)


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--n-pos0", type=int, default=63)
    ap.add_argument("--n-gold0", type=int, default=509)
    ap.add_argument("--hw0", type=float, default=0.088, help="observed pooled paired-Δ CI half-width (device-cluster bootstrap) at n_pos0")
    ap.add_argument("--hw-manual", type=float, default=0.142)
    ap.add_argument("--hw-engaged", type=float, default=0.115)
    ap.add_argument("--n-pos-manual", type=int, default=31)
    ap.add_argument("--n-pos-engaged", type=int, default=32)
    ap.add_argument("--a-head", type=float, default=0.712)
    ap.add_argument("--a-dms", type=float, default=0.637)
    ap.add_argument("--n-pos", default="63,75,88,113,133,195,260")
    ap.add_argument("--deltas", default="0.10,0.075,0.05,0.03")
    ap.add_argument("--pos-rates", default="0.124,0.25", help="expected warranted share among newly labelled clips (gold rate; active-strata expectation)")
    a = ap.parse_args()
    t0 = time.time()
    stamp = time.strftime("%Y%m%dT%H%M%SZ", time.gmtime())
    OUT.mkdir(parents=True, exist_ok=True)
    n_neg0 = a.n_gold0 - a.n_pos0
    r = fit_r(a.a_head, a.a_dms, a.n_pos0, n_neg0, a.hw0)
    n_list = [int(x) for x in a.n_pos.split(",")]
    deltas = [float(x) for x in a.deltas.split(",")]
    rates = [float(x) for x in a.pos_rates.split(",")]
    rows = []
    for n_pos in n_list:
        n_neg = int(round(n_pos * n_neg0 / a.n_pos0))
        hw_cal = a.hw0 * math.sqrt(a.n_pos0 / n_pos)
        hw_hm = paired_half_width(a.a_head, a.a_dms, n_pos, n_neg, r)
        rec = dict(n_pos=n_pos, n_neg_prop=n_neg, half_width_calibrated=hw_cal, half_width_hanley_mcneil=hw_hm, hw_manual_scaled=a.hw_manual * math.sqrt(a.n_pos_manual / max(n_pos * a.n_pos_manual / a.n_pos0, 1)),
                   hw_engaged_scaled=a.hw_engaged * math.sqrt(a.n_pos_engaged / max(n_pos * a.n_pos_engaged / a.n_pos0, 1)), smallest_declarable_delta=hw_cal)
        for d in deltas:
            rec[f"declares_{d:g}"] = bool(d > hw_cal)
        for rate in rates:
            extra_pos = max(n_pos - a.n_pos0, 0)
            rec[f"new_clips_at_{rate:g}"] = int(math.ceil(extra_pos / rate)) if extra_pos else 0
        rows.append(rec)
    for d in deltas:
        n_need = int(math.ceil(a.n_pos0 * (a.hw0 / d) ** 2))
        rows.append(dict(n_pos=n_need, n_neg_prop=int(round(n_need * n_neg0 / a.n_pos0)), half_width_calibrated=a.hw0 * math.sqrt(a.n_pos0 / n_need), half_width_hanley_mcneil=paired_half_width(a.a_head, a.a_dms, n_need, int(round(n_need * n_neg0 / a.n_pos0)), r),
                         hw_manual_scaled=np.nan, hw_engaged_scaled=np.nan, smallest_declarable_delta=d, note=f"n_pos needed to declare Δ = {d:g}", **{f"declares_{dd:g}": bool(dd >= d) for dd in deltas},
                         **{f"new_clips_at_{rate:g}": int(math.ceil(max(n_need - a.n_pos0, 0) / rate)) for rate in rates}))
    df = pd.DataFrame(rows)
    df.to_csv(OUT / f"power_gold_{stamp}.csv", index=False)
    L = [f"## Power of the clean-test gold for the head-vs-DMS paired AUROC claim (eval/power_gold.py, run {stamp})", "",
         f"Anchor: the observed device-cluster paired-Δ CI half-width on the current gold ({a.n_gold0} rows / {a.n_pos0} positives) is {a.hw0:.3f} pooled "
         f"({a.hw_manual:.3f} manual on {a.n_pos_manual} positives, {a.hw_engaged:.3f} engaged on {a.n_pos_engaged}; `warning_head_oof.md`).  Hanley–McNeil positive-limited scaling: half-width(n_pos) = "
         f"{a.hw0:.3f}·√({a.n_pos0}/n_pos); the uncalibrated formula (A_head {a.a_head}, A_DMS {a.a_dms}, negatives in the gold proportion, between-score r fitted to the anchor = {r:+.3f}) is shown for reference. "
         f"A paired Δ is declarable when |Δ| exceeds the half-width.  `new clips` = labelled clips needed for the extra positives at a warranted share of {rates[0]:.1%} (gold rate) / {rates[1]:.0%} (active-strata expectation).", "",
         "| n_pos | half-width (calibrated) | half-width (Hanley–McNeil) | manual / engaged half-width (scaled) | smallest declarable Δ | " + " | ".join(f"Δ={d:g}" for d in deltas) + f" | new clips @{rates[0]:.1%} | new clips @{rates[1]:.0%} | note |",
         "|---|---|---|---|---|" + "---|" * len(deltas) + "---|---|---|"]
    for i, r_ in df.iterrows():
        me = f"{r_.hw_manual_scaled:.3f} / {r_.hw_engaged_scaled:.3f}" if np.isfinite(r_.hw_manual_scaled) else "—"
        note = r_["note"] if "note" in df and isinstance(r_["note"], str) else ""
        L.append(f"| {int(r_.n_pos)} | {r_.half_width_calibrated:.3f} | {r_.half_width_hanley_mcneil:.3f} | {me} | {r_.smallest_declarable_delta:.3f} | "
                 + " | ".join("yes" if bool(r_[f"declares_{d:g}"]) else "no" for d in deltas)
                 + f" | {int(r_[f'new_clips_at_{rates[0]:g}'])} | {int(r_[f'new_clips_at_{rates[1]:g}'])} | {note} |")
    L += ["", f"Reading: today's pooled +0.075 (HGB current vs DMS att@3s) needs n_pos ≥ {int(math.ceil(a.n_pos0 * (a.hw0 / 0.075) ** 2))} to be declared; +0.05 needs ≥ {int(math.ceil(a.n_pos0 * (a.hw0 / 0.05) ** 2))}. "
          f"The engaged stratum alone (32 positives, half-width {a.hw_engaged:.3f}) cannot declare anything below ≈ 0.11 without new engaged positives.  The primary claim stays on the 509-row gold; "
          "the pooled 509 + round-2 stratum (`clean_test_r2`, eval/gold_extended.py) is pre-registered as secondary.", ""]
    md = "\n".join(L) + "\n"
    if not TABLE.exists():
        TABLE.write_text("# Gold power for the warning-head claims (eval/power_gold.py; sections appended per run)\n\n")
    with open(TABLE, "a") as f:
        f.write(md)
    print(md)
    write_run_json(OUT, "power_gold", [Path(__file__)], [], seed=0, t_start=t0, input_listing=False,
                   extra={"anchor": dict(n_pos0=a.n_pos0, n_gold0=a.n_gold0, hw0=a.hw0, hw_manual=a.hw_manual, hw_engaged=a.hw_engaged, a_head=a.a_head, a_dms=a.a_dms, r_fitted=r), "rows": rows, "dms_fields_used": []})


if __name__ == "__main__":
    main()

Xet Storage Details

Size:: 8.89 kB
Xet hash:: cada71dc81259f6a0a9c5baf483a61765b191fffc4d3bc751a3882f5f8117106

Xet efficiently stores files, intelligently splitting them into unique chunks and accelerating uploads and downloads. More info.

	"""cwm3/eval/power_gold.py — how many clean-test POSITIVES the warning-head-vs-DMS claim needs (roadmap §B3 item 4; CPU, seconds).

	Hanley–McNeil (1982) variance of an AUROC A with n_pos positives and n_neg negatives:
	Var(A) = [A(1−A) + (n_pos−1)(Q1 − A²) + (n_neg−1)(Q2 − A²)] / (n_pos·n_neg), Q1 = A/(2−A), Q2 = 2A²/(1+A)
	and for the PAIRED difference of two AUROCs on the same rows Var(Δ) = Var(A1) + Var(A2) − 2·r·SD(A1)·SD(A2) with r the between-score correlation.
	The observed device-cluster paired CI of the head − DMS Δ on the 509-row / 63-positive gold has half-width 0.088 (warning_head_oof.md, pooled; manual
	0.142, engaged 0.115); the empirical resamples include the device clustering that the formula does not, so the table below is CALIBRATED to that anchor:
	half_width(n_pos) = 0.088 · sqrt(63 / n_pos) (the Hanley–McNeil positive-limited scaling; n_neg grows in the same proportion and matters little)
	and reported next to the uncalibrated formula (r fitted so that the formula reproduces 0.088 at 63 positives, A1 = 0.712 / A2 = 0.637, n_neg = 446).
	A paired Δ is declared at the 95 % level when \|Δ\| > half_width. Extending the gold with round-2 labels: n_pos = 63 + n_new_pos, where round-2 clips
	are drawn from the active strata (results/mining/round2_draw_v1.csv), so the warranted share is expected ABOVE the gold's 12.4 % (63 / 509); the table
	gives the number of extra labelled clips needed at 12.4 % (gold rate) and at 25 % (active-strata expectation) for each n_pos.
	python eval/power_gold.py -> results/tables/power_gold.md (append-only section) + results/power_gold/power_gold_<stamp>.csv
	No openpilot DMS field is read.
	"""
	from __future__ import annotations

	import argparse
	import math
	import sys
	import time
	from pathlib import Path

	import numpy as np
	import pandas as pd

	CWM3 = Path(__file__).resolve().parents[1]
	sys.path.insert(0, str(CWM3 / "common"))
	from runmeta import write_run_json # noqa: E402

	TABLE = CWM3 / "results" / "tables" / "power_gold.md"
	OUT = CWM3 / "results" / "power_gold"


	def hm_var(A, n_pos, n_neg):
	Q1, Q2 = A / (2 - A), 2 * A * A / (1 + A)
	return (A * (1 - A) + (n_pos - 1) * (Q1 - A * A) + (n_neg - 1) * (Q2 - A * A)) / (n_pos * n_neg)


	def paired_half_width(A1, A2, n_pos, n_neg, r):
	v = hm_var(A1, n_pos, n_neg) + hm_var(A2, n_pos, n_neg) - 2 * r * math.sqrt(hm_var(A1, n_pos, n_neg) * hm_var(A2, n_pos, n_neg))
	return 1.96 * math.sqrt(max(v, 0.0))


	def fit_r(A1, A2, n_pos, n_neg, target):
	"""between-score correlation that makes the Hanley–McNeil paired half-width equal the observed one."""
	lo, hi = -0.99, 0.99
	for _ in range(60):
	mid = 0.5 * (lo + hi)
	if paired_half_width(A1, A2, n_pos, n_neg, mid) > target:
	lo = mid
	else:
	hi = mid
	return 0.5 * (lo + hi)


	def main():
	ap = argparse.ArgumentParser()
	ap.add_argument("--n-pos0", type=int, default=63)
	ap.add_argument("--n-gold0", type=int, default=509)
	ap.add_argument("--hw0", type=float, default=0.088, help="observed pooled paired-Δ CI half-width (device-cluster bootstrap) at n_pos0")
	ap.add_argument("--hw-manual", type=float, default=0.142)
	ap.add_argument("--hw-engaged", type=float, default=0.115)
	ap.add_argument("--n-pos-manual", type=int, default=31)
	ap.add_argument("--n-pos-engaged", type=int, default=32)
	ap.add_argument("--a-head", type=float, default=0.712)
	ap.add_argument("--a-dms", type=float, default=0.637)
	ap.add_argument("--n-pos", default="63,75,88,113,133,195,260")
	ap.add_argument("--deltas", default="0.10,0.075,0.05,0.03")
	ap.add_argument("--pos-rates", default="0.124,0.25", help="expected warranted share among newly labelled clips (gold rate; active-strata expectation)")
	a = ap.parse_args()
	t0 = time.time()
	stamp = time.strftime("%Y%m%dT%H%M%SZ", time.gmtime())
	OUT.mkdir(parents=True, exist_ok=True)
	n_neg0 = a.n_gold0 - a.n_pos0
	r = fit_r(a.a_head, a.a_dms, a.n_pos0, n_neg0, a.hw0)
	n_list = [int(x) for x in a.n_pos.split(",")]
	deltas = [float(x) for x in a.deltas.split(",")]
	rates = [float(x) for x in a.pos_rates.split(",")]
	rows = []
	for n_pos in n_list:
	n_neg = int(round(n_pos * n_neg0 / a.n_pos0))
	hw_cal = a.hw0 * math.sqrt(a.n_pos0 / n_pos)
	hw_hm = paired_half_width(a.a_head, a.a_dms, n_pos, n_neg, r)
	rec = dict(n_pos=n_pos, n_neg_prop=n_neg, half_width_calibrated=hw_cal, half_width_hanley_mcneil=hw_hm, hw_manual_scaled=a.hw_manual * math.sqrt(a.n_pos_manual / max(n_pos * a.n_pos_manual / a.n_pos0, 1)),
	hw_engaged_scaled=a.hw_engaged * math.sqrt(a.n_pos_engaged / max(n_pos * a.n_pos_engaged / a.n_pos0, 1)), smallest_declarable_delta=hw_cal)
	for d in deltas:
	rec[f"declares_{d:g}"] = bool(d > hw_cal)
	for rate in rates:
	extra_pos = max(n_pos - a.n_pos0, 0)
	rec[f"new_clips_at_{rate:g}"] = int(math.ceil(extra_pos / rate)) if extra_pos else 0
	rows.append(rec)
	for d in deltas:
	n_need = int(math.ceil(a.n_pos0 * (a.hw0 / d) ** 2))
	rows.append(dict(n_pos=n_need, n_neg_prop=int(round(n_need * n_neg0 / a.n_pos0)), half_width_calibrated=a.hw0 * math.sqrt(a.n_pos0 / n_need), half_width_hanley_mcneil=paired_half_width(a.a_head, a.a_dms, n_need, int(round(n_need * n_neg0 / a.n_pos0)), r),
	hw_manual_scaled=np.nan, hw_engaged_scaled=np.nan, smallest_declarable_delta=d, note=f"n_pos needed to declare Δ = {d:g}", **{f"declares_{dd:g}": bool(dd >= d) for dd in deltas},
	**{f"new_clips_at_{rate:g}": int(math.ceil(max(n_need - a.n_pos0, 0) / rate)) for rate in rates}))
	df = pd.DataFrame(rows)
	df.to_csv(OUT / f"power_gold_{stamp}.csv", index=False)
	L = [f"## Power of the clean-test gold for the head-vs-DMS paired AUROC claim (eval/power_gold.py, run {stamp})", "",
	f"Anchor: the observed device-cluster paired-Δ CI half-width on the current gold ({a.n_gold0} rows / {a.n_pos0} positives) is {a.hw0:.3f} pooled "
	f"({a.hw_manual:.3f} manual on {a.n_pos_manual} positives, {a.hw_engaged:.3f} engaged on {a.n_pos_engaged}; `warning_head_oof.md`). Hanley–McNeil positive-limited scaling: half-width(n_pos) = "
	f"{a.hw0:.3f}·√({a.n_pos0}/n_pos); the uncalibrated formula (A_head {a.a_head}, A_DMS {a.a_dms}, negatives in the gold proportion, between-score r fitted to the anchor = {r:+.3f}) is shown for reference. "
	f"A paired Δ is declarable when \|Δ\| exceeds the half-width. `new clips` = labelled clips needed for the extra positives at a warranted share of {rates[0]:.1%} (gold rate) / {rates[1]:.0%} (active-strata expectation).", "",
	"\| n_pos \| half-width (calibrated) \| half-width (Hanley–McNeil) \| manual / engaged half-width (scaled) \| smallest declarable Δ \| " + " \| ".join(f"Δ={d:g}" for d in deltas) + f" \| new clips @{rates[0]:.1%} \| new clips @{rates[1]:.0%} \| note \|",
	"\|---\|---\|---\|---\|---\|" + "---\|" * len(deltas) + "---\|---\|---\|"]
	for i, r_ in df.iterrows():
	me = f"{r_.hw_manual_scaled:.3f} / {r_.hw_engaged_scaled:.3f}" if np.isfinite(r_.hw_manual_scaled) else "—"
	note = r_["note"] if "note" in df and isinstance(r_["note"], str) else ""
	L.append(f"\| {int(r_.n_pos)} \| {r_.half_width_calibrated:.3f} \| {r_.half_width_hanley_mcneil:.3f} \| {me} \| {r_.smallest_declarable_delta:.3f} \| "
	+ " \| ".join("yes" if bool(r_[f"declares_{d:g}"]) else "no" for d in deltas)
	+ f" \| {int(r_[f'new_clips_at_{rates[0]:g}'])} \| {int(r_[f'new_clips_at_{rates[1]:g}'])} \| {note} \|")
	L += ["", f"Reading: today's pooled +0.075 (HGB current vs DMS att@3s) needs n_pos ≥ {int(math.ceil(a.n_pos0 * (a.hw0 / 0.075) ** 2))} to be declared; +0.05 needs ≥ {int(math.ceil(a.n_pos0 * (a.hw0 / 0.05) ** 2))}. "
	f"The engaged stratum alone (32 positives, half-width {a.hw_engaged:.3f}) cannot declare anything below ≈ 0.11 without new engaged positives. The primary claim stays on the 509-row gold; "
	"the pooled 509 + round-2 stratum (`clean_test_r2`, eval/gold_extended.py) is pre-registered as secondary.", ""]
	md = "\n".join(L) + "\n"
	if not TABLE.exists():
	TABLE.write_text("# Gold power for the warning-head claims (eval/power_gold.py; sections appended per run)\n\n")
	with open(TABLE, "a") as f:
	f.write(md)
	print(md)
	write_run_json(OUT, "power_gold", [Path(__file__)], [], seed=0, t_start=t0, input_listing=False,
	extra={"anchor": dict(n_pos0=a.n_pos0, n_gold0=a.n_gold0, hw0=a.hw0, hw_manual=a.hw_manual, hw_engaged=a.hw_engaged, a_head=a.a_head, a_dms=a.a_dms, r_fitted=r), "rows": rows, "dms_fields_used": []})


	if __name__ == "__main__":
	main()