"""Nested cohort validation with training-fold masks and paired query subjects. Starts from counts supplied by the user. It cannot undo prior global depth processing. Draw percentiles quantify reference-set variation, not population CIs. """ import hashlib import json from pathlib import Path import numpy as np import pandas as pd from sklearn.linear_model import LogisticRegression from sklearn.metrics import roc_auc_score from sklearn.preprocessing import StandardScaler from .estimator import CFREF def training_mask(short, long, train, min_median_total=100): """Learn feature eligibility from training rows only.""" mask = np.median((short + long)[train], axis=0) >= min_median_total if not mask.any(): raise ValueError('Training fold retained no features.') return mask def _metrics(scores, y, threshold): return dict(auc=float(roc_auc_score(y, scores)), specificity=float(np.mean(scores[y == 0] <= threshold)), sensitivity=float(np.mean(scores[y == 1] > threshold))) def _reference_draws(y, n_ref, draws, seed, min_query_controls): ctrl = np.flatnonzero(y == 0) if len(ctrl) < n_ref + min_query_controls: return rng = np.random.RandomState(seed) for draw in range(draws): ref = np.sort(rng.choice(ctrl, n_ref, replace=False)) query = np.setdiff1d(np.arange(len(y)), ref) yield draw, ref, query def _cfref_summary(model, X, y, names, n_ref, draws, seed, min_query_controls): spec = [] for _, ref, query in _reference_draws(y, n_ref, draws, seed, min_query_controls): scores = model.decision_function(X[query], X[ref], feature_names=names) spec.append(_metrics(scores, y[query], model.threshold_)['specificity']) if not spec: raise ValueError('Insufficient inner target controls for reference draws.') return float(np.median(spec)) def nested_lopo(short_counts, long_counts, y, cohorts, sample_ids, *, feature_names=None, outer_cohorts=None, seeds=(0, 1, 2), capacities=((64, 24), (128, 48), (256, 96)), episodes=3000, reference_sizes=(5, 10, 20, 30), draws=60, threshold_draws=40, inner_references=5, outer_threshold_references=20, target_specificity=0.95, min_median_total=100, min_query_controls=3, output_dir=None, progress=None, input_preprocessing='Unspecified upstream count processing'): """Return paired per-draw results and nested-selection audit tables. Parameters use the archived reference/threshold conventions by default. Sample IDs must be globally unique; all rows of one cohort are held out together. Each inner fold independently learns its count-coverage mask and scaler. Target labels are used ONLY for reference identification/evaluation, never for training, threshold construction or outer capacity selection. Methods: nested cfREF, source-locked logistic regression, locally recalibrated logistic regression. All use identical target query subjects for each draw. Outputs include draw-level paired differences; these are not independent patient replicates. No capacity is selected using outer performance. """ S = np.asarray(short_counts, dtype=float); L = np.asarray(long_counts, dtype=float) y = np.asarray(y); coh = np.asarray(cohorts); ids = np.asarray(sample_ids).astype(str) if S.ndim != 2 or S.shape != L.shape or min(S.shape) == 0: raise ValueError('Counts must be matching nonempty samples × bins arrays.') if not np.isfinite(S).all() or not np.isfinite(L).all() or (S < 0).any() or (L < 0).any(): raise ValueError('Counts must be finite and nonnegative.') if y.shape != (len(S),) or coh.shape != y.shape or ids.shape != y.shape: raise ValueError('Labels, cohorts and IDs must match sample count.') if not np.isin(y, [0, 1]).all() or len(np.unique(y)) != 2: raise ValueError('Binary labels with both classes are required.') if len(set(ids)) != len(ids) or any(not i.strip() for i in ids): raise ValueError('Subject IDs must be globally unique nonempty strings.') if not all(isinstance(c, (str, np.str_)) and c.strip() for c in coh): raise ValueError('Cohort IDs must be nonempty strings.') if len(np.unique(coh)) < 3: raise ValueError('Nested cohort validation requires at least three cohorts.') if not capacities or len(set(tuple(c) for c in capacities)) != len(capacities): raise ValueError('Supply unique candidate capacities.') if not seeds or len(set(seeds)) != len(seeds) or not all(isinstance(s,int) and 0 <= s < 2**32 for s in seeds): raise ValueError('Supply unique nonnegative integer seeds.') for v in [draws, threshold_draws, episodes, min_query_controls]: if not isinstance(v, int) or v < 1:raise ValueError('Draws, episodes and minimum controls must be positive integers.') if not reference_sizes or len(set(reference_sizes)) != len(reference_sizes) or not all(isinstance(v,int) and v>=2 for v in reference_sizes): raise ValueError('Reference sizes must be unique integers >=2.') if not np.isfinite(min_median_total) or min_median_total < 0:raise ValueError('Invalid coverage threshold.') names = [f'bin_{j}' for j in range(S.shape[1])] if feature_names is None else list(feature_names) if len(names) != S.shape[1] or len(set(names)) != len(names) or not all(isinstance(n,str) and n for n in names): raise ValueError('Feature names must be unique ordered strings matching bins.') outer = list(np.unique(coh)) if outer_cohorts is None else list(outer_cohorts) if len(set(outer)) != len(outer) or any(c not in coh for c in outer):raise ValueError('Invalid outer cohorts.') ratio = S / np.maximum(S + L, 1) result_rows=[]; inner_rows=[]; picks=[]; assignments=[]; masks=[]; skips=[] cfg=dict(seeds=list(seeds),capacities=list(capacities),episodes=episodes, reference_sizes=list(reference_sizes),draws=draws,threshold_draws=threshold_draws, inner_references=inner_references,outer_threshold_references=outer_threshold_references, target_specificity=target_specificity,min_median_total=min_median_total, min_query_controls=min_query_controls,input_preprocessing=input_preprocessing, quantile='linear; in-sample source controls',mask_policy='training rows of each inner/outer fold', selection='mean over inner cohorts of abs(median draw specificity - target); capacity tuple order breaks ties', patient_independent_confidence_intervals=False) def log(msg): if progress is not None:progress(msg) def save_partial(): if output_dir is not None: p=Path(output_dir);p.mkdir(parents=True,exist_ok=True) for name, rows in [('draw_metrics',result_rows),('inner_selection',inner_rows),('capacity_picks',picks),('skipped',skips)]: pd.DataFrame(rows).to_csv(p/(name+'.csv'),index=False) (p/'fold_masks.json').write_text(json.dumps(masks,indent=2)) (p/'protocol.json').write_text(json.dumps(cfg,indent=2)) for ho in outer: tr=coh!=ho;te=~tr; train_cohorts=np.unique(coh[tr]); yy=y[te] if (yy==0).sum()