Make IPF benchmarking fail closed and runnable as one coherent problem

Replace the old order-dependent sequential IPF flow with a validated one-call path. The converter now keeps only closed categorical systems, derives complements only from authored parent totals, reports explicit drop reasons for non-runnable target families, and surfaces incompatible-total failures through structured diagnostics instead of silently chaining margin blocks.

Update the benchmark harness to export retained-authored IPF scoring artifacts, validate the external-IPF input contract, and let L0/GREG opt into scoring on the same retained-authored subset for apples-to-apples comparisons. Refresh the walkthrough notebook, README, example manifest, and add regression coverage for closure logic, runner behavior, export validation, and the new scoring contract.

Author: juaristi22

paper-l0/benchmarking/README.md

@@ -44,7 +44,7 @@ The core workflow is:
 2. export a shared benchmark bundle from a saved calibration package
 3. auto-convert the bundle to IPF inputs when needed
 4. run `L0`, `GREG`, or `IPF`
-5. score all fitted weights against the same shared target matrix
+5. score each method against the target set that matches its benchmark contract
 
 ## Layout
 
@@ -113,30 +113,65 @@ lightweight.
 
 ### IPF inputs
 
-The exporter now auto-generates IPF inputs when the manifest includes `ipf`
-and no external overrides are supplied. It reconstructs an IPF microdata table
-from:
+The exporter auto-generates IPF inputs when the manifest includes `ipf` and no
+external overrides are supplied. It reconstructs an IPF microdata table from:
 
 - the saved calibration package
 - the package metadata's `dataset_path`
 - the package metadata's `db_path`
 - the selected count-like targets and their stratum constraints
 
-The generated `unit_metadata.csv` is currently built for `person_count` and
+The generated `unit_metadata.csv` is built for `person_count` and
 `household_count` targets. It expands cloned households to a person-level table
-when person targets are present, carries a repeated household `unit_index`, and
-adds one derived indicator column per selected target. The generated
-`ipf_target_metadata.csv` then references those indicator columns as numerical
-IPF totals.
+when person targets are present, carries a repeated household `unit_index` so
+per-person weights collapse cleanly back to per-household, and adds one
+string-valued derived category column per declared bucket schema (e.g.
+`age_bracket`, `agi_bracket_district`, `snap_positive`).
+
+The generated `ipf_target_metadata.csv` contains one `categorical_margin` row
+per retained IPF cell after validation. That means:
+
+- authored cells that belong to a closed categorical system are kept
+- binary subset families may gain exactly-derived complement cells when an
+  authored parent total exists on the exact reduced key
+- open subset families are dropped rather than emitted as 1-cell margins
+
+The exporter also writes:
+
+- `ipf_scoring_target_metadata.csv`
+- `ipf_scoring_X_targets_by_units.mtx`
+
+These score IPF on its retained authored targets only. Derived complements are
+recorded for transparency in `ipf_conversion_diagnostics.json`, but they are
+not part of the main benchmark metric set.
+
+When comparing `L0` or `GREG` against that same subset, pass:
+
+```bash
+python paper-l0/benchmarking/benchmark_cli.py run \
+  --method l0 \
+  --run-dir <bundle> \
+  --score-on ipf_retained_authored
+```
 
 External CSVs are still supported through `external_inputs.*` and override the
-automatic conversion path when provided.
+automatic conversion path when provided. The external-IPF contract is strict:
+
+- `external_inputs.ipf_unit_metadata_csv`
+- `external_inputs.ipf_target_metadata_csv`
+- `external_inputs.ipf_scoring_target_metadata_csv`
+- `external_inputs.ipf_scoring_matrix_mtx`
+
+must be provided together. An optional
+`external_inputs.ipf_conversion_diagnostics_json` can also be supplied and will
+be copied through for reporting. External CSVs must also follow the
+`categorical_margin` schema below; the runner rejects `numeric_total` rows.
 
 ### IPF conversion step by step
 
 The IPF conversion is implemented in
-[ipf_conversion.py](/Users/movil1/Desktop/PYTHONJOBS/PolicyEngine/policyengine-us-data/paper-l0/benchmarking/ipf_conversion.py)
-and runs during `benchmark_cli.py export`.
+[ipf_conversion.py](./ipf_conversion.py) and runs during
+`benchmark_cli.py export`.
 
 1. Load the saved calibration package and apply the manifest target filters.
 2. Read `dataset_path`, `db_path`, and `n_clones` from the package metadata.
