Merge fork rift_O4d updates

.github/workflows/ci.yml

@@ -107,6 +107,43 @@ jobs:
       - name: Run import check
         run: python .travis/test-all-mod.py
 
+  dependency-compat-check:
+    needs: install
+    runs-on: ubuntu-latest
+    # Focused upstream-library smoke test for scipy/h5py/numba/matplotlib plus
+    # RIFT's portable precision dtype. This complements the broader import and
+    # test-run jobs with a compact version matrix and small functional probes.
+    strategy:
+      fail-fast: false
+      matrix:
+        include:
+          - lane: legacy
+            python-version: '3.9'
+            numpy-pin: 'numpy==1.24.4'
+          - lane: modern
+            python-version: '3.12'
+            numpy-pin: 'numpy>=2.0,<3.0'
+    name: dependency-compat-check (${{ matrix.lane }} py${{ matrix.python-version }})
+    steps:
+      - uses: actions/checkout@v4
+      - uses: actions/setup-python@v5
+        with:
+          python-version: ${{ matrix.python-version }}
+          cache: 'pip'
+          cache-dependency-path: requirements.txt
+      - name: Enable symlink
+        run: sudo ln -sf $(which python3) /usr/bin/python
+      - name: Install dependencies
+        run: |
+          python -m pip install --upgrade pip --break-system-packages
+          python -m pip install -r requirements.txt --break-system-packages
+          # Pin numpy AFTER requirements.txt so it overrides the unpinned
+          # 'numpy' line in requirements.txt without changing the file.
+          python -m pip install '${{ matrix.numpy-pin }}' --break-system-packages
+          python -m pip install --editable . --break-system-packages
+      - name: Run dependency compatibility check
+        run: python .travis/test-dependency-compat.py
+
   sim-manager-check:
     needs: install
     runs-on: ubuntu-latest

.gitlab-ci.yml

@@ -92,6 +92,11 @@ import_check:
     - python -m pytest -q MonteCarloMarginalizeCode/Code/test/test_lisa_operational_synthetic.py
     - python -m pytest -q MonteCarloMarginalizeCode/Code/test/test_lisa_helper_contract.py
 
+dependency_compat_check:
+  stage: system tests
+  script:
+    - python .travis/test-dependency-compat.py
+
 sim_manager_check:
   stage: unit tests
   script:

.travis/test-dependency-compat.py

