AV: fix the high-SNR empty-live-volume crash and report the collapse; add the extrinsic zoom box

.travis/test-integrate.sh

@@ -1,5 +1,13 @@
 #! /bin/bash
 
+# Unit/regression tests for the sampler helpers themselves (fast, no data needed).
+# The extrinsic "zoom box" limits under the cosine samplers live here: they are pure
+# coordinate-transform + prior-mass identities, so they belong with the integrator gate
+# rather than with the end-to-end run tests.  The adaptive-volume empty-live-volume
+# regression is here for the same reason (it is a threshold identity, not a run test).
+python -m pytest -q MonteCarloMarginalizeCode/Code/test/test_limit_cosine_samplers.py
+python -m pytest -q MonteCarloMarginalizeCode/Code/test/test_av_empty_live_volume.py
+
 python MonteCarloMarginalizeCode/Code/test/test_mcsamplerEnsemble_extended.py --as-test --n-max 100000
 
 python MonteCarloMarginalizeCode/Code/test/test_mcsamplerEnsemble_extended.py --as-test --n-max 100000 --use-lnL

MonteCarloMarginalizeCode/Code/RIFT/integrators/mcsampler.py

@@ -876,6 +876,162 @@ def dec_samp_cdf_inv_vector(p):
     return numpy.arccos(2*p-1) - numpy.pi/2  # target from -pi/2 to pi/2
 
 