@@ -152,18 +187,25 @@ and runs during `benchmark_cli.py export`.
    household-clone `unit_index`.
 7. Calculate the needed source variables from the dataset and attach them to
    the IPF unit table.
-8. For each selected target, evaluate its original stratum logic row by row and
-   materialize the result as a derived indicator column such as
-   `ipf_indicator_00000`.
-9. Write `ipf_target_metadata.csv` so each selected target becomes a
-   `numeric_total` IPF constraint over one of those derived indicator columns.
-10. Run `surveysd::ipf` on the generated unit table and target metadata.
-11. Collapse the fitted IPF row weights back to one weight per shared benchmark
-   `unit_index`, so the fitted result can be scored against the same sparse
-   calibration matrix used by `L0` and `GREG`.
-
-This means the benchmark uses one common scoring space even though `IPF`
-requires a richer input representation than `L0` and `GREG`.
+8. Materialize the string-valued derived category columns the margins cover
+   (e.g. `age_bracket`, `snap_positive`) on that unit table.
+9. Group the resolved targets into margin families, validate them against the
+   observed unit-table support, and keep only families that are already closed
+   or can be closed exactly from authored parent totals.
+10. Emit one `categorical_margin` row per retained authored or exactly-derived
+    cell, sharing a `margin_id` within each family.
+11. Write diagnostics (`dropped_targets`, retained-authored counts, derived
+    complements, and any coherence issues) to
+    `inputs/ipf_conversion_diagnostics.json`.
+12. Run `surveysd::ipf` once on the generated unit table and full validated
+    IPF target metadata.
+13. Collapse the fitted IPF row weights back to one weight per shared
+    benchmark `unit_index`, so the fitted result can be scored against the
+    retained-authored sparse target subset used for the IPF benchmark.
+
+This means the benchmark keeps a shared requested target space for the export,
+but an IPF-specific retained-authored scoring space for the actual IPF
+comparison.
 
 ### Why the IPF conversion exists
 
@@ -182,28 +224,21 @@ matrix directly into `surveysd::ipf`.
 
 ### IPF target metadata schema
 
-`ipf_runner.R` supports two target metadata encodings:
-
-- `numeric_total`
-  One row per target with:
-  - `scope`: `person` or `household`
-  - `target_type`: `numeric_total`
-  - `value_column`: unit-data column to calibrate
-  - `variables`: grouping variables used to wrap the numeric total in a one-cell
-    or multi-cell array
-  - `cell`: pipe-separated assignments for the target cell
-  - `target_value`: numeric total
-- `categorical_margin`
-  One row per margin cell with:
-  - `scope`: `person` or `household`
-  - `target_type`: `categorical_margin`
-  - `margin_id`: identifier for a margin table
-  - `variables`: pipe-separated variable names, e.g. `district_id|age_bin`
-  - `cell`: pipe-separated assignments, e.g.
-    `district_id=0601|age_bin=18_24`
-  - `target_value`: numeric target
-
-The automatic conversion path currently emits `numeric_total` rows.
+`ipf_runner.R` accepts one encoding: `categorical_margin`. One row per
+authored margin cell:
+
+- `scope`: `person` or `household`
+- `target_type`: `categorical_margin`
+- `margin_id`: identifier for a margin block. Rows sharing a `margin_id` are
+  grouped into one `surveysd::ipf` constraint (via `xtabs`).
+- `variables`: pipe-separated variable names, e.g.
+  `congressional_district_geoid|age_bracket`
+- `cell`: pipe-separated assignments, e.g.
+  `congressional_district_geoid=0601|age_bracket=0-4`
+- `target_value`: numeric target
+
+Open subset systems are not exported. If a subset family cannot be closed from
+an authored parent total, it is dropped before the R call.
 
 ## Example Commands
 

paper-l0/benchmarking/benchmark_cli.py

@@ -101,10 +101,51 @@ def _run_greg(run_dir: Path):
     return weights_path, elapsed
 
 
+def _collapse_ipf_rows_to_unit_weights(
+    raw_weights: pd.DataFrame, n_units: int
+) -> np.ndarray:
+    """Validate a per-row IPF output and collapse it to a length-`n_units` vector.
+
+    surveysd::ipf with `meanHH = TRUE` guarantees every row that shares a
+    `unit_index` carries the same fitted weight; the spread check keeps that
+    assumption honest.
+    """
+    if "unit_index" not in raw_weights.columns:
+        raise RuntimeError("IPF runner output must include a unit_index column")
+    if raw_weights["unit_index"].isna().any():
+        raise RuntimeError("IPF runner output contains missing unit_index values")
+
+    raw_weights = raw_weights.copy()
+    raw_weights["unit_index"] = raw_weights["unit_index"].astype(np.int64)
+    if (raw_weights["unit_index"] < 0).any() or (
+        raw_weights["unit_index"] >= n_units
+    ).any():
+        raise RuntimeError("IPF runner output contains out-of-range unit_index values")
+
+    per_unit_spread = raw_weights.groupby("unit_index", sort=True)["fitted_weight"].agg(
+        lambda series: float(series.max() - series.min())
+    )
+    if (per_unit_spread > 1e-9).any():
+        raise RuntimeError(
+            "IPF runner produced inconsistent fitted weights within the same unit_index"
+        )
+
+    weights_by_unit = (
+        raw_weights.groupby("unit_index", sort=True)["fitted_weight"]
+        .first()
+        .reindex(np.arange(n_units, dtype=np.int64))
+    )
+    if weights_by_unit.isna().any():
+        raise RuntimeError(
+            "Aggregated IPF weights do not cover the full benchmark unit range"
+        )
+    return weights_by_unit.to_numpy(dtype=np.float64)
+
+
 def _run_ipf(run_dir: Path):
+    """Run one coherent IPF problem in a single `surveysd::ipf` call."""
     inputs = run_dir / "inputs"
     outputs = run_dir / "outputs"
-    temp_csv = outputs / "_ipf_weights.csv"
 
     with open(inputs / "benchmark_manifest.json") as f:
         manifest = json.load(f)
@@ -116,69 +157,97 @@ def _run_ipf(run_dir: Path):
             "IPF run requires inputs/ipf_target_metadata.csv. "
             "Provide external_inputs.ipf_target_metadata_csv in the manifest."
         )
+    unit_metadata_path = inputs / "unit_metadata.csv"
+    if not unit_metadata_path.exists():
+        raise FileNotFoundError("IPF run requires inputs/unit_metadata.csv.")
+
+    full_targets = pd.read_csv(target_metadata_path)
+    if full_targets.empty:
+        raise RuntimeError("IPF target metadata is empty; nothing to run.")
+    unit_metadata = pd.read_csv(unit_metadata_path)
+    if "unit_index" not in unit_metadata.columns:
+        raise RuntimeError("Unit metadata must include a unit_index column for IPF")
+
+    weight_col = str(options.get("weight_col", "base_weight"))
+    household_id_col = str(options.get("household_id_col", "household_id"))
+
+    initial_weights = np.load(inputs / "initial_weights.npy").astype(np.float64)
+    n_units = len(initial_weights)
+    unit_indices = unit_metadata["unit_index"].astype(np.int64).to_numpy()
+    if unit_indices.min() < 0 or unit_indices.max() >= n_units:
+        raise RuntimeError(
+            "Unit metadata unit_index values fall outside the initial weight vector"
+        )
+    temp_csv = outputs / "_ipf_weights.csv"
+    unit_with_weights = unit_metadata.copy()
+    unit_with_weights[weight_col] = initial_weights[unit_indices]
+    temp_unit_csv = outputs / "_ipf_unit_metadata.csv"
+    unit_with_weights.to_csv(temp_unit_csv, index=False)
 
     cmd = [
         "Rscript",
         str(RUNNERS_DIR / "ipf_runner.R"),
-        str(inputs / "unit_metadata.csv"),
+        str(temp_unit_csv),
         str(target_metadata_path),
         str(inputs / "initial_weights.npy"),
         str(temp_csv),
         str(int(options.get("max_iter", 200))),
         str(float(options.get("bound", 4.0))),
         str(float(options.get("epsP", 1e-6))),
         str(float(options.get("epsH", 1e-2))),
-        str(options.get("household_id_col", "household_id")),
-        str(options.get("weight_col", "base_weight")),
+        household_id_col,
+        weight_col,
     ]
     proc, elapsed = _run_subprocess(cmd)
     if proc.returncode != 0:
         raise RuntimeError(f"IPF runner failed with exit code {proc.returncode}")
 
     raw_weights = pd.read_csv(temp_csv)