@@ -0,0 +1,112 @@
+#!/usr/bin/env python3
+"""Smoke-test upstream scientific-library compatibility for RIFT CI.
+
+This intentionally stays small: the full RIFT import sweep catches broad
+breakage, while this file makes the dependency surface in issue #17 explicit
+and gives CI logs a compact version matrix for scipy/h5py/numba/matplotlib.
+"""
+
+from __future__ import annotations
+
+import json
+import tempfile
+from pathlib import Path
+
+import h5py
+import matplotlib
+
+matplotlib.use("Agg")
+from matplotlib import pyplot as plt  # noqa: E402
+
+import numba
+import numpy as np
+import scipy
+from scipy import integrate, optimize, stats
+
+from RIFT.precision import RIFT_FLOAT_HIGH_PRECISION, RIFT_FLOAT_NAME, RiftFloat
+
+
+@numba.njit(cache=False)
+def _numba_quadratic(x):
+    return x * x + 2.0 * x + 1.0
+
+
+def _check_scipy() -> None:
+    integral, err = integrate.quad(lambda x: x * x, 0.0, 1.0)
+    if not np.isclose(integral, 1.0 / 3.0, atol=1e-12):
+        raise AssertionError(f"scipy.integrate returned {integral} +/- {err}")
+
+    root = optimize.brentq(lambda x: x * x - 2.0, 0.0, 2.0)
+    if not np.isclose(root, np.sqrt(2.0), atol=1e-12):
+        raise AssertionError(f"scipy.optimize returned root {root}")
+
+    cdf = stats.norm.cdf(0.0)
+    if not np.isclose(cdf, 0.5, atol=1e-15):
+        raise AssertionError(f"scipy.stats returned normal cdf {cdf}")
+
+
+def _check_h5py() -> None:
+    values = np.arange(6, dtype=np.float64).reshape(2, 3)
+    with tempfile.TemporaryDirectory() as tmpdir:
+        h5_path = Path(tmpdir) / "compat.h5"
+        with h5py.File(h5_path, "w") as h5_file:
+            h5_file.create_dataset("values", data=values)
+            h5_file.attrs["library"] = "h5py"
+        with h5py.File(h5_path, "r") as h5_file:
+            roundtrip = h5_file["values"][()]
+            marker = h5_file.attrs["library"]
+    if marker != "h5py" or not np.array_equal(values, roundtrip):
+        raise AssertionError("h5py dataset/attribute round trip failed")
+
+
+def _check_numba() -> None:
+    sample = np.asarray([0.0, 1.0, 2.0])
+    expected = np.asarray([1.0, 4.0, 9.0])
+    actual = _numba_quadratic(sample)
+    if not np.allclose(actual, expected):
+        raise AssertionError(f"numba compiled function returned {actual}")
+
+
+def _check_matplotlib() -> None:
+    with tempfile.TemporaryDirectory() as tmpdir:
+        png_path = Path(tmpdir) / "compat.png"
+        fig, ax = plt.subplots(figsize=(2, 2))
+        ax.plot([0, 1], [0, 1])
+        ax.set_xlabel("x")
+        ax.set_ylabel("y")
+        fig.tight_layout()
+        fig.savefig(png_path)
+        plt.close(fig)
+        if png_path.stat().st_size <= 0:
+            raise AssertionError("matplotlib produced an empty PNG")
+
+
+def _check_rift_precision() -> None:
+    dtype = np.dtype(RiftFloat)
+    if dtype.itemsize < np.dtype(np.float64).itemsize:
+        raise AssertionError(f"RiftFloat unexpectedly narrower than float64: {dtype}")
+    if RIFT_FLOAT_HIGH_PRECISION != (dtype.itemsize > np.dtype(np.float64).itemsize):
+        raise AssertionError("RIFT_FLOAT_HIGH_PRECISION does not match RiftFloat width")
+
+
+def main() -> None:
+    matrix = {
+        "numpy": np.__version__,
+        "scipy": scipy.__version__,
+        "h5py": h5py.__version__,
+        "numba": numba.__version__,
+        "matplotlib": matplotlib.__version__,
+        "rift_float": RIFT_FLOAT_NAME,
+        "rift_float_high_precision": RIFT_FLOAT_HIGH_PRECISION,
+    }
+    print(json.dumps(matrix, indent=2, sort_keys=True))
+
+    _check_rift_precision()
+    _check_scipy()
+    _check_h5py()
+    _check_numba()
+    _check_matplotlib()
+
+
+if __name__ == "__main__":
+    main()

MonteCarloMarginalizeCode/Code/ILE_data_handling_probe.py

@@ -210,7 +210,8 @@ def get_unpinned_params(opts, params):
 NR_template_group=None
 NR_template_param=None
 if opts.nr_group and opts.nr_param:
-    import NRWaveformCatalogManager3 as nrwf
+    from RIFT.physics._nrwf_loader import get_nrwf as _rift_get_nrwf
+    nrwf, _useNR = _rift_get_nrwf()    # prefers nrcatalog.compat_nrwf, falls back to NRWaveformCatalogManager3
     NR_template_group = opts.nr_group
     if nrwf.internal_ParametersAreExpressions[NR_template_group]:
         NR_template_param = eval(opts.nr_param)

MonteCarloMarginalizeCode/Code/RIFT/calmarg/adaptive.py

