From 0b74e00eceb77fa2e7451a7a7e173e8aaa1b8062 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Sat, 30 May 2026 20:10:50 -0400 Subject: [PATCH 01/18] calmarg DAG: fix write_calpilot_sub module + add local-condor DAG test target MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Bug fix: create_event_parameter_pipeline_BasicIteration imports dag_utils_generic (aliased dag_utils), so --calmarg-pilot crashed with AttributeError -- write_calpilot_sub was only added to dag_utils.py. Add it to dag_utils_generic.py too (the module the builder actually uses). Caught by the new local DAG build test below. demo/rift/calmarg/Makefile: new `make dag-build` / `dag-validate` / `dag-run` targets that build a real RIFT DAG with --calmarg-pilot (3-IFO + GPU + AV, mirroring the known-good ILE-GPU-Paper batch_gpu target) and validate the cal-pilot topology without needing condor. Verified the produced DAG: - CALPILOT jobs for it=0 (macroiterationprev=-1) and it=1 (=0), running util_CalPilotStage.py - PARENT unify_0 CHILD CALPILOT_0 (pilot after iteration 0 composite, ∥ CIP_0) - PARENT CALPILOT_0 CHILD (the seed barrier: wide_{N+1} waits on pilot_N) - ILE.sub seeded via --calibration-proposal-breadcrumb cal_consolidated_$(macroiterationprev).npz `make dag-run` submits it to local condor (GPU); select the card with CUDA_VISIBLE_DEVICES. Co-Authored-By: Claude Opus 4.8 --- .../Code/RIFT/misc/dag_utils_generic.py | 54 ++++++++++++ .../Code/demo/rift/calmarg/Makefile | 84 ++++++++++++++++++- 2 files changed, 137 insertions(+), 1 deletion(-) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/misc/dag_utils_generic.py b/MonteCarloMarginalizeCode/Code/RIFT/misc/dag_utils_generic.py index 2f320b16..fcb2b41e 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/misc/dag_utils_generic.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/misc/dag_utils_generic.py @@ -2742,6 +2742,60 @@ def write_consolidate_sub_simple(tag='consolidate', exe=None, base=None,target=N +def write_calpilot_sub(tag='calpilot', exe=None, log_dir=None, universe="vanilla", + working_directory=None, ile_args_file=None, top_fraction=0.05, + max_points=32, request_memory=4096, request_gpu=True, + singularity_image=None, max_runtime_minutes=300, **kwargs): + """Submit file for a calibration PILOT stage (Option C; see + RIFT/calmarg/DESIGN_adaptive_driver.md): harvest top-lnL points from iteration + $(macroiteration)'s composite, run ILE --calibration-dump-responsibilities on them, + fit + consolidate a cal proposal that seeds wide_{N+1}. Runs util_CalPilotStage.py. + + Macros expected at instantiation: macroiteration, macroiterationprev. + Produces /cal_consolidated_$(macroiteration).npz (consumed by wide_{N+1} ILE via + --calibration-proposal-breadcrumb). + """ + exe = exe or which("util_CalPilotStage.py") + wd = working_directory + job = pipeline.CondorDAGJob(universe=universe, executable=exe) + sub_name = tag + '.sub' + job.set_sub_file(sub_name) + + # arguments (per-iteration files via condor macros) + job.add_opt("composite", wd + "/consolidated_$(macroiteration).composite") + job.add_opt("ile-args-file", ile_args_file) + job.add_opt("iteration", "$(macroiteration)") + job.add_opt("output-breadcrumb", wd + "/cal_consolidated_$(macroiteration).npz") + job.add_opt("prev-breadcrumb", wd + "/cal_consolidated_$(macroiterationprev).npz") + job.add_opt("top-fraction", str(top_fraction)) + job.add_opt("max-points", str(max_points)) + job.add_opt("workdir", wd) + + uniq_str = "$(macroiteration)-$(cluster)-$(process)" + job.set_log_file("%s%s-%s.log" % (log_dir, tag, uniq_str)) + job.set_stderr_file("%s%s-%s.err" % (log_dir, tag, uniq_str)) + job.set_stdout_file("%s%s-%s.out" % (log_dir, tag, uniq_str)) + + if default_resolved_env: + job.add_condor_cmd('environment', default_resolved_env) + else: + job.add_condor_cmd('getenv', default_getenv_value) + job.add_condor_cmd('request_memory', str(request_memory) + "M") + if request_gpu: + job.add_condor_cmd('request_GPUs', '1') # the pilot runs ILE (GPU path) + if singularity_image: + job.add_condor_cmd("+SingularityImage", '"' + singularity_image + '"') + try: + job.add_condor_cmd('accounting_group', os.environ['LIGO_ACCOUNTING']) + job.add_condor_cmd('accounting_group_user', os.environ['LIGO_USER_NAME']) + except: + print(" LIGO accounting information not available. Add manually to %s !" % sub_name) + if not (max_runtime_minutes is None): + remove_str = 'JobStatus =?= 2 && (CurrentTime - JobStartDate) > ( {})'.format(60 * max_runtime_minutes) + job.add_condor_cmd('periodic_remove', remove_str) + return job, sub_name + + def write_unify_sub_simple(tag='unify', exe=None, base=None,target=None,universe="vanilla",arg_str=None,log_dir=None, use_eos=False,ncopies=1,no_grid=False, max_runtime_minutes=60,extra_text='',**kwargs): """ Write a submit file for launching a consolidation job diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile b/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile index c2faee73..e6336a6d 100644 --- a/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile +++ b/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile @@ -169,6 +169,88 @@ $(LOWSNR_CACHE): low-snr: lowsnr-inputs $(MAKE) all CACHE=$(LOWSNR_CACHE) SIM_XML=$(CURDIR)/mdc_lowsnr.xml.gz +# --- Full DAG test (Option C adaptive cal pilot) on local condor ---------------- +# Builds and (optionally) submits a real RIFT DAG that exercises the per-iteration +# calibration PILOT wiring: calpilot_N runs in parallel with CIP_N, learns a cal +# proposal, and SEEDS the wide ILE jobs of iteration N+1 via a breadcrumb. Mirrors +# the known-good .travis/ILE-GPU-Paper test_workflow_batch_gpu_lowlatency target +# (calls create_event_parameter_pipeline_BasicIteration directly), 3-IFO + GPU + AV, +# with --calmarg-pilot added. Designed to run on a single machine with condor + GPU +# (e.g. cardassia); select the card with CUDA_VISIBLE_DEVICES if needed. +RUN_DAG := $(CURDIR)/rundir_dag +N_IT_DAG ?= 2 # >=2 so iteration 1 consumes the seed from calpilot_0 +NPTS_GRID_DAG ?= 80 # initial intrinsic grid size +NPTS_IT_DAG ?= 12 # intrinsic points / ILE jobs per iteration +NCAL_DAG ?= 20 +NMAX_DAG ?= 40000 +NEFF_DAG ?= 15 +NCHUNK_DAG ?= 10000 + +DAG_ILE_ARGS := --n-chunk $(NCHUNK_DAG) --time-marginalization --sim-xml overlap-grid.xml.gz \ + --reference-freq 100.0 --adapt-weight-exponent 0.1 --event-time $(EVENT_TIME) --save-P 0.1 \ + --cache-file $(CACHE) --fmin-template 10 --n-max $(NMAX_DAG) --fmax 1700.0 --save-deltalnL inf \ + --l-max 2 --n-eff $(NEFF_DAG) --approximant SEOBNRv4 --adapt-floor-level 0.1 --d-max $(DMAX) \ + --psd-file H1=$(PSD) --psd-file L1=$(PSD) --psd-file V1=$(PSD) \ + --channel-name H1=FAKE-STRAIN --channel-name L1=FAKE-STRAIN --channel-name V1=FAKE-STRAIN \ + --inclination-cosine-sampler --declination-cosine-sampler --data-start-time 1000000008 \ + --data-end-time 1000000016 --inv-spec-trunc-time 0 --no-adapt-after-first --no-adapt-distance \ + --srate 4096 --sampler-method AV --vectorized --gpu \ + --calibration-envelope-directory $(CURDIR)/cal_env --calibration-n-realizations $(NCAL_DAG) \ + --calibration-spline-count 10 +DAG_CIP_ARGS := --mc-range [23,35] --eta-range [0.20,0.24999] --parameter mc --parameter-implied eta --parameter-nofit delta_mc --fit-method gp --lnL-offset 120 --cap-points 12000 --n-output-samples 5000 --no-plots --n-eff 5000 +DAG_TEST_ARGS := --always-succeed --method lame --parameter m1 +DAG_PLOT_ARGS := --parameter m1 --parameter m2 +DAG_FILE := marginalize_intrinsic_parameters_BasicIterationWorkflow.dag + +.PHONY: dag dag-build dag-validate dag-run + +# Build the DAG (no condor submission) and validate the cal-pilot wiring is present. +dag-build: cal_env/H1.txt + rm -rf $(RUN_DAG); mkdir -p $(RUN_DAG) + cd $(RUN_DAG) && $(ENV) util_ManualOverlapGrid.py --parameter mc --parameter-range '[23,35]' \ + --parameter delta_mc --parameter-range '[0.0,0.5]' --grid-cartesian-npts $(NPTS_GRID_DAG) \ + --skip-overlap --approx SEOBNRv4 + cd $(RUN_DAG) && cp overlap-grid.xml.gz overlap-grid-0.xml.gz + cd $(RUN_DAG) && printf 'X %s\n' "$(DAG_CIP_ARGS)" > args_cip.txt + cd $(RUN_DAG) && printf 'X %s\n' "$(DAG_TEST_ARGS)" > args_test.txt + cd $(RUN_DAG) && printf 'X %s\n' "$(DAG_PLOT_ARGS)" > args_plot.txt + # args_ile.txt: the wide ILE config + the per-iteration seed breadcrumb. The condor + # macro $(macroiterationprev) is emitted literally ($$ -> $ for Make; single-quoted + # for the shell), exactly as util_RIFT_pseudo_pipe.py would inject it. + cd $(RUN_DAG) && { printf 'X %s' "$(DAG_ILE_ARGS)"; \ + printf ' --calibration-proposal-breadcrumb %s/cal_consolidated_$$(macroiterationprev).npz\n' "$(RUN_DAG)"; } > args_ile.txt + cd $(RUN_DAG) && $(ENV) create_event_parameter_pipeline_BasicIteration --request-gpu-ILE \ + --ile-n-events-to-analyze $(NPTS_IT_DAG) --input-grid overlap-grid-0.xml.gz \ + --ile-exe $(ILE) \ + --ile-args $(RUN_DAG)/args_ile.txt --cip-args $(RUN_DAG)/args_cip.txt \ + --test-args $(RUN_DAG)/args_test.txt --plot-args $(RUN_DAG)/args_plot.txt \ + --request-memory-CIP 2048 --request-memory-ILE 4096 \ + --input-grid $(RUN_DAG)/overlap-grid.xml.gz --n-samples-per-job $(NPTS_IT_DAG) \ + --working-directory $(RUN_DAG) --n-iterations $(N_IT_DAG) \ + --calmarg-pilot --calmarg-pilot-cadence 1 --calmarg-pilot-max-it 1 + $(MAKE) dag-validate + +# Validate the produced DAG carries the cal-pilot topology (no condor needed). +dag-validate: + @test -s "$(RUN_DAG)/CALPILOT.sub" || (echo "FAIL: no CALPILOT.sub" && false) + @grep -q "util_CalPilotStage.py" "$(RUN_DAG)/CALPILOT.sub" || (echo "FAIL: CALPILOT.sub does not run util_CalPilotStage.py" && false) + @grep -q "CALPILOT" $(RUN_DAG)/$(DAG_FILE) || (echo "FAIL: no CALPILOT job in DAG" && false) + @grep -q -- "--calibration-proposal-breadcrumb" "$(RUN_DAG)/ILE.sub" || (echo "FAIL: wide ILE.sub missing seed breadcrumb" && false) + @grep -q "macroiterationprev" "$(RUN_DAG)/ILE.sub" || (echo "FAIL: ILE.sub missing macroiterationprev" && false) + @grep -q "cal_consolidated" "$(RUN_DAG)/CALPILOT.sub" || (echo "FAIL: CALPILOT.sub missing consolidated output" && false) + @echo "OK: DAG built with cal-pilot wiring -- CALPILOT job present, runs util_CalPilotStage.py," + @echo " wide ILE seeded via --calibration-proposal-breadcrumb cal_consolidated_\$$(macroiterationprev).npz." + @echo " Submit it with: make dag-run (or: cd $(RUN_DAG) && condor_submit_dag $(DAG_FILE))" + +# Submit the DAG to the local condor (GPU). Select the card with CUDA_VISIBLE_DEVICES. +dag-run: + cd $(RUN_DAG) && $(ENV) condor_submit_dag -f $(DAG_FILE) + @echo "Submitted. Watch with: condor_q ; tail -f $(RUN_DAG)/$(DAG_FILE).dagman.out" + @echo "After it finishes, the pilot products are $(RUN_DAG)/cal_consolidated_*.npz and" + @echo "$(RUN_DAG)/cal_pilot_resp_*.npz ; the seeded wide-ILE logs report 'SEEDED from proposal breadcrumb'." + +dag: dag-build + clean: rm -f distance_marg.npz out_*.dat *.xml.gz lowsnr.cache snr_*.dat snr-report.txt - rm -rf cal_env lowsnr_mdc + rm -rf cal_env lowsnr_mdc rundir_dag From b39001183ba3e2be1f612d6c936db76cdcde40b5 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Sat, 30 May 2026 21:10:25 -0400 Subject: [PATCH 02/18] Fix stdout pollution corrupting composites (precession warning -> stderr) + robust harvest Local condor smoke test of the cal-pilot DAG failed at the CIP_0 and CALPILOT_0 nodes, both with parse errors on the iteration-0 composite ("could not convert 'Import' to float" / "number of columns changed from 13 to 10"). Root cause: lalsimutils.py printed the missing-`precession` warning to STDOUT, and since many RIFT tools build .dat/.composite files from stdout, that line corrupted the data (one ragged 10-col row among 13-col rows). Not a calmarg bug, but it breaks any pipeline run on a precession-less environment. - RIFT/lalsimutils.py: print the precession ImportError warning to sys.stderr, not stdout. - bin/util_CalHarvestGrid.py: read defensively (genfromtxt usecols=indx..lnL, invalid_raise=False, drop NaN rows) so a ragged/dirty composite degrades gracefully instead of failing the pilot. Verified on the actual corrupted composite (6/27 pts). Caught by `make dag-run` in demo/rift/calmarg (local condor + GPU on cardassia). Co-Authored-By: Claude Opus 4.8 --- MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py | 5 ++++- .../Code/bin/util_CalHarvestGrid.py | 10 +++++++++- 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py index c48fc1d4..c86870b7 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py @@ -51,7 +51,10 @@ try: import precession except ImportError: - print('Import Error - module missing. Please install the module "precession."') + # MUST go to stderr: many RIFT tools write data to stdout (composite/.dat files are + # built from it), and a stray stdout line here corrupts those files -> downstream + # parsers (CIP, util_CalHarvestGrid) choke on ragged/non-numeric rows. + print('Import Error - module missing. Please install the module "precession."', file=sys.stderr) def safe_int(mystr): try: return int(mystr) diff --git a/MonteCarloMarginalizeCode/Code/bin/util_CalHarvestGrid.py b/MonteCarloMarginalizeCode/Code/bin/util_CalHarvestGrid.py index 761cc0c2..aa28cffe 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_CalHarvestGrid.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_CalHarvestGrid.py @@ -67,7 +67,15 @@ def main(argv=None): help="column index of m1 (m2..spins follow)") opts = p.parse_args(argv) - dat = np.loadtxt(opts.composite) + # Robust read: composites are sometimes ragged (variable trailing columns, or a + # stray non-numeric line). Read only the columns we need (indx..lnL) and skip rows + # that don't parse, rather than failing the whole pilot. + n_need = max(opts.lnL_col, opts.mass_col + 7) + 1 + dat = np.genfromtxt(opts.composite, usecols=range(n_need), invalid_raise=False) + dat = np.atleast_2d(dat) + dat = dat[~np.isnan(dat).any(axis=1)] # drop any rows that failed to parse + if len(dat) == 0: + raise ValueError("no valid numeric rows in composite %s" % opts.composite) rows = harvest_top_rows(dat, lnL_col=opts.lnL_col, top_fraction=opts.top_fraction, max_points=opts.max_points) P_list = rows_to_P_list(rows, mass_col=opts.mass_col) From a5d893975f99b5fdd3d7a7d583ed6a10965bcd68 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Sun, 31 May 2026 04:14:44 -0400 Subject: [PATCH 03/18] calmarg pilot: relax first CIP stage sigma-cut for the prior-cal cold start Pipeline blocker found by the local cal-pilot DAG run: iteration 0 draws cal realizations from the broad PRIOR (no pilot proposal learned yet), so its per-point Monte-Carlo error is large (~0.7-0.9 here). CIP's default --sigma-cut 0.6 then strips EVERY point ("Stripped size (0,)"), and the DAG dies. This is generic to any calmarg run, not the demo. Long-term fix (as suggested -- at the helper level, using the iteration structure): - helper_LDG_Events.py: new --calmarg-first-cip-sigma-cut; when set, relax ONLY the FIRST CIP stage's --sigma-cut (that stage runs the cold-start/prior-cal iterations). Later stages run on pilot-seeded iterations with normal errors and keep the default cut. - util_RIFT_pseudo_pipe.py: --calmarg-first-cip-sigma-cut (default 100); passed to the helper automatically when --calmarg-pilot is enabled. Short-term (demo): demo/rift/calmarg Makefile DAG_CIP_ARGS gets --sigma-cut 5 so the local smoke test flows past CIP_0 into the seeded iteration 1. Co-Authored-By: Claude Opus 4.8 --- MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py | 9 +++++++++ .../Code/bin/util_RIFT_pseudo_pipe.py | 5 +++++ .../Code/demo/rift/calmarg/Makefile | 6 +++++- 3 files changed, 19 insertions(+), 1 deletion(-) diff --git a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py index 900d43bd..21b71b69 100755 --- a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py +++ b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py @@ -210,6 +210,7 @@ def get_observing_run(t): parser.add_argument("--assume-volumetric-spin",action='store_true',help="If present, the code will assume a volumetric spin prior in its last iterations. If *not* present, the code will adopt a uniform magnitude spin prior in its last iterations. If not present, generally more iterations are taken.") parser.add_argument("--assume-highq",action='store_true',help="If present, the code will adopt a strategy that drops spin2. Also the precessing strategy will allow perpendicular spin to play a role early on (rather than as a subdominant parameter later)") parser.add_argument("--assume-well-placed",action='store_true',help="If present, the code will adopt a strategy that assumes the initial grid is very well placed, and will minimize the number of early iterations performed. Not as extrme as --propose-flat-strategy") +parser.add_argument("--calmarg-first-cip-sigma-cut",default=None,type=float,help="In-loop calibration marginalization: relax the CIP --sigma-cut on the FIRST cip stage to this value. The cold-start iterations draw cal realizations from the broad PRIOR (no pilot proposal learned yet), so their per-point Monte-Carlo error is large and the default --sigma-cut 0.6 would strip every point. util_RIFT_pseudo_pipe.py passes this automatically when --calmarg-pilot is enabled.") parser.add_argument("--propose-ile-convergence-options",action='store_true',help="If present, the code will try to adjust the adaptation options, Nmax, etc based on experience") parser.add_argument("--internal-propose-ile-convergence-freezeadapt",action='store_true',help="If present, uses the --no-adapt-after-first --no-adapt-distance options (at one point default)") parser.add_argument("--internal-propose-ile-adapt-log",action='store_true',help="If present, uses the --adapt-log argument. Useful for very loud signals. Note only lnL information is used for adapting, not prior, so samples will be *uniform* in prior range if lnL is low") @@ -1741,6 +1742,14 @@ def lambda_m_estimate(m): helper_cip_arg_list += [helper_cip_last_it] +# In-loop calmarg cold start: the first CIP stage runs on iterations whose cal draws come +# from the broad PRIOR (the pilot has not yet learned/seeded a proposal), so their MC error +# is large. Relax that stage's --sigma-cut so the cold-start points are not all stripped +# (CIP default 0.6). Only the FIRST stage is relaxed; later stages run on pilot-seeded +# iterations with normal errors and keep the default cut. +if opts.calmarg_first_cip_sigma_cut is not None and len(helper_cip_arg_list) > 0: + helper_cip_arg_list[0] += " --sigma-cut {} ".format(opts.calmarg_first_cip_sigma_cut) + with open("helper_cip_arg_list.txt",'w+') as f: f.write("\n".join(helper_cip_arg_list)) diff --git a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py index 5f8f1eba..e8d0800f 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py @@ -165,6 +165,7 @@ def unsafe_parse_arg_string_dict(my_argstr): parser.add_argument("--calmarg-pilot-max-it",default=3,type=int,help="Stop launching cal pilots after this iteration (cal is boring; freeze once learned). Default 3.") parser.add_argument("--calmarg-pilot-top-fraction",default=0.05,type=float,help="Fraction of highest-lnL composite points the pilot harvests. Default 0.05.") parser.add_argument("--calmarg-pilot-max-points",default=32,type=int,help="Cap on harvested pilot points per iteration. Default 32.") +parser.add_argument("--calmarg-first-cip-sigma-cut",default=100.0,type=float,help="With --calmarg-pilot: relax the first CIP stage's --sigma-cut to this value, so cold-start (prior-cal) iteration-0 points -- which have large MC error -- are not all stripped by CIP's default 0.6. Threaded to helper_LDG_Events.py. Default 100 (effectively keep all cold-start points).") parser.add_argument("--distance-reweighting",action='store_true',help="Option to add job to DAG to reweight posterior samples due to different distance prior (LVK prod prior)") parser.add_argument("--extra-args-helper",action=None, help="Filename with arguments for the helper. Use to provide alternative channel names and other advanced configuration (--channel-name, data type)!") parser.add_argument("--manual-postfix",default='',type=str) @@ -822,6 +823,10 @@ def unsafe_parse_arg_string_dict(my_argstr): cmd += " --data-LI-seglen "+str(opts.data_LI_seglen) if opts.assume_well_placed: cmd += " --assume-well-placed " +if opts.calmarg_pilot: + # cold-start cal pilots draw cal from the broad PRIOR -> large MC error on iteration 0; + # relax the first CIP stage's sigma-cut so those points are not all stripped. + cmd += " --calmarg-first-cip-sigma-cut {} ".format(opts.calmarg_first_cip_sigma_cut) #if is_event_bns and not opts.no_matter: # cmd += " --assume-matter " # npts_it = 1000 diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile b/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile index e6336a6d..43f7743b 100644 --- a/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile +++ b/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile @@ -197,7 +197,11 @@ DAG_ILE_ARGS := --n-chunk $(NCHUNK_DAG) --time-marginalization --sim-xml overlap --srate 4096 --sampler-method AV --vectorized --gpu \ --calibration-envelope-directory $(CURDIR)/cal_env --calibration-n-realizations $(NCAL_DAG) \ --calibration-spline-count 10 -DAG_CIP_ARGS := --mc-range [23,35] --eta-range [0.20,0.24999] --parameter mc --parameter-implied eta --parameter-nofit delta_mc --fit-method gp --lnL-offset 120 --cap-points 12000 --n-output-samples 5000 --no-plots --n-eff 5000 +# --sigma-cut 5: iteration 0 uses PRIOR cal draws (no pilot seed yet) -> low neff_cal -> +# high per-point MC error (~0.7-0.9 here), which the CIP default --sigma-cut 0.6 would +# strip entirely. Relax it so the cold-start points survive (cf. the long-term +# helper_LDG_Events.py fix that relaxes the first CIP stage when cal pilots are enabled). +DAG_CIP_ARGS := --mc-range [23,35] --eta-range [0.20,0.24999] --parameter mc --parameter-implied eta --parameter-nofit delta_mc --fit-method gp --lnL-offset 120 --sigma-cut 5 --cap-points 12000 --n-output-samples 5000 --no-plots --n-eff 5000 DAG_TEST_ARGS := --always-succeed --method lame --parameter m1 DAG_PLOT_ARGS := --parameter m1 --parameter m2 DAG_FILE := marginalize_intrinsic_parameters_BasicIterationWorkflow.dag From 5ece22dd594660a1a6532a4aac31d9c1d90e2f31 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Sun, 31 May 2026 09:10:18 -0400 Subject: [PATCH 04/18] demo/calmarg: parameterize DAG target (PILOT/FUSED) + vanilla fused-calmarg path The cal-pilot DAG run exposed pilot starvation (a 60-dim cal proposal fit from only ~20 realizations -> near-degenerate covariance -> pathological seeded iteration-1 lnL). Per review priority, add the cleaner production path first: VANILLA in-loop calmarg, no pilots. Makefile dag-build/dag-validate now take PILOT={0,1} and FUSED={0,1} (parse-time ifeq): - PILOT=0 omits --calmarg-pilot and the proposal-breadcrumb seed (every iteration uses prior cal draws -- correct for the no-pilot production path). - FUSED=1 adds --calibration-fused-kernel (Option C fused GPU path). - dag-validate branches: PILOT=1 asserts CALPILOT+breadcrumb; PILOT=0 asserts envelope + (optional) fused flag and NO CALPILOT. Priority demo: make dag-build PILOT=0 FUSED=1 && make dag-run PILOT=0 FUSED=1 (verified the build: ILE.sub carries --calibration-fused-kernel + envelope, no breadcrumb, no CALPILOT.) Pilot starvation backstop (larger pilot n_cal + prior-shrinkage covariance) is tracked separately. Co-Authored-By: Claude Opus 4.8 --- .../Code/demo/rift/calmarg/Makefile | 62 +++++++++++++------ 1 file changed, 44 insertions(+), 18 deletions(-) diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile b/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile index 43f7743b..6f4f3d36 100644 --- a/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile +++ b/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile @@ -178,7 +178,13 @@ low-snr: lowsnr-inputs # with --calmarg-pilot added. Designed to run on a single machine with condor + GPU # (e.g. cardassia); select the card with CUDA_VISIBLE_DEVICES if needed. RUN_DAG := $(CURDIR)/rundir_dag -N_IT_DAG ?= 2 # >=2 so iteration 1 consumes the seed from calpilot_0 +# Toggles: PILOT=1 adds the adaptive cal pilots (seed wide_{N+1}); PILOT=0 is VANILLA +# in-loop calmarg (no pilots -- the clean production path, best for faint sources). +# FUSED=1 uses the fused GPU kernel (Option C). E.g. the priority vanilla-fused demo: +# make dag-build PILOT=0 FUSED=1 && make dag-run PILOT=0 FUSED=1 +PILOT ?= 1 +FUSED ?= 0 +N_IT_DAG ?= 2 # >=2 so iteration 1 consumes the seed from calpilot_0 (PILOT=1) NPTS_GRID_DAG ?= 80 # initial intrinsic grid size NPTS_IT_DAG ?= 12 # intrinsic points / ILE jobs per iteration NCAL_DAG ?= 20 @@ -186,6 +192,17 @@ NMAX_DAG ?= 40000 NEFF_DAG ?= 15 NCHUNK_DAG ?= 10000 +ifeq ($(FUSED),1) + DAG_FUSED_FLAG := --calibration-fused-kernel +else + DAG_FUSED_FLAG := +endif +ifeq ($(PILOT),1) + DAG_PILOT_CEPP := --calmarg-pilot --calmarg-pilot-cadence 1 --calmarg-pilot-max-it 1 +else + DAG_PILOT_CEPP := +endif + DAG_ILE_ARGS := --n-chunk $(NCHUNK_DAG) --time-marginalization --sim-xml overlap-grid.xml.gz \ --reference-freq 100.0 --adapt-weight-exponent 0.1 --event-time $(EVENT_TIME) --save-P 0.1 \ --cache-file $(CACHE) --fmin-template 10 --n-max $(NMAX_DAG) --fmax 1700.0 --save-deltalnL inf \ @@ -196,11 +213,11 @@ DAG_ILE_ARGS := --n-chunk $(NCHUNK_DAG) --time-marginalization --sim-xml overlap --data-end-time 1000000016 --inv-spec-trunc-time 0 --no-adapt-after-first --no-adapt-distance \ --srate 4096 --sampler-method AV --vectorized --gpu \ --calibration-envelope-directory $(CURDIR)/cal_env --calibration-n-realizations $(NCAL_DAG) \ - --calibration-spline-count 10 -# --sigma-cut 5: iteration 0 uses PRIOR cal draws (no pilot seed yet) -> low neff_cal -> -# high per-point MC error (~0.7-0.9 here), which the CIP default --sigma-cut 0.6 would -# strip entirely. Relax it so the cold-start points survive (cf. the long-term -# helper_LDG_Events.py fix that relaxes the first CIP stage when cal pilots are enabled). + --calibration-spline-count 10 $(DAG_FUSED_FLAG) +# --sigma-cut 5: cal draws (PRIOR for vanilla every iteration; PRIOR at iteration 0 for the +# pilot cold start) have high per-point MC error (~0.7-0.9 here), which the CIP default +# --sigma-cut 0.6 would strip entirely. Relax it (cf. the helper_LDG_Events.py fix that +# relaxes the first CIP stage automatically when --calmarg-pilot is enabled). DAG_CIP_ARGS := --mc-range [23,35] --eta-range [0.20,0.24999] --parameter mc --parameter-implied eta --parameter-nofit delta_mc --fit-method gp --lnL-offset 120 --sigma-cut 5 --cap-points 12000 --n-output-samples 5000 --no-plots --n-eff 5000 DAG_TEST_ARGS := --always-succeed --method lame --parameter m1 DAG_PLOT_ARGS := --parameter m1 --parameter m2 @@ -208,7 +225,7 @@ DAG_FILE := marginalize_intrinsic_parameters_BasicIterationWorkflow.dag .PHONY: dag dag-build dag-validate dag-run -# Build the DAG (no condor submission) and validate the cal-pilot wiring is present. +# Build the DAG (no condor submission) and validate the wiring. PILOT/FUSED toggles above. dag-build: cal_env/H1.txt rm -rf $(RUN_DAG); mkdir -p $(RUN_DAG) cd $(RUN_DAG) && $(ENV) util_ManualOverlapGrid.py --parameter mc --parameter-range '[23,35]' \ @@ -218,11 +235,12 @@ dag-build: cal_env/H1.txt cd $(RUN_DAG) && printf 'X %s\n' "$(DAG_CIP_ARGS)" > args_cip.txt cd $(RUN_DAG) && printf 'X %s\n' "$(DAG_TEST_ARGS)" > args_test.txt cd $(RUN_DAG) && printf 'X %s\n' "$(DAG_PLOT_ARGS)" > args_plot.txt - # args_ile.txt: the wide ILE config + the per-iteration seed breadcrumb. The condor - # macro $(macroiterationprev) is emitted literally ($$ -> $ for Make; single-quoted - # for the shell), exactly as util_RIFT_pseudo_pipe.py would inject it. + # args_ile.txt: the wide ILE config; for PILOT=1 also the per-iteration seed breadcrumb + # (the condor macro $(macroiterationprev) is emitted literally: $$ -> $ for Make, then + # single-quoted for the shell -- exactly as util_RIFT_pseudo_pipe.py injects it). cd $(RUN_DAG) && { printf 'X %s' "$(DAG_ILE_ARGS)"; \ - printf ' --calibration-proposal-breadcrumb %s/cal_consolidated_$$(macroiterationprev).npz\n' "$(RUN_DAG)"; } > args_ile.txt + [ "$(PILOT)" = "1" ] && printf ' --calibration-proposal-breadcrumb %s/cal_consolidated_$$(macroiterationprev).npz' "$(RUN_DAG)"; \ + printf '\n'; } > args_ile.txt cd $(RUN_DAG) && $(ENV) create_event_parameter_pipeline_BasicIteration --request-gpu-ILE \ --ile-n-events-to-analyze $(NPTS_IT_DAG) --input-grid overlap-grid-0.xml.gz \ --ile-exe $(ILE) \ @@ -231,20 +249,28 @@ dag-build: cal_env/H1.txt --request-memory-CIP 2048 --request-memory-ILE 4096 \ --input-grid $(RUN_DAG)/overlap-grid.xml.gz --n-samples-per-job $(NPTS_IT_DAG) \ --working-directory $(RUN_DAG) --n-iterations $(N_IT_DAG) \ - --calmarg-pilot --calmarg-pilot-cadence 1 --calmarg-pilot-max-it 1 - $(MAKE) dag-validate + $(DAG_PILOT_CEPP) + $(MAKE) dag-validate PILOT=$(PILOT) FUSED=$(FUSED) -# Validate the produced DAG carries the cal-pilot topology (no condor needed). +# Validate the produced DAG (no condor needed). Branch on PILOT. dag-validate: + @grep -q -- "--calibration-envelope-directory" "$(RUN_DAG)/ILE.sub" || (echo "FAIL: wide ILE.sub missing calibration envelope" && false) +ifeq ($(FUSED),1) + @grep -q -- "--calibration-fused-kernel" "$(RUN_DAG)/ILE.sub" || (echo "FAIL: FUSED=1 but ILE.sub missing --calibration-fused-kernel" && false) +endif +ifeq ($(PILOT),1) @test -s "$(RUN_DAG)/CALPILOT.sub" || (echo "FAIL: no CALPILOT.sub" && false) @grep -q "util_CalPilotStage.py" "$(RUN_DAG)/CALPILOT.sub" || (echo "FAIL: CALPILOT.sub does not run util_CalPilotStage.py" && false) @grep -q "CALPILOT" $(RUN_DAG)/$(DAG_FILE) || (echo "FAIL: no CALPILOT job in DAG" && false) @grep -q -- "--calibration-proposal-breadcrumb" "$(RUN_DAG)/ILE.sub" || (echo "FAIL: wide ILE.sub missing seed breadcrumb" && false) @grep -q "macroiterationprev" "$(RUN_DAG)/ILE.sub" || (echo "FAIL: ILE.sub missing macroiterationprev" && false) - @grep -q "cal_consolidated" "$(RUN_DAG)/CALPILOT.sub" || (echo "FAIL: CALPILOT.sub missing consolidated output" && false) - @echo "OK: DAG built with cal-pilot wiring -- CALPILOT job present, runs util_CalPilotStage.py," - @echo " wide ILE seeded via --calibration-proposal-breadcrumb cal_consolidated_\$$(macroiterationprev).npz." - @echo " Submit it with: make dag-run (or: cd $(RUN_DAG) && condor_submit_dag $(DAG_FILE))" + @echo "OK: cal-PILOT DAG built -- CALPILOT runs util_CalPilotStage.py; wide ILE seeded via breadcrumb." +else + @! test -s "$(RUN_DAG)/CALPILOT.sub" || (echo "FAIL: PILOT=0 but a CALPILOT.sub was produced" && false) + @! grep -q -- "--calibration-proposal-breadcrumb" "$(RUN_DAG)/ILE.sub" || (echo "FAIL: PILOT=0 but ILE.sub has a seed breadcrumb" && false) + @echo "OK: VANILLA calmarg DAG built (no pilots; FUSED=$(FUSED)) -- ILE.sub carries the cal envelope, no CALPILOT." +endif + @echo " Submit with: make dag-run PILOT=$(PILOT) FUSED=$(FUSED)" # Submit the DAG to the local condor (GPU). Select the card with CUDA_VISIBLE_DEVICES. dag-run: From ec2226d73ec9458c7ba7ac0b207ae21815613255 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Sun, 31 May 2026 11:59:30 -0400 Subject: [PATCH 05/18] calmarg pilot backstop: shrink fitted proposal toward prior (fixes starved-cov collapse) Corroborated the seeded-iteration pathology: the pilot fit a 60-D cal Gaussian from only ~20 realizations, so the covariance was rank-deficient and its ~40 uninformed directions collapsed to the cov_floor (1e-8) -- a near-delta proposal -> seeded log(prior/proposal) blew up -> iteration-1 lnL collapsed (+109 -> -131). Fix (robustness, the real backstop): adaptive.fit_proposal gains prior_sigma + shrink. When prior_sigma is given it shrinks the fitted covariance toward diag(prior_sigma**2) with weight rho = (dim+1)/(dim+1+neff): ~all-prior when starved (neff << dim), ->all-data when neff >> dim. Uninformed directions now keep ~prior width (log_w ~ 0) instead of collapsing. Verified on a 60-D/neff=1 starved fit: min cov eigenvalue 1e-8 -> 0.98 (prior var). Threaded prior_sigma through util_CalPilotFit and adaptive.adaptive_cal. Existing tests unchanged: pilot brute-vs-seeded |dlogZ|=0.01 x254; adaptive_cal neff 140/300. Also (demo, per review): N_COPIES_DAG default 1 (drop builder's default-2 redundancy, ~2x turnaround on a single-card test); threaded --n-copies into the demo DAG build. Co-Authored-By: Claude Opus 4.8 --- .../Code/RIFT/calmarg/adaptive.py | 30 +++++++++++++++++-- .../Code/bin/util_CalPilotFit.py | 6 +++- .../Code/demo/rift/calmarg/Makefile | 4 ++- 3 files changed, 35 insertions(+), 5 deletions(-) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/calmarg/adaptive.py b/MonteCarloMarginalizeCode/Code/RIFT/calmarg/adaptive.py index feebf8fb..3e4904df 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/calmarg/adaptive.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/calmarg/adaptive.py @@ -94,10 +94,24 @@ def nodes_to_cal_factors(amp_nodes, phase_nodes, log_f_nodes, T_segment, dT, fmi # --------------------------------------------------------------------------- # Tempered proposal fit + diagnostics # --------------------------------------------------------------------------- -def fit_proposal(nodes, log_resp, beta, cov_floor=1e-8, cov_inflate=1.0): +def fit_proposal(nodes, log_resp, beta, cov_floor=1e-8, cov_inflate=1.0, + prior_sigma=None, shrink=None): """Tempered weighted-Gaussian fit. Weights = softmax(beta * log_resp). beta in (0,1]: small -> broad (many samples), 1 -> full responsibility weighting. + + prior_sigma : if given (length-dim 1-sigma of the diagonal prior), SHRINK the fitted + covariance toward diag(prior_sigma**2). This is essential when the fit is starved + -- a weighted sample covariance from ~neff effective points cannot constrain the + dim*(dim+1)/2 entries of a dim-dimensional covariance (cal node space is ~60-D), + so the UNINFORMED directions otherwise collapse to ~0 variance. A near-zero + proposal variance is a near-delta: seeded draws are pinned and the importance + weights log(prior/proposal) blow up, producing the pathological seeded likelihoods + we saw. Shrinking keeps uninformed directions at ~prior width (log_w ~ 0 there). + shrink : explicit shrinkage weight rho in [0,1] toward the prior; default auto = + (dim+1)/(dim+1+neff), i.e. ~1 (all prior) when starved, ->0 (all data) when + neff >> dim. + Returns (mean, cov).""" lw = beta * log_resp lw = lw - logsumexp(lw) @@ -105,7 +119,16 @@ def fit_proposal(nodes, log_resp, beta, cov_floor=1e-8, cov_inflate=1.0): mean = w @ nodes d = nodes - mean cov = (w[:, None] * d).T @ d - cov = cov_inflate * cov + cov_floor * np.eye(mean.shape[0]) + dim = mean.shape[0] + if prior_sigma is not None: + prior_sigma = np.asarray(prior_sigma, dtype=float) + neff = neff_from_logweights(beta * log_resp) + rho = shrink if shrink is not None else (dim + 1.0) / (dim + 1.0 + neff) + rho = float(min(max(rho, 0.0), 1.0)) + cov = (1.0 - rho) * cov_inflate * cov + rho * np.diag(prior_sigma ** 2) + else: + cov = cov_inflate * cov + cov = cov + cov_floor * np.eye(dim) return mean, cov @@ -156,7 +179,8 @@ def adaptive_cal(evaluate, prior_mean, prior_sigma, n_nodes_amp, # next proposal from tempered posterior responsibilities; inflate the covariance # early (while tempering is on) to keep exploring, relax as beta -> 1. beta = float(betas[min(it, len(betas) - 1)]) - mean, cov = fit_proposal(nodes, log_resp, beta, cov_inflate=1.0 + (1.0 - beta)) + mean, cov = fit_proposal(nodes, log_resp, beta, cov_inflate=1.0 + (1.0 - beta), + prior_sigma=prior_sigma) neff_resp = neff_from_logweights(log_resp) neff_w = neff_from_logweights(log_w) history.append(dict(iter=it, beta=beta, neff_resp=neff_resp, neff_w=neff_w)) diff --git a/MonteCarloMarginalizeCode/Code/bin/util_CalPilotFit.py b/MonteCarloMarginalizeCode/Code/bin/util_CalPilotFit.py index 56e1fa3f..2afc0de9 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_CalPilotFit.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_CalPilotFit.py @@ -84,7 +84,11 @@ def main(argv=None): print(" Auto-tempering: beta=%.3f (target neff>=%.0f; tempered neff=%.1f)" % (beta, target, adaptive.neff_from_logweights(beta * log_resp))) - mean, cov = adaptive.fit_proposal(nodes, log_resp, beta, cov_inflate=opts.cov_inflate) + # shrink toward the prior diagonal: in ~60-D cal node space a pilot with neff~tens + # cannot constrain the full covariance, so uninformed directions must default to + # ~prior width (else the proposal collapses to a near-delta and seeded weights blow up) + mean, cov = adaptive.fit_proposal(nodes, log_resp, beta, cov_inflate=opts.cov_inflate, + prior_sigma=dumps[0]["prior_sigma"]) cal = dict(proposal_mean=mean, proposal_cov=cov, prior_mean=dumps[0]["prior_mean"], prior_sigma=dumps[0]["prior_sigma"], diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile b/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile index 6f4f3d36..04da9628 100644 --- a/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile +++ b/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile @@ -191,6 +191,8 @@ NCAL_DAG ?= 20 NMAX_DAG ?= 40000 NEFF_DAG ?= 15 NCHUNK_DAG ?= 10000 +N_COPIES_DAG ?= 1 # ILE job duplicates (builder default 2); 1 -> no redundancy, + # ~2x faster turnaround on a single-card test ifeq ($(FUSED),1) DAG_FUSED_FLAG := --calibration-fused-kernel @@ -246,7 +248,7 @@ dag-build: cal_env/H1.txt --ile-exe $(ILE) \ --ile-args $(RUN_DAG)/args_ile.txt --cip-args $(RUN_DAG)/args_cip.txt \ --test-args $(RUN_DAG)/args_test.txt --plot-args $(RUN_DAG)/args_plot.txt \ - --request-memory-CIP 2048 --request-memory-ILE 4096 \ + --request-memory-CIP 2048 --request-memory-ILE 4096 --n-copies $(N_COPIES_DAG) \ --input-grid $(RUN_DAG)/overlap-grid.xml.gz --n-samples-per-job $(NPTS_IT_DAG) \ --working-directory $(RUN_DAG) --n-iterations $(N_IT_DAG) \ $(DAG_PILOT_CEPP) From 785e545571ae7f07b868a12bd3554fdee5f5f908 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Sun, 31 May 2026 12:01:36 -0400 Subject: [PATCH 06/18] calmarg: frame zero-cal burn-in of the extrinsic sampler (design) Design framing (per review) for an ILE-level burn-in in analyze_event: cal marginalization -- especially cold-start PRIOR cal on iteration 0 -- makes the extrinsic integral hard to converge (low n_eff), so we fail to seed the grid. The extrinsic posterior is ~cal- independent, so burn the AV sampler in on the cheap ZERO-CAL (n_cal=1) likelihood to a target n_eff first, then switch to the full cal-marg likelihood reusing the adapted proposal. Two mechanisms documented (two-phase integrate; robust warm-start via the existing update_sampling_prior/oracle path) + proposed --calibration-burn-in-neff flag. Composes with the cal pilots (this seeds the EXTRINSIC proposal in-job; pilots seed the CAL proposal across iterations). Implementation tracked separately (needs a GPU test). Co-Authored-By: Claude Opus 4.8 --- .../RIFT/calmarg/DESIGN_adaptive_driver.md | 47 +++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/calmarg/DESIGN_adaptive_driver.md b/MonteCarloMarginalizeCode/Code/RIFT/calmarg/DESIGN_adaptive_driver.md index 4922384c..2330ab7f 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/calmarg/DESIGN_adaptive_driver.md +++ b/MonteCarloMarginalizeCode/Code/RIFT/calmarg/DESIGN_adaptive_driver.md @@ -2,6 +2,53 @@ Status: **planning** (do not implement the multi-stage loop until we pick a path). +## Zero-cal burn-in of the extrinsic sampler (proposed; ILE-level, analyze_event) + +**Problem.** Calibration marginalization makes the extrinsic integral much harder to +converge: with cal drawn from the broad PRIOR the per-time lnL has a large dynamic range, +so the adaptive sampler (AV) struggles to reach a useful `n_eff` -- and on the FIRST +iteration there is no learned cal proposal yet, so every cold-start ILE job is in this +regime. Empirically (local DAG run) iteration-0 points come out with `sigma ~ 0.7-0.9` +and few effective samples; high-SNR sources are hard even WITHOUT cal, and cal makes it +worse. We are effectively failing to seed the *intrinsic* grid because the *extrinsic* +sampler never gets going. + +**Idea (O'Shaughnessy).** Burn the sampler in on a *different, cheaper* likelihood first +-- the ZERO-CAL (n_cal=1) baseline -- until it reaches a minimal `n_eff`, so the +extrinsic sampling proposal is "roughly right", THEN switch to the full cal-marginalized +likelihood for the production estimate. The extrinsic posterior shape is nearly the same +with and without cal (cal mostly rescales / mildly shifts lnL), so the burned-in proposal +is an excellent warm start -- and the zero-cal evaluations are ~`n_cal`x cheaper. + +**Where it lives.** `analyze_event` in `integrate_likelihood_extrinsic_batchmode`. The +likelihood closures already read the module-scope `n_cal_for_likelihood`; the sampler also +already supports a warm start via `sampler.update_sampling_prior(..., external_rvs=...)` +(the existing `oracleRS` path, ~line 2078). Two viable mechanisms: + + 1. **Two-phase integrate (simplest).** Set `n_cal_for_likelihood = 1`, call + `sampler.integrate(likelihood_function, ..., n_eff=burn_in_neff)` (the closures now + evaluate the fast zero-cal baseline), then restore `n_cal_for_likelihood` and call + the production `sampler.integrate(...)` WITHOUT resetting -- reusing the adapted + proposal. Risk: AV's reset semantics across two integrate() calls in one + analyze_event are unverified (it "always resets every iteration" between DAG + iterations; need to confirm it does NOT reset at the start of integrate()). + 2. **Warm-start via update_sampling_prior (robust).** Run the zero-cal burn-in, + harvest its drawn extrinsic samples + lnL, and feed them to + `sampler.update_sampling_prior(external_rvs=...)` exactly like the oracle path, then + run the production integrate. Survives regardless of integrate()'s reset behavior. + +**Proposed flag.** `--calibration-burn-in-neff ` (0/None = off): target n_eff for +the zero-cal burn-in (capped by a fraction of n-max). Default off; opt-in. + +**Relation to the bigger seeding plan.** This is the in-job version of the same idea as +the *zero-cal pilots* that seed the intrinsic/extrinsic grids for high-SNR sources: burn +in cheaply, then pay for cal only once the sampling is on-target. The pilot (`pilot.py`, +util_CalPilotStage) seeds the CAL proposal across iterations; this burn-in seeds the +EXTRINSIC proposal within a single job. They compose. + +Status: framed; needs implementation behind the flag + a GPU smoke test to confirm the +sampler reuses/seeds adaptation across the burn-in -> production handoff. + ## Where we are - In-loop calmarg works and is validated (loop == fused == reference ~1e-14; CPU+GPU; From 4296c55bf9f890b6aa537087ef6c616ff2cbb647 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Sun, 31 May 2026 14:37:47 -0400 Subject: [PATCH 07/18] demo/calmarg: fix wide-ILE settings that crippled n_eff (single-event tuning) Per review, tuned settings on a single-event ILE at the injection (new `tune-single` target) before trusting a full DAG. Found and fixed three issues in the DAG wide-ILE args (all inherited from the ILE-GPU-Paper STANDARD_ILE_OPTS, wrong for calmarg): 1. --no-adapt-distance WITH no distance marginalization -> distance neither adapted nor marginalized -> catastrophic n_eff. Removed the freeze flags; added a DMARG_DAG toggle (default ON: --distance-marginalization + lookup table). 2. --adapt-weight-exponent 0.1 over-tempers the adaptive proposal so it cannot concentrate on the high-dynamic-range cal-marginalized peak. Dropped it -> ILE default 1.0. 3. (separately) the single-event tune must run on the real GPU: cardassia has ONE card (device 0, NVS 510 sm_30); CUDA_VISIBLE_DEVICES must be 0 (3 -> no device -> cupy off -> --gpu silently disabled -> non-GPU path ignores distmarg -> 6-arg likelihood crash). New `make tune-single` runs one ILE at the injection with the exact wide args (lnL col -4, n_eff col -1) -- the fast way to validate settings before an hour-long DAG. Also added --d-min 1 to match the working demo COMMON. N_COPIES default 1, DMARG_DAG default 1. Remaining: even fixed, this SNR~17.5 source + cal gives low single-event n_eff (~3 at 100k) -- high-SNR+cal is intrinsically hard, motivating the zero-cal burn-in (task #24). Co-Authored-By: Claude Opus 4.8 --- .../Code/demo/rift/calmarg/Makefile | 39 ++++++++++++++++--- 1 file changed, 33 insertions(+), 6 deletions(-) diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile b/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile index 04da9628..3c7cc252 100644 --- a/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile +++ b/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile @@ -193,6 +193,17 @@ NEFF_DAG ?= 15 NCHUNK_DAG ?= 10000 N_COPIES_DAG ?= 1 # ILE job duplicates (builder default 2); 1 -> no redundancy, # ~2x faster turnaround on a single-card test +DMARG_DAG ?= 1 # distance marginalization in the wide ILE (default ON). Without + # it (and the removed --no-adapt-distance) distance was neither + # adapted nor marginalized -> catastrophic n_eff. + +ifeq ($(DMARG_DAG),1) + DAG_DMARG_FLAG := --distance-marginalization --distance-marginalization-lookup-table $(CURDIR)/distance_marg.npz + DAG_DMARG_DEP := distance_marg.npz +else + DAG_DMARG_FLAG := + DAG_DMARG_DEP := +endif ifeq ($(FUSED),1) DAG_FUSED_FLAG := --calibration-fused-kernel @@ -205,15 +216,19 @@ else DAG_PILOT_CEPP := endif +# NOTE: adapt-weight-exponent left at the ILE DEFAULT (1.0). Forcing 0.1 (as the +# ILE-GPU-Paper STANDARD_ILE_OPTS does) over-tempers the adaptive proposal -> it cannot +# concentrate on the high-dynamic-range cal-marginalized peak -> n_eff crawls (~2 at +# 100k). At the default 1.0 the same single-event run converges normally. DAG_ILE_ARGS := --n-chunk $(NCHUNK_DAG) --time-marginalization --sim-xml overlap-grid.xml.gz \ - --reference-freq 100.0 --adapt-weight-exponent 0.1 --event-time $(EVENT_TIME) --save-P 0.1 \ + --reference-freq 100.0 --event-time $(EVENT_TIME) --save-P 0.1 \ --cache-file $(CACHE) --fmin-template 10 --n-max $(NMAX_DAG) --fmax 1700.0 --save-deltalnL inf \ - --l-max 2 --n-eff $(NEFF_DAG) --approximant SEOBNRv4 --adapt-floor-level 0.1 --d-max $(DMAX) \ + --l-max 2 --n-eff $(NEFF_DAG) --approximant SEOBNRv4 --d-min 1 --d-max $(DMAX) \ --psd-file H1=$(PSD) --psd-file L1=$(PSD) --psd-file V1=$(PSD) \ --channel-name H1=FAKE-STRAIN --channel-name L1=FAKE-STRAIN --channel-name V1=FAKE-STRAIN \ --inclination-cosine-sampler --declination-cosine-sampler --data-start-time 1000000008 \ - --data-end-time 1000000016 --inv-spec-trunc-time 0 --no-adapt-after-first --no-adapt-distance \ - --srate 4096 --sampler-method AV --vectorized --gpu \ + --data-end-time 1000000016 --inv-spec-trunc-time 0 \ + --srate 4096 --sampler-method AV --vectorized --gpu $(DAG_DMARG_FLAG) \ --calibration-envelope-directory $(CURDIR)/cal_env --calibration-n-realizations $(NCAL_DAG) \ --calibration-spline-count 10 $(DAG_FUSED_FLAG) # --sigma-cut 5: cal draws (PRIOR for vanilla every iteration; PRIOR at iteration 0 for the @@ -225,10 +240,22 @@ DAG_TEST_ARGS := --always-succeed --method lame --parameter m1 DAG_PLOT_ARGS := --parameter m1 --parameter m2 DAG_FILE := marginalize_intrinsic_parameters_BasicIterationWorkflow.dag -.PHONY: dag dag-build dag-validate dag-run +.PHONY: dag dag-build dag-validate dag-run tune-single + +# Single-event ILE at the INJECTION (command-single style) using the SAME wide-ILE args +# the DAG will use -- the fast way to tune settings (n_eff!) for the event before paying +# for a full multi-iteration DAG. Reports lnL (col -4) and n_eff (col -1). +# make tune-single FUSED=1 DMARG_DAG=1 NMAX_DAG=200000 NEFF_DAG=200 +tune-single: cal_env/H1.txt $(DAG_DMARG_DEP) + rm -rf $(CURDIR)/rundir_tune; mkdir -p $(CURDIR)/rundir_tune + cd $(CURDIR)/rundir_tune && $(ENV) $(ILE) $(DAG_ILE_ARGS) \ + --sim-xml $(SIM_XML) --n-events-to-analyze 1 --output-file tune_out + @echo "=== tune-single result (lnL +- err, n_eff) ===" + @cat $(CURDIR)/rundir_tune/tune_out*.dat + @echo " (n_eff is the last column; aim for a healthy value before running the DAG)" # Build the DAG (no condor submission) and validate the wiring. PILOT/FUSED toggles above. -dag-build: cal_env/H1.txt +dag-build: cal_env/H1.txt $(DAG_DMARG_DEP) rm -rf $(RUN_DAG); mkdir -p $(RUN_DAG) cd $(RUN_DAG) && $(ENV) util_ManualOverlapGrid.py --parameter mc --parameter-range '[23,35]' \ --parameter delta_mc --parameter-range '[0.0,0.5]' --grid-cartesian-npts $(NPTS_GRID_DAG) \ From e333ff926148382ca4da19505126b462ca19ea7d Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Sun, 31 May 2026 15:25:15 -0400 Subject: [PATCH 08/18] demo/calmarg: correct misleading tuning comments (record accurate facts) Correcting mischaracterizations from the previous commit (per review), to avoid misleading future readers: - --no-adapt-distance is NOT a freeze: distance is then UNIFORMLY sampled (fine for low SNR with sensible d-min/d-max), and is moot once distance marginalization is on (the default here). It did not "cripple n_eff". - --adapt-weight-exponent is a NO-OP for the AV sampler (it matters for GMM/portfolio -- not dropped globally, just not set in the demo). It did not change convergence here. - SNR ~17.5 is moderate, not loud (loud = 40+); 100k samples is SHORT for RIFT (production runs use millions and let AV creep n_eff up). n_eff ~3 at 100k is early, not pathology. The valid finding stands: baseline (no cal) and calmarg give ~equal n_eff, so cal is not the bottleneck. DMARG_DAG=1 (analytic distance) remains the clean default. Co-Authored-By: Claude Opus 4.8 --- .../Code/demo/rift/calmarg/Makefile | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile b/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile index 3c7cc252..73a37df5 100644 --- a/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile +++ b/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile @@ -193,9 +193,11 @@ NEFF_DAG ?= 15 NCHUNK_DAG ?= 10000 N_COPIES_DAG ?= 1 # ILE job duplicates (builder default 2); 1 -> no redundancy, # ~2x faster turnaround on a single-card test -DMARG_DAG ?= 1 # distance marginalization in the wide ILE (default ON). Without - # it (and the removed --no-adapt-distance) distance was neither - # adapted nor marginalized -> catastrophic n_eff. +DMARG_DAG ?= 1 # distance marginalization in the wide ILE (default ON). With it, + # distance is integrated analytically (not sampled), the clean + # production path. (Aside: --no-adapt-distance is NOT a freeze -- + # distance is then UNIFORMLY sampled, fine for low-SNR with a + # sensible d-min/d-max; it is simply moot once distmarg is on.) ifeq ($(DMARG_DAG),1) DAG_DMARG_FLAG := --distance-marginalization --distance-marginalization-lookup-table $(CURDIR)/distance_marg.npz @@ -216,10 +218,12 @@ else DAG_PILOT_CEPP := endif -# NOTE: adapt-weight-exponent left at the ILE DEFAULT (1.0). Forcing 0.1 (as the -# ILE-GPU-Paper STANDARD_ILE_OPTS does) over-tempers the adaptive proposal -> it cannot -# concentrate on the high-dynamic-range cal-marginalized peak -> n_eff crawls (~2 at -# 100k). At the default 1.0 the same single-event run converges normally. +# NOTE: adapt-weight-exponent left at the ILE DEFAULT (1.0). For the AV sampler used here +# this option is a NO-OP (it matters for the OTHER integrators -- GMM/portfolio -- so do +# not drop it globally). n_eff is naturally low at small n-max (~3 at 100k) regardless: +# 100k samples is SHORT by RIFT standards; production runs use millions and let the AV +# integrator creep n_eff upward. Baseline (no cal) and calmarg give ~equal n_eff here, so +# cal is not the bottleneck -- the moderate SNR (~17.5; "loud" is 40+) x sample count is. DAG_ILE_ARGS := --n-chunk $(NCHUNK_DAG) --time-marginalization --sim-xml overlap-grid.xml.gz \ --reference-freq 100.0 --event-time $(EVENT_TIME) --save-P 0.1 \ --cache-file $(CACHE) --fmin-template 10 --n-max $(NMAX_DAG) --fmax 1700.0 --save-deltalnL inf \ From cf16b030738bbc2ba8a6c7eea580a9b5d0cbd170 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Sun, 31 May 2026 15:29:53 -0400 Subject: [PATCH 09/18] calmarg: zero-cal burn-in of the extrinsic sampler (--calibration-burn-in-neff) Implements task #24 (the generally-useful extrinsic-seeding idea). Cal marginalization -- especially cold-start prior cal -- makes the extrinsic integral slow to converge; the extrinsic posterior is ~cal-independent, so burn the sampler in on the cheap ZERO-CAL (n_cal=1) likelihood first, then run the full cal-marginalized integral reusing the adapted proposal. ILE (analyze_event): --calibration-burn-in-neff (+ --calibration-burn-in-nmax cap). When set and calmarg active, before the production sampler.integrate it toggles the analyze_event-local n_cal_for_likelihood to 1 (the likelihood closures read it, so the SAME like_to_integrate now evaluates the fast baseline), runs a burn-in integrate to the target n_eff with a capped nmax, then restores n_cal_for_likelihood and runs the full cal-marg integral. Correctness is preserved regardless of whether the AV sampler retains adaptation across the two integrate() calls -- worst case the burn-in is wasted; the production integral is always the full cal one. Default off. Threaded --calmarg-burn-in-neff through util_RIFT_pseudo_pipe.py (-> args_ile.txt) and a BURN_IN_NEFF toggle in the demo Makefile (tune-single/DAG). Compiles; flags register; demo emits the flag. NEEDS a GPU smoke test of the burn-in->production handoff (does AV reuse/seed adaptation across the two integrates; else add the update_sampling_prior warm-start documented in DESIGN_adaptive_driver.md). Co-Authored-By: Claude Opus 4.8 --- .../integrate_likelihood_extrinsic_batchmode | 23 +++++++++++++++++++ .../Code/bin/util_RIFT_pseudo_pipe.py | 3 +++ .../Code/demo/rift/calmarg/Makefile | 9 +++++++- 3 files changed, 34 insertions(+), 1 deletion(-) diff --git a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode index e34fd484..6fa2f33f 100755 --- a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode +++ b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode @@ -241,6 +241,8 @@ optp.add_option("--calibration-fused-kernel", action="store_true", default=False optp.add_option("--calibration-proposal-breadcrumb",default=None, help="Opt-in (Option C / adaptive pilot): path to a breadcrumb .npz (RIFT.calmarg.breadcrumbs) carrying a LEARNED Gaussian proposal over cal spline nodes. When set, the cal realizations are drawn from that proposal instead of the broad prior, and the marginalization carries Phase-0 importance weights log(prior/proposal) so it stays unbiased. Requires --calibration-envelope-directory (the prior).") optp.add_option("--calibration-dump-responsibilities",default=None, help="Opt-in (Option C / adaptive pilot): path to write per-cal-realization log-responsibilities (length n_cal), accumulated over the evaluated grid, plus the cal node draws. This is the pilot's output, fitted into a proposal by util_CalPilotFit.py. No effect on the returned likelihood.") optp.add_option("--calibration-pilot-extrinsic",default=256,type=int, help="Pilot only: number of uniform-prior extrinsic samples used to extrinsic-marginalize the per-realization cal responsibility at each intrinsic point. Cal is ~extrinsic-independent, so a modest batch suffices.") +optp.add_option("--calibration-burn-in-neff",default=None,type=float, help="Opt-in (generally useful): before the production cal-marginalized integration, BURN IN the extrinsic sampler on the cheap ZERO-CAL (n_cal=1) likelihood until this effective sample count, then switch to the full cal-marginalized likelihood reusing the adapted proposal. The extrinsic posterior is ~cal-independent, so this warm-starts the (expensive) cal-marg phase. No effect unless calibration marginalization is active.") +optp.add_option("--calibration-burn-in-nmax",default=None,type=int, help="Cap on the number of samples drawn during the zero-cal burn-in (default: the run's --n-max). Keeps the burn-in bounded if it cannot reach --calibration-burn-in-neff.") optp.add_option("--vectorized", action="store_true", help="Perform manipulations of lm and timeseries using numpy arrays, not LAL data structures. (Combine with --gpu to enable GPU use, where available)") optp.add_option("--gpu", action="store_true", help="Perform manipulations of lm and timeseries using numpy arrays, CONVERTING TO GPU when available. You MUST use this option with --vectorized (otherwise it is a no-op). You MUST have a suitable version of cupy installed, your cuda operational, etc") optp.add_option("--force-gpu-only", action="store_true", help="Hard fail if no GPU present (assessed by cupy not loading)") @@ -2086,6 +2088,27 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts,inv_spec_t 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) + # Optional ZERO-CAL burn-in (generally useful; see RIFT/calmarg/DESIGN_adaptive_driver.md). + # Adapt the extrinsic sampler cheaply on the n_cal=1 baseline first, then run the full + # cal-marginalized integration reusing the adapted proposal. The likelihood closures + # read the enclosing n_cal_for_likelihood, so toggling it to 1 makes this same + # like_to_integrate evaluate the fast baseline. Correctness is preserved regardless of + # whether the sampler retains adaptation across the two integrate() calls -- worst case + # the burn-in is simply wasted; the production integral below is always the full cal one. + if opts.calibration_burn_in_neff and calibration_marginalization and n_cal_for_likelihood > 1: + _ncal_full = n_cal_for_likelihood + n_cal_for_likelihood = 1 + _burn_pinned = dict(pinned_params) + _burn_pinned['neff'] = float(opts.calibration_burn_in_neff) + _burn_pinned['nmax'] = int(opts.calibration_burn_in_nmax) if opts.calibration_burn_in_nmax else int(pinned_params.get('nmax', n_max)) + print(" [calmarg burn-in] adapting extrinsic sampler on ZERO-CAL likelihood -> neff>={:g} (nmax cap {})".format(_burn_pinned['neff'], _burn_pinned['nmax'])) + try: + _b = sampler.integrate(like_to_integrate, *unpinned_params, **_burn_pinned) + print(" [calmarg burn-in] done (burn-in neff ~ {}); switching to full cal marginalization".format(_b[2] if _b and len(_b) > 2 else '?')) + except Exception as _eb: + print(" [calmarg burn-in] integrate failed ({}); proceeding to production".format(_eb)) + n_cal_for_likelihood = _ncal_full # restore: production uses the full cal set + res, var, neff, dict_return = sampler.integrate(like_to_integrate, *unpinned_params, **pinned_params) if not(res): # no resut diff --git a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py index e8d0800f..c3ca93ce 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py @@ -166,6 +166,7 @@ def unsafe_parse_arg_string_dict(my_argstr): parser.add_argument("--calmarg-pilot-top-fraction",default=0.05,type=float,help="Fraction of highest-lnL composite points the pilot harvests. Default 0.05.") parser.add_argument("--calmarg-pilot-max-points",default=32,type=int,help="Cap on harvested pilot points per iteration. Default 32.") parser.add_argument("--calmarg-first-cip-sigma-cut",default=100.0,type=float,help="With --calmarg-pilot: relax the first CIP stage's --sigma-cut to this value, so cold-start (prior-cal) iteration-0 points -- which have large MC error -- are not all stripped by CIP's default 0.6. Threaded to helper_LDG_Events.py. Default 100 (effectively keep all cold-start points).") +parser.add_argument("--calmarg-burn-in-neff",default=None,type=float,help="In-loop calmarg: burn the extrinsic sampler in on the cheap zero-cal likelihood to this n_eff before the full cal-marginalized integration (warm start; the extrinsic posterior is ~cal-independent). Threaded to ILE as --calibration-burn-in-neff.") parser.add_argument("--distance-reweighting",action='store_true',help="Option to add job to DAG to reweight posterior samples due to different distance prior (LVK prod prior)") parser.add_argument("--extra-args-helper",action=None, help="Filename with arguments for the helper. Use to provide alternative channel names and other advanced configuration (--channel-name, data type)!") parser.add_argument("--manual-postfix",default='',type=str) @@ -1045,6 +1046,8 @@ def unsafe_parse_arg_string_dict(my_argstr): line += " --calibration-envelope-directory {} --calibration-n-realizations {} --calibration-spline-count {} ".format(cal_dir, opts.calmarg_n_realizations, opts.calmarg_spline_count) if opts.calmarg_fused_kernel: line += " --calibration-fused-kernel " + if opts.calmarg_burn_in_neff: + line += " --calibration-burn-in-neff {} ".format(opts.calmarg_burn_in_neff) if opts.calmarg_pilot: # Option C: wide ILE jobs are SEEDED from the previous iteration's consolidated cal # proposal. The $(macroiterationprev) condor macro resolves per node; ILE falls diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile b/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile index 73a37df5..8c4acf1b 100644 --- a/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile +++ b/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile @@ -212,6 +212,13 @@ ifeq ($(FUSED),1) else DAG_FUSED_FLAG := endif +# Optional zero-cal extrinsic burn-in (task #24): set BURN_IN_NEFF to a target n_eff to +# adapt the sampler on the cheap zero-cal likelihood before the full cal-marg integral. +ifneq ($(BURN_IN_NEFF),) + DAG_BURNIN_FLAG := --calibration-burn-in-neff $(BURN_IN_NEFF) +else + DAG_BURNIN_FLAG := +endif ifeq ($(PILOT),1) DAG_PILOT_CEPP := --calmarg-pilot --calmarg-pilot-cadence 1 --calmarg-pilot-max-it 1 else @@ -234,7 +241,7 @@ DAG_ILE_ARGS := --n-chunk $(NCHUNK_DAG) --time-marginalization --sim-xml overlap --data-end-time 1000000016 --inv-spec-trunc-time 0 \ --srate 4096 --sampler-method AV --vectorized --gpu $(DAG_DMARG_FLAG) \ --calibration-envelope-directory $(CURDIR)/cal_env --calibration-n-realizations $(NCAL_DAG) \ - --calibration-spline-count 10 $(DAG_FUSED_FLAG) + --calibration-spline-count 10 $(DAG_FUSED_FLAG) $(DAG_BURNIN_FLAG) # --sigma-cut 5: cal draws (PRIOR for vanilla every iteration; PRIOR at iteration 0 for the # pilot cold start) have high per-point MC error (~0.7-0.9 here), which the CIP default # --sigma-cut 0.6 would strip entirely. Relax it (cf. the helper_LDG_Events.py fix that From b75ae0c3d75790a797f339bd062d1e0477f50838 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Sun, 31 May 2026 15:43:18 -0400 Subject: [PATCH 10/18] calmarg burn-in: breadcrumb AV reset limitation + future seedable/flexible-AV work Per review: AV (the default, most-efficient extrinsic sampler) COMPLETELY RESETS between integrate() calls -- no seedable AV exists -- so the zero-cal burn-in gives AV no speedup (correctness-safe overhead only). Re-seeding AV is also dangerous: AV can only CONTRACT its volume, never EXPAND or SHIFT boundaries, so an off/too-tight warm start would trap the production phase. GMM/portfolio can reuse sampling models (update_sampling_prior/gmm_dict) but are less efficient. Breadcrumbed in DESIGN_adaptive_driver.md + the --calibration-burn-in-neff help: the burn-in is parked (gated, harmless, ready) pending future work that is broadly useful beyond calmarg -- (1) a seedable AV, (2) a boundary-shifting AV that can expand/translate its volume, not only contract. Until then the cal PILOT (across-iteration proposal learning) + the prior-shrinkage backstop are the load-bearing cal pieces. Co-Authored-By: Claude Opus 4.8 --- .../RIFT/calmarg/DESIGN_adaptive_driver.md | 25 +++++++++++++++++-- .../integrate_likelihood_extrinsic_batchmode | 2 +- 2 files changed, 24 insertions(+), 3 deletions(-) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/calmarg/DESIGN_adaptive_driver.md b/MonteCarloMarginalizeCode/Code/RIFT/calmarg/DESIGN_adaptive_driver.md index 2330ab7f..a1f1db98 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/calmarg/DESIGN_adaptive_driver.md +++ b/MonteCarloMarginalizeCode/Code/RIFT/calmarg/DESIGN_adaptive_driver.md @@ -46,8 +46,29 @@ in cheaply, then pay for cal only once the sampling is on-target. The pilot (`p util_CalPilotStage) seeds the CAL proposal across iterations; this burn-in seeds the EXTRINSIC proposal within a single job. They compose. -Status: framed; needs implementation behind the flag + a GPU smoke test to confirm the -sampler reuses/seeds adaptation across the burn-in -> production handoff. +Status: implemented behind `--calibration-burn-in-neff` (two-phase integrate, toggling +`n_cal_for_likelihood`), correctness-safe (the production integral is always the full cal +one). BUT see the sampler limitation below. + +### Sampler limitation (measured / per review) -- BREADCRUMB for future work +The default and most efficient extrinsic sampler, **AV (mcsamplerAdaptiveVolume), +COMPLETELY RESETS between `integrate()` calls** -- there is currently no seedable AV. So +the two-phase burn-in gives AV **no speedup** (the adapted proposal is thrown away); it is +only ever correctness-safe overhead there. Moreover, **re-seeding AV is inherently +dangerous**: AV can only *contract* its volume, never *expand* or *shift* its boundaries, +so a burn-in proposal that is slightly off / too tight would trap the production phase in +the wrong region. The other integrators (**GMM**, **portfolio**) CAN reuse sampling +models (their `update_sampling_prior` / `gmm_dict` hooks) but are less efficient overall. + +Therefore the zero-cal burn-in only pays off once we have: + 1. a **seedable AV** (warm-start AV from external samples / a prior proposal), and/or + 2. a **more flexible, boundary-shifting AV** that can EXPAND and TRANSLATE its sampling + volume (not only contract) -- so a warm start can't trap it. +Both would be broadly useful well beyond calmarg (every iteration, every warm start). +Until then: keep the burn-in flag (harmless, gated, ready) but do NOT rely on it for AV; +for GMM/portfolio it can warm-start via their model-reuse hooks (untested). The cal PILOT +(across-iteration proposal learning) and the prior-shrinkage backstop remain the load- +bearing pieces for cal; extrinsic seeding waits on seedable AV. ## Where we are diff --git a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode index 6fa2f33f..91e1fb45 100755 --- a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode +++ b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode @@ -241,7 +241,7 @@ optp.add_option("--calibration-fused-kernel", action="store_true", default=False optp.add_option("--calibration-proposal-breadcrumb",default=None, help="Opt-in (Option C / adaptive pilot): path to a breadcrumb .npz (RIFT.calmarg.breadcrumbs) carrying a LEARNED Gaussian proposal over cal spline nodes. When set, the cal realizations are drawn from that proposal instead of the broad prior, and the marginalization carries Phase-0 importance weights log(prior/proposal) so it stays unbiased. Requires --calibration-envelope-directory (the prior).") optp.add_option("--calibration-dump-responsibilities",default=None, help="Opt-in (Option C / adaptive pilot): path to write per-cal-realization log-responsibilities (length n_cal), accumulated over the evaluated grid, plus the cal node draws. This is the pilot's output, fitted into a proposal by util_CalPilotFit.py. No effect on the returned likelihood.") optp.add_option("--calibration-pilot-extrinsic",default=256,type=int, help="Pilot only: number of uniform-prior extrinsic samples used to extrinsic-marginalize the per-realization cal responsibility at each intrinsic point. Cal is ~extrinsic-independent, so a modest batch suffices.") -optp.add_option("--calibration-burn-in-neff",default=None,type=float, help="Opt-in (generally useful): before the production cal-marginalized integration, BURN IN the extrinsic sampler on the cheap ZERO-CAL (n_cal=1) likelihood until this effective sample count, then switch to the full cal-marginalized likelihood reusing the adapted proposal. The extrinsic posterior is ~cal-independent, so this warm-starts the (expensive) cal-marg phase. No effect unless calibration marginalization is active.") +optp.add_option("--calibration-burn-in-neff",default=None,type=float, help="Opt-in: before the production cal-marginalized integration, BURN IN the extrinsic sampler on the cheap ZERO-CAL (n_cal=1) likelihood until this effective sample count, then switch to the full cal-marginalized likelihood. The extrinsic posterior is ~cal-independent. CAVEAT: the AV sampler RESETS between integrate() calls (no seedable AV yet), so this gives AV no speedup (correctness-safe only). It can warm-start GMM/portfolio (model reuse). Awaiting a seedable / boundary-shifting AV; see DESIGN_adaptive_driver.md. No effect unless calmarg is active.") optp.add_option("--calibration-burn-in-nmax",default=None,type=int, help="Cap on the number of samples drawn during the zero-cal burn-in (default: the run's --n-max). Keeps the burn-in bounded if it cannot reach --calibration-burn-in-neff.") optp.add_option("--vectorized", action="store_true", help="Perform manipulations of lm and timeseries using numpy arrays, not LAL data structures. (Combine with --gpu to enable GPU use, where available)") optp.add_option("--gpu", action="store_true", help="Perform manipulations of lm and timeseries using numpy arrays, CONVERTING TO GPU when available. You MUST use this option with --vectorized (otherwise it is a no-op). You MUST have a suitable version of cupy installed, your cuda operational, etc") From d0497a815f68c1e96048db079665ec6f2dfc44fa Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Sun, 31 May 2026 19:08:25 -0400 Subject: [PATCH 11/18] demo/calmarg: tune-condor target (single condor job, unbuffered) for n_eff creep Per review: a full DAG is overkill for single-event settings tuning, and an interactive background launch buffers stdout (can't watch n_eff creep) and is subject to the NVS 510 display watchdog. tune-condor writes ONE submit file running `python -u --sim-xml ` with request_GPUs=1 + getenv (cupy), persistent output/error/log, initialdir in rundir_tune. Unbuffered -> live progress; condor -> watchdog-isolated + restartable. Verified: submits, runs on the GPU, streams the iteration/Neff table live. conda run -n rift_gpu2 make tune-condor FUSED=1 DMARG_DAG=1 NMAX_DAG=4000000 tail -f rundir_tune/tune_condor.out ; cat rundir_tune/tune_out*.dat (lnL col -4, neff -1) Co-Authored-By: Claude Opus 4.8 --- .../Code/demo/rift/calmarg/Makefile | 22 ++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile b/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile index 8c4acf1b..5b6c3bb2 100644 --- a/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile +++ b/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile @@ -251,7 +251,27 @@ DAG_TEST_ARGS := --always-succeed --method lame --parameter m1 DAG_PLOT_ARGS := --parameter m1 --parameter m2 DAG_FILE := marginalize_intrinsic_parameters_BasicIterationWorkflow.dag -.PHONY: dag dag-build dag-validate dag-run tune-single +.PHONY: dag dag-build dag-validate dag-run tune-single tune-condor + +# Single-event ILE at the injection as ONE condor job (a full DAG is overkill for tuning): +# persistent logs, watchdog-isolated (condor, not an interactive launch), restartable, and +# UNBUFFERED via `python -u` so n_eff creep is watchable live. Submit from inside the env +# (getenv captures cupy): conda run -n rift_gpu2 make tune-condor FUSED=1 NMAX_DAG=4000000 +# watch: tail -f rundir_tune/tune_condor.out ; cat rundir_tune/tune_out*.dat +tune-condor: cal_env/H1.txt $(DAG_DMARG_DEP) + rm -rf $(CURDIR)/rundir_tune; mkdir -p $(CURDIR)/rundir_tune + @PYEXE=$$(command -v python); \ + { echo "universe = vanilla"; \ + echo "executable = $$PYEXE"; \ + echo "arguments = -u $(ILE) $(DAG_ILE_ARGS) --sim-xml $(SIM_XML) --n-events-to-analyze 1 --output-file tune_out"; \ + echo "request_GPUs = 1"; \ + echo "getenv = True"; \ + echo "initialdir = $(CURDIR)/rundir_tune"; \ + echo "output = tune_condor.out"; echo "error = tune_condor.err"; echo "log = tune_condor.log"; \ + echo "queue"; } > $(CURDIR)/rundir_tune/tune_condor.sub + cd $(CURDIR)/rundir_tune && $(ENV) condor_submit tune_condor.sub + @echo "Submitted. Watch live: tail -f $(CURDIR)/rundir_tune/tune_condor.out" + @echo "Result when done: cat $(CURDIR)/rundir_tune/tune_out*.dat (lnL col -4, n_eff col -1)" # Single-event ILE at the INJECTION (command-single style) using the SAME wide-ILE args # the DAG will use -- the fast way to tune settings (n_eff!) for the event before paying From 702569fe3bb9e6308f0bc69672413f467c7c1adc Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Sun, 31 May 2026 19:11:04 -0400 Subject: [PATCH 12/18] Use portable RIFT float dtype in likelihood tools Thread RiftFloat through the ILE executables, related waveform and posterior utilities, and likelihood sampler tests so platforms without numpy.float128 use the centralized float64 fallback. --- .../Code/bin/integrate_likelihood_extrinsic | 11 ++++++----- .../Code/bin/integrate_likelihood_extrinsic_batchmode | 11 ++++++----- ..._ConstructIntrinsicPosterior_GenericCoordinates.py | 4 ++-- .../Code/bin/util_GenerateMaxlnLWaveform.py | 5 +++-- .../bin/util_GenerateMaxlnLWaveform_NRFromIndex.py | 5 +++-- .../Code/test/test_like_and_samp.py | 7 ++++--- .../Code/test/test_like_and_samp_simplified.py | 4 ++-- 7 files changed, 26 insertions(+), 21 deletions(-) diff --git a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic index 347e0bce..e55072b3 100755 --- a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic +++ b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic @@ -39,6 +39,7 @@ import glue.lal #import pylal import RIFT.lalsimutils as lalsimutils +from RIFT.precision import RiftFloat import RIFT.likelihood.factored_likelihood as factored_likelihood import RIFT.integrators.mcsampler as mcsampler import RIFT.misc.sky_rotations as sky_rotations @@ -1047,7 +1048,7 @@ if not opts.time_marginalization: incl = numpy.arccos(incl) # use EXTREMELY many bits - lnL = numpy.zeros(right_ascension.shape,dtype=numpy.float128) + lnL = numpy.zeros(right_ascension.shape,dtype=RiftFloat) i = 0 for ph, th, tr, phr, ic, ps, di in zip(right_ascension, dec, t_ref, phi_orb, incl, psi, distance): @@ -1078,7 +1079,7 @@ else: # Sum over time for every point in other extrinsic params incl = numpy.arccos(incl) # use EXTREMELY many bits - lnL = numpy.zeros(right_ascension.shape,dtype=numpy.float128) + lnL = numpy.zeros(right_ascension.shape,dtype=RiftFloat) i = 0 tvals = numpy.linspace(-t_ref_wind,t_ref_wind,int((t_ref_wind)*2/P.deltaT)) # choose an array at the target sampling rate. P is inherited globally @@ -1140,7 +1141,7 @@ else: # Sum over time for every point in other extrinsic params incl = numpy.arccos(incl) # use EXTREMELY many bits - lnL = numpy.zeros(right_ascension.shape,dtype=numpy.float128) + lnL = numpy.zeros(right_ascension.shape,dtype=RiftFloat) P.phi = right_ascension.astype(float) # cast to float P.theta = dec #declination.astype(float) P.tref = float(fiducial_epoch) @@ -1278,7 +1279,7 @@ else: # Sum over time for every point in other extrinsic params if opts.inclination_cosine_sampler: incl = numpy.arccos(incl) - lnL = numpy.zeros(len(right_ascension),dtype=numpy.float128) + lnL = numpy.zeros(len(right_ascension),dtype=RiftFloat) i = 0 tvals = numpy.linspace(-t_ref_wind,t_ref_wind,int((t_ref_wind)*2/P.deltaT)) # choose an array at the target sampling rate. P is inherited globally @@ -1378,7 +1379,7 @@ if opts.maximize_only and opts.output_file: sampler._rvs["inclination"][indx_guess]/(numpy.pi),\ sampler._rvs["psi"][indx_guess]/numpy.pi,\ sampler._rvs["distance"][indx_guess]/dmax\ - ],dtype=numpy.float128) + ],dtype=RiftFloat) x0 = numpy.fmod(x0,numpy.ones(len(x0))) # had BETTER be defined on this range! # Pick the best starting time. BRUTE FORCE METHOD: use grid def fn_scaled_t(t,x0): diff --git a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode index 64c9e88d..59f13d94 100755 --- a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode +++ b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode @@ -48,6 +48,7 @@ from igwn_ligolw import utils, ligolw import glue.lal import RIFT.lalsimutils as lalsimutils +from RIFT.precision import RiftFloat import RIFT.integrators.mcsampler as mcsampler import RIFT.misc.sky_rotations as sky_rotations try: @@ -1574,7 +1575,7 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts,inv_spec_t distance = redshift_to_distance(distance) # use EXTREMELY many bits - lnL = numpy.zeros(right_ascension.shape,dtype=numpy.float128) + lnL = numpy.zeros(right_ascension.shape,dtype=RiftFloat) i = 0 for ph, th, tr, phr, ic, ps, di in zip(right_ascension, dec, t_ref, phi_orb, incl, psi, distance): @@ -1607,7 +1608,7 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts,inv_spec_t distance = redshift_to_distance(distance) # use EXTREMELY many bits - lnL = numpy.zeros(right_ascension.shape,dtype=numpy.float128) + lnL = numpy.zeros(right_ascension.shape,dtype=RiftFloat) i = 0 tvals = numpy.linspace(-t_ref_wind,t_ref_wind,int((t_ref_wind)*2/P.deltaT)) # choose an array at the target sampling rate. P is inherited globally @@ -1650,7 +1651,7 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts,inv_spec_t distance = redshift_to_distance(distance) # use EXTREMELY many bits - lnL = numpy.zeros(right_ascension.shape,dtype=numpy.float128) + lnL = numpy.zeros(right_ascension.shape,dtype=RiftFloat) P.phi = right_ascension.astype(float) # cast to float P.theta = dec #declination.astype(float) P.tref = float(fiducial_epoch) @@ -1889,7 +1890,7 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts,inv_spec_t if opts.d_prior_redshift: distance = redshift_to_distance(distance) - lnL = numpy.zeros(len(right_ascension),dtype=numpy.float128) + lnL = numpy.zeros(len(right_ascension),dtype=RiftFloat) # i = 0 tvals = numpy.linspace(-t_ref_wind,t_ref_wind,int((t_ref_wind)*2/P.deltaT)) # choose an array at the target sampling rate. P is inherited globally @@ -2360,7 +2361,7 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts,inv_spec_t sampler._rvs["inclination"][indx_guess]/(numpy.pi),\ sampler._rvs["psi"][indx_guess]/numpy.pi,\ sampler._rvs["distance"][indx_guess]/dmax\ - ],dtype=numpy.float128) + ],dtype=RiftFloat) x0 = numpy.fmod(x0,numpy.ones(len(x0))) # had BETTER be defined on this range! # Pick the best starting time. BRUTE FORCE METHOD: use grid def fn_scaled_t(t,x0): diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py index aaa43560..a12eb885 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py @@ -16,6 +16,7 @@ import RIFT.interpolators.BayesianLeastSquares as BayesianLeastSquares +from RIFT.precision import RiftFloat import argparse import sys @@ -3115,7 +3116,7 @@ def parse_corr_params(my_str): # Note also downselects NOT applied: no range cuts, unless applied as part of aligned_prior, etc. # - use for Bayes factors with GREAT CARE for this reason; should correct for with indx_ok log_res_reweighted = lnLmax + np.log(np.mean(weights)) -sigma_reweighted= np.std(weights,dtype=np.float128)/np.mean(weights) +sigma_reweighted= np.std(weights,dtype=RiftFloat)/np.mean(weights) neff_reweighted = np.sum(weights)/np.max(weights) np.savetxt(opts.fname_output_integral+"_withpriorchange.dat", [log_res_reweighted]) # should agree with the usual result, if no prior changes with open(opts.fname_output_integral+"_withpriorchange+annotation.dat", 'w') as file_out: @@ -3693,4 +3694,3 @@ def parse_corr_params(my_str): sys.exit(0) - diff --git a/MonteCarloMarginalizeCode/Code/bin/util_GenerateMaxlnLWaveform.py b/MonteCarloMarginalizeCode/Code/bin/util_GenerateMaxlnLWaveform.py index b95eb26d..49cfbee6 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_GenerateMaxlnLWaveform.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_GenerateMaxlnLWaveform.py @@ -15,10 +15,12 @@ import subprocess from shutil import copyfile +from RIFT.precision import RiftFloat + parser = argparse.ArgumentParser() parser.add_argument("--l-max", type=int, default=2, help="Include all (l,m) modes with l less than or equal to this value.") parser.add_argument("--run-dir",default=None, help="directory code was run in") -parser.add_argument("--event-time", type=np.float128, default=None,help="Event time. If not present, will use event.log") +parser.add_argument("--event-time", type=RiftFloat, default=None,help="Event time. If not present, will use event.log") parser.add_argument("--save-plots", default=False, action='store_true',help="saves waveform plot and hoft data") parser.add_argument("--use-NR", default=False, action='store_true', help="generate plots using NR data") parser.add_argument("--no-memory",action='store_true') @@ -250,4 +252,3 @@ os.chdir(opts.run_dir) print(cmd) os.system(cmd) - diff --git a/MonteCarloMarginalizeCode/Code/bin/util_GenerateMaxlnLWaveform_NRFromIndex.py b/MonteCarloMarginalizeCode/Code/bin/util_GenerateMaxlnLWaveform_NRFromIndex.py index 6424445a..2dcd4061 100644 --- a/MonteCarloMarginalizeCode/Code/bin/util_GenerateMaxlnLWaveform_NRFromIndex.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_GenerateMaxlnLWaveform_NRFromIndex.py @@ -15,6 +15,8 @@ import subprocess from shutil import copyfile +from RIFT.precision import RiftFloat + import NRWaveformCatalogManager3 as nrwf import lal import RIFT.lalsimutils as lalsimutils @@ -26,7 +28,7 @@ parser.add_argument("--nr-param",default=None,help="param") parser.add_argument("--fname-indexed",default=None,help="File name of *.indexed file") parser.add_argument("--n-max",default=1e4,type=int,help="Override settings in command-single.sh") -parser.add_argument("--event-time", type=np.float128, default=None,help="Event time. If not present, will use event.log") +parser.add_argument("--event-time", type=RiftFloat, default=None,help="Event time. If not present, will use event.log") parser.add_argument("--save-plots", default=False, action='store_true',help="saves waveform plot and hoft data") parser.add_argument("--use-NR", default=False, action='store_true', help="generate plots using NR data") parser.add_argument("--no-memory",action='store_true') @@ -180,4 +182,3 @@ os.chdir(opts.run_dir) print(cmd) os.system(cmd) - diff --git a/MonteCarloMarginalizeCode/Code/test/test_like_and_samp.py b/MonteCarloMarginalizeCode/Code/test/test_like_and_samp.py index 895a764e..d8eb6809 100755 --- a/MonteCarloMarginalizeCode/Code/test/test_like_and_samp.py +++ b/MonteCarloMarginalizeCode/Code/test/test_like_and_samp.py @@ -81,6 +81,7 @@ import pickle import numpy as np +from RIFT.precision import RiftFloat try: import cupy @@ -846,7 +847,7 @@ def likelihood_function(right_ascension, declination,t_ref, phi_orb, inclination global nEvals global lnLOffsetValue # use EXTREMELY many bits - lnL = np.zeros(right_ascension.shape,dtype=np.float128) + lnL = np.zeros(right_ascension.shape,dtype=RiftFloat) i = 0 # if opts.rotate_sky_coordinates: # print " -Sky ring width ", np.std(declination), " note contribution from floor is of order p_floor*(pi)/sqrt(12) ~ 0.9 pfloor" @@ -887,7 +888,7 @@ def likelihood_function(right_ascension, declination,t_ref, phi_orb, inclination global nEvals global lnLOffsetValue # use EXTREMELY many bits - lnL = np.zeros(right_ascension.shape,dtype=np.float128) + lnL = np.zeros(right_ascension.shape,dtype=RiftFloat) i = 0 # if opts.rotate_sky_coordinates: # print " -Sky ring width ", np.std(declination), " note contribution from floor is of order p_floor*(pi)/sqrt(12) ~ 0.9 pfloor" @@ -964,7 +965,7 @@ def likelihood_function(right_ascension, declination,t_ref, phi_orb, inclination global nEvals global lnLOffsetValue # use EXTREMELY many bits - lnL = np.zeros(right_ascension.shape,dtype=np.float128) + lnL = np.zeros(right_ascension.shape,dtype=RiftFloat) i = 0 # if opts.rotate_sky_coordinates: # print " -Sky ring width ", np.std(declination), " note contribution from floor is of order p_floor*(pi)/sqrt(12) ~ 0.9 pfloor" diff --git a/MonteCarloMarginalizeCode/Code/test/test_like_and_samp_simplified.py b/MonteCarloMarginalizeCode/Code/test/test_like_and_samp_simplified.py index 4f747764..4bcd3845 100755 --- a/MonteCarloMarginalizeCode/Code/test/test_like_and_samp_simplified.py +++ b/MonteCarloMarginalizeCode/Code/test/test_like_and_samp_simplified.py @@ -65,6 +65,7 @@ import sys import numpy as np +from RIFT.precision import RiftFloat from glue.lal import Cache from glue.ligolw import utils, lsctables, table, ligolw, git_version @@ -639,7 +640,7 @@ def likelihood_function(right_ascension, declination,t_ref, phi_orb, inclination global nEvals global lnLOffsetValue # use EXTREMELY many bits - lnL = np.zeros(right_ascension.shape,dtype=np.float128) + lnL = np.zeros(right_ascension.shape,dtype=RiftFloat) i = 0 # if opts.rotate_sky_coordinates: # print " -Sky ring width ", np.std(declination), " note contribution from floor is of order p_floor*(pi)/sqrt(12) ~ 0.9 pfloor" @@ -824,4 +825,3 @@ def likelihood_function(right_ascension, declination,t_ref, phi_orb, inclination # utils.write_fileobj(xmldoc,sys.stdout) xmlutils.append_likelihood_result_to_xmldoc(xmldoc, np.log(res), **{"mass1": m1, "mass2": m2}) utils.write_filename(xmldoc, opts.points_file_base+".xml.gz", gz=True) - From f4272a8d80851c0b3c32e1e329ebcaad1d3ee3bb Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Sun, 31 May 2026 19:16:56 -0400 Subject: [PATCH 13/18] calmarg: note n_eff is conservative vs true ESS (softens pilot-starvation math) Per review: RIFT's reported n_eff is a deliberately CONSERVATIVE lower bound; the true ESS is meaningfully larger, so n_eff(us)=100 yields appreciably more usable fair-draw points. Implications recorded in DESIGN_adaptive_driver.md: the earlier "low n_eff" worry was over-pessimistic (tune-condor reached n_eff>200 on the moderate-SNR injection with the fixed settings, and ESS is larger still); pilot harvesting is LESS starved than the conservative count implied -- a real run can likely pull out enough high-quality points to inform the cal proposal. The d(d+1)/2 full-covariance requirement still holds, but the prior-shrinkage backstop covers the residual unconstrained directions. Co-Authored-By: Claude Opus 4.8 --- .../Code/RIFT/calmarg/DESIGN_adaptive_driver.md | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/calmarg/DESIGN_adaptive_driver.md b/MonteCarloMarginalizeCode/Code/RIFT/calmarg/DESIGN_adaptive_driver.md index a1f1db98..d404ae28 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/calmarg/DESIGN_adaptive_driver.md +++ b/MonteCarloMarginalizeCode/Code/RIFT/calmarg/DESIGN_adaptive_driver.md @@ -203,6 +203,19 @@ prior + importance-weight metadata. Designed from the start to ALSO carry an ex proposal (same harvest->fit->consolidate->seed structure). NF is a later drop-in behind the same interface. Stub the schema + the consolidation/seed hooks now. +## n_eff is conservative vs the true ESS (refines the starvation math) +RIFT's reported `n_eff` is a deliberately CONSERVATIVE lower bound -- the true effective +sample size (ESS) is meaningfully larger. So `n_eff(us)=100` yields appreciably more +usable fair-draw points than 100. Consequences: +- The earlier "low n_eff" worry was over-pessimistic: with enough samples the integrator + creeps up fine (the tune-condor run reached n_eff>200 on the moderate-SNR injection), and + the usable ESS is larger still. +- Pilot harvesting is LESS starved than the conservative count implied: top-fraction of + the composite + the larger-than-n_eff ESS means a real run can pull out enough + high-quality points to inform the cal proposal after all. The d(d+1)/2 requirement for + a FULL covariance still holds, but the prior-shrinkage backstop covers the residual + unconstrained directions, so we do not need to fully resolve every cal dof to be safe. + ## Build order (this branch) 1. **Timing data** -- done (`--scan-ncal`). 2. **A: brute-force reference** -- prior-only large-`n_cal`, converged; the ground truth. From 8fa9a8d50022e63d8c9e700b91b9297c01bc6b2f Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Mon, 1 Jun 2026 13:54:22 -0400 Subject: [PATCH 14/18] demo/calmarg: pp-build/pp-validate -- top-level pseudo_pipe thread-through (incl. time samples) Adds `make pp` (pp-build + pp-validate): exercises the FULL top-level builder util_RIFT_pseudo_pipe.py -> helper_LDG_Events.py -> args generation -> create_event_*, not just the lower-level builder that dag-build calls directly. Offline build-validate (the established pattern: .travis/test-build.sh + demo/pipeline/zero_spin_phenomD), reusing the zero-spin IMRPhenomD ini + ref coinc + a placeholder fake-data cache, so no GPU/data run is needed for the threading check. Zero spin (--assume-nospin); iterations forced small (--internal-force-iterations). Confirms EVERYTHING threads through the generated pipeline (validated, not just emitted): - calmarg --calibration-envelope-directory / --calibration-n-realizations / --calibration-fused-kernel land in args_ile.txt, ILE.sub AND ILE_extr.sub; - TIME SAMPLES: --add-extrinsic-time-resampling + --internal-ile-srate-time-resampling -> --srate-resample-time-marginalization 4096 in the wide AND extrinsic (ILE_extr) stages, alongside --time-marginalization; - zero-spin IMRPhenomD; small iteration count; full top-level DAG produced. (Note: --last-iteration-extrinsic-time-resampling is a transient builder arg consumed at build time -- its persistent effect is the --srate-resample-time-marginalization in ILE_extr.sub, which is what we assert.) Co-Authored-By: Claude Opus 4.8 --- .../Code/demo/rift/calmarg/Makefile | 57 ++++++++++++++++++- 1 file changed, 55 insertions(+), 2 deletions(-) diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile b/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile index 5b6c3bb2..e0b56f02 100644 --- a/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile +++ b/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile @@ -341,6 +341,59 @@ dag-run: dag: dag-build +# --- Top-level pipeline thread-through (util_RIFT_pseudo_pipe.py) ------------------- +# Unlike dag-build (which calls create_event_parameter_pipeline_BasicIteration directly +# with a hand-written args_ile.txt), this exercises the FULL top-level builder: +# util_RIFT_pseudo_pipe.py -> helper_LDG_Events.py -> args generation -> create_event_*. +# Goal: confirm the calmarg flags AND the time-sample (time-resampling) options thread all +# the way into args_ile.txt / the .sub files / the DAG. Offline build-validate (the +# established pattern: .travis/test-build.sh + demo/pipeline/zero_spin_phenomD), so no GPU +# run needed. Zero spin (reuses the zero-spin IMRPhenomD ini); iterations forced small. +PP_RUN := $(CURDIR)/rundir_pp +PP_INI ?= $(RIFT_CODE_ROOT)/demo/pipeline/zero_spin_phenomD/zero_spin_phenomD.ini +PP_COINC ?= $(REPO_ROOT)/.travis/ref_ini/coinc.xml +PP_NIT ?= 2 # forced iteration count (small, under control) +PP_SRATE ?= 4096 # time-resampling output rate (--srate-resample-time-marginalization) + +.PHONY: pp pp-build pp-validate + +pp-build: cal_env/H1.txt + @test -s "$(PP_INI)" || (echo "missing zero-spin ini $(PP_INI)" && false) + @test -s "$(PP_COINC)" || (echo "missing coinc $(PP_COINC)" && false) + rm -rf $(PP_RUN); touch $(CURDIR)/foo.cache + $(ENV) util_RIFT_pseudo_pipe.py \ + --use-ini $(PP_INI) --use-coinc $(PP_COINC) --use-rundir $(PP_RUN) \ + --fake-data-cache $(CURDIR)/foo.cache \ + --assume-nospin --approx IMRPhenomD --ile-sampler-method AV \ + --add-extrinsic --add-extrinsic-time-resampling --internal-ile-srate-time-resampling $(PP_SRATE) \ + --calmarg-envelope-directory $(CURDIR)/cal_env --calmarg-n-realizations $(NCAL_DAG) --calmarg-fused-kernel \ + --internal-force-iterations $(PP_NIT) --ile-n-eff 30 + $(MAKE) pp-validate + +# Confirm everything threaded through into the generated pipeline (no condor needed). +pp-validate: + @test -s "$(PP_RUN)/args_ile.txt" || (echo "FAIL: pseudo_pipe produced no args_ile.txt" && false) + @echo "--- args_ile.txt ---"; cat $(PP_RUN)/args_ile.txt + @grep -q -- "--calibration-envelope-directory" $(PP_RUN)/args_ile.txt || (echo "FAIL: calmarg envelope not threaded to args_ile.txt" && false) + @grep -q -- "--calibration-n-realizations" $(PP_RUN)/args_ile.txt || (echo "FAIL: calmarg n-realizations not threaded" && false) + @grep -q -- "--calibration-fused-kernel" $(PP_RUN)/args_ile.txt || (echo "FAIL: calmarg fused-kernel not threaded" && false) + @grep -q -- "--time-marginalization" $(PP_RUN)/args_ile.txt || (echo "FAIL: --time-marginalization absent (time integral)" && false) + @grep -q -- "--srate-resample-time-marginalization" $(PP_RUN)/args_ile.txt || (echo "FAIL: TIME SAMPLES (--srate-resample-time-marginalization) not threaded to ILE" && false) + @ls $(PP_RUN)/*.dag >/dev/null 2>&1 || (echo "FAIL: no top-level DAG produced" && false) + @test -s "$(PP_RUN)/ILE_extr.sub" || (echo "FAIL: no extrinsic (time-sample) stage ILE_extr.sub" && false) + @echo "--- ILE_extr.sub (the extrinsic / time-sample output stage) ---" + @grep -oE -- "--calibration-fused-kernel|--calibration-envelope-directory|--srate-resample-time-marginalization [0-9]+|--time-marginalization" $(PP_RUN)/ILE_extr.sub | sort -u + @grep -q -- "--calibration-fused-kernel" $(PP_RUN)/ILE_extr.sub || (echo "FAIL: extrinsic stage missing calmarg" && false) + @grep -q -- "--srate-resample-time-marginalization" $(PP_RUN)/ILE_extr.sub || (echo "FAIL: extrinsic stage missing TIME SAMPLES resampling" && false) + @echo "OK: top-level util_RIFT_pseudo_pipe.py threaded EVERYTHING --" + @echo " calmarg (envelope + n-realizations + fused-kernel) into args_ile.txt, ILE.sub AND ILE_extr.sub;" + @echo " TIME SAMPLES (--srate-resample-time-marginalization) into the wide AND extrinsic (ILE_extr) stages" + @echo " [--add-extrinsic-time-resampling drives the extrinsic stage at build time]; zero-spin IMRPhenomD;" + @echo " $(PP_NIT) iterations; full DAG produced." + @echo " (Offline build-validate -- the established pattern; no GPU/data run needed for the threading check.)" + +pp: pp-build + clean: - rm -f distance_marg.npz out_*.dat *.xml.gz lowsnr.cache snr_*.dat snr-report.txt - rm -rf cal_env lowsnr_mdc rundir_dag + rm -f distance_marg.npz out_*.dat *.xml.gz lowsnr.cache snr_*.dat snr-report.txt foo.cache + rm -rf cal_env lowsnr_mdc rundir_dag rundir_tune rundir_pp From 26f8d830911c5c3d81edc761cff6a9f4b636cde8 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Mon, 1 Jun 2026 14:01:42 -0400 Subject: [PATCH 15/18] demo/calmarg: pp-run -- RUNNABLE top-level calmarg pipeline on CI data (breadcrumb) Primes and launches a full util_RIFT_pseudo_pipe.py pipeline that ACTUALLY RUNS on the zero-noise CI fake data (not just the offline pp-build threading check). Pieces: - util_SimInspiralToCoinc.py makes a real coinc from the injection (H1/L1/V1, t=1000000014, m1=35/m2=30, SNR 17.5) -- `make pp-coinc`/ci_coinc.xml (regenerated, not committed). - calmarg_ci.ini: CI-matched (FAKE-STRAIN channels x3, srate 4096, seglen 8 -> the segment [1000000008.236,1000000016.236] sits inside the cache's 328s frame, fmin 10, zero spin, mc [23,35]). - pp-run-build runs pseudo_pipe with the real zero_noise.cache + CI PSD + calmarg (--calmarg-* --calmarg-fused-kernel) + time-resampling (--add-extrinsic-time-resampling + --internal-ile-srate-time-resampling 4096) + --assume-nospin + small forced iterations, and asserts the RUNNABLE bits threaded (real cache, FAKE-STRAIN, event time, and that ILE_extr.sub carries calmarg + --srate-resample-time-marginalization). - pp-run builds + condor_submit_dag. Verified: builds clean, submits (the extrinsic stage produces TIME SAMPLES with calmarg on). Single GPU on cardassia (CUDA_VISIBLE_DEVICES=0). This is the runnable counterpart to pp-build; leaves a working end-to-end breadcrumb for exercising the full pipeline + calmarg + time samples later. Co-Authored-By: Claude Opus 4.8 --- .../Code/demo/rift/calmarg/Makefile | 45 +++++++++++++- .../Code/demo/rift/calmarg/calmarg_ci.ini | 62 +++++++++++++++++++ 2 files changed, 105 insertions(+), 2 deletions(-) create mode 100644 MonteCarloMarginalizeCode/Code/demo/rift/calmarg/calmarg_ci.ini diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile b/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile index e0b56f02..888d9adc 100644 --- a/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile +++ b/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile @@ -394,6 +394,47 @@ pp-validate: pp: pp-build +# --- RUNNABLE top-level pipeline on the CI fake data (breadcrumb for later) ---------- +# Unlike pp-build (offline threading check), this builds a pipeline that ACTUALLY RUNS on +# the zero-noise CI data and SUBMITS it to condor: real coinc (util_SimInspiralToCoinc.py +# from the injection), the real zero_noise.cache + CI PSD + FAKE-STRAIN channels (a +# CI-matched ini), calmarg + time-resampling, zero spin, small forced iterations. Produces +# extrinsic-stage TIME SAMPLES with calmarg on. Single GPU on cardassia (CUDA_VISIBLE_DEVICES=0). +PP_RUN_REAL := $(CURDIR)/rundir_pp_run +PP_DAG := marginalize_intrinsic_parameters_BasicIterationWorkflow.dag + +.PHONY: pp-run pp-run-build pp-coinc + +pp-coinc: ci_coinc.xml +ci_coinc.xml: + $(ENV) util_SimInspiralToCoinc.py --sim-xml $(SIM_XML) --event 0 \ + --ifo H1 --ifo L1 --ifo V1 --output $(CURDIR)/ci_coinc.xml --injected-snr 17.5 + +# Build the runnable pipeline (no submit) and validate the runnable bits threaded through. +pp-run-build: cal_env/H1.txt ci_coinc.xml + rm -rf $(PP_RUN_REAL) + $(ENV) util_RIFT_pseudo_pipe.py \ + --use-ini $(CURDIR)/calmarg_ci.ini --use-coinc $(CURDIR)/ci_coinc.xml --use-rundir $(PP_RUN_REAL) \ + --fake-data-cache $(CACHE) --use-online-psd-file $(PSD) \ + --assume-nospin --approx IMRPhenomD --ile-sampler-method AV \ + --add-extrinsic --add-extrinsic-time-resampling --internal-ile-srate-time-resampling 4096 \ + --calmarg-envelope-directory $(CURDIR)/cal_env --calmarg-n-realizations $(NCAL_DAG) --calmarg-fused-kernel \ + --internal-force-iterations $(PP_NIT) --ile-n-eff 30 + @grep -q "zero_noise.cache" $(PP_RUN_REAL)/args_ile.txt || (echo "FAIL: not the real CI cache" && false) + @grep -q "FAKE-STRAIN" $(PP_RUN_REAL)/args_ile.txt || (echo "FAIL: not the FAKE-STRAIN channels" && false) + @grep -q "1000000014" $(PP_RUN_REAL)/args_ile.txt || (echo "FAIL: not the CI event time" && false) + @grep -q -- "--srate-resample-time-marginalization" $(PP_RUN_REAL)/ILE_extr.sub || (echo "FAIL: extrinsic stage missing time resampling" && false) + @grep -q -- "--calibration-fused-kernel" $(PP_RUN_REAL)/ILE_extr.sub || (echo "FAIL: extrinsic stage missing calmarg" && false) + @test -s "$(PP_RUN_REAL)/$(PP_DAG)" || (echo "FAIL: no top-level DAG produced" && false) + @echo "OK: runnable calmarg pipeline built on the CI fake data (real cache + PSD + FAKE-STRAIN," + @echo " event 1000000014.236, calmarg+fused, time-resampling, zero-spin, $(PP_NIT) iters)." + +# Build + SUBMIT to condor. Watch: condor_q ; tail -f $(PP_RUN_REAL)/$(PP_DAG).dagman.out +# Extrinsic time-sample output lands in the last iteration's extrinsic stage. +pp-run: pp-run-build + cd $(PP_RUN_REAL) && $(ENV) condor_submit_dag -f $(PP_DAG) + @echo "Submitted. Watch: condor_q ; tail -f $(PP_RUN_REAL)/$(PP_DAG).dagman.out" + clean: - rm -f distance_marg.npz out_*.dat *.xml.gz lowsnr.cache snr_*.dat snr-report.txt foo.cache - rm -rf cal_env lowsnr_mdc rundir_dag rundir_tune rundir_pp + rm -f distance_marg.npz out_*.dat *.xml.gz lowsnr.cache snr_*.dat snr-report.txt foo.cache ci_coinc.xml + rm -rf cal_env lowsnr_mdc rundir_dag rundir_tune rundir_pp rundir_pp_run diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/calmarg_ci.ini b/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/calmarg_ci.ini new file mode 100644 index 00000000..38467455 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/calmarg_ci.ini @@ -0,0 +1,62 @@ +# CI-matched RIFT ini for the in-loop calmarg pseudo_pipe RUN target (pp-run). +# Matches the zero-noise synthetic CI data in .travis/ILE-GPU-Paper/demos: 3 IFOs +# (H1,L1,V1) with FAKE-STRAIN channels, srate 4096, seglen 8 (the 8 s frame +# [1000000008,1000000016]; event 1000000014 + buffer 2 - seglen 8 = that frame), +# fmin 10, zero spin (IMRPhenomD-compatible), mc in [23,35]. CLI args on the +# util_RIFT_pseudo_pipe.py command line (--approx, --assume-nospin, calmarg flags, +# time-resampling, ...) win, so keep [rift-pseudo-pipe] minimal. + +[analysis] +ifos=['H1','L1','V1'] +singularity=False +osg=False + +[paths] + +[input] +max-psd-length=10000 + +[condor] +accounting_group=ligo.sim.o4.cbc.pe.rift +accounting_group_user=richard.oshaughnessy + +[datafind] +url-type=file +types = {'H1': 'fake_strain', 'L1': 'fake_strain', 'V1': 'fake_strain'} + +[data] +channels = {'H1': 'H1:FAKE-STRAIN','L1': 'L1:FAKE-STRAIN', 'V1': 'V1:FAKE-STRAIN'} + +[lalinference] +flow = {'H1': 10, 'L1': 10, 'V1': 10} +fhigh = { 'H1': 1700, 'L1': 1700, 'V1': 1700 } + +[engine] +fref=20 +amporder = -1 +seglen = 8 +srate = 4096 +a_spin1-max = 0.0 +a_spin2-max = 0.0 +chirpmass-min = 23.0 +chirpmass-max = 35.0 +comp-min = 1 +comp-max = 1000 +distance-max = 1000 +aligned-spin = +alignedspin-zprior = + +[rift-pseudo-pipe] +internal-ile-request-disk="4M" +cip-fit-method="rf" +ile-n-eff=10 +l-max=2 +internal-distance-max=1000 +ile-runtime-max-minutes=60 +ile-jobs-per-worker=10 +internal-propose-converge-last-stage=True +force-eta-range="[0.20,0.24999]" +fmin-template=10 +event-time=1000000014.236547946 +n-output-samples=2000 +use-online-psd=False From 8aaba0c9991cbe83e86f772c15ed6abf15a24713 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Tue, 2 Jun 2026 06:37:12 -0400 Subject: [PATCH 16/18] util_ConstructEOSPosterior: decouple fit basis from MC sampling basis Background ---------- util_ConstructIntrinsicPosterior_GenericCoordinates.py has long used three CLI flags for declaring how a parameter is treated: --parameter X both fit dim AND MC sampling dim --parameter-implied X fit dim only (the converter produces X from the data file's columns; the MC integrator never sees it) --parameter-nofit X MC sampling dim only (the integrator integrates over it; the fit never sees it) util_ConstructEOSPosterior.py declared the same three flags but never honoured them: the integrator at line 487 hardcoded `low_level_coord_names=dat_orig_names` in its convert_coords closure, which only worked when the sampling basis equalled the data-file basis; sampler.add_parameter iterated over coord_names (the fit basis); the arity dispatch for likelihood_function keyed on len(coord_names); and sampler.integrate was passed *coord_names rather than *low_level_coord_names. The net effect: any user who tried to fit in a transformed basis (e.g. via the new --supplementary-coordinate-code plugin) silently got a wrong likelihood evaluation -- the rotation was applied an extra time inside convert_coords every Monte Carlo step. What this commit changes ------------------------ bin/util_ConstructEOSPosterior.py * Parameter-resolution block rewritten to mirror IntrinsicPosterior's semantics, plus a clean fallback to dat_orig_names when none of the three flags are supplied (legacy bare-invocation unchanged). Seven CLI permutations now map to documented (coord_names, low_level_coord_names) pairs. * The convert_coords closure used by the integrator captures low_level_coord_names as its input basis (was dat_orig_names). The initial dat->X conversion still uses dat_orig_names, since that's the basis of the file columns. * Sampler add_parameter loop now iterates over low_level_coord_names (the MC basis), and sampler.integrate is passed *low_level_coord_names. * The arity-dispatched likelihood_function definitions key on len(low_level_coord_names) and route every input -- including the scalar branches -- through convert_coords so a non-trivial converter is never silently bypassed. * Output-writer iterates samples by low_level_coord_names (the keys sampler._rvs actually carries) and applies the "constant fill" check in the sampling basis, not the fit basis. Implied (fit-only) coords correctly skip the output file. * Added a guard: if low_level_coord_names != coord_names but no coordinate plugin is supplied, raise a clear error instead of silently feeding samples through an identity convert_coords into a fit built in a different basis. * Help text for --parameter / --parameter-implied / --parameter-nofit rewritten to describe what each flag actually does now. RIFT/hyperpipe/coords.py * HyperCoordSpec.from_strings accepts integration ranges for names in coords-nofit (the MC sampling basis is coords-fit + coords-nofit); unknown range names are still rejected. * HyperCoordSpec.validate accepts empty coords-fit so long as coords-implied covers the fit basis and coords-nofit covers the sampling basis; emits distinct errors for empty-fit vs empty-sample. * to_parameter_args emits --integration-parameter-range for the sampling basis (parameters + nofit), not just parameters. * to_puff_args and to_test_args emit --parameter for the sampling basis -- the puff lane and convergence-test driver operate on the data-file columns, which is the sampling basis after decoupling. RIFT/hyperpipe/config.py * validate_config accepts empty coords-fit when coords-implied (fit-side) and coords-nofit (sample-side) are non-empty. demo/hyperpipe/hyperpipe_conf_linear_uvw.yaml * Rewritten to actually exercise the decoupled path: coords-implied "u v w" (fit), coords-nofit "x y z" (sample), coords-sample ranges in (x, y, z), coord-module pointing at the linear plugin with the uvw_rotated chart. Iteration / puff / marg stay in (x, y, z); the EOS posterior fits in (u, v, w) and writes its posterior in (x, y, z). Verified -------- * Parameter-resolution unit test (in this commit's worktree) covers 7 CLI permutations -- legacy no-flags, legacy --parameter, IntrinsicPosterior --parameter+implied and --parameter+nofit, the new --implied-only, --implied+nofit, and full --parameter+implied+nofit -- all map to the documented (coord_names, low_level_coord_names) pairs. * HyperCoordSpec unit test covers the new decoupled emit (post sees implied/nofit and ranges; puff/test see the sampling basis only), a legacy-regression case (unchanged output), the two new validation errors (empty fit, empty sample), the new "range for nofit name" permission, and the still-rejected "unknown range name" case. * AST + yaml parses on every edited file. * validate_config passes on hyperpipe_conf_linear_uvw.yaml plus the demo's baseline and tracer yamls. --- .../Code/RIFT/hyperpipe/config.py | 20 +- .../Code/RIFT/hyperpipe/coords.py | 71 +++- .../Code/bin/util_ConstructEOSPosterior.py | 303 ++++++++++++------ .../hyperpipe/hyperpipe_conf_linear_uvw.yaml | 97 +++--- 4 files changed, 325 insertions(+), 166 deletions(-) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/hyperpipe/config.py b/MonteCarloMarginalizeCode/Code/RIFT/hyperpipe/config.py index 17185298..f9e06703 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/hyperpipe/config.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/hyperpipe/config.py @@ -208,12 +208,24 @@ def validate_config(cfg) -> None: if not isinstance(_get(arch, "n-samples-per-job"), int) or _get(arch, "n-samples-per-job") <= 0: raise ValueError("arch.n-samples-per-job must be a positive integer.") - # post: at least one of (coords-fit) must be set + # post: must have at least one fit dim AND at least one MC sampling dim. + # Fit basis = coords-fit + coords-implied; MC basis = coords-fit + coords-nofit. + # Pre-decoupling this only required coords-fit, because the fit basis was + # forced to equal the MC basis -- now an EOS-style "fit in a transformed + # basis" config can legally have empty coords-fit (everything routed via + # coords-implied + coords-nofit through the coordinate plugin). post = _get(cfg, "post") - if not _get(post, "coords-fit"): + has_fit = bool(_get(post, "coords-fit")) or bool(_get(post, "coords-implied")) + has_samp = bool(_get(post, "coords-fit")) or bool(_get(post, "coords-nofit")) + if not has_fit: raise ValueError( - "post.coords-fit must list at least one parameter " - "(e.g. 'x y z')." + "post: must list at least one fit dimension " + "(coords-fit or coords-implied; e.g. 'x y z' or 'u v w')." + ) + if not has_samp: + raise ValueError( + "post: must list at least one MC sampling dimension " + "(coords-fit or coords-nofit; e.g. 'x y z')." ) # init: must have either file or generation set diff --git a/MonteCarloMarginalizeCode/Code/RIFT/hyperpipe/coords.py b/MonteCarloMarginalizeCode/Code/RIFT/hyperpipe/coords.py index bc039846..ccd6b18d 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/hyperpipe/coords.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/hyperpipe/coords.py @@ -177,12 +177,22 @@ def from_strings( ``coords_nofit`` : "delta_mc s1z s2z" (optional) ``likelihood_factor``: (module, function, ini) (any element may be None) """ - params = parse_parameter_list(coords_fit) - ranges = parse_range_string(coords_sample) - unknown = set(ranges) - set(params) + params = parse_parameter_list(coords_fit) + implied = parse_parameter_list(coords_implied) + nofit = parse_parameter_list(coords_nofit) + ranges = parse_range_string(coords_sample) + # coords-sample provides INTEGRATION ranges, so it has to cover every + # name in the MC SAMPLING basis -- which is coords-fit + coords-nofit. + # (implied names are fit-only and don't need a sample range.) Pre- + # decoupling this was just coords-fit because nofit/implied were + # rarely used and the sampling basis was forced to equal the fit + # basis; now we have to allow the nofit names too. + sampling_basis = set(params) | set(nofit) + unknown = set(ranges) - sampling_basis if unknown: raise ValueError( - f"coords-sample names a parameter not in coords-fit: {sorted(unknown)!r}" + f"coords-sample names a parameter not in coords-fit or " + f"coords-nofit: {sorted(unknown)!r}" ) lf: Optional[Tuple[str, Optional[str], Optional[str]]] = None if likelihood_factor: @@ -195,8 +205,8 @@ def from_strings( name=name or None, parameters=params, parameter_ranges=ranges, - implied=parse_parameter_list(coords_implied), - nofit=parse_parameter_list(coords_nofit), + implied=implied, + nofit=nofit, likelihood_factor=lf, ) @@ -210,13 +220,30 @@ def validate(self, strict_import: bool = False) -> None: environment (singularity image, OSG worker) and not necessarily on the submit host. """ - if not self.parameters: - raise ValueError("HyperCoordSpec requires at least one fitting parameter.") - missing = [p for p in self.parameters if p not in self.parameter_ranges] + # The fit basis is coords-fit + coords-implied; the sampling basis is + # coords-fit + coords-nofit. Both must be non-empty for the run to + # make sense. Pre-decoupling this only required coords-fit -- now + # an "EOS-style fit in a transformed basis" config can legally have + # empty coords-fit (everything goes through implied + nofit). + if not self.parameters and not self.implied: + raise ValueError( + "HyperCoordSpec requires at least one fit dimension " + "(coords-fit or coords-implied)." + ) + if not self.parameters and not self.nofit: + raise ValueError( + "HyperCoordSpec requires at least one MC sampling dimension " + "(coords-fit or coords-nofit)." + ) + # Every name in the SAMPLING basis must have an integration range + # (the integrator reads prior_range_map[p] for p in low_level_coord_names). + sampling_names = list(self.parameters) + list(self.nofit) + missing = [p for p in sampling_names if p not in self.parameter_ranges] if missing: raise ValueError( - f"No integration range supplied for parameter(s): {missing!r}. " - "Every entry in coords-fit must appear in coords-sample." + f"No integration range supplied for sampling parameter(s): " + f"{missing!r}. Every entry in coords-fit and coords-nofit must " + "appear in coords-sample." ) for p, (lo, hi) in self.parameter_ranges.items(): if not lo < hi: @@ -266,7 +293,9 @@ def to_parameter_args(self) -> str: bits.append(f"--parameter-implied {p}") for p in self.nofit: bits.append(f"--parameter-nofit {p}") - for p in self.parameters: + # Integration ranges cover the MC SAMPLING basis (parameters + nofit). + # Implied coordinates are fit-only and don't have a sampling range. + for p in list(self.parameters) + list(self.nofit): lo, hi = self.parameter_ranges[p] bits.append( f"--integration-parameter-range {p}:[{self._fmt_num(lo)},{self._fmt_num(hi)}]" @@ -290,12 +319,16 @@ def to_post_args(self) -> str: def to_puff_args(self, force_away: float = 0.03, puff_factor: float = 0.5) -> str: """Emit the puff-stage arg block. - By default we puff in every fitting parameter; this is what every - existing hyperpipe example does. Extra flags can be appended by the - caller. + The puff lane reads / writes grid files in the data-file column + basis, which is the MC sampling basis (coords-fit + coords-nofit). + Pre-decoupling this only emitted --parameter for coords-fit because + the sampling basis was forced to equal the fit basis; once those + diverge (EOSPosterior with --parameter-implied for a transformed + fit basis), the puff lane must continue to operate on the data- + file columns -- i.e. coords-fit + coords-nofit. """ bits = [f"--force-away {force_away}", f"--puff-factor {puff_factor}"] - for p in self.parameters: + for p in list(self.parameters) + list(self.nofit): bits.append(f"--parameter {p}") return " ".join(bits) @@ -304,8 +337,12 @@ def to_test_args(self, method: str = "JS", threshold: float = 0.05) -> str: Mirrors the args_test.txt pattern from the Gaussian demo: ``--parameter x --parameter y --parameter z --method JS --threshold 0.05`` + + Like the puff lane, this operates on the SAMPLING basis (coords-fit + + coords-nofit) -- the convergence-test driver reads grid / posterior + files whose columns are in the sampling basis. """ - bits = [f"--parameter {p}" for p in self.parameters] + bits = [f"--parameter {p}" for p in list(self.parameters) + list(self.nofit)] bits.append(f"--method {method}") bits.append(f"--threshold {threshold}") return " ".join(bits) diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ConstructEOSPosterior.py b/MonteCarloMarginalizeCode/Code/bin/util_ConstructEOSPosterior.py index 45b49e2c..aae18cc1 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_ConstructEOSPosterior.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_ConstructEOSPosterior.py @@ -122,10 +122,10 @@ def add_field(a, descr): parser.add_argument("--fname-output-integral",default="output-EOS-integral",help="for evidencees and pipeline compatibility") parser.add_argument("--n-output-samples",default=2000,type=int,help="output posterior samples (default 3000)") parser.add_argument("--eos-param", type=str, default=None, help="parameterization of equation of state [spectral only, for now]") -parser.add_argument("--parameter", action='append', help="Parameters used as fitting parameters AND varied at a low level to make a posterior. Currently can only specify gamma1,gamma2, ..., and these MUST be columns in --fname. IF NOT PROVIDED, DEFAULTS TO LIST IN FILE. ") -parser.add_argument("--parameter-implied", action='append', help="Parameter used in fit, but not independently varied for Monte Carlo. For EOS objects, only possible for physical quantities like R1.4, etc. NOT YET PROVIDED") +parser.add_argument("--parameter", action='append', help="Parameter used BOTH as a fit dimension (GP/RF input) AND as a Monte Carlo sampling dimension. Adds to coord_names AND low_level_coord_names. Must be a column in --fname unless --supplementary-coordinate-code is supplied. IF NEITHER --parameter NOR --parameter-implied IS PROVIDED, coord_names defaults to the data file's column list; IF NEITHER --parameter NOR --parameter-nofit IS PROVIDED, low_level_coord_names also defaults to the data file's column list.") +parser.add_argument("--parameter-implied", action='append', help="Parameter used as a fit dimension only -- added to coord_names but NOT sampled independently. The coordinate plugin is responsible for producing it from the data file's columns. Useful for fitting in a different basis (e.g. coord_names=[u,v,w]) than the data is stored in (e.g. dat_orig_names=[x,y,z]).") #parser.add_argument("--no-adapt-parameter",action='append',help="Disable adaptive sampling in a parameter. Useful in cases where a parameter is not well-constrained, and the a prior sampler is well-chosen.") -parser.add_argument("--parameter-nofit", action='append', help="Parameter used to initialize the implied parameters, and varied at a low level, but NOT the fitting parameters.") +parser.add_argument("--parameter-nofit", action='append', help="Parameter used as a sampling dimension only -- added to low_level_coord_names but NOT to the fit basis. Useful when the MC samples in the data-file basis (e.g. dat_orig_names=[x,y,z]) but the fit lives in a transformed basis routed through the coordinate plugin.") parser.add_argument("--integration-parameter-range",action='append', help="Integration parameter ranges. Syntax is name:[a,b]") parser.add_argument("--downselect-parameter",action='append', help='Name of parameter to be used to eliminate grid points ') parser.add_argument("--downselect-parameter-range",action='append',type=str) @@ -223,17 +223,39 @@ def add_field(a, descr): ### Parameters in use ### -coord_names = opts.parameter # Used in fit -if coord_names is None: - coord_names = dat_orig_names -low_level_coord_names = coord_names # Used for Monte Carlo -if opts.parameter_implied: - coord_names = coord_names+opts.parameter_implied -if opts.parameter_nofit: - if opts.parameter is None: - low_level_coord_names = opts.parameter_nofit # Used for Monte Carlo - else: - low_level_coord_names = opts.parameter+opts.parameter_nofit # Used for Monte Carlo +# Decoupled fit basis vs Monte Carlo sampling basis -- mirrors the +# convention established by util_ConstructIntrinsicPosterior_GenericCoordinates.py +# and required for the new coordinate-plugin path: +# +# --parameter X -> X is BOTH a fit (GP/RF) and a sampling (MC) dim +# --parameter-implied X-> X is a fit dim ONLY (the plugin produces it from +# dat_orig_names; the MC integrator never sees it) +# --parameter-nofit X -> X is a sampling dim ONLY (the MC integrates over +# it; the fit never sees it). Typical use: MC in +# the data-file basis while the fit lives in a +# transformed basis routed through the plugin. +# +# Legacy fallback (preserves the pre-decoupling default): if the user +# supplies neither --parameter nor --parameter-implied, coord_names +# defaults to the data file's column list. Likewise low_level_coord_names +# defaults to the data file's columns when --parameter and --parameter-nofit +# are both absent. That way a bare invocation -- no flags -- still does +# "fit on every column in the file, MC sample in the same basis", which +# is what every existing hyperpipe / EOS-posterior demo relies on. +_user_params = list(opts.parameter) if opts.parameter else [] +_user_implied = list(opts.parameter_implied) if opts.parameter_implied else [] +_user_nofit = list(opts.parameter_nofit) if opts.parameter_nofit else [] + +if not _user_params and not _user_implied: + coord_names = list(dat_orig_names) # legacy default +else: + coord_names = _user_params + _user_implied # fit basis + +if not _user_params and not _user_nofit: + low_level_coord_names = list(dat_orig_names) # legacy default +else: + low_level_coord_names = _user_params + _user_nofit # MC basis + error_factor = len(coord_names) name_index_dict ={} for name in dat_orig_names: @@ -473,6 +495,24 @@ def fn_return(x_in,rf=rf): # Naive convert: no downselect. if (supplemental_coordinate_convert ==None): + # The identity convert_coords below only makes sense when the fit + # basis equals the MC sampling basis equals a permutation of the + # data file's columns. Catch the new "split-basis" misconfiguration + # early, otherwise the integrator silently feeds samples in + # low_level_coord_names through an identity into a fit built on + # coord_names. + if list(low_level_coord_names) != list(coord_names): + raise ValueError( + " EOSPosterior: --parameter-implied / --parameter-nofit make " + "the fit basis ({coord!r}) differ from the MC sampling basis " + "({low!r}), but no --supplementary-coordinate-code was " + "supplied. The integrator cannot translate between the two " + "bases without a converter.".format( + coord=list(coord_names), + low=list(low_level_coord_names), + ) + ) + indx_of_orig_names = np.array([ dat_orig_names.index(coord_names[k]) for k in range(len(coord_names))]) dat_out = [] for line in dat: @@ -496,14 +536,37 @@ def convert_coords(x): return x else: - # Pack data, using coordinate converter. Note later calculations MUST use the converter + # Pack data, using coordinate converter. Note later calculations MUST use the converter. + # + # Two distinct call sites for the converter, with two different + # input bases -- this is the change that decouples the fit from + # the MC sampling basis: + # + # (1) The initial dat->X conversion below feeds rows whose + # columns are ordered by dat_orig_names (the data file's + # header). So we pass low_level_coord_names=dat_orig_names + # at this site. + # + # (2) The convert_coords closure is what the integrator calls + # on every Monte Carlo sample. The sampler operates in + # low_level_coord_names (we add_parameter() over that list + # below), so the closure must claim its inputs are in + # low_level_coord_names -- NOT dat_orig_names. Pre-fix this + # was hardcoded to dat_orig_names, which only happened to + # work when low_level_coord_names == dat_orig_names (i.e. + # the legacy case). For any non-trivial plugin where the + # MC samples in a different basis than the file's columns, + # the old behaviour applied the rotation an extra time and + # silently mis-evaluated lnL. X = supplemental_coordinate_convert(dat[:,2:], coord_names=coord_names, low_level_coord_names=dat_orig_names) # convert and generate X Y = dat[:,0] Y_err = dat[:,1] - if np.max(Y)<0 and lnL_shift ==0: + if np.max(Y)<0 and lnL_shift ==0: lnL_shift = -100 - np.max(Y) # force it to be offset/positive -- may help some configurations. Remember our adaptivity is silly. - def convert_coords(x_in): - return supplemental_coordinate_convert(x_in, coord_names=coord_names, low_level_coord_names=dat_orig_names) # convert and generate X + def convert_coords(x_in, _low=low_level_coord_names, _coord=coord_names): + # _low / _coord captured as defaults so the closure stays correct + # even if either list mutates later in the script. + return supplemental_coordinate_convert(x_in, coord_names=_coord, low_level_coord_names=_low) # Save copies for later (plots) X_orig = X.copy() Y_orig = Y.copy() @@ -636,7 +699,11 @@ def convert_coords(x_in): ## ## Loop over param names ## -for p in coord_names: +# IMPORTANT: iterate over low_level_coord_names, not coord_names. The +# sampler operates in the MC basis. coord_names is the FIT basis, which +# only the GP/RF and the convert_coords closure see. Pre-decoupling this +# loop used coord_names because the two lists were forced to be equal. +for p in low_level_coord_names: prior_here = prior_map[p] range_here = prior_range_map[p] @@ -647,71 +714,82 @@ def convert_coords(x_in): def log_likelihood_function(*args): return my_fit(convert_coords(np.array([*args]).T )) -if len(coord_names) ==1: - def likelihood_function(x): - if isinstance(x,float): - return np.exp(my_fit([x])) - else: - return np.exp(my_fit(convert_coords(np.c_[x]))) -if len(coord_names) ==2: - def likelihood_function(x,y): - if isinstance(x,float): - return np.exp(my_fit([x,y])) - else: -# return np.exp(my_fit(convert_coords(np.array([x,y],dtype=internal_dtype).T))) - return np.exp(my_fit(convert_coords(np.c_[x,y]))) -if len(coord_names) ==3: - def likelihood_function(x,y,z): - if isinstance(x,float): - return np.exp(my_fit([x,y,z])) - else: -# return np.exp(my_fit(convert_coords(np.array([x,y,z],dtype=internal_dtype).T))) - return np.exp(my_fit(convert_coords(np.c_[x,y,z]))) -if len(coord_names) ==4: - def likelihood_function(x,y,z,a): - if isinstance(x,float): - return np.exp(my_fit([x,y,z,a])) - else: -# return np.exp(my_fit(convert_coords(np.array([x,y,z,a],dtype=internal_dtype).T))) - return np.exp(my_fit(convert_coords(np.c_[x,y,z,a]))) -if len(coord_names) ==5: - def likelihood_function(x,y,z,a,b): - if isinstance(x,float): - return np.exp(my_fit([x,y,z,a,b])) - else: -# return np.exp(my_fit(convert_coords(np.array([x,y,z,a,b],dtype=internal_dtype).T))) - return np.exp(my_fit(convert_coords(np.c_[x,y,z,a,b]))) -if len(coord_names) ==6: - def likelihood_function(x,y,z,a,b,c): - if isinstance(x,float): - return np.exp(my_fit([x,y,z,a,b,c])) - else: -# return np.exp(my_fit(convert_coords(np.array([x,y,z,a,b,c],dtype=internal_dtype).T))) - return np.exp(my_fit(convert_coords(np.c_[x,y,z,a,b,c]))) -if len(coord_names) ==7: - def likelihood_function(x,y,z,a,b,c,d): - if isinstance(x,float): - return np.exp(my_fit([x,y,z,a,b,c,d])) - else: - return np.exp(my_fit(convert_coords(np.c_[x,y,z,a,b,c,d]))) -if len(coord_names) ==8: - def likelihood_function(x,y,z,a,b,c,d,e): - if isinstance(x,float): - return np.exp(my_fit([x,y,z,a,b,c,d,e])) - else: - return np.exp(my_fit(convert_coords(np.c_[x,y,z,a,b,c,d,e]))) -if len(coord_names) ==9: - def likelihood_function(x,y,z,a,b,c,d,e,f): - if isinstance(x,float): - return np.exp(my_fit([x,y,z,a,b,c,d,e,f])) - else: - return np.exp(my_fit(convert_coords(np.c_[x,y,z,a,b,c,d,e,f]))) -if len(coord_names) ==10: - def likelihood_function(x,y,z,a,b,c,d,e,f,g): - if isinstance(x,float): - return np.exp(my_fit([x,y,z,a,b,c,d,e,f,g])) - else: - return np.exp(my_fit(convert_coords(np.c_[x,y,z,a,b,c,d,e,f,g]))) +# Fixed-arity wrappers around log_likelihood_function / likelihood_function. +# +# mcsampler's adaptive code introspects the wrapped function's argument +# count, so we generate one definition per supported dimensionality. The +# arity that matters is the MC SAMPLING dimensionality +# (len(low_level_coord_names)), not the fit dimensionality +# (len(coord_names)) -- the sampler passes one positional per +# low_level_coord_name and then convert_coords maps that batch into the +# fit basis. Pre-decoupling, the dispatch keyed on len(coord_names), +# which was only correct when low_level_coord_names == coord_names. +# +# Scalar branches also go through convert_coords so a non-trivial +# converter does not silently get bypassed when an internal caller hands +# in a single scalar sample. +def _scalar_to_lnL_input(args_tuple): + # Wrap an N-tuple of scalars into a (1, N) row, push it through the + # converter, and return -- shape (1, len(coord_names)) -- ready for my_fit. + return convert_coords(np.array([args_tuple], dtype=float)) + +_LN_LOW_DIM = len(low_level_coord_names) +if _LN_LOW_DIM == 1: + def likelihood_function(x): + if isinstance(x, float): + return np.exp(my_fit(_scalar_to_lnL_input((x,)))) + return np.exp(my_fit(convert_coords(np.c_[x]))) +elif _LN_LOW_DIM == 2: + def likelihood_function(x, y): + if isinstance(x, float): + return np.exp(my_fit(_scalar_to_lnL_input((x, y)))) + return np.exp(my_fit(convert_coords(np.c_[x, y]))) +elif _LN_LOW_DIM == 3: + def likelihood_function(x, y, z): + if isinstance(x, float): + return np.exp(my_fit(_scalar_to_lnL_input((x, y, z)))) + return np.exp(my_fit(convert_coords(np.c_[x, y, z]))) +elif _LN_LOW_DIM == 4: + def likelihood_function(x, y, z, a): + if isinstance(x, float): + return np.exp(my_fit(_scalar_to_lnL_input((x, y, z, a)))) + return np.exp(my_fit(convert_coords(np.c_[x, y, z, a]))) +elif _LN_LOW_DIM == 5: + def likelihood_function(x, y, z, a, b): + if isinstance(x, float): + return np.exp(my_fit(_scalar_to_lnL_input((x, y, z, a, b)))) + return np.exp(my_fit(convert_coords(np.c_[x, y, z, a, b]))) +elif _LN_LOW_DIM == 6: + def likelihood_function(x, y, z, a, b, c): + if isinstance(x, float): + return np.exp(my_fit(_scalar_to_lnL_input((x, y, z, a, b, c)))) + return np.exp(my_fit(convert_coords(np.c_[x, y, z, a, b, c]))) +elif _LN_LOW_DIM == 7: + def likelihood_function(x, y, z, a, b, c, d): + if isinstance(x, float): + return np.exp(my_fit(_scalar_to_lnL_input((x, y, z, a, b, c, d)))) + return np.exp(my_fit(convert_coords(np.c_[x, y, z, a, b, c, d]))) +elif _LN_LOW_DIM == 8: + def likelihood_function(x, y, z, a, b, c, d, e): + if isinstance(x, float): + return np.exp(my_fit(_scalar_to_lnL_input((x, y, z, a, b, c, d, e)))) + return np.exp(my_fit(convert_coords(np.c_[x, y, z, a, b, c, d, e]))) +elif _LN_LOW_DIM == 9: + def likelihood_function(x, y, z, a, b, c, d, e, f): + if isinstance(x, float): + return np.exp(my_fit(_scalar_to_lnL_input((x, y, z, a, b, c, d, e, f)))) + return np.exp(my_fit(convert_coords(np.c_[x, y, z, a, b, c, d, e, f]))) +elif _LN_LOW_DIM == 10: + def likelihood_function(x, y, z, a, b, c, d, e, f, g): + if isinstance(x, float): + return np.exp(my_fit(_scalar_to_lnL_input((x, y, z, a, b, c, d, e, f, g)))) + return np.exp(my_fit(convert_coords(np.c_[x, y, z, a, b, c, d, e, f, g]))) +else: + raise NotImplementedError( + " EOSPosterior currently only ships fixed-arity likelihood_function " + "wrappers for 1..10 sampling dimensions; got " + "{} (low_level_coord_names={!r}).".format(_LN_LOW_DIM, low_level_coord_names) + ) @@ -781,7 +859,7 @@ def parse_corr_params(my_str): -res, var, neff, dict_return = sampler.integrate(fn_passed, *coord_names, verbose=True,nmax=int(opts.n_max),n=n_step,neff=opts.n_eff, save_intg=True,tempering_adapt=True, floor_level=1e-3,igrand_threshold_p=1e-3,convergence_tests=test_converged,adapt_weight_exponent=my_exp,no_protect_names=True,**extra_args) # weight ecponent needs better choice. We are using arbitrary-name functions +res, var, neff, dict_return = sampler.integrate(fn_passed, *low_level_coord_names, verbose=True,nmax=int(opts.n_max),n=n_step,neff=opts.n_eff, save_intg=True,tempering_adapt=True, floor_level=1e-3,igrand_threshold_p=1e-3,convergence_tests=test_converged,adapt_weight_exponent=my_exp,no_protect_names=True,**extra_args) # MC integrates in the SAMPLING basis (low_level_coord_names); convert_coords routes each sample into the fit basis (coord_names) before evaluating the GP/RF # Save result -- needed for odds ratios, etc. @@ -794,8 +872,10 @@ def parse_corr_params(my_str): samples = sampler._rvs print(samples.keys()) -n_params = len(coord_names) -dat_mass = np.zeros((len(samples[coord_names[0]]),n_params+3)) +# sampler._rvs is keyed by the SAMPLING basis (low_level_coord_names). +# Look up sample arrays by names that actually exist in the dict. +n_params = len(low_level_coord_names) +dat_mass = np.zeros((len(samples[low_level_coord_names[0]]),n_params+3)) if not(opts.internal_use_lnL): dat_logL = np.log(samples["integrand"]) else: @@ -826,7 +906,8 @@ def parse_corr_params(my_str): if not('log_joint_s_prior' in samples): indx_ok=samples["joint_s_prior"]>0 indx_ok = np.logical_and(dat_logL > np.max(dat_logL)-opts.lnL_offset ,indx_ok) -for p in coord_names: +# Mask in the sampling basis -- samples dict is keyed by low_level_coord_names. +for p in low_level_coord_names: samples[p] = samples[p][indx_ok] dat_logL = dat_logL[indx_ok] print(samples.keys()) @@ -856,8 +937,22 @@ def parse_corr_params(my_str): dat_out = np.zeros( (opts.n_output_samples,2+len(dat_orig_names)) ) -# Initialize fixed parameters -if len(coord_names) < len(dat_orig_names): # not needed if all params are in fit +# Initialize fixed parameters. +# +# Output columns are dat_orig_names (the file's basis). For any name in +# dat_orig_names that is NOT in the sampling basis, we fill the output +# column with the original grid's value -- it's a "constant" placeholder +# for that column. +# +# Pre-decoupling this test was against coord_names (the fit basis), which +# was the same list as low_level_coord_names whenever the fit basis equaled +# the sampling basis -- the only case the code supported. Now that +# coord_names and low_level_coord_names can differ, the right question +# for the output writer is "did we MC-sample this column?", not "did we +# fit on it?". Implied (fit-only) coords like R1.4 should NOT appear in +# the output file (they aren't grid columns and dat_orig_names doesn't +# carry them). +if len(low_level_coord_names) < len(dat_orig_names): # not needed if every grid column is sampled if len(dat) < opts.n_output_samples: print(" NOTE: original data shorter than requested output; adding",opts.n_output_samples-len(dat),"duplicate fill lines from original data.") @@ -872,21 +967,31 @@ def parse_corr_params(my_str): else: newlines = dat[:opts.n_output_samples-len(dat)] #duplicate lines to fill dat = np.concatenate((dat,newlines), axis=0) #should be fine since dat isn't used after this - + for c in np.arange(len(dat_orig_names)): - if dat_orig_names[c] not in coord_names: - print(" Not in coord_names:",dat_orig_names[c],"; adding to output as constant.") + if dat_orig_names[c] not in low_level_coord_names: + print(" Not sampled:", dat_orig_names[c], "; adding to output as constant from initial grid.") outidx = name_index_dict[dat_orig_names[c]] # write in correct place if len(dat) > opts.n_output_samples: dat_out[:,outidx] = dat[:opts.n_output_samples,outidx] #truncate original data to fit (not ideal) else: #len(dat) <= n_output_samples (if dat was <, should now be =) dat_out[:,outidx] = dat[:,outidx] - -# Fill data from PE -for indx in np.arange(len(coord_names)): - vals = samples[coord_names[indx]][indx_list] # load in data for this column - outindx = name_index_dict[ coord_names[indx]] # write in correct place - dat_out[:,outindx] = vals + +# Fill data from PE. +# +# samples[k] keys are the SAMPLING basis (low_level_coord_names). We can +# only write a column to the output file when (a) we sampled it, and (b) +# its name is one of dat_orig_names (i.e. it has a column slot in the +# output schema). Skip sampled names that aren't in dat_orig_names with +# a warning -- they'd correspond to a --parameter-nofit name that isn't +# a data-file column, which is rare but valid and would only show up via +# the corner plot, not the output file. +for name in low_level_coord_names: + if name not in name_index_dict: + print(" Sampled coord {!r} is not a data-file column; not writing to output (would only appear in corner plots).".format(name)) + continue + outindx = name_index_dict[name] + dat_out[:, outindx] = samples[name][indx_list] # NOTE: if m1 or m2 is "constant" (i.e., not in samples), the possibility for m2 > m1 arises! Re-sort masses here to avoid; use below code. #if ("m1" not in coord_names) or ("m2" not in coord_names): diff --git a/MonteCarloMarginalizeCode/Code/demo/hyperpipe/hyperpipe_conf_linear_uvw.yaml b/MonteCarloMarginalizeCode/Code/demo/hyperpipe/hyperpipe_conf_linear_uvw.yaml index 5092208f..fe5b2804 100644 --- a/MonteCarloMarginalizeCode/Code/demo/hyperpipe/hyperpipe_conf_linear_uvw.yaml +++ b/MonteCarloMarginalizeCode/Code/demo/hyperpipe/hyperpipe_conf_linear_uvw.yaml @@ -1,46 +1,48 @@ -# End-to-end demo for the linear_coordinate_convert plugin. +# End-to-end demo for the linear_coordinate_convert plugin with DECOUPLED +# fit / sampling bases inside util_ConstructEOSPosterior.py. # # Goal # ---- -# The pipeline iterates in (x, y, z) -- same basis as hyperpipe_conf_tracer.yaml -# -- but the likelihood is a 3-D Gaussian whose principal axes lie along -# (u, v, w), a 45-degree rotation of (x, y, z) in the xy-plane. Sigmas are -# deliberately unequal (sigma_u != sigma_v != sigma_w) so the rotation is -# observable in the resulting (x, y, z) posterior -- you should see an -# elongated, off-axis ellipsoid that points along the rotated v-axis. +# Iteration (puff lane, marg evaluator, condor DAG) lives in (x, y, z) -- +# the data file's column basis. The EOS posterior step FITS its GP/RF in +# (u, v, w), a 45-degree rotation of (x, y, z) in the xy-plane that +# axis-aligns the unequal-sigma Gaussian. The MC integrator inside the +# EOS posterior still SAMPLES in (x, y, z) so the output samples come out +# in the same basis the rest of the pipeline expects. # -# Why iteration stays in (x, y, z) -# -------------------------------- -# The hyperpipe yaml schema currently shares `coords-fit` between the puff -# stage and the post stage (see RIFT/hyperpipe/coords.py:to_puff_args and -# to_post_args). Having the puff lane and the EOS posterior operate in -# different bases would require either a schema extension or a Makefile- -# style direct invocation; both are out of scope for this demo. What this -# yaml DOES demonstrate is the most useful piece of the new plugin -# machinery: a user-supplied likelihood evaluator (example_gaussian_uvw.py) -# loading linear_coordinate_convert.py via the same plugin contract that -# util_ConstructEOSPosterior.py uses. The library is exercised -# end-to-end; the post-stage's `--supplementary-coordinate-code` path is -# covered by the unit tests in RIFT/misc/coordinate_plugin.py. +# This is the configuration the IntrinsicPosterior-style coordinate +# argument mechanism (--parameter-implied / --parameter-nofit) enables +# inside the EOS posterior: implied = fit dim only, nofit = sample dim +# only, the coordinate plugin routes between them. # # Reading this file # ----------------- # Differences from hyperpipe_conf_tracer.yaml: -# * `marg-list[0].exe` -> example_gaussian_uvw.py (uses linear plugin) -# * `marg-list[0].args` ships the absolute paths to the plugin / ini so -# the marg job finds them regardless of the condor scratch cwd. -# * `general.rundir` renamed to keep the two tracer rundirs side-by-side. -# Everything else mirrors hyperpipe_conf_tracer.yaml exactly. +# * post.coords-fit is EMPTY -- no parameter is both fit and sampled. +# * post.coords-implied: "u v w" -- the FIT basis. The coordinate +# plugin produces these from the data file's (x, y, z) columns; +# the MC integrator never sees them. +# * post.coords-nofit: "x y z" -- the MC SAMPLING basis. These are the +# data file's columns; the integrator samples them, the puff lane +# puffs them, and the output posterior is in this basis. +# * post.coord-module: linear_coordinate_convert.py with the +# uvw_rotated chart, plus the ini that carries A and b. +# * marg-list[0].exe: example_gaussian_uvw.py, which uses the same +# plugin library to define its Gaussian likelihood in (u, v, w). # # Run with: # util_RIFT_hyperpipe.py --config ./hyperpipe_conf_linear_uvw.yaml # -# Then plot with: +# Then plot the (x, y, z) posterior from rundir_linear_uvw/: # (cd rundir_linear_uvw; plot_posterior_corner.py \ # --posterior-file posterior-4.dat \ # --composite-file all.marg_net --composite-file-has-labels \ # --parameter x --parameter y --parameter z \ # --lnL-cut 15 --use-all-composite-but-grayscale) +# +# The recovered ellipsoid in (x, y, z) should point along the rotated +# v-axis ((-1, +1, 0)/sqrt(2)) with the narrowest extent along the +# u-axis ((+1, +1, 0)/sqrt(2)). arch: method: default @@ -50,32 +52,39 @@ arch: start-iteration: 0 post: - # Iteration basis (== puff basis == data-file column basis). The post - # stage's GP/RF fit happens here too; the rotated (u, v, w) basis only - # appears inside the likelihood evaluator below. - coords-fit: "x y z" - # Integration ranges in (x, y, z). Sized to comfortably enclose the - # initial grid plus a margin for the Gaussian tails after rotation: - # the Gaussian peaks at (u, v, w) = (0, 4.95, 3.5), which back-projects - # to (x, y, z) ~= (-3.5, 3.5, 3.5), and the broadest sigma (sigma_w = - # 1.5 along z) needs ~+/-5 in z to contain ~3 sigma. - coords-sample: "x:[-7,1] y:[0,7] z:[-1,8]" + # FIT basis -- only seen by the GP/RF. Plugin produces these from the + # data-file columns. Nothing about implied coords is sampled by the MC, + # so they don't need integration ranges. + coords-implied: "u v w" + # MC SAMPLING basis -- same as the data file's columns. The integrator + # samples these and the output posterior is written in this basis. + coords-nofit: "x y z" + # Sampling-basis integration ranges sized to enclose ~3 sigma of the + # rotated Gaussian (centred at (u, v, w) = (0, 4.95, 3.5), which back- + # projects to (x, y, z) ~= (-3.5, 3.5, 3.5); broadest extent is along + # z with sigma_w = 1.5). + coords-sample: "x:[-7,1] y:[0,7] z:[-1,8]" + # The coordinate plugin sits in this directory. ${hydra:runtime.cwd} + # resolves to where you invoked util_RIFT_hyperpipe.py from -- run it + # from demo/hyperpipe/ for the path below to work. + coord-module: ${hydra:runtime.cwd}/linear_coordinate_convert.py + # --supplementary-coordinate-{ini,chart,function} land here. function + # defaults to convert_coordinates so we don't have to set it. + extra-args: "--supplementary-coordinate-ini ${hydra:runtime.cwd}/linear_coordinate_convert.ini --supplementary-coordinate-chart uvw_rotated" settings: fit-method: rf marg-list: - name: Gaussian_uvw exe: example_gaussian_uvw.py - # The marg job lives inside a condor scratch dir, so the plugin / - # ini paths need to be absolute. ${hydra:runtime.cwd} is Hydra's - # original-invocation cwd resolver -- this assumes you run - # util_RIFT_hyperpipe.py from inside demo/hyperpipe/. If you invoke - # from elsewhere, edit the next two args to point at wherever you - # put linear_coordinate_convert.{py,ini}. args: "--outdir Gaussian_uvw_example --conforming-output-name --coord-plugin ${hydra:runtime.cwd}/linear_coordinate_convert.py --coord-ini ${hydra:runtime.cwd}/linear_coordinate_convert.ini --mu-uvw 0.0,4.95,3.5 --sigma-uvw 0.5,1.0,1.5" n-chunk: 100 puff: + # Tracer drop-in. Reads all.marg_net (which has columns (x, y, z)) and + # produces grid-{k+1}.dat in (x, y, z). Sees --parameter x y z because + # to_puff_args emits the SAMPLING basis (coords-fit + coords-nofit), + # not the fit basis. exe: util_HyperparameterTracerUpdate.py input-source: marg_net puff-factor: 0.5 @@ -89,10 +98,6 @@ puff: rng-seed: null init: - # Reuse the same initial grid as the baseline tracer demo. Columns are - # (x, y, z); the Gaussian centre at (u, v, w) = (0, 4.95, 3.5) lies near - # the centroid of this grid after rotation, so the iteration has - # something to lock onto right from iteration zero. file: blind_gaussian_3d_xy_plus.dat general: From d1995aa30e6d5a5c3c78ea05a604efcc9948f020 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Tue, 2 Jun 2026 06:48:00 -0400 Subject: [PATCH 17/18] demo/hyperpipe: README tour of the four yaml configs --- .../Code/demo/hyperpipe/README.md | 187 +++++++++++++++++- 1 file changed, 180 insertions(+), 7 deletions(-) diff --git a/MonteCarloMarginalizeCode/Code/demo/hyperpipe/README.md b/MonteCarloMarginalizeCode/Code/demo/hyperpipe/README.md index 89e0c1bb..2002cae5 100644 --- a/MonteCarloMarginalizeCode/Code/demo/hyperpipe/README.md +++ b/MonteCarloMarginalizeCode/Code/demo/hyperpipe/README.md @@ -1,14 +1,187 @@ -README for HyperPipe code +# `demo/hyperpipe` — HyperPipe example pipelines -The Makefile contains commands for running this HyperPipe code for a sample Gaussian distribution, and instructions for how the makefile can be modified for a new user supplied code run. -For more details, see the auto-generated documentation [here](https://rift-documentation.readthedocs.io/en/latest/hyperpipe.html) for a description of hyperpipe. +This directory holds runnable demos for the RIFT HyperPipe pipeline. Each +configuration drives the iterative `MARG → CON → UNIFY → EOS_POST → PUFF → +TEST → next iteration` loop on the same 3-D Gaussian toy, but each one +exercises a different slice of the pipeline so you can pick the one closest +to what you are testing. +The auto-generated documentation at + +covers the pipeline design. `technical_doc.txt` in this directory has a +pedagogical writeup of the procedure and implementation. -Similar content is also available in the supplementary document 'technical_doc.txt' for a pedagogical overview of the proceedure and implementation of the pipeline. -# Hyperpipe hydra pipeline writer -New tool +## The driver + +All hydra-style demos in this directory are launched with the same tool: + +```bash +util_RIFT_hyperpipe.py --config ./ +``` + +`util_RIFT_hyperpipe.py` reads the yaml, emits `args_*.txt` files into the +rundir specified by `general.rundir`, calls `create_eos_posterior_pipeline` +to assemble the condor DAG (`marginalize_hyperparameters.dag`), and submits +it. Inspect any rundir's `.sub` files to see exactly what every executable +was invoked with. + +Configs in this directory: + +| Config | Rundir | What it exercises | +|-------------------------------------|-----------------------|-------------------------------------------------------------------------------------| +| `hyperpipe_conf.yaml` | `rundir/` | Baseline pipeline: posterior-only resampling (no puff lane). | +| `hyperpipe_conf_tracer.yaml` | `rundir_tracer/` | Parsimonious-placement (tracer) variant — MARG_PUFF lane suppressed. | +| `hyperpipe_conf_osg.yaml` | `rundir_osg/` | OSG / IGWN submit-host setup (containers, OSG attributes). | +| `hyperpipe_conf_linear_uvw.yaml` | `rundir_linear_uvw/` | Coordinate transformation: fit in `(u,v,w)` while sampling in `(x,y,z)`. | + + +## `hyperpipe_conf.yaml` — baseline + +**What it exercises.** The simplest end-to-end loop: MARG evaluates +likelihoods on the initial grid, the EOS posterior fits a GP/RF, the next +iteration draws candidate grid points by resampling the posterior. There +is no PUFF lane (no `puff.exe` set) and no tracer placement — the +posterior alone seeds the next iteration's grid. + +**Build.** + +```bash +util_RIFT_hyperpipe.py --config ./hyperpipe_conf.yaml +# equivalent to: make rundir +``` + +**Inspect.** After the DAG finishes: + +```bash +(cd rundir; plot_posterior_corner.py \ + --posterior-file posterior-2.dat \ + --composite-file all.marg_net --composite-file-has-labels \ + --parameter x --parameter y --parameter z \ + --lnL-cut 15 --use-all-composite-but-grayscale) +``` + +You should see the recovered isotropic Gaussian centred where +`example_gaussian.py` placed it (default `[-5, 0, 0]`). + + +## `hyperpipe_conf_tracer.yaml` — parsimonious placement + +**What it exercises.** The tracer pathway (`util_HyperparameterTracerUpdate.py` +as `puff.exe`). This consumes `all.marg_net` directly and writes +`grid-{k+1}.dat`, suppressing the MARG_PUFF lane entirely. Per the comment +in the config file, the saving is about 1.7–1.8× for `N=5–6` iterations on +this toy — MARG still runs every iteration, but no separate MARG_PUFF +lane is built. + +The yaml exposes the SMC-MALA / birth-death sampler hyperparameters in +`puff.settings` so you can twiddle them without escaping into +`extra-args`. + +**Build.** + +```bash +util_RIFT_hyperpipe.py --config ./hyperpipe_conf_tracer.yaml +# equivalent to: make rundir_tracer +``` + +**Inspect.** Same plot command as above with `--posterior-file +rundir_tracer/posterior-3.dat`, or use the shell helper: + +```bash +make rundir_tracer_plots ``` -util_RIFT_hyperpipe.py --config-name hyperpipe_conf_tracer.yaml + + +## `hyperpipe_conf_osg.yaml` — OSG / IGWN setup + +**What it exercises.** Submit-host configuration for a condor pool with +OSG attributes and an IGWN-prefixed condor-local nonworker setup. Same +3-D Gaussian toy as the baseline, but `general.use-osg`, +`general.condor-local-nonworker`, and +`general.condor-local-nonworker-igwn-prefix` are all enabled. Use this +as a starting point when adapting one of the other demos for OSG or LIGO +clusters. + +**Build.** + +```bash +util_RIFT_hyperpipe.py --config ./hyperpipe_conf_osg.yaml ``` +> **Note.** As of writing the OSG yaml has a couple of stale keys +> (`explode marg jobs`, `puff factor`) that hydra would reject because +> they use spaces where the schema expects hyphens. Fix those locally +> before running (`explode-marg-jobs`, `puff-factor`) or copy from the +> baseline yaml. + + +## `hyperpipe_conf_linear_uvw.yaml` — coordinate transformation + +**What it exercises.** The decoupled-bases path in +`util_ConstructEOSPosterior.py`: the iteration (puff, marg evaluator, +convergence test) operates in the data-file basis `(x, y, z)`, but the +EOS posterior step *fits its GP/RF in a transformed basis `(u, v, w)`* +routed through the `linear_coordinate_convert.py` plugin. The +likelihood evaluator (`example_gaussian_uvw.py`) uses the same plugin +library to define a Gaussian whose principal axes lie along `(u, v, w)` +with deliberately unequal sigmas — so the rotation is observable in the +recovered posterior. + +Two pieces of the new machinery are exercised together: + +1. **Coordinate plugin contract.** `post.coord-module: + linear_coordinate_convert.py` (plus `--supplementary-coordinate-ini` + and `--supplementary-coordinate-chart uvw_rotated` via + `post.extra-args`) walks the loader at + `RIFT.misc.coordinate_plugin.load_coordinate_converter`. + +2. **`--parameter-implied` / `--parameter-nofit` semantics in + EOSPosterior.** `post.coords-implied: u v w` is the fit basis, + `post.coords-nofit: x y z` is the MC sampling basis, and + `post.coords-fit` is empty. This is the IntrinsicPosterior-style + coord-arg mechanism, now wired through to EOSPosterior. + +**Build.** + +```bash +util_RIFT_hyperpipe.py --config ./hyperpipe_conf_linear_uvw.yaml +``` + +**Inspect.** The recovered `(x, y, z)` posterior should be an +elongated ellipsoid pointing along the rotated `v`-axis +(`(-1, +1, 0) / sqrt(2)`), with the narrowest extent along the +`u`-axis (`(+1, +1, 0) / sqrt(2)`): + +```bash +(cd rundir_linear_uvw; plot_posterior_corner.py \ + --posterior-file posterior-4.dat \ + --composite-file all.marg_net --composite-file-has-labels \ + --parameter x --parameter y --parameter z \ + --lnL-cut 15 --use-all-composite-but-grayscale) +``` + + +## Pre-hydra Makefile-style demos + +The Makefile also carries two targets that bypass `util_RIFT_hyperpipe.py` +and drive `create_eos_posterior_pipeline` directly: + +- `make Gaussian_adaptive_unimodal` — single-event adaptive run with the + baseline puffball. Useful when you want to inspect the + `args_*.txt` files by hand or when hydra is in the way. +- `make Gaussian_adaptive_bimodal` — two events (`example_gaussian.py` + + `example_gaussian2.py`) into the same posterior, illustrating + multi-event heterogeneous-driver mode. + +These predate the hydra wrapper; the hydra yamls above are the +recommended entry point for new work. + + +## Initial grid + +Every demo above seeds from `blind_gaussian_3d_xy_plus.dat` (1000 uniform +points covering a corner of the `[-7, 7]^3` cube). Regenerate it with +`make blind_gaussian_3d_xy_plus.dat` if you delete it. The Gaussian's +centre is configurable in `example_gaussian.py` — `make +change_center_location` rewrites it via `sed`. From 5643386932c83a8b7febe8aa8524fbfb18f419a8 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Tue, 2 Jun 2026 11:54:12 -0400 Subject: [PATCH 18/18] plot_posterior_corner: additive coordinate-plugin hook Adds --supplementary-coordinate-{code,function,ini,chart} plus the two input/output parameter list flags to plot_posterior_corner.py, mirroring the surface already in util_ConstructEOSPosterior.py. When the plugin flag is set, _materialize_plugin_columns runs once per loaded posterior and once per loaded composite file *after* the existing RIFT postprocessing -- it computes the plugin's output columns from existing record-array fields and splices them in via add_field. Critically, the hook is strictly ADDITIVE. Any output name already present in samples.dtype.names is skipped, so the hardcoded extract_combination_from_LI and the per-file postprocess loops (mc / eta / chi_eff / LambdaTilde / chi1_perp / ...) always win. Legacy invocations with no --supplementary-coordinate-code flag are byte- identical to the pre-plugin tool -- the helper returns input unchanged when the converter is None. CLI surface ----------- --supplementary-coordinate-code SPEC 'rift_default' | filesystem path to a .py | dotted module name. --supplementary-coordinate-function NAME Entry-point callable. Defaults to 'convert_coordinates'. --supplementary-coordinate-ini PATH Optional; parsed and handed to prepare(). --supplementary-coordinate-chart NAME Required only when the plugin defines multiple charts. --supplementary-coordinate-input-parameter NAME (action='append') Override the plugin-declared INPUT_PARAMETERS / chart's input_parameters list. --supplementary-coordinate-output-parameter NAME (action='append') Override the plugin-declared OUTPUT_PARAMETERS / chart's parameters list. When the input / output name lists aren't given on the CLI they're resolved from CHARTS[chart] (input_parameters, parameters) and then from the module-level INPUT_PARAMETERS / OUTPUT_PARAMETERS attributes. Verified -------- Five synthetic cases on an (m1, x, y, z) record array with the linear plugin requesting (u, v, w): * happy path -- (x, y, z) -> (u, v, w) values match u=(x+y)/sqrt(2), v=(y-x)/sqrt(2), w=z on every row; (m1, x, y, z) untouched. * output-reorder -- works regardless of the order u/v/w are listed. * name-collision -- samples pre-seeded with u=99; the plugin leaves u alone (RIFT path wins) and still adds v, w. * missing-input -- samples without an x column; helper logs a skip, returns input unchanged, no crash. * no-plugin -- helper is identity (out is samples). --- .../Code/bin/plot_posterior_corner.py | 165 +++++++++++++++++- 1 file changed, 164 insertions(+), 1 deletion(-) diff --git a/MonteCarloMarginalizeCode/Code/bin/plot_posterior_corner.py b/MonteCarloMarginalizeCode/Code/bin/plot_posterior_corner.py index 66af83e3..f6c2b466 100755 --- a/MonteCarloMarginalizeCode/Code/bin/plot_posterior_corner.py +++ b/MonteCarloMarginalizeCode/Code/bin/plot_posterior_corner.py @@ -254,6 +254,35 @@ def render_coordinates(coord_names,logparams=[]): parser.add_argument("--no-mod-psi",action="store_true",help="Default is to take psi mod pi. If present, does not do this") parser.add_argument("--downselect-parameter",action='append', help='Name of parameter to be used to eliminate grid points ') parser.add_argument("--downselect-parameter-range",action='append',type=str) +# ---- User-supplied coordinate-convert plugin (additive; the hardcoded RIFT +# conversion path is untouched). When --supplementary-coordinate-code is +# omitted these flags are a no-op and plotting behaves byte-identically +# to the legacy tool. When supplied, the plugin runs AFTER the existing +# per-file postprocessing -- it can only ADD columns to the record array, +# never override or shadow a name the RIFT path already produced. See +# RIFT.misc.coordinate_plugin for the plugin contract. +parser.add_argument("--supplementary-coordinate-code", default=None, type=str, + help="Coordinate conversion plugin spec. Accepts the literal 'rift_default', " + "a filesystem path to a .py file, or an importable dotted module name. " + "When set, the plugin is loaded and used to materialize additional " + "named columns on every loaded posterior / composite file -- columns " + "already produced by the RIFT hardcoded path (extract_combination_from_LI " + "etc.) are NOT overridden.") +parser.add_argument("--supplementary-coordinate-function", default=None, type=str, + help="Entry-point callable inside the plugin module. Defaults to 'convert_coordinates'.") +parser.add_argument("--supplementary-coordinate-ini", default=None, type=str, + help="Optional ini file parsed and handed to the plugin's prepare() hook.") +parser.add_argument("--supplementary-coordinate-chart", default=None, type=str, + help="Which chart (coordinate system) defined by the plugin to use. Required only " + "when the plugin defines multiple charts; otherwise auto-resolved.") +parser.add_argument("--supplementary-coordinate-input-parameter", action='append', default=None, + help="Existing posterior/composite column name to feed to the plugin as input. " + "Repeat for each input column. If omitted, the plugin's CHARTS[chart] " + "input_parameters / INPUT_PARAMETERS attribute is used.") +parser.add_argument("--supplementary-coordinate-output-parameter", action='append', default=None, + help="Name of a column the plugin should produce and add to the record arrays. " + "Repeat for each output column. If omitted, the plugin's CHARTS[chart] " + "parameters / OUTPUT_PARAMETERS attribute is used.") parser.add_argument("--verbose",action='store_true',help='print matplotlibrc data') opts= parser.parse_args() @@ -304,7 +333,131 @@ def render_coordinates(coord_names,logparams=[]): downselect_dict[dlist[indx]] = dlist_ranges[indx] if opts.downselect_parameter: print("Parameter downselect " , downselect_dict) - + + +# --------------------------------------------------------------------------- +# Optional coordinate-convert plugin. +# +# Loaded ONCE up front; the converter and the input/output name lists are +# captured as module-level closures so the per-file materialize step below +# is cheap. Skipped entirely when --supplementary-coordinate-code is unset, +# so the legacy invocation is byte-identical to the pre-plugin tool. +# +# Design note (don't break the hardcoded RIFT path): the plugin only ever +# ADDS columns to a record array. The `_materialize_plugin_columns` +# helper skips any output name that already exists in samples.dtype.names, +# so the converter cannot shadow a parameter produced by +# `extract_combination_from_LI` or by the file-load hot loop's +# `add_field`-based postprocessing. The standard `samples[param] ... else +# extract_combination_from_LI(...)` fallback at the plot-data extraction +# sites just transparently finds the new columns when the user asks for +# them. +# --------------------------------------------------------------------------- +_coord_plugin_converter = None +_coord_plugin_in_names = [] +_coord_plugin_out_names = [] +if opts.supplementary_coordinate_code: + from RIFT.misc.coordinate_plugin import load_coordinate_converter + + # Reuse the same loader util_ConstructEOSPosterior.py uses, so the + # plugin contract (CHARTS / prepare / register_priors / etc.) is + # exercised consistently across the codebase. prior_map/prior_range_map + # are unused here (plotting doesn't sample anything) -- pass None so the + # loader doesn't try to install priors we'd ignore. + _coord_plugin_converter, _coord_plugin_module = load_coordinate_converter( + spec=opts.supplementary_coordinate_code, + function_name=opts.supplementary_coordinate_function, + ini_path=opts.supplementary_coordinate_ini, + coord_names=opts.supplementary_coordinate_output_parameter, + low_level_coord_names=opts.supplementary_coordinate_input_parameter, + chart=opts.supplementary_coordinate_chart, + opts=opts, + prior_map=None, + prior_range_map=None, + ) + + # Resolve the input/output name lists. Priority: explicit CLI override + # > selected chart's declared names > module-level INPUT/OUTPUT_PARAMETERS. + _chart_spec = ( + getattr(_coord_plugin_module, "CHARTS", {}).get(opts.supplementary_coordinate_chart) + if opts.supplementary_coordinate_chart + else None + ) + if _chart_spec is None: + _charts = getattr(_coord_plugin_module, "CHARTS", None) or {} + if len(_charts) == 1: + _chart_spec = next(iter(_charts.values())) + + _coord_plugin_out_names = list( + opts.supplementary_coordinate_output_parameter + or (_chart_spec.get("parameters") if _chart_spec else None) + or getattr(_coord_plugin_module, "OUTPUT_PARAMETERS", []) + ) + _coord_plugin_in_names = list( + opts.supplementary_coordinate_input_parameter + or (_chart_spec.get("input_parameters") if _chart_spec else None) + or getattr(_coord_plugin_module, "INPUT_PARAMETERS", []) + ) + if not _coord_plugin_out_names or not _coord_plugin_in_names: + raise ValueError( + "Coordinate plugin loaded but the input/output column names " + "could not be determined from CHARTS / INPUT_PARAMETERS / " + "OUTPUT_PARAMETERS. Pass --supplementary-coordinate-input-parameter " + "and --supplementary-coordinate-output-parameter explicitly." + ) + print( + " plot_posterior_corner: coordinate plugin will materialize {} " + "from {} on every loaded posterior / composite file.".format( + _coord_plugin_out_names, _coord_plugin_in_names + ) + ) + + +def _materialize_plugin_columns(samples, source_label=""): + """Return ``samples`` with the plugin's output columns spliced in. + + Pure no-op when no plugin is loaded, when the plugin's required input + columns are missing from ``samples``, or when every output name is + already present (the legacy RIFT path wins by virtue of running first). + """ + if _coord_plugin_converter is None: + return samples + missing_in = [n for n in _coord_plugin_in_names if n not in samples.dtype.names] + if missing_in: + print( + " coordinate plugin: input column(s) {!r} missing from {!r}; " + "skipping conversion for this file.".format(missing_in, source_label or "samples") + ) + return samples + # Names already in samples are NOT overridden -- the RIFT hardcoded path + # (and any per-file postprocessing) wins. This is the only line that + # protects the legacy code path from a misconfigured plugin. + new_names = [n for n in _coord_plugin_out_names if n not in samples.dtype.names] + if not new_names: + return samples + x_in = np.column_stack( + [np.asarray(samples[n], dtype=float) for n in _coord_plugin_in_names] + ) + y = _coord_plugin_converter( + x_in, + coord_names=_coord_plugin_out_names, + low_level_coord_names=_coord_plugin_in_names, + ) + y = np.asarray(y, dtype=float) + if y.shape != (len(samples), len(_coord_plugin_out_names)): + raise ValueError( + "coordinate plugin returned shape {!r}, expected {!r}".format( + y.shape, (len(samples), len(_coord_plugin_out_names)) + ) + ) + out = samples + for name in new_names: + col_indx = _coord_plugin_out_names.index(name) + out = add_field(out, [(name, float)]) + out[name] = y[:, col_indx] + return out + + truth_P_list = None P_ref = None truth_dat = None @@ -481,6 +634,11 @@ def render_coordinates(coord_names,logparams=[]): samples = samples[indx_ok] + # Splice in any user-supplied coordinate-plugin output columns + # AFTER the legacy RIFT postprocessing has finished, never before. + # No-op when --supplementary-coordinate-code is unset. + samples = _materialize_plugin_columns(samples, source_label=fname) + # Save samples posterior_list.append(samples) @@ -658,6 +816,11 @@ def render_coordinates(coord_names,logparams=[]): samples = samples[indx_ok] + # Same plugin hook as the posterior loop above: only ADDS columns, + # never overrides what the composite-file postprocessing produced. + samples = _materialize_plugin_columns(samples, source_label=fname) + samples_orig = _materialize_plugin_columns(samples_orig, source_label=fname + " (pre-lnL-cut)") + composite_list.append(samples) composite_full_list.append(samples_orig)