diff --git a/.travis/test-lisa.sh b/.travis/test-lisa.sh index e821b5739..91021697f 100644 --- a/.travis/test-lisa.sh +++ b/.travis/test-lisa.sh @@ -15,4 +15,11 @@ fi MonteCarloMarginalizeCode/Code/test/test_lisa_helper_contract.py \ MonteCarloMarginalizeCode/Code/test/test_lisa_pseudo_pipe_contract.py \ MonteCarloMarginalizeCode/Code/test/test_lisa_pp_surface.py \ - MonteCarloMarginalizeCode/Code/test/test_lisa_synthetic_demo.py + MonteCarloMarginalizeCode/Code/test/test_lisa_synthetic_demo.py \ + MonteCarloMarginalizeCode/Code/test/test_lisa_fairdraw_weights.py \ + MonteCarloMarginalizeCode/Code/test/test_lisa_l0_rescue.py \ + MonteCarloMarginalizeCode/Code/test/test_lisa_sampler_plumbing.py \ + MonteCarloMarginalizeCode/Code/test/test_lisa_av_state.py \ + MonteCarloMarginalizeCode/Code/test/test_lisa_use_lnL_branches.py \ + MonteCarloMarginalizeCode/Code/test/test_lisa_portfolio_method_integrity.py \ + MonteCarloMarginalizeCode/Code/test/test_lisa_driver_drift.py diff --git a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode_lisa b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode_lisa index 3754e5452..aca690a95 100755 --- a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode_lisa +++ b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode_lisa @@ -306,7 +306,39 @@ integration_params.add_option("--internal-use-lnL",action='store_true',help="lik integration_params.add_option("--sampler-method",default="adaptive_cartesian_gpu",help="adaptive_cartesian|GMM|adaptive_cartesian_gpu") integration_params.add_option("--sampler-portfolio",default=None,action='append',type=str,help="comma-separated strings, matching sampler methods other than portfolio") integration_params.add_option("--sampler-portfolio-args",default=None, action='append', type=str, help='eval-able dictionaryo to be passed to that sampler') +# Portfolio freeze/allocation policy. Pure pass-through to the shared portfolio sampler. +# Definitions copied verbatim from bin/integrate_likelihood_extrinsic_batchmode; pinned by +# test_lisa_sampler_plumbing.py. +integration_params.add_option("--portfolio-adaptive-alloc",action='store_true',default=False,help="Portfolio: ENABLE (opt-in) adaptive-probe draw allocation -- concentrate draws on the best per-chunk-n_ess member. Good on strongly-correlated targets; NOT recommended for AV-favorable high-SNR events (it starves the slow-contracting AV workhorse). Off by default (legacy n_ess reweighting).") +integration_params.add_option("--portfolio-alloc-exponent",default=None,type=float,help="Portfolio: adaptive allocation ~ member_quality^exponent. Higher concentrates harder on the winner. Sampler default 1.0.") +integration_params.add_option("--portfolio-freeze-wt",default=None,type=float,help="Portfolio: a member whose balance weight is below this stops updating its proposal (subject to grace/revive/VARAHA-exemption). Sampler default 0.05.") +integration_params.add_option("--portfolio-grace-iters",default=None,type=int,help="Portfolio: never freeze ANY member during the first N integration chunks (let slow starters contract). Sampler default 25.") +integration_params.add_option("--portfolio-probe-period",default=None,type=int,help="Portfolio: round-robin probe one member at a raised draw share every N chunks (breaks the under-observation trap). 0 disables probing. Sampler default 4.") +integration_params.add_option("--portfolio-quality-signal",default=None,type=str,help="Portfolio adaptive allocation: which per-member quality signal to rank members by. 'global' (default) = marginal gain in POOLED n_eff per sample (credits weight mass, debits weight variance); 'credit' = q_mix-native MIS credit assignment, sum_i [frac_m q_m/q_mix]_i * w_i per drawn sample (credits a member for COVERING where the integrand is, even if it drew few samples there); 'ness' = legacy per-member Kish n_ess (scale-invariant, misranks a slow-contracting AV -- see DESIGN_portfolio_freeze_policy.md).") +integration_params.add_option("--portfolio-revive-period",default=None,type=int,help="Portfolio: every N chunks, update even a frozen member one step so it can recover. 0 disables. Sampler default 8.") +integration_params.add_option("--portfolio-varaha-can-freeze",action='store_true',default=False,help="Portfolio: DISABLE the VARAHA freeze-exemption, so VARAHA/AV members obey the grace/revive/weight freeze schedule like other members. Use only if a VARAHA member is a known-bad fit and you want to save its selfish-draw eval cycles.") +integration_params.add_option("--portfolio-varaha-max-frac",default=None,type=float,help="Portfolio: CAP the combined DRAW fraction of VARAHA/AV members (0/unset = no cap). Use WITH --portfolio-varaha-min-frac to constrain the VARAHA share to a BAND. Rationale: a floor alone stops the mixture degenerating to peaked-member-only (which strips q_mix of its broad backstop, so a missed mode goes uncovered and lnZ is silently low while n_eff looks GOOD), but the share can then run away the OTHER way to ~1 and the mixture degenerates to VARAHA-only instead. A band (e.g. 0.25/0.75) keeps q_mix genuinely mixed by construction. Unbiased either way (balance heuristic), so it costs at most draws, never correctness.") +integration_params.add_option("--portfolio-varaha-min-frac",default=None,type=float,help="Portfolio: reserve this combined DRAW fraction for VARAHA/AV members (0/unset = off). never-freeze keeps a VARAHA member UPDATING, but both allocation rules score by per-chunk n_ess, which sits at ~1 during VARAHA's slow cumulative contraction -- so a member that looks instantly good can take nearly the whole budget (measured on S250114ax post-#33: GMM took ~0.84 and the portfolio collapsed to n_eff ~2 vs ~100 for standalone AV). Unbiased for any allocation (q_mix); trades efficiency only.") +integration_params.add_option("--portfolio-varaha-never-freeze",action='store_true',default=False,help="Portfolio: VARAHA/AV members always update every chunk past their breakpoint (freeze-exempt). This is the sampler default; the flag is here for explicitness/pipe pass-through.") +integration_params.add_option("--portfolio-weight-clip",default=None,type=float,help="Portfolio: OPT-IN truncated importance sampling applied to the PROPOSAL-FIT INPUT ONLY. Caps the weights fed to member.update_sampling_prior (the GMM covariance fit) at tau = C*sqrt(n)*mean(w) (0/unset = off; C~1 is the standard Ionides choice), so one enormous weight cannot make that fit degenerate. The estimator (ln Z, n_eff), the n_ess report, and the allocation signal all use the TRUE unclipped weights, so they stay exactly unbiased and undistorted. Do NOT clip the estimator (measured on S250114ax: n_eff=100 2x faster than AV but ln Z biased -11.5 nats) or the n_ess report (clipping inflates the clipped member's n_ess and starves the AV workhorse). The withheld tail mass is tracked and reported as a diagnostic. NOTE: if huge weights come from q_mix UNDERFLOW (watch for the warning) they are a numerical artifact, not tail mass.") integration_params.add_option("--sampler-xpy",default=None,help="numpy|cupy if the adaptive_cartesian_gpu sampler is active, use that.") +# AV live-volume state, per-axis bin allocation, and the collapse gate. Copied verbatim +# from bin/integrate_likelihood_extrinsic_batchmode; pinned by test_lisa_av_state.py. +integration_params.add_option("--sampler-save-state",default=None,help="AV only: after integration, write the adapted live-volume state (.npz) for reuse by later instances/iterations. Point --sampler-load-state at the same file across a grid to warm-start each point from the previous one.") +integration_params.add_option("--sampler-load-state",default=None,help="AV only: load a saved live-volume state (.npz from --sampler-save-state) to warm-start this integration. Overrides --sampler-warmstart-samples.") +integration_params.add_option("--sampler-anisotropic-bins",action="store_true",help="AV only: give each extrinsic axis a DIFFERENT number of bins during contraction -- fine where the live points cluster tightly (phase/polarization/sky), coarse where they are broad (distance/inclination) -- instead of the default equal split. Keeps the same total bin budget, so the estimator is unchanged; helps AV wrap a correlated/degenerate posterior more tightly.") +optp.add_option("--reject-collapsed-live-volume",action='store_true',default=False, help="DROP an event whose adaptive-volume live volume degenerated (see the [AV COLLAPSE] report) instead of exporting it: the integration is treated as a failure, so no likelihood row, XML or posterior samples are written for it. Such a run's lnZ and samples describe a single mode of the integrand and are NOT a fair posterior draw, and nothing downstream can distinguish them from a converged export. Default off, because dropping the event silently THINS the posterior in an SNR-dependent way -- that was the pre-fix behaviour, when this case crashed. Left off, the event is exported but announces itself loudly and (with --mc-error-replicas>0) triggers replication. Turn it on when a contaminated point is worse than a missing one.") +# L0 auto-rescue. Ported from bin/integrate_likelihood_extrinsic_batchmode; defaults and help +# text kept IDENTICAL there and here on purpose -- see test_lisa_l0_rescue.py, which pins them. +integration_params.add_option("--sampler-warmstart-retry-neff",type=float,default=None,help="AV or portfolio (L0 auto-rescue): if a pass finishes below this n_eff (i.e. it stalled on a very sharp / high-amplitude peak), automatically re-run a second pass warm-started from THIS point's own highest-likelihood samples. Same-problem reuse in the sense that the seed provably contains the peak the cold pass found -- but NOT that every mode is represented, so the warm pass can be biased low if the seed missed one. The rescue still runs as before; its result is rejected in favour of the cold pass only on positive evidence of lost mass (see --sampler-l0-rescue-reject-dlnZ). A portfolio is unaffected: its GMM member carries a defensive component. Directly targets the high-SNR n_eff LOTTERY (a large fraction of independent runs collapse to n_eff~1 by contracting onto the wrong spot); the rescue re-seeds a collapsed run from the peak it did find. Recommended for high-SNR events; e.g. 5.") +integration_params.add_option("--sampler-l0-rescue-reject-dlnZ", type=float, default=3.0, help="Evidence threshold (nats) for rejecting the L0 rescue's warm pass: reject when the full-support cold pass reports lnZ this much HIGHER, which would indicate the seed missed mass. Larger = more permissive. DEFAULT RAISED 0.5 -> 3.0 ON MEASUREMENT (see test/expensive_before_merging/integrators/L0_REJECT_DLNZ_MEASUREMENT.md): across 160 known-lnZ passes the gate caught 0 of 55 genuinely truncated warm passes at EVERY threshold, while at 0.5 it binned 25% of GOOD portfolio warm passes. 0.5 was therefore strictly dominated -- it bought no detection and cost one good pass in four. 3.0 keeps a safety net for a genuinely large discrepancy at ~0% false-positive rate. This gate is NOT a working truncation detector; do not rely on it as one.") +integration_params.add_option("--sampler-l0-rescue-accept-truncated", action='store_true', default=False, help="Report the L0 rescue's warm pass even when it lands well below the full-support cold pass (see --sampler-l0-rescue-reject-dlnZ). Default OFF: on that evidence the cold result is kept instead, since the warm pass is confined to the seeded peak and may be missing a mode. The rescue itself still runs either way.") +integration_params.add_option("--sampler-l0-rescue-puff-scale", type='choice', choices=['fixed','auto'], default='auto', help="How wide to puff the L0 rescue's seed when it is rank-deficient in the adaptive dimensions. 'auto' (default) measures the posterior scale AND correlations from every finite lnL the collapsed pass already drew; 'fixed' uses --sampler-l0-rescue-puff-width-frac of each parameter's prior range, which is the historical behaviour and knows nothing about the posterior (which narrows as 1/rho). 'auto' falls back to 'fixed' when there are too few finite points to estimate a covariance.") +integration_params.add_option("--sampler-l0-rescue-puff-width-frac", type=float, default=0.005, help="Isotropic puff width for the L0 rescue's rank-deficient seed, as a fraction of each parameter's prior range. Used by --sampler-l0-rescue-puff-scale fixed, and as the 'auto' fallback. Default 0.005 = the historical hardcoded 1/200.") +integration_params.add_option("--sampler-l0-rescue-puff-factor", type=float, default=2.0, help="Multiply the L0 rescue's puff width by this factor. Default 2 is the measured optimum on a known-lnZ 6-D target (mean lnZ error +0.08 nats, ESS 52); BOTH tails are wrong, so do not treat wide as free -- x0.5 truncates (-8.5 nats), x6 biases high (+3.0) and costs efficiency, x12 is a cold start in all but name and re-collapses (-30).") +# Also consumed by the rescue (it is the lnL window build_warm_seed keeps), which is why it +# lands in this pass rather than with the sequential warm start it is named for. +integration_params.add_option("--sampler-sequential-warmstart-deltalnL",type=float,default=15.0,help="Keep previous-point samples within this lnL of the max as the warm seed for the next point. Default 15.") integration_params.add_option("--supplementary-likelihood-factor-code", default=None,type=str,help="Import a module (in your pythonpath!) containing a supplementary factor for the likelihood. Used to impose supplementary external priors of arbitrary complexity and external dependence (e.g., EM observations). EXPERTS-ONLY") integration_params.add_option("--supplementary-likelihood-factor-function", default=None,type=str,help="With above option, specifies the specific function used as an external prior. EXPERTS ONLY") integration_params.add_option("--supplementary-likelihood-factor-ini", default=None,type=str,help="With above option, specifies an ini file that is parsed (here) and passed to the preparation code, called when the module is first loaded, to configure the module. EXPERTS ONLY") @@ -790,6 +822,7 @@ params = {} sampler = mcsampler.MCSampler() xpy_asarray_already = functools.partial(xpy_default.asarray,dtype=np.float64) +use_gmm_member=False # set when a portfolio carries a GMM member (see the portfolio setup loop) if opts.sampler_method == "adaptive_cartesian_gpu": print(" ILE: {}".format(opts.sampler_method)) sampler = mcsamplerGPU.MCSampler() @@ -844,13 +877,32 @@ elif opts.sampler_method == "portfolio": for name in sampler_types: if name =='AV': sampler = mcsamplerAdaptiveVolume.MCSampler(n_chunk=opts.n_chunk) # enforce now, so provided for setup phase - if name =='GMM': + elif name =='GMM': sampler = mcsamplerEnsemble.MCSampler() - # following override means sampler_method is CHANGED, so THIS MUST BE LAST, and can't condition on portfolio - opts.sampler_method = 'GMM' # this will force the creation/parsing of GMM-specific arguments below, so they are properly passed - if name == "adaptive_cartesian_gpu": + # A GMM member needs the GMM-specific argument blocks below to run so its config is + # forwarded. This used to CLOBBER opts.sampler_method='GMM', which silently broke + # every downstream `sampler_method == "portfolio"` test -- most importantly the L0 + # auto-rescue gate, which then NEVER FIRED for a portfolio carrying a GMM member -- + # and made a portfolio take GMM-only branches (e.g. return_lnI). Flag it + # non-destructively instead: sampler_method stays 'portfolio', and the GMM blocks + # below key off `use_gmm_args` = standalone GMM OR a portfolio with a GMM member. + # Ported from bin/integrate_likelihood_extrinsic_batchmode, which fixed this. + use_gmm_member = True + elif name == "adaptive_cartesian_gpu" or name == 'AC': sampler = mcsamplerGPU.MCSampler() mcsampler = mcsamplerGPU # force use of routines in that file, for properly configured GPU-accelerated code as needed + elif name in mcsamplerPortfolio.known_pipelines: # everything else, including nflow + sampler = mcsamplerPortfolio.known_pipelines[name]() + else: + # No else clause here meant an unrecognized name left `sampler` bound to its + # previous value -- the plain MCSampler built before this chain, or, on the second + # and later iterations, the PREVIOUS member -- and appended it silently. The + # portfolio then ran with a member the user never asked for, and a typo in + # --sampler-portfolio produced a duplicate rather than an error. Ported from + # bin/integrate_likelihood_extrinsic_batchmode. (The chain above is now elif for + # the same reason: with plain `if`, a name matching no branch fell through every + # test and reused whatever `sampler` still held.) + raise Exception(" --sampler-portfolio: unknown member '{}'. Known: AV, GMM, AC/adaptive_cartesian_gpu, {}".format(name, sorted(mcsamplerPortfolio.known_pipelines))) print('PORTFOLIO: adding {} '.format(name)) # enable xpy for low level sampler as needed if hasattr(sampler, 'xpy'): @@ -1197,10 +1249,175 @@ if opts.sampler_method=="GMM" and opts.internal_use_lnL: if opts.sampler_method =="adaptive_cartesian_gpu" and opts.internal_use_lnL: return_lnL=True pinned_params.update({"use_lnL":True}) +if opts.sampler_method =="AV" and opts.internal_use_lnL: + # AV integrates in log space natively (integrate() is a thin wrapper over integrate_log); + # without this, --internal-use-lnL --sampler-method AV passed the ok_lnL_methods check but + # silently did nothing, so exp(lnL) overflowed at high SNR when no logarithm offset was set. + # + # PORTED FROM THE MAIN DRIVER, where this branch already exists. It was missing here, and + # the drift audit could not see it: a missing `if` branch is not a FUNC/OPTION/CONST/ATTR, + # so it produces no gap item. High-SNR is the LISA MBHB regime, which is exactly the case + # the main driver's comment describes. + return_lnL=True + pinned_params.update({"use_lnL":True}) if opts.sampler_method =="portfolio": return_lnL=True pinned_params.update({"use_lnL":True}) -if opts.sampler_method == "GMM": + +# What the sampler will actually STORE in _rvs['integrand'], derived from the pinned params +# above rather than from the CLI. This is not the same predicate as opts.internal_use_lnL: +# that option is accepted for adaptive_cartesian_gpu and portfolio too (see the branches +# directly above), which set use_lnL WITHOUT return_lnI and therefore still store linear L. +# Keying the weight helpers off the option would compute L + ln p - ln p_s for those, which +# is the failure the main driver documents at ln_weights_from_rvs. +rvs_integrand_is_lnL = bool(pinned_params.get("return_lnI", False)) + + +# --------------------------------------------------------------------------------------- +# Fair-draw weighting helpers. Ported from bin/integrate_likelihood_extrinsic_batchmode +# (PR #87); see test/expensive_before_merging/integrators/RVS_FAIRDRAW_AUDIT.md. +# +# WHY THESE ARE HERE, given this driver has no .dgrid/.dslice/proposal-breadcrumb exports +# (the three consumers whose double-weighting PR #87 actually fixed): this driver DOES set +# igrand_fairdraw_samples from --fairdraw-extrinsic-output, so its _rvs can be a fair draw, +# and every shared sampler already sets the provenance marker at its rebind. The marker was +# arriving here and nothing was reading it. The helpers are the correct thing for the next +# person to reach for, which is the whole argument of the audit's Recommendation 1. +# +# KEEP IN STEP WITH THE MAIN DRIVER. These are deliberate copies, not an import, because the +# two drivers are a deliberate fork; audit_lisa_driver_drift.py is what makes the copy visible. +# --------------------------------------------------------------------------------------- +def _rvs_lnL_convention(use_lnL=None): + """Resolve the stored-'integrand' convention for a helper call. + + Returns the explicit argument when given, else the run's `rvs_integrand_is_lnL`. Falls + back to False (the historical linear reading) when that global is absent, which is what + happens when these helpers are lifted out of the driver by the unit tests. Read through + globals() rather than by name so a missing global cannot become a NameError swallowed by + a caller's bare `except Exception`. + """ + if use_lnL is not None: + return bool(use_lnL) + return bool(globals().get('rvs_integrand_is_lnL', False)) + + +def ln_weights_from_rvs(rvs, convert=None, use_lnL=False): + """THE importance log-weight of an _rvs record: lnL + ln(prior) - ln(sampling_prior). + + ONE definition, because the alternative has already cost us. A stored 'log_weights' + column does not mean the same thing in every sampler: mcsamplerPortfolio stores the true + importance weight, but mcsamplerGPU stores tempering_exp*lnL + ln p - ln p_s -- the + ADAPTATION weight, with --adapt-weight-exponent baked in. That exponent is not 1 in + production and --no-adapt drives it to 0, removing the likelihood from the column + entirely. A consumer preferring that cache silently reweights its output by L^(e-1). + + So the cache is never read here: the weight is DERIVED from the canonical components -- + log form first, then the linear (mcsamplerEnsemble) form, out-of-support rows -inf. + Raises when neither set is present: an explicit failure beats a plausible wrong number. + + `use_lnL` is REQUIRED to read the linear form correctly, because mcsamplerEnsemble reuses + 'integrand' for BOTH conventions (it stores lnL when given return_lnI). Taking log() of + lnL compresses tens of nats into log(tens), leaving an almost flat weight vector, and the + positivity cut is wrong in that mode too: non-positive means a low-likelihood point, not + a rejected one, so `ig > 0` would discard every sample with lnL <= 0. + + PASS THE STORED CONVENTION, NOT THE CLI OPTION -- `rvs_integrand_is_lnL`, not + `opts.internal_use_lnL`. In this driver the two genuinely differ: --internal-use-lnL is + also accepted for adaptive_cartesian_gpu and portfolio, which set use_lnL without + return_lnI and still store linear L. + """ + conv = convert if convert is not None else (lambda x: x) + if all(k in rvs for k in ('log_integrand', 'log_joint_prior', 'log_joint_s_prior')): + return (numpy.asarray(conv(rvs['log_integrand']), dtype=float) + + numpy.asarray(conv(rvs['log_joint_prior']), dtype=float) + - numpy.asarray(conv(rvs['log_joint_s_prior']), dtype=float)) + if all(k in rvs for k in ('integrand', 'joint_prior', 'joint_s_prior')): + ig = numpy.asarray(conv(rvs['integrand']), dtype=float) + jp = numpy.asarray(conv(rvs['joint_prior']), dtype=float) + js = numpy.asarray(conv(rvs['joint_s_prior']), dtype=float) + out = numpy.full(len(ig), -numpy.inf) + if use_lnL: + # 'integrand' already holds lnL: do not log it again, do not cut on its sign. + keep = numpy.isfinite(ig) & (jp > 0) & (js > 0) + out[keep] = ig[keep] + numpy.log(jp[keep]) - numpy.log(js[keep]) + else: + keep = (ig > 0) & (jp > 0) & (js > 0) + out[keep] = numpy.log(ig[keep]) + numpy.log(jp[keep]) - numpy.log(js[keep]) + return out + raise Exception("cannot build importance weights from sampler._rvs (keys={})".format( + sorted(rvs.keys()))) + + +def _rvs_len(rvs): + for v in rvs.values(): + try: + return len(numpy.atleast_1d(numpy.asarray(v)).ravel()) + except Exception: + continue + return 0 + + +def _rvs_is_export_resample(sampler): + """True when the ROWS of _rvs were drawn in proportion to weight. + + Set by the samplers at the rebind itself, so it means "the draw FIRED", which is NOT the + same predicate as `opts.fairdraw_extrinsic_output`: the draw is skipped when it would not + shrink the record (n_extr >= len(_rvs)), and then the rows are still the retained set + carrying real importance weights. Keying off the CLI flag would flatten those -- the same + class of error in the other direction. + + SURVIVES POOLING by design. This driver does not pool replicas today; the distinction is + kept anyway so that adding --mc-error-replicas here later cannot quietly get it wrong. + """ + return bool(getattr(sampler, '_rvs_is_fairdraw', False)) + + +def _rvs_is_equal_weight(sampler): + """True when EVERY row of _rvs carries the same posterior weight. + + Two properties, deliberately not one flag: + + rows resampled -- each row drawn proportional to w (per-BLOCK property) + equal weight -- the record as a whole is uniform (property of the WHOLE record) + + A single fair draw has both. A POOLED record has the first and not the second, because + pooling weights block k by the replica evidence Z_k/K. Conflating them broke two things + in opposite directions in the main driver (audit Finding 6), which is why the split is + carried over here even though this driver has no pooling yet. + """ + return (bool(getattr(sampler, '_rvs_is_fairdraw', False)) + and not bool(getattr(sampler, '_rvs_is_pooled', False))) + + +def ln_weights_for_posterior(rvs, sampler, convert=None, use_lnL=None): + """The weights to use when treating an _rvs record as a POSTERIOR SAMPLE SET. + + NOT the same question as `ln_weights_from_rvs`, which answers "what is the importance + weight of this record" and is always right about that. The question here is "how should + these rows be weighted to represent the posterior", and the answer depends on whether the + fair draw already did it. + + A fair-drawn record was resampled WITH REPLACEMENT proportional to w, so its rows are + already an equal-weight draw from the posterior. Weighting them by w again applies w^2 + and over-concentrates the result -- measured at a 13% shift in the posterior mean of a + weight-correlated coordinate (verify_skew.py). + + So: uniform (zero log-weight) for a fair-drawn record, the derived importance weight + otherwise. Returns a float array the length of the record. + + `use_lnL` is passed THROUGH UNRESOLVED, exactly as in the main driver: a caller that + omits it gets the linear reading, not the run's convention. That is a trap in both + drivers, and it is deliberately reproduced rather than fixed here -- a helper of the + same name behaving differently in the two forked drivers would be a worse defect than + the one it fixes. Callers must pass `use_lnL=rvs_integrand_is_lnL`, or route through + `_rvs_lnL_convention` first, the way the main driver's call sites do. + """ + if _rvs_is_equal_weight(sampler): + return numpy.zeros(_rvs_len(rvs), dtype=float) + return numpy.asarray(ln_weights_from_rvs(rvs, convert=convert, use_lnL=use_lnL), + dtype=float) +use_gmm_args = (opts.sampler_method == "GMM") or use_gmm_member +if use_gmm_args: # standalone GMM, or a portfolio carrying a GMM member n_step =pinned_params["n"] n_max_blocks = ((1.0*int(opts.n_max))/n_step) # pairing coordinates for adaptive integration: see definition of order below @@ -1283,13 +1500,491 @@ if use_portfolio: if not(isinstance(opts.sampler_portfolio_args[indx], dict)): print(indx,opts.sampler_portfolio_args[indx]) print(" ARGS ", opts.sampler_portfolio_args) - sampler.setup(portfolio_args=opts.sampler_portfolio_args, **pinned_params) # directly pass all parameters set above to low-level portfolios. In particular, GMM setup + # Assemble freeze-policy overrides from the CLI. Only include options the user actually + # set (None = unset) so the sampler keeps its built-in defaults otherwise. The two VARAHA + # flags are mutually exclusive; --portfolio-varaha-can-freeze wins if both are given. + _freeze_policy_kwargs = {} + if opts.portfolio_grace_iters is not None: + _freeze_policy_kwargs['portfolio_grace_iters'] = opts.portfolio_grace_iters + if opts.portfolio_revive_period is not None: + _freeze_policy_kwargs['portfolio_revive_period'] = opts.portfolio_revive_period + if opts.portfolio_freeze_wt is not None: + _freeze_policy_kwargs['portfolio_freeze_wt'] = opts.portfolio_freeze_wt + if opts.portfolio_varaha_can_freeze: + _freeze_policy_kwargs['portfolio_varaha_never_freeze'] = False + elif opts.portfolio_varaha_never_freeze: + _freeze_policy_kwargs['portfolio_varaha_never_freeze'] = True + # adaptive-probe draw allocation (OPT-IN; off by default in the sampler) + if opts.portfolio_adaptive_alloc: + _freeze_policy_kwargs['portfolio_adaptive_alloc'] = True + if opts.portfolio_varaha_min_frac is not None: + _freeze_policy_kwargs['portfolio_varaha_min_frac'] = opts.portfolio_varaha_min_frac + if opts.portfolio_varaha_max_frac is not None: + _freeze_policy_kwargs['portfolio_varaha_max_frac'] = opts.portfolio_varaha_max_frac + if opts.portfolio_weight_clip is not None: + _freeze_policy_kwargs['portfolio_weight_clip'] = opts.portfolio_weight_clip + if opts.portfolio_quality_signal is not None: + _freeze_policy_kwargs['portfolio_quality_signal'] = opts.portfolio_quality_signal + if opts.portfolio_alloc_exponent is not None: + _freeze_policy_kwargs['portfolio_alloc_exponent'] = opts.portfolio_alloc_exponent + if opts.portfolio_probe_period is not None: + _freeze_policy_kwargs['portfolio_probe_period'] = opts.portfolio_probe_period + print(" PORTFOLIO freeze-policy overrides: ", _freeze_policy_kwargs) + sampler.setup(portfolio_args=opts.sampler_portfolio_args, **_freeze_policy_kwargs, **pinned_params) # directly pass all parameters set above to low-level portfolios. In particular, GMM setup # initialize sampler, before we call integrate, so we can seed it if opts.sampler_method == 'adaptive_cartesian_gpu' and opts.skymap_file: sampler.setup() +# --------------------------------------------------------------------------------------- +# L0 auto-rescue. Ported from bin/integrate_likelihood_extrinsic_batchmode (PR #79/#84/#87); +# see test/expensive_before_merging/integrators/RVS_FAIRDRAW_AUDIT.md Findings 1 and 5. +# +# ONE DELIBERATE STRUCTURAL DIVERGENCE FROM THE MAIN DRIVER. There the rescue is inline in +# the single analyze_event. This driver has TWO -- analyze_event_LISA (used with --LISA) and +# analyze_event (the non-LISA fallback) -- each with its own integrate call and export block, +# already ~50% duplicated. Inlining the rescue twice would create a third copy to keep in +# step, which is the failure mode this whole exercise exists to prevent. So the block lives +# in _maybe_l0_rescue below and both call it. The helpers are byte-identical to main's and +# are pinned that way by test_lisa_l0_rescue.py. +# --------------------------------------------------------------------------------------- +def _lnZ_of_rvs(rvs, already_pooled=True, use_lnL=None): + """log of the evidence implied by an _rvs record. + + For a POOLED record the weights already carry their 1/(K n_k) factor, so the estimate is the + plain sum; for a single run it is the mean. Returns None when the weights cannot be rebuilt. + """ + try: + try: + lw = ln_weights_from_rvs(rvs, use_lnL=_rvs_lnL_convention(use_lnL)) + except Exception: + return None + lw = lw[numpy.isfinite(lw)] + if lw.size == 0: + return None + m = numpy.max(lw) + tot = m + numpy.log(numpy.sum(numpy.exp(lw - m))) + return float(tot if already_pooled else tot - numpy.log(lw.size)) + except Exception: + return None + + +def _kish_neff_of_rvs(rvs, use_lnL=None): + """Kish effective sample size of an _rvs record, or None if the weights are not reconstructible.""" + try: + try: + lw = ln_weights_from_rvs(rvs, use_lnL=_rvs_lnL_convention(use_lnL)) + except Exception: + return None + lw = lw[numpy.isfinite(lw)] + if lw.size == 0: + return None + lw = lw - numpy.max(lw) + w = numpy.exp(lw) + return float(numpy.sum(w) ** 2 / numpy.sum(w ** 2)) + except Exception: + return None + + +def _lnZ_of_reserve_or_rvs(sampler, rvs, reserve=None): + """lnZ of a completed pass, from the points it RETAINED where that is available. + + The rescue's reject gate compares the warm pass's lnZ against the cold pass's, and both + were read out of _rvs -- which the fair draw has already replaced with + min(n_extr, 1.5*eff_samp, 1.5*neff) rows resampled WITH REPLACEMENT, proportional to + weight. That is not a smaller unbiased sample of the same estimator, it is a DIFFERENT + and biased one: _lnZ_of_rvs forms logsumexp(w)/n, so drawing n rows proportional to w + returns something near max(w) rather than mean(w), high by roughly + + log(n_retained / eff_samp) + + and the two passes are drawn at wildly different n and eff_samp. The gate was reading a + multi-nat artifact of its own two subsample sizes as evidence that the warm seed had + missed mass. + + Falls back to the old _rvs reading when no reserve was kept, so the comparison degrades to + the previous behaviour rather than to no gate at all. + + Returns (lnZ, source) -- the caller MUST check that both sides came from the same source, + because the two readings are not interchangeable. + """ + _res = reserve if reserve is not None else getattr(sampler, '_warm_seed_reserve', None) + if isinstance(_res, dict) and 'log_joint_prior' in _res and 'log_joint_s_prior' in _res: + try: + # NOT _lnZ_of_rvs: it averages over the rows it is handed, and the reserve is + # neither the draw set nor a uniform sample of it. lnZ_from_reserve restores the + # original proposal-draw normalization from n_finite/n_retained; without it a + # PORTFOLIO reading is high by ~log(n_retained/n_finite), and the error does NOT + # cancel in the gate because the two passes have different finite fractions. + _v = mcsamplerAdaptiveVolume.lnZ_from_reserve(_res) + if _v is not None and numpy.isfinite(_v): + return _v, 'retained' + except Exception: + pass + return _lnZ_of_rvs(rvs, already_pooled=False), 'fairdraw' + + +def _snapshot_pass_state(sampler, res, var, neff, dict_return, rvs=None): + """Everything that must move TOGETHER when a completed pass is put back -> dict. + + THE POINT IS THE WORD "everything". A pass is described by more than its samples, and the + reject path used to restore only some of it: `_rvs`, the estimate and `dict_return` went + back to the cold pass while `_warm_seed_reserve` was left holding the REJECTED warm cloud. + In the main driver that stayed latent until --sampler-sequential-warmstart began seeding + the next intrinsic point from the reserve, at which point a rejected, truncated warm pass + became the seed for the next point -- the exact failure the reject gate exists to prevent, + reintroduced one attribute over. THAT OPTION DOES NOT EXIST IN THIS DRIVER YET, so the + reserve restore is pre-emptive here; it is also what makes porting the capture safe, which + is why it lands first. + + So the snapshot carries the reserve and the fair-draw marker as well, including the + per-member reserves: `_warm_seed_reserve_for` falls through to `portfolio_realizations`, + so restoring only the aggregate would leave that fallback pointing at the warm pass. + """ + return dict( + rvs=(dict(sampler._rvs) if rvs is None else rvs), + res=res, var=var, neff=neff, dict_return=dict_return, + warm_seed_reserve=getattr(sampler, '_warm_seed_reserve', None), + rvs_is_fairdraw=bool(getattr(sampler, '_rvs_is_fairdraw', False)), + rvs_is_pooled=bool(getattr(sampler, '_rvs_is_pooled', False)), + member_reserves=[getattr(_m, '_warm_seed_reserve', None) + for _m in list(getattr(sampler, 'portfolio_realizations', []) or [])], + ) + + +def _restore_pass_state(sampler, state): + """Undo of _snapshot_pass_state -> (res, var, neff, dict_return). + + Both callers (the reject path and the exception handler) go through here, so the set of + attributes that travels with a restored pass cannot drift between them. + """ + sampler._rvs = state['rvs'] + sampler._warm_seed_reserve = state['warm_seed_reserve'] + sampler._rvs_is_fairdraw = state['rvs_is_fairdraw'] + sampler._rvs_is_pooled = state['rvs_is_pooled'] + _members = list(getattr(sampler, 'portfolio_realizations', []) or []) + for _m, _r in zip(_members, state.get('member_reserves', [])): + _m._warm_seed_reserve = _r + return state['res'], state['var'], state['neff'], state['dict_return'] + + +def _warm_seed_reserve_for(sampler): + """The retained-sample reserve a completed pass left behind, or None. + + THE ONE LOOKUP FOR SEED CONSUMERS. In the main driver there are two -- the L0 auto-rescue + (which re-seeds a collapsed pass from its own peak) and the --sampler-sequential-warmstart + capture (which seeds the NEXT intrinsic point). ONLY THE RESCUE EXISTS IN THIS DRIVER; the + shared lookup is kept so the pair cannot drift once the capture is ported. Both otherwise + fall back to sampler._rvs, which by then has been rebound to a fair-draw subset + taken WITH REPLACEMENT -- and the whole point of the reserve is that on the collapsed pass + a warm start exists for, that subset is a handful of rows several of which are the same + point twice. + + A PORTFOLIO keeps the reserve on the aggregate, not on its members, but a bare AV member + can be the one that has it; check the sampler first, then its realizations. + + COLUMN ORDER MUST MATCH or the seed is scrambled: the reserve stores X in the column order + of the sampler that built it, and a seed handed to bootstrap_from_samples is read + positionally against params_ordered. A mismatch is silent and produces a seed in the + wrong coordinates, so decline the reserve rather than use it. + + NOT SHARED WITH `_lnZ_of_reserve_or_rvs` above, deliberately. That one reads only lnL and + the two prior columns, never X, so a column-order mismatch is harmless to it and declining + would throw away a good lnZ reading and silently downgrade the gate to its fallback. + """ + _res = getattr(sampler, '_warm_seed_reserve', None) + if _res is None: + for _m in list(getattr(sampler, 'portfolio_realizations', []) or []): + _res = getattr(_m, '_warm_seed_reserve', None) + if _res is not None: + break + if _res is not None and list(_res.get('params_ordered', [])) != list(sampler.params_ordered): + return None + return _res + + +def _warm_seed_geometry(sampler): + """Which columns a warm seed must span, and the box it must lie in -> (axes, lo, hi). + + The seed is judged on the ADAPTIVE axes, because those are the only ones the live-volume + grid resolves (the rest get a single bin), and that is the set the [AV COLLAPSE] report + counts against. Ask the sampler that will consume the seed rather than assuming all + dimensions: with --force-adapt-all they coincide, without it a rank test over every column + would demand a seed span directions the grid cannot resolve and puff for nothing. + + A PORTFOLIO has no adaptive axes of its own -- they live on its AV-style members -- so fall + through to the first member that can answer. If nobody can, every column it is. + """ + _lo = np.array([sampler.llim[p] for p in sampler.params_ordered], dtype=float) + _hi = np.array([sampler.rlim[p] for p in sampler.params_ordered], dtype=float) + for _s in [sampler] + list(getattr(sampler, 'portfolio_realizations', []) or []): + if hasattr(_s, 'warm_seed_axes'): + try: + return list(_s.warm_seed_axes()), _lo, _hi + except Exception: + pass + return list(range(len(sampler.params_ordered))), _lo, _hi + + +def _clear_warm_state(sampler): + """Clear a warm-start seed AND any grid it installed, reaching PORTFOLIO MEMBERS too. + + `sampler._warm = None` alone is not enough for mcsamplerPortfolio: `_warm` and the + contracted AV grid live on each MEMBER, and portfolio.integrate_log() does not rerun each + member's setup(), so the next point would silently draw from the PREVIOUS point's + contracted live volume. If the new point's support falls outside it, lnZ is biased low + with a healthy-looking n_eff and no error. Portfolio exposes clear_warm_state(); + everything else keeps the old behaviour. + """ + # Deliberately NOT wrapped in try/except. A reset that quietly did not happen leaves the + # next point drawing from the previous point's contracted grid -- the exact silent bias + # this guards against -- so a failure must abort the point rather than degrade to a log + # line nobody reads. + if hasattr(sampler, 'clear_warm_state'): + sampler.clear_warm_state() + else: + sampler._warm = None + sampler._warm_applied = False + + +def _maybe_load_av_state(sampler): + """Warm-start this integration from a saved AV live-volume state (--sampler-load-state).""" + try: + if opts.sampler_load_state and hasattr(sampler, 'load_state'): + print(" warm-start: loading saved sampler state from", opts.sampler_load_state) + sampler.load_state(opts.sampler_load_state) + except Exception as _e_ls: + print(" AV state load skipped (", _e_ls, ")") + + +def _maybe_save_av_state(sampler): + """Persist the adapted live-volume state for reuse by later instances/iterations.""" + if opts.sampler_method == 'AV' and opts.sampler_save_state and hasattr(sampler, 'save_state'): + if not getattr(sampler, '_av_state_reuse_safe', True): + print(" AV: not saving live-volume state from a rejected/failed rescue") + return + try: + sampler.save_state(opts.sampler_save_state) + print(" AV: saved live-volume state to", opts.sampler_save_state) + except Exception as _e_ss: + print(" AV: could not save state (", _e_ss, ")") + + +def _maybe_enable_anisotropic_bins(sampler): + """Opt-in per-axis bin allocation, on the AV sampler AND any AV portfolio members.""" + if getattr(opts, 'sampler_anisotropic_bins', False): + _aniso_targets = [sampler] + list(getattr(sampler, 'portfolio_realizations', [])) + for _t in _aniso_targets: + if hasattr(_t, 'anisotropic_bins'): + _t.anisotropic_bins = True + print(" AV: anisotropic per-axis bin allocation ENABLED") + + +def _reject_if_collapsed(dd, stage): + """Apply --reject-collapsed-live-volume to whatever the CURRENT verdict is. + + In the main driver this is called TWICE -- once on the first run and again on the + replica pool, because replication can turn a healthy first run into a collapsed POOL. + This driver has no replica pooling yet, so only the first call exists here; the second + call site must be added WITH --mc-error-replicas, or the flag is silently bypassed for + exactly the case pooling introduces. + """ + if not opts.reject_collapsed_live_volume: + return + if not (isinstance(dd, dict) and dd.get('live_volume_collapsed', False)): + return + # Route through the ordinary failure path, so the caller skips this binary and writes no + # result row -- the pre-fix outcome, but for a stated reason. + _exc = mcsamplerAdaptiveVolume.LiveVolumeCollapse if mcsampler_AV_ok else RuntimeError + raise _exc( + "extrinsic integration collapsed ({}): live volume degenerated ({}); " + "--reject-collapsed-live-volume is set, so this event is being dropped " + "rather than exported".format(stage, dd.get('collapse_reason', ''))) + + +def _report_and_gate_collapse(dict_return, stage="first run"): + """Announce a collapsed live volume, then apply the rejection gate.""" + _collapsed = bool(dict_return.get('live_volume_collapsed', False)) if isinstance(dict_return, dict) else False + _collapse_reason = (dict_return or {}).get('collapse_reason', '') if isinstance(dict_return, dict) else '' + if _collapsed: + print(" [mc error] *** LIVE VOLUME COLLAPSED *** {}".format(_collapse_reason)) + print(" [mc error] this event's lnZ and exported samples are NOT a fair draw from the posterior.") + _reject_if_collapsed(dict_return, stage) + + +def _maybe_l0_rescue(sampler, res, var, neff, dict_return, + like_to_integrate, unpinned_params, pinned_params, + lnL_offset=0.0): + """Run the L0 auto-rescue if this pass stalled -> (res, var, neff, dict_return). + + Returns its arguments unchanged when the rescue does not apply, so the call site is a + single unconditional assignment. + + MUST BE CALLED BEFORE the `if not(res): raise` guard. A degenerate early termination + (mcsamplerPortfolio/AV returning (None,None,None,None) when the live volume never found + finite in-volume samples) is the STRONGEST rescue trigger, not a reason to skip -- such a + pass still populated _rvs, so the peak seed is available. In the main driver that guard + sits ~200 lines further down and the ordering is implicit; here it is immediately after + integrate, so the ordering is stated and pinned by a test. + + `lnL_offset` is this event's manual_avoid_overflow_logarithm, used only to print absolute + lnZ values. It is a local of the caller in both analyze_event variants, hence a parameter. + """ + # Reset per event before ANY early return. The sampler is reused: a rejected rescue on + # one event must not suppress saving a later healthy event that needs no rescue at all. + sampler._av_state_reuse_safe = True + # APPLICABILITY FIRST, then n_eff. The main driver evaluates + # _neff_val = None if neff is None else float(sampler.identity_convert(neff)) + # BEFORE its guard, which is safe there only by luck: identity_convert comes from + # MCSamplerGeneric, and RIFT.integrators.mcsampler.MCSampler -- the object this driver + # keeps for --sampler-method adaptive_cartesian -- does NOT inherit it. Evaluating it + # unconditionally therefore raises AttributeError on EVERY adaptive_cartesian event, at + # the end of a completed integration and before --output-file is written, losing the + # whole point's compute. The rescue is AV/portfolio-only regardless, so nothing is lost + # by asking whether it applies before touching the sampler's conversion helpers. + # + # DELIBERATE DIVERGENCE from the main driver, which has the same latent defect on the + # line above its own guard and should take the same reordering. + if not (opts.sampler_method in ('AV', 'portfolio') and opts.sampler_warmstart_retry_neff + and hasattr(sampler, 'bootstrap_from_samples')): + return res, var, neff, dict_return + # A DEGENERATE EARLY TERMINATION (neff None) counts as below threshold, not as "skip". + _neff_val = None if neff is None else float(sampler.identity_convert(neff)) + _needs_l0_rescue = (_neff_val is None) or (_neff_val < float(opts.sampler_warmstart_retry_neff or 0)) + if not _needs_l0_rescue: + return res, var, neff, dict_return + + # Cold state to fall back on, captured only once the warm pass is actually about to run. + # `None` means nothing has been disturbed yet, so the handler must not "restore". + _cold_state_l0 = None + try: + # SEED FROM THE POINTS THE PASS RETAINED, not from what survived the fair draw. + # sampler._rvs has by now been REBOUND to a fair-draw subset taken WITH REPLACEMENT -- + # a resample built for EXPORT. On the collapsed pass this rescue exists for the + # effective sample size is ~1, so _rvs can be a single row, or a handful several of + # which are the same point twice. The live set held a thousand. + _res_l0 = _warm_seed_reserve_for(sampler) + if _res_l0 is not None: + _cols = np.asarray(_res_l0['X'], dtype=float) + _lnv = np.asarray(_res_l0['lnL'], dtype=float).ravel() + print(" [L0 auto-rescue] seeding from {} retained sample(s) of {} (fair draw left {} in _rvs)".format( + len(_lnv), _res_l0.get('n_retained', '?'), + len(np.asarray(sampler.identity_convert(sampler._rvs['log_integrand'])).ravel()) + if 'log_integrand' in sampler._rvs else '?')) + else: + _lnkey = 'log_integrand' if 'log_integrand' in sampler._rvs else ('integrand' if 'integrand' in sampler._rvs else None) + _lnv = np.asarray(sampler.identity_convert(sampler._rvs[_lnkey]), dtype=float).ravel() if _lnkey else np.array([]) + _cols = (np.vstack([np.asarray(sampler.identity_convert(sampler._rvs[p]), dtype=float).ravel() + for p in sampler.params_ordered]).T if _lnv.size else np.zeros((0, len(sampler.params_ordered)))) + if _lnv.size >= 1 and np.any(np.isfinite(_lnv)): + # RANK, not count, decides whether this seed can define a live volume. A 2-to-5 + # point seed passes a count test and is still rank-deficient in 6 adaptive + # dimensions, so the warm start contracts onto a degenerate subspace and reports a + # healthy n_eff over a sliver of the support. build_warm_seed applies the rank + # test through the SAME seed_affine_rank the grid builder uses, and puffs to full + # rank when it is short. + _ax_l0, _lo_l0, _hi_l0 = _warm_seed_geometry(sampler) + _seed, _seed_info = mcsamplerAdaptiveVolume.build_warm_seed( + _cols, _lnv, _lo_l0, _hi_l0, _ax_l0, + deltalnL=opts.sampler_sequential_warmstart_deltalnL, + puff_scale=opts.sampler_l0_rescue_puff_scale, + puff_width_frac=opts.sampler_l0_rescue_puff_width_frac, + puff_factor=opts.sampler_l0_rescue_puff_factor) + print(" [L0 auto-rescue] cold n_eff {} < {}; re-running warm from this point's peak ({} pts)".format( + "DEGENERATE (early termination)" if _neff_val is None else "{:.1f}".format(_neff_val), + opts.sampler_warmstart_retry_neff, len(_seed))) + if _seed_info['puffed']: + print(" [L0 auto-rescue] seed of {} point(s) had affine rank {}/{}: PUFFED to rank" + " {}/{} with {} points ({} scale, x{:g}), keeping the original point(s)".format( + _seed_info['n_core'], _seed_info['rank_core'], _seed_info['dim'], + _seed_info['rank_final'], _seed_info['dim'], _seed_info['n_puff'], + _seed_info['puff_scale'], opts.sampler_l0_rescue_puff_factor)) + if _seed_info['rank_final'] < _seed_info['dim']: + print(" [L0 auto-rescue] *** the puffed seed is STILL rank-deficient" + " ({}/{}); the warm pass will be reported as collapsed.".format( + _seed_info['rank_final'], _seed_info['dim'])) + # The warm pass is an estimate over TRUNCATED support: the seeded box provably + # contains the peak the cold pass found, and says nothing about what that pass did + # not reach, so it is biased low by any missed mode. The rescue still runs, because + # it exists to fix the high-SNR n_eff lottery and removing it by default would be a + # certain production regression traded against a possible bias. What the gate below + # changes is only the case where there is POSITIVE EVIDENCE of lost mass. + # + # SNAPSHOT, not an alias: integrate_log repopulates sampler._rvs IN PLACE, so + # `_cold_rvs = sampler._rvs` would be holding the warm samples by the time the + # restore ran -- i.e. the reject path would report the cold lnZ while exporting the + # warm cloud, exactly what it exists to prevent. + _cold_rvs = dict(sampler._rvs) + # Snapshot the RESERVE for the same reason and at the same moment: the warm pass's + # integrate_log clears and rewrites it, so reading it after the fact would compare + # the warm pass against itself. + _cold_reserve_l0 = getattr(sampler, '_warm_seed_reserve', None) + _cold_lnZ, _cold_src = _lnZ_of_reserve_or_rvs(sampler, _cold_rvs, + reserve=_cold_reserve_l0) + # dict_return too: khat, block scatter, ESS and the confidence interval all read it, + # so keeping the warm pass's diagnostics beside a restored cold result would describe + # a run we did not report. And the reserve and the fair-draw marker, one level out. + _cold_state_l0 = _snapshot_pass_state(sampler, res, var, neff, dict_return, + rvs=_cold_rvs) + sampler.bootstrap_from_samples(_seed, cover_frac=0.0) + res, var, neff, dict_return = sampler.integrate(like_to_integrate, *unpinned_params, **pinned_params) + _warm_lnZ, _warm_src = _lnZ_of_reserve_or_rvs(sampler, sampler._rvs) + # BOTH SIDES FROM THE SAME READING, or the difference is not a difference. A + # fair-drawn lnZ sits ~log(n_retained/eff_samp) above a retained-set one, so a mixed + # comparison manufactures a gap of several nats in whichever direction the mismatch + # happens to fall. If the two passes did not produce the same kind of estimate, read + # BOTH from _rvs -- the old behaviour, at least self-consistent. + if _cold_src != _warm_src: + print(" [L0 auto-rescue] lnZ provenance differs (cold={}, warm={});" + " re-reading both from the fair-draw record so the comparison is" + " like-for-like.".format(_cold_src, _warm_src)) + _cold_lnZ = _lnZ_of_rvs(_cold_rvs, already_pooled=False) + _warm_lnZ = _lnZ_of_rvs(sampler._rvs, already_pooled=False) + _cold_src = _warm_src = 'fairdraw' + _evidence_of_loss = ( + (_cold_lnZ is not None) and (_warm_lnZ is not None) + and numpy.isfinite(_cold_lnZ) and numpy.isfinite(_warm_lnZ) + and (_cold_lnZ - _warm_lnZ) > float(opts.sampler_l0_rescue_reject_dlnZ)) + if _evidence_of_loss: + print(" [L0 auto-rescue] *** REJECTING the warm pass *** its lnZ {:.3f} is" + " {:.3f} nats BELOW the full-support cold pass ({:.3f}), which is evidence" + " the seed missed mass the cold pass reached.".format( + _warm_lnZ + lnL_offset, _cold_lnZ - _warm_lnZ, _cold_lnZ + lnL_offset)) + if opts.sampler_l0_rescue_accept_truncated: + print(" [L0 auto-rescue] --sampler-l0-rescue-accept-truncated set:" + " reporting the warm pass anyway (may be biased LOW).") + else: + print(" [L0 auto-rescue] keeping the COLD (full-support) result; its n_eff" + " is lower but it is not missing mass. A portfolio avoids this" + " trade entirely -- its GMM member carries a defensive component.") + # The RESERVE goes back too. Once --sampler-sequential-warmstart is + # ported here, omitting this would seed the next intrinsic point from the + # warm cloud this gate just rejected: _warm_seed_reserve_for would return + # the warm pass's record while _rvs, the estimate and the diagnostics all + # describe the cold one. Nothing reads it in this driver today. + res, var, neff, dict_return = _restore_pass_state(sampler, _cold_state_l0) + sampler._av_state_reuse_safe = False + _clear_warm_state(sampler) + except Exception as _e_l0: + # "skipped" is only true if the warm pass never started. If it raised PARTWAY THROUGH + # sampler.integrate(), the assignment never completed, so res/var/neff/dict_return still + # hold the COLD pass -- while sampler._rvs was repopulated in place and now holds the + # WARM samples. Reporting cold k-hat / ESS / lnZ beside a warm export describes a run + # that was never made, and it did so silently for a whole campaign. + print(" [L0 auto-rescue] *** FAILED *** (", _e_l0, ")") + import traceback as _tb_l0 + _tb_l0.print_exc() + if _cold_state_l0 is not None: + print(" [L0 auto-rescue] the warm pass may already have replaced the stored" + " samples; restoring the COLD pass so the reported diagnostics and the" + " exported samples describe the same integral.") + res, var, neff, dict_return = _restore_pass_state(sampler, _cold_state_l0) + sampler._av_state_reuse_safe = False + _clear_warm_state(sampler) + return res, var, neff, dict_return + + def resample_samples_LISA(my_samples, rholms, cross_terms, right_ascension, declination, P, modes, reference_distance): """This function takes in extrinsic samples and for each sample samples a time shift. This is done by generating a likelihood time series at an extrinsic sample and then weighted sampling in time.""" # How many time samples? Same as the extrinsic samples being passed @@ -1487,7 +2182,7 @@ def analyze_event_LISA(P_list, indx_event, data_dict, psd_dict, fmax, opts, inv_ if 'distance' in sampler.params: sampler.reset_sampling('distance') sampler.reset_sampling('inclination') - elif opts.sampler_method == "GMM": + elif use_gmm_args: # standalone GMM or a portfolio with a GMM member (gmm_dict exists in both) if 'distance' in sampler.params: pair_d_incl = sampler_param_tuple(sampler, ['distance','inclination']) if pair_d_incl in gmm_dict: @@ -1517,11 +2212,31 @@ def analyze_event_LISA(P_list, indx_event, data_dict, psd_dict, fmax, opts, inv_ lnL_oracles = np.zeros(opts.n_chunk) sampler.update_sampling_prior(lnL_oracles, opts.n_chunk, external_rvs=rvs_train,log_scale_weights=True,floor_integrated_probability=opts.adapt_floor_level) + _maybe_load_av_state(sampler) + _maybe_enable_anisotropic_bins(sampler) res, var, neff, dict_return = sampler.integrate(like_to_integrate, *unpinned_params, **pinned_params) + # L0 auto-rescue: on a very sharply-peaked (high-amplitude) point a cold AV can stall at + # n_eff ~ 1 because it never draws near the tiny peak. If so, seed a SECOND pass from this + # same point's own highest-likelihood samples and re-run. Opt-in via + # --sampler-warmstart-retry-neff. MUST run BEFORE the not(res) guard below: a degenerate + # early termination returns (None,None,None,None) and is the strongest rescue trigger, so + # raising on it first would skip exactly the case the rescue exists for. + res, var, neff, dict_return = _maybe_l0_rescue( + sampler, res, var, neff, dict_return, + like_to_integrate, unpinned_params, pinned_params, + lnL_offset=manual_avoid_overflow_logarithm) + if not(res): # no resut raise ValueError(" No integral result returned") + # Collapse gate AFTER the result check, matching the main driver's ordering. + _report_and_gate_collapse(dict_return, "first run") + # Persist only a result we are actually willing to report. In particular, never write + # a collapsed grid that --reject-collapsed-live-volume just rejected, nor a warm grid + # whose rescue result was rejected/failed and replaced by the cold estimate. + _maybe_save_av_state(sampler) + if not(opts.internal_use_lnL): log_res = numpy.log(res) sqrt_var_over_res = numpy.sqrt(var)/res @@ -2191,7 +2906,7 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts, inv_spec_ if 'distance' in sampler.params: sampler.reset_sampling('distance') sampler.reset_sampling('inclination') - elif opts.sampler_method == "GMM": + elif use_gmm_args: # standalone GMM or a portfolio with a GMM member (gmm_dict exists in both) if 'distance' in sampler.params: pair_d_incl = sampler_param_tuple(sampler, ['distance','inclination']) if pair_d_incl in gmm_dict: @@ -2222,11 +2937,29 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts, inv_spec_ lnL_oracles = np.zeros(opts.n_chunk) sampler.update_sampling_prior(lnL_oracles, opts.n_chunk, external_rvs=rvs_train,log_scale_weights=True,floor_integrated_probability=opts.adapt_floor_level) + _maybe_load_av_state(sampler) + _maybe_enable_anisotropic_bins(sampler) res, var, neff, dict_return = sampler.integrate(like_to_integrate, *unpinned_params, **pinned_params) + # L0 auto-rescue: on a very sharply-peaked (high-amplitude) point a cold AV can stall at + # n_eff ~ 1 because it never draws near the tiny peak. If so, seed a SECOND pass from this + # same point's own highest-likelihood samples and re-run. Opt-in via + # --sampler-warmstart-retry-neff. MUST run BEFORE the not(res) guard below: a degenerate + # early termination returns (None,None,None,None) and is the strongest rescue trigger, so + # raising on it first would skip exactly the case the rescue exists for. + res, var, neff, dict_return = _maybe_l0_rescue( + sampler, res, var, neff, dict_return, + like_to_integrate, unpinned_params, pinned_params, + lnL_offset=manual_avoid_overflow_logarithm) + if not(res): # no resut raise ValueError(" No integral result returned") + # Collapse gate AFTER the result check, matching the main driver's ordering. + _report_and_gate_collapse(dict_return, "first run") + # See the LISA variant above: only persist accepted, reusable AV state. + _maybe_save_av_state(sampler) + if not(opts.internal_use_lnL): log_res = numpy.log(res) sqrt_var_over_res = numpy.sqrt(var)/res @@ -2468,7 +3201,7 @@ for indx in numpy.arange(len(P_list)): if opts.sampler_method == "adaptive_cartesian_gpu": for name in sampler.params: sampler.reset_sampling(name) - elif opts.sampler_method == "GMM": + elif use_gmm_args: # standalone GMM or a portfolio with a GMM member # reset the GMM dictionary for component in gmm_dict: gmm_dict[component] = None diff --git a/MonteCarloMarginalizeCode/Code/test/expensive_before_merging/integrators/LISA_DRIVER_DRIFT.md b/MonteCarloMarginalizeCode/Code/test/expensive_before_merging/integrators/LISA_DRIVER_DRIFT.md new file mode 100644 index 000000000..12ad51dcb --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/expensive_before_merging/integrators/LISA_DRIVER_DRIFT.md @@ -0,0 +1,179 @@ +# The LISA ILE driver, against the main one + +The two drivers are a **deliberate fork** (RO, 2026-08-13: *"It is super annoying we have to +have two of them, but the overhead of one ring to rule them all is too high."*). Nothing here +argues for merging them. The purpose is to make the *consequence* of the fork -- drift -- +mechanically visible, so it stays a choice. + + bin/integrate_likelihood_extrinsic_batchmode 4,883 lines moves fast + bin/integrate_likelihood_extrinsic_batchmode_lisa 2,526 lines lags + +Both import the SAME integrators and expose the SAME `ok_lnL_methods` +(`GMM, adaptive_cartesian, adaptive_cartesian_gpu, AV, portfolio` -- verified identical), so +anything landed in `RIFT/integrators/` already reaches LISA. **All measured drift is in the +driver.** + +## How to regenerate this + +Do not trust the numbers below; they are a snapshot. The tooling is the authority. + +``` +python3 audit_lisa_driver_drift.py --summary # counts per category and decision +python3 audit_lisa_driver_drift.py --undecided # what nobody has classified yet +python3 audit_lisa_driver_drift.py --check # the CI gate +python3 make_lisa_drift_ledger.py # regenerate lisa_drift_ledger.json +``` + +`audit_lisa_driver_drift.py` extracts four categories from both drivers by AST and diffs them: +`FUNC` (def names, qualified by enclosing function), `OPTION` (`--foo` literals given to +`add_option`/`add_argument`), `CONST` (module-level `UPPER_CASE`), `ATTR` (sampler provenance +markers -- `_rvs_is_*`, `_warm_seed*`, including the `getattr(obj, 'name', default)` form, +which is how the readers actually access them). + +The judgements live in `make_lisa_drift_ledger.py` as ordered +(pattern -> decision + reason) rules, first match wins, so a whole family is decided once. +An item matching no rule is reported and left out, which fails `--check`. That is the +intended path for newly-drifted code: **a person has to classify it.** + +## What this audit CANNOT see + +Stated plainly, because an adversarial review defeated the gate with four realistic drifts +and the honest answer is that some of them are out of scope by construction rather than by +oversight. + +**It is a NAME-PRESENCE set difference.** It answers "does the LISA driver have a thing +called X". It does not compare behaviour. So all of these produce **zero** gap items: + +* a **changed default** on an option present in both drivers (`--adapt-floor-level` going + 0.1 -> 0.9 is invisible here); +* **changed help text**; +* a **changed body** of a same-named function -- the anti-drift tests in + `test/test_lisa_*.py` cover this for the specific helpers that were ported, and nothing + covers it for anything else; +* a **missing `if` branch or `pinned_params` key**, which is not a FUNC/OPTION/CONST/ATTR at + all. A real example is below. + +**Option names built at runtime evade the extractor.** `add_option(_name_var, ...)`, +options added in a `for` loop, and `"--evade-" + "concat"` are all missed, because the +extractor reads string LITERALS out of the AST. Since `OPTION` is the large majority of the +gap, this is the biggest hole. Neither driver does any of this today. + +**`ATTR` is presence-anywhere.** A marker READ but never WRITTEN counts as present, so a +reader-ported/writer-missing port looks closed. `_rvs_is_pooled` is exactly that today, and +its ledger entry says so. + +The gate is worth having anyway -- it catches the ordinary case, which is a helper or an +option appearing in the main driver and nobody asking the LISA question. It is not a proof +of equivalence, and it should not be described as one. + +## The gate + +`test/test_lisa_driver_drift.py`, wired into the `lisa-check` CI job via +`.travis/test-lisa.sh`. It does **not** assert the gap is empty or that anything was ported. +It asserts that every gap item carries one of `PORT` / `PORTED` / `NA` / `PHYSICS` **with a +reason**, that no item claims `PORTED` while still absent, and that the ledger holds no +entries for items that have left the gap. + +*"Does not apply to LISA" is a fine answer; silence is not.* + +## Snapshot, 2026-08-15 (junior/rift_O4d @ 364a22fd) + +132 items before this pass; 8 ported here, leaving 124. + +| decision | n | meaning | +|---|---|---| +| `PORT` | 70 | belongs in LISA, not there yet -- open work | +| `NA` | 43 | does not apply, with the reason | +| `PHYSICS` | 11 | blocked on a physics decision, with the question | +| `PORTED` | 8 | carried across in this pass | + +### Ported in this pass -- the fair-draw correctness family (PR #87) + +`ln_weights_from_rvs`, `ln_weights_for_posterior`, `_rvs_is_export_resample`, +`_rvs_is_equal_weight`, `_rvs_len`, `_rvs_lnL_convention`, and reads of the `_rvs_is_fairdraw` +/ `_rvs_is_pooled` markers. + +The three consumers whose double-weighting PR #87 actually fixed -- the +`--extrinsic-proposal-output` breadcrumb, the `.dgrid` exporter, the `.dslice` reweight core +-- **do not exist in the LISA driver**, so there was no live `w^2` bug there. What existed was +the hazard: the LISA driver sets `igrand_fairdraw_samples` from `--fairdraw-extrinsic-output`, +so its `_rvs` can be a fair draw, and all seven shared rebind sites already set +`_rvs_is_fairdraw`. **The marker was arriving and nothing read it.** This port is preventive, +and it is the "correct thing to reach for" that the audit's Recommendation 1 asks for. + +Tests: `test/test_lisa_fairdraw_weights.py` (29), revert-checked -- each fix broken in turn, +the named test confirmed failing, the file restored and verified byte-identical. + +Two things deliberately NOT done: + +* `ln_weights_for_posterior` passes `use_lnL` **through unresolved**, exactly as the main + driver does, so a caller that omits it gets the linear reading rather than the run's + convention. That is a latent trap **in both drivers**; reproducing it beats having a + same-named helper behave differently in the two forks. Worth fixing in both, together. +* `_truthy_option` was initially classified with this family and moved out: its only caller in + the main driver is the `--interpolate-time` normalizer, so porting it here would have added + dead code. + +### `NA` -- does not apply to LISA (43) + +| family | n | why | +|---|---|---| +| `--calibration-*` + 4 helpers | 19 | LIGO/Virgo **spline calibration envelopes**. The LISA driver models no instrument calibration: no envelope directory, no cal nodes, response applied analytically by `factored_likelihood_LISA`. | +| `.dslice` / `.dgrid` distance export | 11 | Data products for a downstream LIGO CIP distance workflow the LISA pipeline does not run. No consumer. | +| `--freqresponse*` | 3 | Finite light-travel-time across the arms for **3G ground** detectors, on `lalsimulation` geometry with an arm length in metres. LISA's finite-size response is not an add-on -- it is the TDI response the driver already applies. | +| `--rotation-*` | 3 | Sidereal time-dependence of an **Earth-based** antenna pattern. The constellation's motion is already in the LISA response; this would apply Earth rotation to a heliocentric detector. | +| data/waveform io | 6 | LISA has its own equivalents under different names -- `--data-integration-window-half` for the storage window, `--internal-waveform-*` fd/L-frame passthroughs, h5 frames instead of gwpy, rate from the frame rather than `--srate-internal`. | +| `--e-freq`, `--save-meanPerAno` | 2 | Ground-based eccentric-waveform path (TEOBResumS); LISA's own export is `--save-eccentricity`. | + +### `PHYSICS` -- needs a decision before it can be answered (11) + +These are the ones that need you, not more code reading. + +1. **`--d-prior-redshift`, `dLofz`, `dVdz`** (4 items incl. constants) — *which cosmology and + which redshift range should a LISA distance prior use?* Arguably **more** important for + LISA than for ground-based work, since MBHB sit at z~1-20 where a Euclidean `d^2` prior is + badly wrong -- but the main driver's helpers were built and gridded for the ground-based + range. +2. **`--internal-reparam-dl-incl`, `_reparam_A_of_incl`, `_REPARAM_*`** (5 items) — *does the + quadrupole amplitude `A(iota)=sqrt(((1+cos^2 i)/2)^2+cos^2 i)` remain the right axis to + reparameterize distance against under the LISA TDI response?* It is a pure l=|m|=2 + statement; LISA MBHB are strongly higher-mode and TDI mixes the polarizations differently, + so the degeneracy it straightens may not be the degeneracy LISA has. +3. **`--limit-right-ascension`, `--limit-declination`** — *what should a sky zoom box mean for + LISA?* The driver reuses the key names `right_ascension`/`declination` for its sampled sky + pair, but the values are ecliptic and may be further rotated by + `--internal-sky-network-coordinates`. LISA already has + `--ecliptic-latitude`/`--ecliptic-longitude`/`--lisa-fixed-sky`, which may be the intended + mechanism. (`--limit-psi`/`--limit-inclination` have no such ambiguity and are `PORT`.) +4. **`--sampler-warmstart-samples`** — *what frame are the named columns of a LISA pilot file + in?* Same key-names-different-meaning problem; needs a stated convention, or a pilot + written by the LISA driver itself. + +### `PORT` -- open work, highest value first (70) + +Nothing here is blocked on physics; all of it is sampler-agnostic or pure plumbing. + +| family | n | note | +|---|---|---| +| L0 rescue + warm-start state | 15 | **Highest value.** Triggers on low `n_eff`; LISA MBHB are high-SNR, the regime that stalls. `_snapshot_pass_state`/`_restore_pass_state` must port **as a set** -- Finding 5 was a rejected warm pass restoring `_rvs` but not the reserve. | +| portfolio tuning | 12 | Reachable today via LISA's `--sampler-portfolio-args` eval-dict; porting is pipeline parity. | +| GMM tuning | 7 | Pure pass-through to `mcsamplerEnsemble`. | +| MC-error replicas + pooling | 7 | Includes `_pool_replica_rvs`; port the **per-replica sequence** form, not the boolean (Finding 6). | +| extrinsic proposal handoff | 6 | `--extrinsic-proposal-output` is a Finding-2 site: port it **on top of** `ln_weights_for_posterior`, never with a bare `w`. | +| lnZ / n_eff helpers | 3 | `_lnZ_of_rvs`, `_kish_neff_of_rvs`, `_lnZ_of_reserve_or_rvs` -- needed by the two families above. | +| AV state + binning | 3 | `--sampler-save/load-state`, `--sampler-anisotropic-bins`. | +| misc plumbing | 17 | `--limit-psi`/`--limit-inclination` (port the **post-#58** form, incl. the `cos(iota)` endpoint swap), `--check-good-enough`, `--random-event`, `--fairdraw-extrinsic-output-n-max`, interpolate-time normalizer, etc. | + +**One trap recorded against `--fairdraw-extrinsic-output-n-max`:** the LISA driver currently +hardcodes the cap to `opts.n_eff`, while main's default for the flag is **5**. Adopting main's +default verbatim would silently shrink every LISA export by orders of magnitude. Port the flag +with LISA's present behaviour as its default. + +## Note on CI + +The LISA driver is **not** uncovered -- the `lisa-check` job runs nine test files. But all nine +are import / contract / smoke level: they check the driver loads, exposes its CLI surface and +runs a synthetic demo. None asserts anything about integrator weighting or fair-draw +correctness, which is how 2,357 lines of drift accumulated with CI green. That is the gap the +drift gate closes -- not by testing the physics, but by refusing to let a new item through +without a recorded human decision. diff --git a/MonteCarloMarginalizeCode/Code/test/expensive_before_merging/integrators/audit_lisa_driver_drift.py b/MonteCarloMarginalizeCode/Code/test/expensive_before_merging/integrators/audit_lisa_driver_drift.py new file mode 100644 index 000000000..94ced56c9 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/expensive_before_merging/integrators/audit_lisa_driver_drift.py @@ -0,0 +1,272 @@ +#!/usr/bin/env python3 +""" +Audit: what the main ILE driver has that the LISA ILE driver does not. + +The two drivers are a DELIBERATE FORK (RO, 2026-08-13: "the overhead of one ring to +rule them all is too high"). This script does not argue with that. It makes the +consequence -- drift -- mechanically visible, so the fork stays a choice rather than +an accident. + + bin/integrate_likelihood_extrinsic_batchmode <- main, moves fast + bin/integrate_likelihood_extrinsic_batchmode_lisa <- LISA, lags + +Both import the SAME integrators (``mcsampler``, ``mcsamplerEnsemble``, ``mcsamplerGPU``, +``mcsamplerAdaptiveVolume``, ``mcsamplerPortfolio``), so anything landed in +``RIFT/integrators/`` already reaches LISA. The drift measured here is entirely in the +driver: helpers, CLI options, module constants and sampler provenance markers. + +WHAT IS EXTRACTED +----------------- +``FUNC`` ``def`` names, qualified by enclosing function (``analyze_event._foo``), so a + nested helper is not confused with a top-level one of the same name. +``OPTION`` ``--foo`` literals passed to ``add_option``/``add_argument``. These drivers + use ``optparse``; both call forms are scanned so a future port to argparse + does not silently empty this category. +``CONST`` module-level ``UPPER_CASE`` assignments -- the sentinels (``_SEQ_WS_PENDING``) + and tuning constants that travel with a feature. +``ATTR`` provenance markers set/read on the sampler object (``_rvs_is_fairdraw``, + ``_warm_seed_reserve``, ...). These are the fair-draw correctness family + from PR #87 and are the reason this audit exists. + +THE LEDGER +---------- +Every gap item needs a recorded decision in ``lisa_drift_ledger.json``: + +``PORT`` belongs in LISA and is not there yet -- an open work item. +``PORTED`` carried across; the item should have disappeared from the gap, so a + ``PORTED`` entry still showing up in the gap is itself an error. +``NA`` does not apply to LISA, WITH A REASON. "Does not apply" is a fine answer; + silence is not. +``PHYSICS`` needs a physics decision before it can be answered, with the question + recorded verbatim. + +USAGE +----- + python3 audit_lisa_driver_drift.py # human-readable gap report + python3 audit_lisa_driver_drift.py --summary # counts per category and decision + python3 audit_lisa_driver_drift.py --json # machine-readable + python3 audit_lisa_driver_drift.py --undecided # only items with no ledger entry + python3 audit_lisa_driver_drift.py --check # CI gate: exit 1 on an undecided item + +``--check`` is the CI form, and it is deliberately weak about physics: it does not assert +that the gap is empty, or that any particular item was ported. Closing the gap is not the +goal -- the fork is intentional. It asserts only that no item drifted in unnoticed. A new +helper or option in the main driver fails the build until a person classifies it, which is +the property we want and the one that was missing when 2,357 lines accumulated. + +Keyed by NAME, not by source hash (the fair-draw audit next door keys by hash because it +tracks reads of one attribute, which move). Names are the stable identity here: renaming a +helper in the main driver SHOULD invalidate its verdict, since the thing being tracked is +"does LISA have this", and a rename means nobody has answered that about the new name. + +Needs Python >= 3.8. +""" +import argparse +import ast +import json +import os +import sys + +HERE = os.path.dirname(os.path.abspath(__file__)) +CODE_ROOT = os.path.abspath(os.path.join(HERE, "..", "..", "..")) + +MAIN = "bin/integrate_likelihood_extrinsic_batchmode" +LISA = "bin/integrate_likelihood_extrinsic_batchmode_lisa" + +LEDGER = os.path.join(HERE, "lisa_drift_ledger.json") + +DECISIONS = ("PORT", "PORTED", "NA", "PHYSICS") + +# Sampler attributes worth tracking as provenance markers. Prefix-matched. Kept narrow +# on purpose: every one of these is a boolean or a record describing MUTABLE SHARED STATE, +# which is the shape that produced six defects in PR #87 (see RVS_FAIRDRAW_AUDIT.md). +ATTR_PREFIXES = ("_rvs_is", "_warm_seed", "_retained", "_export_") + + +def _is_str(node): + return isinstance(node, ast.Constant) and isinstance(node.value, str) + + +class _Collector(ast.NodeVisitor): + def __init__(self): + self.funcs = {} # qualified name -> lineno + self.options = {} # "--foo" -> lineno + self.consts = {} # NAME -> lineno + self.attrs = {} # attr name -> lineno + self._stack = [] + + def visit_FunctionDef(self, node): + qual = ".".join(self._stack + [node.name]) + self.funcs.setdefault(qual, node.lineno) + self._stack.append(node.name) + self.generic_visit(node) + self._stack.pop() + + visit_AsyncFunctionDef = visit_FunctionDef + + def visit_Call(self, node): + func = node.func + if isinstance(func, ast.Attribute) and func.attr in ("add_option", "add_argument"): + for arg in node.args: + if _is_str(arg) and arg.value.startswith("--"): + self.options.setdefault(arg.value, node.lineno) + # getattr(sampler, '_rvs_is_fairdraw', False) names an attribute just as much as + # sampler._rvs_is_fairdraw does, and the defensive getattr form is the one the + # provenance READERS use. Missing it would let a real port look like a no-op. + if isinstance(func, ast.Name) and func.id in ("getattr", "setattr", "hasattr"): + for arg in node.args[1:2]: + if _is_str(arg) and any(arg.value.startswith(p) for p in ATTR_PREFIXES): + self.attrs.setdefault(arg.value, node.lineno) + self.generic_visit(node) + + def visit_Assign(self, node): + if not self._stack: + for tgt in node.targets: + name = getattr(tgt, "id", None) + if name and name.upper() == name and any(c.isalpha() for c in name): + self.consts.setdefault(name, node.lineno) + self.generic_visit(node) + + def visit_Attribute(self, node): + if any(node.attr.startswith(p) for p in ATTR_PREFIXES): + self.attrs.setdefault(node.attr, node.lineno) + self.generic_visit(node) + + +def collect(relpath): + path = os.path.join(CODE_ROOT, relpath) + with open(path) as fh: + tree = ast.parse(fh.read(), filename=path) + c = _Collector() + c.visit(tree) + return {"FUNC": c.funcs, "OPTION": c.options, "CONST": c.consts, "ATTR": c.attrs} + + +def compute_gap(): + """Items present in the main driver and absent from the LISA driver. + + Returns (gap, extras) where gap is a list of dicts and extras lists LISA-only + items -- reported but never gated, since LISA is allowed its own surface. + """ + main = collect(MAIN) + lisa = collect(LISA) + gap, extras = [], [] + # A FUNC is satisfied by its BARE name as well as its qualified one. The main driver has + # ONE analyze_event and nests helpers inside it; this driver has TWO (analyze_event_LISA + # and analyze_event), so a helper ported here must be hoisted to module level or else + # duplicated -- and duplicating is the failure mode this audit exists to prevent. Without + # this, every correctly-hoisted port would sit in the gap forever as a false positive, + # which is how a gate gets trained out of people. + _lisa_bare = {n.rsplit(".", 1)[-1] for n in lisa["FUNC"]} + for cat in ("FUNC", "OPTION", "CONST", "ATTR"): + for name in sorted(set(main[cat]) - set(lisa[cat])): + if cat == "FUNC" and name.rsplit(".", 1)[-1] in _lisa_bare: + continue + gap.append({"category": cat, "name": name, + "key": "%s:%s" % (cat, name), "main_line": main[cat][name]}) + for name in sorted(set(lisa[cat]) - set(main[cat])): + extras.append({"category": cat, "name": name, "lisa_line": lisa[cat][name]}) + return gap, extras + + +def load_ledger(): + """Missing ledger => empty => --check reports every item, which is the safe direction.""" + if not os.path.exists(LEDGER): + return {} + with open(LEDGER) as fh: + raw = json.load(fh) + return raw.get("entries", raw) + + +def annotate(gap, ledger): + for item in gap: + entry = ledger.get(item["key"]) + item["decision"] = entry.get("decision") if entry else None + item["reason"] = entry.get("reason") if entry else None + return gap + + +def main(): + ap = argparse.ArgumentParser( + description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--json", action="store_true") + ap.add_argument("--summary", action="store_true") + ap.add_argument("--undecided", action="store_true") + ap.add_argument("--check", action="store_true") + args = ap.parse_args() + + gap, extras = compute_gap() + ledger = load_ledger() + gap = annotate(gap, ledger) + + undecided = [g for g in gap if g["decision"] is None] + # A PORTED item that is still missing from LISA means the ledger is lying about the + # tree -- either the port was reverted or it never landed. Louder than undecided. + stale = [g for g in gap if g["decision"] == "PORTED"] + # A ledger entry naming an item no longer in the gap is spent: either it was ported + # (good) or the main driver dropped it (also fine). Not a failure, but worth showing + # so the ledger does not accumulate fiction. + gap_keys = {g["key"] for g in gap} + spent = sorted(k for k in ledger if k not in gap_keys) + + if args.json: + json.dump({"gap": gap, "lisa_only": extras, "undecided": len(undecided), + "stale_ported": [s["key"] for s in stale], "spent_entries": spent}, + sys.stdout, indent=2, sort_keys=True) + print() + return 0 + + if args.summary: + print("LISA driver drift: %d items in main and absent from lisa" % len(gap)) + for cat in ("FUNC", "OPTION", "CONST", "ATTR"): + rows = [g for g in gap if g["category"] == cat] + if not rows: + continue + counts = {} + for r in rows: + counts[r["decision"] or "UNDECIDED"] = counts.get(r["decision"] or "UNDECIDED", 0) + 1 + detail = " ".join("%s=%d" % (k, counts[k]) for k in sorted(counts)) + print(" %-7s %3d %s" % (cat, len(rows), detail)) + print(" LISA-only surface (never gated): %d" % len(extras)) + if spent: + print(" spent ledger entries (no longer in gap): %d" % len(spent)) + return 0 + + rows = undecided if args.undecided else gap + if args.undecided and not rows: + print("no undecided items: every gap item carries a recorded decision") + for cat in ("FUNC", "OPTION", "CONST", "ATTR"): + sel = [g for g in rows if g["category"] == cat] + if not sel: + continue + print("=== %s (%d)" % (cat, len(sel))) + for g in sel: + print(" %-12s %-52s main:%d" % (g["decision"] or "UNDECIDED", g["name"], g["main_line"])) + if g["reason"]: + print(" %s" % g["reason"]) + print() + + if args.check: + rc = 0 + if undecided: + print("FAIL: %d gap item(s) carry no decision in %s" % ( + len(undecided), os.path.basename(LEDGER)), file=sys.stderr) + for g in undecided: + print(" %s (main:%d)" % (g["key"], g["main_line"]), file=sys.stderr) + print("\nClassify each as PORT / PORTED / NA / PHYSICS with a reason.", + file=sys.stderr) + rc = 1 + if stale: + print("FAIL: %d item(s) marked PORTED are still absent from the LISA driver:" + % len(stale), file=sys.stderr) + for g in stale: + print(" %s" % g["key"], file=sys.stderr) + rc = 1 + if rc == 0: + print("OK: all %d gap items carry a recorded decision" % len(gap)) + return rc + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/MonteCarloMarginalizeCode/Code/test/expensive_before_merging/integrators/lisa_drift_ledger.json b/MonteCarloMarginalizeCode/Code/test/expensive_before_merging/integrators/lisa_drift_ledger.json new file mode 100644 index 000000000..945ff8d11 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/expensive_before_merging/integrators/lisa_drift_ledger.json @@ -0,0 +1,373 @@ +{ + "_comment": "GENERATED by make_lisa_drift_ledger.py -- edit the RULES there, not this file.", + "entries": { + "CONST:_REPARAM_A_MAX": { + "decision": "PHYSICS", + "reason": "Tuning constants for --internal-reparam-dl-incl." + }, + "CONST:_REPARAM_A_MIN": { + "decision": "PHYSICS", + "reason": "Tuning constants for --internal-reparam-dl-incl." + }, + "CONST:_REPARAM_LNF": { + "decision": "PHYSICS", + "reason": "Tuning constants for --internal-reparam-dl-incl." + }, + "CONST:_SEQ_WS_PENDING": { + "decision": "PORT", + "reason": "Sentinel for the deferred sequential warm-start capture; ports with --sampler-sequential-warmstart." + }, + "CONST:_TI_LEGACY_BOOLEAN": { + "decision": "PORT", + "reason": "Legacy-boolean vocabulary for --interpolate-time. Main (PR #97) now accepts STENCIL NAMES there -- nearest/cubic/sinc -- normalizing into opts._noloop_time_interp, with this tuple for back-compat and an explicit typo guard so a misspelling is not absorbed as falsey. LISA still passes the raw --interpolate-time value straight to the likelihood, so porting means normalizing it AND teaching the LISA time path the stencil name; it travels with _normalize_interpolate_time_argv and _truthy_option." + }, + "FUNC:_cal_setup_prior_with_nodes": { + "decision": "NA", + "reason": "Calibration-envelope internals; see the --calibration-* reason." + }, + "FUNC:_draw_more_calibration_draws": { + "decision": "NA", + "reason": "Calibration-envelope internals; see the --calibration-* reason." + }, + "FUNC:_normalize_interpolate_time_argv": { + "decision": "PORT", + "reason": "Normalizes --interpolate-time argv forms. LISA exposes --interpolate-time, so the same normalization applies." + }, + "FUNC:_pool_replica_rvs": { + "decision": "PORT", + "reason": "Pools replica records by evidence. Ports with --mc-error-replicas. NOTE its per-replica already_resampled sequence (Finding 6): a single global boolean is wrong near the n_extr boundary, so port the sequence form, not the boolean." + }, + "FUNC:_pool_replica_rvs._block_resampled": { + "decision": "PORT", + "reason": "Pools replica records by evidence. Ports with --mc-error-replicas. NOTE its per-replica already_resampled sequence (Finding 6): a single global boolean is wrong near the n_extr boundary, so port the sequence form, not the boolean." + }, + "FUNC:_reparam_A_of_incl": { + "decision": "PHYSICS", + "reason": "Implementation of --internal-reparam-dl-incl." + }, + "FUNC:_truthy_option": { + "decision": "PORT", + "reason": "Tolerant truthiness for optparse values that may arrive as strings from the pipe. Belongs with _normalize_interpolate_time_argv, its ONLY caller in the main driver (opts._noloop_time_interp), not with the fair-draw family -- porting it alongside those helpers would have added dead code to the LISA driver." + }, + "FUNC:analyze_event._cal_error_probe": { + "decision": "NA", + "reason": "Calibration Monte-Carlo error probe; see the --calibration-* reason." + }, + "FUNC:analyze_event._cal_error_probe._draw_dist": { + "decision": "NA", + "reason": "Calibration Monte-Carlo error probe; see the --calibration-* reason." + }, + "FUNC:analyze_event._extract_mc_diag": { + "decision": "PORT", + "reason": "Diagnostics for the replica triggers." + }, + "FUNC:dLofz": { + "decision": "PHYSICS", + "reason": "Cosmology helpers behind --d-prior-redshift. Same question: the interpolation range has to be re-chosen for MBHB redshifts." + }, + "FUNC:dVdz": { + "decision": "PHYSICS", + "reason": "Cosmology helpers behind --d-prior-redshift. Same question: the interpolation range has to be re-chosen for MBHB redshifts." + }, + "OPTION:--calibration-burn-in-neff": { + "decision": "NA", + "reason": "LIGO/Virgo spline calibration-envelope marginalization. The LISA driver models no instrument calibration: it takes no envelope directory, has no cal nodes, and its response is applied analytically by factored_likelihood_LISA. LISA calibration, if it is ever modelled, will not have this data product or this spline parameterization, so porting the LIGO machinery would be actively misleading." + }, + "OPTION:--calibration-burn-in-nmax": { + "decision": "NA", + "reason": "LIGO/Virgo spline calibration-envelope marginalization. The LISA driver models no instrument calibration: it takes no envelope directory, has no cal nodes, and its response is applied analytically by factored_likelihood_LISA. LISA calibration, if it is ever modelled, will not have this data product or this spline parameterization, so porting the LIGO machinery would be actively misleading." + }, + "OPTION:--calibration-conjugate-phase": { + "decision": "NA", + "reason": "LIGO/Virgo spline calibration-envelope marginalization. The LISA driver models no instrument calibration: it takes no envelope directory, has no cal nodes, and its response is applied analytically by factored_likelihood_LISA. LISA calibration, if it is ever modelled, will not have this data product or this spline parameterization, so porting the LIGO machinery would be actively misleading." + }, + "OPTION:--calibration-dump-responsibilities": { + "decision": "NA", + "reason": "LIGO/Virgo spline calibration-envelope marginalization. The LISA driver models no instrument calibration: it takes no envelope directory, has no cal nodes, and its response is applied analytically by factored_likelihood_LISA. LISA calibration, if it is ever modelled, will not have this data product or this spline parameterization, so porting the LIGO machinery would be actively misleading." + }, + "OPTION:--calibration-envelope-directory": { + "decision": "NA", + "reason": "LIGO/Virgo spline calibration-envelope marginalization. The LISA driver models no instrument calibration: it takes no envelope directory, has no cal nodes, and its response is applied analytically by factored_likelihood_LISA. LISA calibration, if it is ever modelled, will not have this data product or this spline parameterization, so porting the LIGO machinery would be actively misleading." + }, + "OPTION:--calibration-export-posterior": { + "decision": "NA", + "reason": "LIGO/Virgo spline calibration-envelope marginalization. The LISA driver models no instrument calibration: it takes no envelope directory, has no cal nodes, and its response is applied analytically by factored_likelihood_LISA. LISA calibration, if it is ever modelled, will not have this data product or this spline parameterization, so porting the LIGO machinery would be actively misleading." + }, + "OPTION:--calibration-fused-kernel": { + "decision": "NA", + "reason": "LIGO/Virgo spline calibration-envelope marginalization. The LISA driver models no instrument calibration: it takes no envelope directory, has no cal nodes, and its response is applied analytically by factored_likelihood_LISA. LISA calibration, if it is ever modelled, will not have this data product or this spline parameterization, so porting the LIGO machinery would be actively misleading." + }, + "OPTION:--calibration-global-norm": { + "decision": "NA", + "reason": "LIGO/Virgo spline calibration-envelope marginalization. The LISA driver models no instrument calibration: it takes no envelope directory, has no cal nodes, and its response is applied analytically by factored_likelihood_LISA. LISA calibration, if it is ever modelled, will not have this data product or this spline parameterization, so porting the LIGO machinery would be actively misleading." + }, + "OPTION:--calibration-mc-error-extrinsic": { + "decision": "NA", + "reason": "LIGO/Virgo spline calibration-envelope marginalization. The LISA driver models no instrument calibration: it takes no envelope directory, has no cal nodes, and its response is applied analytically by factored_likelihood_LISA. LISA calibration, if it is ever modelled, will not have this data product or this spline parameterization, so porting the LIGO machinery would be actively misleading." + }, + "OPTION:--calibration-n-realizations": { + "decision": "NA", + "reason": "LIGO/Virgo spline calibration-envelope marginalization. The LISA driver models no instrument calibration: it takes no envelope directory, has no cal nodes, and its response is applied analytically by factored_likelihood_LISA. LISA calibration, if it is ever modelled, will not have this data product or this spline parameterization, so porting the LIGO machinery would be actively misleading." + }, + "OPTION:--calibration-n-realizations-max": { + "decision": "NA", + "reason": "LIGO/Virgo spline calibration-envelope marginalization. The LISA driver models no instrument calibration: it takes no envelope directory, has no cal nodes, and its response is applied analytically by factored_likelihood_LISA. LISA calibration, if it is ever modelled, will not have this data product or this spline parameterization, so porting the LIGO machinery would be actively misleading." + }, + "OPTION:--calibration-neff-cal-target": { + "decision": "NA", + "reason": "LIGO/Virgo spline calibration-envelope marginalization. The LISA driver models no instrument calibration: it takes no envelope directory, has no cal nodes, and its response is applied analytically by factored_likelihood_LISA. LISA calibration, if it is ever modelled, will not have this data product or this spline parameterization, so porting the LIGO machinery would be actively misleading." + }, + "OPTION:--calibration-pilot-extrinsic": { + "decision": "NA", + "reason": "LIGO/Virgo spline calibration-envelope marginalization. The LISA driver models no instrument calibration: it takes no envelope directory, has no cal nodes, and its response is applied analytically by factored_likelihood_LISA. LISA calibration, if it is ever modelled, will not have this data product or this spline parameterization, so porting the LIGO machinery would be actively misleading." + }, + "OPTION:--calibration-proposal-breadcrumb": { + "decision": "NA", + "reason": "LIGO/Virgo spline calibration-envelope marginalization. The LISA driver models no instrument calibration: it takes no envelope directory, has no cal nodes, and its response is applied analytically by factored_likelihood_LISA. LISA calibration, if it is ever modelled, will not have this data product or this spline parameterization, so porting the LIGO machinery would be actively misleading." + }, + "OPTION:--calibration-spline-count": { + "decision": "NA", + "reason": "LIGO/Virgo spline calibration-envelope marginalization. The LISA driver models no instrument calibration: it takes no envelope directory, has no cal nodes, and its response is applied analytically by factored_likelihood_LISA. LISA calibration, if it is ever modelled, will not have this data product or this spline parameterization, so porting the LIGO machinery would be actively misleading." + }, + "OPTION:--check-good-enough": { + "decision": "PORT", + "reason": "Early-exit when the pipeline has written an 'ile_good_enough' sentinel. Pipeline plumbing, detector-agnostic." + }, + "OPTION:--d-prior-redshift": { + "decision": "PHYSICS", + "reason": "QUESTION: which cosmology and which redshift range should a LISA distance prior use? This is arguably MORE important for LISA than for ground-based work -- MBHB sit at z~1-20 where a Euclidean d^2 prior is badly wrong -- but the main driver's helper was built and gridded for the ground-based range. Needs a stated cosmology and a z ceiling before porting." + }, + "OPTION:--distance-slice-all-fresh": { + "decision": "NA", + "reason": "The .dslice export and its placement/tuning knobs. This is a data product for a downstream LIGO CIP distance workflow that the LISA pipeline does not run; there is no consumer. If a LISA distance workflow is ever built, note that the .dslice reweight core was the third Finding-2 site and must not be revived in its pre-#87 form." + }, + "OPTION:--distance-slice-chunk": { + "decision": "NA", + "reason": "The .dslice export and its placement/tuning knobs. This is a data product for a downstream LIGO CIP distance workflow that the LISA pipeline does not run; there is no consumer. If a LISA distance workflow is ever built, note that the .dslice reweight core was the third Finding-2 site and must not be revived in its pre-#87 form." + }, + "OPTION:--distance-slice-randomize": { + "decision": "NA", + "reason": "The .dslice export and its placement/tuning knobs. This is a data product for a downstream LIGO CIP distance workflow that the LISA pipeline does not run; there is no consumer. If a LISA distance workflow is ever built, note that the .dslice reweight core was the third Finding-2 site and must not be revived in its pre-#87 form." + }, + "OPTION:--distance-slice-skip-threshold": { + "decision": "NA", + "reason": "The .dslice export and its placement/tuning knobs. This is a data product for a downstream LIGO CIP distance workflow that the LISA pipeline does not run; there is no consumer. If a LISA distance workflow is ever built, note that the .dslice reweight core was the third Finding-2 site and must not be revived in its pre-#87 form." + }, + "OPTION:--distance-slice-wing-delta-lnL": { + "decision": "NA", + "reason": "The .dslice export and its placement/tuning knobs. This is a data product for a downstream LIGO CIP distance workflow that the LISA pipeline does not run; there is no consumer. If a LISA distance workflow is ever built, note that the .dslice reweight core was the third Finding-2 site and must not be revived in its pre-#87 form." + }, + "OPTION:--distance-slice-wing-neff": { + "decision": "NA", + "reason": "The .dslice export and its placement/tuning knobs. This is a data product for a downstream LIGO CIP distance workflow that the LISA pipeline does not run; there is no consumer. If a LISA distance workflow is ever built, note that the .dslice reweight core was the third Finding-2 site and must not be revived in its pre-#87 form." + }, + "OPTION:--distance-slice-wing-nmax": { + "decision": "NA", + "reason": "The .dslice export and its placement/tuning knobs. This is a data product for a downstream LIGO CIP distance workflow that the LISA pipeline does not run; there is no consumer. If a LISA distance workflow is ever built, note that the .dslice reweight core was the third Finding-2 site and must not be revived in its pre-#87 form." + }, + "OPTION:--e-freq": { + "decision": "NA", + "reason": "TEOBResumS eccentric-frequency convention. Tied to a ground-based eccentric waveform path the LISA driver does not offer (it takes --modes / h5 frames)." + }, + "OPTION:--export-distance-slices": { + "decision": "NA", + "reason": "The .dslice export and its placement/tuning knobs. This is a data product for a downstream LIGO CIP distance workflow that the LISA pipeline does not run; there is no consumer. If a LISA distance workflow is ever built, note that the .dslice reweight core was the third Finding-2 site and must not be revived in its pre-#87 form." + }, + "OPTION:--export-marginal-distance-grid": { + "decision": "NA", + "reason": "The .dgrid export. Same absent consumer as .dslice, and the second Finding-2 double-weighting site." + }, + "OPTION:--extrinsic-proposal-adapt": { + "decision": "PORT", + "reason": "Consumes the breadcrumb above. Ports with it." + }, + "OPTION:--extrinsic-proposal-breadcrumb": { + "decision": "PORT", + "reason": "Consumes the breadcrumb above. Ports with it." + }, + "OPTION:--extrinsic-proposal-field": { + "decision": "PORT", + "reason": "AV proposal-field handoff, built by util_BuildProposalField.py from a previous ILE iteration. Sampler-agnostic; blocked only on the LISA pipeline growing that stage, so it is a work item rather than an exclusion." + }, + "OPTION:--extrinsic-proposal-field-cover-frac": { + "decision": "PORT", + "reason": "AV proposal-field handoff, built by util_BuildProposalField.py from a previous ILE iteration. Sampler-agnostic; blocked only on the LISA pipeline growing that stage, so it is a work item rather than an exclusion." + }, + "OPTION:--extrinsic-proposal-field-inflate": { + "decision": "PORT", + "reason": "AV proposal-field handoff, built by util_BuildProposalField.py from a previous ILE iteration. Sampler-agnostic; blocked only on the LISA pipeline growing that stage, so it is a work item rather than an exclusion." + }, + "OPTION:--extrinsic-proposal-output": { + "decision": "PORT", + "reason": "Fits the run's extrinsic posterior to a GMM and writes it as a breadcrumb. This is one of the three Finding-2 double-weighting sites, so it MUST be ported on top of ln_weights_for_posterior (done here) and never with a bare w." + }, + "OPTION:--fairdraw-extrinsic-output-n-max": { + "decision": "PORT", + "reason": "Caps rows per fair-draw export. LISA currently hardcodes this to opts.n_eff at the igrand_fairdraw_samples_max call site. WARNING for the port: main's default is 5, so adopting main's default verbatim would silently shrink every LISA export by orders of magnitude. Port the flag with LISA's present behaviour as its default." + }, + "OPTION:--freqresponse": { + "decision": "NA", + "reason": "Finite light-travel-time transfer across the arms for 3G ground detectors (CE/ET), built on lalsimulation detector geometry and an arm-length override in metres. LISA's finite-size response is not an add-on: it is the whole point of the TDI response the LISA driver already applies." + }, + "OPTION:--freqresponse-arm-length": { + "decision": "NA", + "reason": "Finite light-travel-time transfer across the arms for 3G ground detectors (CE/ET), built on lalsimulation detector geometry and an arm-length override in metres. LISA's finite-size response is not an add-on: it is the whole point of the TDI response the LISA driver already applies." + }, + "OPTION:--freqresponse-qmax": { + "decision": "NA", + "reason": "Finite light-travel-time transfer across the arms for 3G ground detectors (CE/ET), built on lalsimulation detector geometry and an arm-length override in metres. LISA's finite-size response is not an add-on: it is the whole point of the TDI response the LISA driver already applies." + }, + "OPTION:--internal-data-storage-window-half": { + "decision": "NA", + "reason": "Half-width of the main driver's internal precompute storage window. The LISA driver has its own equivalent under a different name, --data-integration-window-half, which it passes straight into PrecomputeAlignedSpinLISA. Same role, already present." + }, + "OPTION:--internal-gmm-adaptive-components": { + "decision": "PORT", + "reason": "mcsamplerEnsemble (GMM) tuning. LISA wires that sampler and exposes the same 'GMM' method string, so these knobs are reachable physics-wise but simply not plumbed. Pure pass-through. Caveat for whoever ports --internal-gmm-sky-components: the default grouping is (sky)(distance,inclination)(psi,phi), and 'sky' for LISA is the ecliptic pair -- the grouping still makes sense, the docstring does not." + }, + "OPTION:--internal-gmm-correlate-all": { + "decision": "PORT", + "reason": "mcsamplerEnsemble (GMM) tuning. LISA wires that sampler and exposes the same 'GMM' method string, so these knobs are reachable physics-wise but simply not plumbed. Pure pass-through. Caveat for whoever ports --internal-gmm-sky-components: the default grouping is (sky)(distance,inclination)(psi,phi), and 'sky' for LISA is the ecliptic pair -- the grouping still makes sense, the docstring does not." + }, + "OPTION:--internal-gmm-defensive-frac": { + "decision": "PORT", + "reason": "mcsamplerEnsemble (GMM) tuning. LISA wires that sampler and exposes the same 'GMM' method string, so these knobs are reachable physics-wise but simply not plumbed. Pure pass-through. Caveat for whoever ports --internal-gmm-sky-components: the default grouping is (sky)(distance,inclination)(psi,phi), and 'sky' for LISA is the ecliptic pair -- the grouping still makes sense, the docstring does not." + }, + "OPTION:--internal-gmm-inflate": { + "decision": "PORT", + "reason": "mcsamplerEnsemble (GMM) tuning. LISA wires that sampler and exposes the same 'GMM' method string, so these knobs are reachable physics-wise but simply not plumbed. Pure pass-through. Caveat for whoever ports --internal-gmm-sky-components: the default grouping is (sky)(distance,inclination)(psi,phi), and 'sky' for LISA is the ecliptic pair -- the grouping still makes sense, the docstring does not." + }, + "OPTION:--internal-gmm-max-components": { + "decision": "PORT", + "reason": "mcsamplerEnsemble (GMM) tuning. LISA wires that sampler and exposes the same 'GMM' method string, so these knobs are reachable physics-wise but simply not plumbed. Pure pass-through. Caveat for whoever ports --internal-gmm-sky-components: the default grouping is (sky)(distance,inclination)(psi,phi), and 'sky' for LISA is the ecliptic pair -- the grouping still makes sense, the docstring does not." + }, + "OPTION:--internal-gmm-phase-components": { + "decision": "PORT", + "reason": "mcsamplerEnsemble (GMM) tuning. LISA wires that sampler and exposes the same 'GMM' method string, so these knobs are reachable physics-wise but simply not plumbed. Pure pass-through. Caveat for whoever ports --internal-gmm-sky-components: the default grouping is (sky)(distance,inclination)(psi,phi), and 'sky' for LISA is the ecliptic pair -- the grouping still makes sense, the docstring does not." + }, + "OPTION:--internal-gmm-sky-components": { + "decision": "PORT", + "reason": "mcsamplerEnsemble (GMM) tuning. LISA wires that sampler and exposes the same 'GMM' method string, so these knobs are reachable physics-wise but simply not plumbed. Pure pass-through. Caveat for whoever ports --internal-gmm-sky-components: the default grouping is (sky)(distance,inclination)(psi,phi), and 'sky' for LISA is the ecliptic pair -- the grouping still makes sense, the docstring does not." + }, + "OPTION:--internal-precompute-ignore-threshold": { + "decision": "PORT", + "reason": "Drops negligible modes during precompute. LISA is mode-heavy (--modes, --restricted-mode-list-file) and pays more per mode than a ground-based run, so if anything this matters more there. No LIGO-specific assumption." + }, + "OPTION:--internal-reparam-dl-incl": { + "decision": "PHYSICS", + "reason": "QUESTION: does the quadrupole amplitude A(iota)=sqrt(((1+cos^2 i)/2)^2+cos^2 i) remain the right axis to reparameterize distance against under the LISA TDI response? The reparameterization is a pure l=|m|=2 statement; LISA MBHB are strongly higher-mode and the TDI channels mix the two polarizations differently, so the degeneracy it straightens may not be the degeneracy LISA has." + }, + "OPTION:--internal-use-gwpy": { + "decision": "NA", + "reason": "gwpy low-level frame io. The LISA driver reads its data from h5 frames (--h5-frame/--h5-frame-FD), not from GWF via gwpy." + }, + "OPTION:--internal-waveform-extra-kwargs": { + "decision": "NA", + "reason": "lalsimulation taper / extra-kwargs passthrough for the ground-based waveform path. The LISA driver has its own passthroughs for the generator it uses (--internal-waveform-extra-lalsuite-args, --internal-waveform-fd-L-frame, --internal-waveform-fd-no-condition)." + }, + "OPTION:--internal-waveform-taper": { + "decision": "NA", + "reason": "lalsimulation taper / extra-kwargs passthrough for the ground-based waveform path. The LISA driver has its own passthroughs for the generator it uses (--internal-waveform-extra-lalsuite-args, --internal-waveform-fd-L-frame, --internal-waveform-fd-no-condition)." + }, + "OPTION:--limit-declination": { + "decision": "PHYSICS", + "reason": "QUESTION: what should a sky zoom box mean for LISA? The LISA driver reuses the KEY NAMES right_ascension/declination for its sampled sky pair, but the values are ecliptic (lambda,beta) and may be further rotated by --internal-sky-network-coordinates. A box is therefore well-defined only once it is stated which frame the user is quoting -- and LISA already has --ecliptic-latitude/--ecliptic-longitude/--lisa-fixed-sky, which may already be the intended mechanism." + }, + "OPTION:--limit-inclination": { + "decision": "PORT", + "reason": "Zoom-box limits on psi and inclination. These parameters mean the same thing in both drivers and LISA exposes --inclination-cosine-sampler, which is exactly the case junior PR #58 found silently ignored -- so port the POST-#58 form, including the cos(iota) endpoint swap." + }, + "OPTION:--limit-psi": { + "decision": "PORT", + "reason": "Zoom-box limits on psi and inclination. These parameters mean the same thing in both drivers and LISA exposes --inclination-cosine-sampler, which is exactly the case junior PR #58 found silently ignored -- so port the POST-#58 form, including the cos(iota) endpoint swap." + }, + "OPTION:--limit-right-ascension": { + "decision": "PHYSICS", + "reason": "QUESTION: what should a sky zoom box mean for LISA? The LISA driver reuses the KEY NAMES right_ascension/declination for its sampled sky pair, but the values are ecliptic (lambda,beta) and may be further rotated by --internal-sky-network-coordinates. A box is therefore well-defined only once it is stated which frame the user is quoting -- and LISA already has --ecliptic-latitude/--ecliptic-longitude/--lisa-fixed-sky, which may already be the intended mechanism." + }, + "OPTION:--mc-error-ess-trigger": { + "decision": "PORT", + "reason": "Replica-based lnL error stabilization. Triggers on weight-tail diagnostics of the run's own weights; nothing detector-specific. Valuable for LISA for the same reason as for high-SNR ground events: the reported sigma is the thing downstream CIP trusts." + }, + "OPTION:--mc-error-khat-trigger": { + "decision": "PORT", + "reason": "Replica-based lnL error stabilization. Triggers on weight-tail diagnostics of the run's own weights; nothing detector-specific. Valuable for LISA for the same reason as for high-SNR ground events: the reported sigma is the thing downstream CIP trusts." + }, + "OPTION:--mc-error-replicas": { + "decision": "PORT", + "reason": "Replica-based lnL error stabilization. Triggers on weight-tail diagnostics of the run's own weights; nothing detector-specific. Valuable for LISA for the same reason as for high-SNR ground events: the reported sigma is the thing downstream CIP trusts." + }, + "OPTION:--mc-error-sigma-trigger": { + "decision": "PORT", + "reason": "Replica-based lnL error stabilization. Triggers on weight-tail diagnostics of the run's own weights; nothing detector-specific. Valuable for LISA for the same reason as for high-SNR ground events: the reported sigma is the thing downstream CIP trusts." + }, + "OPTION:--n-distance-slice-core": { + "decision": "NA", + "reason": "The .dslice export and its placement/tuning knobs. This is a data product for a downstream LIGO CIP distance workflow that the LISA pipeline does not run; there is no consumer. If a LISA distance workflow is ever built, note that the .dslice reweight core was the third Finding-2 site and must not be revived in its pre-#87 form." + }, + "OPTION:--n-distance-slice-wing": { + "decision": "NA", + "reason": "The .dslice export and its placement/tuning knobs. This is a data product for a downstream LIGO CIP distance workflow that the LISA pipeline does not run; there is no consumer. If a LISA distance workflow is ever built, note that the .dslice reweight core was the third Finding-2 site and must not be revived in its pre-#87 form." + }, + "OPTION:--nf-flow-load": { + "decision": "PORT", + "reason": "Normalizing-flow persistence is detector-agnostic, but the LISA portfolio factory currently constructs only AV, GMM, and adaptive_cartesian_gpu members. Port the NF member construction and route load/save to that member before exposing these flags; hooks on the portfolio aggregate are a silent no-op because it has no flow API." + }, + "OPTION:--nf-flow-save": { + "decision": "PORT", + "reason": "Normalizing-flow persistence is detector-agnostic, but the LISA portfolio factory currently constructs only AV, GMM, and adaptive_cartesian_gpu members. Port the NF member construction and route load/save to that member before exposing these flags; hooks on the portfolio aggregate are a silent no-op because it has no flow API." + }, + "OPTION:--random-event": { + "decision": "PORT", + "reason": "Pick a random event from the input file. Detector-agnostic; flagged dangerous in its own help text for oversampling reasons that apply equally to LISA." + }, + "OPTION:--rotation-n-harmonics": { + "decision": "NA", + "reason": "Sidereal time-dependence of an EARTH-BASED antenna pattern F(t). The LISA constellation's motion is already carried by the LISA response itself (factored_likelihood_LISA + the h5/TDI frames), so this correction is both unnecessary and wrong there -- it would apply Earth rotation to a heliocentric detector." + }, + "OPTION:--rotation-p-max": { + "decision": "NA", + "reason": "Sidereal time-dependence of an EARTH-BASED antenna pattern F(t). The LISA constellation's motion is already carried by the LISA response itself (factored_likelihood_LISA + the h5/TDI frames), so this correction is both unnecessary and wrong there -- it would apply Earth rotation to a heliocentric detector." + }, + "OPTION:--rotation-slow": { + "decision": "NA", + "reason": "Sidereal time-dependence of an EARTH-BASED antenna pattern F(t). The LISA constellation's motion is already carried by the LISA response itself (factored_likelihood_LISA + the h5/TDI frames), so this correction is both unnecessary and wrong there -- it would apply Earth rotation to a heliocentric detector." + }, + "OPTION:--sampler-sequential-warmstart": { + "decision": "PORT", + "reason": "Warm-start each intrinsic point from the previous one's cloud. Applies whenever --n-events-to-analyze>1, which LISA supports. Its snapshot/restore prerequisites (Finding 5) already landed with the L0 rescue, so this is now capture + the event-loop wiring only." + }, + "OPTION:--sampler-sequential-warmstart-cover-frac": { + "decision": "PORT", + "reason": "Coverage floor for the above; meaningless without it, so they travel together." + }, + "OPTION:--sampler-warmstart-cover-frac": { + "decision": "PORT", + "reason": "Coverage floor and inflation for a handed-off seed. Pure geometry on the sampled unit cube." + }, + "OPTION:--sampler-warmstart-inflate": { + "decision": "PORT", + "reason": "Coverage floor and inflation for a handed-off seed. Pure geometry on the sampled unit cube." + }, + "OPTION:--sampler-warmstart-samples": { + "decision": "PHYSICS", + "reason": "QUESTION: what frame are the named columns of a LISA pilot file in? The reader expects right_ascension/declination/inclination/psi/phi_orb/distance, and the LISA driver does use those KEY NAMES internally -- but they carry ecliptic (and, with --internal-sky-network-coordinates, rotated) values, so a file is only meaningful if the writer and reader agree on the convention. Needs a stated convention before it can be ported, or a pilot written by the LISA driver itself." + }, + "OPTION:--save-meanPerAno": { + "decision": "NA", + "reason": "Exports the eccentric mean anomaly. Tied to the ground-based eccentric waveform path (see --e-freq); the LISA driver's own eccentricity export is --save-eccentricity." + }, + "OPTION:--save-samples-process-params": { + "decision": "PORT", + "reason": "Retain the process_params table in the XML output. Pure output plumbing." + }, + "OPTION:--srate-internal": { + "decision": "NA", + "reason": "Separate internal sampling rate for the ground-based precompute. LISA's precompute takes its rate from the h5 frame and P.deltaT; there is no second internal rate to set." + }, + "OPTION:--srate-resample-time-marginalization": { + "decision": "PORT", + "reason": "Interpolate the lnL time series onto a finer grid before time resampling. LISA already has --resample-time-marginalization and its own time-resampling block, so this is the matching resolution knob and applies directly." + } + } +} diff --git a/MonteCarloMarginalizeCode/Code/test/expensive_before_merging/integrators/make_lisa_drift_ledger.py b/MonteCarloMarginalizeCode/Code/test/expensive_before_merging/integrators/make_lisa_drift_ledger.py new file mode 100644 index 000000000..9ccc03994 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/expensive_before_merging/integrators/make_lisa_drift_ledger.py @@ -0,0 +1,393 @@ +#!/usr/bin/env python3 +""" +Regenerate ``lisa_drift_ledger.json`` -- the recorded decision for every item the main +ILE driver has and the LISA ILE driver does not. + + python3 make_lisa_drift_ledger.py # rewrite the ledger + python3 make_lisa_drift_ledger.py --dry-run # show what would change + python3 audit_lisa_driver_drift.py --check # CI gate over the result + +The gap itself is computed by ``audit_lisa_driver_drift.py``; this file holds only the +JUDGEMENTS, as ordered (pattern -> decision + reason) rules so a whole family is decided +once. First match wins, so put specific items above their family. + +DECISIONS + PORT belongs in LISA, not there yet. An open work item. + PORTED carried across. The audit re-checks these: a PORTED item still missing from + the LISA driver fails the build. + NA does not apply to LISA, with the reason. + PHYSICS cannot be answered without a physics decision, with the question. + +An item matching NO rule is reported and left out of the ledger, so ``--check`` fails on +it. That is the intended path for newly-drifted code: it must be classified by a person. + +WHY THESE DECISIONS LOOK THE WAY THEY DO +The two drivers import the SAME integrators and expose the SAME ``ok_lnL_methods`` +(``GMM, adaptive_cartesian, adaptive_cartesian_gpu, AV, portfolio``, verified identical +2026-08-15). So anything that is pure sampler plumbing applies to LISA by construction and +is PORT; the NA items are the ones tied to a ground-based detector, to LIGO/Virgo +calibration envelopes, or to a downstream pipeline stage LISA does not run. +""" +import argparse +import json +import os +import re +import sys + +import audit_lisa_driver_drift as audit + +HERE = os.path.dirname(os.path.abspath(__file__)) +OUT = os.path.join(HERE, "lisa_drift_ledger.json") + +# --------------------------------------------------------------------------------------- +# Ordered rules. (regex over the audit key "CATEGORY:name", decision, reason) +# First match wins. +# --------------------------------------------------------------------------------------- +RULES = [ + + # ---------------------------------------------------------------- the fair-draw family + # PORTED in this pass. These are the PR #87 correctness helpers. They are pure + # functions of the _rvs record plus the sampler's own provenance markers, and the + # markers are already set by the shared integrators at all seven rebind sites, so they + # already arrive on LISA's sampler objects at runtime -- only the driver-side readers + # were missing. + (r"^FUNC:ln_weights_from_rvs$", "PORTED", + "Importance weight of an _rvs record. Pure function of the record; no extrinsic " + "coordinate assumptions. LISA sets igrand_fairdraw_samples, so its records can be " + "fair draws and need the same answer."), + (r"^FUNC:ln_weights_for_posterior$", "PORTED", + "How rows should be weighted to REPRESENT THE POSTERIOR, as distinct from their " + "importance weight. Returns zeros on an equal-weight record. This is the helper " + "that makes the w^2 double-weighting defect unrepresentable."), + (r"^FUNC:_rvs_is_export_resample$", "PORTED", + "Predicate: rows were drawn proportional to w (survives pooling). Reads the shared " + "marker the integrators already set."), + (r"^FUNC:_rvs_is_equal_weight$", "PORTED", + "Predicate: record is globally equal-weight (fairdraw and not pooled). Finding 6 " + "split this from _rvs_is_export_resample; porting one without the other rebuilds " + "the flag-answering-two-questions bug."), + (r"^FUNC:_rvs_len$", "PORTED", + "Row count of an _rvs record, tolerant of the tuple-keyed sky column. Support " + "helper for the above."), + (r"^ATTR:_rvs_is_fairdraw$", "PORTED", + "Set by all seven shared rebind sites in RIFT/integrators/, so it already reaches " + "LISA at runtime; the LISA driver simply never read it."), + (r"^ATTR:_rvs_is_pooled$", "PORTED", + "READER ONLY, deliberately. The marker is read by _rvs_is_equal_weight and carried " + "by the pass snapshot/restore; nothing in this driver ever SETS it, because there " + "is no replica pooling here yet. Main's reset-on-entry (Finding 7: the marker " + "outliving a FAILED event) is therefore NOT ported and MUST come with " + "--mc-error-replicas -- without it the first pooled record would leave the marker " + "set on the next event. Note this is also the ATTR category's blind spot: a name " + "read anywhere counts as present, so reader-ported/writer-missing looks closed."), + + # ---------------------------------------------------------------------- lnZ / n_eff + (r"^FUNC:_lnZ_of_rvs$", "PORTED", + "Evidence of an _rvs record with the already_pooled/fairdraw correction. Landed " + "with the L0 rescue gate, which is its first consumer here."), + (r"^FUNC:_kish_neff_of_rvs$", "PORTED", + "Kish n_eff of a record. Landed with _lnZ_of_rvs; its own consumer (replica " + "pooling) arrives in the MC-error pass."), + (r"^FUNC:_lnZ_of_reserve_or_rvs$", "PORTED", + "Reads a pass's lnZ from the points it RETAINED where available, so the reject " + "gate is not comparing two differently-sized fair-draw artifacts."), + (r"^FUNC:(_snapshot_pass_state|_restore_pass_state)$", "PORTED", + "Snapshot/restore of everything that must travel with a put-back pass -- the " + "reserve and the fair-draw marker included (Finding 5). Ported as a SET with the " + "rescue; either one alone rebuilds the defect."), + (r"^FUNC:(_warm_seed_reserve_for|_warm_seed_geometry|_clear_warm_state)$", "PORTED", + "Shared reserve lookup (with the column-order guard), adaptive-axis geometry for " + "the rank test, and the warm-state clear that reaches portfolio MEMBERS."), + (r"^ATTR:_warm_seed_reserve$", "PORTED", + "The retained-sample reserve the rescue seeds from and the snapshot carries."), + (r"^OPTION:--sampler-warmstart-retry-neff$", "PORTED", + "The L0 rescue trigger. High value for LISA: MBHB are high-SNR, which is the " + "regime that stalls at n_eff~1."), + (r"^OPTION:--sampler-l0-rescue-", "PORTED", + "L0 rescue tuning, defaults and help text kept identical to the main driver " + "(including reject-dlnZ 3.0, the measured value -- see " + "L0_REJECT_DLNZ_MEASUREMENT.md). Pinned by test_lisa_l0_rescue.py."), + (r"^OPTION:--sampler-sequential-warmstart-deltalnL$", "PORTED", + "The lnL window build_warm_seed keeps. Consumed by the L0 rescue, so it landed " + "with that pass rather than with the sequential warm start it is named for."), + + # --------------------------------------------------------------- L0 rescue / warm start + (r"^OPTION:--reject-collapsed-live-volume$", "PORTED", + "AV live-volume collapse rejection. AV is wired in the LISA driver identically. " + "NOTE the main driver calls its gate TWICE -- first run and replica pool -- and only " + "the first call exists here, because there is no pooling yet; the second MUST be " + "added with --mc-error-replicas or the flag is bypassed for the case pooling creates."), + (r"^FUNC:analyze_event\._reject_if_collapsed$", "PORTED", + "Hoisted to module level rather than nested, because this driver has TWO " + "analyze_event variants. The audit matches FUNC items on the bare name for exactly " + "this reason."), + (r"^OPTION:--sampler-sequential-warmstart$", "PORT", + "Warm-start each intrinsic point from the previous one's cloud. Applies whenever " + "--n-events-to-analyze>1, which LISA supports. Its snapshot/restore prerequisites " + "(Finding 5) already landed with the L0 rescue, so this is now capture + the " + "event-loop wiring only."), + (r"^OPTION:--sampler-sequential-warmstart-cover-frac$", "PORT", + "Coverage floor for the above; meaningless without it, so they travel together."), + (r"^OPTION:--sampler-anisotropic-bins$", "PORTED", + "AV per-axis bin counts during contraction. AV is wired in LISA, and the argument " + "for it is if anything stronger there: the LISA extrinsic axes are no more " + "isotropic than the ground-based ones, and a sky pair that localizes tightly " + "while distance stays broad is the exact case this exists for."), + (r"^OPTION:--sampler-(save|load)-state$", "PORTED", + "AV live-volume state serialization. AV is wired in LISA; the state is the " + "sampler's own internal grid, so it carries no LIGO-specific convention."), + (r"^OPTION:--sampler-warmstart-(cover-frac|inflate)$", "PORT", + "Coverage floor and inflation for a handed-off seed. Pure geometry on the " + "sampled unit cube."), + (r"^OPTION:--sampler-warmstart-samples$", "PHYSICS", + "QUESTION: what frame are the named columns of a LISA pilot file in? The reader " + "expects right_ascension/declination/inclination/psi/phi_orb/distance, and the " + "LISA driver does use those KEY NAMES internally -- but they carry ecliptic " + "(and, with --internal-sky-network-coordinates, rotated) values, so a file is only " + "meaningful if the writer and reader agree on the convention. Needs a stated " + "convention before it can be ported, or a pilot written by the LISA driver itself."), + + # --------------------------------------------------------------------- MC error replicas + (r"^OPTION:--mc-error-(replicas|sigma-trigger|ess-trigger|khat-trigger)$", "PORT", + "Replica-based lnL error stabilization. Triggers on weight-tail diagnostics of the " + "run's own weights; nothing detector-specific. Valuable for LISA for the same " + "reason as for high-SNR ground events: the reported sigma is the thing downstream " + "CIP trusts."), + (r"^FUNC:_pool_replica_rvs(\._block_resampled)?$", "PORT", + "Pools replica records by evidence. Ports with --mc-error-replicas. NOTE its " + "per-replica already_resampled sequence (Finding 6): a single global boolean is " + "wrong near the n_extr boundary, so port the sequence form, not the boolean."), + (r"^FUNC:analyze_event\._extract_mc_diag$", "PORT", "Diagnostics for the replica triggers."), + + # ------------------------------------------------------------------------ GMM plumbing + (r"^OPTION:--internal-gmm-", "PORT", + "mcsamplerEnsemble (GMM) tuning. LISA wires that sampler and exposes the same " + "'GMM' method string, so these knobs are reachable physics-wise but simply not " + "plumbed. Pure pass-through. Caveat for whoever ports --internal-gmm-sky-components: " + "the default grouping is (sky)(distance,inclination)(psi,phi), and 'sky' for LISA " + "is the ecliptic pair -- the grouping still makes sense, the docstring does not."), + + # ------------------------------------------------------------------ portfolio plumbing + (r"^OPTION:--portfolio-", "PORTED", + "mcsamplerPortfolio freeze/allocation policy. Definitions copied verbatim and the " + "_freeze_policy_kwargs assembly is textually identical to the main driver's, so " + "unset options (None) stay out of the dict and the sampler keeps its own defaults. " + "--portfolio-varaha-can-freeze wins over --portfolio-varaha-never-freeze, as there."), + + # ------------------------------------------------------------------- NF flow plumbing + (r"^OPTION:--nf-flow-(load|save)$", "PORT", + "Normalizing-flow persistence is detector-agnostic, but the LISA portfolio factory " + "currently constructs only AV, GMM, and adaptive_cartesian_gpu members. Port the NF " + "member construction and route load/save to that member before exposing these flags; " + "hooks on the portfolio aggregate are a silent no-op because it has no flow API."), + + # --------------------------------------------------------- extrinsic proposal handoff + (r"^OPTION:--extrinsic-proposal-output$", "PORT", + "Fits the run's extrinsic posterior to a GMM and writes it as a breadcrumb. This " + "is one of the three Finding-2 double-weighting sites, so it MUST be ported on top " + "of ln_weights_for_posterior (done here) and never with a bare w."), + (r"^OPTION:--extrinsic-proposal-(breadcrumb|adapt)$", "PORT", + "Consumes the breadcrumb above. Ports with it."), + (r"^OPTION:--extrinsic-proposal-field(-cover-frac|-inflate)?$", "PORT", + "AV proposal-field handoff, built by util_BuildProposalField.py from a previous " + "ILE iteration. Sampler-agnostic; blocked only on the LISA pipeline growing that " + "stage, so it is a work item rather than an exclusion."), + + # ------------------------------------------------------------------------ fair-draw size + (r"^OPTION:--fairdraw-extrinsic-output-n-max$", "PORT", + "Caps rows per fair-draw export. LISA currently hardcodes this to opts.n_eff at " + "the igrand_fairdraw_samples_max call site. WARNING for the port: main's default " + "is 5, so adopting main's default verbatim would silently shrink every LISA " + "export by orders of magnitude. Port the flag with LISA's present behaviour as " + "its default."), + + # ------------------------------------------------------- LIGO/Virgo calibration envelopes + (r"^OPTION:--calibration-", "NA", + "LIGO/Virgo spline calibration-envelope marginalization. The LISA driver models no " + "instrument calibration: it takes no envelope directory, has no cal nodes, and its " + "response is applied analytically by factored_likelihood_LISA. LISA calibration, if " + "it is ever modelled, will not have this data product or this spline parameterization, " + "so porting the LIGO machinery would be actively misleading."), + (r"^FUNC:(_cal_setup_prior_with_nodes|_draw_more_calibration_draws)$", "NA", + "Calibration-envelope internals; see the --calibration-* reason."), + (r"^FUNC:analyze_event\._cal_error_probe(\._draw_dist)?$", "NA", + "Calibration Monte-Carlo error probe; see the --calibration-* reason."), + + # ------------------------------------------------------- ground-based detector geometry + (r"^OPTION:--rotation-(slow|n-harmonics|p-max)$", "NA", + "Sidereal time-dependence of an EARTH-BASED antenna pattern F(t). The LISA " + "constellation's motion is already carried by the LISA response itself " + "(factored_likelihood_LISA + the h5/TDI frames), so this correction is both " + "unnecessary and wrong there -- it would apply Earth rotation to a heliocentric " + "detector."), + (r"^OPTION:--freqresponse(-arm-length|-qmax)?$", "NA", + "Finite light-travel-time transfer across the arms for 3G ground detectors " + "(CE/ET), built on lalsimulation detector geometry and an arm-length override in " + "metres. LISA's finite-size response is not an add-on: it is the whole point of " + "the TDI response the LISA driver already applies."), + (r"^OPTION:--e-freq$", "NA", + "TEOBResumS eccentric-frequency convention. Tied to a ground-based eccentric " + "waveform path the LISA driver does not offer (it takes --modes / h5 frames)."), + + # ---------------------------------------------------------- distance slice / grid export + (r"^OPTION:--(export-distance-slices|distance-slice-|n-distance-slice-)", "NA", + "The .dslice export and its placement/tuning knobs. This is a data product for a " + "downstream LIGO CIP distance workflow that the LISA pipeline does not run; there " + "is no consumer. If a LISA distance workflow is ever built, note that the .dslice " + "reweight core was the third Finding-2 site and must not be revived in its " + "pre-#87 form."), + (r"^OPTION:--export-marginal-distance-grid$", "NA", + "The .dgrid export. Same absent consumer as .dslice, and the second Finding-2 " + "double-weighting site."), + + # ----------------------------------------------------------------- cosmology / d prior + (r"^OPTION:--d-prior-redshift$", "PHYSICS", + "QUESTION: which cosmology and which redshift range should a LISA distance prior " + "use? This is arguably MORE important for LISA than for ground-based work -- MBHB " + "sit at z~1-20 where a Euclidean d^2 prior is badly wrong -- but the main driver's " + "helper was built and gridded for the ground-based range. Needs a stated " + "cosmology and a z ceiling before porting."), + (r"^FUNC:(dLofz|dVdz)$", "PHYSICS", + "Cosmology helpers behind --d-prior-redshift. Same question: the interpolation " + "range has to be re-chosen for MBHB redshifts."), + + # -------------------------------------------------------------- distance/incl reparam + (r"^OPTION:--internal-reparam-dl-incl$", "PHYSICS", + "QUESTION: does the quadrupole amplitude A(iota)=sqrt(((1+cos^2 i)/2)^2+cos^2 i) " + "remain the right axis to reparameterize distance against under the LISA TDI " + "response? The reparameterization is a pure l=|m|=2 statement; LISA MBHB are " + "strongly higher-mode and the TDI channels mix the two polarizations differently, " + "so the degeneracy it straightens may not be the degeneracy LISA has."), + (r"^FUNC:_reparam_A_of_incl$", "PHYSICS", "Implementation of --internal-reparam-dl-incl."), + (r"^CONST:_REPARAM_", "PHYSICS", "Tuning constants for --internal-reparam-dl-incl."), + + # ---------------------------------------------------------------------- extrinsic boxes + (r"^OPTION:--limit-(psi|inclination)$", "PORT", + "Zoom-box limits on psi and inclination. These parameters mean the same thing in " + "both drivers and LISA exposes --inclination-cosine-sampler, which is exactly the " + "case junior PR #58 found silently ignored -- so port the POST-#58 form, including " + "the cos(iota) endpoint swap."), + (r"^OPTION:--limit-(right-ascension|declination)$", "PHYSICS", + "QUESTION: what should a sky zoom box mean for LISA? The LISA driver reuses the " + "KEY NAMES right_ascension/declination for its sampled sky pair, but the values " + "are ecliptic (lambda,beta) and may be further rotated by " + "--internal-sky-network-coordinates. A box is therefore well-defined only once it " + "is stated which frame the user is quoting -- and LISA already has " + "--ecliptic-latitude/--ecliptic-longitude/--lisa-fixed-sky, which may already be " + "the intended mechanism."), + + # --------------------------------------------------------------------- data / waveform io + (r"^OPTION:--internal-data-storage-window-half$", "NA", + "Half-width of the main driver's internal precompute storage window. The LISA " + "driver has its own equivalent under a different name, --data-integration-window-half, " + "which it passes straight into PrecomputeAlignedSpinLISA. Same role, already present."), + (r"^OPTION:--internal-use-gwpy$", "NA", + "gwpy low-level frame io. The LISA driver reads its data from h5 frames " + "(--h5-frame/--h5-frame-FD), not from GWF via gwpy."), + (r"^OPTION:--internal-waveform-(taper|extra-kwargs)$", "NA", + "lalsimulation taper / extra-kwargs passthrough for the ground-based waveform " + "path. The LISA driver has its own passthroughs for the generator it uses " + "(--internal-waveform-extra-lalsuite-args, --internal-waveform-fd-L-frame, " + "--internal-waveform-fd-no-condition)."), + (r"^OPTION:--srate-internal$", "NA", + "Separate internal sampling rate for the ground-based precompute. LISA's " + "precompute takes its rate from the h5 frame and P.deltaT; there is no second " + "internal rate to set."), + (r"^OPTION:--srate-resample-time-marginalization$", "PORT", + "Interpolate the lnL time series onto a finer grid before time resampling. LISA " + "already has --resample-time-marginalization and its own time-resampling block, " + "so this is the matching resolution knob and applies directly."), + (r"^CONST:_TI_LEGACY_BOOLEAN$", "PORT", + "Legacy-boolean vocabulary for --interpolate-time. Main (PR #97) now accepts STENCIL " + "NAMES there -- nearest/cubic/sinc -- normalizing into opts._noloop_time_interp, with " + "this tuple for back-compat and an explicit typo guard so a misspelling is not " + "absorbed as falsey. LISA still passes the raw --interpolate-time value straight to " + "the likelihood, so porting means normalizing it AND teaching the LISA time path the " + "stencil name; it travels with _normalize_interpolate_time_argv and _truthy_option."), + (r"^FUNC:_normalize_interpolate_time_argv$", "PORT", + "Normalizes --interpolate-time argv forms. LISA exposes --interpolate-time, so " + "the same normalization applies."), + (r"^OPTION:--internal-precompute-ignore-threshold$", "PORT", + "Drops negligible modes during precompute. LISA is mode-heavy (--modes, " + "--restricted-mode-list-file) and pays more per mode than a ground-based run, so " + "if anything this matters more there. No LIGO-specific assumption."), + + # ------------------------------------------------------------------------------- misc + (r"^OPTION:--check-good-enough$", "PORT", + "Early-exit when the pipeline has written an 'ile_good_enough' sentinel. Pipeline " + "plumbing, detector-agnostic."), + (r"^OPTION:--random-event$", "PORT", + "Pick a random event from the input file. Detector-agnostic; flagged dangerous in " + "its own help text for oversampling reasons that apply equally to LISA."), + (r"^OPTION:--save-samples-process-params$", "PORT", + "Retain the process_params table in the XML output. Pure output plumbing."), + (r"^OPTION:--save-meanPerAno$", "NA", + "Exports the eccentric mean anomaly. Tied to the ground-based eccentric waveform " + "path (see --e-freq); the LISA driver's own eccentricity export is " + "--save-eccentricity."), + (r"^OPTION:--calibration-spline-count$", "NA", "See the --calibration-* reason."), + (r"^CONST:_SEQ_WS_PENDING$", "PORT", + "Sentinel for the deferred sequential warm-start capture; ports with " + "--sampler-sequential-warmstart."), + (r"^FUNC:_truthy_option$", "PORT", + "Tolerant truthiness for optparse values that may arrive as strings from the pipe. " + "Belongs with _normalize_interpolate_time_argv, its ONLY caller in the main driver " + "(opts._noloop_time_interp), not with the fair-draw family -- porting it alongside " + "those helpers would have added dead code to the LISA driver."), + (r"^FUNC:_rvs_lnL_convention$", "PORTED", + "Resolves the stored-integrand convention from the run's rvs_integrand_is_lnL. " + "Ported alongside the weight helpers because it is how a caller is SUPPOSED to " + "obtain use_lnL: ln_weights_for_posterior passes the argument through unresolved " + "in both drivers, so omitting it silently yields the linear reading."), +] + + +def classify(key): + for pat, decision, reason in RULES: + if re.search(pat, key): + return decision, reason + return None, None + + +def main(): + ap = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--dry-run", action="store_true") + args = ap.parse_args() + + gap, _extras = audit.compute_gap() + entries, unmatched = {}, [] + for item in gap: + decision, reason = classify(item["key"]) + if decision is None: + unmatched.append(item) + continue + entries[item["key"]] = {"decision": decision, "reason": reason} + + counts = {} + for e in entries.values(): + counts[e["decision"]] = counts.get(e["decision"], 0) + 1 + print("classified %d/%d gap items: %s" % ( + len(entries), len(gap), " ".join("%s=%d" % kv for kv in sorted(counts.items())))) + + if unmatched: + print("\n%d item(s) match NO rule -- add one, or they fail --check:" % len(unmatched)) + for item in unmatched: + print(" %-58s main:%d" % (item["key"], item["main_line"])) + + if args.dry_run: + return 1 if unmatched else 0 + + payload = { + "_comment": "GENERATED by make_lisa_drift_ledger.py -- edit the RULES there, not this file.", + "entries": entries, + } + with open(OUT, "w") as fh: + json.dump(payload, fh, indent=2, sort_keys=True) + fh.write("\n") + print("wrote %s" % os.path.relpath(OUT, HERE)) + return 1 if unmatched else 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/MonteCarloMarginalizeCode/Code/test/expensive_before_merging/integrators/make_rvs_fairdraw_ledger.py b/MonteCarloMarginalizeCode/Code/test/expensive_before_merging/integrators/make_rvs_fairdraw_ledger.py index 5dec706af..3ab925189 100644 --- a/MonteCarloMarginalizeCode/Code/test/expensive_before_merging/integrators/make_rvs_fairdraw_ledger.py +++ b/MonteCarloMarginalizeCode/Code/test/expensive_before_merging/integrators/make_rvs_fairdraw_ledger.py @@ -73,7 +73,15 @@ def verdict(h): "so, rather than reporting a plausible wrong number.") # --- the L0 rescue and the sequential warm start --------------------------------- - if f == "bin/integrate_likelihood_extrinsic_batchmode": + # BOTH ILE drivers. The LISA driver now carries a ported copy of the L0 rescue, and the + # `_rvs` reads inside it are byte-identical to the ones here -- these rules match on + # source TEXT, so the same text earns the same verdict. (It lives in a module-level + # _maybe_l0_rescue there rather than inlined in analyze_event, because that driver has TWO + # analyze_event variants; the enclosing function name is not part of the match.) Rules in + # this block naming things the LISA driver does not have -- _rep_rvs, extrinsic_handoff, + # the sequential-warm-start seeds -- simply never match for it. + if f in ("bin/integrate_likelihood_extrinsic_batchmode", + "bin/integrate_likelihood_extrinsic_batchmode_lisa"): if "_lnZ_of_reserve_or_rvs" in s: return ("FIXED", "PR #79, re-landed as #86. sampler._rvs is passed as the FALLBACK " diff --git a/MonteCarloMarginalizeCode/Code/test/expensive_before_merging/integrators/rvs_fairdraw_verdicts.json b/MonteCarloMarginalizeCode/Code/test/expensive_before_merging/integrators/rvs_fairdraw_verdicts.json index bdecbd3d8..a5253302f 100644 --- a/MonteCarloMarginalizeCode/Code/test/expensive_before_merging/integrators/rvs_fairdraw_verdicts.json +++ b/MonteCarloMarginalizeCode/Code/test/expensive_before_merging/integrators/rvs_fairdraw_verdicts.json @@ -419,6 +419,51 @@ "verdict": "BENIGN", "why": "MAP-seed read: argmax over the record, then the SAME row's coordinates. The fair draw picks rows proportional to weight, so the retained argmax is a genuinely-drawn high-likelihood point; it seeds a local search and a printed diagnostic, and no downstream number is a statistic of it. Degraded (it may not be the global MAP), not wrong." }, + "bin/integrate_likelihood_extrinsic_batchmode_lisa:_maybe_l0_rescue:0ab512f38a": { + "source": "_warm_lnZ = _lnZ_of_rvs(sampler._rvs, already_pooled=False)", + "verdict": "BROKEN", + "why": "PR #79's CROSS-SOURCE FALLBACK: fires only when the cold and warm passes produced different reading sources (one had a reserve, the other did not), and then re-reads BOTH sides from the fair-draw record. That is self-consistent, which is what #79 claims for it, but it is not unbiased: the two passes sit at different n_eff, so the log(n/n_eff) artifact does not cancel and this branch is back in the regime measured at +3.48 nats / 100% rejection at the 0.5 default. A known, documented, BOUNDED residual -- not a defect anyone introduced. Closing it needs a retained-set reading on both sides, i.e. a reserve for the samplers that keep none. Follow-up, not a regression." + }, + "bin/integrate_likelihood_extrinsic_batchmode_lisa:_maybe_l0_rescue:1331dd70ea": { + "source": "len(np.asarray(sampler.identity_convert(sampler._rvs['log_integrand'])).ravel())", + "verdict": "BENIGN", + "why": "Reports how many rows the fair draw left, for the log line that contrasts it with the retained count. Reading the resample's size is the POINT here." + }, + "bin/integrate_likelihood_extrinsic_batchmode_lisa:_maybe_l0_rescue:32c64bcd73": { + "source": "_lnv = np.asarray(sampler.identity_convert(sampler._rvs[_lnkey]), dtype=float).ravel() if _lnkey else np.array([])", + "verdict": "BENIGN", + "why": "The DELIBERATE no-reserve fallback for a sampler that keeps none: reads _rvs so the feature degrades to its previous behaviour rather than to no seed at all. The rank hazard is still handled -- build_warm_seed puffs the result to full rank -- so what remains is only 'fewer points', which cannot be improved without a reserve." + }, + "bin/integrate_likelihood_extrinsic_batchmode_lisa:_maybe_l0_rescue:36e58e97fb": { + "source": "if 'log_integrand' in sampler._rvs else '?'))", + "verdict": "PER_ROW", + "why": "Key-presence / column reference, not a population statistic." + }, + "bin/integrate_likelihood_extrinsic_batchmode_lisa:_maybe_l0_rescue:7a4e161eb0": { + "source": "_cold_rvs = dict(sampler._rvs)", + "verdict": "PER_ROW", + "why": "Snapshots the cold record so the reject path can restore it. A dict copy of whatever rows exist; makes no claim about their statistics." + }, + "bin/integrate_likelihood_extrinsic_batchmode_lisa:_maybe_l0_rescue:a6ea5209f8": { + "source": "_cols = (np.vstack([np.asarray(sampler.identity_convert(sampler._rvs[p]), dtype=float).ravel()", + "verdict": "BENIGN", + "why": "The DELIBERATE no-reserve fallback for a sampler that keeps none: reads _rvs so the feature degrades to its previous behaviour rather than to no seed at all. The rank hazard is still handled -- build_warm_seed puffs the result to full rank -- so what remains is only 'fewer points', which cannot be improved without a reserve." + }, + "bin/integrate_likelihood_extrinsic_batchmode_lisa:_maybe_l0_rescue:cfe2518e26": { + "source": "_lnkey = 'log_integrand' if 'log_integrand' in sampler._rvs else ('integrand' if 'integrand' in sampler._rvs else None)", + "verdict": "PER_ROW", + "why": "Key-presence / column reference, not a population statistic." + }, + "bin/integrate_likelihood_extrinsic_batchmode_lisa:_maybe_l0_rescue:d2514d9663": { + "source": "_warm_lnZ, _warm_src = _lnZ_of_reserve_or_rvs(sampler, sampler._rvs)", + "verdict": "FIXED", + "why": "PR #79, re-landed as #86. sampler._rvs is passed as the FALLBACK argument only: the helper prefers the retained reserve via lnZ_from_reserve, and the gate refuses to compare across sources (_cold_src != _warm_src forces BOTH back to the fair-draw reading, which is at least self-consistent). Measured before #79: two passes with identical true lnZ at n_eff 1.8 vs 53 produced a +3.48 nat gap and rejected the good warm pass 100% of the time at the 0.5 default." + }, + "bin/integrate_likelihood_extrinsic_batchmode_lisa:_snapshot_pass_state:02594db9f0": { + "source": "rvs=(dict(sampler._rvs) if rvs is None else rvs),", + "verdict": "PER_ROW", + "why": "_snapshot_pass_state takes a dict copy of whatever rows are present so a rejected warm pass can be undone. It makes no claim about their statistics, and it snapshots the fair-draw MARKER alongside them so the restored record and the marker describing it cannot disagree." + }, "bin/integrate_likelihood_extrinsic_batchmode_lisa:analyze_event:16e8b48c86": { "source": "(sampler._rvs[\"phi_orb\"][indx_guess]/(2*numpy.pi)), \\", "verdict": "BENIGN", diff --git a/MonteCarloMarginalizeCode/Code/test/test_lisa_av_state.py b/MonteCarloMarginalizeCode/Code/test/test_lisa_av_state.py new file mode 100644 index 000000000..e578b242f --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_lisa_av_state.py @@ -0,0 +1,328 @@ +#!/usr/bin/env python +""" +Tests for the AV live-volume state, per-axis bin allocation and collapse gate ported into +the LISA ILE driver (bin/integrate_likelihood_extrinsic_batchmode_lisa). + +Four options, all sampler-agnostic: --sampler-save-state / --sampler-load-state (the AV +grid, which carries no detector convention), --sampler-anisotropic-bins, and +--reject-collapsed-live-volume. + +THE ONE THING TO KNOW. The main driver calls its collapse gate TWICE -- once on the first +run, and again on the replica POOL, because replication can turn a healthy first run into a +collapsed pool. This driver has no replica pooling yet, so only the first call exists here. +When --mc-error-replicas is ported the second call MUST come with it, or the flag is +silently bypassed for exactly the case pooling introduces. That is recorded at the helper, +in the drift ledger, and asserted below. +""" + +import ast +import os + +import pytest + +_HERE = os.path.dirname(os.path.abspath(__file__)) +_LISA = os.path.join(_HERE, '..', 'bin', 'integrate_likelihood_extrinsic_batchmode_lisa') +_MAIN = os.path.join(_HERE, '..', 'bin', 'integrate_likelihood_extrinsic_batchmode') + +OPTS = ["--sampler-save-state", "--sampler-load-state", + "--sampler-anisotropic-bins", "--reject-collapsed-live-volume"] + +HELPERS = ['_maybe_load_av_state', '_maybe_save_av_state', + '_maybe_enable_anisotropic_bins', '_reject_if_collapsed', + '_report_and_gate_collapse'] + + +def _src(path): + with open(path) as fh: + return fh.read() + + +def _option_nodes(path): + out = {} + for n in ast.walk(ast.parse(_src(path), filename=path)): + if (isinstance(n, ast.Call) and isinstance(n.func, ast.Attribute) + and n.func.attr in ("add_option", "add_argument")): + names = [a.value for a in n.args + if isinstance(a, ast.Constant) and isinstance(a.value, str)] + if names and names[0].startswith("--"): + out[names[0]] = n + return out + + +def _kwargs_of(node): + out = {} + for kw in node.keywords: + try: + out[kw.arg] = ast.literal_eval(kw.value) + except Exception: + out[kw.arg] = ast.dump(kw.value) + return out + + +class _Collapse(Exception): + pass + + +class _AVModule(object): + LiveVolumeCollapse = _Collapse + + +def _load(**optkw): + base = {"sampler_load_state": None, "sampler_save_state": None, + "sampler_anisotropic_bins": False, "reject_collapsed_live_volume": False, + "sampler_method": "AV"} + base.update(optkw) + defs = {n.name: n for n in ast.parse(_src(_LISA)).body + if isinstance(n, ast.FunctionDef) and n.name in HELPERS} + missing = sorted(set(HELPERS) - set(defs)) + assert not missing, "LISA driver is missing: %s" % missing + mod = ast.Module(body=[defs[n] for n in HELPERS], type_ignores=[]) + ns = {"opts": type("O", (), base)(), + "mcsamplerAdaptiveVolume": _AVModule, "mcsampler_AV_ok": True} + exec(compile(ast.fix_missing_locations(mod), "av_state", "exec"), ns) + return ns + + +# ------------------------------------------------------------------------------- options +@pytest.mark.parametrize("opt", OPTS) +def test_option_present_and_matches_the_main_driver(opt): + a, b = _kwargs_of(_option_nodes(_LISA)[opt]), _kwargs_of(_option_nodes(_MAIN)[opt]) + for key in ("default", "type", "action", "choices"): + assert a.get(key) == b.get(key), "%s: %s differs" % (opt, key) + + +# ---------------------------------------------------------------------------- load / save +class _AV(object): + def __init__(self): + self.loaded = self.saved = None + + def load_state(self, p): + self.loaded = p + + def save_state(self, p): + self.saved = p + + +class _NoState(object): + pass + + +def test_state_hooks_are_noops_when_unset(): + ns = _load() + s = _AV() + ns['_maybe_load_av_state'](s) + ns['_maybe_save_av_state'](s) + assert s.loaded is None and s.saved is None + + +def test_state_round_trip_reaches_the_sampler(): + ns = _load(sampler_load_state="/in.npz", sampler_save_state="/out.npz") + s = _AV() + ns['_maybe_load_av_state'](s) + ns['_maybe_save_av_state'](s) + assert s.loaded == "/in.npz" and s.saved == "/out.npz" + + +def test_save_state_is_restricted_to_the_AV_method(): + """Main gates the save on sampler_method == 'AV'; a portfolio's aggregate has no such grid.""" + ns = _load(sampler_save_state="/out.npz", sampler_method="portfolio") + s = _AV() + ns['_maybe_save_av_state'](s) + assert s.saved is None + + +def test_rejected_or_failed_rescue_state_is_not_saved(capsys): + """The warm grid may outlive restoration of the cold result; never persist that mismatch.""" + ns = _load(sampler_save_state="/out.npz") + s = _AV() + s._av_state_reuse_safe = False + ns['_maybe_save_av_state'](s) + assert s.saved is None + assert "not saving" in capsys.readouterr().out + + +def test_state_hooks_tolerate_a_sampler_without_state_support(): + ns = _load(sampler_load_state="/in.npz", sampler_save_state="/out.npz") + ns['_maybe_load_av_state'](_NoState()) + ns['_maybe_save_av_state'](_NoState()) + + +def test_a_bad_state_file_degrades_to_a_cold_run(): + """A missing/corrupt state must not kill the point.""" + ns = _load(sampler_load_state="/in.npz", sampler_save_state="/out.npz") + + class _Boom(object): + def load_state(self, p): + raise IOError("nope") + + def save_state(self, p): + raise IOError("read-only") + + ns['_maybe_load_av_state'](_Boom()) + ns['_maybe_save_av_state'](_Boom()) + + +# --------------------------------------------------------------------------- anisotropic +class _Binned(object): + anisotropic_bins = False + + +def test_anisotropic_bins_is_opt_in(): + ns = _load() + s = _Binned() + ns['_maybe_enable_anisotropic_bins'](s) + assert s.anisotropic_bins is False + + +def test_anisotropic_bins_reaches_portfolio_members_too(): + """The grid lives on the MEMBERS; setting it only on the aggregate would do nothing.""" + ns = _load(sampler_anisotropic_bins=True) + m1, m2 = _Binned(), _Binned() + s = _Binned() + s.portfolio_realizations = [m1, m2] + ns['_maybe_enable_anisotropic_bins'](s) + assert s.anisotropic_bins and m1.anisotropic_bins and m2.anisotropic_bins + + +def test_anisotropic_bins_skips_members_that_do_not_support_it(): + ns = _load(sampler_anisotropic_bins=True) + s = _Binned() + s.portfolio_realizations = [_NoState()] + ns['_maybe_enable_anisotropic_bins'](s) # must not raise + assert s.anisotropic_bins is True + + +# -------------------------------------------------------------------------- collapse gate +COLLAPSED = {'live_volume_collapsed': True, 'collapse_reason': 'zero volume'} +HEALTHY = {'live_volume_collapsed': False} + + +def test_gate_is_inert_when_the_flag_is_off(): + _load()['_reject_if_collapsed'](COLLAPSED, "first run") + + +def test_gate_is_inert_on_a_healthy_run(): + _load(reject_collapsed_live_volume=True)['_reject_if_collapsed'](HEALTHY, "first run") + + +def test_gate_raises_when_flag_set_and_run_collapsed(): + with pytest.raises(_Collapse): + _load(reject_collapsed_live_volume=True)['_reject_if_collapsed'](COLLAPSED, "first run") + + +def test_gate_message_names_the_stage_and_reason(): + """The stage is in the message because the main driver calls this at two stages.""" + with pytest.raises(_Collapse) as e: + _load(reject_collapsed_live_volume=True)['_reject_if_collapsed'](COLLAPSED, "pooled") + assert "pooled" in str(e.value) and "zero volume" in str(e.value) + + +@pytest.mark.parametrize("dd", [None, "not a dict", {}]) +def test_gate_tolerates_a_missing_or_malformed_dict_return(dd): + _load(reject_collapsed_live_volume=True)['_reject_if_collapsed'](dd, "first run") + + +def test_report_announces_a_collapse_even_when_the_gate_is_off(capsys): + """Not rejecting is not the same as not telling anyone.""" + ns = _load() + capsys.readouterr() + ns['_report_and_gate_collapse'](COLLAPSED) + out = capsys.readouterr().out + assert "LIVE VOLUME COLLAPSED" in out and "NOT a fair draw" in out + + +def test_report_says_nothing_on_a_healthy_run(capsys): + ns = _load() + capsys.readouterr() + ns['_report_and_gate_collapse'](HEALTHY) + assert "COLLAPSED" not in capsys.readouterr().out + + +def test_report_still_raises_when_gated(): + with pytest.raises(_Collapse): + _load(reject_collapsed_live_volume=True)['_report_and_gate_collapse'](COLLAPSED) + + +def test_gate_falls_back_to_RuntimeError_without_AV(): + """mcsampler_AV_ok False -> the AV exception class is unavailable.""" + defs = {n.name: n for n in ast.parse(_src(_LISA)).body + if isinstance(n, ast.FunctionDef) and n.name == '_reject_if_collapsed'} + mod = ast.Module(body=[defs['_reject_if_collapsed']], type_ignores=[]) + ns = {"opts": type("O", (), {"reject_collapsed_live_volume": True})(), + "mcsamplerAdaptiveVolume": None, "mcsampler_AV_ok": False} + exec(compile(ast.fix_missing_locations(mod), "av_state", "exec"), ns) + with pytest.raises(RuntimeError): + ns['_reject_if_collapsed'](COLLAPSED, "first run") + + +# ------------------------------------------------------------------------- call-site wiring +def test_both_analyze_event_variants_get_every_hook(): + tree = ast.parse(_src(_LISA)) + fns = {n.name: n for n in tree.body + if isinstance(n, ast.FunctionDef) and n.name in ('analyze_event', 'analyze_event_LISA')} + assert set(fns) == {'analyze_event', 'analyze_event_LISA'} + for name, node in fns.items(): + called = {c.func.id for c in ast.walk(node) + if isinstance(c, ast.Call) and isinstance(c.func, ast.Name)} + for hook in ('_maybe_load_av_state', '_maybe_save_av_state', + '_maybe_enable_anisotropic_bins', '_report_and_gate_collapse'): + assert hook in called, "%s does not call %s" % (name, hook) + + +def test_hook_ordering_at_both_call_sites(): + """Only a nonempty, collapse-approved result may persist its live-volume state.""" + src = _src(_LISA) + pos = 0 + for _ in range(2): + load = src.index("_maybe_load_av_state(sampler)", pos) + aniso = src.index("_maybe_enable_anisotropic_bins(sampler)", load) + integ = src.index("sampler.integrate(like_to_integrate", aniso) + guard = src.index("if not(res): # no resut", integ) + gate = src.index("_report_and_gate_collapse(dict_return", guard) + save = src.index("_maybe_save_av_state(sampler)", gate) + assert load < aniso < integ < guard < gate < save + pos = save + 1 + + +def test_the_second_gate_call_site_is_recorded_as_missing(): + """Main gates twice; this driver gates once because it has no replica pooling yet. + + If someone ports --mc-error-replicas without adding the second call, the flag is + silently bypassed for the case pooling creates. This asserts the warning is still + written down where that person will be working. + """ + src = _src(_LISA) + fn = src[src.index("def _reject_if_collapsed"):] + fn = fn[:fn.index("\ndef ")] + assert "mc-error-replicas" in fn and "TWICE" in fn + # Count CALLS, not textual occurrences: the `def` line matches the same substring. + calls = [c for c in ast.walk(ast.parse(src)) + if isinstance(c, ast.Call) and isinstance(c.func, ast.Name) + and c.func.id == '_report_and_gate_collapse'] + assert len(calls) == 2, \ + "expected exactly one gate call per analyze_event variant, found %d" % len(calls) + + +# ---------------------------------------------------------------- anti-drift vs the main driver +def _named(path, name): + for n in ast.walk(ast.parse(_src(path))): + if isinstance(n, ast.FunctionDef) and n.name == name: + return n + raise AssertionError("%s not found in %s" % (name, os.path.basename(path))) + + +def _normalized(fn): + node = ast.parse(ast.unparse(fn)).body[0] if hasattr(ast, "unparse") else fn + body = list(node.body) + if (body and isinstance(body[0], ast.Expr) + and isinstance(getattr(body[0], "value", None), ast.Constant) + and isinstance(body[0].value.value, str)): + body = body[1:] + return ast.dump(ast.fix_missing_locations(ast.Module(body=body, type_ignores=[]))) + + +def test_reject_if_collapsed_body_is_identical_to_the_main_drivers(): + """Hoisted out of analyze_event here, but the body must not have changed with it.""" + assert (_normalized(_named(_LISA, '_reject_if_collapsed')) + == _normalized(_named(_MAIN, '_reject_if_collapsed'))), \ + "_reject_if_collapsed has drifted between the two drivers (docstrings excluded)" diff --git a/MonteCarloMarginalizeCode/Code/test/test_lisa_driver_drift.py b/MonteCarloMarginalizeCode/Code/test/test_lisa_driver_drift.py new file mode 100644 index 000000000..52fe0d0be --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_lisa_driver_drift.py @@ -0,0 +1,158 @@ +#!/usr/bin/env python +""" +Drift gate for the LISA ILE driver. + +The two ILE drivers are a DELIBERATE fork: + + bin/integrate_likelihood_extrinsic_batchmode <- main, moves fast + bin/integrate_likelihood_extrinsic_batchmode_lisa <- LISA, lags + +RO, 2026-08-13: "It is super annoying we have to have two of them, but the overhead of one +ring to rule them all is too high." So this gate does NOT try to close the gap, and does +not assert that any particular item was ported. Closing the gap is not the goal. + +What it asserts is that nothing drifts in UNNOTICED: every helper, CLI option, module +constant and sampler provenance marker present in the main driver and absent from the LISA +one carries a recorded decision -- PORT / PORTED / NA / PHYSICS -- with a reason. "Does not +apply to LISA" is a fine answer; silence is not. + +When this fails, the fix is to classify the new item, not to delete the test: + + cd test/expensive_before_merging/integrators + python3 audit_lisa_driver_drift.py --undecided # what is unclassified + $EDITOR make_lisa_drift_ledger.py # add a rule, with a reason + python3 make_lisa_drift_ledger.py # regenerate the ledger + +This exists because 2,357 lines of drift accumulated while the LISA driver's nine CI tests +(all import/contract/smoke level) stayed green. +""" + +import os +import sys + +import pytest + +_HERE = os.path.dirname(os.path.abspath(__file__)) +_AUDIT_DIR = os.path.join(_HERE, 'expensive_before_merging', 'integrators') + +if _AUDIT_DIR not in sys.path: + sys.path.insert(0, _AUDIT_DIR) + +audit = pytest.importorskip("audit_lisa_driver_drift", + reason="LISA drift auditor not present") + + +@pytest.fixture(scope="module") +def state(): + gap, extras = audit.compute_gap() + ledger = audit.load_ledger() + return audit.annotate(gap, ledger), extras, ledger + + +def test_the_gap_is_non_empty_so_the_audit_is_actually_looking(state): + """Guard against a silently broken extractor reporting a clean tree.""" + gap, _extras, _ledger = state + assert len(gap) > 0, "the audit found no drift at all, which almost certainly means " \ + "the extractor broke rather than that the drivers converged" + + +def test_every_gap_item_carries_a_recorded_decision(state): + gap, _extras, _ledger = state + undecided = [g for g in gap if g["decision"] is None] + assert not undecided, ( + "%d item(s) drifted into the main ILE driver with no recorded decision about the " + "LISA driver:\n%s\n\nClassify each as PORT / PORTED / NA / PHYSICS with a reason " + "in make_lisa_drift_ledger.py, then regenerate the ledger." + % (len(undecided), "\n".join(" %s (main:%d)" % (g["key"], g["main_line"]) + for g in undecided))) + + +def test_no_item_claims_to_be_ported_while_still_missing(state): + """A PORTED verdict is a claim about the tree, so the tree gets to contradict it. + + This is the regression direction: if a ported helper is later deleted from the LISA + driver, the item reappears in the gap still marked PORTED, and this fails. + """ + gap, _extras, _ledger = state + stale = [g for g in gap if g["decision"] == "PORTED"] + assert not stale, ( + "marked PORTED but absent from the LISA driver: %s" + % ", ".join(g["key"] for g in stale)) + + +def test_every_decision_is_a_known_verdict(state): + gap, _extras, _ledger = state + bad = sorted({g["decision"] for g in gap + if g["decision"] is not None and g["decision"] not in audit.DECISIONS}) + assert not bad, "unknown decision value(s) in the ledger: %s" % bad + + +def test_every_decision_carries_a_reason(state): + """A verdict without a reason is silence with extra steps.""" + gap, _extras, _ledger = state + thin = [g["key"] for g in gap + if g["decision"] is not None and len((g["reason"] or "").strip()) < 20] + assert not thin, "decision recorded with no usable reason: %s" % ", ".join(thin) + + +def test_ledger_has_no_entries_for_items_outside_the_gap(state): + """Spent entries are not a failure, but they should not pile up as fiction. + + An entry naming something no longer in the gap means it was ported or the main driver + dropped it; regenerating the ledger clears it. + """ + gap, _extras, ledger = state + gap_keys = {g["key"] for g in gap} + spent = sorted(k for k in ledger if k not in gap_keys) + assert not spent, ("ledger describes %d item(s) that are no longer in the gap: %s\n" + "Regenerate with make_lisa_drift_ledger.py." + % (len(spent), ", ".join(spent))) + + +def test_the_committed_ledger_matches_what_its_generator_produces(): + """The ledger is GENERATED. Nothing enforced that until this test. + + An adversarial audit added an option to the main driver and hand-wrote a + ``{"decision": "NA", "reason": "..."}`` entry straight into the JSON: the whole gate + passed while make_lisa_drift_ledger.py still reported the item as matching no rule. + The stated property -- that a person has to classify new drift AS A RULE, with a reason + -- was silenceable by a one-line JSON edit. + + So regenerate in memory and compare. This also catches a ledger left stale after the + main driver moved. + """ + gen = pytest.importorskip("make_lisa_drift_ledger", + reason="LISA drift ledger generator not present") + gap, _extras = audit.compute_gap() + expected, unmatched = {}, [] + for item in gap: + decision, reason = gen.classify(item["key"]) + if decision is None: + unmatched.append(item["key"]) + else: + expected[item["key"]] = {"decision": decision, "reason": reason} + + assert not unmatched, ( + "%d gap item(s) match no rule in make_lisa_drift_ledger.py: %s\n" + "Add a rule with a reason -- do not hand-edit the JSON." + % (len(unmatched), ", ".join(unmatched))) + + committed = audit.load_ledger() + assert committed == expected, ( + "lisa_drift_ledger.json does not match make_lisa_drift_ledger.py.\n" + "Regenerate it (python3 make_lisa_drift_ledger.py) rather than editing the JSON:\n" + " only in committed: %s\n only in generated: %s\n differing: %s" + % (sorted(set(committed) - set(expected)), + sorted(set(expected) - set(committed)), + sorted(k for k in set(committed) & set(expected) if committed[k] != expected[k]))) + + +def test_the_fairdraw_helpers_ported_in_this_pass_are_present_in_lisa(): + """Belt and braces: name them, so deleting one fails here as well as via the ledger.""" + lisa = audit.collect(audit.LISA) + for name in ('ln_weights_from_rvs', 'ln_weights_for_posterior', + '_rvs_is_export_resample', '_rvs_is_equal_weight', '_rvs_len', + '_rvs_lnL_convention'): + assert name in lisa["FUNC"], "%s is missing from the LISA driver" % name + for marker in ('_rvs_is_fairdraw', '_rvs_is_pooled'): + assert marker in lisa["ATTR"], "%s is no longer read by the LISA driver" % marker diff --git a/MonteCarloMarginalizeCode/Code/test/test_lisa_fairdraw_weights.py b/MonteCarloMarginalizeCode/Code/test/test_lisa_fairdraw_weights.py new file mode 100644 index 000000000..8d8f947cc --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_lisa_fairdraw_weights.py @@ -0,0 +1,305 @@ +#!/usr/bin/env python +""" +Tests for the fair-draw weighting helpers ported into the LISA ILE driver +(bin/integrate_likelihood_extrinsic_batchmode_lisa) from the main driver, PR #87. + +WHY THE LISA DRIVER NEEDS THEM AT ALL. The three consumers whose double-weighting PR #87 +fixed -- the `--extrinsic-proposal-output` breadcrumb, the `.dgrid` exporter and the +`.dslice` reweight core -- do not exist in the LISA driver, so there is no live w^2 bug +there today. What DOES exist is the hazard: the LISA driver sets +`igrand_fairdraw_samples` from `--fairdraw-extrinsic-output`, so its `_rvs` can be a fair +draw, and every shared sampler in RIFT/integrators/ already sets `_rvs_is_fairdraw` at its +rebind. The marker was arriving and nothing read it. These tests pin the readers. + +TWO DISTINCT PROPERTIES, deliberately not one flag (audit Finding 6): + + rows resampled -- each row drawn proportional to w (per-BLOCK property) + equal weight -- the record as a whole is uniform (property of the WHOLE record) + +and the anti-drift test at the bottom pins the LISA copies to the main driver's, because +these are deliberate COPIES in a deliberate fork, not an import. + +Conventions follow test_fairdraw_double_weighting.py and test_l0_rescue_seed.py: the driver +scripts are not importable (they parse argv at import), so the helpers are exec'd out. +""" + +import ast +import os + +import numpy as np +import pytest + +_HERE = os.path.dirname(os.path.abspath(__file__)) +_LISA = os.path.join(_HERE, '..', 'bin', 'integrate_likelihood_extrinsic_batchmode_lisa') +_MAIN = os.path.join(_HERE, '..', 'bin', 'integrate_likelihood_extrinsic_batchmode') + +# The helpers ported in this pass. Named explicitly: if a future edit drops one, the +# extraction below fails loudly rather than silently testing a smaller surface. +PORTED = ['_rvs_lnL_convention', 'ln_weights_from_rvs', '_rvs_len', + '_rvs_is_export_resample', '_rvs_is_equal_weight', 'ln_weights_for_posterior'] + + +def _extract(path, names): + """Return {name: ast.FunctionDef} for top-level defs, by name.""" + with open(path) as fh: + tree = ast.parse(fh.read(), filename=path) + found = {n.name: n for n in tree.body + if isinstance(n, ast.FunctionDef) and n.name in names} + missing = sorted(set(names) - set(found)) + assert not missing, "%s is missing ported helper(s): %s" % (os.path.basename(path), missing) + return found + + +def _load(path, names=PORTED): + """Exec the named helpers out of a driver script into a namespace.""" + defs = _extract(path, names) + mod = ast.Module(body=[defs[n] for n in names], type_ignores=[]) + ns = {"numpy": np, "np": np} + exec(compile(ast.fix_missing_locations(mod), "lisa_weight_helpers", "exec"), ns) + return ns + + +@pytest.fixture(scope="module") +def H(): + return _load(_LISA) + + +# --------------------------------------------------------------------------- record builders +def _log_record(n=6, seed=0): + rng = np.random.default_rng(seed) + return {'log_integrand': rng.normal(size=n) * 3.0, + 'log_joint_prior': rng.normal(size=n), + 'log_joint_s_prior': rng.normal(size=n), + 'right_ascension': rng.uniform(0, 2 * np.pi, size=n)} + + +def _linear_record(n=6, seed=1, lnL=False): + rng = np.random.default_rng(seed) + ig = (rng.normal(size=n) * 3.0) if lnL else rng.uniform(0.1, 5.0, size=n) + return {'integrand': ig, + 'joint_prior': rng.uniform(0.1, 2.0, size=n), + 'joint_s_prior': rng.uniform(0.1, 2.0, size=n), + 'psi': rng.uniform(0, np.pi, size=n)} + + +class _FakeSampler(object): + def __init__(self, fairdraw=None, pooled=None): + if fairdraw is not None: + self._rvs_is_fairdraw = fairdraw + if pooled is not None: + self._rvs_is_pooled = pooled + + +# ------------------------------------------------------------------- ln_weights_from_rvs +def test_log_form_is_the_canonical_combination(H): + r = _log_record() + got = H['ln_weights_from_rvs'](r) + want = r['log_integrand'] + r['log_joint_prior'] - r['log_joint_s_prior'] + assert np.allclose(got, want) + + +def test_log_form_preferred_over_linear_when_both_present(H): + """The log columns win. A record carrying both must not be read the linear way.""" + r = _log_record() + r.update({'integrand': np.full(len(r['log_integrand']), 1.0), + 'joint_prior': np.full(len(r['log_integrand']), 1.0), + 'joint_s_prior': np.full(len(r['log_integrand']), 1.0)}) + got = H['ln_weights_from_rvs'](r) + want = r['log_integrand'] + r['log_joint_prior'] - r['log_joint_s_prior'] + assert np.allclose(got, want), "linear columns shadowed the canonical log ones" + + +def test_linear_form_linear_convention(H): + r = _linear_record(lnL=False) + got = H['ln_weights_from_rvs'](r, use_lnL=False) + want = np.log(r['integrand']) + np.log(r['joint_prior']) - np.log(r['joint_s_prior']) + assert np.allclose(got, want) + + +def test_linear_form_out_of_support_rows_are_minus_inf(H): + r = _linear_record(lnL=False) + r['joint_prior'][2] = 0.0 # zero prior -> out of support + r['integrand'][4] = 0.0 # zero L -> out of support + got = H['ln_weights_from_rvs'](r, use_lnL=False) + assert got[2] == -np.inf and got[4] == -np.inf + assert np.isfinite(got[[0, 1, 3, 5]]).all() + + +def test_lnL_convention_does_not_log_twice_and_keeps_negative_lnL(H): + """The bug this argument exists for. + + mcsamplerEnsemble reuses 'integrand' for BOTH conventions. Under return_lnI it holds + lnL, so the linear reading would (a) take log() of it, compressing tens of nats into + log(tens), and (b) apply `ig > 0`, silently discarding every sample with lnL <= 0. + """ + r = _linear_record(lnL=True) + r['integrand'][0] = -12.5 # a perfectly good low-likelihood point + got = H['ln_weights_from_rvs'](r, use_lnL=True) + want = r['integrand'] + np.log(r['joint_prior']) - np.log(r['joint_s_prior']) + assert np.allclose(got, want) + assert np.isfinite(got[0]), "a negative lnL row was discarded as out-of-support" + + wrong = H['ln_weights_from_rvs'](r, use_lnL=False) + assert not np.allclose(np.nan_to_num(wrong, neginf=-1e9), got), \ + "the two conventions agree, so this test cannot detect reading lnL as L" + + +def test_raises_when_neither_component_set_is_present(H): + """An explicit failure beats a plausible wrong number.""" + with pytest.raises(Exception): + H['ln_weights_from_rvs']({'psi': np.zeros(4), 'log_weights': np.zeros(4)}) + + +def test_cached_log_weights_column_is_never_read(H): + """mcsamplerGPU stores the ADAPTATION weight there, with adapt-weight-exponent baked in.""" + r = _log_record() + r['log_weights'] = np.full(len(r['log_integrand']), 999.0) + got = H['ln_weights_from_rvs'](r) + assert not np.allclose(got, 999.0) + + +# ------------------------------------------------------------------------ the two predicates +@pytest.mark.parametrize("fairdraw,pooled,resample,equal", [ + (None, None, False, False), # markers absent entirely -> both False, no AttributeError + (False, False, False, False), + (True, False, True, True), # a plain fair draw has BOTH properties + (True, True, True, False), # pooled: rows resampled, record NOT globally uniform + (False, True, False, False), +]) +def test_predicate_truth_table(H, fairdraw, pooled, resample, equal): + s = _FakeSampler(fairdraw, pooled) + assert H['_rvs_is_export_resample'](s) is resample + assert H['_rvs_is_equal_weight'](s) is equal + + +def test_predicates_differ_on_a_pooled_record(H): + """The Finding-6 property: one flag cannot answer both questions.""" + s = _FakeSampler(fairdraw=True, pooled=True) + assert H['_rvs_is_export_resample'](s) != H['_rvs_is_equal_weight'](s) + + +# --------------------------------------------------------------- ln_weights_for_posterior +def test_fair_drawn_record_gets_uniform_posterior_weights(H): + """The anti-double-weighting property: rows already ~w must not be weighted by w again.""" + r = _log_record() + w = H['ln_weights_for_posterior'](r, _FakeSampler(fairdraw=True, pooled=False)) + assert w.shape == (len(r['log_integrand']),) + assert np.allclose(w, 0.0) + + +def test_non_fairdrawn_record_gets_the_importance_weights(H): + r = _log_record() + s = _FakeSampler(fairdraw=False, pooled=False) + assert np.allclose(H['ln_weights_for_posterior'](r, s), H['ln_weights_from_rvs'](r)) + + +def test_pooled_record_keeps_its_between_block_weights(H): + """Pooling weights block k by the replica evidence: uniform here would discard that.""" + r = _log_record() + w = H['ln_weights_for_posterior'](r, _FakeSampler(fairdraw=True, pooled=True)) + assert not np.allclose(w, 0.0) + assert np.allclose(w, H['ln_weights_from_rvs'](r)) + + +def test_double_weighting_would_shift_a_posterior_mean(H): + """Why it matters, not just that it differs. + + Build a record whose weight correlates with a coordinate, fair-draw it, then compare the + mean under the correct (uniform) weights against the mean under a second application of + w. The second application concentrates toward high-w rows and moves the answer. + """ + rng = np.random.default_rng(7) + n = 4000 + x = rng.uniform(0.0, 1.0, size=n) + lnw = 4.0 * x # weight correlated with the coordinate + w = np.exp(lnw - lnw.max()) + idx = rng.choice(n, size=n, replace=True, p=w / w.sum()) # the fair draw + rec = {'log_integrand': lnw[idx], 'log_joint_prior': np.zeros(n), + 'log_joint_s_prior': np.zeros(n), 'x': x[idx]} + + correct = H['ln_weights_for_posterior'](rec, _FakeSampler(fairdraw=True, pooled=False)) + assert np.allclose(correct, 0.0) + mean_correct = np.average(rec['x'], weights=np.exp(correct - correct.max())) + + doubled = H['ln_weights_from_rvs'](rec) # what the pre-fix consumers did + mean_doubled = np.average(rec['x'], weights=np.exp(doubled - doubled.max())) + + shift = abs(mean_doubled - mean_correct) / abs(mean_correct) + assert shift > 0.05, ("double weighting should move the posterior mean materially; " + "got %.3f%%" % (100 * shift)) + + +# ----------------------------------------------------------------------------- _rvs_len +def test_rvs_len_counts_rows(H): + assert H['_rvs_len'](_log_record(n=9)) == 9 + + +def test_rvs_len_survives_an_unsized_entry(H): + r = _log_record(n=5) + r['not_an_array'] = None + assert H['_rvs_len'](r) == 5 + + +# ------------------------------------------------------------------ the convention resolver +def test_lnL_convention_prefers_the_explicit_argument(H): + assert H['_rvs_lnL_convention'](True) is True + assert H['_rvs_lnL_convention'](False) is False + + +def test_lnL_convention_falls_back_to_linear_outside_the_driver(H): + """No `rvs_integrand_is_lnL` in scope (which is the case in these tests) -> False.""" + assert H['_rvs_lnL_convention'](None) is False + + +# ------------------------------------------------------------- source-level wiring in LISA +def _lisa_src(): + with open(_LISA) as fh: + return fh.read() + + +def test_lisa_derives_the_convention_from_pinned_params_not_the_cli_option(): + """The trap this port had to avoid. + + --internal-use-lnL is ALSO accepted for adaptive_cartesian_gpu and portfolio, and those + branches set use_lnL WITHOUT return_lnI -- they still store linear L. Deriving the + stored convention from the option would read those records as lnL. + """ + src = _lisa_src() + assert 'rvs_integrand_is_lnL = bool(pinned_params.get("return_lnI", False))' in src, \ + "the stored-integrand convention is not derived from pinned_params['return_lnI']" + assert 'rvs_integrand_is_lnL = bool(opts.internal_use_lnL' not in src, \ + "the convention is keyed off the CLI option, which is a different predicate" + + +def test_lisa_still_requests_the_fair_draw(): + """If this ever stops being set, the helpers become dead code and should be revisited.""" + assert '"igrand_fairdraw_samples": opts.fairdraw_extrinsic_output' in _lisa_src() + + +# ------------------------------------------------------------------ anti-drift vs the main driver +def _normalized(fn): + """AST dump of a function with its docstring stripped. + + Docstrings are deliberately allowed to differ -- the LISA copies carry LISA-specific + notes. Everything the interpreter runs must match. + """ + node = ast.parse(ast.unparse(fn)).body[0] if hasattr(ast, "unparse") else fn + body = list(node.body) + if (body and isinstance(body[0], ast.Expr) + and isinstance(getattr(body[0], "value", None), ast.Constant) + and isinstance(body[0].value.value, str)): + body = body[1:] + stripped = ast.Module(body=body, type_ignores=[]) + return ast.dump(ast.fix_missing_locations(stripped)) + + +@pytest.mark.parametrize("name", PORTED) +def test_ported_helper_is_identical_to_the_main_driver(name): + """These are COPIES in a deliberate fork. A copy that quietly changes is the whole risk. + + If you intend to change one, change both -- or record the divergence explicitly. + """ + lisa = _extract(_LISA, [name])[name] + main = _extract(_MAIN, [name])[name] + assert _normalized(lisa) == _normalized(main), ( + "%s has drifted between the two drivers (docstrings excluded)" % name) diff --git a/MonteCarloMarginalizeCode/Code/test/test_lisa_l0_rescue.py b/MonteCarloMarginalizeCode/Code/test/test_lisa_l0_rescue.py new file mode 100644 index 000000000..ca3793ccf --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_lisa_l0_rescue.py @@ -0,0 +1,609 @@ +#!/usr/bin/env python +""" +Tests for the L0 auto-rescue ported into the LISA ILE driver +(bin/integrate_likelihood_extrinsic_batchmode_lisa) from the main driver. + +WHY IT BELONGS IN LISA. The rescue targets the high-SNR n_eff LOTTERY: a large fraction of +independent AV/portfolio runs collapse to n_eff ~ 1 by contracting onto the wrong spot, and +the rescue re-seeds such a run from the peak it did find. LISA MBHB are high-SNR by +construction, so this is the regime, not an edge case. The sampler-side machinery +(`build_warm_seed`, `lnZ_from_reserve`, the reserve itself) already lives in +RIFT/integrators/ and therefore already reached LISA; only the driver-side wiring was missing. + +ONE DELIBERATE STRUCTURAL DIVERGENCE. The main driver inlines the rescue in its single +`analyze_event`. This driver has TWO -- `analyze_event_LISA` (with --LISA) and +`analyze_event` (the fallback) -- so the block was lifted into `_maybe_l0_rescue` and both +call it. That is a divergence in SHAPE, not behaviour, and it buys something main does not +have: the reject gate becomes unit-testable. The audit notes that in main these call sites +"cannot be exercised from a unit test" because analyze_event needs data, PSDs and a waveform. +Here the gate is a function of its arguments, so the tests below drive it directly. + +ORDERING IS LOAD-BEARING (see test_rescue_runs_before_the_no_result_guard). In main the +`if not(res): raise` guard sits ~200 lines below the integrate call and the ordering is +implicit. In this driver it is immediately after, so the rescue had to be inserted BETWEEN +them: a degenerate early termination returns (None,None,None,None) and is the STRONGEST +rescue trigger, so raising on it first would skip exactly the case the rescue exists for. +""" + +import ast +import os + +import numpy as np +import pytest + +_HERE = os.path.dirname(os.path.abspath(__file__)) +_LISA = os.path.join(_HERE, '..', 'bin', 'integrate_likelihood_extrinsic_batchmode_lisa') +_MAIN = os.path.join(_HERE, '..', 'bin', 'integrate_likelihood_extrinsic_batchmode') + +# Helpers ported verbatim from the main driver. _maybe_l0_rescue is NOT in this list: it is +# the LISA-only wrapper, and has no counterpart to be identical to. +PORTED = ['_lnZ_of_rvs', '_kish_neff_of_rvs', '_lnZ_of_reserve_or_rvs', + '_snapshot_pass_state', '_restore_pass_state', + '_warm_seed_reserve_for', '_warm_seed_geometry', '_clear_warm_state'] + +# Everything the exec'd namespace needs, in dependency order. +_DEPS = ['_rvs_lnL_convention', 'ln_weights_from_rvs', '_rvs_len', + '_rvs_is_export_resample', '_rvs_is_equal_weight', 'ln_weights_for_posterior'] + + +def _defs(path, names): + with open(path) as fh: + tree = ast.parse(fh.read(), filename=path) + found = {n.name: n for n in tree.body + if isinstance(n, ast.FunctionDef) and n.name in names} + missing = sorted(set(names) - set(found)) + assert not missing, "%s is missing: %s" % (os.path.basename(path), missing) + return found + + +class _FakeAV(object): + """Stand-in for RIFT.integrators.mcsamplerAdaptiveVolume inside the helpers.""" + lnZ_value = None + seed_info = {'puffed': False, 'n_core': 3, 'rank_core': 3, 'dim': 3, + 'rank_final': 3, 'n_puff': 0, 'puff_scale': 'auto'} + + @classmethod + def lnZ_from_reserve(cls, reserve): + return cls.lnZ_value + + @classmethod + def build_warm_seed(cls, cols, lnL, lo, hi, axes, **kw): + return np.asarray(cols, dtype=float), dict(cls.seed_info) + + +class _Opts(object): + sampler_method = 'AV' + sampler_warmstart_retry_neff = 5.0 + sampler_l0_rescue_reject_dlnZ = 3.0 + sampler_l0_rescue_accept_truncated = False + sampler_l0_rescue_puff_scale = 'auto' + sampler_l0_rescue_puff_width_frac = 0.005 + sampler_l0_rescue_puff_factor = 2.0 + sampler_sequential_warmstart_deltalnL = 15.0 + + def __init__(self, **kw): + for k, v in kw.items(): + setattr(self, k, v) + + +def _load(opts=None, av=None): + """Exec the rescue helpers out of the LISA driver with injected globals.""" + names = _DEPS + PORTED + ['_maybe_l0_rescue'] + defs = _defs(_LISA, names) + mod = ast.Module(body=[defs[n] for n in names], type_ignores=[]) + ns = {"numpy": np, "np": np, + "opts": opts if opts is not None else _Opts(), + "mcsamplerAdaptiveVolume": av if av is not None else _FakeAV} + exec(compile(ast.fix_missing_locations(mod), "lisa_l0_helpers", "exec"), ns) + return ns + + +@pytest.fixture +def H(): + return _load() + + +# ------------------------------------------------------------------------------ fake sampler +class _Sampler(object): + def __init__(self, rvs=None, reserve=None, params=('a', 'b'), members=None, + integrate_result=None, raise_in_integrate=False): + self._rvs = rvs if rvs is not None else {} + self._warm_seed_reserve = reserve + self.params_ordered = list(params) + self.llim = {p: 0.0 for p in self.params_ordered} + self.rlim = {p: 1.0 for p in self.params_ordered} + self.portfolio_realizations = members or [] + self._warm = "stale" + self._warm_applied = True + self._integrate_result = integrate_result + self._raise_in_integrate = raise_in_integrate + self.bootstrapped = None + self.warm_rvs = None + + def identity_convert(self, x): + return x + + def bootstrap_from_samples(self, seed, cover_frac=0.0): + self.bootstrapped = (np.asarray(seed), cover_frac) + + def integrate(self, fn, *a, **kw): + if self._raise_in_integrate: + # Repopulate _rvs IN PLACE first, then raise: this is the dangerous shape -- + # the assignment at the call site never completes, so res/var/neff still hold + # the COLD pass while _rvs holds the WARM samples. + self._rvs = dict(self.warm_rvs or {}) + raise RuntimeError("warm pass exploded") + if self.warm_rvs is not None: + self._rvs = dict(self.warm_rvs) + return self._integrate_result + + +def _rec(lnL, n=None): + lnL = np.asarray(lnL, dtype=float) + n = len(lnL) if n is None else n + return {'log_integrand': lnL, + 'log_joint_prior': np.zeros(n), + 'log_joint_s_prior': np.zeros(n), + 'a': np.linspace(0.1, 0.9, n), 'b': np.linspace(0.2, 0.8, n)} + + +# ------------------------------------------------------------------------------- lnZ helpers +def test_lnZ_pooled_is_the_sum_and_unpooled_is_the_mean(H): + r = _rec([0.0, 0.0, 0.0, 0.0]) + pooled = H['_lnZ_of_rvs'](r, already_pooled=True) + single = H['_lnZ_of_rvs'](r, already_pooled=False) + assert np.isclose(pooled, np.log(4.0)) + assert np.isclose(single, 0.0) + assert np.isclose(pooled - single, np.log(4.0)) + + +def test_lnZ_returns_none_when_weights_cannot_be_rebuilt(H): + assert H['_lnZ_of_rvs']({'a': np.zeros(3)}) is None + + +def test_lnZ_ignores_non_finite_rows(H): + r = _rec([0.0, -np.inf, 0.0]) + assert np.isclose(H['_lnZ_of_rvs'](r, already_pooled=True), np.log(2.0)) + + +def test_kish_neff_of_equal_weights_is_the_row_count(H): + assert np.isclose(H['_kish_neff_of_rvs'](_rec(np.zeros(7))), 7.0) + + +def test_kish_neff_collapses_on_one_dominant_row(H): + neff = H['_kish_neff_of_rvs'](_rec([0.0, -50.0, -50.0, -50.0])) + assert 1.0 <= neff < 1.01 + + +# -------------------------------------------------------------- reserve-vs-fairdraw provenance +def test_lnZ_prefers_the_retained_reserve_and_says_so(H): + _FakeAV.lnZ_value = -1.25 + s = _Sampler(reserve={'log_joint_prior': np.zeros(3), 'log_joint_s_prior': np.zeros(3)}) + val, src = H['_lnZ_of_reserve_or_rvs'](s, _rec([0.0, 0.0])) + assert src == 'retained' and np.isclose(val, -1.25) + + +def test_lnZ_falls_back_to_the_fairdraw_record_when_no_reserve(H): + s = _Sampler(reserve=None) + val, src = H['_lnZ_of_reserve_or_rvs'](s, _rec([0.0, 0.0])) + assert src == 'fairdraw' and np.isclose(val, 0.0) + + +def test_lnZ_falls_back_when_lnZ_from_reserve_is_not_finite(H): + """Degrade to the previous behaviour, not to no gate at all.""" + _FakeAV.lnZ_value = np.nan + s = _Sampler(reserve={'log_joint_prior': np.zeros(3), 'log_joint_s_prior': np.zeros(3)}) + _val, src = H['_lnZ_of_reserve_or_rvs'](s, _rec([0.0, 0.0])) + assert src == 'fairdraw' + + +# ----------------------------------------------------------------- snapshot / restore (Finding 5) +def test_snapshot_restore_round_trips_the_whole_pass(H): + member = _Sampler(params=('a', 'b')) + member._warm_seed_reserve = {'tag': 'cold-member'} + s = _Sampler(rvs=_rec([1.0, 2.0]), reserve={'tag': 'cold'}, members=[member]) + s._rvs_is_fairdraw, s._rvs_is_pooled = True, False + + snap = H['_snapshot_pass_state'](s, 'RES', 'VAR', 'NEFF', {'d': 1}) + + # the warm pass overwrites everything + s._rvs = _rec([9.0]) + s._warm_seed_reserve = {'tag': 'WARM'} + s._rvs_is_fairdraw, s._rvs_is_pooled = False, True + member._warm_seed_reserve = {'tag': 'WARM-member'} + + out = H['_restore_pass_state'](s, snap) + assert out == ('RES', 'VAR', 'NEFF', {'d': 1}) + assert s._warm_seed_reserve == {'tag': 'cold'}, "the RESERVE did not come back (Finding 5)" + assert member._warm_seed_reserve == {'tag': 'cold-member'}, "per-member reserve did not come back" + assert s._rvs_is_fairdraw is True and s._rvs_is_pooled is False + assert np.allclose(s._rvs['log_integrand'], [1.0, 2.0]) + + +def test_snapshot_takes_a_copy_not_an_alias(H): + """integrate_log repopulates _rvs IN PLACE, so an alias would hold the warm samples.""" + s = _Sampler(rvs=_rec([1.0, 2.0])) + snap = H['_snapshot_pass_state'](s, 1, 2, 3, {}) + s._rvs['log_integrand'] = np.array([99.0, 99.0]) + assert snap['rvs'] is not s._rvs + + +# -------------------------------------------------------------------- the reserve lookup guard +def test_reserve_lookup_declines_a_column_order_mismatch(H): + """A silent mismatch produces a seed in the wrong coordinates, so decline it.""" + s = _Sampler(reserve={'params_ordered': ['b', 'a']}, params=('a', 'b')) + assert H['_warm_seed_reserve_for'](s) is None + + +def test_reserve_lookup_accepts_matching_column_order(H): + res = {'params_ordered': ['a', 'b']} + assert H['_warm_seed_reserve_for'](_Sampler(reserve=res, params=('a', 'b'))) is res + + +def test_reserve_lookup_falls_through_to_a_portfolio_member(H): + member = _Sampler(params=('a', 'b')) + member._warm_seed_reserve = {'params_ordered': ['a', 'b'], 'tag': 'member'} + s = _Sampler(reserve=None, params=('a', 'b'), members=[member]) + assert H['_warm_seed_reserve_for'](s)['tag'] == 'member' + + +# ------------------------------------------------------------------------------- seed geometry +def test_geometry_uses_the_samplers_adaptive_axes_when_it_has_them(H): + s = _Sampler(params=('a', 'b')) + s.warm_seed_axes = lambda: [1] + axes, lo, hi = H['_warm_seed_geometry'](s) + assert axes == [1] and np.allclose(lo, [0, 0]) and np.allclose(hi, [1, 1]) + + +def test_geometry_defaults_to_every_column(H): + axes, _lo, _hi = H['_warm_seed_geometry'](_Sampler(params=('a', 'b', 'c'))) + assert axes == [0, 1, 2] + + +def test_geometry_falls_through_to_a_portfolio_member(H): + member = _Sampler(params=('a', 'b')) + member.warm_seed_axes = lambda: [0] + s = _Sampler(params=('a', 'b'), members=[member]) + assert H['_warm_seed_geometry'](s)[0] == [0] + + +# ---------------------------------------------------------------------------- clearing warm state +def test_clear_warm_state_prefers_the_portfolio_hook(H): + s = _Sampler() + calls = [] + s.clear_warm_state = lambda: calls.append(1) + H['_clear_warm_state'](s) + assert calls == [1], "portfolio members would keep the previous point's contracted grid" + + +def test_clear_warm_state_falls_back_to_the_attributes(H): + s = _Sampler() + H['_clear_warm_state'](s) + assert s._warm is None and s._warm_applied is False + + +def test_clear_warm_state_does_not_swallow_failures(H): + """A reset that quietly did not happen is the silent bias this guards against.""" + s = _Sampler() + + def _boom(): + raise RuntimeError("no") + s.clear_warm_state = _boom + with pytest.raises(RuntimeError): + H['_clear_warm_state'](s) + + +# ------------------------------------------------------------------------- the rescue itself +def _run(H, sampler, res=1.0, var=0.1, neff=1.0, dict_return=None): + return H['_maybe_l0_rescue'](sampler, res, var, neff, dict_return or {'cold': True}, + lambda *a, **k: None, (), {}) + + +def _assert_declined(H, sampler, capsys, **runkw): + """The rescue must DECLINE silently -- not run and get rescued by its own except. + + Asserting only the return value is not enough, and an earlier version of these tests + made exactly that mistake: with a guard removed the rescue starts, throws somewhere + inside, and `except Exception` returns the inputs unchanged -- so the return value is + identical either way. The observable difference is that a declining rescue says + NOTHING and never touches the sampler. + """ + capsys.readouterr() + out_vals = _run(H, sampler, **runkw) + printed = capsys.readouterr().out + assert "[L0 auto-rescue]" not in printed, \ + "the rescue engaged when it should have declined: %r" % printed + assert getattr(sampler, 'bootstrapped', None) is None + return out_vals + + +def test_rescue_is_a_noop_when_the_option_is_off(capsys): + """Uses a DEGENERATE neff, so the option guard is the only thing declining. + + With neff=None, `_needs_l0_rescue` is True on its own; only the + `opts.sampler_warmstart_retry_neff` conjunct can stop the rescue here. A healthy neff + would make this test pass with that conjunct deleted. + """ + H = _load(opts=_Opts(sampler_warmstart_retry_neff=None)) + s = _Sampler(rvs=_rec([1.0, 2.0, 3.0]), integrate_result=('R2', 'V2', 42.0, {'warm': True})) + assert _assert_declined(H, s, capsys, neff=None) == (1.0, 0.1, None, {'cold': True}) + + +def test_rescue_is_a_noop_for_a_sampler_method_it_does_not_apply_to(capsys): + """AV/portfolio only. Every other conjunct is satisfied here.""" + H = _load(opts=_Opts(sampler_method='GMM')) + s = _Sampler(rvs=_rec([1.0, 2.0, 3.0]), integrate_result=('R2', 'V2', 42.0, {'warm': True})) + _assert_declined(H, s, capsys, neff=1.0) + + +def test_rescue_is_a_noop_for_a_sampler_that_cannot_warm_start(capsys): + """mcsampler/GMM have no bootstrap_from_samples; the rescue must decline, not crash.""" + H = _load() + + class _NoBootstrap(object): + def __init__(self): + self._rvs = _rec([1.0]) + self.params_ordered = ['a', 'b'] + + def identity_convert(self, x): + return x + + s = _NoBootstrap() + assert not hasattr(s, 'bootstrap_from_samples') + _assert_declined(H, s, capsys, neff=1.0) + + +def test_rescue_does_not_touch_identity_convert_before_deciding_it_applies(capsys): + """Regression: RIFT.integrators.mcsampler.MCSampler has NO identity_convert. + + That is the object this driver keeps for --sampler-method adaptive_cartesian. The main + driver evaluates `sampler.identity_convert(neff)` BEFORE its guard, so porting it + verbatim made every adaptive_cartesian event die with AttributeError at the end of a + completed integration, before --output-file was written. The applicability guards must + run first. + """ + H = _load(opts=_Opts(sampler_method='adaptive_cartesian')) + + class _NoConvert(object): + """Exactly mcsampler.MCSampler's relevant shape: no identity_convert.""" + def __init__(self): + self._rvs = _rec([1.0]) + self.params_ordered = ['a', 'b'] + + s = _NoConvert() + assert not hasattr(s, 'identity_convert') + capsys.readouterr() + assert _run(H, s, neff=1.0) == (1.0, 0.1, 1.0, {'cold': True}) + + +def test_rescue_is_a_noop_when_neff_is_healthy(capsys): + H = _load() + s = _Sampler(rvs=_rec([1.0]), integrate_result=('R2', 'V2', 42.0, {'warm': True})) + assert _assert_declined(H, s, capsys, neff=500.0)[3] == {'cold': True} + + +def test_degenerate_early_termination_triggers_the_rescue(): + """neff=None is the STRONGEST trigger, not a reason to skip.""" + H = _load() + s = _Sampler(rvs=_rec([1.0, 2.0, 3.0]), integrate_result=('R2', 'V2', 42.0, {'warm': True})) + s.warm_rvs = _rec([5.0, 5.0, 5.0]) + out = _run(H, s, neff=None) + assert s.bootstrapped is not None, "a degenerate pass did not trigger the rescue" + assert out[2] == 42.0 + + +def test_accepted_warm_pass_replaces_the_cold_result(): + H = _load() + s = _Sampler(rvs=_rec([0.0, 0.0]), integrate_result=('R2', 'V2', 42.0, {'warm': True})) + s.warm_rvs = _rec([0.0, 0.0]) # same lnZ -> no evidence of loss + out = _run(H, s) + assert out == ('R2', 'V2', 42.0, {'warm': True}) + assert s._av_state_reuse_safe is True + + +def test_warm_pass_far_below_cold_is_rejected_and_cold_is_restored(): + """The gate: positive evidence of lost mass keeps the full-support cold pass.""" + H = _load() + cold = _rec([0.0, 0.0, 0.0, 0.0]) # lnZ = 0 + s = _Sampler(rvs=cold, reserve={'tag': 'cold'}, + integrate_result=('R2', 'V2', 42.0, {'warm': True})) + s._rvs_is_fairdraw = True + s.warm_rvs = _rec([-20.0, -20.0, -20.0, -20.0]) # lnZ = -20, far below + out = _run(H, s, res='R1', var='V1', neff=1.0, dict_return={'cold': True}) + assert out == ('R1', 'V1', 1.0, {'cold': True}), "the warm pass was not rejected" + assert s._warm_seed_reserve == {'tag': 'cold'}, "the reserve did not come back (Finding 5)" + assert np.allclose(s._rvs['log_integrand'], cold['log_integrand']) + assert s._av_state_reuse_safe is False, "the rejected warm grid could be persisted" + + +def test_a_later_healthy_event_resets_the_state_save_veto(): + """Sampler objects are reused; an earlier rejection must not poison later state saves.""" + H = _load() + s = _Sampler(rvs=_rec([0.0, 0.0]), integrate_result=('R2', 'V2', 42.0, {})) + s.warm_rvs = _rec([-20.0, -20.0]) + _run(H, s) # rejected warm pass + assert s._av_state_reuse_safe is False + _run(H, s, neff=42.0) # healthy next event; returns before attempting a rescue + assert s._av_state_reuse_safe is True + + +def test_reject_message_reports_lnZ_on_the_events_offset_scale(capsys): + """lnL_offset is this event's manual_avoid_overflow_logarithm. + + It exists so the *** REJECTING *** line quotes absolute lnZ rather than the internally + offset value. Nothing else reads it, so dropping it at the call sites is invisible + unless a test drives it at a NON-ZERO value -- which is what made it possible to delete + `lnL_offset=manual_avoid_overflow_logarithm` from both call sites with 81 tests green. + """ + H = _load() + s = _Sampler(rvs=_rec([0.0, 0.0]), integrate_result=('R2', 'V2', 42.0, {'warm': True})) + s.warm_rvs = _rec([-20.0, -20.0]) + capsys.readouterr() + H['_maybe_l0_rescue'](s, 'R1', 'V1', 1.0, {'cold': True}, + lambda *a, **k: None, (), {}, lnL_offset=1000.0) + out = capsys.readouterr().out + assert "REJECTING" in out + # cold lnZ 0.0 and warm lnZ -20.0, both shifted by +1000 in the report + assert "1000.000" in out and "980.000" in out, \ + "the reject message did not quote lnZ on the event's offset scale: %r" % out + + +def test_both_call_sites_pass_the_events_offset(): + """Source-level, because the value comes from a local of each analyze_event.""" + src = _src() + assert src.count("lnL_offset=manual_avoid_overflow_logarithm") == 2, \ + "a call site dropped the event's lnL offset, so its reject message would quote " \ + "the internally-offset lnZ instead of the absolute one" + + +def test_accept_truncated_reports_the_warm_pass_anyway(): + H = _load(opts=_Opts(sampler_l0_rescue_accept_truncated=True)) + s = _Sampler(rvs=_rec([0.0, 0.0]), integrate_result=('R2', 'V2', 42.0, {'warm': True})) + s.warm_rvs = _rec([-20.0, -20.0]) + assert _run(H, s)[0] == 'R2' + + +def test_reject_threshold_is_respected(): + """A shortfall smaller than the threshold is not evidence of loss.""" + H = _load(opts=_Opts(sampler_l0_rescue_reject_dlnZ=50.0)) + s = _Sampler(rvs=_rec([0.0, 0.0]), integrate_result=('R2', 'V2', 42.0, {'warm': True})) + s.warm_rvs = _rec([-20.0, -20.0]) + assert _run(H, s)[0] == 'R2', "a 20-nat drop was rejected against a 50-nat threshold" + + +def test_a_raising_warm_pass_restores_the_cold_state(): + """The silent-for-a-campaign shape: _rvs holds warm samples, res/neff still hold cold.""" + H = _load() + cold = _rec([0.0, 0.0]) + s = _Sampler(rvs=cold, reserve={'tag': 'cold'}, raise_in_integrate=True) + s.warm_rvs = _rec([7.0, 7.0]) + out = _run(H, s, res='R1', var='V1', neff=1.0, dict_return={'cold': True}) + assert out == ('R1', 'V1', 1.0, {'cold': True}) + assert np.allclose(s._rvs['log_integrand'], cold['log_integrand']), \ + "cold diagnostics were reported beside a warm export" + assert s._warm_seed_reserve == {'tag': 'cold'} + assert s._av_state_reuse_safe is False, "the failed warm grid could be persisted" + + +def test_rescue_clears_warm_state_afterwards(): + H = _load() + s = _Sampler(rvs=_rec([0.0, 0.0]), integrate_result=('R2', 'V2', 42.0, {})) + s.warm_rvs = _rec([0.0, 0.0]) + _run(H, s) + assert s._warm is None, "the next point would draw from this point's contracted grid" + + +def test_mixed_lnZ_provenance_falls_back_to_a_like_for_like_comparison(): + """Cold read from the reserve, warm from the fair draw, is not a difference. + + The two readings differ by ~log(n_retained/eff_samp), so a mixed comparison manufactures + a gap of several nats out of nothing. The numbers here are chosen so the two paths + DISAGREE about the outcome -- an earlier version of this test used values where both + accepted, and it passed with the guard disabled. + + mixed (broken): cold 'retained' +10.0 vs warm 'fairdraw' 0.0 -> 10 nats -> REJECT + like-for-like : both re-read from _rvs, 0.0 vs 0.0 -> 0 nats -> ACCEPT + """ + class _AV(_FakeAV): + calls = {'n': 0} + + @classmethod + def lnZ_from_reserve(cls, reserve): + # available for the cold read, gone for the warm one + cls.calls['n'] += 1 + return 10.0 if cls.calls['n'] == 1 else None + _AV.calls['n'] = 0 + H = _load(av=_AV) + s = _Sampler(rvs=_rec([0.0, 0.0]), + reserve={'log_joint_prior': np.zeros(2), 'log_joint_s_prior': np.zeros(2)}, + integrate_result=('R2', 'V2', 42.0, {'warm': True})) + s.warm_rvs = _rec([0.0, 0.0]) + out = _run(H, s, res='R1', var='V1', neff=1.0, dict_return={'cold': True}) + assert out[0] == 'R2', ("a like-for-like lnZ comparison found no evidence of loss, so the " + "warm pass must stand; rejecting it means the gate compared a " + "'retained' reading against a 'fairdraw' one") + + +# ---------------------------------------------------------------------- source-level wiring +def _src(): + with open(_LISA) as fh: + return fh.read() + + +def test_both_analyze_event_variants_call_the_rescue(): + """This driver has two; a rescue wired into only one is a silent half-port.""" + tree = ast.parse(_src()) + fns = {n.name: n for n in tree.body + if isinstance(n, ast.FunctionDef) and n.name in ('analyze_event', 'analyze_event_LISA')} + assert set(fns) == {'analyze_event', 'analyze_event_LISA'} + for name, node in fns.items(): + called = any(isinstance(c, ast.Call) and isinstance(c.func, ast.Name) + and c.func.id == '_maybe_l0_rescue' for c in ast.walk(node)) + assert called, "%s does not call _maybe_l0_rescue" % name + + +def test_rescue_runs_before_the_no_result_guard(): + """Ordering is load-bearing. + + A degenerate early termination returns (None,None,None,None); `if not(res): raise` would + abort on it, skipping the strongest rescue trigger. In the main driver that guard sits + ~200 lines below the integrate call so the ordering is implicit -- here it is adjacent, + so it is pinned. + """ + src = _src() + guard = "if not(res): # no resut" + assert src.count(guard) == 2, "expected the guard in both analyze_event variants" + pos = 0 + for _ in range(2): + g = src.index(guard, pos) + call = src.rindex("_maybe_l0_rescue(", 0, g) + integ = src.rindex("sampler.integrate(like_to_integrate", 0, call) + assert integ < call < g, "the rescue must sit between integrate and the not(res) guard" + pos = g + 1 + + +def test_rescue_is_not_hidden_behind_the_LISA_flag(): + """Both variants get it; nothing keys the rescue off opts.LISA.""" + tree = ast.parse(_src()) + fn = [n for n in tree.body if isinstance(n, ast.FunctionDef) and n.name == '_maybe_l0_rescue'][0] + body = ast.dump(fn) + assert "'LISA'" not in body and 'attr=\'LISA\'' not in body + + +@pytest.mark.parametrize("opt,default", [ + ("--sampler-l0-rescue-reject-dlnZ", "default=3.0"), + ("--sampler-l0-rescue-puff-width-frac", "default=0.005"), + ("--sampler-l0-rescue-puff-factor", "default=2.0"), + ("--sampler-sequential-warmstart-deltalnL", "default=15.0"), +]) +def test_option_defaults_match_the_main_driver(opt, default): + """A knob that means something different in the two drivers is worse than a missing one. + + reject-dlnZ 3.0 in particular is a MEASURED value (L0_REJECT_DLNZ_MEASUREMENT.md); the + old 0.5 binned 25% of good portfolio warm passes while catching 0 of 55 truncated ones. + """ + for path in (_LISA, _MAIN): + with open(path) as fh: + src = fh.read() + i = src.index('"%s"' % opt) + line = src[i:src.index("\n", i)] + assert default.replace(" ", "") in line.replace(" ", ""), \ + "%s: %s does not carry %s" % (os.path.basename(path), opt, default) + + +# ------------------------------------------------------------------ anti-drift vs the main driver +def _normalized(fn): + node = ast.parse(ast.unparse(fn)).body[0] if hasattr(ast, "unparse") else fn + body = list(node.body) + if (body and isinstance(body[0], ast.Expr) + and isinstance(getattr(body[0], "value", None), ast.Constant) + and isinstance(body[0].value.value, str)): + body = body[1:] + return ast.dump(ast.fix_missing_locations(ast.Module(body=body, type_ignores=[]))) + + +@pytest.mark.parametrize("name", PORTED) +def test_ported_helper_is_identical_to_the_main_driver(name): + """Deliberate COPIES in a deliberate fork. Change one, change both.""" + assert _normalized(_defs(_LISA, [name])[name]) == _normalized(_defs(_MAIN, [name])[name]), \ + "%s has drifted between the two drivers (docstrings excluded)" % name diff --git a/MonteCarloMarginalizeCode/Code/test/test_lisa_portfolio_method_integrity.py b/MonteCarloMarginalizeCode/Code/test/test_lisa_portfolio_method_integrity.py new file mode 100644 index 000000000..f247ae2a3 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_lisa_portfolio_method_integrity.py @@ -0,0 +1,331 @@ +#!/usr/bin/env python +""" +`opts.sampler_method` must survive portfolio construction. + +THE DEFECT. Building a portfolio that carries a GMM member used to CLOBBER +`opts.sampler_method = 'GMM'`, so the GMM-specific argument blocks further down would run +and forward that member's config. It worked for that, and silently broke everything else +that asks "what sampler is this run using", because by then the honest answer -- 'portfolio' +-- had been overwritten. + +The consequence that matters here: the **L0 auto-rescue never fired for a portfolio**. Its +guard is + + opts.sampler_method in ('AV', 'portfolio') + +so for the single most common portfolio configuration -- one carrying a GMM member -- the +rescue silently declined, on a driver where the rescue had just been ported specifically +because LISA MBHB are high-SNR and that is the regime that stalls. No error, no log line; +the feature was simply absent. + +A portfolio also took GMM-only branches, `return_lnI` among them, which feeds +`rvs_integrand_is_lnL` and therefore how `ln_weights_from_rvs` reads the record. + +THE FIX, ported from the main driver, which had already made it: flag the member +non-destructively. `opts.sampler_method` stays 'portfolio'; the GMM blocks key off +`use_gmm_args = (sampler_method == "GMM") or use_gmm_member`. + +WHY THIS FILE IS SHAPED AS AN INVARIANT. A mutation of a shared option is not a FUNC, +OPTION, CONST or ATTR, so the drift audit produces zero gap items for it -- the same blind +spot that hid the missing AV/`use_lnL` branch. The first test below is therefore the +general rule ("nothing assigns opts.sampler_method") rather than a check on this one site, +because the next such clobber will be somewhere else. +""" + +import ast +import os + +import pytest + +_HERE = os.path.dirname(os.path.abspath(__file__)) +_LISA = os.path.join(_HERE, '..', 'bin', 'integrate_likelihood_extrinsic_batchmode_lisa') +_MAIN = os.path.join(_HERE, '..', 'bin', 'integrate_likelihood_extrinsic_batchmode') + + +def _src(path): + with open(path) as fh: + return fh.read() + + +def _assignments_to(path, attr): + """Line numbers where `opts.` is assigned (=, augmented, or walrus-ish).""" + tree = ast.parse(_src(path), filename=path) + hits = [] + for node in ast.walk(tree): + targets = [] + if isinstance(node, ast.Assign): + targets = node.targets + elif isinstance(node, ast.AugAssign): + targets = [node.target] + for t in targets: + if (isinstance(t, ast.Attribute) and t.attr == attr + and isinstance(t.value, ast.Name) and t.value.id == 'opts'): + hits.append(node.lineno) + return hits + + +# ------------------------------------------------------------------------- the invariant +@pytest.mark.parametrize("path,label", [(_LISA, 'lisa'), (_MAIN, 'main')]) +def test_nothing_assigns_opts_sampler_method(path, label): + """The general rule, in BOTH drivers. + + `opts.sampler_method` is read by the L0 rescue gate, the AV state save, the use_lnL + branch table and several per-event resets. Any code that reassigns it makes every one + of those answer a question about a sampler the run is not using. + """ + hits = _assignments_to(path, 'sampler_method') + assert not hits, ( + "%s driver assigns opts.sampler_method at line(s) %s. Flag the condition " + "non-destructively (see use_gmm_member) instead of overwriting the run's identity." + % (label, hits)) + + +def test_the_rescue_guard_still_reads_sampler_method(): + """If the guard stops reading it, the invariant above protects nothing. + + Pins the two together so neither can be quietly relaxed on its own. + """ + assert "opts.sampler_method in ('AV', 'portfolio')" in _src(_LISA) + + +# ------------------------------------------------------------------- the replacement flag +def test_portfolio_loop_flags_a_GMM_member_without_clobbering(): + src = _src(_LISA) + assert 'use_gmm_member = True' in src, "the GMM member is not flagged at all" + assert "opts.sampler_method = 'GMM'" not in src, "the clobber is back" + assert 'use_gmm_member=False' in src, "the flag is never initialised" + + +def test_use_gmm_args_is_standalone_GMM_or_a_portfolio_member(): + assert 'use_gmm_args = (opts.sampler_method == "GMM") or use_gmm_member' in _src(_LISA) + + +def test_use_gmm_args_is_defined_before_every_use(): + src = _src(_LISA) + define = src.index('use_gmm_args = (opts.sampler_method') + first_use = src.index('if use_gmm_args:') + assert define < first_use + tree = ast.parse(src) + define_line = min(n.lineno for n in ast.walk(tree) + if isinstance(n, ast.Assign) and len(n.targets) == 1 + and getattr(n.targets[0], 'id', None) == 'use_gmm_args') + module_level_uses = [n.lineno for n in ast.walk(tree) + if isinstance(n, ast.Name) and n.id == 'use_gmm_args' + and isinstance(n.ctx, ast.Load)] + # uses inside analyze_event run later regardless; only module-level order can break. + assert min(module_level_uses) >= define_line + + +def test_the_GMM_setup_block_runs_for_a_portfolio_member(): + """This is what the clobber existed to achieve, now achieved honestly.""" + src = _src(_LISA) + i = src.index('use_gmm_args = (opts.sampler_method') + block = src[i:i + 400] + assert 'if use_gmm_args:' in block, "the GMM setup block no longer runs for a portfolio member" + + +def test_per_event_gmm_resets_key_off_use_gmm_args(): + """gmm_dict exists for a portfolio-with-GMM too, so the resets must reach it. + + Two analyze_event variants plus the --force-reset-all block: three sites. + """ + src = _src(_LISA) + assert src.count('elif use_gmm_args:') == 3, \ + "expected the two per-event resets and --force-reset-all to key off use_gmm_args" + + +def test_return_lnI_still_keys_on_the_method_not_the_member(): + """A portfolio must NOT take the GMM lnL branch. + + This is the other half of the clobber's damage: with sampler_method overwritten, a + portfolio run set return_lnI, which flips rvs_integrand_is_lnL and changes how + ln_weights_from_rvs reads the record. + """ + src = _src(_LISA) + assert 'if opts.sampler_method=="GMM" and opts.internal_use_lnL:' in src + i = src.index('if opts.sampler_method=="GMM" and opts.internal_use_lnL:') + assert 'return_lnI' in src[i:i + 300] + # and it must not have been widened to the member flag + assert 'use_gmm_args' not in src[i:i + 300], \ + "the return_lnI branch was widened to portfolios carrying a GMM member" + + +# ------------------------------------------------------------- the rescue actually fires +def _load_rescue(sampler_method): + """Exec the rescue with a chosen opts.sampler_method.""" + names = ['_rvs_lnL_convention', 'ln_weights_from_rvs', '_rvs_len', + '_rvs_is_export_resample', '_rvs_is_equal_weight', 'ln_weights_for_posterior', + '_lnZ_of_rvs', '_kish_neff_of_rvs', '_lnZ_of_reserve_or_rvs', + '_snapshot_pass_state', '_restore_pass_state', '_warm_seed_reserve_for', + '_warm_seed_geometry', '_clear_warm_state', '_maybe_l0_rescue'] + import numpy as np + defs = {n.name: n for n in ast.parse(_src(_LISA)).body + if isinstance(n, ast.FunctionDef) and n.name in names} + mod = ast.Module(body=[defs[n] for n in names], type_ignores=[]) + + class _AV(object): + @staticmethod + def lnZ_from_reserve(r): + return None + + @staticmethod + def build_warm_seed(cols, lnL, lo, hi, axes, **kw): + return np.asarray(cols, dtype=float), {'puffed': False, 'n_core': 3, + 'rank_core': 3, 'dim': 3, + 'rank_final': 3, 'n_puff': 0, + 'puff_scale': 'auto'} + + opts = type('O', (), { + 'sampler_method': sampler_method, 'sampler_warmstart_retry_neff': 5.0, + 'sampler_l0_rescue_reject_dlnZ': 3.0, 'sampler_l0_rescue_accept_truncated': False, + 'sampler_l0_rescue_puff_scale': 'auto', 'sampler_l0_rescue_puff_width_frac': 0.005, + 'sampler_l0_rescue_puff_factor': 2.0, + 'sampler_sequential_warmstart_deltalnL': 15.0})() + ns = {"numpy": np, "np": np, "opts": opts, "mcsamplerAdaptiveVolume": _AV} + exec(compile(ast.fix_missing_locations(mod), "rescue", "exec"), ns) + return ns + + +class _Sampler(object): + def __init__(self): + import numpy as np + n = 3 + self._rvs = {'log_integrand': np.zeros(n), 'log_joint_prior': np.zeros(n), + 'log_joint_s_prior': np.zeros(n), + 'a': np.linspace(0.1, 0.9, n), 'b': np.linspace(0.2, 0.8, n)} + self._warm_seed_reserve = None + self.params_ordered = ['a', 'b'] + self.llim = {'a': 0.0, 'b': 0.0} + self.rlim = {'a': 1.0, 'b': 1.0} + self.portfolio_realizations = [] + self._warm = None + self._warm_applied = False + self.bootstrapped = None + + def identity_convert(self, x): + return x + + def bootstrap_from_samples(self, seed, cover_frac=0.0): + self.bootstrapped = seed + + def integrate(self, fn, *a, **k): + return ('R2', 'V2', 42.0, {'warm': True}) + + +@pytest.mark.parametrize("method", ['portfolio', 'AV']) +def test_rescue_fires_for_both_eligible_methods(method): + """The end the whole fix serves. + + With the clobber, a portfolio carrying a GMM member arrived here as 'GMM' and this + returned untouched -- no bootstrap, no warm pass, no message. + """ + ns = _load_rescue(method) + s = _Sampler() + out = ns['_maybe_l0_rescue'](s, 'R1', 'V1', 1.0, {'cold': True}, + lambda *a, **k: None, (), {}) + assert s.bootstrapped is not None, "the rescue did not fire for %s" % method + assert out[2] == 42.0 + + +def test_rescue_declines_for_a_clobbered_method(): + """The failure mode itself, pinned: if the method ever reads 'GMM', the rescue is off. + + Not an argument that declining for standalone GMM is wrong -- it is correct, GMM has no + bootstrap_from_samples in practice. It documents that the guard is exactly what the + clobber defeated, so the invariant above is what protects it. + """ + ns = _load_rescue('GMM') + s = _Sampler() + ns['_maybe_l0_rescue'](s, 'R1', 'V1', 1.0, {'cold': True}, lambda *a, **k: None, (), {}) + assert s.bootstrapped is None + + +# --------------------------------------------------- the member-dispatch chain itself +def _member_loop(path): + """The `for name in sampler_types:` loop body, as AST.""" + for node in ast.walk(ast.parse(_src(path), filename=path)): + if (isinstance(node, ast.For) and isinstance(node.target, ast.Name) + and node.target.id == 'name' + and isinstance(node.iter, ast.Name) and node.iter.id == 'sampler_types'): + return node + raise AssertionError("no `for name in sampler_types` loop in %s" % os.path.basename(path)) + + +@pytest.mark.parametrize("path,label", [(_LISA, 'lisa'), (_MAIN, 'main')]) +def test_member_dispatch_is_a_single_elif_chain(path, label): + """A chain of separate `if`s reuses the previous member on an unmatched name. + + With `if name == 'AV': ... ; if name == 'GMM': ...` a name matching NOTHING falls + through every test and leaves `sampler` bound to whatever it last held -- the plain + MCSampler built before the chain, or on later iterations the PREVIOUS member -- which is + then appended. A typo in --sampler-portfolio silently produced a DUPLICATE member + rather than an error. + """ + loop = _member_loop(path) + # Only the statements that DISPATCH ON THE MEMBER NAME. The loop body also holds an + # `if hasattr(sampler, 'xpy')` after the chain in both drivers, which is not part of it. + dispatch = [st for st in loop.body + if isinstance(st, ast.If) + and any(isinstance(n, ast.Name) and n.id == 'name' + for n in ast.walk(st.test))] + assert len(dispatch) == 1, ( + "%s: member dispatch is %d separate `if` statements, not one elif chain; an " + "unmatched name reuses the previous member" % (label, len(dispatch))) + + +@pytest.mark.parametrize("path,label", [(_LISA, 'lisa'), (_MAIN, 'main')]) +def test_an_unknown_member_name_raises(path, label): + """The chain must END in an else that raises, not fall off silently.""" + node = [st for st in _member_loop(path).body + if isinstance(st, ast.If) + and any(isinstance(n, ast.Name) and n.id == 'name' + for n in ast.walk(st.test))][0] + while isinstance(node, ast.If): + tail = node.orelse + if len(tail) == 1 and isinstance(tail[0], ast.If): + node = tail[0] + continue + break + assert tail, "%s: the member dispatch chain has no else clause" % label + assert any(isinstance(st, ast.Raise) for st in tail), ( + "%s: the else clause does not raise, so an unknown --sampler-portfolio member is " + "accepted silently" % label) + + +def test_plugin_pipelines_are_dispatched_before_the_error(): + """A plugin member (nflow, ...) must CONSTRUCT, not fall through to the raise. + + Checking that the string "known_pipelines" merely appears is not enough: it also appears + in the error message, so deleting the whole dispatch branch left that check green. Walk + the chain and require a branch that both TESTS and SUBSCRIPTS known_pipelines. + """ + node = [st for st in _member_loop(_LISA).body + if isinstance(st, ast.If) + and any(isinstance(n, ast.Name) and n.id == 'name' + for n in ast.walk(st.test))][0] + found = False + while isinstance(node, ast.If): + tests_it = any(isinstance(a, ast.Attribute) and a.attr == 'known_pipelines' + for a in ast.walk(node.test)) + builds_it = any(isinstance(sub, ast.Subscript) + and any(isinstance(a, ast.Attribute) and a.attr == 'known_pipelines' + for a in ast.walk(sub.value)) + for sub in ast.walk(ast.Module(body=node.body, type_ignores=[]))) + if tests_it and builds_it: + found = True + break + node = node.orelse[0] if (len(node.orelse) == 1 + and isinstance(node.orelse[0], ast.If)) else None + if node is None: + break + assert found, ("no branch dispatches to mcsamplerPortfolio.known_pipelines, so a plugin " + "member falls through to the unknown-member error") + + +def test_the_unknown_member_error_names_what_is_known(): + assert "--sampler-portfolio: unknown member" in _src(_LISA) + + +def test_AC_is_accepted_as_an_alias(): + """main accepts 'AC' alongside 'adaptive_cartesian_gpu'; a portfolio spec is shared.""" + assert "name == 'AC'" in _src(_LISA) diff --git a/MonteCarloMarginalizeCode/Code/test/test_lisa_sampler_plumbing.py b/MonteCarloMarginalizeCode/Code/test/test_lisa_sampler_plumbing.py new file mode 100644 index 000000000..976d17169 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_lisa_sampler_plumbing.py @@ -0,0 +1,223 @@ +#!/usr/bin/env python +""" +Tests for the portfolio freeze/allocation policy ported into the LISA ILE driver +(bin/integrate_likelihood_extrinsic_batchmode_lisa). + +This is pure PASS-THROUGH plumbing to samplers the LISA driver already wires -- it exposes +the same ``ok_lnL_methods`` as the main driver (``GMM, adaptive_cartesian, +adaptive_cartesian_gpu, AV, portfolio``, verified identical) and builds mcsamplerPortfolio +the same way. Before this port the knobs were reachable only through +``--sampler-portfolio-args``, an eval-able dict; the pipeline passes the named flags. + +WHAT CAN ACTUALLY GO WRONG HERE, and is therefore what these tests check: + + * a default that differs between the two drivers. Worse than a missing option: the same + command line then means two different things depending on which driver ran it. + * an option that is UNSET leaking into the kwargs as ``None`` and overriding the sampler's + own default with nothing. The assembly's whole shape -- ``if opts.x is not None`` -- + exists for that, and a single dropped guard is invisible until a run behaves oddly. + * the two mutually-exclusive VARAHA flags resolving the wrong way round. + +The freeze-policy assembly is inline in both drivers (not a function), so it is exercised +here by extracting the block and exec'ing it against a fake ``opts``. That tests the real +source, not a paraphrase of it. +""" + +import ast +import os +import re +import textwrap + +import pytest + +_HERE = os.path.dirname(os.path.abspath(__file__)) +_LISA = os.path.join(_HERE, '..', 'bin', 'integrate_likelihood_extrinsic_batchmode_lisa') +_MAIN = os.path.join(_HERE, '..', 'bin', 'integrate_likelihood_extrinsic_batchmode') + +PORTFOLIO_OPTS = [ + "--portfolio-adaptive-alloc", "--portfolio-alloc-exponent", "--portfolio-freeze-wt", + "--portfolio-grace-iters", "--portfolio-probe-period", "--portfolio-quality-signal", + "--portfolio-revive-period", "--portfolio-varaha-can-freeze", + "--portfolio-varaha-max-frac", "--portfolio-varaha-min-frac", + "--portfolio-varaha-never-freeze", "--portfolio-weight-clip", +] + + +def _src(path): + with open(path) as fh: + return fh.read() + + +def _option_nodes(path): + """{'--foo': ast.Call} for every add_option in a driver.""" + out = {} + for n in ast.walk(ast.parse(_src(path), filename=path)): + if (isinstance(n, ast.Call) and isinstance(n.func, ast.Attribute) + and n.func.attr in ("add_option", "add_argument")): + names = [a.value for a in n.args + if isinstance(a, ast.Constant) and isinstance(a.value, str)] + if names and names[0].startswith("--"): + out[names[0]] = n + return out + + +def _kwargs_of(node): + out = {} + for kw in node.keywords: + try: + out[kw.arg] = ast.literal_eval(kw.value) + except Exception: + out[kw.arg] = ast.dump(kw.value) + return out + + +@pytest.fixture(scope="module") +def opts_lisa(): + return _option_nodes(_LISA) + + +@pytest.fixture(scope="module") +def opts_main(): + return _option_nodes(_MAIN) + + +# ------------------------------------------------------------------------------ presence +@pytest.mark.parametrize("opt", PORTFOLIO_OPTS) +def test_option_is_present_in_the_lisa_driver(opt, opts_lisa): + assert opt in opts_lisa + + +# ------------------------------------------------------------------------------- defaults +@pytest.mark.parametrize("opt", PORTFOLIO_OPTS) +def test_option_signature_matches_the_main_driver(opt, opts_lisa, opts_main): + """Same default, same type, same action. + + A knob that means something different in the two drivers is worse than a missing one: + the same pipeline command line would then produce two different integrations. + """ + a, b = _kwargs_of(opts_lisa[opt]), _kwargs_of(opts_main[opt]) + for key in ("default", "type", "action", "choices"): + assert a.get(key) == b.get(key), ( + "%s: %s differs (lisa=%r, main=%r)" % (opt, key, a.get(key), b.get(key))) + + +@pytest.mark.parametrize("opt", [ + "--portfolio-alloc-exponent", "--portfolio-freeze-wt", "--portfolio-grace-iters", + "--portfolio-probe-period", "--portfolio-quality-signal", "--portfolio-revive-period", + "--portfolio-varaha-max-frac", "--portfolio-varaha-min-frac", "--portfolio-weight-clip", +]) +def test_tuning_options_default_to_none_so_the_sampler_keeps_its_own(opt, opts_lisa): + """None is the sentinel the assembly keys on. A default of 0/0.0 would silently + override the sampler's built-in value for every run that never set the flag.""" + assert _kwargs_of(opts_lisa[opt]).get("default") is None + + +@pytest.mark.parametrize("opt", [ + "--portfolio-adaptive-alloc", "--portfolio-varaha-can-freeze", + "--portfolio-varaha-never-freeze", +]) +def test_flags_are_store_true_and_default_false(opt, opts_lisa): + kw = _kwargs_of(opts_lisa[opt]) + assert kw.get("action") == "store_true" and kw.get("default") is False + + +# --------------------------------------------------------- the assembly block, executed +_START = "_freeze_policy_kwargs = {}" +_END = 'print(" PORTFOLIO freeze-policy overrides: "' + + +def _assembly_block(path): + """The inline freeze-policy assembly, dedented so it can be exec'd on its own. + + Slice from the START OF THE LINE holding the sentinel, not from the sentinel itself: + otherwise the first line carries no indentation while the rest do, and dedent finds no + common prefix. + """ + src = _src(path) + i = src.rindex("\n", 0, src.index(_START)) + 1 + j = src.index(_END, i) + j = src.rindex("\n", i, j) + 1 + return textwrap.dedent(src[i:j]) + + +class _Opts(object): + """Every portfolio option at its documented default.""" + portfolio_grace_iters = None + portfolio_revive_period = None + portfolio_freeze_wt = None + portfolio_varaha_can_freeze = False + portfolio_varaha_never_freeze = False + portfolio_adaptive_alloc = False + portfolio_varaha_min_frac = None + portfolio_varaha_max_frac = None + portfolio_weight_clip = None + portfolio_quality_signal = None + portfolio_alloc_exponent = None + portfolio_probe_period = None + + def __init__(self, **kw): + for k, v in kw.items(): + assert hasattr(type(self), k), "unknown option %s" % k + setattr(self, k, v) + + +def _assemble(**kw): + ns = {"opts": _Opts(**kw)} + exec(compile(_assembly_block(_LISA), "freeze_policy", "exec"), ns) + return ns["_freeze_policy_kwargs"] + + +def test_nothing_set_means_nothing_overridden(): + """The important one: an all-defaults run must not touch the sampler's policy at all.""" + assert _assemble() == {} + + +def test_each_tuning_option_passes_through_when_set(): + got = _assemble(portfolio_grace_iters=7, portfolio_revive_period=3, + portfolio_freeze_wt=0.25, portfolio_varaha_min_frac=0.2, + portfolio_varaha_max_frac=0.8, portfolio_weight_clip=1.0, + portfolio_quality_signal='credit', portfolio_alloc_exponent=2.0, + portfolio_probe_period=5) + assert got == {'portfolio_grace_iters': 7, 'portfolio_revive_period': 3, + 'portfolio_freeze_wt': 0.25, 'portfolio_varaha_min_frac': 0.2, + 'portfolio_varaha_max_frac': 0.8, 'portfolio_weight_clip': 1.0, + 'portfolio_quality_signal': 'credit', 'portfolio_alloc_exponent': 2.0, + 'portfolio_probe_period': 5} + + +def test_zero_is_passed_through_not_treated_as_unset(): + """0 disables probing/reviving and is a REAL value; `if x:` would drop it.""" + got = _assemble(portfolio_probe_period=0, portfolio_revive_period=0) + assert got == {'portfolio_probe_period': 0, 'portfolio_revive_period': 0} + + +def test_varaha_never_freeze_sets_true(): + assert _assemble(portfolio_varaha_never_freeze=True) == {'portfolio_varaha_never_freeze': True} + + +def test_varaha_can_freeze_sets_false(): + assert _assemble(portfolio_varaha_can_freeze=True) == {'portfolio_varaha_never_freeze': False} + + +def test_can_freeze_wins_when_both_are_given(): + """Documented precedence; the two flags are mutually exclusive.""" + got = _assemble(portfolio_varaha_can_freeze=True, portfolio_varaha_never_freeze=True) + assert got == {'portfolio_varaha_never_freeze': False} + + +def test_adaptive_alloc_is_opt_in_only(): + assert 'portfolio_adaptive_alloc' not in _assemble() + assert _assemble(portfolio_adaptive_alloc=True) == {'portfolio_adaptive_alloc': True} + + +def test_assembly_block_is_identical_to_the_main_drivers(): + """Deliberate copies in a deliberate fork. Change one, change both.""" + def norm(s): + return re.sub(r"\s+", " ", s).strip() + assert norm(_assembly_block(_LISA)) == norm(_assembly_block(_MAIN)) + + +def test_assembly_result_is_actually_handed_to_setup(): + """Building the dict and not passing it would be a silent no-op.""" + src = _src(_LISA) + assert "sampler.setup(portfolio_args=opts.sampler_portfolio_args, **_freeze_policy_kwargs" in src diff --git a/MonteCarloMarginalizeCode/Code/test/test_lisa_use_lnL_branches.py b/MonteCarloMarginalizeCode/Code/test/test_lisa_use_lnL_branches.py new file mode 100644 index 000000000..894881634 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_lisa_use_lnL_branches.py @@ -0,0 +1,193 @@ +#!/usr/bin/env python +""" +The per-sampler `use_lnL` / `return_lnI` branches in the LISA ILE driver. + +FOUND BY ADVERSARIAL AUDIT, NOT BY THE DRIFT GATE. The main driver has + + if opts.sampler_method == "AV" and opts.internal_use_lnL: + return_lnL = True + pinned_params.update({"use_lnL": True}) + +with the comment: *"without this, --internal-use-lnL --sampler-method AV passed the +ok_lnL_methods check but silently did nothing, so exp(lnL) overflowed at high SNR when no +logarithm offset was set."* The LISA driver had branches for GMM, adaptive_cartesian_gpu +and portfolio -- and none for AV. + +High SNR is the LISA MBHB regime, so this is the case, not an edge. + +WHY THE DRIFT AUDIT MISSED IT, and why this file exists. A missing `if` branch is not a +FUNC, OPTION, CONST or ATTR, so it produces zero gap items. The audit is a name-presence +set difference; behaviour behind a shared name is invisible to it. These tests close that +specific hole by pinning the branch TABLE in both drivers against each other. +""" + +import ast +import os + +import pytest + +_HERE = os.path.dirname(os.path.abspath(__file__)) +_LISA = os.path.join(_HERE, '..', 'bin', 'integrate_likelihood_extrinsic_batchmode_lisa') +_MAIN = os.path.join(_HERE, '..', 'bin', 'integrate_likelihood_extrinsic_batchmode') + +# Samplers both drivers accept. Verified identical in both ok_lnL_methods lists. +METHODS = ['GMM', 'adaptive_cartesian', 'adaptive_cartesian_gpu', 'AV', 'portfolio'] + + +def _src(path): + with open(path) as fh: + return fh.read() + + +def _pinned_updates(path): + """{method: {key: value}} for every `pinned_params.update({...})` guarded by a method test. + + Walks module-level `if` statements, works out which sampler method each one is about + from the string constants in its test, and records the pinned_params keys it sets. + """ + tree = ast.parse(_src(path), filename=path) + out = {} + for node in tree.body: + if not isinstance(node, ast.If): + continue + methods = {c.value for c in ast.walk(node.test) + if isinstance(c, ast.Constant) and c.value in METHODS} + if not methods: + continue + uses_lnL_opt = any(isinstance(a, ast.Attribute) and a.attr == 'internal_use_lnL' + for a in ast.walk(node.test)) + keys = {} + for call in ast.walk(node): + if (isinstance(call, ast.Call) and isinstance(call.func, ast.Attribute) + and call.func.attr == 'update' + and isinstance(call.func.value, ast.Name) + and call.func.value.id == 'pinned_params'): + for arg in call.args: + if isinstance(arg, ast.Dict): + for k, v in zip(arg.keys, arg.values): + if isinstance(k, ast.Constant): + try: + keys[k.value] = ast.literal_eval(v) + except Exception: + keys[k.value] = '' + if keys: + for m in methods: + rec = out.setdefault(m, {"keys": {}, "gated_on_internal_use_lnL": False}) + rec["keys"].update(keys) + rec["gated_on_internal_use_lnL"] |= uses_lnL_opt + return out + + +@pytest.fixture(scope="module") +def lisa(): + return _pinned_updates(_LISA) + + +@pytest.fixture(scope="module") +def main(): + return _pinned_updates(_MAIN) + + +def test_AV_sets_use_lnL_under_internal_use_lnL(lisa): + """The regression this file exists for. + + Without it, --sampler-method AV --internal-use-lnL is a SILENT no-op: the option passes + the ok_lnL_methods check and changes nothing, so the integrand stays linear and exp(lnL) + overflows at high SNR unless a manual logarithm offset happens to be set. + """ + assert 'AV' in lisa, "no AV branch sets pinned_params at all" + assert lisa['AV']['keys'].get('use_lnL') is True, \ + "--sampler-method AV --internal-use-lnL does not set use_lnL: silent no-op" + assert lisa['AV']['gated_on_internal_use_lnL'], \ + "the AV branch must be gated on --internal-use-lnL, not unconditional" + + +@pytest.mark.parametrize("method", ['GMM', 'adaptive_cartesian_gpu', 'AV']) +def test_branch_table_matches_the_main_driver(method, lisa, main): + """Same method -> same pinned_params keys in both drivers. + + This is the check that would have caught the missing AV branch, and it is the shape the + name-based drift audit cannot express. + """ + assert method in main, "the main driver has no %s branch to compare against" % method + assert method in lisa, "the LISA driver has no %s branch" % method + assert lisa[method]['keys'] == main[method]['keys'], ( + "%s: pinned_params differ (lisa=%r, main=%r)" + % (method, lisa[method]['keys'], main[method]['keys'])) + + +def test_portfolio_differs_from_main_only_by_the_deferred_GMM_forwarding(lisa, main): + """portfolio is the ONE branch still divergent, and only in a known, recorded way. + + The main driver's portfolio branch also forwards the --internal-gmm-* knobs to its GMM + member (gmm_adaptive / gmm_defensive_frac / gmm_inflate). Those options are deliberately + deferred: main wires them through its group-pairing setup, and this driver's GMM block is + structured differently, so they need their own pass. + + Asserting the delta EXACTLY -- rather than skipping portfolio -- means any OTHER + divergence in this branch still fails, and this test tightens on its own once the GMM + pass lands. + """ + deferred = {'gmm_adaptive', 'gmm_defensive_frac', 'gmm_inflate'} + lk, mk = lisa['portfolio']['keys'], main['portfolio']['keys'] + assert set(mk) - set(lk) == deferred, ( + "portfolio branch diverges beyond the deferred GMM forwarding: missing here = %s" + % sorted(set(mk) - set(lk))) + assert not set(lk) - set(mk), "the LISA portfolio branch sets keys main does not: %s" \ + % sorted(set(lk) - set(mk)) + for k in set(lk) & set(mk): + assert lk[k] == mk[k], "portfolio: %s differs (lisa=%r, main=%r)" % (k, lk[k], mk[k]) + + +def test_only_GMM_requests_return_lnI(lisa): + """return_lnI is what makes 'integrand' hold lnL, and it drives rvs_integrand_is_lnL. + + If another sampler gains it, the stored-convention derivation has to be revisited -- + ln_weights_from_rvs reads that convention to decide whether to log the integrand. + """ + with_lnI = {m for m, rec in lisa.items() if rec['keys'].get('return_lnI') is True} + assert with_lnI == {'GMM'}, "unexpected return_lnI set: %s" % sorted(with_lnI) + + +def test_adaptive_cartesian_has_no_use_lnL_branch(lisa): + """Plain adaptive_cartesian (RIFT.integrators.mcsampler) has no use_lnL handling at all. + + It always stores linear L. A branch here would make ln_weights_from_rvs read its + records as lnL, which is the failure the helper's docstring warns about. + """ + assert 'adaptive_cartesian' not in lisa or \ + 'use_lnL' not in lisa['adaptive_cartesian']['keys'] + + +def test_the_convention_is_still_derived_from_pinned_params(): + """Adding a branch must not tempt anyone back to the CLI option.""" + src = _src(_LISA) + assert 'rvs_integrand_is_lnL = bool(pinned_params.get("return_lnI", False))' in src + + +def test_the_convention_is_derived_after_every_branch_that_could_set_return_lnI(): + """Ordering: pinned_params must be final where the convention is read off it. + + The main driver derives it "where pinned_params is final". If a later update ever + carried return_lnI, deriving it early would silently pick the wrong convention. + """ + src = _src(_LISA) + tree = ast.parse(src) + derive_line = None + for node in ast.walk(tree): + if (isinstance(node, ast.Assign) and len(node.targets) == 1 + and getattr(node.targets[0], 'id', None) == 'rvs_integrand_is_lnL'): + derive_line = node.lineno + assert derive_line is not None, "rvs_integrand_is_lnL is never assigned" + # AST, not a text search: 'return_lnI' also appears in docstrings that DESCRIBE the + # convention, and an earlier version of this test matched those and failed on prose. + later = [c.lineno for c in ast.walk(tree) + if isinstance(c, ast.Call) and isinstance(c.func, ast.Attribute) + and c.func.attr == 'update' + and isinstance(c.func.value, ast.Name) and c.func.value.id == 'pinned_params' + and c.lineno > derive_line + and any(isinstance(k, ast.Constant) and k.value == 'return_lnI' + for a in c.args if isinstance(a, ast.Dict) for k in a.keys)] + assert not later, \ + "pinned_params gains return_lnI at line(s) %s, AFTER the stored convention is " \ + "derived from it at line %d" % (later, derive_line)