+###
+### Truncated isotropic angle samplers ("zoom box" support)
+###
+# RIFT can sample sky location / orientation either in the ANGLE itself (dec,
+# iota) or -- with --declination-cosine-sampler / --inclination-cosine-sampler --
+# in the cosine variable that makes the isotropic prior flat.  The two
+# conventions, read off the conversions applied in
+# bin/integrate_likelihood_extrinsic_batchmode, are
+#
+#   declination:  dec  = pi/2 - arccos(z)  <=>  z = sin(dec),   z in [-1,1]
+#                 sin() is INCREASING on [-pi/2, pi/2], so a box [lo,hi] maps to
+#                 [sin(lo), sin(hi)]: the order of the limits is PRESERVED.
+#   inclination:  iota = arccos(z)         <=>  z = cos(iota),  z in [-1,1]
+#                 cos() is DECREASING on [0, pi], so a box [lo,hi] maps to
+#                 [cos(hi), cos(lo)]: the order of the limits SWAPS.
+#
+# In both cases the physical isotropic prior is p(z) dz = dz/2, i.e. a CONSTANT
+# density 1/2 in the cosine coordinate.  That constant is deliberately NOT
+# renormalized over a restricted box, so that restricting the box reduces the
+# prior mass (and hence lnZ) by exactly the same factor as in the angle
+# coordinate, where the prior density is 0.5*cos(dec) resp. 0.5*sin(iota).
+
+_COSINE_SAMPLER_CONVENTIONS = {
+    # name: (angle_min, angle_max, angle->cosine-coordinate map, reverses_order)
+    'declination': (-numpy.pi/2, numpy.pi/2, numpy.sin, False),
+    'inclination': (0.0,         numpy.pi,   numpy.cos, True),
+}
+
+
+def infer_array_module(x, xpy=None):
+    """Return the array module (numpy, cupy, ...) that should be used to operate on `x`.
+
+    The truncated samplers below are handed to whichever backend the run selected:
+    mcsampler and mcsamplerAdaptiveVolume call them with host (numpy) arrays, while
+    mcsamplerGPU.draw_simplified() calls `self.pdf[param](samples)` / `self.cdf_inv[param](p)`
+    with a SINGLE positional argument holding a cupy array.  Defaulting to numpy would then
+    hit `numpy.asarray(cupy_array)`, which raises, so the backend is inferred from the
+    argument instead of assumed.  An explicit `xpy=` always wins.
+
+    The lookup is by the array type's top-level module, taken from `sys.modules` (a cupy
+    array cannot exist unless cupy is already imported), so this file keeps its numpy-only
+    import list and works for any duck-typed backend exposing the numpy API.
+    """
+    if xpy is not None:
+        return xpy
+    mod_name = type(x).__module__.split('.')[0]
+    if mod_name in ('numpy', 'builtins'):
+        return numpy
+    mod = sys.modules.get(mod_name)
+    if mod is not None and all(hasattr(mod, _attr) for _attr in ('asarray', 'where', 'clip')):
+        return mod
+    return numpy
+
+
+def clip_angle_limits(lo, hi, kind):
+    """Clip an angular range [lo,hi] (radians) to the physical domain of `kind`
+    ('declination' -> [-pi/2,pi/2], 'inclination' -> [0,pi]).
+
+    Raises ValueError if the requested range is empty/inverted, or if it does
+    not overlap the physical domain.  Returns (lo, hi) as floats with lo < hi.
+    """
+    if kind not in _COSINE_SAMPLER_CONVENTIONS:
+        raise ValueError("clip_angle_limits: unknown angle '{}' (expected one of {})".format(kind, sorted(_COSINE_SAMPLER_CONVENTIONS)))
+    angle_min, angle_max, _, _ = _COSINE_SAMPLER_CONVENTIONS[kind]
+    lo = float(lo)
+    hi = float(hi)
+    if not numpy.isfinite(lo) or not numpy.isfinite(hi):
+        raise ValueError("clip_angle_limits: non-finite {} range [{}, {}]".format(kind, lo, hi))
+    if not (hi > lo):
+        raise ValueError("clip_angle_limits: empty or inverted {} range [{}, {}] (need LO < HI, in radians)".format(kind, lo, hi))
+    lo_c = min(max(lo, angle_min), angle_max)
+    hi_c = min(max(hi, angle_min), angle_max)
+    if not (hi_c > lo_c):
+        raise ValueError("clip_angle_limits: {} range [{}, {}] does not overlap the physical domain [{}, {}]".format(kind, lo, hi, angle_min, angle_max))
+    return lo_c, hi_c
+
+
+def cosine_sampler_limits(lo, hi, kind):
+    """Map an angular range [lo,hi] (radians) into the coordinate actually sampled
+    by RIFT's 'cosine' sky/orientation samplers.
+
+    kind='declination': sampled variable is z = sin(dec); sin is increasing, so
+        the limit order is preserved:  [lo,hi] -> [sin(lo), sin(hi)].
+    kind='inclination': sampled variable is z = cos(iota); cos is DECREASING, so
+        the limit order SWAPS:         [lo,hi] -> [cos(hi), cos(lo)].
+
+    The range is clipped to the physical angular domain first, and the result is
+    clipped to the sampler domain [-1,1].  Raises ValueError on an empty or
+    inverted request.  Returns (z_lo, z_hi) with z_lo < z_hi.
+    """
+    lo_c, hi_c = clip_angle_limits(lo, hi, kind)
+    _, _, fn, reverses = _COSINE_SAMPLER_CONVENTIONS[kind]
+    z_a = float(fn(lo_c))
+    z_b = float(fn(hi_c))
+    z_lo, z_hi = (z_b, z_a) if reverses else (z_a, z_b)
+    z_lo = max(z_lo, -1.0)
+    z_hi = min(z_hi,  1.0)
+    if not (z_hi > z_lo):
+        raise ValueError("cosine_sampler_limits: {} range [{}, {}] maps to an empty sampling interval [{}, {}]".format(kind, lo, hi, z_lo, z_hi))
+    return z_lo, z_hi
+
+
+def ret_dec_samp_vector(dec_lo, dec_hi):
+    """Sampling pdf in DECLINATION for a uniform-in-sin(dec) draw truncated to
+    [dec_lo, dec_hi].  Normalized to unity over that box (the samplers that use
+    an explicit cdf_inv do not renormalize the pdf themselves).  Reduces to
+    dec_samp_vector for the full range."""
+    z_lo, z_hi = cosine_sampler_limits(dec_lo, dec_hi, 'declination')
+    lo_c, hi_c = clip_angle_limits(dec_lo, dec_hi, 'declination')
+    norm = z_hi - z_lo
+    def _pdf(x, xpy=None):
+        xpy = infer_array_module(x, xpy)
+        x = xpy.asarray(x, dtype=numpy.float64)
+        vals = xpy.cos(x)/norm
+        return xpy.where((x >= lo_c) & (x <= hi_c), vals, xpy.zeros_like(vals))
+    return _pdf
+
+
+def ret_dec_samp_cdf_inv_vector(dec_lo, dec_hi):
+    """Inverse CDF (p in [0,1] -> declination) for uniform-in-sin(dec) truncated
+    to [dec_lo, dec_hi].  Monotonically increasing in p."""
+    z_lo, z_hi = cosine_sampler_limits(dec_lo, dec_hi, 'declination')
+    def _cdf_inv(p, xpy=None):
+        xpy = infer_array_module(p, xpy)
+        p = xpy.asarray(p, dtype=numpy.float64)
+        return xpy.arcsin(xpy.clip(z_lo + p*(z_hi - z_lo), -1.0, 1.0))
+    return _cdf_inv
+
+
+def ret_cos_samp_vector(incl_lo, incl_hi):
+    """Sampling pdf in INCLINATION for a uniform-in-cos(iota) draw truncated to
+    [incl_lo, incl_hi].  Normalized to unity over that box.  Reduces to
+    cos_samp_vector for the full range."""
+    z_lo, z_hi = cosine_sampler_limits(incl_lo, incl_hi, 'inclination')
+    lo_c, hi_c = clip_angle_limits(incl_lo, incl_hi, 'inclination')
+    norm = z_hi - z_lo
+    def _pdf(x, xpy=None):
+        xpy = infer_array_module(x, xpy)
+        x = xpy.asarray(x, dtype=numpy.float64)
+        vals = xpy.sin(x)/norm
+        return xpy.where((x >= lo_c) & (x <= hi_c), vals, xpy.zeros_like(vals))
+    return _pdf
+
+
+def ret_cos_samp_cdf_inv_vector(incl_lo, incl_hi):
+    """Inverse CDF (p in [0,1] -> inclination) for uniform-in-cos(iota)
+    truncated to [incl_lo, incl_hi].  Monotonically increasing in p: p=0 gives
+    incl_lo (which is arccos of the UPPER cosine limit -- note the swap)."""
+    z_lo, z_hi = cosine_sampler_limits(incl_lo, incl_hi, 'inclination')
+    def _cdf_inv(p, xpy=None):
+        xpy = infer_array_module(p, xpy)
+        p = xpy.asarray(p, dtype=numpy.float64)
+        return xpy.arccos(xpy.clip(z_hi - p*(z_hi - z_lo), -1.0, 1.0))
+    return _cdf_inv
+
+
 def pseudo_dist_samp(r0,r):
         return r*r*numpy.exp( - (r0/r)*(r0/r)/2. + r0/r)+0.01  # put a floor on probability, so we converge. Note this floor only cuts out NEARBY distances
 