@@ -170,11 +170,23 @@ def cal_mc_error_from_components(comp, cal_log_weights=None, sample_log_weights=
     logw = np.zeros(n_cal) if cal_log_weights is None else np.asarray(cal_log_weights, dtype=float)
     lc = comp + logw[None, :]                       # log( w_c L_jc )
     lnL_marg = logsumexp(lc, axis=1)                # per-sample log sum_c w_c L_jc (norm cancels)
-    log_r = lc - lnL_marg[:, None]                  # responsibilities r_jc, sum_c r_jc = 1
+    # Guard total underflow: a sample whose likelihood is -inf across ALL cal draws has
+    # lnL_marg = -inf, so lc - lnL_marg = -inf - -inf = NaN, which would poison log_a (and
+    # hence a_c / sigma / neff, and the adaptive-doubling stop test) for the whole batch.
+    # Such samples carry zero posterior weight, so force their responsibilities to -inf
+    # (contribute nothing) instead of NaN.  A pathological extrinsic point (e.g. a sky
+    # position with ~zero antenna response) can produce this on a real run.
+    finite = np.isfinite(lnL_marg)
+    log_r = np.full_like(lc, -np.inf)               # responsibilities r_jc, sum_c r_jc = 1
+    if np.any(finite):
+        log_r[finite] = lc[finite] - lnL_marg[finite, None]
     if sample_log_weights is None:
-        slw = lnL_marg                              # prior-drawn batch -> weight by marginal L
+        slw = np.where(finite, lnL_marg, -np.inf)   # prior-drawn batch -> weight by marginal L
     else:
         slw = np.asarray(sample_log_weights, dtype=float)
+    slw = np.where(np.isfinite(slw), slw, -np.inf)
+    if not np.any(np.isfinite(slw)):                # fully-degenerate batch: nothing to weight
+        return float('inf'), 1.0, np.full(n_cal, 1.0 / n_cal)
     slw = slw - logsumexp(slw)                      # sum_j W_j = 1
     log_a = logsumexp(slw[:, None] + log_r, axis=0) # a_c = sum_j W_j r_jc
     a_c = np.exp(log_a - logsumexp(log_a))          # exact renormalization

MonteCarloMarginalizeCode/Code/RIFT/calmarg/breadcrumbs.py

@@ -54,14 +54,14 @@ def save(path, cal=None, extrinsic=None, kind="gaussian", meta=None):
             cal_prior_sigma=np.asarray(cal["prior_sigma"], dtype=float),
             cal_node_log_f=np.asarray(cal["node_log_f"], dtype=float),
             cal_n_nodes_amp=np.int64(cal["n_nodes_amp"]),
-            cal_dets=np.array(list(cal["dets"]), dtype=object),
+            cal_dets=np.array([str(x) for x in cal["dets"]]),   # string dtype (NOT object): no
         )
     if extrinsic is not None:
         groups = extrinsic["groups"]
         d["ext_kind"] = str(extrinsic.get("kind", "gmm"))
         d["ext_n_groups"] = np.int64(len(groups))
         for i, g in enumerate(groups):
-            d["ext_g%d_params" % i] = np.array(list(g["params"]), dtype=object)
+            d["ext_g%d_params" % i] = np.array([str(x) for x in g["params"]])   # string, not object
             d["ext_g%d_means" % i] = np.asarray(g["means"], dtype=float)
             d["ext_g%d_covs" % i] = np.asarray(g["covariances"], dtype=float)
             d["ext_g%d_weights" % i] = np.asarray(g["weights"], dtype=float)
@@ -113,6 +113,12 @@ def load(path):
     assert g["kind"] == "gaussian" and g["cal"]["dets"] == ["H1", "L1", "V1"]
     assert np.allclose(g["cal"]["proposal_mean"], cal["proposal_mean"])
     assert g["meta"]["iteration"] == 2
