Add Bilby-compatible noise evidence output

MonteCarloMarginalizeCode/Code/RIFT/likelihood/noise_evidence.py

@@ -0,0 +1,61 @@
+"""Bilby-compatible fixed-PSD noise-evidence utilities.
+
+For the Gaussian transient likelihood used by Bilby and RIFT, the noise
+hypothesis has no sampled parameters.  Its log evidence is therefore the
+zero-signal log likelihood, ``-0.5 * sum_k (d_k | d_k)``.  The Gaussian
+determinant term is omitted, matching Bilby's GW likelihood convention.
+"""
+
+import math
+
+
+def compute_network_log_noise_evidence(
+        data_dict, psd_dict, fmin, fmax, fnyq,
+        inv_spec_trunc_Q=False, T_spec=0.0, inner_product_factory=None):
+    """Compute ``-0.5 * sum_k (d_k | d_k)`` for a detector network.
+
+    Parameters use the same conditioned data, PSDs, and frequency bounds as
+    ILE.  ``inner_product_factory`` is injectable to keep the bookkeeping
+    independently testable; production calls use :class:`RIFT.lalsimutils.ComplexIP`.
+
+    Returns
+    -------
+    total : float
+        Network log noise evidence.
+    per_detector : dict
+        Per-detector ``d_inner_d`` and ``log_noise_evidence`` values.
+    """
+    data_detectors = set(data_dict)
+    psd_detectors = set(psd_dict)
+    if not data_detectors:
+        raise ValueError("Cannot compute noise evidence without detector data")
+    if data_detectors != psd_detectors:
+        raise ValueError(
+            "Detector mismatch between data ({}) and PSDs ({})".format(
+                sorted(data_detectors), sorted(psd_detectors)))
+
+    if inner_product_factory is None:
+        from RIFT.lalsimutils import ComplexIP
+        inner_product_factory = ComplexIP
+
+    per_detector = {}
+    total = 0.0
+    for detector in sorted(data_detectors):
+        data = data_dict[detector]
+        inner_product = inner_product_factory(
+            fLow=fmin, fMax=fmax, fNyq=fnyq, deltaF=data.deltaF,
+            psd=psd_dict[detector], analyticPSD_Q=False,
+            inv_spec_trunc_Q=inv_spec_trunc_Q, T_spec=T_spec)
+        d_inner_d = float(inner_product.ip(data, data).real)
+        if not math.isfinite(d_inner_d) or d_inner_d < 0:
+            raise ValueError(
+                "Invalid (d|d)={} for detector {}".format(
+                    d_inner_d, detector))
+        detector_log_evidence = -0.5 * d_inner_d
+        per_detector[detector] = {
+            "d_inner_d": d_inner_d,
+            "log_noise_evidence": detector_log_evidence,
+        }
+        total += detector_log_evidence
+
+    return total, per_detector

MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode

@@ -8,6 +8,7 @@ Integrate the extrinsic parameters of the prefactored likelihood function.
 """
 
 import sys
+import json
 import functools
 from optparse import OptionParser, OptionGroup
 
@@ -48,6 +49,7 @@ from igwn_ligolw import utils,  ligolw
 import glue.lal
 
 import RIFT.lalsimutils as lalsimutils
+from RIFT.likelihood.noise_evidence import compute_network_log_noise_evidence
 import RIFT.integrators.mcsampler as mcsampler
 import RIFT.misc.sky_rotations as sky_rotations
 try:
@@ -177,6 +179,8 @@ optp.add_option( "--zero-likelihood", action='store_true', help="Run with exactl
 optp.add_option("-c", "--cache-file", default=None, help="LIGO cache file containing all data needed.")
 optp.add_option("-C", "--channel-name", action="append", help="instrument=channel-name, e.g. H1=FAKE-STRAIN. Can be given multiple times for different instruments.")
 optp.add_option("-p", "--psd-file", action="append", help="instrument=psd-file, e.g. H1=H1_PSD.xml.gz. Can be given multiple times for different instruments.")
+optp.add_option("--log-noise-evidence-output", default=None, help="Write Bilby-compatible fixed-PSD log noise evidence and per-detector (d|d) values to this JSON file. Uses the same conditioned data and PSD as ILE.")
+optp.add_option("--log-noise-evidence-only", action="store_true", help="Exit after writing --log-noise-evidence-output, before waveform generation and Monte Carlo integration. If no output path is supplied, use log_noise_evidence.json.")
 optp.add_option("-k", "--skymap-file", help="Use skymap stored in given FITS file.")
 optp.add_option("-x", "--coinc-xml", help="gstlal_inspiral XML file containing coincidence information.")
 optp.add_option("-I", "--sim-xml", help="XML file containing mass grid to be evaluated")
@@ -332,6 +336,8 @@ for pin_param in LIKELIHOOD_PINNABLE_PARAMS:
 optp.add_option_group(pinnable)
 
 opts, args = optp.parse_args()
+if opts.log_noise_evidence_only and opts.log_noise_evidence_output is None:
+  opts.log_noise_evidence_output = "log_noise_evidence.json"
 
 # cosmo d prior tools for interpolation:  not used normally, but set if needed
 final_scipy_interpolate=None
@@ -733,6 +739,37 @@ for Psig in P_list:
   Psig.deltaf = P.deltaF
 
 
+# Bilby's fixed-PSD noise hypothesis has no sampled parameters: its evidence is
+# the zero-signal likelihood, -1/2 sum_k (d_k|d_k).  Compute it only after all
+# ILE data/PSD conditioning so the normalization exactly matches this run.
+if opts.log_noise_evidence_output:
+  log_noise_evidence, noise_evidence_by_detector = compute_network_log_noise_evidence(
+      data_dict, psd_dict, fmin=P.fmin, fmax=fmax, fnyq=0.5/P.deltaT,
+      inv_spec_trunc_Q=inv_spec_trunc_Q, T_spec=T_spec)
+  noise_evidence_payload = {
+      "schema": "rift.log-noise-evidence/v1",
+      "convention": "bilby-fixed-psd-no-gaussian-determinant",
+      "log_noise_evidence": log_noise_evidence,
+      "per_detector": noise_evidence_by_detector,
+      "analysis": {
+          "fmin": float(P.fmin),
+          "fmax": float(fmax),
+          "fnyq": float(0.5/P.deltaT),
+          "inverse_spectrum_truncation_time": float(T_spec),
+          "data_start_time": float(start_time),
+          "data_end_time": float(end_time),
+      },
+  }
+  with open(opts.log_noise_evidence_output, "w") as noise_evidence_file:
+    json.dump(noise_evidence_payload, noise_evidence_file,
+              indent=2, sort_keys=True)
+    noise_evidence_file.write("\n")
+  print("Wrote log noise evidence {} to {}".format(
+      log_noise_evidence, opts.log_noise_evidence_output))
+  if opts.log_noise_evidence_only:
+    sys.exit(0)
+
+
 
 #
 # Set up parameters and bounds

MonteCarloMarginalizeCode/Code/test/test_noise_evidence.py

@@ -0,0 +1,57 @@
+import unittest
+
+from RIFT.likelihood.noise_evidence import compute_network_log_noise_evidence
+
+
+class _Data:
+    def __init__(self, value, delta_f):
+        self.value = value
+        self.deltaF = delta_f
+
+
+class _InnerProduct:
+    calls = []
+
+    def __init__(self, **kwargs):
+        self.kwargs = kwargs
+        self.__class__.calls.append(kwargs)
+
+    def ip(self, first, second):
+        return complex(first.value * second.value * self.kwargs["psd"])
+
+
+class TestNoiseEvidence(unittest.TestCase):
+    def setUp(self):
+        _InnerProduct.calls = []
+
+    def test_network_sum_and_bilby_sign(self):
+        data = {"L1": _Data(2.0, 0.25), "H1": _Data(3.0, 0.25)}
+        psds = {"H1": 2.0, "L1": 4.0}
+
+        total, per_detector = compute_network_log_noise_evidence(
+            data, psds, fmin=20.0, fmax=1024.0, fnyq=2048.0,
+            inv_spec_trunc_Q=True, T_spec=8.0,
+            inner_product_factory=_InnerProduct)
+
+        self.assertEqual(per_detector["H1"]["d_inner_d"], 18.0)
+        self.assertEqual(per_detector["L1"]["d_inner_d"], 16.0)
+        self.assertEqual(total, -17.0)
+        self.assertEqual([call["psd"] for call in _InnerProduct.calls], [2.0, 4.0])
+        for call in _InnerProduct.calls:
+            self.assertEqual(call["fLow"], 20.0)
+            self.assertEqual(call["fMax"], 1024.0)
+            self.assertEqual(call["fNyq"], 2048.0)
+            self.assertEqual(call["deltaF"], 0.25)
+            self.assertTrue(call["inv_spec_trunc_Q"])
+            self.assertEqual(call["T_spec"], 8.0)
+
+    def test_detector_mismatch_is_rejected(self):
+        with self.assertRaisesRegex(ValueError, "Detector mismatch"):
+            compute_network_log_noise_evidence(
+                {"H1": _Data(1.0, 0.25)}, {"L1": 1.0},
+                fmin=20.0, fmax=1024.0, fnyq=2048.0,
+                inner_product_factory=_InnerProduct)
+
+
+if __name__ == "__main__":
+    unittest.main()