distance slices: cupy-safe export (GPU runs were silently skipping every binary)

On a GPU run (--gpu/--force-xpy) the sampler stores CUPY arrays in _rvs and the
cached likelihood returns cupy. The .dslice export did plain numpy on them
(np.array(sampler._rvs["distance"]) etc.), so cupy raised "Implicit conversion to
a NumPy array is not allowed". That throw happens inside analyze_event, so ILE
caught it, printed the MISLEADING generic hint "Probable reasons: SEOB nyquist
..." + "Skipping the following binary!", and skipped EVERY binary at the extrinsic
stage -> all_dslice.dat empty (and consolidate --allow-empty masked it). Intrinsic
iterations never run the export, so they were clean -- the symptom was "only the
extrinsic ILE barfed" on the same binaries.

Fix: convert cupy -> numpy at the boundary.
- ILE dslice block: read _rvs / pdf outputs through the module-level
  identity_convert (= cupy.asnumpy, or no-op without cupy): distance, log_weights,
  log_integrand, joint_prior/joint_s_prior (+ log_ variants), integrand, the
  sampler pdf/prior_pdf evaluations, and prior_pdf_d(d_all).
- distance_slices.py: add _to_cpu() helper; apply to every sampler._rvs read in
  _ln_omega_iw_factor and importance_reweight_slices, to the cached-likelihood
  return in importance_reweight_slices and fresh_sample_slices. Pure-numpy callers
  (the CPU path, unit tests) are unaffected.

Affects both core+wing and all-fresh. A CPU smoke test passes even while broken,
so this MUST be validated on GPU (done separately). Re-verified the CPU all-fresh
smoke test still writes a 6-row .dslice.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

Author: oshaughn

MonteCarloMarginalizeCode/Code/RIFT/misc/distance_slices.py

@@ -32,6 +32,21 @@
 import numpy as np
 
 