MonteCarloMarginalizeCode/Code/RIFT/integrators/mcsamplerAdaptiveVolume.py

@@ -94,6 +94,57 @@ def __init__(self, value):
     def __str__(self):
         return repr(self.value)
 
+class LiveVolumeCollapse(Exception):
+    """The adaptive-volume live set is empty (or carries no usable information).
+
+    Raised INSTEAD of the bare numpy/cupy "zero-size array to reduction operation
+    ... which has no identity" that the empty-live-volume path used to produce, so
+    callers and logs can tell a degenerate contraction apart from a waveform
+    generation failure.  See the collapse discussion in get_likelihood_threshold.
+    """
+    pass
+
+
+def ess_from_log_weights(log_wt):
+    """Kish effective sample size (sum w)^2 / sum w^2, from LOG weights."""
+    lw = np.asarray(log_wt, dtype=float)
+    lw = lw[np.isfinite(lw)]
+    if len(lw) == 0:
+        return 0.0
+    lse = special.logsumexp
+    return float(np.exp(2 * lse(lw) - lse(2 * lw)))
+
+
+def live_volume_collapse_verdict(n_live, ndim, ess=None):
+    """Has the adaptive-volume live set degenerated?  -> (collapsed, [reasons])
+
+    A degenerate contraction must be REPORTED rather than silently exported.  With the
+    threshold clamp in place the run no longer crashes, so it now returns a lnZ and a
+    sample cloud -- but at high SNR both can describe a SINGLE mode, and the cloud is
+    then not a fair posterior draw.  Unreported, that turns a crash (which at least
+    kept the point out of the posterior) into a silent contamination, which is worse.
+
+    Thresholds, against the separation measured on zero-noise injections at a fixed
+    intrinsic point (rho_net 51 -> 147):
+        healthy    ESS 16.3-34.3   (rho 51.4, converges cold)
+        collapsed  ESS 1.0-2.0     (rho 103-147)
+      * n_live <= ndim -- a live set no larger than the dimension cannot span the
+        space, let alone describe a posterior in it.  Geometric, not tuned.
+      * ESS < 2 -- fewer than two effective samples IS one sample.
+    The gap between the regimes is an order of magnitude wide, so these sit far from
+    both sides of it.  Pareto k-hat is deliberately NOT used: it exceeds its nominal
+    0.70 "unresolved tail" threshold even in the healthy runs on this problem (the 12
+    converged rho=51.4 replicates measure 0.819-1.605), so gating on it would have
+    false-flagged 12 of 12 good exports.
+    """
+    reasons = []
+    if n_live <= ndim:
+        reasons.append("final live volume holds {} sample(s) in {} dimensions".format(n_live, ndim))
+    if ess is not None and ess < 2.0:
+        reasons.append("ESS={:.2f}".format(ess))
+    return bool(reasons), reasons
+
+
 ### V. Tiwari routines
 
 def get_likelihood_threshold(lkl, lkl_thr, nsel, discard_prob,xpy_here=xpy_default):
