diff --git a/ML/Calibration/calibration_runner.py b/ML/Calibration/calibration_runner.py index ccaf131e..ccdbb485 100644 --- a/ML/Calibration/calibration_runner.py +++ b/ML/Calibration/calibration_runner.py @@ -93,15 +93,19 @@ def _list_jobs_and_exit(cfg, args): # UID splitting # -------------------------------------------------------------------------------- -def _uid_c2st_intervals(job): - """Build the UID splitter and the 'c2st'/'c2st_val' bucket intervals for a job. - - Calibration always evaluates the split held out from BIT training (the - 'c2st_train'/'c2st_val' buckets), matching the convention used for C2ST checks. +def _uid_split_interval(split_cfg: dict, *split_names: str): + """Build the UID splitter and the merged bucket interval for one or more named splits. + + `split_cfg` is a `splitting:` block (`{enabled, uid_fields, seed, n_buckets, scheme}`), + either a job's or `defaults.splitting` directly -- both have the same shape. + Multiple split names are merged into one interval by concatenating their [lo,hi) + bucket ranges in scheme order; this is only a contiguous range when the named + splits are adjacent in `scheme`, which holds for ('c2st_train','c2st_val') and + for any standalone split such as ('final_eval',). """ - split_cfg = job.get("splitting") or {} + split_cfg = split_cfg or {} if not bool(split_cfg.get("enabled", False)): - raise RuntimeError("Calibration requires job.splitting.enabled=True (UID splitting) to avoid data leakage.") + raise RuntimeError("Requires splitting.enabled=True (UID splitting) to avoid data leakage.") uid_fields = split_cfg.get("uid_fields", ["run", "luminosityBlock", "event"]) uid_seed = int(split_cfg.get("seed", 0)) @@ -121,17 +125,27 @@ def _uid_c2st_intervals(job): uid_intervals[k] = (lo, lo + int(sz)) lo += int(sz) - if "c2st_train" not in uid_intervals or "c2st_val" not in uid_intervals: - raise RuntimeError("splitting.scheme must define 'c2st_train' and 'c2st_val'.") - - # assumes they're always in the order c2st_train, c2st_valid - c2st_interval = (uid_intervals["c2st_train"][0], uid_intervals["c2st_val"][1]) + missing = [nm for nm in split_names if nm not in uid_intervals] + if missing: + raise RuntimeError(f"splitting.scheme must define {missing}.") + + merged = (min(uid_intervals[nm][0] for nm in split_names), + max(uid_intervals[nm][1] for nm in split_names)) logger.info("[UID] fields=%s seed=%d n_buckets=%d", uid_fields, uid_seed, uid_n_buckets) logger.info("[UID] scheme intervals: %s", uid_intervals) - logger.info("[UID] BIT train split 'c2st_train'+'c2st_val' -> %s", c2st_interval) + logger.info("[UID] split %s -> %s", split_names, merged) + + return uid_splitter, list(uid_fields), merged - return uid_splitter, list(uid_fields), c2st_interval + +def _uid_c2st_intervals(job): + """Build the UID splitter and the 'c2st_train'+'c2st_val' bucket interval for a job. + + Calibration always evaluates the split held out from BIT training (the + 'c2st_train'/'c2st_val' buckets), matching the convention used for C2ST checks. + """ + return _uid_split_interval(job.get("splitting") or {}, "c2st_train", "c2st_val") # -------------------------------------------------------------------------------- diff --git a/ML/Calibration/eft_calibration.py b/ML/Calibration/eft_calibration.py index f2d7d83c..16d88335 100644 --- a/ML/Calibration/eft_calibration.py +++ b/ML/Calibration/eft_calibration.py @@ -20,7 +20,7 @@ sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), "..", ".."))) sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) -from common.derivative_providers import EFTDerivativeProvider +from common.derivative_providers import build_derivative_provider import calibration_runner as cr @@ -32,7 +32,7 @@ def main(): if not job.get("eft"): raise RuntimeError(f"Job '{job['id']}' has no 'eft' block; use pdf_calibration.py for PDF jobs.") - provider = EFTDerivativeProvider(job["eft"].get("parameters", [])) + provider = build_derivative_provider(job) cr.run_calibration(cfg, job, samples_mod, args, provider) diff --git a/ML/Calibration/pdf_calibration.py b/ML/Calibration/pdf_calibration.py index 06bb6de6..3e396255 100644 --- a/ML/Calibration/pdf_calibration.py +++ b/ML/Calibration/pdf_calibration.py @@ -21,7 +21,7 @@ sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), "..", ".."))) sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) -from common.derivative_providers import PDFDerivativeProvider +from common.derivative_providers import build_derivative_provider import calibration_runner as cr @@ -33,7 +33,7 @@ def main(): if not job.get("pdf"): raise RuntimeError(f"Job '{job['id']}' has no 'pdf' block; use eft_calibration.py for EFT jobs.") - provider = PDFDerivativeProvider(job["pdf"]) + provider = build_derivative_provider(job) cr.run_calibration(cfg, job, samples_mod, args, provider) diff --git a/common/derivative_providers.py b/common/derivative_providers.py index 924a8f8d..52ebcbb1 100644 --- a/common/derivative_providers.py +++ b/common/derivative_providers.py @@ -80,3 +80,12 @@ def __init__(self, parameters): def truth_weight_matrix(self, G, w, observer_names): return self._eft.make_weight_matrix(G, observer_names, w) + + +def build_derivative_provider(job): + """Build the PDF or EFT derivative provider for a BIT job, dispatching on its config block.""" + if job.get("pdf"): + return PDFDerivativeProvider(job["pdf"]) + if job.get("eft"): + return EFTDerivativeProvider(job["eft"].get("parameters", [])) + raise RuntimeError(f"Job '{job['id']}' has neither a 'pdf' nor an 'eft' block.") diff --git a/configs/Eta_unbinned/toys_example.yaml b/configs/Eta_unbinned/toys_example.yaml new file mode 100644 index 00000000..fcb314bb --- /dev/null +++ b/configs/Eta_unbinned/toys_example.yaml @@ -0,0 +1,28 @@ +# Example toy spec for fit/ToyGenerator.py, for configs/Eta_unbinned/Eta_unbinned_2016.yaml +# (the multi-process reference region: SR_2016 classes are TTLep_pow_2016, EtaS_2016, +# EtaP_2016 -- see the "Multi-process handling" section of +# user/ricardo/claude/plans/toy-generation-for-pseudo-experiments.md, verification 6-7). +# +# python fit/ToyGenerator.py configs/Eta_unbinned/Eta_unbinned_2016.yaml \ +# --toySpec configs/Eta_unbinned/toys_example.yaml --toyPoint etaS_up \ +# --seeds 0-499 --outputDir + +source: truth +split: [c2st_train, c2st_val] # null makes cache-mode-vs-truth-mode pull *widths* comparable (verification 6) +throw_nuisances: true +allow_negative_weights: false + +points: + # nominal truth-mode toy: no modifiers -- the surrogate-correct limit. Every + # class gets an implicit ProcessSource at scale_factor=1.0 from its own sample. + - name: nominal + + # EtaS_2016 has a rate_shift POI (a model parameter, so this specific case is + # also expressible in cache mode via hypothesis: {rate_shift_EtaS: ...}); this + # truth-mode version instead probes whether the ICH/ICPH-free unbinned + # surrogates absorb a 20% scale mismodelling correctly. + - name: etaS_up + injection: + SR_2016: + EtaS_2016: + scale_factor: 1.2 diff --git a/configs/unbinned_v6/toys_example.yaml b/configs/unbinned_v6/toys_example.yaml new file mode 100644 index 00000000..2bccc280 --- /dev/null +++ b/configs/unbinned_v6/toys_example.yaml @@ -0,0 +1,53 @@ +# Example toy spec for fit/ToyGenerator.py, for configs/unbinned_v6/unbinned_2016.yaml. +# +# python fit/ToyGenerator.py configs/unbinned_v6/unbinned_2016.yaml \ +# --toySpec configs/unbinned_v6/toys_example.yaml --toyPoint c0_0p5_bit \ +# --seeds 0-499 --outputDir +# +# The spec defines the ensemble (source, split, which points); the CLI only +# selects which point and which seed(s) to generate, and where to write the +# result. See user/ricardo/claude/plans/toy-generation-for-pseudo-experiments.md +# for the full design. + +source: truth # or cache +split: [c2st_train, c2st_val] # truth mode only; null for the whole dataset +throw_nuisances: false +allow_negative_weights: true + +points: + # nominal truth-mode toy: no modifiers, no hypothesis -- the surrogate-correct + # limit, useful as the mean-pull baseline for verification 6. + - name: nominal + + # exact PDF weight injection at c0=0.5, generated from generator-level info. + - name: c0_0p5 + injection: + SR_2016: + TTLep_pow_2016: + coefficients: {c0: 0.5} + + # exact JES NP uncertainty injection via alternative sample + - name: CMS_scale_j_FlavorPureQuark_up + injection: + SR_2016: + TTLep_pow_2016: + sample_name: TTLep_pow_2016_CMS_scale_j_FlavorPureQuark_up + + # exact renormalization scale uncertainty injection via additional weight + - name: mu_ren_up + injection: + SR_2016: + TTLep_pow_2016: + weight_branches: [scale_ren2p0_fac1p0] + + # the same point injected through the trained BIT instead (surrogate route), + # on the *same* truth-mode split -- the control for c0_0p5 (see the plan's + # "two POI injection routes" section). Do not also set 'coefficients' for the + # same POI on the same source -- the two routes are mutually exclusive. + - name: c0_0p5_bit + hypothesis: {c0: 0.5} + +# For a per-process rate mismodelling with no rate nuisance to absorb it (a +# background scaled by e.g. 20%, everything else nominal), see the multi-process +# example in configs/Eta_unbinned/toys_example.yaml, which has more than one class +# per region -- this config's SR_2016 region has only TTLep_pow_2016. diff --git a/configs/unbinned_v7/toys_cache_example.yaml b/configs/unbinned_v7/toys_cache_example.yaml new file mode 100644 index 00000000..e6aad7c0 --- /dev/null +++ b/configs/unbinned_v7/toys_cache_example.yaml @@ -0,0 +1,19 @@ +# Cache-mode toy spec for configs/unbinned_v7/unbinned_2018.yaml, used for the +# first real (non-synthetic) validation of fit/ToyGenerator.py cache-mode +# generation. See user/ricardo/claude/plans/toy-generation-status.md. +# +# python fit/ToyGenerator.py configs/unbinned_v7/unbinned_2018.yaml \ +# --toySpec configs/unbinned_v7/toys_cache_example.yaml --toyPoint nominal \ +# --seeds 0-9 --outputDir + +# TTLep_pow_2018/SingleTop_2018/DrellYan_LO_HTbinned_2018/TTSemi_pow_2018 are all +# NLO-generated samples with ~0.33% negative-weight MC events (routine, not a +# modelling issue), so allow_negative_weights is needed even at the nominal point. +source: cache +throw_nuisances: true +allow_negative_weights: true + +points: + - name: nominal + - name: c0_0p5 + hypothesis: {c0: 0.5} diff --git a/fit/Likelihood.py b/fit/Likelihood.py index 4c050ba1..06e2f1da 100644 --- a/fit/Likelihood.py +++ b/fit/Likelihood.py @@ -1339,18 +1339,26 @@ def _assemble_nuA_groups(self, rid: str, hypothesis) -> Dict[str, list[tuple[dic nuA_per_group[cid] = groups return nuA_per_group - def _compute_T_chunk(self, rid: str, cA_per_class, nuA_per_group, ln_bias_map, rate_shift_map, start: int, stop: int) -> np.ndarray: + def _compute_T_from_columns(self, rid: str, columns_by_class, cA_per_class, nuA_per_group, + ln_bias_map, rate_shift_map, start: int, stop: int) -> np.ndarray: """ - Compute T(x; c, ν) on [start:stop) for a single region rid, summing over classes. + Compute T(x; c, ν) on [start:stop) for a single region rid, summing over classes, + from already-materialized per-class columns. + + `columns_by_class[cid]` may be either an open HDF5 group (the cache) or an + in-memory dict of numpy arrays (an observed/toy dataset's by-class block) -- + both support `col[start:stop]` slicing and `dset_name in col`, which is all + this needs, so one implementation serves Asimov, observation and toy generation. + T_i = Σ_p g_p(x_i) * [ (c⋅R_p)(x_i) * e^{Σ_s ν_B Δ_{p,B}(x_i)} + (e^{...} - 1) ]. """ M = stop - start T = np.zeros(M, dtype=np.float64) for cid in self._class_ids_by_region[rid]: - f = self._h5[(rid, cid)] - g_slice = f['g'][start:stop] # (M,) - R_slice = f['R'][start:stop, :] # (M, nA) + f = columns_by_class[cid] + g_slice = np.asarray(f['g'][start:stop], dtype=np.float64) # (M,) + R_slice = np.asarray(f['R'][start:stop, :], dtype=np.float64) # (M, nA) cA = cA_per_class[cid] if R_slice.shape[1] == 0: @@ -1370,7 +1378,9 @@ def _compute_T_chunk(self, rid: str, cA_per_class, nuA_per_group, ln_bias_map, r expo = np.zeros_like(g_slice) for gm, nuA in nuA_per_group[cid]: dset = gm.get("dset", f"Delta::{gm['id']}") - dA = f[dset][start:stop, :] # (M, nB) + if dset not in f: + raise RuntimeError(f"[N2LL] Missing '{dset}' for {rid}/{cid}.") + dA = np.asarray(f[dset][start:stop, :], dtype=np.float64) # (M, nB) if dA.shape[1] != nuA.shape[0]: raise RuntimeError(f"[N2LL] Δ dim {dA.shape[1]} != ν_A dim {nuA.shape[0]} for {rid}/{cid}/{gm['id']}") expo = expo + (dA @ nuA) # (M,) @@ -1381,6 +1391,91 @@ def _compute_T_chunk(self, rid: str, cA_per_class, nuA_per_group, ln_bias_map, r return T + def _compute_T_chunk(self, rid: str, cA_per_class, nuA_per_group, ln_bias_map, rate_shift_map, start: int, stop: int) -> np.ndarray: + """Compute T(x; c, ν) on [start:stop) for region rid from the cached HDF5 columns.""" + columns_by_class = {cid: self._h5[(rid, cid)] for cid in self._class_ids_by_region[rid]} + return self._compute_T_from_columns(rid, columns_by_class, cA_per_class, nuA_per_group, + ln_bias_map, rate_shift_map, start, stop) + + def _eval_region_surrogates(self, rid: str, X: np.ndarray, feature_names: list[str]) -> dict: + """ + Compute, in memory, the per-class arrays needed to build T(x; c,nu): + - classifier probs g(x) if a classifier is configured (else ones) + - BIT basis R_A(x) + - each PNN group's Δ_B(x) matrix + Returns a by-class dict mirroring the on-disk cache layout: + { cid: {'g': (N,), 'R': (N,nA), 'Delta::': (N,nB), ...}, ... } + + A class with no configured BIT predictor gets `R` of shape (N,0) rather than + raising -- needed for classes with no `POI` block at all, and for `rate_shift` + POI blocks that carry no `predictor`. This is permissive by necessity for + setObservation/toy generation, at a real cost: any *other* caller relying on + this method (e.g. evaluate_ratio) that expects a BIT predictor to always be + present will now silently get c_dot_R == 0 for such a class instead of a + loud failure. + """ + # resolve the region cfg and class list + region = next((R for R in self.regions if R['id'] == rid), None) + if region is None: + raise RuntimeError(f"[_eval_region_surrogates] Unknown region id '{rid}'.") + classes = list(region.get('classes', []) or []) + n_proc = len(classes) + + # input features + mask utility + feat_names = list(feature_names or []) + if not feat_names: + raise RuntimeError("[_eval_region_surrogates] feature_names must be provided.") + X = np.asarray(X, dtype=np.float64, order='C') + if X.ndim != 2: + raise RuntimeError(f"[_eval_region_surrogates] X must be 2D, got shape {X.shape}.") + N = X.shape[0] + + # classifier g(x) + clf = region.get('_classifier_predictor', None) + if clf is None or n_proc <= 1: + g_all = np.ones((N, n_proc), dtype=np.float64) + else: + if not hasattr(clf, "feature_names"): + raise RuntimeError("[_eval_region_surrogates] classifier predictor lacks feature_names.") + mask = self.make_column_mask(feat_names, list(clf.feature_names)) + g_all = _predict_classifier(clf, X[:, mask]) # (N, n_proc) + if g_all.shape[1] != n_proc: + raise RuntimeError(f"[_eval_region_surrogates] classifier outputs {g_all.shape[1]} != {n_proc} classes for region '{rid}'.") + + # per-class outputs + by_class: dict[str, dict] = {} + + for C in classes: + cid = C['id'] + comp = {} + + # g for this process + p_index = class_index(classes, cid) + comp['g'] = np.asarray(g_all[:, p_index], dtype=np.float64, order='C') # (N,) + + # BIT R_A(x); a class with no predictor contributes no POI ratio. + poi_predictor = (C.get('POI') or {}).get('predictor') + if poi_predictor is None: + R_A = np.empty((N, 0), dtype=np.float64) + else: + mask_bit = self.make_column_mask(feat_names, list(poi_predictor.feature_names)) + R_A = predict_bit_ratio(poi_predictor, X[:, mask_bit]) # (N, nA) + comp['R'] = np.asarray(R_A, dtype=np.float64, order='C') + + # PNN Δ groups + for S in C.get('_pnn_systs', []): + sid = S['id'] + pnn = S.get('predictor', None) + if pnn is None: + raise RuntimeError(f"[_eval_region_surrogates] Missing PNN predictor for {rid}/{cid}/{sid}.") + mask_pnn = self.make_column_mask(feat_names, list(pnn.feature_names)) + dA = predict_pnn_deltaA(pnn, X[:, mask_pnn]) # (N, nB) + comp[f"Delta::{sid}"] = np.asarray(dA, dtype=np.float64, order='C') + + by_class[cid] = comp + + return by_class + def setAsimov(self, hypothesis=None) -> None: """ Set an off-nominal Asimov hypothesis (c', ν') and precompute T'(x; c', ν') @@ -1484,27 +1579,11 @@ def setObservation(self, if not self._runtime_prepared: raise RuntimeError("[N2LL.setObservation] Call prepare_runtime() before setting observation.") - # You can’t mix observed-data mode with Asimov in the same evaluation flow. - # We allow switching, but make it explicit and clear. - if getattr(self, "_asimov_hypothesis_set", False) and getattr(self, "_asimov_active", False): - print("[N2LL.setObservation] An Asimov hypothesis had been set; disabling it in favor of observed-data mode.") - if not ignore_weights: print("[N2LL.setObservation] Using weighted sample in setObservation.") - self._asimov_hypothesis_set = False - self._asimov_active = False - self._asimov_hyp = None - - # hypothesis for non-central Asimov - self._asimov_T.clear() - self._binned_asimov_lambda.clear() - - # Reset observation containers - self._obs_unbinned = {} - self._obs_binned_counts = {} - # Flag we’re now in observed-data mode - self._observation_set = True + unbinned_blocks: dict = {} + binned_counts: dict = {} # ------------- UNBINNED OBSERVATION ------------- if unbinned_loaders: @@ -1540,80 +1619,18 @@ def setObservation(self, X_all = np.empty((0, len(getattr(loader, "feature_names", []))), dtype=np.float64) w_all = np.empty((0,), dtype=np.float64) - self._obs_unbinned[rid] = {'X': X_all, 'w': w_all} - # pre-evaluating surrogates on observed events # to avoid evaluating everytime n2ll(hyp) is called - # NB: in binned, evaluation is done directly in n2ll(hyp), # because ICH and ICPH already give the ratio for nominal vs. alternative + by_class = self._eval_region_surrogates(rid, X_all, loader.feature_names) + unbinned_blocks[rid] = {'X': X_all, 'w': w_all, 'by_class': by_class} - self._obs_unbinned[rid]['by_class'] = {} - - # load region info from likelihood object list - # does not protect against two regions with - # the same name - is this allowed ? - # currently, if there's two regions with the same name, - # keeps the last one checked - region_info = {} - for region in self.regions: - if region['id']==rid: - region_info = region - - if region_info == {}: - raise ValueError(f"Did not find information in config corresponding to region {rid}") - + region_info = next(R for R in self.regions if R['id'] == rid) n_classes = len(region_info['classes']) - - class_predictor = region_info.get('_classifier_predictor', None) - - n_classifiers = 0 - if class_predictor is None or n_classes <= 1: - g_obs = np.ones((len(w_all), n_classes), dtype=np.float64) - else: - # in the case where class_predictor uses - # subset of the features in dataloader - class_predictor_column_mask = self.make_column_mask(loader.feature_names, class_predictor.feature_names) - - g_obs = _predict_classifier(class_predictor, X_all[:, class_predictor_column_mask]) # (n_events, n_classes) - if g_obs.shape[1] != n_classes: - raise RuntimeError(f"[N2LL] Classifier outputs {g_obs.shape[1]} != {n_classes} classes in region '{rid}'") - n_classifiers+=1 - - n_poi_predictors = 0 - n_syst_predictors = 0 - for class_info in region_info['classes']: - cid = class_info['id'] - i_class = class_index(region_info['classes'], cid) - - self._obs_unbinned[rid]['by_class'][cid] = {} - - # classifier output per classfrom global classifier - self._obs_unbinned[rid]['by_class'][cid]['g'] = g_obs[:,i_class] - - # POI predictor for each class - poi_predictor = (class_info.get('POI') or {}).get('predictor') - - if poi_predictor is None: - R_A = np.empty((len(w_all), 0), dtype=np.float64) - else: - poi_predictor_column_mask = self.make_column_mask(loader.feature_names, poi_predictor.feature_names) - R_A = predict_bit_ratio(poi_predictor, X_all[:,poi_predictor_column_mask]) - n_poi_predictors +=1 - - self._obs_unbinned[rid]['by_class'][cid]['R'] = R_A - - # syst predictors for each class - for syst_info in class_info['_pnn_systs']: - - syst_column_mask = self.make_column_mask(loader.feature_names, syst_info['predictor'].feature_names) - pnn = syst_info['predictor'] - syst_id = syst_info['id'] - dA = predict_pnn_deltaA(pnn, X_all[:, syst_column_mask]) # (N_events, nB) - dset = f'Delta::{syst_id}' - - self._obs_unbinned[rid]['by_class'][cid][dset] = dA - n_syst_predictors += 1 + n_classifiers = 1 if (region_info.get('_classifier_predictor') is not None and n_classes > 1) else 0 + n_poi_predictors = sum(1 for C in region_info['classes'] if (C.get('POI') or {}).get('predictor') is not None) + n_syst_predictors = sum(len(C.get('_pnn_systs', []) or []) for C in region_info['classes']) print(f"[setObservation] Unbinned region '{rid}': loaded {X_all.shape[0]:,} events " f"({'unit weights' if ignore_weights else 'with weights'}).") @@ -1672,10 +1689,87 @@ def setObservation(self, counts2d += H.astype(np.float64) flat_counts = counts if len(edges) == 1 else counts2d.reshape(-1) - self._obs_binned_counts[rid] = flat_counts + binned_counts[rid] = flat_counts print(f"[setObservation] Binned region '{rid}': filled {flat_counts.size} bins " f"({'unit weights' if ignore_weights else 'with weights'}).") + self.setObservationArrays(unbinned_blocks=unbinned_blocks or None, binned_counts=binned_counts or None) + + def setObservationArrays(self, unbinned_blocks: dict | None = None, binned_counts: dict | None = None) -> None: + """ + Register an observed/toy dataset from already-evaluated arrays (no loaders). + + Parameters + ---------- + unbinned_blocks : dict or None + Mapping {region_id -> {'X': features (N,d) or None, 'w': weights (N,), + 'by_class': {cid: {'g': (N,), 'R': (N,nA), 'Delta::': (N,nB), ...}}}}. + 'X' may be None (e.g. cache-mode toys, which have no materialized X). + binned_counts : dict or None + Mapping {region_id -> counts_flat (Nflat,)}. + + Effects mirror setObservation: sets `self._observation_set = True`, disables + any previously set Asimov bias, and populates `self._obs_unbinned` / + `self._obs_binned_counts`. `setObservation` ends by calling this, so this is + the single place that owns the observed/Asimov mode switch. + """ + if not self._runtime_prepared: + raise RuntimeError("[N2LL.setObservationArrays] Call prepare_runtime() before setting observation.") + + # You can’t mix observed-data mode with Asimov in the same evaluation flow. + # We allow switching, but make it explicit and clear. + if getattr(self, "_asimov_hypothesis_set", False) and getattr(self, "_asimov_active", False): + print("[N2LL.setObservationArrays] An Asimov hypothesis had been set; disabling it in favor of observed-data mode.") + + self._asimov_hypothesis_set = False + self._asimov_active = False + self._asimov_hyp = None + + # hypothesis for non-central Asimov + self._asimov_T.clear() + self._binned_asimov_lambda.clear() + + # Reset observation containers + self._obs_unbinned = {} + self._obs_binned_counts = {} + + # Flag we’re now in observed-data mode + self._observation_set = True + + known_region_ids = {R['id'] for R in self.regions} + if unbinned_blocks: + for rid, block in unbinned_blocks.items(): + if rid not in known_region_ids: + raise RuntimeError(f"[setObservationArrays:unbinned] Unknown region id '{rid}'.") + w = np.asarray(block['w'], dtype=np.float64) + for cid, comp in block.get('by_class', {}).items(): + for col_name, col in comp.items(): + if len(col) != len(w): + raise RuntimeError( + f"[setObservationArrays:unbinned:{rid}/{cid}] " + f"len({col_name})={len(col)} != len(w)={len(w)}.") + self._obs_unbinned[rid] = block + + if binned_counts: + for rid, counts in binned_counts.items(): + if rid not in self._binned_unroll: + raise RuntimeError(f"[setObservationArrays:binned] Region '{rid}' has no binned definition in current likelihood.") + self._obs_binned_counts[rid] = np.asarray(counts, dtype=np.float64) + + def setToy(self, toy: dict, hypothesis) -> None: + """ + Register a toy (from fit/ToyGenerator.py's generate_toy/load_toy) as the + observation, and apply its thrown constraint centres to `hypothesis`. + """ + unbinned_blocks = { + rid: {'X': block.get('X'), 'w': block['w'], 'by_class': block['by_class']} + for rid, block in toy.get('unbinned_blocks', {}).items() + } + self.setObservationArrays(unbinned_blocks=unbinned_blocks or None, + binned_counts=toy.get('binned_counts') or None) + base = hypothesis._base if hasattr(hypothesis, '_base') else hypothesis + base.set_constraint_centers(toy.get('constraint_centers', {})) + def __call__(self, hypothesis) -> float: """ Evaluate -2 log L for either: @@ -1711,7 +1805,6 @@ def __call__(self, hypothesis) -> float: byc = block['by_class'] W = np.asarray(block['w'], dtype=np.float64) N = len(W) - T = np.zeros(N, dtype=np.float64) # current hypothesis A-basis and ν_A groups cA_per_class = self._assemble_cA_per_class(rid, hypothesis._base) @@ -1745,32 +1838,7 @@ def __call__(self, hypothesis) -> float: # 1: get T evaluated on the observed events, summing over predictions of surrogates for each class # 2: make weighted sum of log1p(T), for observed weights can be one or not # depends on whether setObservation was called with ignoreWeights=True - for cid, comp in byc.items(): - g_slice = np.asarray(comp['g'], dtype=np.float64) # (N,) - R_slice = np.asarray(comp['R'], dtype=np.float64) # (N, nA) - cA = cA_per_class[cid] # (nA,) - if R_slice.shape[1] != cA.shape[0]: - raise RuntimeError(f"[N2LL:obs:unbinned] BIT dim {R_slice.shape[1]} != |A| {cA.shape[0]} for {rid}/{cid}") - c_dot_R = R_slice @ cA # (N,) - - # Adding rate shift - rs = rate_shift.get(cid, 0.0) - if rs != 0.0: - c_dot_R = c_dot_R + rs - - expo = np.zeros_like(g_slice) - for gm, nuA in nuA_per_group[cid]: - dset = gm.get("dset", f"Delta::{gm['id']}") - if dset not in comp: - raise RuntimeError(f"[N2LL:obs:unbinned] Missing '{dset}' for {rid}/{cid}.") - dA = np.asarray(comp[dset], dtype=np.float64) # (N, nB) - if dA.shape[1] != nuA.shape[0]: - raise RuntimeError(f"[N2LL:obs:unbinned] Δ dim {dA.shape[1]} != ν_A dim {nuA.shape[0]} for {rid}/{cid}/{gm['id']}") - expo = expo + (dA @ nuA) - - # include per-class lnN bias additively in exponent - expo += ln_bias[cid] # (M,) - T += g_slice * (c_dot_R * np.exp(expo) + np.expm1(expo)) + T = self._compute_T_from_columns(rid, byc, cA_per_class, nuA_per_group, ln_bias, rate_shift, 0, N) total_unbinned += np.dot(W,np.log1p(T)) @@ -2353,6 +2421,9 @@ def plot_correlation_root(out_dir, base, rotated, names, corr, suffix=""): p.add_argument("--no_syst", action="store_true", help="Disable all nuisances (freeze to 0).") p.add_argument("--syst_only", action="store_true", help="Disable all POIs (freeze to 0).") p.add_argument("--data", action="store_true", help="Fits to data defined in config.") + p.add_argument("--toyFile", default=None, + help="Path to a toy HDF5 file generated by fit/ToyGenerator.py; fits to it instead of " + "data/Asimov. Generation-only: no spec parsing, seeding or point selection happens here.") p.add_argument("--asimov", nargs="+", default=None, metavar=("PAR", "VAL"), help="Set an off-nominal Asimov hypothesis via pairs: --asimov par1 val1 par2 val2 ...") p.add_argument("--shuffle", nargs="+", default=None, help="Shuffle these features") @@ -2462,6 +2533,12 @@ def plot_correlation_root(out_dir, base, rotated, names, corr, suffix=""): if args.data: suffix += "_data" print("Fitting to data!") + if args.toyFile: + with h5py.File(args.toyFile, "r") as _toy_meta_f: + _toy_point = str(_toy_meta_f["meta"].attrs.get("point", "")) or "toy" + _toy_seed = int(_toy_meta_f["meta"].attrs["seed"]) + suffix += f"_{_toy_point}_toy{_toy_seed}" + print(f"Fitting to toy '{_toy_point}' seed {_toy_seed} from {args.toyFile}") os.makedirs(user.output_directory, exist_ok=True) out_path = os.path.join(user.output_directory, f"{base}_{version}{suffix}_fit.json") @@ -2599,8 +2676,16 @@ def plot_correlation_root(out_dir, base, rotated, names, corr, suffix=""): if (region['data'].get('ignore_weights', True) is False) and ignore_weights: ignore_weights = False + # ---- toy: generated separately by fit/ToyGenerator.py, loaded (never generated) here ---- + if args.toyFile: + if args.data or args.asimov is not None: + raise RuntimeError("--toyFile cannot be combined with --data or --asimov.") + from fit.ToyGenerator import load_toy + toy = load_toy(args.toyFile, n2ll) + n2ll.setToy(toy, hyp_for_fit) + # data - if args.data: + elif args.data: if (not unbinned_dataloaders) and (not binned_dataloaders): raise ValueError("You asked for a data fit, but did not define any dataset in your config. Exiting!") n2ll.setObservation(unbinned_dataloaders, binned_dataloaders, ignore_weights=ignore_weights) diff --git a/fit/Modeling.py b/fit/Modeling.py index ed6e134e..6f4fd351 100644 --- a/fit/Modeling.py +++ b/fit/Modeling.py @@ -36,12 +36,13 @@ class ModelParameter: Lightweight model parameter representation with convenient helpers. """ def __init__(self, name, val=0.0, *, isPOI=False, - isFrozen=False, isPenalized=False): - self.name = str(name) - self.val = _coerce_parameter_value(val) - self.isPOI = bool(isPOI) - self.isFrozen = bool(isFrozen) - self.isPenalized = bool(isPenalized) + isFrozen=False, isPenalized=False, constraint_center=0.0): + self.name = str(name) + self.val = _coerce_parameter_value(val) + self.isPOI = bool(isPOI) + self.isFrozen = bool(isFrozen) + self.isPenalized = bool(isPenalized) + self.constraint_center = float(constraint_center) def __repr__(self): tags = [] @@ -49,7 +50,7 @@ def __repr__(self): else: tags.append("Nuis.") if self.isFrozen: tags.append("frozen") if self.isPenalized: tags.append("pen.") - return f"<{self.name}({','.join(tags)})={self.val:.6e}>" + return f"<{self.name}({','.join(tags)})={float(getval(self.val)):.6e}>" def __str__(self): return self.__repr__().lstrip('<').rstrip('>') @@ -203,8 +204,18 @@ def penalized(self): # ---------- penalty ---------- def penalty(self): - """Compute the penalty (sum v**2) from all penalized nuisances.""" - return sum(p.val**2 for p in self.parameters if p.isPenalized) + """Compute the penalty (sum (v - center)**2) from all penalized nuisances.""" + return sum((p.val - p.constraint_center)**2 for p in self.parameters if p.isPenalized) + + def set_constraint_centers(self, mapping: dict) -> None: + """Set constraint centers by parameter name (e.g. a toy's thrown nu_obs).""" + for name, center in mapping.items(): + if name not in self: + raise KeyError(f"[Hypothesis.set_constraint_centers] Unknown parameter '{name}'.") + p = self[name] + if not p.isPenalized: + raise RuntimeError(f"[Hypothesis.set_constraint_centers] Parameter '{name}' is not penalized.") + p.constraint_center = float(center) # ---------- mutators ---------- def modify(self, **kwargs): @@ -338,7 +349,8 @@ def __init__(self, base: Hypothesis, json_filename: str, name: str | None = None ModelParameter(name=nm, val=float(getattr(self._base, nm).val), isPOI=False, - isPenalized=bool(getattr(self._base, nm).isPenalized)) + isPenalized=bool(getattr(self._base, nm).isPenalized), + constraint_center=float(getattr(self._base, nm).constraint_center)) for nm in self._nuis_names ] all_params = d_params + nuis_params diff --git a/fit/N2LLExtensions.py b/fit/N2LLExtensions.py index 34d09073..bba2cc69 100644 --- a/fit/N2LLExtensions.py +++ b/fit/N2LLExtensions.py @@ -467,86 +467,8 @@ def _assemble_nuA_groups_from_cfg(self, rid: str, hypothesis) -> dict[str, list[ out[cid] = groups return out - def _eval_region_surrogates(self, rid: str, X: np.ndarray, feature_names: list[str]) -> dict: - """ - Compute, in memory, the per-class arrays needed to build T(x; c,nu): - - classifier probs g(x) if a classifier is configured (else ones) - - BIT basis R_A(x) - - each PNN group's Δ_B(x) matrix - Returns a by-class dict mirroring the shape used elsewhere: - { cid: {'g': (N,), 'R': (N,nA), 'Delta::': (N,nB), ...}, ... } - """ - import numpy as np - - # resolve the region cfg and class list - region = next((R for R in self.regions if R['id'] == rid), None) - if region is None: - raise RuntimeError(f"[evaluate_ratio] Unknown region id '{rid}'.") - classes = list(region.get('classes', []) or []) - n_proc = len(classes) - - # input features + mask utility - feat_names = list(feature_names or []) - if not feat_names: - raise RuntimeError("[evaluate_ratio] feature_names must be provided.") - X = np.asarray(X, dtype=np.float64, order='C') - if X.ndim != 2: - raise RuntimeError(f"[evaluate_ratio] X must be 2D, got shape {X.shape}.") - N = X.shape[0] - - # classifier g(x) - g_all = None - clf = region.get('_classifier_predictor', None) - if clf is None or n_proc <= 1: - g_all = np.ones((N, n_proc), dtype=np.float64) - else: - if not hasattr(clf, "feature_names"): - raise RuntimeError("[evaluate_ratio] classifier predictor lacks feature_names.") - mask = self.make_column_mask(feat_names, list(clf.feature_names)) - g_all = _predict_classifier(clf, X[:, mask]) # (N, n_proc) - if g_all.shape[1] != n_proc: - raise RuntimeError(f"[evaluate_ratio] classifier outputs {g_all.shape[1]} != {n_proc} classes for region '{rid}'.") - - # per-class outputs - by_class: dict[str, dict] = {} - - for C in classes: - cid = C['id'] - comp = {} - - # g for this process - p_index = class_index(classes, cid) - comp['g'] = np.asarray(g_all[:, p_index], dtype=np.float64, order='C') # (N,) - - # BIT R_A(x) - poi = (C.get('POI') or {}) - bit = poi.get('predictor', None) - if bit is None: - raise RuntimeError(f"[evaluate_ratio] Missing BIT predictor for {rid}/{cid}.") - if not hasattr(bit, "feature_names"): - raise RuntimeError(f"[evaluate_ratio] BIT predictor lacks feature_names for {rid}/{cid}.") - mask_bit = self.make_column_mask(feat_names, list(bit.feature_names)) - R_A = predict_bit_ratio(bit, X[:, mask_bit]) # (N, nA) - comp['R'] = np.asarray(R_A, dtype=np.float64, order='C') - - # PNN Δ groups - for S in C.get('_pnn_systs', []): - sid = S['id'] - pnn = S.get('predictor', None) - if pnn is None: - raise RuntimeError(f"[evaluate_ratio] Missing PNN predictor for {rid}/{cid}/{sid}.") - # Feature selection for PNN: prefer predictor.feature_names if present - if hasattr(pnn, "feature_names") and pnn.feature_names: - mask_pnn = self.make_column_mask(feat_names, list(pnn.feature_names)) - dA = predict_pnn_deltaA(pnn, X[:, mask_pnn]) # (N, nB) - else: - # Fallback: assume pnn.deltaA expects the provided X ordering. - dA = predict_pnn_deltaA(pnn, X) - comp[f"Delta::{sid}"] = np.asarray(dA, dtype=np.float64, order='C') - - by_class[cid] = comp - - return by_class + # _eval_region_surrogates moved to the base N2LL class (fit/Likelihood.py), + # since setObservation/toy generation need it too. def evaluate_ratio(self, rid: str, X: np.ndarray, feature_names: list[str], hypothesis, *, cached: bool = True, return_T: bool = False, diff --git a/fit/ToyGenerator.py b/fit/ToyGenerator.py new file mode 100644 index 00000000..8b9b5590 --- /dev/null +++ b/fit/ToyGenerator.py @@ -0,0 +1,909 @@ +"""Toy dataset generation for pseudo-experiments. + +Two independent axes: `source` (this module's "cache" vs "truth", see below) and +`mode` (N2LL's existing Asimov/observation distinction in fit/Likelihood.py, +untouched here). A toy is a plain dict, produced by `generate_toy`, persisted with +`save_toy`/`load_toy`, and applied to an N2LL via `N2LL.setToy`. + +- `source="cache"`: the model (trained surrogates) is taken as truth. A toy is a + Poisson draw on the cached MC's own intensity `w0 * (1 + T(hypothesis))`. Answers + "given the surrogates are correct, are my intervals calibrated?" +- `source="truth"`: pseudo-data is generated from an exact reweighting, a + per-process scale factor, or an alternate sample (`ProcessSource`), while the fit + still uses the trained surrogates. Probes surrogate mismodelling as a bias in the + fitted coefficients. + +Known caveat (documented, not fixed): toys are drawn from the same cached/streamed +MC that supplies the expected-yield term, so toy and template fluctuations are +correlated. Sub-percent effect on interval widths at current MC statistics; see +user/ricardo/claude/plans/toy-generation-design-decisions.md for why this is not +corrected. +""" +from __future__ import annotations + +import os +import sys +import json +import hashlib +import logging +import importlib +from dataclasses import dataclass +from typing import Callable, Optional + +import numpy as np +import h5py + +# project root (this file lives in fit/) + ML/Calibration (for the shared UID-split helper) +_REPO_ROOT = os.path.abspath(os.path.join(os.path.dirname(__file__), "..")) +sys.path.insert(0, _REPO_ROOT) +sys.path.insert(0, os.path.join(_REPO_ROOT, "ML", "Calibration")) + +from fit.Likelihood import N2LL, expand_pois_linear_quadratic, build_hypothesis_from_likelihood +from fit.Modeling import Hypothesis +import calibration_runner as cr # _uid_split_interval + +# common.derivative_providers is imported lazily inside _materialize_truth_weights +# (see below): it pulls in data/samples_eft.py, which eagerly constructs RDataLoaders +# for every declared EFT sample at import time, so importing it at module level would +# make cache-mode generation depend on the EFT sample data being reachable even though +# cache mode never touches it. + +logger = logging.getLogger(__name__) + + +# ============================================================================ +# Multi-process truth sources +# ============================================================================ + +@dataclass +class ProcessSource: + """One process's contribution to a truth-mode toy in one region. + + w_truth = w_coefficients * scale_factor * prod(weight_branches) * weight_function(X) + reducing to the nominal event weight when `coefficients` is None. + """ + class_id: str + sample_name: Optional[str] = None + coefficients: Optional[dict] = None + scale_factor: float = 1.0 + weight_branches: Optional[list] = None + weight_function: Optional[Callable] = None + + +# ============================================================================ +# RNG: stable per-(region[,class]) spawning +# ============================================================================ + +def _stable_str_to_int(s: str) -> int: + """Process-stable string hash (Python's built-in hash() is salted per-process).""" + return int.from_bytes(hashlib.blake2b(s.encode("utf-8"), digest_size=8).digest(), "little") + + +def _spawn_rng(seed: int, *parts: str) -> np.random.Generator: + """RNG for a single (seed, region_id[, class_id]) tuple, via a stable hash of the + parts -- not sequential spawning off one Generator -- so adding a region or + process does not perturb another's draws.""" + entropy = [int(seed) & 0xFFFFFFFF] + [_stable_str_to_int(p) & 0xFFFFFFFF for p in parts] + return np.random.default_rng(np.random.SeedSequence(entropy)) + + +# ============================================================================ +# constraint centers +# ============================================================================ + +def throw_constraint_centers(hypothesis, rng: np.random.Generator) -> dict: + """Throw nu_obs ~ Normal(nu_true, 1) for every penalized nuisance; {name: nu_obs}.""" + return {p.name: float(rng.normal(p.val, 1.0)) for p in hypothesis.parameters if p.isPenalized} + + +# ============================================================================ +# cache-mode unbinned +# ============================================================================ + +def _compute_cache_intensity(n2ll: N2LL, region_id: str, hypothesis) -> np.ndarray: + """w0_i * (1 + T_i(hypothesis)) over the whole cached MC for one region -- the + model's expected (Asimov) per-event intensity, unclipped. Shared by + `generate_unbinned_toy_from_cache` (which draws Poisson toys from it) and the + diagnostic Asimov reference histogram in plot/toys (which sums it as-is).""" + N = n2ll._N_region.get(region_id, 0) + if N == 0: + return np.empty(0, dtype=np.float64) + + cA_per_class = n2ll._assemble_cA_per_class(region_id, hypothesis) + nuA_per_group = n2ll._assemble_nuA_groups(region_id, hypothesis) + nu_vals = {p.name: p.val for p in getattr(hypothesis, "parameters", []) if not p.isPOI} + ln_bias = { + cid: sum(a * nu_vals.get(nm, 0.0) for nm, a in n2ll._lnN_by_class.get((region_id, cid), [])) + for cid in n2ll._class_ids_by_region[region_id] + } + rate_shift = n2ll._assemble_rate_shift_per_class(region_id, hypothesis) + + class_ids = n2ll._class_ids_by_region[region_id] + w0 = np.asarray(n2ll._h5[(region_id, class_ids[0])]["w0"], dtype=np.float64) + T = n2ll._compute_T_chunk(region_id, cA_per_class, nuA_per_group, ln_bias, rate_shift, 0, N) + return w0 * (1.0 + T) + + +def generate_unbinned_toy_from_cache(n2ll: N2LL, region_id: str, hypothesis, rng: np.random.Generator, + *, allow_negative_weights: bool = False) -> dict: + """Cache-mode unbinned toy: n_i ~ Poisson(w0_i * (1 + T_i(hypothesis))) over the + whole cached MC (no UID split -- see the plan's "UID splitting" section for why + cache mode uses everything). + + Negative intensity has two distinct causes, both real: at hypothesis=0, T is + identically 0 for every event, so a negative w0_i * (1+T_i) there can only come + from a negative-weight MC event (routine for NLO generators like aMC@NLO/POWHEG, + not a modelling problem); away from hypothesis=0 it can also come from T < -1 + (the model extrapolating outside its valid range). Both are reported and treated + the same way as truth mode's negative-weight handling: crash by default with the + offending fraction, clip to zero under allow_negative_weights=True. + + Returns row indices and multiplicities for n_i > 0 (not sliced columns -- the + surrogate columns are verbatim cache slices, rehydrated by `load_toy`). + """ + N = n2ll._N_region.get(region_id, 0) + if N == 0: + return {"indices": np.empty(0, dtype=np.int64), "n": np.empty(0, dtype=np.float64)} + + lam = _compute_cache_intensity(n2ll, region_id, hypothesis) + + neg = lam < 0.0 + n_neg = int(np.sum(neg)) + if n_neg > 0: + neg_weight_sum = float(np.sum(lam[neg])) + if not allow_negative_weights: + raise RuntimeError( + f"[toy:cache:{region_id}] {n_neg}/{N} events ({n_neg / N:.4%}) have negative intensity " + f"w0*(1+T) (summed negative intensity {neg_weight_sum:.4g}); pass " + f"allow_negative_weights=True to clip them to zero." + ) + logger.info("[toy:cache:%s] clipping %d/%d negative-intensity events (summed %.4g) to zero.", + region_id, n_neg, N, neg_weight_sum) + lam = np.where(neg, 0.0, lam) + + n_i = rng.poisson(lam) + keep = n_i > 0 + return {"indices": np.nonzero(keep)[0].astype(np.int64), "n": n_i[keep].astype(np.float64)} + + +def _rehydrate_cache_by_class(n2ll: N2LL, region_id: str, indices: np.ndarray) -> dict: + """Rebuild by-class surrogate columns for a cache-mode toy by indexing the live + cache. Structural binding: an out-of-range index (e.g. a retrained/shrunk cache) + raises IndexError rather than silently misaligning columns.""" + class_ids = n2ll._class_ids_by_region[region_id] + by_class = {} + for cid in class_ids: + cols = n2ll._h5[(region_id, cid)] + by_class[cid] = {k: np.asarray(cols[k])[indices] for k in cols if k != "w0"} + return by_class + + +def materialize_cache_region_diagnostics(n2ll: N2LL, region_id: str, hypothesis, feature_names: list, + plot_opts: dict, extra_indices: Optional[np.ndarray] = None): + """Single re-stream pass over a region's Asimov samples, for diagnostic plotting + (see plot/toys/toy_diagnostic_plots.py). Cache-mode toys store only (indices, + surrogate columns), not X, to avoid duplicating the surrogate cache, so this is + the opt-in, costly path back to raw kinematics. + + Returns (asimov_hist, extra_X): + - asimov_hist: {feature_name: weighted histogram of w0*(1+T) at `hypothesis` + over the WHOLE cached MC, binned per `plot_opts`} -- what toys generated at + this hypothesis should fluctuate around, always available since the + surrogate cache is always built regardless of toy source. Negative + intensities are clipped to zero for display (a diagnostic reference isn't + the place to enforce generation's crash-by-default policy), with a logged + count/fraction. + - extra_X: raw feature matrix at `extra_indices` (e.g. the union of + cache-mode toys' kept row indices), or None if `extra_indices` is None. + Piggy-backs on the same pass since computing `asimov_hist` already visits + every cached event once. + + `extra_indices`, if given, must be sorted ascending (true of both + `generate_unbinned_toy_from_cache`'s and `np.unique`'s output). + """ + region = _find_region(n2ll, region_id) + + lam = _compute_cache_intensity(n2ll, region_id, hypothesis) + neg = lam < 0.0 + n_neg = int(np.sum(neg)) + if n_neg > 0: + logger.info( + "[toy:plot:%s] clipping %d/%d negative-intensity events (summed %.4g) to zero " + "for the Asimov reference histogram.", region_id, n_neg, len(lam), float(np.sum(lam[neg])), + ) + lam = np.where(neg, 0.0, lam) + + edges = { + feat: np.linspace(*plot_opts[feat]["binning"][1:], int(plot_opts[feat]["binning"][0]) + 1) + for feat in feature_names if feat in plot_opts + } + asimov_hist = {feat: np.zeros(len(e) - 1, dtype=np.float64) for feat, e in edges.items()} + + extra_indices = None if extra_indices is None else np.asarray(extra_indices, dtype=np.int64) + if extra_indices is not None and not np.all(np.diff(extra_indices) >= 0): + raise ValueError(f"[toy:plot:{region_id}] extra_indices must be sorted ascending.") + extra_X = np.empty((len(extra_indices), len(feature_names)), dtype=np.float64) if extra_indices is not None else None + + col_positions = None + row_ptr = 0 + global_offset = 0 + for feat_names, X, _w0 in n2ll._iter_asimov_batches(region): + if col_positions is None: + pos = {f: i for i, f in enumerate(feat_names)} + missing = [f for f in feature_names if f not in pos] + if missing: + raise KeyError(f"[toy:plot:{region_id}] Feature(s) {missing} not in region samples.") + col_positions = [pos[f] for f in feature_names] + + batch_hi = global_offset + len(X) + batch_lam = lam[global_offset:batch_hi] + for feat, e in edges.items(): + asimov_hist[feat] += np.histogram(X[:, pos[feat]], bins=e, weights=batch_lam)[0] + + if extra_indices is not None: + while row_ptr < len(extra_indices) and extra_indices[row_ptr] < batch_hi: + extra_X[row_ptr] = X[extra_indices[row_ptr] - global_offset, col_positions] + row_ptr += 1 + global_offset = batch_hi + + if extra_indices is not None and row_ptr < len(extra_indices): + raise RuntimeError( + f"[toy:plot:{region_id}] {len(extra_indices) - row_ptr} indices exceeded the re-streamed " + f"event count ({global_offset}); does the toy's cache match this config's samples?" + ) + return asimov_hist, extra_X + + +# ============================================================================ +# truth-mode: resolving sources, feature union, per-source truth weights +# ============================================================================ + +def _find_region(n2ll: N2LL, region_id: str) -> dict: + for region in list(n2ll.regions) + list(n2ll.binned): + if region["id"] == region_id: + return region + raise RuntimeError(f"[toy] Unknown region '{region_id}'.") + + +def _find_class(region: dict, class_id: str) -> dict: + for C in region.get("classes", []) or []: + if C["id"] == class_id: + return C + raise RuntimeError(f"[toy] Unknown class '{class_id}' in region '{region['id']}'.") + + +def _resolve_process_sources(classes: list, explicit_sources) -> dict: + """Merge explicit ProcessSources with implicit ones (scale_factor=1 from the + class's own sample) for every class with no explicit source.""" + by_cid: dict = {} + for s in explicit_sources or []: + if s.class_id in by_cid: + raise RuntimeError(f"[toy] Duplicate ProcessSource for class '{s.class_id}'.") + by_cid[s.class_id] = s + known_cids = {C["id"] for C in classes} + unknown = set(by_cid) - known_cids + if unknown: + raise RuntimeError(f"[toy] Unknown class id(s) in truth sources: {sorted(unknown)}") + for C in classes: + cid = C["id"] + if cid not in by_cid: + by_cid[cid] = ProcessSource(class_id=cid) + return by_cid + + +def _region_feature_union(region: dict) -> list: + """Union of every feature list the region's surrogates consume (classifier + + every class's BIT + every class's PNNs), built from the loaded predictors' + feature_names -- not from any job['features'] list.""" + predictors = [] + clf = region.get("_classifier_predictor") + if clf is not None: + predictors.append(clf) + for C in region.get("classes", []) or []: + poi_pred = (C.get("POI") or {}).get("predictor") + if poi_pred is not None: + predictors.append(poi_pred) + for S in C.get("_pnn_systs", []) or []: + pred = S.get("predictor") + if pred is not None: + predictors.append(pred) + + names, seen = [], set() + for pred in predictors: + for f in getattr(pred, "feature_names", []) or []: + if f not in seen: + seen.add(f) + names.append(f) + return names + + +def _class_bit_job(n2ll: N2LL, region_id: str, class_id: str) -> dict: + region = _find_region(n2ll, region_id) + C = _find_class(region, class_id) + job_id = (C.get("POI") or {}).get("job") + if not job_id: + raise RuntimeError(f"[toy] Class '{region_id}/{class_id}' has no configured BIT job (POI.job).") + jobs_by_id = getattr(n2ll, "_toy_jobs_by_id", None) + if not jobs_by_id: + raise RuntimeError( + "[toy] n2ll has no attached job registry; ToyGenerator.__main__ must set n2ll._toy_jobs_by_id." + ) + job = jobs_by_id.get(job_id) + if job is None: + raise RuntimeError(f"[toy] BIT job '{job_id}' not found in the loaded config's jobs.") + return job + + +def _materialize_truth_weights(n2ll: N2LL, region_id: str, source: ProcessSource, feature_names: list, + split: Optional[str], hypothesis, allow_negative_weights: bool): + """Stream one ProcessSource's sample, reconstruct the exact truth weight, + restrict to `split`, rescale by the retained *weight* fraction, and (if + `hypothesis` is given) multiply by the model's region-level (1+T) at that + hypothesis -- the surrogate-route injection, evaluated on this source's own + kinematics via the same `_eval_region_surrogates`/`_compute_T_from_columns` + N2LL uses for real data (see the plan's "two POI injection routes"). + + Returns (X (M,d), w_truth (M,), diagnostics dict). + """ + region = _find_region(n2ll, region_id) + C = _find_class(region, source.class_id) + class_sample = C.get("sample") + if not class_sample: + raise RuntimeError(f"[toy] Class '{region_id}/{source.class_id}' has no configured 'sample'.") + sample_name = source.sample_name or class_sample + if source.sample_name is not None and source.sample_name != class_sample: + logger.warning( + "[toy:%s/%s] sample_name '%s' differs from the class's configured sample '%s' " + "-- this is a deliberate mismodelling injection.", + region_id, source.class_id, source.sample_name, class_sample, + ) + + provider = None + required_observers: list = [] + if source.coefficients is not None: + if hypothesis is not None: + poi_vals = {p.name: p.val for p in getattr(hypothesis, "POIs", [])} + overlap = [k for k in source.coefficients if poi_vals.get(k, 0.0) != 0.0] + if overlap: + raise ValueError( + f"[toy:{region_id}/{source.class_id}] POI(s) {overlap} injected via both " + f"'coefficients' and a nonzero hypothesis value; pick one route." + ) + from common.derivative_providers import build_derivative_provider + job = _class_bit_job(n2ll, region_id, source.class_id) + provider = build_derivative_provider(job) + expected_poi_order = n2ll._poi_order.get((region_id, source.class_id), []) + if expected_poi_order and list(provider.parameters) != list(expected_poi_order): + raise RuntimeError( + f"[toy:{region_id}/{source.class_id}] provider POI order {provider.parameters} " + f"!= class POI order {expected_poi_order}." + ) + required_observers = list(provider.required_observers) + + splitting_cfg = getattr(n2ll, "_toy_splitting_defaults", None) + uid_fields = list((splitting_cfg or {}).get("uid_fields", ["run", "luminosityBlock", "event"])) + uid_splitter = lo = hi = None + if split is not None: + if not splitting_cfg: + raise RuntimeError( + f"[toy:{region_id}/{source.class_id}] split='{split}' requested but n2ll has no attached " + f"splitting config (ToyGenerator.__main__ sets n2ll._toy_splitting_defaults)." + ) + split_names = [split] if isinstance(split, str) else list(split) + uid_splitter, uid_fields, (lo, hi) = cr._uid_split_interval(splitting_cfg, *split_names) + + observer_names = list(dict.fromkeys(required_observers + uid_fields)) + if source.weight_branches: + for weight_branch in source.weight_branches: + observer_names.append(weight_branch) + loader = n2ll.factory.get(sample_name).clone() + loader.setFeatures(feature_names, observer_names=observer_names) + + Xs, Ws = [], [] + sum_w_all = 0.0 + sum_w_split = 0.0 + for shard in range(int(getattr(loader, "n_split", 1))): + X, G, w = loader.materialize(shard=shard, what="fow") + X = np.asarray(X, dtype=np.float64) + G = np.asarray(G, dtype=np.float64) + w = np.asarray(w, dtype=np.float64) + sum_w_all += float(np.sum(w)) + + if provider is not None: + deriv_w = provider.truth_weight_matrix(G, w, observer_names) # (N, M), col 0 = nominal + w_coef = deriv_w[:, 0] + deriv_w[:, 1:] @ expand_pois_linear_quadratic( + provider.parameters, source.coefficients + ) + else: + w_coef = w + + w_truth = w_coef * source.scale_factor + if source.weight_branches: + branch_mask = N2LL.make_column_mask(observer_names, source.weight_branches) + w_truth = w_truth * np.prod(G[:, branch_mask], axis=1) + if source.weight_function is not None: + w_truth = w_truth * np.asarray(source.weight_function(X, feature_names), dtype=np.float64) + + if split is not None: + on2idx = {n: i for i, n in enumerate(observer_names)} + G_uid = G[:, [on2idx[f] for f in uid_fields]] + m_keep = uid_splitter.mask_from_np(G_uid, uid_fields, lo, hi) + else: + m_keep = np.ones(len(w), dtype=bool) + sum_w_split += float(np.sum(w[m_keep])) + + Xs.append(X[m_keep]) + Ws.append(w_truth[m_keep]) + + if sum_w_all <= 0: + raise RuntimeError(f"[toy:{region_id}/{source.class_id}] sample has zero total weight.") + f_weight = sum_w_split / sum_w_all + if f_weight <= 0: + raise RuntimeError(f"[toy:{region_id}/{source.class_id}] split='{split}' selected zero weight.") + + X_all = np.concatenate(Xs, axis=0) if Xs else np.empty((0, len(feature_names)), dtype=np.float64) + w_truth_all = (np.concatenate(Ws, axis=0) if Ws else np.empty(0, dtype=np.float64)) / f_weight + + nominal_yield = sum_w_all + if hypothesis is not None and len(X_all): + by_class_all = n2ll._eval_region_surrogates(region_id, X_all, feature_names) + cA_per_class = n2ll._assemble_cA_per_class(region_id, hypothesis) + nuA_per_group = n2ll._assemble_nuA_groups(region_id, hypothesis) + nu_vals = {p.name: p.val for p in getattr(hypothesis, "parameters", []) if not p.isPOI} + ln_bias = { + cid: sum(a * nu_vals.get(nm, 0.0) for nm, a in n2ll._lnN_by_class.get((region_id, cid), [])) + for cid in n2ll._class_ids_by_region[region_id] + } + rate_shift = n2ll._assemble_rate_shift_per_class(region_id, hypothesis) + T = n2ll._compute_T_from_columns(region_id, by_class_all, cA_per_class, nuA_per_group, + ln_bias, rate_shift, 0, len(X_all)) + w_truth_all = w_truth_all * (1.0 + T) + + neg = w_truth_all < 0.0 + n_neg = int(np.sum(neg)) + truth_yield = float(np.sum(w_truth_all)) + neg_yield_share = float(np.sum(w_truth_all[neg]) / truth_yield) if (n_neg and truth_yield) else 0.0 + if n_neg > 0: + if not allow_negative_weights: + raise RuntimeError( + f"[toy:{region_id}/{source.class_id}] {n_neg}/{len(w_truth_all)} events have negative " + f"truth weight (share of yield {neg_yield_share:.4f}); pass allow_negative_weights=True " + f"to clip them to zero." + ) + w_truth_all = np.where(neg, 0.0, w_truth_all) + truth_yield = float(np.sum(w_truth_all)) + + diag = { + "f_weight": f_weight, "split": split, + "nominal_yield": nominal_yield, "truth_yield": truth_yield, + "n_negative": n_neg, "negative_yield_share": neg_yield_share, + } + return X_all, w_truth_all, diag + + +# ============================================================================ +# truth-mode unbinned toy (multi-process) +# ============================================================================ + +def generate_unbinned_toy_from_truth(n2ll: N2LL, region_id: str, sources, rng: np.random.Generator, *, + split: tuple[str] | None = ("c2st_train", "c2st_val"), hypothesis=None, + allow_negative_weights: bool = False) -> dict: + """Multi-process truth-mode unbinned toy (see module docstring / plan + "Multi-process handling"). Per source (one per class): reconstruct the truth + weight, throw n_i ~ Poisson(w_truth_i), keep n_i > 0, record an origin label + (diagnostics only). Concatenate across sources and evaluate the region's + surrogates once on the whole unlabeled set, exactly as setObservation does for + real data. + """ + region = _find_region(n2ll, region_id) + classes = region.get("classes", []) or [] + by_cid = _resolve_process_sources(classes, sources) + feature_names = _region_feature_union(region) + + Xs, Ns, origins, diagnostics = [], [], [], {} + for i_class, (cid, source) in enumerate(by_cid.items()): + X, w_truth, diag = _materialize_truth_weights( + n2ll, region_id, source, feature_names, split, hypothesis, allow_negative_weights + ) + n_i = rng.poisson(np.clip(w_truth, 0.0, None)) + keep = n_i > 0 + n_kept = n_i[keep] + dup = n_kept > 1 + diag.update({ + "drawn_yield": float(n_kept.sum()), + "n_multiplicity_gt1": int(dup.sum()), + "multiplicity_gt1_yield_share": float(n_kept[dup].sum() / n_kept.sum()) if n_kept.sum() else 0.0, + }) + diagnostics[cid] = diag + logger.info("[toy:%s/%s] %s", region_id, cid, diag) + + Xs.append(X[keep]) + Ns.append(n_kept.astype(np.float64)) + origins.append(np.full(int(keep.sum()), i_class, dtype=np.int64)) + + X_toy = np.concatenate(Xs, axis=0) if Xs else np.empty((0, len(feature_names)), dtype=np.float64) + n_toy = np.concatenate(Ns, axis=0) if Ns else np.empty(0, dtype=np.float64) + origin_toy = np.concatenate(origins, axis=0) if origins else np.empty(0, dtype=np.int64) + + by_class = n2ll._eval_region_surrogates(region_id, X_toy, feature_names) + + return { + "X": X_toy, "w": n_toy, "by_class": by_class, "origin": origin_toy, + "feature_names": feature_names, "diagnostics": diagnostics, + } + + +# ============================================================================ +# binned toy +# ============================================================================ + +def generate_binned_toy(n2ll: N2LL, region_id: str, rng: np.random.Generator, *, + hypothesis=None, sources=None, split: tuple[str] | None = ("c2st_train", "c2st_val"), + allow_negative_weights: bool = False) -> dict: + """Binned toy: N_obs ~ Poisson(lambda) per bin, thrown once per bin (see the + plan's "Binned" section for why histogram-first is exact and cheaper than + per-event Poisson-then-histogram). + + Cache mode (sources is None): lambda = n2ll._compute_lambda_binned(...), the + model's expected counts. + Truth mode (sources given): lambda = the truth-weighted histogram summed over + sources' events, into the same bin edges -- so it differs from the cache-mode + lambda by exactly the ICH/ICPH mismodelling. + """ + un = n2ll._binned_unroll[region_id] + edges = un["edges"] + Nflat = len(un["flat_bins"]) + + if sources is None: + if hypothesis is None: + raise ValueError("[toy:binned] cache mode (sources=None) requires a hypothesis.") + lam = n2ll._compute_lambda_binned(region_id, hypothesis) + else: + region = _find_region(n2ll, region_id) + classes = region.get("classes", []) or [] + by_cid = _resolve_process_sources(classes, sources) + axes = un["axes"] + lam = np.zeros(Nflat, dtype=np.float64) + for cid, source in by_cid.items(): + X, w_truth, diag = _materialize_truth_weights( + n2ll, region_id, source, axes, split, hypothesis, allow_negative_weights + ) + logger.info("[toy:%s/%s] %s", region_id, cid, diag) + if len(edges) == 1: + H, _ = np.histogram(X[:, 0], bins=edges[0], weights=w_truth) + else: + H, _, _ = np.histogram2d(X[:, 0], X[:, 1], bins=[edges[0], edges[1]], weights=w_truth) + lam = lam + H.reshape(-1) + + neg = lam < 0.0 + n_neg = int(np.sum(neg)) + if n_neg > 0: + if not allow_negative_weights: + raise RuntimeError( + f"[toy:binned:{region_id}] negative expected counts in {n_neg}/{Nflat} bins " + f"(summed {float(np.sum(lam[neg])):.4g}); pass allow_negative_weights=True to clip them to zero." + ) + lam = np.where(neg, 0.0, lam) + + counts = rng.poisson(lam) + return {"counts": counts.astype(np.float64)} + + +# ============================================================================ +# top-level +# ============================================================================ + +def generate_toy(n2ll: N2LL, seed: int, *, source: str, hypothesis=None, truth_sources=None, + split: tuple[str] | None = ("c2st_train", "c2st_val"), throw_nuisances: bool = False, + allow_negative_weights: bool = False) -> dict: + """Generate one full toy (every unbinned + binned region of n2ll's likelihood). + + `source`: "cache" (model is truth) or "truth" (exact reweighting/scale/sample). + `hypothesis`: the injection point. In cache mode, defaults to nominal + (build_hypothesis_from_likelihood). In truth mode, None means no surrogate-route + (1+T) multiplication -- ProcessSource modifiers are still applied. + """ + if source not in ("cache", "truth"): + raise ValueError(f"source must be 'cache' or 'truth', got {source!r}") + if source == "cache" and truth_sources: + raise ValueError("truth_sources is only valid with source='truth'.") + + gen_hypothesis = hypothesis if hypothesis is not None else build_hypothesis_from_likelihood(n2ll.lk, name="toy") + constraint_centers = throw_constraint_centers(gen_hypothesis, _spawn_rng(seed, "__constraints__")) if throw_nuisances else {} + + unbinned_blocks: dict = {} + binned_counts: dict = {} + diagnostics: dict = {"unbinned": {}, "binned": {}} + + for region in n2ll.regions: + rid = region["id"] + region_rng = _spawn_rng(seed, rid) + if source == "cache": + toy = generate_unbinned_toy_from_cache(n2ll, rid, gen_hypothesis, region_rng, + allow_negative_weights=allow_negative_weights) + unbinned_blocks[rid] = { + "X": None, "w": toy["n"], "indices": toy["indices"], + "by_class": _rehydrate_cache_by_class(n2ll, rid, toy["indices"]), + } + else: + sources_for_region = (truth_sources or {}).get(rid, []) + toy = generate_unbinned_toy_from_truth( + n2ll, rid, sources_for_region, region_rng, split=split, hypothesis=hypothesis, + allow_negative_weights=allow_negative_weights, + ) + unbinned_blocks[rid] = {"X": toy["X"], "w": toy["w"], "by_class": toy["by_class"], "origin": toy["origin"]} + diagnostics["unbinned"][rid] = toy["diagnostics"] + + for region in n2ll.binned: + rid = region["id"] + region_rng = _spawn_rng(seed, rid) + if source == "cache": + toy = generate_binned_toy(n2ll, rid, region_rng, hypothesis=gen_hypothesis, + allow_negative_weights=allow_negative_weights) + else: + sources_for_region = (truth_sources or {}).get(rid, []) + toy = generate_binned_toy( + n2ll, rid, region_rng, sources=sources_for_region, hypothesis=hypothesis, + split=split, allow_negative_weights=allow_negative_weights, + ) + binned_counts[rid] = toy["counts"] + + return { + "seed": int(seed), + "source": source, + "hypothesis": {p.name: float(p.val) for p in gen_hypothesis.parameters}, + "constraint_centers": constraint_centers, + "unbinned_blocks": unbinned_blocks, + "binned_counts": binned_counts, + "diagnostics": diagnostics, + "config_version": getattr(n2ll, "version", None), + } + + +# ============================================================================ +# persistence +# ============================================================================ + +def save_toy(path: str, toy: dict) -> None: + """Persist a toy to HDF5. Cache-mode regions store only (indices, n) -- the + surrogate columns are verbatim cache slices, rehydrated by load_toy rather than + duplicated on disk. Truth-mode regions store the materialized X/by_class/origin. + """ + out_dir = os.path.dirname(os.path.abspath(path)) + os.makedirs(out_dir, exist_ok=True) + with h5py.File(path, "w") as f: + meta = f.create_group("meta") + meta.attrs["seed"] = int(toy["seed"]) + meta.attrs["source"] = toy["source"] + meta.attrs["hypothesis"] = json.dumps(toy["hypothesis"]) + meta.attrs["config_version"] = str(toy.get("config_version") or "") + meta.attrs["point"] = toy.get("point", "") + + cc = f.create_group("constraint_centers") + for name, val in (toy.get("constraint_centers") or {}).items(): + cc.attrs[name] = float(val) + + ub = f.create_group("unbinned") + for rid, block in toy["unbinned_blocks"].items(): + g = ub.create_group(rid) + g.create_dataset("n", data=np.asarray(block["w"], dtype=np.float64)) + if block.get("indices") is not None: + g.create_dataset("indices", data=np.asarray(block["indices"], dtype=np.int64)) + else: + g.create_dataset("X", data=np.asarray(block["X"], dtype=np.float64)) + g.create_dataset("origin", data=np.asarray(block["origin"], dtype=np.int64)) + bc = g.create_group("by_class") + for cid, comp in block["by_class"].items(): + cg = bc.create_group(cid) + for col_name, col in comp.items(): + cg.create_dataset(col_name, data=np.asarray(col, dtype=np.float64)) + + bn = f.create_group("binned") + for rid, counts in toy["binned_counts"].items(): + bn.create_dataset(rid, data=np.asarray(counts, dtype=np.float64)) + + logger.info("[toy] Saved %s", path) + + +def load_toy(path: str, n2ll: N2LL) -> dict: + """Load a toy saved by save_toy. Cache-mode regions are rehydrated by indexing + n2ll's live cache -- a shrunk/retrained cache raises IndexError (or a + _rehydrate_cache_by_class KeyError for a missing class) rather than silently + misaligning columns.""" + with h5py.File(path, "r") as f: + meta = f["meta"] + seed = int(meta.attrs["seed"]) + source = str(meta.attrs["source"]) + hypothesis = json.loads(meta.attrs["hypothesis"]) + config_version = str(meta.attrs.get("config_version", "")) + point = str(meta.attrs.get("point", "")) + + constraint_centers = {name: float(val) for name, val in f["constraint_centers"].attrs.items()} + + unbinned_blocks = {} + for rid in f["unbinned"]: + g = f["unbinned"][rid] + n = np.asarray(g["n"]) + if "indices" in g: + indices = np.asarray(g["indices"], dtype=np.int64) + unbinned_blocks[rid] = { + "X": None, "w": n, "indices": indices, + "by_class": _rehydrate_cache_by_class(n2ll, rid, indices), + } + else: + by_class = {} + for cid in g["by_class"]: + cg = g["by_class"][cid] + by_class[cid] = {col_name: np.asarray(cg[col_name]) for col_name in cg} + unbinned_blocks[rid] = { + "X": np.asarray(g["X"]), "w": n, "by_class": by_class, + "origin": np.asarray(g["origin"]), + } + + binned_counts = {rid: np.asarray(f["binned"][rid]) for rid in f["binned"]} + + return { + "seed": seed, "source": source, "hypothesis": hypothesis, + "constraint_centers": constraint_centers, + "unbinned_blocks": unbinned_blocks, "binned_counts": binned_counts, + "config_version": config_version, "point": point, + } + + +# ============================================================================ +# spec file parsing + CLI +# ============================================================================ + +def _resolve_weight_function(dotted_path: str) -> Callable: + module_name, func_name = dotted_path.rsplit(".", 1) + mod = importlib.import_module(module_name) + return getattr(mod, func_name) + + +def _parse_injection(injection: dict) -> dict: + truth_sources = {} + for region_id, classes in (injection or {}).items(): + sources = [] + for class_id, fields in (classes or {}).items(): + fields = dict(fields or {}) + wf = fields.pop("weight_function", None) + sources.append(ProcessSource( + class_id=class_id, + sample_name=fields.pop("sample_name", None), + coefficients=fields.pop("coefficients", None), + scale_factor=float(fields.pop("scale_factor", 1.0)), + weight_branches=fields.pop("weight_branches", None), + weight_function=_resolve_weight_function(wf) if wf else None, + )) + if fields: + raise RuntimeError(f"Unknown ProcessSource field(s) {list(fields)} for {region_id}/{class_id}.") + truth_sources[region_id] = sources + return truth_sources + + +def _hypothesis_from_point(n2ll: N2LL, point_hyp) -> Optional[Hypothesis]: + if not point_hyp: + return None + base = build_hypothesis_from_likelihood(n2ll.lk, name="toy_point") + for name, val in point_hyp.items(): + if name not in base: + raise KeyError(f"Unknown parameter '{name}' in point hypothesis.") + base[name].val = float(val) + return base + + +def _parse_seeds(spec: str) -> list: + seeds = [] + for part in spec.split(","): + part = part.strip() + if "-" in part: + lo, hi = part.split("-", 1) + seeds.extend(range(int(lo), int(hi) + 1)) + else: + seeds.append(int(part)) + return seeds + + +if __name__ == "__main__": + import argparse + import common.yaml_loader as yaml_loader + from fit.Likelihood import load_likelihood + + logging.basicConfig(level=logging.INFO, format="%(message)s") + + p = argparse.ArgumentParser(description="Toy dataset generation for pseudo-experiments.") + p.add_argument("configs", nargs="+", help="Path to one or more global YAML configs") + p.add_argument("--toySpec", required=True, help="Path to the toy spec YAML") + p.add_argument("--toyPoint", required=True, help="Name of the point in the spec to generate") + p.add_argument("--seeds", required=True, help="Seed or range, e.g. '0-499' or '3,7,12'") + p.add_argument("--outputDir", required=True, help="Directory to write _toy.h5 files") + p.add_argument("--overwrite", nargs="?", const="all", default=None, choices=["fit", "all"], + help="Overwrite the surrogate cache ('all') before generating.") + p.add_argument("--plot", action="store_true", + help="After generating, write per-feature diagnostic plots (config's " + "default_features, overlaid across the requested seeds) under " + "/toys///. Cache-mode toys pay the " + "cost of re-streaming the raw samples once to recover kinematics.") + args = p.parse_args() + + list_configs = [] + for config_path in args.configs: + aux_cfg = yaml_loader.load_yaml(config_path) + yaml_loader.print_summary(aux_cfg, config_path, yaml_loader._INCLUDE_TRACE) + yaml_loader.load_surrogates(aux_cfg, config_path, overwrite=False) + list_configs.append(aux_cfg) + cfg = yaml_loader.combine_configs(list_configs) + + like_info = load_likelihood(cfg) + + samples_mod = importlib.import_module(cfg["defaults"]["module_samples"]) + from common.yaml_loader import _resolve_features_list + default_features = cfg["defaults"].get("default_features", None) + features = _resolve_features_list(default_features) if default_features else None + factory = samples_mod.Factory( + features=features, + selection=cfg["defaults"].get("default_selection", None), + selection_features=cfg["defaults"].get("default_selection_features", None), + ) + + base = "_".join(os.path.splitext(os.path.basename(c))[0] for c in args.configs) + n2ll = N2LL( + like_info, factory=factory, + cache_subdir=os.path.join("NN2LCache", base, cfg["version"]), + cache_root=None, + overwrite=(args.overwrite == "all"), + ) + + spec = yaml_loader.load_yaml(args.toySpec) + spec_source = spec.get("source", "cache") + spec_split = spec.get("split", ("c2st_train","c2st_val")) + if isinstance(spec_split, str): + logger.info("spec_split is a string") + spec_split = [part.strip() for part in spec_split.split(",") if part.strip()] + else: + logger.info("spec_split is a list") + spec_throw_nuisances = bool(spec.get("throw_nuisances", False)) + spec_allow_negative = bool(spec.get("allow_negative_weights", False)) + + n2ll.shuffle_features = None + if spec_source == "cache": + n2ll.build_cache() + n2ll.prepare_runtime() + else: + logger.warning("Running with truth mode, skipping build_cache() and prepare_runtime(). If cache doesn't exist, running with --plot will crash.") + + n2ll.version = cfg.get("version") + n2ll._toy_splitting_defaults = (cfg.get("defaults") or {}).get("splitting") + n2ll._toy_jobs_by_id = {j["id"]: j for j in (cfg.get("jobs") or []) if j.get("id")} + + point = next((pt for pt in (spec.get("points") or []) if pt.get("name") == args.toyPoint), None) + if point is None: + available = [pt.get("name") for pt in (spec.get("points") or [])] + raise RuntimeError(f"Point '{args.toyPoint}' not found in {args.toySpec}. Available: {available}") + + truth_sources = _parse_injection(point.get("injection")) if spec_source == "truth" else None + hypothesis = _hypothesis_from_point(n2ll, point.get("hypothesis")) + + seeds = _parse_seeds(args.seeds) + os.makedirs(args.outputDir, exist_ok=True) + generated_toys = [] + for seed in seeds: + toy = generate_toy( + n2ll, seed, source=spec_source, hypothesis=hypothesis, truth_sources=truth_sources, + split=spec_split, throw_nuisances=spec_throw_nuisances, allow_negative_weights=spec_allow_negative, + ) + toy["point"] = args.toyPoint + out_path = os.path.join(args.outputDir, f"{args.toyPoint}_toy{seed}.h5") + save_toy(out_path, toy) + generated_toys.append(toy) + + if args.plot: + if not features: + raise RuntimeError("--plot requires defaults.default_features to be set in the config.") + import common.user as user + import common.syncer as syncer + from plot.toys.toy_diagnostic_plots import plot_toy_feature_distributions + + plot_dir = os.path.join(user.plot_directory, "toys", str(cfg.get("version")), args.toyPoint) + plot_toy_feature_distributions(n2ll, generated_toys, features, plot_dir) + syncer.sync() diff --git a/plot/bit/eft_modification_plot.py b/plot/bit/eft_modification_plot.py index 58994b74..ee9373ce 100644 --- a/plot/bit/eft_modification_plot.py +++ b/plot/bit/eft_modification_plot.py @@ -31,7 +31,7 @@ sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) import common.syncer as syncer -from common.derivative_providers import EFTDerivativeProvider +from common.derivative_providers import build_derivative_provider import modification_plotter as mp @@ -43,7 +43,7 @@ def main(): if not job.get("eft"): raise RuntimeError(f"Job '{job['id']}' has no 'eft' block; use pdf_modification_plot.py for PDF jobs.") - provider = EFTDerivativeProvider(job["eft"].get("parameters", [])) + provider = build_derivative_provider(job) mp.make_modification_plots(cfg, job, samples_mod, args, provider) diff --git a/plot/bit/pdf_modification_plot.py b/plot/bit/pdf_modification_plot.py index 7436126f..9345bac7 100644 --- a/plot/bit/pdf_modification_plot.py +++ b/plot/bit/pdf_modification_plot.py @@ -32,7 +32,7 @@ sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) import common.syncer as syncer -from common.derivative_providers import PDFDerivativeProvider +from common.derivative_providers import build_derivative_provider import modification_plotter as mp @@ -44,7 +44,7 @@ def main(): if not job.get("pdf"): raise RuntimeError(f"Job '{job['id']}' has no 'pdf' block; use eft_modification_plot.py for EFT jobs.") - provider = PDFDerivativeProvider(job["pdf"]) + provider = build_derivative_provider(job) mp.make_modification_plots(cfg, job, samples_mod, args, provider) diff --git a/plot/toys/toy_diagnostic_plots.py b/plot/toys/toy_diagnostic_plots.py new file mode 100644 index 00000000..51d0927c --- /dev/null +++ b/plot/toys/toy_diagnostic_plots.py @@ -0,0 +1,178 @@ +"""Diagnostic plots for generated toys: per-feature distributions overlaid across +toy seeds, for each unbinned region, using the config's default_features and +data/plot_options binning. Also draws the full Asimov-from-cache distribution +(the model's expected w0*(1+T) at the toys' generation hypothesis, summed over +the WHOLE cached MC, not just the sampled toy events) as a reference each toy +should fluctuate around -- always available, regardless of toy source, since the +surrogate cache is always built. + +For cache-mode toys (model is truth), raw kinematics are not part of the toy file +-- only cache row indices are, to avoid duplicating the surrogate cache -- so +features are recovered by re-streaming the region's samples once, at the union of +all requested toys' kept indices, in the same pass that also accumulates the +Asimov reference histogram (see fit.ToyGenerator.materialize_cache_region_diagnostics). +This is costly by design and meant for occasional diagnostic use, not routine toy +generation. + +For truth-mode toys, raw kinematics (X) are already part of the toy file, but only +for the region's own surrogate feature union (whatever the classifier/BIT/PNNs +consume) -- default_features not in that union are skipped with a warning. + +Entry point: fit/ToyGenerator.py --plot (see its __main__). +""" +from __future__ import annotations + +import os +import sys +import logging + +import numpy as np +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from matplotlib.patches import Patch +import mplhep as hep + +sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), "..", ".."))) + +import common.user as user +from data.plot_options import plot_options as DEFAULT_PLOT_OPTS +from fit.Likelihood import build_hypothesis_from_likelihood +from fit.ToyGenerator import _find_region, _region_feature_union, materialize_cache_region_diagnostics + +# Same as data.colors.cmap_petroff10_mpl, inlined to avoid importing that module's +# ROOT dependency here -- ROOT's cling JIT has been observed to segfault when +# initialized a second time in the same process after fit/ToyGenerator.py has +# already streamed samples (RDataLoader/PyROOT), so this purely-matplotlib +# diagnostic script stays ROOT-free. +COLOR_HEX = ['#3f90da', '#ffa90e', '#bd1f01', '#94a4a2', '#832db6', + '#a96b59', '#e76300', '#b9ac70', '#717581', '#92dadd'] + +logger = logging.getLogger(__name__) + +hep.style.use("CMS") + +from common.helpers import copyIndexPHP + +def _hypothesis_from_toy(n2ll, toy: dict): + """Rebuild the generation Hypothesis from a toy's stored {name: value} dict + (all parameters, not just nonzero ones -- see `generate_toy`'s return).""" + hyp = build_hypothesis_from_likelihood(n2ll.lk, name="toy_diagnostics") + for name, val in toy["hypothesis"].items(): + hyp[name].val = float(val) + return hyp + + +def _region_diagnostics(n2ll, region_id: str, toys: list, feature_names: list, plot_opts: dict): + """Gathers everything plot_toy_feature_distributions needs for one region: + (available_features, by_seed={seed: (X, w)}, asimov_hist={feature: histogram}). + + Cache-mode toys can recover any raw feature by name (via re-streaming); truth- + mode toys are limited to whatever the region's own surrogates consume. The + Asimov reference histogram is always computed, regardless of toy source, in + the same re-streaming pass that recovers cache-mode toys' raw kinematics. + """ + blocks = [(t["seed"], t["unbinned_blocks"][region_id]) for t in toys] + cache_blocks = [(seed, b) for seed, b in blocks if b.get("indices") is not None] + + region = _find_region(n2ll, region_id) + if cache_blocks: + available = list(feature_names) + else: + union = set(_region_feature_union(region)) + available = [f for f in feature_names if f in union] + skipped = [f for f in feature_names if f not in union] + if skipped: + logger.warning( + "[toy:plot:%s] Skipping feature(s) %s -- not part of this truth-mode toy's " + "stored kinematics (region surrogates don't consume them).", region_id, skipped, + ) + + union_indices = None + if cache_blocks: + union_indices = np.unique(np.concatenate([np.asarray(b["indices"]) for _, b in cache_blocks])) + + hypothesis = _hypothesis_from_toy(n2ll, toys[0]) + asimov_hist, union_X = materialize_cache_region_diagnostics( + n2ll, region_id, hypothesis, available, plot_opts, extra_indices=union_indices, + ) + + cache_X_by_seed = {} + if union_indices is not None: + pos_of = {int(idx): i for i, idx in enumerate(union_indices)} + for seed, b in cache_blocks: + rows = [pos_of[int(i)] for i in b["indices"]] + cache_X_by_seed[seed] = union_X[rows] + + truth_feature_union = _region_feature_union(region) + by_seed = {} + for seed, b in blocks: + w = np.asarray(b["w"], dtype=np.float64) + if seed in cache_X_by_seed: + X = cache_X_by_seed[seed] + else: + full_X = np.asarray(b["X"], dtype=np.float64) + col_positions = [truth_feature_union.index(f) for f in available] + X = full_X[:, col_positions] + by_seed[seed] = (X, w) + return available, by_seed, asimov_hist + + +def plot_toy_feature_distributions(n2ll, toys: list, feature_names: list, out_dir: str, + plot_opts: dict = DEFAULT_PLOT_OPTS) -> None: + """One figure per (unbinned region, feature): the full Asimov-from-cache + distribution (grey filled reference, at the toys' generation hypothesis, summed + over the whole cached MC) behind each requested toy seed's weighted step + histogram, binned per data/plot_options. Binned regions have no per-event + kinematics and are skipped.""" + os.makedirs(out_dir, exist_ok=True) + seeds = sorted(t["seed"] for t in toys) + colors = {seed: COLOR_HEX[i % len(COLOR_HEX)] for i, seed in enumerate(seeds)} + + for region in n2ll.regions: + rid = region["id"] + plottable = [f for f in feature_names if f in plot_opts] + skipped_opts = [f for f in feature_names if f not in plot_opts] + if skipped_opts: + logger.warning("[toy:plot:%s] No plot_options entry for %s, skipping.", rid, skipped_opts) + if not plottable: + continue + available, by_seed, asimov_hist = _region_diagnostics(n2ll, rid, toys, plottable, plot_opts) + if not available: + continue + + region_dir = os.path.join(out_dir, rid) + copyIndexPHP(region_dir) + os.makedirs(region_dir, exist_ok=True) + for i_feat, feat in enumerate(available): + n_bins, x_lo, x_hi = plot_opts[feat]["binning"] + edges = np.linspace(x_lo, x_hi, int(n_bins) + 1, dtype=np.float64) + tex = "$" + plot_opts[feat]["tex"].replace("#", "\\") + "$" + + fig, ax = plt.subplots(figsize=(8.0, 6.5)) + asimov = asimov_hist[feat] + ax.fill_between(edges[:-1], 0.0, asimov, step="post", color="0.88", zorder=0) + ax.stairs(asimov, edges, color="0.35", linewidth=1.2, zorder=1) + for seed, (X, w) in by_seed.items(): + hist, _ = np.histogram(X[:, i_feat], bins=edges, weights=w) + ax.stairs(hist, edges, color=colors[seed], linewidth=1.6, label=f"toy {seed}", zorder=2) + + if plot_opts[feat].get("logY", False): + ax.set_yscale("log") + ax.set_xlabel(tex) + ax.set_ylabel("Events") + handles, labels = ax.get_legend_handles_labels() + handles = [Patch(facecolor="0.88", edgecolor="0.35", label="Asimov (cache)")] + handles + labels = ["Asimov (cache)"] + labels + n_col = 1 if len(seeds) <= 8 else 2 + ax.legend(handles, labels, frameon=False, fontsize=8, ncol=n_col, loc="upper right") + + hep.cms.label("Internal", data=False, ax=ax, loc=0) + fig.tight_layout() + stub = os.path.join(region_dir, feat) + plt.savefig(stub + ".png", dpi=200) + plt.savefig(stub + ".pdf") + plt.close(fig) + + logger.info("[toy:plot] Wrote diagnostic plots to %s", out_dir) diff --git a/user/ricardo/claude/plans/toy-generation-design-decisions.md b/user/ricardo/claude/plans/toy-generation-design-decisions.md new file mode 100644 index 00000000..ade7cfbb --- /dev/null +++ b/user/ricardo/claude/plans/toy-generation-design-decisions.md @@ -0,0 +1,125 @@ +# Toy generation: design decisions + +Rejected alternatives for the work planned in +`toy-generation-for-pseudo-experiments.md`. Recorded for auditability. Several are also noted inline where they bear on a +specific section; this consolidates them with the reasoning. Design settled +2026-07-27. + +**Statistical / physics choices** + +- *Generating toys from an MC subset statistically independent of the one + supplying the expected-yield term.* Rejected. It would remove the + toy/template correlation, but halves effective MC statistics and requires + reweighting the two halves. The correlation is a sub-percent effect on the + widths at current MC statistics. Revisit only if MC statistics grow enough + that the cost disappears, or if the widths are ever shown to be biased at the + percent level. +- *Fixing nuisances at their true values instead of throwing global + observables.* Rejected. Coverage would be optimistic wherever systematics + matter, which is the regime of interest. Cost is the `constraint_center` + addition to `fit/Modeling.py`, which is small and reduces exactly to current + behaviour at centre zero. +- *Applying a UID split in cache mode.* Rejected on both grounds: it would + require UID columns in the HDF5 cache, a format change forcing every cache to + be rebuilt, and it is unnecessary because the model *is* the surrogate by + definition, so training leakage cannot bias coverage of the model against + itself. Only MC granularity is affected, which argues for using everything. +- *Rescaling a split by the nominal bucket fraction (0.20) rather than the + measured weight fraction.* Rejected. UID hashing distributes buckets, not + weights, so the nominal fraction injects a normalization bias that would + masquerade as exactly the rate bias under measurement. +- *Per-event Poisson then histogramming, for binned truth toys.* Rejected in + favour of one draw per bin. The two are distributionally identical — a sum of + independent Poissons is Poisson with the summed mean, and regions are disjoint + selections so there is no shared-event correlation — so histogram-first is + strictly cheaper. +- *Re-reading ROOT and re-evaluating surrogates per toy in cache mode.* + Rejected. Physically identical to slicing the cached columns, but far slower, + since `g`, `R` and `Delta::*` are already stored per event. The genuinely + different capability this seemed to offer became truth mode instead. + +**Scope and interface choices** + +- *A `toys:` block in the analysis YAML config.* Rejected in favour of a + standalone `--toySpec` file. The injection is per-invocation, not part of the + statistical model — hundreds of runs share one model — whereas + `classifier.asimov` is in YAML precisely because it defines the model's + template MC. A new top-level block would also need validation in + `_apply_defaults_and_checks`, a merge rule in `combine_configs()` and a + `print_summary` entry, all in `common/yaml_loader.py`, which every training + entry point depends on. Since CLAUDE.md forbids compatibility shims, a later + schema change would mean editing every config that adopted it. Precedent for + the standalone file: `--rotate`. Promoting the spec into the config later is + mechanical if it proves to be model after all. +- *A toy-fit driver, Slurm fan-out, and merge/plot step.* Deliberately out of + scope. The generator plus `fit/ToyGenerator.py.__main__` is enough to produce + and refit toys; the fan-out pattern already exists in `fit/slurm_utils.py` + when wanted. +- *CLI grammar `--toyScale REGION:CLASS=FACTOR`, `--toySample`, + `--toyCoefficients`, `--toyHypothesis`, `--toySource`, `--throwNuisances`.* + Rejected in favour of named points inside the spec file. The ad-hoc + `REGION:CLASS=` parsing was error-prone, and nine flags collapsed to four once + the spec defines the ensemble and the CLI only selects a member of it. Point + names also land in output filenames for free. +- *A `ToyDataset` in-memory loader class.* Rejected. Once `setObservationArrays` + exists nothing consumes the loader contract; `data/toy_data.py:_ArrayLoader` + does not expose the required `n_split`; and routing through `setObservation` + would hit its `ignore_weights=True` default and overwrite the multiplicities + with ones. +- *An `apply_model_hypothesis` flag on truth-mode generation.* Rejected as + redundant — the presence of a `hypothesis` block on a spec point is the + switch. +- *The name `setObservationFromBlocks`.* Rejected for `setObservationArrays`. + "Block" is real precedent ([Likelihood.py:1710](fit/Likelihood.py#L1710)) but + undocumented jargon at a public signature; the meaningful distinction from + `setObservation` is that it takes already-evaluated arrays rather than + loaders. + +**Simplicity pass, 2026-07-28** + +Added after a cold-read review against the "Code simplicity" section of +CLAUDE.md, which postdates most of the plan. + +- *Toy generation in `fit/Likelihood.py.__main__` as well as `ToyGenerator.py`.* + Rejected. Once standalone generation exists — justified by generate-once / + fit-many-ways, which avoids comparing fit variants across two different random + ensembles — a second generating entry point duplicates spec parsing, point + selection and seeding with no failure of its own to justify it, and + contradicts the plan's "no toy-fit driver" scope line. The fitter takes + `--toyFile` only, which also keeps `importlib` resolution of a user-supplied + dotted path out of the fit entry point. +- *Persisting `g`/`R`/`Delta::*` per cache-mode toy.* Rejected. Those rows are + verbatim slices of `n2ll._h5`, so 500 toys would write roughly 500 copies of + the cache. Storing `(indices, multiplicities)` and rehydrating at load is + strictly smaller and makes column misalignment impossible by construction + rather than by validation. +- *A `/meta` binding payload of region and class ids, `nA`/`nB` dims and + delta-group ids, cross-checked on load.* Rejected as mostly redundant once + cache-mode toys index the live cache. Truth mode asserts shapes against the + loaded predictors instead — checking the real runtime objects rather than a + schema that must be kept in sync with them. +- *`Toy` as a class.* Rejected for a plain dict plus `save_toy` / `load_toy`. + No state outlives a call and there is no behaviour beyond serialization. +- *A `derivative_job` override on `ProcessSource`.* Rejected. Justified only by + a hypothetical ("e.g. a different truncation") with no study behind it and no + verification item. The provider comes from the class's configured BIT job. + +**Considered in the same pass and deliberately kept** + +Recorded because a cold reader will propose cutting these again — they look +unjustified without the context below. + +- *`sample_name`, `weight_branches`, `weight_function` on `ProcessSource`.* All + three are explicit user requirements: the three mechanisms chosen for + supplying truth density (an alternative sample, extra weight branches, an + arbitrary callable). "No named use case in the plan" is true and beside the + point — the requirement predates the plan. +- *`allow_negative_weights`.* Kept. Crash-and-report is the default; the flag + exists so a study can proceed knowingly at large `|c|` where the truncated + quadratic goes negative. +- *`split: null`.* Kept. Verification 6 needs it to make pull *widths* + comparable between truth and cache mode; without it that comparison is + mean-only. +- *Named `points:` in the spec.* Kept. One file per injection point turns an + injection scan into a directory of near-duplicate files; point names also land + in output filenames for free. diff --git a/user/ricardo/claude/plans/toy-generation-for-pseudo-experiments.md b/user/ricardo/claude/plans/toy-generation-for-pseudo-experiments.md new file mode 100644 index 00000000..b7ed9153 --- /dev/null +++ b/user/ricardo/claude/plans/toy-generation-for-pseudo-experiments.md @@ -0,0 +1,760 @@ +# Toy dataset generation for pseudo-experiments + +## Context + +GOLLUM currently supports exactly two data modes in `N2LL`: Asimov +(`setAsimov`, [Likelihood.py:1384](fit/Likelihood.py#L1384)) and observed data +(`setObservation`, [Likelihood.py:1448](fit/Likelihood.py#L1448)). There is no +toy/pseudo-experiment generation anywhere in `fit/` — the only `np.random` call +is the unseeded feature shuffle used by the `--shuffle` diagnostic +([Likelihood.py:872](fit/Likelihood.py#L872)). + +Without toys we cannot (a) measure bias of the extracted PDF coefficients, +(b) measure coverage of the reported intervals, or (c) test whether the +`-2dlnL` test statistic follows the asymptotic chi-square distribution, which +is not guaranteed for a POD-basis fit with many weakly-constrained directions. + +Two **sources** are needed — the spec key is `source`, and "mode" is reserved +throughout for `N2LL`'s existing Asimov/observation distinction. They are **not** +interchangeable: + +1. **`source: cache`** — the model is taken as truth. Fast, and answers "given + that the surrogates are correct, are my intervals calibrated?" +2. **`source: truth`** — pseudo-data is generated from an exact reweighting, a + per-process rescaling, or a different sample, while the model still fits + with the trained surrogates. Probes surrogate mismodelling as a bias in the + fitted coefficients. + +Below, "cache mode" and "truth mode" are shorthand for these two `source` +values. + +Decision rule: **anything expressible as a parameter value belongs in cache +mode.** A normalization scaling is already a model parameter — the floating +`nu_norm_ttbar` lnN in `configs/unbinned_v6_rate/unbinned_2016_rate.yaml:52-53`, or +a `rate_shift` POI ([Likelihood.py:654](fit/Likelihood.py#L654)) — so +a spec point carrying `hypothesis: {nu_norm_ttbar: 1.0}` generates it exactly, +with no MC re-reading. +Truth mode is for shifts the model has **no parameter to absorb**, e.g. scaling +a background that carries no rate nuisance and measuring the induced bias on +`c0..c5`. + +One deliberate exception: a surrogate-route point used as the *control* for a +truth point must run in truth mode too, so both draw from the same events. See +§"The two POI injection routes". + +Scope: the **generator** plus the minimum `N2LL` / `Modeling` surface needed to +fit a toy. No toy-fit driver and no Slurm fan-out. + +## Statistical content + +`rho(x) = dN/dx` denotes the intensity of the unbinned Poisson point process — +the differential expected yield, whose integral is the total expected yield. +It is not a normalized PDF. `rho_0` is the intensity at `c = 0, nu = 0`. + +### Unbinned, cache mode + +The model factorizes as `rho(x) = rho_0(x) * (1 + T(x; c, nu))`, with `T` from +`_compute_T_chunk` ([Likelihood.py:1342](fit/Likelihood.py#L1342)); `1 + T` is +the differential cross-section ratio the BIT/PNN surrogates learn. Writing out +the extended unbinned NLL, + + -lnL = int rho - sum_obs ln rho(x_i) + = int rho_0 + int rho_0*T - sum_obs ln rho_0 - sum_obs log1p(T) + +the first and third terms are parameter-independent and drop, leaving +`int rho_0*T - sum_obs log1p(T)`. The first is estimated as `sum_i w0_i * T_i` +over the cached MC — literally [Likelihood.py:1742](fit/Likelihood.py#L1742), +with the per-event term at [Likelihood.py:1775](fit/Likelihood.py#L1775). +**`rho_0` is never evaluated pointwise**; the cache is only its MC estimator, +`sum_i w0_i delta(x - x_i)`. + +Symbols used throughout: `w0` is the MC event weight; `g` the classifier +probability for a class; `R` the BIT ratio basis; `Delta::` a PNN group's +basis matrix. All four are the per-event columns the cache stores. + +A sum of independent Poissons placed on those MC points is exactly a Poisson +point process with that intensity, so a toy is + + n_i ~ Poisson( w0_i * (1 + T_i(c', nu')) ) + +Rows with `n_i > 0` are the toy events, carrying observed weight `n_i`. Since +the likelihood consumes `sum_i w_i log1p(T_i)`, multiplicities can be carried +as weights — no `np.repeat`. No surrogate re-evaluation is needed, since `g`, +`R` and `Delta::*` are already cached per event. + +### Unbinned, truth mode + +Drawn from an externally supplied weight instead of the surrogate prediction: +`n_i ~ Poisson(w_truth_i)`, with `g`, `R`, `Delta::*` then evaluated by the +trained surrogates on the selected events. The expected-yield term is unchanged, +so both the shape and the normalization of any mismodelling propagate into the +fit — which is the point. + +### Binned + +In cache mode, `_compute_lambda_binned(rid, hyp)` +([Likelihood.py:1236](fit/Likelihood.py#L1236)) already returns the expected +count vector, so `N_obs ~ Poisson(lambda(c', nu'))` into +`_obs_binned_counts[rid]`. Needed because the v6 configs mix an unbinned +`SR_` with a binned `CR_` in one fit. + +In truth mode the binned regions must follow the same truth as the unbinned +ones. The substantive choice is **which lambda is thrown on**, not where the +Poisson is applied: + +- `_compute_lambda_binned` returns the **model's** expected counts, from the + ICH/ICPH surrogates. +- Histogramming the truth-weighted events into `_binned_unroll[rid]['edges']` + (the histogramming `setObservation` already does) gives the **truth's** + expected counts, `lambda_truth_b = sum_{i in b} w_truth_i`. + +These differ by exactly the ICH/ICPH mismodelling, which is what truth mode +exists to measure, so truth mode throws on `lambda_truth`. Using the model's +lambda would leave binned regions model-truth while unbinned regions are +reweighted, making the resulting bias uninterpretable. + +Throw **once per bin**, not per event. Per-event `n_i ~ Poisson(w_truth_i)` +followed by histogramming is a sum of independent Poissons, hence exactly +`Poisson(lambda_truth_b)` — distributionally identical, and regions are disjoint +selections so there is no shared-event correlation. Histogram-first is therefore +strictly cheaper: one draw per bin rather than one per MC event, with no +multiplicity array to materialize. + +### Nuisance throwing + +Each toy throws `nu_obs ~ Normal(nu_true, 1)` per penalized nuisance, with the +constraint becoming `(nu - nu_obs)^2`. `Hypothesis.penalty()` +([Modeling.py:205](fit/Modeling.py#L205)) hardcodes `sum(p.val**2)`, so this +needs the change below. Floating nuisances get `isPenalized=False` +([Likelihood.py:338-341](fit/Likelihood.py#L338-L341)), so free rate parameters +such as `nu_norm_ttbar` are correctly never thrown. + +### UID splitting: which events generate toys + +`build_cache` applies **no** UID filtering — `_iter_asimov_batches` +([Likelihood.py:717](fit/Likelihood.py#L717)) materializes every shard, and +nothing under `fit/` consumes `defaults.splitting`. Meanwhile the configs +reserve `final_eval: {fraction: 0.20}`, which **no script consumes at all** — +repo-wide it appears in exactly two places, both comments +([pnn_training_closure_mpl.py:320](ML/PNN/pnn_training_closure_mpl.py#L320), +[modification_plotter.py:347](plot/bit/modification_plotter.py#L347)). + +Verified that it is genuinely held out: the BIT trains on the PNN's buckets — +`bit_train_key = "pnn_train"`, `bit_val_key = "pnn_val"` at +[pdf_bit_training.py:158-159](ML/BIT/pdf_bit_training.py#L158-L159) and +[eft_bit_training.py:164-165](ML/BIT/eft_bit_training.py#L164-L165) — and the +TFMC does the same at +[tfmc_training.py:184-185](ML/TFMC/tfmc_training.py#L184-L185). No trainer of any +surrogate the likelihood uses touches `final_eval`. + +The right choice differs by mode: + +**Cache mode — full dataset, no split.** The model *is* the trained surrogate by +definition, and the toy is a draw from the model's own intensity +`w0 * (1 + T)`. Whether an event was in training does not change what the model +is; it only changes which points carry the MC estimator of `rho_0`. Coverage of +the model against itself is therefore unaffected by overfitting, while MC +granularity argues for using everything. Memorized structure in `T` on training +events appears identically in the toy and in the expected-yield term and +cancels — that belongs to the MC-reuse caveat below, not to a bias. + +**Truth mode — `final_eval` by default.** Here the events carry the truth +density and the study measures the gap between exact `w(c)` and surrogate +`R(c)`. On training events the surrogates fit better than on fresh events, so +the measured bias is optimistic in precisely the quantity being measured. +Costless to apply, since truth mode re-reads from loaders: add `uid_fields` to +`required_observers` and mask with `mask_from_np`, the pattern at +[calibration_runner.py:168-169](ML/Calibration/calibration_runner.py#L168-L169) +and [calibration_runner.py:201](ML/Calibration/calibration_runner.py#L201). + +Bucket intervals are built from the scheme the way `_uid_c2st_intervals` +([calibration_runner.py:96](ML/Calibration/calibration_runner.py#L96)) does, +but that helper hardcodes the `c2st_train`/`c2st_val` keys. Generalize it to +take the split name as an argument and change its one existing call site over, +rather than duplicating the interval arithmetic. The scheme is read from the +merged config's `defaults.splitting`, which is what the fit path already +carries — not from a per-job block, since toy generation is not job-scoped. + +**Normalization must be restored.** Generating from 20% of events while the +expected-yield term `sum w0 * T` runs over the full cache would put the toy +yield at a fifth of the expectation and pull the fit violently. Toy weights are +rescaled by `1 / f_weight` with + + f_weight = sum(w over split) / sum(w over all) + +the **weight** fraction, not the nominal bucket fraction `0.20` — UID hashing +distributes buckets, not weights, and using `0.20` would inject a normalization +bias masquerading as the rate bias under measurement. This needs no extra pass: +truth mode already streams every event to evaluate the mask, so both sums +accumulate in that pass while non-split events are discarded. + +The cost is granularity: rescaling by ~5x means each retained event stands for +five times more data events, so `n_i ~ Poisson(5 * w_i)` yields many more +duplicates — negligible for `TTLep_pow_*`, potentially severe for a small +background. The `n_i > 1` diagnostic must therefore be reported **after** +rescaling. + +Spec key `split` (default `final_eval`, `null` for the whole dataset) keeps the +escape hatch for a purely statistical study where leakage is irrelevant and +granularity dominates. + +Applying a split in cache mode would require UID columns in the HDF5 cache, a +format change forcing every cache to be rebuilt. Not done here; noted as the +extension path if it is ever wanted. + +### Known caveat (documented, not fixed) + +Toys are drawn from the same cached MC that supplies the expected-yield term, so +toy and template fluctuations are correlated. Sub-percent on the widths at +current MC statistics. Record in the module docstring; add no machinery. + +## Multi-process handling (truth mode) + +`build_cache` streams **all** `classifier.asimov` samples into one concatenated +event set and writes per-class `g`/`R`/`Delta::*` columns for the *same rows*, +sharing one `w0` (per-class columns written together in `build_cache`, +[Likelihood.py:786](fit/Likelihood.py#L786)). The cache +therefore represents the **total** nominal intensity, and the process +decomposition lives entirely in the classifier output `g_p(x)` +(`T = sum_p g_p * (...)`, [Likelihood.py:1364](fit/Likelihood.py#L1364)). The +likelihood never sees a process label. Reference multi-class region: +`configs/Eta_unbinned/Eta_unbinned_2016.yaml:24`, +`asimov: [TTLep_pow_2016, EtaS_2016, EtaP_2016]`. + +So a truth toy is built **per process** — that is where per-process weights and +truth modifications live — then concatenated and pushed through the surrogates +**once, as a single unlabeled event set**, exactly as real data is treated. + +```python +@dataclass +class ProcessSource: + """One process's contribution to a truth-mode toy in one region.""" + class_id: str # must match a region class id + sample_name: str | None = None # default: the class's configured `sample` + coefficients: dict[str, float] | None = None # PDF/EFT truth injection + scale_factor: float = 1.0 + weight_branches: list[str] | None = None + weight_function: Callable | None = None # f(X, feature_names) -> (N,) +``` + +Modifiers are optional and multiplicative, so they compose: + + w_truth = w_coefficients * scale_factor * prod(weight_branches) * weight_function(X) + +reducing to the nominal weight when `coefficients` is None. + +### Exact PDF/EFT weights (`coefficients`) + +The BIT is trained only on *derivatives*, so no branch on any sample carries the +full modified weight `w(c)` and no alternative sample provides it. It has to be +reconstructed — but the machinery already exists and must be **reused, not +re-derived**. + +`provider.truth_weight_matrix(G, w, observer_names)` in +`common/derivative_providers.py` returns an `(N, M)` matrix aligned to +`provider.combinations = [(), (p,)..., (p,q) for p<=q]`, already in weight +units, with column 0 the nominal weight +([derivative_providers.py:69](common/derivative_providers.py#L69) returns +`deriv * w`; [EFTWeightInterface.py:59](eft/EFTWeightInterface.py#L59) sets +`out[0] = nominal_weight`). Columns `1:` are the true first and symmetrized +second derivatives. The Taylor sum is then exactly `expand_pois_linear_quadratic` +([Likelihood.py:383](fit/Likelihood.py#L383)): + +```python +w_coefficients = deriv_w[:, 0] + deriv_w[:, 1:] @ expand_pois_linear_quadratic( + provider.parameters, coefficients) +``` + +The orderings match term for term: `expand_pois_linear_quadratic` emits linears +then `i<=j` quadratics in `combinations_with_replacement` order, identical to +`provider.combinations` after dropping column 0, and its `0.5 if i==j else 1` +factor is the `(1/2) sum_{a,b} -> sum_{a<=b}` rewrite explained in the comment at +[Likelihood.py:388-393](fit/Likelihood.py#L388-L393). `PODBasis.derivatives` +returns true second derivatives (`H = outer + outer.T`, and the parametrization +is linear in `c`, so `d2f/dc_a dc_b = 0` individually). + +**The convention is inherited, not re-derived** — verified: BIT training builds +its targets from the identical calls, `pdf.derivatives(x1=..., x2=..., id1=..., +id2=..., Q=...)` at [pdf_bit_training.py:251](ML/BIT/pdf_bit_training.py#L251) +and `eft.make_weight_matrix(G, obs_names, w)` at +[eft_bit_training.py:266](ML/BIT/eft_bit_training.py#L266), which are exactly +what `PDFDerivativeProvider` / `EFTDerivativeProvider` wrap. Truth weight and +model ratio `R = 1 + c_A R_A` therefore agree by construction rather than by +getting an algebra factor right independently — and this holds for EFT without +needing to know whether the `der_*` branches carry a factor of one half. + +This is the sharpest available surrogate test: generate from the exact `w(c)`, +fit with the BIT-predicted `R(c)`; residual bias in `c0..c5` is purely BIT +mismodelling plus MC statistics. + +Four requirements this imposes: + +1. **Observers must be requested.** For PDF this is nearly free: the five + `Generator_*` branches *and* the UID fields are already in + `observables.OBSERVERS` ([observables.py:103](data/observables.py#L103)), + which the samples set as `observer_names` + ([samples_RunII.py:98](data/samples_RunII.py#L98)). EFT needs + `EFTWeight_SM` + `der_*`, which are not in that default set. Either way the + generator must `clone()` the loader + ([RDataLoader.py:638](data/RDataLoader.py#L638)) and call + `setFeatures(..., observer_names=provider.required_observers + uid_fields)` + explicitly, following + [calibration_runner.py:171](ML/Calibration/calibration_runner.py#L171) — + cloning so the shared factory loader is never mutated, and explicit so a + sample lacking a branch fails loudly via `strict_branches` rather than + silently. Materialize with `what="fow"`. +2. **POI order must be asserted.** `self._poi_order[(rid, cid)]` must equal + `provider.parameters`, or the c-vector contracts against the wrong columns + silently. Hard assert, no fallback. +3. **Negative truth weights are physical here.** A truncated quadratic goes + negative at large `|c|` — truncation breaking down, not a bug. Report the + fraction and the negative yield share rather than clipping silently. +4. **The provider comes from the class's configured BIT job**, so truth and + model share a parametrization by construction. There is no override — a + study that wants a deliberately different parametrization can point + `sample_name` elsewhere or edit the config, and no current study needs it. + +Cost note: the PDF provider evaluates the parametrization per event +(`PODBasis.evaluate`), so this is materially slower than cache mode — once per +source per toy. + +`truth_sources` is `dict[region_id, list[ProcessSource]]`. The `injection:` block +of a spec point (§5) deserializes into exactly this — region id, then class id, +then the `ProcessSource` fields — so `injection`, `truth_sources` and +`ProcessSource` are three views of one object. **Any class with no +explicit source gets an implicit one at `scale_factor=1.0` from its own +configured sample**, so the common case is one line: + +```python +truth_sources = {"SR_2016": [ProcessSource(class_id="DY_2016", scale_factor=1.2)]} +``` + +Per-region algorithm: + +1. Resolve the class list from `n2ll.regions` (unbinned) or `n2ll.binned` + (binned); merge explicit sources with implicit defaults. Reject unknown + `class_id`s and duplicates. + The loaders are given the **union** of every feature list the region's + surrogates consume — classifier, BIT and every PNN — because + `_eval_region_surrogates` calls `make_column_mask` for each, and each mask + must resolve against the toy `X`. Build the union from the loaded + predictors' `feature_names`, not from `job["features"]`. +2. Per source: `factory.get(sample_name).clone()`, then + `setFeatures(..., observer_names=provider.required_observers + uid_fields)`. + Stream shards with `materialize(what="fow")`, accumulating `sum(w)` over all + events and `sum(w)` over the `split` mask; apply the modifiers to the masked + events, rescale by `1 / f_weight`, throw `n_i ~ Poisson(w_truth_i)`, keep + `n_i > 0`, record a per-event **origin label** (diagnostics and plotting + only, never used in the likelihood). +3. Assert all sources share identical `feature_names` — the same check + `_iter_asimov_batches` makes at [Likelihood.py:731](fit/Likelihood.py#L731). +4. Concatenate to `X_toy (M,d)`, `n_toy (M,)`, `origin_toy (M,)`. +5. `by_class = n2ll._eval_region_surrogates(rid, X_toy, feature_names)`, one + call on the whole set. +6. Emit `{'X': X_toy, 'w': n_toy, 'by_class': by_class}` — the `_obs_unbinned` + layout. + +A `sample_name` differing from the class's configured sample is legal but warns: +the surrogates were trained on the configured sample, so this is deliberately a +mismodelling injection. + +Per-source diagnostics, printed and stored: `f_weight` and the split used, +nominal yield, truth yield, drawn +yield, negative-weight fraction, and **the count of rows with `n_i > 1` plus +their share of the yield**. The last matters for low-statistics backgrounds — a +DY event with `w = 8` yields ~8 duplicates at one phase-space point, and a lumpy +toy is an MC-stat artifact, not a fluctuation to read as coverage failure. + +## Changes + +### 1. `fit/Modeling.py` — constraint centers + +- `ModelParameter.__init__`: add keyword `constraint_center=0.0`, stored as a + plain float (keeps `penalty()` differentiable in `run_autograd_fit`). +- `Hypothesis.penalty()` becomes + `sum((p.val - p.constraint_center)**2 for p in self.parameters if p.isPenalized)`. + With the default centre this is bit-identical to today. +- `Hypothesis.set_constraint_centers(mapping)`: set by name, raising on an + unknown name or a non-penalized target. +- `Rotated`: mirrored nuisance parameters + ([Modeling.py:338-342](fit/Modeling.py#L338-L342)) copy `constraint_center` + for correct printing. `Rotated.penalty()` already forwards to `_base` + ([Modeling.py:390](fit/Modeling.py#L390)), so nothing else changes. + +### 2. `fit/Likelihood.py` — collapse the two copies of `T`, add an array setter + +`T` is implemented twice: `_compute_T_chunk` reading `self._h5` +([Likelihood.py:1342](fit/Likelihood.py#L1342)) and again inline in observation +mode reading `byc` ([Likelihood.py:1748-1773](fit/Likelihood.py#L1748-L1773)). + +- Extract `_compute_T_from_columns(rid, columns_by_class, cA_per_class, + nuA_per_group, ln_bias, rate_shift, start, stop)`; `_compute_T_chunk` calls it + with the cache, observation mode calls it with `byc`. One implementation for + Asimov, observation and toy generation, and it makes the surrogate-route + injection in truth mode nearly free. + **The merge must keep `_compute_T_chunk`'s tolerance of `R_slice.shape[1] == 0`** + ([Likelihood.py:1358-1360](fit/Likelihood.py#L1358-L1360)); the inline + observation copy raises unconditionally on a dim mismatch + ([Likelihood.py:1752](fit/Likelihood.py#L1752)), which would break + `rate_shift`-only and POI-less classes. +- Move `_eval_region_surrogates` **up** from the subclass `N2LLExtensions` + ([N2LLExtensions.py:470](fit/N2LLExtensions.py#L470); `class + N2LLExtensions(N2LL)` at [N2LLExtensions.py:12](fit/N2LLExtensions.py#L12)) + into the base `N2LL`, which is what + [Likelihood.py:2556](fit/Likelihood.py#L2556) instantiates — that is why the + move is necessary at all. Rewrite `setObservation` to call it, deleting the + duplicated inline block at + [Likelihood.py:1551-1616](fit/Likelihood.py#L1551-L1616). + **The merged version must take `setObservation`'s permissive branch for a + missing BIT predictor** — the `(N,0)` array at + [Likelihood.py:1598](fit/Likelihood.py#L1598), not the hard raise at + [N2LLExtensions.py:525](fit/N2LLExtensions.py#L525). Two cases need it: the + `Eta_unbinned` `TTLep_pow_2016` class has no `POI` block at all, and the + `rate_shift` POI blocks that *are* present + (`Eta_unbinned_2016.yaml:160,299`) carry no `predictor`, so + `poi.get('predictor')` is `None` either way. + **What this costs:** `evaluate_ratio` currently fails loudly on a missing BIT + predictor via that raise; after the merge it will silently return + `c_dot_R == 0` instead. That is correct for the observation and toy paths, but + it is a real weakening for `evaluate_ratio` — note it in the docstring rather + than letting a future reader discover it. +- Add `setObservationArrays(unbinned_blocks=None, binned_counts=None)`: + populates `_obs_unbinned` / `_obs_binned_counts` from prepared arrays and + flips the mode flags exactly as + [Likelihood.py:1494-1507](fit/Likelihood.py#L1494-L1507). `setObservation` is + refactored to end in a call to it, so one place owns the mode switch. + Validates known region ids and that `len(w)` matches every `by_class` column. + Named for what it takes — already-evaluated arrays, as against the loaders + `setObservation` takes — rather than after the `block` local at + [Likelihood.py:1710](fit/Likelihood.py#L1710), which is real precedent but + undocumented jargon at a public signature. +- Add `setToy(toy, hypothesis)`: calls `setObservationArrays` and applies the + toy's thrown constraint centres to `hypothesis`. + +### 3. `common/derivative_providers.py` — shared provider factory + +`PDFDerivativeProvider(job["pdf"])` / `EFTDerivativeProvider(job["eft"].get("parameters", []))` +is constructed identically at four call sites: +[pdf_calibration.py:36](ML/Calibration/pdf_calibration.py#L36), +[eft_calibration.py:35](ML/Calibration/eft_calibration.py#L35), +[pdf_modification_plot.py:47](plot/bit/pdf_modification_plot.py#L47), +[eft_modification_plot.py:46](plot/bit/eft_modification_plot.py#L46). The toy +generator would be a fifth. Four instances warrants the helper: add +`build_derivative_provider(job)` dispatching on the presence of a `pdf` or `eft` +block, and change all four existing call sites over (no compatibility shim). + +### 4. `fit/ToyGenerator.py` (new) + +Pure generation — no plotting, no fit driving. + +No `ToyDataset` loader class: every generator below emits `_obs_unbinned`-shaped +arrays consumed by `setObservationArrays`, so nothing ever needs the loader +contract. (A loader would also have to carry `n_split`, which +`data/toy_data.py:_ArrayLoader` does not expose, and would hit +`setObservation`'s `ignore_weights=True` default and overwrite the +multiplicities with ones.) + +- `class ProcessSource` — the only normative definition is in §Multi-process + above. +No `Toy` class either: a toy is a plain dict, with `save_toy(path, toy)` and +`load_toy(path, n2ll)` as functions. There is no state outliving a call and no +behaviour beyond serialization, so per CLAUDE.md this is a function until proven +otherwise. The dict holds per-region unbinned blocks, binned counts, thrown +`constraint_centers`, generating hypothesis values, seed, `source` and +per-source diagnostics. + +- `generate_unbinned_toy_from_cache(n2ll, region_id, hypothesis, rng)` — reads + `n2ll._h5[(rid, cid)]`, computes `T` via `_compute_T_from_columns` (reusing + `_assemble_cA_per_class`, `_assemble_nuA_groups`, + `_assemble_rate_shift_per_class`, `_lnN_by_class`), throws multiplicities, and + returns **row indices plus multiplicities** for `n_i > 0` — not sliced + columns. The columns are verbatim slices of the cache, so copying them would + write roughly one cache per toy; `load_toy` rehydrates them by indexing + `n2ll._h5`. No `X`, which the cache does not store and `__call__` never uses. +- `generate_unbinned_toy_from_truth(n2ll, region_id, sources, rng, *, split="final_eval", hypothesis=None)` — + the multi-process algorithm above. A non-None `hypothesis` multiplies + `w_truth` by `1 + T`, which is the only injection route for nuisances and the + surrogate alternative to `coefficients` for POIs. +- `generate_binned_toy(n2ll, region_id, rng, *, hypothesis=None, sources=None, split="final_eval")` — + with `sources` (truth) it resolves the class list from `n2ll.binned`, + histograms truth weights into `_binned_unroll[rid]['edges']` and throws on + `lambda_truth`; without it (cache) it throws on `_compute_lambda_binned`. The + `sources` argument is required for the truth path — a signature taking only a + hypothesis could only ever produce the model lambda that §Binned rules out. +- `throw_constraint_centers(hypothesis, rng)` — for every `p.isPenalized`. +- `generate_toy(n2ll, seed, *, source, hypothesis=None, truth_sources=None, split="final_eval", throw_nuisances=False, allow_negative_weights=False)` — + top-level. Unbinned regions come from `n2ll.regions`; binned regions from + `n2ll.binned` (**not** `n2ll.regions`, which holds unbinned only — + [Likelihood.py:569](fit/Likelihood.py#L569) vs + [Likelihood.py:579](fit/Likelihood.py#L579); `_binned_regions_ids` is a list + of ids, not configs, so it cannot supply a class list). + +RNG: one `np.random.SeedSequence(seed)`, spawned per `(region_id, class_id)` +via a stable hash, so adding a region or process does not perturb another's +draws. + +Negative MC weights: Poisson is undefined for a negative mean. Crash by default +reporting the offending fraction; `allow_negative_weights=True` clips to zero +and logs the clipped weight sum, so the size of the approximation is always +visible. + +Persistence (HDF5, one file per toy), mirroring the `NN2LCache` layout. Cache +mode stores indices; truth mode has no cache to index into, so it stores the +evaluated columns: + + /meta attrs: seed, source, hypothesis json, config version + /unbinned//n (M,) multiplicities = observed weights + /unbinned//indices (M,) cache mode only — rows of _h5 + /unbinned//X (M,d) truth mode only + /unbinned//origin (M,) truth mode only, class index + /unbinned//by_class//g (M,) truth mode only + /unbinned//by_class//R (M,nA) truth mode only + /unbinned//by_class//Delta:: (M,nB) truth mode only + /binned//counts (Nflat,) + /constraint_centers (names + values) + +`/meta` carries no ids or dimensions. For cache mode the binding is structural — +you are indexing the live cache, so misalignment cannot occur. For truth mode +`setObservationArrays` already validates region ids and that `len(w)` matches +every `by_class` column (§2); extend that to assert `nA`/`nB` against the loaded +predictors, which checks the real runtime objects rather than a schema that has +to be kept in sync with them. + +### 5. The toy spec file + +A standalone YAML/JSON passed by path, **deliberately not a block in the +analysis config**. This mirrors the existing `--rotate` precedent +([Likelihood.py:2350](fit/Likelihood.py#L2350), loaded at +[Likelihood.py:2500](fit/Likelihood.py#L2500)): model-adjacent, versioned, +reusable across invocations, but outside the analysis YAML — so +`common/yaml_loader.py`, `combine_configs()` and `print_summary` are untouched +and no schema is locked in that other configs depend on. + +Governing split: **the spec defines the ensemble; the CLI selects which member +of it and where the result goes.** Everything that defines the pseudo-experiment +lives in the file; only the seed and output paths are flags. + +The spec declares **named points**, so one file covers a whole study and each +point's name lands in the output filename: + +```yaml +source: truth # or cache +split: final_eval # truth mode only; null for the whole dataset +throw_nuisances: true +allow_negative_weights: false +points: + - name: nominal + - name: c0_0p5 + injection: + SR_2018: + TTLep_pow_2018: + coefficients: {c0: 0.5} # exact PDF/EFT weights + - name: c0_0p5_bit + hypothesis: {c0: 0.5} # same point via the BIT surrogate + - name: dy_up + injection: + SR_2018: + DY_2018: {scale_factor: 1.2} +``` + +Per-point keys map onto `ProcessSource` fields (`sample_name`, `coefficients`, +`scale_factor`, `weight_branches`, `weight_function`) plus a +`hypothesis` block for surrogate-route injection. `weight_function` is a dotted +import path resolved with the `importlib` pattern already used for +`defaults.module_samples` ([Likelihood.py:2481](fit/Likelihood.py#L2481)). + +The presence of a `hypothesis` block **is** the switch for multiplying by +`1 + T`, so no separate `apply_model_hypothesis` flag is needed. For nuisances +it is the only injection route; for POIs it is the surrogate alternative to +`coefficients`. + +### The two POI injection routes + +`c0_0p5` and `c0_0p5_bit` above inject the same value by different routes, and +both are *fitted* with the BIT: + +- `c0_0p5` weights each event by the exact parametrization evaluated on its + **generator-level** info (`Generator_x1/x2/id1/id2/scalePDF`). No network + makes the data. +- `c0_0p5_bit` weights by `w0 * (1 + T(c0=0.5))`, with `T` evaluated on the + event's **reconstructed features**. + +`c0_0p5_bit` is therefore the control — generated and fitted by the same object, +unbiased by construction — and the difference in mean fitted `c0` is the BIT +mismodelling. + +The exact weight depends on generator quantities absent from the feature vector, +so the BIT cannot reproduce it event by event and should not: the most it can +learn is `E[w(c)/w(0) | x_reco]`, which is exactly what the likelihood needs, +since only the feature-space density ratio enters. The truth toy's intensity is +`rho_0(x) * R_true(x)` and the BIT toy's is `rho_0(x) * R_BIT(x)`, so **a perfect +BIT makes the two densities identical** — which is what makes this a clean null +test. + +They are not equally noisy, however. The truth toy carries extra variance from +the spread of `w_exact / w_0` at fixed `x`, a real MC fluctuation present even +with a perfect BIT. **The comparison statistic is the mean fitted POI, not the +width**; reading a width difference as mismodelling would be an error. + +Both points must live under the same `source: truth` and `split`, so they run +over the same events and the difference is apples-to-apples. A cache-mode +control would use all events and break that. + +Injecting POIs through both routes at once (`coefficients` and `hypothesis` +naming the same POI) is rejected — they exist precisely so their results can be +compared, as in verification item 9. + +### 6. Entry points + +**`fit/ToyGenerator.py.__main__` — standalone generation.** Nothing in +generation needs the fitter; it needs the config, loaded surrogates, +`build_cache()` and `prepare_runtime()`, i.e. everything +`Likelihood.py.__main__` does before it reaches `run_autograd_fit`. Generating +separately is worth more than convenience: + +- **Generate once, fit many ways.** Comparing `--no_syst` against full, or + different `--rotate` bases, or POI truncations, is far less noisy on + *identical* pseudo-data; regenerating per variant compares two ensembles. +- **Split the resource profile.** Generation is memory-heavy (whole cache in + RAM) and done once; fitting is CPU-heavy and fans out. +- **Inspect before spending.** Check per-source diagnostics and plot toy + distributions before committing hundreds of fits. + +``` +python fit/ToyGenerator.py --toySpec PATH --toyPoint NAME \ + --seeds 0-499 --outputDir DIR +``` + +Files are written as `/_toy.h5`, the same convention +`--toyFile` expects, so a Slurm array task computes its own path from +`$SLURM_ARRAY_TASK_ID` without an index file. + +**`fit/Likelihood.py.__main__` gains exactly one flag: `--toyFile PATH`.** It +loads a toy and never generates one. Generation lives solely in +`ToyGenerator.py`, so the fit entry point needs no spec parsing, no point +selection, no seeding, and — importantly — no `importlib` resolution of a +user-supplied dotted path. Two flags is one too many here: once standalone +generation exists for the reasons above, a second generating entry point is a +duplicate with no failure of its own to justify it, and it would contradict this +plan's own "no toy-fit driver" scope line. + +Rejecting `--data` and `--asimov` alongside it is not sufficient: nominal Asimov +is the *fall-through default* (`elif args.asimov is None:` at +[Likelihood.py:2609](fit/Likelihood.py#L2609)), so the toy branch must be +inserted **before** it, not merely guarded against an explicit flag. + +Applied via `n2ll.setToy(toy, hyp)` where `--data` calls `setObservation` +([Likelihood.py:2603](fit/Likelihood.py#L2603)). The point name and seed are +read from the toy's `/meta`, and the output suffix gains `_{POINT}_toy{SEED}` so +`serialize_result` writes one JSON per toy without collision. + +## Verification + +Run from the repo root against `configs/unbinned_v6/unbinned_2018.yaml`, and +`configs/Eta_unbinned/Eta_unbinned_2016.yaml` for the multi-process path. + +1. **No regression** — rerun a known Asimov fit; `fval` and covariance must be + unchanged, since constraint centres default to 0 and the `T` refactor is + pure extraction: + `python fit/Likelihood.py configs/unbinned_v6/unbinned_2018.yaml --overwrite fit` +2. **Deterministic injection check** — build a "toy" whose multiplicities are + the raw `w0` instead of a Poisson draw, inject via `setObservationArrays`, + and confirm `n2ll(hyp)` **equals** the nominal Asimov value for several + hypotheses. Equality is exact, not up to an offset: with `w == w0` the + observation branch `-2(sum w log1p(T) - sum w0 T)` + ([Likelihood.py:1742](fit/Likelihood.py#L1742), + [Likelihood.py:1775](fit/Likelihood.py#L1775)) is algebraically the Asimov + branch `sum w0 (log1p(T) - T)` + ([Likelihood.py:1829](fit/Likelihood.py#L1829)). Validates the whole cache + path with no random numbers involved. +3. **`T` refactor equivalence** — capture `_compute_T_chunk` output **before** + the refactor and assert `_compute_T_from_columns` reproduces it elementwise. + Comparing the two afterwards is vacuous, since `_compute_T_chunk` will + delegate to it. +4. **Yield closure** — over ~200 cache-mode toys at nominal, the mean of + `sum n_i` must equal `sum w0_i * (1 + T)` within its Poisson error per + region; likewise binned counts against `lambda`. +5. **Pulls and coverage** — ~500 toys at a point carrying + `hypothesis: {c0: 0.5}` with `throw_nuisances: true`; the pull + `(c_fit - c_true)/sigma_fit` should centre at 0 with unit width. Departures + of the width from 1 are the non-asymptotic effect this exists to expose. +6. **Truth mode closure** — `source: truth` with a point carrying no modifiers, + on the multi-class Eta config, is the surrogate-correct limit, so its pull + **mean** must agree with cache mode. **Compare means only, not widths.** + Truth mode defaults to `split: final_eval` (20% of events, weights rescaled + ~5x) while cache mode uses the full dataset, so the two ensembles have + different MC granularity by construction and their pull widths will differ + for reasons unrelated to the surrogates. A mean shift is a mismodelling + signal; a width difference is not. (Setting `split: null` for this one point + makes the widths comparable too, if that is wanted.) +7. **Multi-process injection** — a point with `EtaS_2016: {scale_factor: 1.2}`; + confirm the drawn yield of that source alone moves by 20% in the diagnostics + while the others are unchanged, and inspect the induced shift in the fitted + POIs. +8. **Exact-weight reconstruction check** — at `c = 0`, `w_coefficients` must + equal the nominal weight elementwise. Then verify the reconstruction against + an independent route: for a single-coefficient injection, compare + `w_coefficients` to a direct finite-difference of the parametrization, and + confirm the ratio `w_coefficients / w_nominal` matches the BIT's `1 + c_A R_A` + to within the BIT's known calibration accuracy (the quantity + `ML/Calibration/pdf_calibration.py` already measures). This is the check that + catches a mis-ordered or mis-normalized contraction. +9. **Surrogate bias measurement** — the headline use, and why the two injection + routes both exist. Fit the `c0_0p5` point (exact weights) and the + `c0_0p5_bit` point (surrogate weights) from the example spec above. The + difference in mean fitted `c0` is the BIT mismodelling bias, with MC + statistics common to both. +10. **Split normalization** — the failure mode here is silent. Generate the same + truth-mode point with `split: final_eval` and `split: null`; the mean drawn + yield per region must agree within its Poisson error, and both must match + `sum w0` over the cache. Confirm `f_weight` is computed from weights and + differs from the nominal bucket fraction `0.20` by a small but nonzero + amount, and check that the `n_i > 1` share is reported after rescaling and + is roughly five times worse on the split. +11. **Negative-weight reporting** — push `|c|` up until the truncated quadratic + turns negative and confirm the generator reports the negative fraction and + refuses to proceed unless `allow_negative_weights` is set. +12. **Persistence round-trip and binding** — generate with + `fit/ToyGenerator.py`, refit from the saved file, confirm identical `fval`. + Then load the same toy against a config with a different `version` or a + retrained BIT and confirm `load_toy` crashes rather than misaligning columns + silently. + +## Implementation order + +Work in a **git worktree**, one commit per piece of functionality: + +| # | Commit | Contents | Gate | +|---|---|---|---| +| 1 | provider factory | Change 3: `build_derivative_provider(job)` + the four call sites | existing calibration/plot scripts still run | +| 2 | `T` and surrogate-eval refactor | Change 2: `_compute_T_from_columns`, `_eval_region_surrogates` moved to base `N2LL`, `setObservationArrays`, `setToy` | **verification 1-3, exact equality** | +| 3 | constraint centers | Change 1: `fit/Modeling.py` | verification 1 again (centres default to 0, so `fval` unchanged) | +| 4 | the generator | Change 4: `fit/ToyGenerator.py`; Change 5: spec file; Change 6: `--toyFile` | verification 4-12 | + +Commits 1 and 3 are independent of everything else and can land in any order. +Commit 2 is the risky one and is gated hardest. Commit 4 is the only one that +adds new behaviour rather than moving existing behaviour. + + +**Do Changes 2 first, as its own commit, and run verification 1-3 before +touching anything else.** The `T` refactor and the `_eval_region_surrogates` +move are the highest-risk, lowest-creativity part of this work: they sit next to +numba-jitted code (`_weighted_sum_log1p_minus_x`, +[Likelihood.py:1829](fit/Likelihood.py#L1829)) and must preserve numerics +exactly. Verification 2 is an exact-equality check, so it either passes or it +does not. Everything else builds on that foundation and is far cheaper to debug +once it is known good. + +**If verification 2 fails and the cause is not obvious within a couple of +attempts, stop and escalate rather than patch.** Silent numerical wrongness is +the failure mode throughout this feature — a toy that is subtly wrong still +produces plausible pulls, plausible coverage and plausible plots. There is no +downstream check that will catch it. Treat a failing exact-equality assertion as +a hard stop, not as a tolerance to loosen. + +## Rejected alternatives + +Recorded separately in `toy-generation-design-decisions.md`, alongside this +file. Those decisions are settled — read that file before reopening any of +them. Rejections that bear on a specific section are also noted inline where +they apply, so implementing from this file alone will not reintroduce one. diff --git a/user/ricardo/claude/plans/toy-generation-status.md b/user/ricardo/claude/plans/toy-generation-status.md new file mode 100644 index 00000000..8daf963e --- /dev/null +++ b/user/ricardo/claude/plans/toy-generation-status.md @@ -0,0 +1,313 @@ +# Toy generation: implementation status + +Living status file for finishing the work in +`toy-generation-for-pseudo-experiments.md` (design decisions in +`toy-generation-design-decisions.md`). Update this as you debug/extend on your +end so the history doesn't only live in chat. + +Branch: `toy-generation-for-pseudo-experiments`, worktree at +`/users/ricardo.barrue/nsbi_gluon_pdf/GOLLUM-toygen`. + +## Landed (2026-07-29) + +- `f265664` -- `--plot` flag on `fit/ToyGenerator.py`: after generating a + point's toys, writes one figure per (unbinned region, default_feature), + overlaying each requested seed's weighted distribution (binning from + `data/plot_options`). New `plot/toys/toy_diagnostic_plots.py`; new + `fit.ToyGenerator.materialize_cache_region_features` (re-streams a region's + Asimov samples once, at the union of requested seeds' cache indices, to + recover raw kinematics for cache-mode toys, which otherwise only store + surrogate columns + row indices -- costly by design, opt-in only). Truth-mode + toys already carry X, but only for the region's own surrogate feature union; + default_features outside that union are skipped with a logged warning. + **Real-run verified**: `unbinned_2018.yaml` + `toys_cache_example.yaml`, + seeds 0-2, nominal point -- ran end to end, no crashes, exactly the 16 + expected `TOP_KINEMATICS + LEPTON_KINEMATICS + ASYMMETRY` features plotted + for `SR_2018` (32 files: png+pdf). One bug found and fixed along the way: + the plotting module originally imported `data.colors.cmap_petroff10_mpl` + (matching `plot/bit/modification_plotter.py`'s convention) for the toy-seed + color palette, but that module imports ROOT, and ROOT's cling JIT segfaulted + when initialized a second time in-process after `fit/ToyGenerator.py` had + already streamed samples via PyROOT/RDataLoader. Fixed by inlining the same + 10-color hex list directly instead of importing `data.colors`, keeping this + purely-matplotlib diagnostic script ROOT-free. (Separately, CERN EOS sync + via `common.syncer` fails in this environment -- no valid Kerberos ticket -- + unrelated to the plotting code itself.) + +- `f336ff1` -- full Asimov-from-cache reference band added to the diagnostic + plots: `materialize_cache_region_features` generalized into + `materialize_cache_region_diagnostics`, which accumulates a weighted + (`w0*(1+T)` at the toys' generation hypothesis, over the WHOLE cached MC) + histogram in the same re-streaming pass that recovers cache-mode toys' raw + kinematics -- no extra pass needed. `_compute_cache_intensity` factored out + of `generate_unbinned_toy_from_cache`, shared with the new diagnostics path. + Always computed regardless of toy source (cache/truth), since the surrogate + cache is always built. **Real-run verified** against + `unbinned_v7/unbinned_2016.yaml` (`SR_2016`, nominal point, seeds 0-2): grey + Asimov band tracks the toy seeds' mean, with visibly larger Poisson scatter + in the low-statistics tail (e.g. `tr_ttbar_mass` above ~2500 GeV), as + expected. + - **Real environment finding (not fixed, not code-related)**: two + consecutive real `--plot` runs (2018 retry, then 2016) crashed with the + same ROOT/cling segfault as the earlier `data.colors` issue -- but this + time triggered by `common.syncer`'s own module-level `import ROOT` + (it monkey-patches `TCanvas.Print`), not by anything in this branch's + code. Confirmed by re-running the plotting step directly against the + already-saved toy files with `syncer` excluded: completed cleanly, no + crash, correct-looking plots. You explicitly decided to keep `syncer` as + invoked (matches house convention; you'll refresh Kerberos manually for + EOS uploads) and treat the crash as pre-existing environment flakiness in + this conda env's cling JIT (mismatched compiler/std headers, see the + repeated `cannot extract standard library include paths` / `assert.h not + found` warnings that appear on every ROOT touch here, crash or not) -- + not something to work around in this branch. If it becomes a recurring + annoyance for the real `--plot` output path (`www/toys/...`), the isolated + fix is dropping `syncer` from `ToyGenerator.py`'s `--plot` branch alone + (this diagnostic script has no real need for EOS auto-sync), not touching + `common/syncer.py` itself. + - Example plots (grey Asimov band + 3 toy seeds, `SR_2016`) are sitting in + `/users/ricardo.barrue/.claude/jobs/cb7ec444/tmp/toys_plot_test_2016_plots_only/SR_2016/` + -- job-scratch, not durable; you said you'd copy them yourself. + +## Landed (2026-07-28) + +Four commits, in plan order: + +1. `44705f9` -- `build_derivative_provider(job)` factory in + `common/derivative_providers.py`, four call sites switched over + (`pdf_calibration.py`, `eft_calibration.py`, `pdf_modification_plot.py`, + `eft_modification_plot.py`). +2. `66b67d7` -- the risky one: `_compute_T_from_columns` unifies + `_compute_T_chunk` (cache) and `setObservation`'s inline by_class loop; + `_eval_region_surrogates` moved from `N2LLExtensions` to base `N2LL`, now + permissive on a missing BIT predictor (documented as a real weakening for + `evaluate_ratio`); `setObservationArrays(unbinned_blocks, binned_counts)` + added, `setObservation` refactored to end by calling it. +3. `2029f58` -- `ModelParameter.constraint_center`, + `Hypothesis.penalty()` centred, `Hypothesis.set_constraint_centers(mapping)`, + `Rotated` mirrors `constraint_center` on nuisance passthrough. +4. `7439354` -- `fit/ToyGenerator.py` (generation), `N2LL.setToy` + + `--toyFile` in `fit/Likelihood.py`, `_uid_split_interval` generalized in + `ML/Calibration/calibration_runner.py`, two example specs + (`configs/unbinned_v6/toys_example.yaml`, + `configs/Eta_unbinned/toys_example.yaml`). + +## What's actually verified vs. not + +No trained surrogates or MC caches exist anywhere in this environment +(`common/yaml_loader.py` reported 0/64 ready artifacts even on the unmodified +`sbi-eft-CMS` branch), so nothing here has run against a real config yet. +Everything below was checked with a synthetic in-memory `N2LL` (fake +classifier/BIT/PNN predictors, fabricated cache arrays) instead: + +- **Verified, bit-for-bit / exact:** + - `_compute_T_chunk` / `_compute_T_from_columns` reproduce the pre-refactor + reference elementwise (verification item 3). + - Deterministic-injection identity: observation branch == Asimov branch when + `w == w0` (verification item 2), to float roundoff (~1e-14). + - `ModelParameter.constraint_center` default 0 leaves `penalty()`/`fval` + unchanged. + - `save_toy` -> `load_toy` -> `setToy` round trip gives a bit-identical + `n2ll(hyp)` for a cache-mode toy. + - Exact PDF/EFT weight reconstruction formula + (`deriv_w[:,0] + deriv_w[:,1:] @ expand_pois_linear_quadratic(...)`) + checked at c=0 (reduces to nominal) and against an independent + finite-difference derivative (verification item 8's core math). +- **Verified, statistically (Poisson closure over ~300 toys):** + - Cache-mode unbinned yield closure (verification item 4, unbinned half). + - Cache-mode binned yield closure, via a monkeypatched `_compute_lambda_binned` + (verification item 4, binned half) -- real ICH/ICPH wiring not exercised. + - `throw_constraint_centers` mean/width. +- **Implemented but NOT exercised at all:** + - The entire truth-mode ROOT-streaming path in + `_materialize_truth_weights`: `RDataLoader.clone()`/`setFeatures`, + `UIDSplitter.mask_from_np`, `f_weight` renormalization, + `ProcessSource.weight_branches`/`weight_function`, negative-weight + reporting, per-source diagnostics, multi-process concatenation + (`generate_unbinned_toy_from_truth`, `generate_binned_toy` truth branch). + - `_class_bit_job` / `n2ll._toy_jobs_by_id` / `n2ll._toy_splitting_defaults` + wiring -- the pragmatic choice (not literally specified in the plan) of + attaching these to `n2ll` from `ToyGenerator.py.__main__` rather than + threading `cfg` through every generator call. Worth a second look once you + can actually run it. + - `fit/ToyGenerator.py.__main__` end to end (config loading, `build_cache`, + `prepare_runtime`, spec parsing, `--seeds` range parsing). + - `--toyFile` in `fit/Likelihood.py.__main__` (loads a toy, calls `setToy`, + suffixes the output path with `_{point}_toy{seed}`). + - Verification items 1, 5, 6, 7, 9, 10, 11, 12 from the plan -- all need a + real config with trained surrogates. + +## Suggested first real-run order + +1. **Verification 1** (no-regression Asimov fit) on + `configs/unbinned_v6/unbinned_2018.yaml` once its surrogates are trained -- + this is the cheapest way to catch anything the synthetic fixture couldn't, + since it exercises the real `_prepare_structure`/`prepare_runtime` path the + synthetic tests bypassed. +2. Generate a handful of cache-mode toys (`source: cache`, small `--seeds` + range) and check verification 4 (yield closure) for real -- fast, no ROOT + streaming involved, isolates whether the cache-mode generator is fine + before touching truth mode. +3. Truth mode: start with the `nominal` point in + `configs/unbinned_v6/toys_example.yaml` (no modifiers) at `--seeds 0-2` and + just confirm it runs and produces sane diagnostics -- this is the first + real test of the whole `RDataLoader`/UID-split/`_eval_region_surrogates` + chain. Watch for `f_weight` being close to the `final_eval` nominal bucket + fraction (~0.2) but not exactly it (verification item 10's "small but + nonzero" check). +4. Then `c0_0p5` / `c0_0p5_bit` (verification 8-9) and the multi-process + config (verification 6-7). + +## Log + +- 2026-07-28: initial implementation landed (see "Landed" above). Not yet run + against a real config. +- 2026-07-28: **verification 1 (no-regression Asimov fit), run for real.** + Config: `configs/unbinned_v7/unbinned_2018.yaml` (multi-process SR_2018: + TTLep_pow_2018, SingleTop_2018, DrellYan_LO_HTbinned_2018, TTSemi_pow_2018; + 73/75 artifacts ready, only the two TFMC classifier jobs missing -> classifier + degrades to `g=1` for all classes, harmless for this bit-identity check). + Method: a temporary detached-HEAD worktree at `ed8532d` (the commit right + before any of these 4 commits) at `/users/ricardo.barrue/nsbi_gluon_pdf/GOLLUM-ref`, + pointed at the same `common/user.py`/`common/directories.py` paths as the + main worktree, ran `python fit/Likelihood.py configs/unbinned_v7/unbinned_2018.yaml + --overwrite fit --verbosity 0` there (pre-refactor baseline) and in + `GOLLUM-toygen` (post-refactor), both against the *same* pre-built cache + (`NN2LCache/unbinned_2018/unbinned_v7`), then diffed the two `_fit.json` + outputs. Result: **fval, edm, niter, all 42 parameter values, all 42 + parameter errors, and the full 42x42 covariance and correlation matrices are + bit-for-bit identical (max abs diff = 0.0 everywhere).** `fval=0` both times, + as expected for a nominal Asimov self-fit. + Two environment issues found and fixed along the way, unrelated to the toy + generation code itself: (1) the worktree had briefly been nested inside + `GOLLUM/`, which broke every `sys.path.insert(0, '..')` in the codebase + (CWD-relative, not script-relative) -- moved back to being a sibling + directory; (2) `fit/Likelihood.py.__main__` unconditionally constructs a + `Data_` loader even in Asimov mode (pre-existing, confirmed via `git + diff ed8532d HEAD -- fit/Likelihood.py` that this code is untouched) -- + needed `common/directories.py`'s `SAMPLES_RUNII_BASE_DIRECTORY` pointed at + `v2-3-2_nJ2p_nB2p_2l` (has `Data_nominal.root` for all four eras) rather than + `v2-3_nJ2p_nB2p_trvalid` (only has a `2016/` subfolder, no Data at all). + The `GOLLUM-ref` worktree is still on disk in case it's useful for further + comparisons -- remove with `git worktree remove` (from `GOLLUM` or + `GOLLUM-toygen`) plus deleting the directory once done with it. + Remaining from "Suggested first real-run order": cache-mode toy generation + (step 2) and truth-mode (steps 3-4) still untested against real data. +- 2026-07-28: **verification 4 (cache-mode yield closure) and the + `--toyFile`/`setToy`/`load_toy` wiring, run for real**, on the same + `configs/unbinned_v7/unbinned_2018.yaml`. + - New spec `configs/unbinned_v7/toys_cache_example.yaml` (`source: cache`, + `throw_nuisances: true`, points `nominal` and `c0_0p5`). + - Bug found and fixed: `generate_unbinned_toy_from_cache` hard-crashed + unconditionally on **any** negative intensity, but 0.33% of events in every + class of this cache have negative `w0` (routine NLO-generator negative + weights -- confirmed directly from the cache: `TTLep_pow_2018.h5['w0']` + has 25683/7760846 negative entries, min -0.11), which triggers even at the + nominal hypothesis where `T` is identically 0 everywhere. Fixed by + extending the already-designed `allow_negative_weights` escape hatch (so + far only wired for truth mode) to `generate_unbinned_toy_from_cache` and + `generate_binned_toy`'s cache branch too, and threading it through + `generate_toy`'s cache-mode calls. Same crash-by-default / + clip-when-allowed semantics as truth mode, just applied to the mundane + "negative MC weight" cause as well as the "T < -1" cause. + - Also found and fixed: `fit/ToyGenerator.py` imported + `common.derivative_providers.build_derivative_provider` at module level, + which eagerly imports `data/samples_eft.py` (constructs RDataLoaders for + every declared EFT sample at import time) -- so cache-mode generation + depended on unrelated EFT sample data being reachable just to *import* the + module. Moved the import inline into `_materialize_truth_weights` (its + only call site, truth-mode only); `import fit.ToyGenerator` now succeeds + standalone with no EFT stub needed. + - Generated 5 cache-mode toys (`nominal` point, seeds 0-4, + `allow_negative_weights: true`) via + `python fit/ToyGenerator.py configs/unbinned_v7/unbinned_2018.yaml --toySpec + configs/unbinned_v7/toys_cache_example.yaml --toyPoint nominal --seeds 0-4 + --outputDir ...`. Metadata round-trips correctly (seed, source, point, 36 + thrown constraint centers). Yield closure: mean drawn yield over 5 toys = + 234660.8 vs expected `sum(clip(w0,0,None))` = 234907.5, diff = 1.14 sigma + (Poisson error) -- consistent. + - `--toyFile` end-to-end: `load_toy` + `setToy` + `n2ll(hyp)` all confirmed + working via a lightweight direct check (no minimizer): finite `-2lnL` at + the nominal hypothesis (35.78, nonzero as expected for a real thrown + fluctuation, unlike the trivial Asimov self-fit), rises sensibly away from + it (1239 at `c0=0.3`, 44.3 at `nu_pu=0.5`), and repeat-evaluation at the + same hypothesis is bit-identical (no hidden state mutation across + `setToy` calls). + - The **actual fit** (`python fit/Likelihood.py ... --toyFile + nominal_toy0.h5 --overwrite fit`) was launched in the background + (2018 toy) and ran for ~3 hours burning 1625% CPU with zero visible log + output (Python fully block-buffers stdout when redirected to a file, so + `prepare_runtime()`'s prints were sitting in an unflushed buffer the whole + time -- not evidence it was stuck). **Killed manually** once the 2016 run + below revealed why it was never going to finish cleanly. See "NaN at first + optimizer step" below. + - **NaN at first optimizer step -- real finding, pre-existing code, not a + toy-generation bug.** Reran the same kind of fit against a fresh 2016 toy + (`configs/unbinned_v7/unbinned_2016.yaml`, same spec) with `--verbosity 2` + and `PYTHONUNBUFFERED=1` for visibility. Two bugs surfaced: + 1. `ModelParameter.__repr__` (`fit/Modeling.py:53`) crashed with + `TypeError: unsupported format string passed to ArrayBox.__format__` + -- `f"{self.val:.6e}"` fails when `self.val` is an autograd `ArrayBox` + (during gradient tracing) rather than a plain float. Pre-existing, + apparently never hit before since `--verbosity 2` (the only level that + calls `Hypothesis.print()` from inside `grad(fcn)`) had never been used + against a real fit. **Fixed and committed (`48c81f0`)**: wrap in + `float(getval(self.val))` before formatting. + 2. After that fix, the fit reran and got further: eval 1 (`f=36.02`, pure + penalty term -- `T` is identically 0 at `x=0` for every event by + construction, so `total_unbinned(x=0)=0` exactly and the whole value + comes from the 36 thrown constraint centers' `sum(offset^2)`), eval 2 + (autograd's internal forward retrace at the same `x=0`, same value, + confirms it), **eval 3 = `f=nan`, crashing with `RuntimeError: NaN + likelihood!`**. Diagnosis: eval 3 is L-BFGS-B's first real step (not a + per-parameter probe -- autograd computes the whole 42-dim gradient in + one reverse-mode pass), and that step pushes `T < -1` for at least one + event, making `log1p(T)` NaN. Likely driven by a large gradient + component from an outlier event (plausibly among the negative-weight / + extreme-kinematics ones found earlier) combined with L-BFGS-B taking an + unscaled first step from `x=0` (the `steps` dict `run_autograd_fit` + computes is never actually passed to `scipy.optimize.minimize` -- + looks like dead code). **This is not toy-generation-specific**: `T(x=0)` + is identically 0 for *any* observed dataset regardless of content + (real data included), so this fragility is inherent to + `run_autograd_fit`'s unscaled first step from the origin, just never + exercised before because the only fits run so far were nominal-Asimov + self-fits that never leave `x=0` (`niter=0` in verification 1). + - Per your instruction: **not fixing this now** (2026-07-28) -- logging it + here as a known blocker for verification 12 (refit-converges), out of + scope for the toy-generation plan itself since it lives in + `run_autograd_fit`'s optimizer setup, not in anything this branch + touched. The 2018 background fit was killed rather than left running, + since it's the same code path against the same kind of toy and was very + likely going to hit the identical wall. + - To pick this up: the fix is probably either (a) guard `log1p`/clip `T` + at `-1+eps` before the log (cheap, but masks the real issue), or (b) + actually use the computed `steps` dict to scale the initial L-BFGS-B + step (`options={"eps": ...}` or a proper trust-region method), or + (c) identify and understand the specific outlier event(s) driving the + huge gradient component first. Reproduce with: + `PYTHONUNBUFFERED=1 python -u fit/Likelihood.py + configs/unbinned_v7/unbinned_2016.yaml --toyFile + /users/ricardo.barrue/.claude/jobs/cb7ec444/tmp/toys_cache_2016/nominal_toy0.h5 + --overwrite fit --verbosity 2` (toy file may need regenerating if the + job tmp dir was cleaned). + - **(2026-07-29) Observation from you**: this NaN does *not* happen when + fitting the same toy with `run_iminuit_fit` (Minuit/MIGRAD) -- only + `run_autograd_fit` (L-BFGS-B) hits it. This is consistent with the + diagnosis above and narrows it further: MIGRAD's first step is scaled + by its own internal error/step estimates rather than taking an + unscaled step from `x=0`, so it doesn't overshoot into `T < -1`. + Strengthens candidate fix (b) (use the computed `steps` dict to scale + L-BFGS-B's first step) over (a) (clip/guard, masks the issue) as the + more likely real fix, though (c) (identify the outlier event driving + the huge gradient component) is still worth doing first to confirm. + - Remaining: verification 12 (refit-converges) blocked on the above; truth- + mode generation (steps 3-4 of the suggested order) still untested against + real data. + - 2026-07-29: generation of toy in truth mode tested for c0 = 0.5. + - added examples for truth-mode toy generation from an alternative sample + or additional weights (e.g. QCD scales). The latter required adding the + additional weight branches to the list of observers in the toy generator. + - Minuit fits with truth-mode toys are still crashing. Will debug later. + +