+    # PORTABILITY: the file must contain NO pickled objects, so a breadcrumb written by
+    # one numpy (e.g. the container's 2.x) loads under any other (e.g. the host's 1.x).
+    # np.load(allow_pickle=False) raises if any array is object-dtype -- guards the dets/
+    # params regression where dtype=object silently pickled and broke cross-version load.
+    with np.load(p, allow_pickle=False) as _z:
+        assert _z["cal_dets"].dtype.kind in ("U", "S"), _z["cal_dets"].dtype
 
     # extrinsic (GMM) round-trip
     ext = dict(kind="gmm", groups=[
@@ -133,4 +139,7 @@ def load(path):
     assert np.allclose(g2["extrinsic"]["groups"][0]["means"], ext["groups"][0]["means"])
     assert np.allclose(g2["extrinsic"]["groups"][1]["covariances"], ext["groups"][1]["covariances"])
     assert g2["cal"] is not None   # cal + extrinsic coexist in one breadcrumb
-    print("PASS: breadcrumb save/load round-trips (cal Gaussian + extrinsic GMM, schema v%d)." % SCHEMA_VERSION)
+    with np.load(p2, allow_pickle=False) as _z2:   # portability: no pickled object arrays
+        assert _z2["ext_g0_params"].dtype.kind in ("U", "S"), _z2["ext_g0_params"].dtype
+    print("PASS: breadcrumb save/load round-trips (cal Gaussian + extrinsic GMM, schema v%d), "
+          "pickle-free (portable across numpy versions)." % SCHEMA_VERSION)

MonteCarloMarginalizeCode/Code/RIFT/calmarg/test_cal_mc_error.py

@@ -94,10 +94,31 @@ def test_extrinsic_batch_weighting():
     print("test_extrinsic_batch_weighting: OK")
 
 
+def test_total_underflow_guard():
+    # A sample with -inf lnL across ALL cal draws (e.g. a zero-response extrinsic point)
+    # must not poison the batch with NaN; it carries zero weight, so the answer must equal
+    # the finite-samples-only answer.
+    rng = np.random.default_rng(17)
+    n_cal = 200
+    good = rng.normal(0.0, 0.8, size=(8, n_cal))
+    s_good, neff_good, a_good = cal_mc_error_from_components(good)
+    comp = np.vstack([good, np.full((1, n_cal), -np.inf)])   # append a dead sample
+    s, neff, a = cal_mc_error_from_components(comp)
+    assert np.isfinite(s) and np.isfinite(neff), (s, neff)
+    assert abs(a.sum() - 1.0) < 1e-12 and np.all(np.isfinite(a))
+    assert abs(s - s_good) < 1e-9 and abs(neff - neff_good) < 1e-7, (s, s_good, neff, neff_good)
+    # fully-degenerate batch (every sample dead) returns a finite sentinel, not NaN
+    dead = np.full((3, n_cal), -np.inf)
+    s0, neff0, a0 = cal_mc_error_from_components(dead)
+    assert np.isinf(s0) and neff0 == 1.0 and abs(a0.sum() - 1.0) < 1e-12
+    print("test_total_underflow_guard: OK")
+
+
 if __name__ == "__main__":
     test_lognormal_closed_form()
     test_brute_force_scatter()
     test_neff_dominated()
     test_importance_weights_consistency()
     test_extrinsic_batch_weighting()
+    test_total_underflow_guard()
     print("ALL OK")

MonteCarloMarginalizeCode/Code/RIFT/calmarg/test_seed_fallback.py

@@ -0,0 +1,123 @@
+#!/usr/bin/env python3
+"""
+Regression tests for the calmarg PILOT "graceful degradation" contract.
+
+Background
+----------
+In the in-loop-calmarg PILOT pipeline the wide ILE jobs (and the last-iteration
+EXTRINSIC ILE jobs) are SEEDED from the previous iteration's consolidated cal
+proposal, referenced as  cal_consolidated_$(macroiterationprev).npz  and (on OSG
+file transfer) listed in condor transfer_input_files.  A calpilot only produces
+that seed for iterations it<=--calmarg-pilot-max-it on-cadence, so a wide/extrinsic
+iteration whose seed was never produced used to HARD-HOLD condor on a missing
+transfer source (HoldReasonCode 13) and deadlock the whole DAG.
+
+The fix pre-seeds, for EVERY referenced iteration index, a placeholder copy of the
+always-present  cal_consolidated_-1.npz  iteration-0 seed -- which is a VALID prior
+breadcrumb (proposal == prior).  A real calpilot overwrites its placeholder at run
+time, so behavior is unchanged whenever the learned seed IS produced; a missing seed
+now degrades to the PRIOR instead of dead-holding.  This module locks the two
+properties that make that safe:
+
+  1. A prior breadcrumb (the placeholder content) round-trips and, when seeded from,
+     yields UNWEIGHTED prior cal draws (importance log-weights == 0).  I.e. "fall back
+     to the placeholder" is exactly equivalent to "draw from the broad prior".
+  2. The absent/empty-seed guard the ILE binary uses --
+     `os.path.exists(p) and os.path.getsize(p) > 0` -- correctly classifies a missing
+     path and a 0-byte placeholder as "not present yet" (-> prior fallback), and
+     breadcrumbs.load() raises on an empty file (so the ILE except-fallback fires).
+
+Run:  python3 -m RIFT.calmarg.test_seed_fallback
+"""
+import os
+import tempfile
+
+import numpy as np
+
+import RIFT.calmarg.breadcrumbs as breadcrumbs
+import RIFT.calmarg.generate_realizations as genr
+
+
+def _prior_breadcrumb_cal(n_amp=5, dets=("H1", "L1"), fmin=20.0, fmax=500.0,
+                          prior_sigma_val=0.1):
+    """A 'placeholder' cal breadcrumb: proposal == prior (mean, diagonal cov).
+
+    This is the semantic content of cal_consolidated_-1.npz and of the copies the
+    pipeline pre-seeds into iterations a calpilot never produced.
+    """
+    dim = 2 * n_amp * len(dets)
+    prior_mean = np.zeros(dim)
+    prior_sigma = np.full(dim, prior_sigma_val)
+    return dict(
+        proposal_mean=prior_mean.copy(),
+        proposal_cov=np.diag(prior_sigma ** 2),   # proposal == prior
+        prior_mean=prior_mean,
+        prior_sigma=prior_sigma,
+        node_log_f=np.linspace(np.log10(fmin), np.log10(fmax), n_amp),
+        n_nodes_amp=n_amp,
+        dets=list(dets),
+    )
+
+
+def test_prior_placeholder_seeds_as_prior():
+    """Seeding from a prior placeholder == unweighted prior draws (log_w == 0)."""
+    n_amp = 5
+    dets = ("H1", "L1")
+    fmin, fmax = 20.0, 500.0
+    cal = _prior_breadcrumb_cal(n_amp=n_amp, dets=dets, fmin=fmin, fmax=fmax)
+
+    p = os.path.join(tempfile.mkdtemp(), "cal_consolidated_-1.npz")
+    breadcrumbs.save(p, cal=cal, meta=dict(placeholder=True, iteration=-1))
+    bc = breadcrumbs.load(p)
+
+    _dat, cal_log_weights, nodes = genr.seed_realizations_from_breadcrumb(
+        bc, T_segment=4.0, dT=1.0 / 1024, fmin=fmin, fmax=fmax,
+        n_spline_points=n_amp, n_realizations=256,
+        rng=np.random.default_rng(1234))
+
+    assert nodes.shape == (256, 2 * n_amp * len(dets)), nodes.shape
+    # proposal == prior  =>  log(prior/proposal) == 0 for every realization.
+    assert np.allclose(cal_log_weights, 0.0, atol=1e-8), \
+        "prior-placeholder importance weights should be exactly 0; max|w|=%g" \
+        % np.max(np.abs(cal_log_weights))
+    print("test_prior_placeholder_seeds_as_prior: OK "
+          "(max|log_w|=%.2e over %d realizations)"
+          % (np.max(np.abs(cal_log_weights)), len(cal_log_weights)))
+
+
+def _seed_present(path):
+    """Exact guard used by the ILE binary to decide 'seed present vs fall back'."""
+    return os.path.exists(path) and os.path.getsize(path) > 0
+
+
+def test_absent_and_empty_seed_guard():
+    """A missing path and a 0-byte placeholder both classify as 'not present'."""
+    d = tempfile.mkdtemp()
+
+    missing = os.path.join(d, "cal_consolidated_3.npz")   # never produced
+    assert not _seed_present(missing), "missing seed must classify as absent"
+
+    empty = os.path.join(d, "cal_consolidated_.npz")      # 0-byte placeholder
+    open(empty, "a").close()
+    assert not _seed_present(empty), "0-byte seed must classify as absent"
+    # ILE loads inside try/except; an empty .npz must raise so the except-fallback fires.
+    raised = False
+    try:
+        breadcrumbs.load(empty)
+    except Exception:
+        raised = True
+    assert raised, "breadcrumbs.load() must raise on an empty (0-byte) placeholder"
+
+    real = os.path.join(d, "cal_consolidated_0.npz")       # a genuine seed
+    breadcrumbs.save(real, cal=_prior_breadcrumb_cal(), meta=dict(iteration=0))
+    assert _seed_present(real), "a genuine breadcrumb must classify as present"
+    breadcrumbs.load(real)   # must load without raising
+    print("test_absent_and_empty_seed_guard: OK "
+          "(missing -> absent, 0-byte -> absent+raises, real -> present)")
+
+
+if __name__ == "__main__":
+    test_prior_placeholder_seeds_as_prior()
+    test_absent_and_empty_seed_guard()
+    print("ALL OK: calmarg pilot seed absent/placeholder -> graceful prior fallback "
+          "(no transfer hard-hold).")