@@ -104,7 +155,18 @@ def get_likelihood_threshold(lkl, lkl_thr, nsel, discard_prob,xpy_here=xpy_defau
     nsel : integer, has to do with size of array of likelihoods used to evaluate for next array.
     discard_prob: threshold on CDF to throw away an entire bin.  Should be very small
     """
-
+    if len(lkl) == 0:
+        # Caller must not ask for a threshold on an empty live volume: every reduction
+        # below (max, argsort, [0]) is undefined, and the bare backend error for that
+        # ("zero-size array to reduction operation CUPY_CUB_MAX which has no identity")
+        # is what the ILE used to mis-report as a waveform problem.  Name the cause.
+        raise LiveVolumeCollapse(
+            "adaptive-volume live set is empty: the likelihood returned no finite value "
+            "inside the sampled volume.  At high network SNR this is likelihood UNDERFLOW "
+            "(exp() of a lnL more than ~745 nats below the peak returns 0), so a cold "
+            "extrinsic prior yields almost no usable draw.  This is NOT a waveform Nyquist/"
+            "start-frequency/duration problem.  Narrow the extrinsic prior or seed the sampler.")
+
     w = xpy_here.exp(lkl - np.max(lkl))
     npoints = len(w)
     sumw = xpy_here.sum(w)
@@ -121,6 +183,37 @@ def get_likelihood_threshold(lkl, lkl_thr, nsel, discard_prob,xpy_here=xpy_defau
         lkl_stop_thr = lkl_stop_thr[-1]
     lkl_thr = min(lkl_stop_thr, prob_stop_thr)
 
+    # CLAMP (backport of rift_O4d PR #63).  The threshold is applied downstream as a
+    # STRICT `lkl > thr`, so a threshold at or above max(lkl) discards the ENTIRE live
+    # volume; the next cycle then reduces over an empty array and raises
+    #   "zero-size array to reduction operation CUPY_CUB_MAX which has no identity",
+    # which the ILE driver misreports as a waveform problem.
+    #
+    # It happens whenever the live set is small AND one weight dominates -- i.e. at high
+    # network SNR, where a double-precision exp() of a lnL more than ~745 nats below the
+    # peak underflows to zero, so almost every cold draw is -inf and the live set starts
+    # with a handful of members.  Then prob_stop_thr saturates at max(lkl) (every other
+    # weight is 0, so the discard_prob quantile IS the top sample) while lkl_stop_thr
+    # falls back to the smallest.  Each cycle discards at least one point regardless of
+    # merit and the live set ratchets to zero.
+    #
+    # Back the threshold off to the largest value strictly below the maximum, so the peak
+    # always survives.  In a healthy run len(lkl) >> nsel and lkl_stop_thr is the nsel-th
+    # largest, far below the max, so this clamp never engages and results are unchanged.
+    #
+    # Reduce on the ACTIVE backend and move only the scalar: identity_convert(lkl) here
+    # would copy the whole live set device->host every cycle, on the healthy path too.
+    lkl_max = float(identity_convert(xpy_here.max(lkl)))
+    if not (float(identity_convert(lkl_thr)) < lkl_max):
+        lkl_host = identity_convert(lkl)   # rare branch: the live set is tiny by construction
+        below = lkl_host[lkl_host < lkl_max]
+        if len(below):
+            lkl_thr = np.max(below)   # keep only the maximum: maximal, but safe, contraction
+        else:
+            # every live point is at the maximum: no threshold can separate them, so
+            # take one below all of them and leave the live volume intact.
+            lkl_thr = np.nextafter(lkl_max, -np.inf)
+
     truncp = xpy_here.sum(w[lkl < lkl_thr]) / sumw
 
     return identity_convert(lkl_thr), identity_convert(truncp)  # send both to CPU as needed
@@ -686,6 +779,29 @@ def integrate_log(self, lnF, *args, xpy=xpy_default,**kwargs):
               self._rvs[name] = identity_convert(self._rvs[name])   # this is trivial if xpy_default is numpy, and a conversion otherwise
 
         dict_return = {}
+
+        # ------------------------------------------------------------------
+        # LIVE-VOLUME COLLAPSE VERDICT.  The threshold clamp in get_likelihood_threshold
+        # stops the empty-reduction CRASH, but the run that used to crash now completes and
+        # exports -- sometimes from a single surviving sample.  Silently exporting that is
+        # worse than the crash was: a crashed export writes nothing and is visibly missing,
+        # whereas a degenerate one enters downstream posterior assembly looking ordinary.
+        # So say so.  Thresholds and the measured regimes are documented on
+        # live_volume_collapse_verdict.  Purely additive: no sampling decision reads this.
+        n_live_final = int(len(log_wt))
+        _ess = ess_from_log_weights(log_wt)
+        collapsed, _reasons = live_volume_collapse_verdict(n_live_final, ndim, ess=_ess)
+        dict_return['n_ESS'] = _ess
+        dict_return['n_live_final'] = n_live_final
+        dict_return['live_volume_collapsed'] = collapsed
+        if collapsed:
+            dict_return['collapse_reason'] = "; ".join(_reasons)
+            print(" [AV COLLAPSE] the live volume degenerated: " + dict_return['collapse_reason'] + ".")
+            print(" [AV COLLAPSE] lnZ and the exported samples describe a SINGLE mode of the integrand and are")
+            print(" [AV COLLAPSE] NOT a fair draw from the posterior.  Do not use this export unweighted.")
+            print(" [AV COLLAPSE] At high network SNR this is likelihood underflow over a cold extrinsic prior;")
+            print(" [AV COLLAPSE] narrow the extrinsic prior (--limit-* zoom box) or seed the sampler.")
+
         return log_int, np.log(rel_var)  +2*log_int, eff_samp, dict_return
 
         # if outvals: