diff --git a/MonteCarloMarginalizeCode/Code/RIFT/calmarg/DESIGN_adaptive_driver.md b/MonteCarloMarginalizeCode/Code/RIFT/calmarg/DESIGN_adaptive_driver.md index 4922384c..d404ae28 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/calmarg/DESIGN_adaptive_driver.md +++ b/MonteCarloMarginalizeCode/Code/RIFT/calmarg/DESIGN_adaptive_driver.md @@ -2,6 +2,74 @@ 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: 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 - In-loop calmarg works and is validated (loop == fused == reference ~1e-14; CPU+GPU; @@ -135,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. 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/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/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/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/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/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 e34fd484..d74f1402 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: @@ -241,6 +242,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: 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") optp.add_option("--force-gpu-only", action="store_true", help="Hard fail if no GPU present (assessed by cupy not loading)") @@ -1678,7 +1681,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): @@ -1711,7 +1714,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 @@ -1754,7 +1757,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) @@ -2014,7 +2017,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 @@ -2086,6 +2089,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 @@ -2485,7 +2509,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/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) 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) 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/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/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/bin/util_RIFT_pseudo_pipe.py b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py index 5f8f1eba..c3ca93ce 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py @@ -165,6 +165,8 @@ 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("--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) @@ -822,6 +824,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 @@ -1040,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/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`. 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: diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile b/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile index c2faee73..888d9adc 100644 --- a/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile +++ b/MonteCarloMarginalizeCode/Code/demo/rift/calmarg/Makefile @@ -169,6 +169,272 @@ $(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 +# 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 +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 +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 + DAG_DMARG_DEP := distance_marg.npz +else + DAG_DMARG_FLAG := + DAG_DMARG_DEP := +endif + +ifeq ($(FUSED),1) + DAG_FUSED_FLAG := --calibration-fused-kernel +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 + DAG_PILOT_CEPP := +endif + +# 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 \ + --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 \ + --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) $(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 +# 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 +DAG_FILE := marginalize_intrinsic_parameters_BasicIterationWorkflow.dag + +.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 +# 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_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) \ + --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; 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)"; \ + [ "$(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) \ + --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 --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) + $(MAKE) dag-validate PILOT=$(PILOT) FUSED=$(FUSED) + +# 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) + @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: + 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 + +# --- 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 + +# --- 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 - rm -rf cal_env lowsnr_mdc + 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 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) -