-    if "unit_index" not in raw_weights.columns:
-        raise RuntimeError("IPF runner output must include a unit_index column")
-    if raw_weights["unit_index"].isna().any():
-        raise RuntimeError("IPF runner output contains missing unit_index values")
-
-    raw_weights["unit_index"] = raw_weights["unit_index"].astype(np.int64)
-    n_units = len(np.load(inputs / "initial_weights.npy"))
-    if (raw_weights["unit_index"] < 0).any() or (
-        raw_weights["unit_index"] >= n_units
-    ).any():
-        raise RuntimeError("IPF runner output contains out-of-range unit_index values")
-
-    per_unit_spread = raw_weights.groupby("unit_index", sort=True)["fitted_weight"].agg(
-        lambda series: float(series.max() - series.min())
-    )
-    inconsistent_units = per_unit_spread[per_unit_spread > 1e-9]
-    if not inconsistent_units.empty:
-        raise RuntimeError(
-            "IPF runner produced inconsistent fitted weights within the same unit_index"
-        )
-
-    weights_by_unit = (
-        raw_weights.groupby("unit_index", sort=True)["fitted_weight"]
-        .first()
-        .reindex(np.arange(n_units, dtype=np.int64))
-    )
-    if weights_by_unit.isna().any():
-        raise RuntimeError(
-            "Aggregated IPF weights do not cover the full benchmark unit range"
-        )
-    weights = weights_by_unit.to_numpy(dtype=np.float64)
+    current_weights = _collapse_ipf_rows_to_unit_weights(raw_weights, n_units)
     weights_path = outputs / "fitted_weights.npy"
-    np.save(weights_path, weights)
+    np.save(weights_path, current_weights)
     temp_csv.unlink(missing_ok=True)
+    temp_unit_csv.unlink(missing_ok=True)
     return weights_path, elapsed
 
 
+def _select_scoring_inputs(
+    run_dir: Path, method: str, score_on: str
+) -> tuple[Path, Path, str]:
+    inputs = run_dir / "inputs"
+    ipf_targets = inputs / "ipf_scoring_target_metadata.csv"
+    ipf_matrix = inputs / "ipf_scoring_X_targets_by_units.mtx"
+    has_ipf_scoring = ipf_targets.exists() and ipf_matrix.exists()
+
+    if score_on == "ipf_retained_authored":
+        if not has_ipf_scoring:
+            raise FileNotFoundError(
+                "Requested score_on=ipf_retained_authored, but "
+                "inputs/ipf_scoring_target_metadata.csv and "
+                "inputs/ipf_scoring_X_targets_by_units.mtx are not both present."
+            )
+        return ipf_targets, ipf_matrix, "ipf_retained_authored"
+
+    if score_on == "auto" and method == "ipf" and has_ipf_scoring:
+        return ipf_targets, ipf_matrix, "ipf_retained_authored"
+    return (
+        inputs / "target_metadata.csv",
+        inputs / "X_targets_by_units.mtx",
+        "shared_requested",
+    )
+
+
 def cmd_run(args):
     run_dir = Path(args.run_dir)
     inputs = run_dir / "inputs"
     outputs = run_dir / "outputs"
     outputs.mkdir(parents=True, exist_ok=True)
-    targets_df = load_targets_csv(inputs / "target_metadata.csv")
+    targets_path, matrix_path, scoring_target_set = _select_scoring_inputs(
+        run_dir,
+        args.method,
+        getattr(args, "score_on", "auto"),
+    )
+    targets_df = load_targets_csv(targets_path)
 
     started = time.time()
     if args.method == "l0":
@@ -195,11 +264,12 @@ def cmd_run(args):
     summary = compute_common_metrics(
         weights=weights,
         targets_df=targets_df,
-        matrix_path=inputs / "X_targets_by_units.mtx",
+        matrix_path=matrix_path,
     )
     summary["method"] = args.method
     summary["run_dir"] = str(run_dir.resolve())
     summary["runtime_seconds"] = elapsed
+    summary["scoring_target_set"] = scoring_target_set
     write_method_summary(summary, outputs / f"{args.method}_summary.json")
     print(json.dumps(summary, indent=2, sort_keys=True))
     return 0
@@ -225,6 +295,16 @@ def build_parser():
     run_parser.add_argument(
         "--run-dir", required=True, help="Exported benchmark bundle directory"
     )
+    run_parser.add_argument(
+        "--score-on",
+        default="auto",
+        choices=["auto", "shared_requested", "ipf_retained_authored"],
+        help=(
+            "Scoring target set. 'auto' uses IPF-retained-authored targets only "
+            "for method=ipf when available; the other methods default to the "
+            "shared requested target set unless explicitly overridden."
+        ),
+    )
     run_parser.set_defaults(func=cmd_run)
 
     return parser