CIP: give the master sample mask a name it keeps

.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

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)

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[<key>] = samples[<key>][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))