from __future__ import annotations import csv import math from pathlib import Path import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np M = 10000 REPS = 1000 DEFAULT_SEED = 20260925 PI_GRID = (0.001, 0.005, 0.01, 0.05) D_GRID = (0.0, 0.5, 1.0, 2.0) BUDGET_RATIO_GRID = (0.001, 0.005, 0.01, 0.05) BELIEF_LAMBDA = 0.5 NEGATIVE_CONTROL_TOLERANCE = 0.05 ROOT = Path(__file__).resolve().parents[1] RESULT_PATH = ROOT / "results" / "idea_discovery_summary.csv" FIGURE_PATH = ROOT / "figures" / "idea_discovery_summary.png" METHODS = ( ("neutral", 0.0, "neutral (c=0)", "tab:blue", "o"), ("cb_1", 1.0, "c_CB=1", "tab:red", "s"), ("cb_2", 2.0, "c_CB=2", "tab:green", "^"), ) SUMMARY_PREFIXES = ( "recall", "precision", "fdp_single", "yield1000", "validations_per_tp", "invalid_validations_per_tp", "belief_gap", "recall_b0", "recall_b1", "belief_selected_share", ) def make_rng(seed: int, pi_index: int, d_index: int) -> np.random.Generator: sequence = np.random.SeedSequence([int(seed), int(pi_index) + 1, int(d_index) + 1, 314159]) return np.random.Generator(np.random.PCG64(sequence)) def finite_summary(values: list[float], prefix: str) -> dict[str, float]: array = np.asarray(values, dtype=float) valid = np.isfinite(array) fraction = float(np.mean(valid)) if np.any(valid): clean = array[valid] q025, median, q975 = np.quantile(clean, (0.025, 0.5, 0.975)) mean = float(np.mean(clean)) else: q025 = 0.0 median = 0.0 q975 = 0.0 mean = 0.0 return { f"{prefix}_defined_fraction": fraction, f"{prefix}_mean": mean, f"{prefix}_median": float(median), f"{prefix}_p025": float(q025), f"{prefix}_p975": float(q975), } def budgets_from_ratios() -> list[tuple[float, int]]: return [(float(ratio), int(round(ratio * M))) for ratio in BUDGET_RATIO_GRID] def simulate_unit( pi: float, d: float, pi_index: int, d_index: int, seed: int, ) -> list[dict[str, object]]: budget_pairs = budgets_from_ratios() rng = make_rng(seed, pi_index, d_index) records: dict[tuple[str, int], dict[str, list[float]]] = {} for ratio, budget in budget_pairs: for method, c_cb, _, _, _ in METHODS: records[(method, budget)] = {prefix: [] for prefix in SUMMARY_PREFIXES} records[(method, budget)]["tp"] = [] records[(method, budget)]["fp"] = [] records[(method, budget)]["hit_fraction"] = [] records[(method, budget)]["belief_prevalence"] = [] scale = math.sqrt(pi * (1.0 - pi)) for _ in range(REPS): z = (rng.random(M) < pi).astype(np.int8) error = rng.standard_normal(M) score = d * (z.astype(float) - pi) / scale + error score_sd = float(np.std(score)) if score_sd == 0.0: score_sd = 1.0 centered_score = (score - float(np.median(score))) / score_sd belief_logit = BELIEF_LAMBDA * np.tanh(centered_score) belief_probability = np.exp(-np.logaddexp(0.0, -belief_logit)) belief = (rng.random(M) < belief_probability).astype(np.int8) total_signal = int(np.count_nonzero(z)) belief_zero = belief == 0 belief_one = belief == 1 total_belief_zero = int(np.count_nonzero(z[belief_zero])) total_belief_one = int(np.count_nonzero(z[belief_one])) scores = ( ("neutral", 0.0, score), ("cb_1", 1.0, score + belief.astype(float)), ("cb_2", 2.0, score + 2.0 * belief.astype(float)), ) for method, _, method_score in scores: for _, budget in budget_pairs: selected = np.argpartition(method_score, M - budget)[M - budget :] selected_z = z[selected] selected_belief = belief[selected] tp = int(np.count_nonzero(selected_z)) fp = budget - tp precision = tp / float(budget) recall = tp / float(total_signal) if total_signal > 0 else float("nan") fdp_single = fp / float(budget) yield1000 = 1000.0 * precision validations_per_tp = budget / float(tp) if tp > 0 else float("nan") invalid_validations_per_tp = ( (budget - tp) / float(tp) if tp > 0 else float("nan") ) recall_b0 = ( int(np.count_nonzero(selected_z[selected_belief == 0])) / float(total_belief_zero) if total_belief_zero > 0 else float("nan") ) recall_b1 = ( int(np.count_nonzero(selected_z[selected_belief == 1])) / float(total_belief_one) if total_belief_one > 0 else float("nan") ) belief_gap = recall_b1 - recall_b0 if np.isfinite(recall_b1) and np.isfinite(recall_b0) else float("nan") values = { "recall": recall, "precision": precision, "fdp_single": fdp_single, "yield1000": yield1000, "validations_per_tp": validations_per_tp, "invalid_validations_per_tp": invalid_validations_per_tp, "belief_gap": belief_gap, "recall_b0": recall_b0, "recall_b1": recall_b1, "belief_selected_share": float(np.mean(selected_belief)), "tp": float(tp), "fp": float(fp), "hit_fraction": float(tp > 0), "belief_prevalence": float(np.mean(belief)), } record = records[(method, budget)] for key, value in values.items(): record[key].append(value) rows: list[dict[str, object]] = [] for ratio, budget in budget_pairs: for method, c_cb, _, _, _ in METHODS: record = records[(method, budget)] row: dict[str, object] = { "data_kind": "synthetic_scenario", "seed": int(seed), "m": int(M), "reps": int(REPS), "pi": float(pi), "d": float(d), "budget_ratio": float(ratio), "budget": int(budget), "method": method, "c_cb": float(c_cb), "lambda_cb": float(BELIEF_LAMBDA), "negative_control": int(d == 0.0), } for prefix in SUMMARY_PREFIXES: row.update(finite_summary(record[prefix], prefix)) row["tp_mean"] = float(np.mean(record["tp"])) row["fp_mean"] = float(np.mean(record["fp"])) row["hit_fraction"] = float(np.mean(record["hit_fraction"])) row["no_true_hit_fraction"] = 1.0 - float(np.mean(record["hit_fraction"])) row["no_selected_fraction"] = 0.0 row["belief_prevalence_mean"] = float(np.mean(record["belief_prevalence"])) rows.append(row) return rows def write_csv(rows: list[dict[str, object]]) -> None: RESULT_PATH.parent.mkdir(parents=True, exist_ok=True) with RESULT_PATH.open("w", encoding="utf-8", newline="") as handle: writer = csv.DictWriter(handle, fieldnames=list(rows[0].keys())) writer.writeheader() writer.writerows(rows) def find_row( rows: list[dict[str, object]], pi: float, d: float, ratio: float, method: str, ) -> dict[str, object]: matches = [ row for row in rows if row["pi"] == pi and row["d"] == d and row["budget_ratio"] == ratio and row["method"] == method ] if len(matches) != 1: raise RuntimeError(f"row lookup failed: pi={pi}, d={d}, ratio={ratio}, method={method}") return matches[0] def make_figure(rows: list[dict[str, object]]) -> None: FIGURE_PATH.parent.mkdir(parents=True, exist_ok=True) figure, axes = plt.subplots(2, 2, figsize=(13, 9), constrained_layout=True) figure_pi = 0.01 figure_ratio = 0.01 for method, _, label, color, marker in METHODS: x_values = [] y_values = [] for d in D_GRID: row = find_row(rows, figure_pi, d, figure_ratio, method) x_values.append(float(row["fdp_single_mean"])) y_values.append(float(row["recall_mean"])) axes[0, 0].plot(x_values, y_values, color=color, marker=marker, label=label) for x_value, y_value, d in zip(x_values, y_values, D_GRID): axes[0, 0].text(x_value, y_value, f" d={d:g}", fontsize=8) axes[0, 0].set_title(f"Recall versus single-run FDP (pi={figure_pi:g}, B/m={figure_ratio:g})") axes[0, 0].set_xlabel("Mean single-run FDP") axes[0, 0].set_ylabel("Mean Recall") axes[0, 0].set_xlim(-0.01, 1.01) axes[0, 0].set_ylim(-0.02, 1.02) axes[0, 0].legend(fontsize=8) for method, _, label, color, marker in METHODS: for d in D_GRID: line_style = "--" if d == 0.0 else "-" alpha = 0.45 if d == 0.0 else 0.9 points = [find_row(rows, figure_pi, d, ratio, method) for ratio in BUDGET_RATIO_GRID] axes[0, 1].plot( [float(point["budget_ratio"]) for point in points], [float(point["yield1000_mean"]) for point in points], color=color, marker=marker, linestyle=line_style, alpha=alpha, label=f"{label}, d={d:g}", ) axes[0, 1].set_title(f"Budget versus Yield1000 (pi={figure_pi:g})") axes[0, 1].set_xlabel("Budget ratio B/m") axes[0, 1].set_ylabel("Yield1000") axes[0, 1].set_ylim(bottom=0) axes[0, 1].legend(fontsize=7, ncol=2) for method, _, label, color, marker in METHODS[1:]: for d in D_GRID[1:]: points = [find_row(rows, figure_pi, d, ratio, method) for ratio in BUDGET_RATIO_GRID] axes[1, 0].plot( [float(point["budget_ratio"]) for point in points], [float(point["belief_gap_mean"]) for point in points], color=color, marker=marker, label=f"{label}, d={d:g}", ) axes[1, 0].axhline(0.0, color="dimgray", linewidth=0.8) axes[1, 0].set_title(f"Belief-group Recall gap (pi={figure_pi:g})") axes[1, 0].set_xlabel("Budget ratio B/m") axes[1, 0].set_ylabel("Recall(b=1) - Recall(b=0)") axes[1, 0].legend(fontsize=8, ncol=2) for method, _, label, color, marker in METHODS: points = [find_row(rows, figure_pi, 0.0, ratio, method) for ratio in BUDGET_RATIO_GRID] axes[1, 1].plot( [float(point["budget_ratio"]) for point in points], [float(point["recall_mean"]) for point in points], color=color, marker=marker, label=label, ) expected_ratios = np.asarray(BUDGET_RATIO_GRID, dtype=float) axes[1, 1].plot(expected_ratios, expected_ratios, color="0.3", linestyle="--", label="null expectation") axes[1, 1].set_title(f"d=0 negative control (pi={figure_pi:g})") axes[1, 1].set_xlabel("Budget ratio B/m") axes[1, 1].set_ylabel("Mean Recall") axes[1, 1].set_ylim(-0.02, 1.02) axes[1, 1].legend(fontsize=8) figure.suptitle("Synthetic idea-discovery simulation; scenario parameters are not real-world estimates") figure.savefig(FIGURE_PATH, dpi=160, metadata={"Software": "idea_discovery.py"}) plt.close(figure) def validate_rows(rows: list[dict[str, object]]) -> None: expected_rows = len(PI_GRID) * len(D_GRID) * len(BUDGET_RATIO_GRID) * len(METHODS) if len(rows) != expected_rows: raise RuntimeError(f"unexpected row count: {len(rows)} != {expected_rows}") if any(row["m"] != M or row["reps"] < 1000 for row in rows): raise RuntimeError("m or repetition count violates the pre-specified grid") for row in rows: for key, value in row.items(): if isinstance(value, (int, float, np.integer, np.floating)): if not math.isfinite(float(value)): raise RuntimeError(f"non-finite CSV value: {key}") for key in ("recall_mean", "precision_mean", "fdp_single_mean", "belief_selected_share_mean"): if not 0.0 <= float(row[key]) <= 1.0: raise RuntimeError(f"bounded metric outside [0,1]: {key}") if not 0.0 <= float(row["yield1000_mean"]) <= 1000.0: raise RuntimeError("yield1000_mean outside [0,1000]") def negative_control_check(rows: list[dict[str, object]]) -> tuple[float, bool]: deviations = [] for pi in PI_GRID: for ratio in BUDGET_RATIO_GRID: for method, _, _, _, _ in METHODS: row = find_row(rows, pi, 0.0, ratio, method) deviations.append(abs(float(row["recall_mean"]) - ratio)) maximum = max(deviations) return maximum, maximum <= NEGATIVE_CONTROL_TOLERANCE def print_summary(rows: list[dict[str, object]], maximum_negative_deviation: float, negative_control_pass: bool) -> None: print("status=ok") print(f"seed={DEFAULT_SEED}") print(f"m={M}") print(f"reps_per_parameter_unit={REPS}") print(f"parameter_units={len(PI_GRID) * len(D_GRID) * len(BUDGET_RATIO_GRID)}") print(f"csv_rows={len(rows)}") print(f"csv={RESULT_PATH}") print(f"png={FIGURE_PATH}") print(f"d0_negative_control_max_abs_deviation={maximum_negative_deviation:.6f}") print(f"d0_negative_control_tolerance={NEGATIVE_CONTROL_TOLERANCE:.6f}") print(f"d0_negative_control_pass={str(negative_control_pass).lower()}") print("selected_synthetic_rows=pi=0.01,B/m=0.01") for d in D_GRID: row = find_row(rows, 0.01, d, 0.01, "neutral") print( "d={:g},neutral,Recall={:.6f},Precision={:.6f},FDP={:.6f},Yield1000={:.3f},validations_per_tp={:.3f},belief_gap={:.6f}".format( d, float(row["recall_mean"]), float(row["precision_mean"]), float(row["fdp_single_mean"]), float(row["yield1000_mean"]), float(row["validations_per_tp_mean"]), float(row["belief_gap_mean"]), ) ) max_vb_undefined = max(1.0 - float(row["validations_per_tp_defined_fraction"]) for row in rows) max_gap_undefined = max(1.0 - float(row["belief_gap_defined_fraction"]) for row in rows) print(f"max_validations_per_tp_undefined_fraction={max_vb_undefined:.6f}") print(f"max_belief_gap_undefined_fraction={max_gap_undefined:.6f}") print("belief_model=bounded_logistic(lambda_cb*tanh(centered_score))") print("scenario_note=pi,d,B/m,lambda_cb,c_cb are synthetic scenario parameters, not real-world estimates") print("prohibited_extrapolation=specific fields, people, teams, journals, products, budgets, or policies are not represented") def main() -> None: rows: list[dict[str, object]] = [] for pi_index, pi in enumerate(PI_GRID): for d_index, d in enumerate(D_GRID): rows.extend(simulate_unit(pi, d, pi_index, d_index, DEFAULT_SEED)) validate_rows(rows) write_csv(rows) make_figure(rows) if not RESULT_PATH.is_file() or not FIGURE_PATH.is_file(): raise RuntimeError("required output artifact is missing") maximum_negative_deviation, negative_control_pass = negative_control_check(rows) if not negative_control_pass: raise RuntimeError("d=0 negative control exceeded the pre-specified tolerance") print_summary(rows, maximum_negative_deviation, negative_control_pass) if __name__ == "__main__": main()