diff --git a/.travis/test-posterior.sh b/.travis/test-posterior.sh index a15fd4e29..b1e384dcd 100755 --- a/.travis/test-posterior.sh +++ b/.travis/test-posterior.sh @@ -1,5 +1,11 @@ #!/bin/bash +# Static guard on the posterior drivers, ahead of the end-to-end arms below: `indx_ok` keeps +# one meaning after the master post-integration mask is applied. Pure ast, no RIFT import, +# ~2 s. The defect is invisible to the runs below, because the rebindings it forbids were +# never read as the master mask -- it is a trap for the next edit, not a live bug. +python -m pytest -q MonteCarloMarginalizeCode/Code/test/test_cip_indx_ok_scope.py + python .travis/make_fake_composite.py # Test default sampler (constant fit) util_ConstructIntrinsicPosterior_GenericCoordinates.py --fname fake.composite --parameter mtot --parameter q --parameter s1z --parameter s2z --use-precessing --no-plots diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py index 4880eea85..f13d89fb8 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py @@ -3726,20 +3726,26 @@ def parse_corr_params(my_str): s2z = samples['chiz_plus'] - samples['chiz_minus'] val1 = np.array(s1z**2+samples["s1y"]**2 + samples["s1x"]**2,dtype=internal_dtype); chi1 = np.sqrt(val1) val2 = np.array(s2z**2+samples["s2y"]**2 + samples["s2x"]**2,dtype=internal_dtype); chi2= np.sqrt(val2) - indx_ok = np.logical_and(chi1<=chi_max , chi2<=chi_small_max) - weights[ np.logical_not(indx_ok)] = 0 # Zero out failing samples. Has effect of fixing prior range! - weights[indx_ok] *= 9.*(chi_max**2 * chi_small_max**2)/(chi1*chi1*chi2*chi2)[indx_ok] + indx_in_range = np.logical_and(chi1<=chi_max , chi2<=chi_small_max) + weights[ np.logical_not(indx_in_range)] = 0 # Zero out failing samples. Has effect of fixing prior range! + weights[indx_in_range] *= 9.*(chi_max**2 * chi_small_max**2)/(chi1*chi1*chi2*chi2)[indx_in_range] +# DEAD BRANCH, KEPT DEAD ON PURPOSE. This guard repeats the one above, so this body never +# runs. It is the alternate-sampling twin -- it divides prior_weight out, the branch above +# does not -- so the `not` looks like a copy-paste slip. Do NOT drop it: the body reads +# samples["s1x"]/["s1y"], which an aligned-spin chiz_plus coordinate set does not have, so +# dropping it crashes a CLI combination that completes today. Reviving this needs the body +# fixed and a physics decision, not a guard edit. Measured in junior PR #357. elif opts.pseudo_uniform_magnitude_prior and 'chiz_plus' in samples.keys() and not opts.pseudo_uniform_magnitude_prior_alternate_sampling: s1z = samples['chiz_plus'] + samples['chiz_minus'] s2z = samples['chiz_plus'] - samples['chiz_minus'] val1 = np.array(s1z**2+samples["s1y"]**2 + samples["s1x"]**2,dtype=internal_dtype); chi1 = np.sqrt(val1) val2 = np.array(s2z**2+samples["s2y"]**2 + samples["s2x"]**2,dtype=internal_dtype); chi2= np.sqrt(val2) - indx_ok = np.logical_and(chi1<=chi_max , chi2<=chi_small_max) - weights[ np.logical_not(indx_ok)] = 0 # Zero out failing samples. Has effect of fixing prior range! + indx_in_range = np.logical_and(chi1<=chi_max , chi2<=chi_small_max) + weights[ np.logical_not(indx_in_range)] = 0 # Zero out failing samples. Has effect of fixing prior range! prior_weight = np.prod([prior_map[x](samples[x]) for x in ['s1x','s1y', 's2x', 's2y','chiz_plus','chiz_minus'] ],axis=0) - indx_ok = np.logical_and(indx_ok, divisible_sampling_density(prior_weight, len(weights))) - weights[ np.logical_not(indx_ok)] = 0 - weights[indx_ok] *= 9.*(chi_max**2 * chi_small_max**2)/(chi1*chi1*chi2*chi2)[indx_ok]/prior_weight[indx_ok] # undo chizplus, chizminus prior + indx_in_range = np.logical_and(indx_in_range, divisible_sampling_density(prior_weight, len(weights))) + weights[ np.logical_not(indx_in_range)] = 0 + weights[indx_in_range] *= 9.*(chi_max**2 * chi_small_max**2)/(chi1*chi1*chi2*chi2)[indx_in_range]/prior_weight[indx_in_range] # undo chizplus, chizminus prior # If we are using alignedspin-zprior AND chiz+, chiz-, then we need to reweight .. that prior cannot be evaluated internally @@ -3749,10 +3755,10 @@ def parse_corr_params(my_str): prior_weight = np.prod([prior_map[x](samples[x]) for x in ['chiz_plus','chiz_minus'] ],axis=0) s1z = samples['chiz_plus'] + samples['chiz_minus'] s2z =samples['chiz_plus'] - samples['chiz_minus'] - indx_ok = np.logical_and(np.abs(s1z)<=chi_max , np.abs(s2z)<=chi_max) - indx_ok = np.logical_and(indx_ok, divisible_sampling_density(prior_weight, len(weights))) - weights[ np.logical_not(indx_ok)] = 0 # Zero out failing samples. Has effect of fixing prior range! - weights[indx_ok] *= s_component_zprior( s1z[indx_ok])*s_component_zprior(s2z[indx_ok])/(prior_weight[indx_ok]) # correct for uniform + indx_in_range = np.logical_and(np.abs(s1z)<=chi_max , np.abs(s2z)<=chi_max) + indx_in_range = np.logical_and(indx_in_range, divisible_sampling_density(prior_weight, len(weights))) + weights[ np.logical_not(indx_in_range)] = 0 # Zero out failing samples. Has effect of fixing prior range! + weights[indx_in_range] *= s_component_zprior( s1z[indx_in_range])*s_component_zprior(s2z[indx_in_range])/(prior_weight[indx_in_range]) # correct for uniform if opts.pseudo_gaussian_mass_prior: # mass normalization (assuming mc, eta limits are bounds - as is invariably the case) @@ -3796,7 +3802,7 @@ def parse_corr_params(my_str): # Integral result v2: using modified prior. # 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 +# - use for Bayes factors with GREAT CARE for this reason; should correct for with the indx_in_range cuts above # Same absolute-scale restoration as for the integral above: lnLmax here is a maximum of the # CENTRED integrand, and this file is documented to agree with integral_result.dat -- so leaving the # plugin's constant out of one and not the other turns a check into a spurious disagreement. @@ -3989,9 +3995,9 @@ def parse_corr_params(my_str): fig_base = overlay_or_warn(dat_out_low_level_coord_names, range_here, low_level_coord_names, "input grid", weights=np.ones(len(X))/len(X), plot_datapoints=True,plot_density=False,plot_contours=False,quantiles=None,fig=fig_base, data_kwargs={'color':my_cmap_values},hist_kwargs={'color':'g', 'linestyle':'dashed'}) # TRUNCATED data set used here - indx_ok = Y > Y.max() - scipy.stats.chi2.isf(0.1,len(low_level_coord_names))/2 # approximate threshold for significant points,from inverse cdf 90% - n_ok = np.sum(indx_ok) - fig_base = overlay_or_warn(dat_out_low_level_coord_names[indx_ok], range_here, low_level_coord_names, "significant grid points", weights=np.ones(n_ok)*1.0/n_ok, plot_datapoints=True,plot_density=False,plot_contours=False,quantiles=None,fig=fig_base, data_kwargs={'color':'b'},hist_kwargs={'color':'b', 'linestyle':'dashed'}) + indx_significant = Y > Y.max() - scipy.stats.chi2.isf(0.1,len(low_level_coord_names))/2 # approximate threshold for significant points,from inverse cdf 90% + n_ok = np.sum(indx_significant) + fig_base = overlay_or_warn(dat_out_low_level_coord_names[indx_significant], range_here, low_level_coord_names, "significant grid points", weights=np.ones(n_ok)*1.0/n_ok, plot_datapoints=True,plot_density=False,plot_contours=False,quantiles=None,fig=fig_base, data_kwargs={'color':'b'},hist_kwargs={'color':'b', 'linestyle':'dashed'}) #except: else: @@ -4374,9 +4380,9 @@ def parse_corr_params(my_str): # BEFORE truncation, note, to highlight region explored. ONLY for this plot fig_base = corner.corner(X_orig, weights=np.ones(len(X_orig))/len(X_orig),plot_datapoints=True,plot_density=False,plot_contours=False,quantiles=None,fig=fig_base, data_kwargs={'color':'g'},hist_kwargs={'color':'g', 'linestyle':'dashed'},range=range_here) # A subset of the truncated data set - indx_ok = Y > Y.max() - scipy.stats.chi2.isf(0.1,len(low_level_coord_names))/2 # approximate threshold for significant points,from inverse cdf 90% - n_ok = np.sum(indx_ok) - fig_base = corner.corner(X[indx_ok],weights=np.ones(n_ok)*1.0/n_ok, plot_datapoints=True,plot_density=False,plot_contours=False,quantiles=None,fig=fig_base, data_kwargs={'color':'r'},hist_kwargs={'color':'b', 'linestyle':'dashed'},range=range_here) + indx_significant = Y > Y.max() - scipy.stats.chi2.isf(0.1,len(low_level_coord_names))/2 # approximate threshold for significant points,from inverse cdf 90% + n_ok = np.sum(indx_significant) + fig_base = corner.corner(X[indx_significant],weights=np.ones(n_ok)*1.0/n_ok, plot_datapoints=True,plot_density=False,plot_contours=False,quantiles=None,fig=fig_base, data_kwargs={'color':'r'},hist_kwargs={'color':'b', 'linestyle':'dashed'},range=range_here) plt.legend(handles=line_handles, bbox_to_anchor=corner_legend_location, prop=corner_legend_prop,loc=4) @@ -4475,10 +4481,10 @@ def parse_corr_params(my_str): print(" Rendering past samples for ", extra_plot_coord_names[indx], " based on ", len(dat_points_here)) fig_base = overlay_or_warn(dat_points_here, range_here, coord_names_here, "input grid", weights=np.ones(len(dat_points_here))*1.0/len(dat_points_here), plot_datapoints=True,plot_density=False,plot_contours=False,quantiles=None,fig=fig_base, data_kwargs={'color':'g'},hist_kwargs={'color':'g', 'linestyle':'dashed'}) # Render points available. Note we use the ORIGINAL data set, and truncate it - indx_ok = Y_orig > Y_orig.max() - scipy.stats.chi2.isf(0.1,len(low_level_coord_names))/2 # approximate threshold for significant points,from inverse cdf 90% - n_ok = np.sum(indx_ok) + indx_significant = Y_orig > Y_orig.max() - scipy.stats.chi2.isf(0.1,len(low_level_coord_names))/2 # approximate threshold for significant points,from inverse cdf 90% + n_ok = np.sum(indx_significant) print(" Adding points for figure ", n_ok, extra_plot_coord_names[indx], " drawn from original ") - fig_base = overlay_or_warn(dat_points_here[indx_ok], range_here, coord_names_here, "significant grid points", weights=np.ones(n_ok)*1.0/n_ok, plot_datapoints=True,plot_density=False,plot_contours=False,quantiles=None,fig=fig_base, data_kwargs={'color':'b'},hist_kwargs={'color':'b', 'linestyle':'dashed'}) + fig_base = overlay_or_warn(dat_points_here[indx_significant], range_here, coord_names_here, "significant grid points", weights=np.ones(n_ok)*1.0/n_ok, plot_datapoints=True,plot_density=False,plot_contours=False,quantiles=None,fig=fig_base, data_kwargs={'color':'b'},hist_kwargs={'color':'b', 'linestyle':'dashed'}) plt.legend(handles=line_handles, bbox_to_anchor=corner_legend_location, prop=corner_legend_prop,loc=4) diff --git a/MonteCarloMarginalizeCode/Code/test/test_cip_indx_ok_scope.py b/MonteCarloMarginalizeCode/Code/test/test_cip_indx_ok_scope.py new file mode 100644 index 000000000..c813b44ed --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_cip_indx_ok_scope.py @@ -0,0 +1,158 @@ +"""`indx_ok` must keep one meaning in the posterior drivers. + +Both drivers build one master post-integration mask and apply it to the whole sample set: + + indx_ok = np.logical_and(dat_logL > lnLmax - opts.lnL_offset, samples["joint_s_prior"] > 0) + ... + dat_logL = dat_logL[indx_ok] + +Downstream blocks used to rebind that same name for unrelated throwaway cuts -- spin-prior +range cuts, and the "significant points" subset in the corner plots. Nothing read the master +mask after its four applications, so this was never a live bug; it is a trap. A later block +that indexes with `indx_ok` when the nearest rebinding is a few lines further up gets the +master mask instead, with the right length and the wrong meaning, and no error. + +These tests are STATIC (ast-based). The drivers are top-level scripts, not importable +modules, so the block cannot be exercised without a full inference job. +""" +import ast +import os + +import pytest + +HERE = os.path.dirname(os.path.abspath(__file__)) +BIN = os.path.abspath(os.path.join(HERE, "..", "bin")) +DRIVERS = [ + "util_ConstructIntrinsicPosterior_GenericCoordinates.py", + "util_ConstructEOSPosterior.py", +] +NAME = "indx_ok" + +SCOPES = (ast.FunctionDef, ast.AsyncFunctionDef, ast.Lambda, ast.ClassDef) + + +def _tree(fname): + path = os.path.join(BIN, fname) + # NOT pytest.skip: a missing driver disarms every assertion below, and a skip exits 0. + assert os.path.exists(path), ( + "%s: driver not found at %s. This guard cannot run, which is a failure, not a pass " + "-- if the driver moved, update BIN/DRIVERS." % (fname, path)) + with open(path) as f: + return ast.parse(f.read(), filename=path) + + +def _refs(node, _top=True): + """Module-level Name nodes for NAME under `node`. + + Nested scopes have their own `indx_ok` locals (the fit wrappers near the top of the CIP + driver), which are fine and are skipped -- INCLUDING when the statement handed in is + itself a def/class, which an earlier version of this descended into. A nested scope that + declares `global indx_ok` is not local, so those are reported (see _global_rebinds). + """ + if _top and isinstance(node, SCOPES): + return + for ch in ast.iter_child_nodes(node): + if isinstance(ch, SCOPES): + continue + if isinstance(ch, ast.Name) and ch.id == NAME: + yield ch + for sub in _refs(ch, _top=False): + yield sub + + +def _global_rebinds(tree): + """Nested scopes that declare `global indx_ok` AND assign it -- a module-level rebind.""" + out = [] + for node in ast.walk(tree): + if not isinstance(node, SCOPES): + continue + declares = any(isinstance(n, ast.Global) and NAME in n.names for n in ast.walk(node)) + if not declares: + continue + for n in ast.walk(node): + if isinstance(n, ast.Name) and n.id == NAME and isinstance(n.ctx, ast.Store): + out.append(n.lineno) + return sorted(set(out)) + + +def _anchor(tree): + """Index of the module statement `dat_logL = dat_logL[indx_ok]`, the mask's defining use.""" + for i, stmt in enumerate(tree.body): + if not isinstance(stmt, ast.Assign) or len(stmt.targets) != 1: + continue + tgt, val = stmt.targets[0], stmt.value + if (isinstance(tgt, ast.Name) and tgt.id == "dat_logL" + and isinstance(val, ast.Subscript) + and isinstance(val.value, ast.Name) and val.value.id == "dat_logL" + and any(_refs(val))): + return i + return None + + +def _is_sample_remask(stmt): + """`samples[] = samples[][indx_ok]` -- the only shape allowed to follow the anchor. + + The object being re-masked must be `samples`. Accepting any `Subscript = Subscript[Name]` + let `prior_range_map['mc'] = weights[indx_ok]` through, which is exactly the consumer this + file exists to catch. + """ + if not (isinstance(stmt, ast.Assign) and len(stmt.targets) == 1): + return False + tgt, val = stmt.targets[0], stmt.value + if not (isinstance(tgt, ast.Subscript) + and isinstance(tgt.value, ast.Name) and tgt.value.id == "samples"): + return False + return (isinstance(val, ast.Subscript) + and isinstance(val.value, ast.Subscript) + and isinstance(val.value.value, ast.Name) and val.value.value.id == "samples" + and isinstance(val.slice, ast.Name) and val.slice.id == NAME + and len(list(_refs(stmt))) == 1) + + +@pytest.mark.parametrize("fname", DRIVERS) +def test_master_mask_is_present(fname): + """The anchor must exist -- otherwise every other assertion here passes vacuously.""" + tree = _tree(fname) + assert _anchor(tree) is not None, ( + "%s: no module-level `dat_logL = dat_logL[%s]`; the master mask was renamed or " + "removed, and this guard no longer guards anything" % (fname, NAME)) + + +@pytest.mark.parametrize("fname", DRIVERS) +def test_master_mask_is_never_rebound(fname): + """No module-level rebinding of `indx_ok` after the master mask has been applied.""" + tree = _tree(fname) + i = _anchor(tree) + stores = [(n.lineno, ast.unparse(stmt).splitlines()[0][:90]) + for stmt in tree.body[i + 1:] + for n in _refs(stmt) if isinstance(n.ctx, ast.Store)] + assert not stores, ( + "%s: `%s` is rebound after the master mask is applied, at %s. Give the local cut its " + "own name (e.g. indx_in_range, indx_significant)." % (fname, NAME, stores)) + + +@pytest.mark.parametrize("fname", DRIVERS) +def test_only_remasking_reads_the_master_mask(fname): + """After the anchor, `indx_ok` may only re-mask another `samples[...]` entry.""" + tree = _tree(fname) + i = _anchor(tree) + bad = [(next(_refs(stmt)).lineno, ast.unparse(stmt).splitlines()[0][:90]) + for stmt in tree.body[i + 1:] + if any(_refs(stmt)) and not _is_sample_remask(stmt)] + assert not bad, ( + "%s: `%s` is read after the master mask is applied, at %s. A later block that indexes " + "with this name gets the master mask, not a local cut." % (fname, NAME, bad)) + + +@pytest.mark.parametrize("fname", DRIVERS) +def test_no_nested_scope_rebinds_the_master_mask(fname): + """A helper may keep a local `indx_ok`; it may not `global indx_ok` and assign it. + + The three tests above read module-level statements only, so a `global` rebind inside a + function is invisible to them while still changing the mask everything downstream uses. + """ + tree = _tree(fname) + lines = _global_rebinds(tree) + assert not lines, ( + "%s: a nested scope declares `global %s` and assigns it, at line(s) %s. That rebinds " + "the master mask from inside a helper." % (fname, NAME, lines))