+def _to_cpu(x):
+    """Host (numpy) view of x, converting cupy arrays if needed.
+
+    On GPU ILE runs the sampler stores CUPY arrays in ``_rvs`` and the cached
+    likelihood returns cupy; numpy refuses implicit conversion of those
+    ("Implicit conversion to a NumPy array is not allowed. Please use .get()").
+    This collapses cupy -> numpy at the boundary. Plain numpy arrays / python
+    scalars pass straight through (cupy is not even imported here).
+    """
+    get = getattr(x, "get", None)
+    if get is not None and type(x).__module__.split(".")[0] == "cupy":
+        return get()
+    return x
+
+
 DISTANCE_SLICE_FIELDS = (
     "lnL",        # extrinsic-marginalized lnL at d=dist (pure likelihood,
                   # i.e. distance sampling prior divided out)
@@ -105,14 +120,14 @@ def _ln_omega_iw_factor(rvs, ln_prior_d_at_samples, ln_proposal_d_at_samples):
     """
     # Pull joint prior / proposal ratio
     if "joint_prior" in rvs and "joint_s_prior" in rvs:
-        jp = np.asarray(rvs["joint_prior"], float)
-        jsp = np.asarray(rvs["joint_s_prior"], float)
+        jp = np.asarray(_to_cpu(rvs["joint_prior"]), float)
+        jsp = np.asarray(_to_cpu(rvs["joint_s_prior"]), float)
         with np.errstate(divide="ignore"):
             ln_pi_over_q_joint = np.log(np.maximum(jp, np.finfo(float).tiny)) \
                                  - np.log(np.maximum(jsp, np.finfo(float).tiny))
     elif "log_joint_prior" in rvs and "log_joint_s_prior" in rvs:
-        ln_pi_over_q_joint = np.asarray(rvs["log_joint_prior"], float) \
-                             - np.asarray(rvs["log_joint_s_prior"], float)
+        ln_pi_over_q_joint = np.asarray(_to_cpu(rvs["log_joint_prior"]), float) \
+                             - np.asarray(_to_cpu(rvs["log_joint_s_prior"]), float)
     else:
         raise KeyError("sampler._rvs missing joint prior/proposal columns")
     return ln_pi_over_q_joint - (np.asarray(ln_prior_d_at_samples, float)
@@ -159,7 +174,7 @@ def importance_reweight_slices(
         if a not in rvs:
             raise KeyError("sampler._rvs missing required column {!r} for "
                             "slice reweighting".format(a))
-        fixed_inputs[a] = np.asarray(rvs[a])
+        fixed_inputs[a] = np.asarray(_to_cpu(rvs[a]))
 
     K = len(d_slices)
     lnL_out = np.empty(K)
@@ -173,7 +188,7 @@ def importance_reweight_slices(
             else:
                 like_inputs.append(fixed_inputs[a])
         lnL_at = like_to_integrate(*like_inputs)
-        lnL_at = np.asarray(lnL_at, dtype=np.float64)
+        lnL_at = np.asarray(_to_cpu(lnL_at), dtype=np.float64)
         if not return_lnL:
             # function returned exp(lnL - overflow); take log
             with np.errstate(divide="ignore"):
@@ -435,7 +450,9 @@ def like_at_pinned_d(**kw):
                 eps = 1e-12 * max(abs(hi - lo), 1.0)
                 full[p] = np.clip(np.asarray(arr, float), lo + eps, hi - eps)
             full["distance"] = d_arr
-            return like_to_integrate(*(full[a] for a in arg_names))
+            # like_to_integrate is the cached ILE likelihood -> returns CUPY on a
+            # GPU run; the fresh AV integrator here is host-side, so bring it back.
+            return _to_cpu(like_to_integrate(*(full[a] for a in arg_names)))
 
         try:
             res = sampler.integrate_log(

MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode

@@ -2632,11 +2632,16 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts,inv_spec_t
           if opts.sampler_method == "GMM" and neff < 50:
             print("  WARNING: --export-distance-slices with --sampler-method GMM at main n_eff={:.1f} (<50). ".format(neff)
                   + "B2-reweight may be biased; prefer --sampler-method AV or raise --n-max.")
-          dL_samp = np.array(sampler._rvs["distance"])
+          # On GPU the sampler stores CUPY arrays in _rvs / pdf outputs; convert
+          # to numpy at the boundary (identity_convert = cupy.asnumpy, or a no-op
+          # without cupy) before ANY numpy math. Otherwise cupy raises "Implicit
+          # conversion to a NumPy array is not allowed", analyze_event throws, and
+          # EVERY binary is skipped -> empty .dslice (silent failure).
+          dL_samp = np.asarray(identity_convert(sampler._rvs["distance"]), float)
           # ln(pi_d) at samples uses the actual sampler prior (volumetric or
           # pseudo-cosmo, whichever was registered).
           prior_pdf_d = sampler.prior_pdf["distance"]
-          pi_d_samp = np.asarray(prior_pdf_d(dL_samp), float)
+          pi_d_samp = np.asarray(identity_convert(prior_pdf_d(dL_samp)), float)
           pi_d_samp = np.where(pi_d_samp > 0, pi_d_samp, np.finfo(float).tiny)
           ln_pi_d_samp = np.log(pi_d_samp)
           # ln(q_d) at samples; for the standard ILE path the proposal is the
@@ -2647,7 +2652,7 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts,inv_spec_t
           # and q_d at the SAMPLES to isolate the Omega-only factor. Compute
           # q_d as a normalized 1-D density on the supported range.
           try:
-            q_d_raw = np.asarray(sampler.pdf["distance"](dL_samp), float)
+            q_d_raw = np.asarray(identity_convert(sampler.pdf["distance"](dL_samp)), float)
           except Exception:
             q_d_raw = np.ones_like(dL_samp)
           q_d_norm = float(getattr(sampler, "_pdf_norm", {}).get("distance", 1.0)) or 1.0
@@ -2658,13 +2663,14 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts,inv_spec_t
           # log-importance weights with the same fallback chain we use for
           # .dgrid.
           rvs = sampler._rvs
+          _rvs = lambda k: np.asarray(identity_convert(rvs[k]), float)  # cupy-safe column read
           if 'log_weights' in rvs:
-            ln_w_full = np.array(rvs['log_weights'])
+            ln_w_full = _rvs('log_weights')
           elif 'log_integrand' in rvs:
-            ln_w_full = np.array(rvs['log_integrand'] + rvs['log_joint_prior'] - rvs['log_joint_s_prior'])
+            ln_w_full = _rvs('log_integrand') + _rvs('log_joint_prior') - _rvs('log_joint_s_prior')
           else:
-            integrand = np.asarray(rvs['integrand'])
-            jp = np.asarray(rvs['joint_prior']); jsp = np.asarray(rvs['joint_s_prior'])
+            integrand = _rvs('integrand')
+            jp = _rvs('joint_prior'); jsp = _rvs('joint_s_prior')
             keep = (integrand > 0) & (jp > 0) & (jsp > 0)
             ln_w_full = np.full(len(integrand), -np.inf)
             ln_w_full[keep] = np.log(integrand[keep]) + np.log(jp[keep]) - np.log(jsp[keep])
@@ -2760,7 +2766,7 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts,inv_spec_t
               np.full(len(d_core), distance_slices.METHOD_REWEIGHT, dtype=int),
               np.full(len(d_wings), distance_slices.METHOD_FRESH, dtype=int),
           ])
-          ln_pi_d_all = np.log(np.maximum(prior_pdf_d(d_all), np.finfo(float).tiny))
+          ln_pi_d_all = np.log(np.maximum(np.asarray(identity_convert(prior_pdf_d(d_all)), float), np.finfo(float).tiny))
           order = np.argsort(d_all)
 
           try: