From cc15a306bfe5083a1dfacee93e0c017208d90225 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Fri, 28 Aug 2026 13:48:34 -0400 Subject: [PATCH 01/16] EOSManager: support reviewed lalsim multibranch families --- CHANGES.rst | 6 + .../Code/RIFT/physics/EOSManager.py | 175 ++++++++++++---- .../Code/RIFT/physics/lalsim_eos_compat.py | 183 +++++++++++++++++ .../Code/test/test_lalsim_eos_compat.py | 186 ++++++++++++++++++ 4 files changed, 517 insertions(+), 33 deletions(-) create mode 100644 MonteCarloMarginalizeCode/Code/RIFT/physics/lalsim_eos_compat.py create mode 100644 MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py diff --git a/CHANGES.rst b/CHANGES.rst index c824160a7..afd21236e 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,3 +1,9 @@ +0.0.17.13 compatibility backstop +-------------------------------- + - EOSManager supports released and reviewed multipart/multibranch LALSimulation + family interfaces while retaining the existing scalar and NuclearMatter-Backend + sequence contracts. See ``docs/eos-interface-contract.md``. + 0.0.17.13 --------- MR https://git.ligo.org/rapidpe-rift/rift/-/merge_requests/55 , for ln(e) parameter access in pipeline diff --git a/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py b/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py index 5a4123d19..d131b5be2 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py @@ -21,6 +21,11 @@ import scipy.interpolate as interp import scipy +from RIFT.physics.lalsim_eos_compat import ( + LALSimNeutronStarFamilyAdapter, + create_family, +) + try: from natsort import natsorted except: @@ -72,16 +77,46 @@ def __init__(self,name=None): self.eos_fam = None return None - def lambda_from_m(self, m): - eos_fam = self.eos_fam + def _set_lalsim_family(self, minimal=True): + """Create and cache a released-or-multibranch LAL family.""" + self._lalsim_family_adapter = create_family( + self.eos, minimal=minimal, lalsim_module=lalsim + ) + self.eos_fam = self._lalsim_family_adapter.family + self.mMaxMsun = self._lalsim_family_adapter.maximum_mass() / lal.MSUN_SI + return self.eos_fam + + def _get_lalsim_family_adapter(self): + adapter = getattr(self, "_lalsim_family_adapter", None) + if adapter is None or adapter.family is not self.eos_fam: + adapter = LALSimNeutronStarFamilyAdapter.from_family( + self.eos_fam, lalsim_module=lalsim + ) + self._lalsim_family_adapter = adapter + return adapter + + def branches_for_m(self, m): + """Return stable LAL family branches available at mass ``m``.""" + if m < 10**15: + m = m * lal.MSUN_SI + return self._get_lalsim_family_adapter().branches_for_mass(m) + + def radius_from_m(self, m, branch_id=None): + """Return radius in metres; require ``branch_id`` for twin stars.""" + if m < 10**15: + m = m * lal.MSUN_SI + return self._get_lalsim_family_adapter().radius(m, branch_id=branch_id) + + def lambda_from_m(self, m, branch_id=None): if m<10**15: m=m*lal.MSUN_SI if m/lal.MSUN_SI > 0.999*self.mMaxMsun: return 1e-8 # not exactly zero, but tiny - k2=lalsim.SimNeutronStarLoveNumberK2(m, eos_fam) - r=lalsim.SimNeutronStarRadius(m, eos_fam) + family = self._get_lalsim_family_adapter() + k2 = family.love_number_k2(m, branch_id=branch_id) + r = family.radius(m, branch_id=branch_id) m=m*lal.G_SI/lal.C_SI**2 lam=2./(3*lal.G_SI)*k2*r**5 @@ -89,6 +124,18 @@ def lambda_from_m(self, m): return dimensionless_lam + def lambda_from_m_vector(self, m, branch_id=None): + # vectorize above. Naive vectorization, will depend on improved swig interface + # alternative is to store meta-information and re-interpolatie + if not(isinstance(m, np.ndarray)): + return self.lambda_from_m(m, branch_id=branch_id) + else: + # horribly inefficient right now! Hopefully vectorized code available soon + out = np.array([ + self.lambda_from_m(m_here, branch_id=branch_id) for m_here in m + ]) + return out + def estimate_baryon_mass_from_mg(self,m): r""" Estimate m_b = m_g + m_g^2/(R_{1.4}/km) based on https://arxiv.org/pdf/1905.03784.pdf Eq. (6) @@ -97,7 +144,7 @@ def estimate_baryon_mass_from_mg(self,m): but lalsuite doesn't provide access to this low-level info !! This function is only for use when LALEOS is created. Use RePrimAnd's baryon_mass_from_mg preferably for most other purposes!! """ - r1p4 =lalsim.SimNeutronStarRadius(1.4*lal.MSUN_SI, self.eos_fam)/1e3 + r1p4 = self.radius_from_m(1.4, branch_id=None) / 1e3 return m + (1./r1p4)*m**2 def pressure_density_on_grid_alternate(self,logrho_grid,enforce_causal=False): @@ -188,7 +235,7 @@ def test_speed_of_sound_causal(self, test_only_under_mmax=True,fast_test=True): hmax = lalsim.SimNeutronStarEOSMaxPseudoEnthalpy(eos) else: try: - pmax = lalsim.SimNeutronStarCentralPressure(m_max_SI,fam) + pmax = self._get_lalsim_family_adapter().central_pressure(m_max_SI) hmax = lalsim.SimNeutronStarEOSPseudoEnthalpyOfPressure(pmax,eos) except: # gatch gsl interpolation errors for example @@ -222,11 +269,43 @@ def __init__(self,name): self.name=name self.eos = lalsim.SimNeutronStarEOSByName(name) - self.eos_fam = lalsim.CreateSimNeutronStarFamily(self.eos) - self.mMaxMsun = lalsim.SimNeutronStarMaximumMass(self.eos_fam) / lal.MSUN_SI + self._set_lalsim_family() return None +class EOSLALSimulationFromFile(EOSConcrete): + """Load a released two-column or reviewed nine-column LAL EOS table. + + The reviewed LALSimulation reader detects clean phase transitions in both + formats and preserves all thermodynamic columns in the new format. Set + ``dirty_phase_transitions`` to request its opt-in correction of numerically + imperfect pressure plateaus. + """ + + def __init__(self, fname, name=None, dirty_phase_transitions=False, + skip_family=False, minimal_family=True): + self.name = name or os.path.basename(fname) + self.fname = fname + self.eos = None + self.eos_fam = None + dirty_reader = getattr( + lalsim, "SimNeutronStarEOSFromFileChoiceDirtyPT", None + ) + if dirty_phase_transitions: + if dirty_reader is None: + raise NotImplementedError( + "dirty phase-transition correction requires the reviewed " + "LALSimulation multipart EOS interface" + ) + self.eos = dirty_reader(fname, 1) + else: + self.eos = lalsim.SimNeutronStarEOSFromFile(fname) + if not skip_family: + self._set_lalsim_family(minimal=minimal_family) + else: + self.mMaxMsun = None + + ### ### SERVICE 2: EOSFromFile @@ -325,7 +404,7 @@ def __init__(self,name=None,eos_data=None,eos_units=None,reject_phase_transition self.eos = lalsim.SimNeutronStarEOSFromFile(eos_fname) self.eos_fam = None if not(skip_family): - self.eos_fam = lalsim.CreateSimNeutronStarFamily(self.eos) + self._set_lalsim_family() return None @@ -380,13 +459,15 @@ def eos_ls(self): eos_fname = "./" +eos_name + "_geom.dat" # assume write acces np.savetxt(eos_fname, np.transpose((press, edens)), delimiter='\t') eos = lalsim.SimNeutronStarEOSFromFile(eos_fname) - fam = lalsim.CreateSimNeutronStarFamily(eos) + family_adapter = create_family(eos) + fam = family_adapter.family else: print(" No such file ", self.fname) sys.exit(0) - self.mMaxMsun = lalsim.SimNeutronStarMaximumMass(fam) / lal.MSUN_SI + self._lalsim_family_adapter = family_adapter + self.mMaxMsun = family_adapter.maximum_mass() / lal.MSUN_SI return eos, fam def p_rho_arrays(self): @@ -486,8 +567,7 @@ def __init__(self,name,param_dict=None): self.mMaxMsun=None self.eos=lalsim.SimNeutronStarEOS4ParameterPiecewisePolytrope(param_dict['logP1'], param_dict['gamma1'], param_dict['gamma2'], param_dict['gamma3']) - self.eos_fam=lalsim.CreateSimNeutronStarFamily(self.eos) - self.mMaxMsun = lalsim.SimNeutronStarMaximumMass(self.eos_fam) / lal.MSUN_SI + self._set_lalsim_family() return None @@ -522,22 +602,16 @@ def __init__(self,name=None,spec_params=None,verbose=False,use_lal_spec_eos=Fals if check_cs_builtin: # this valid = self.test_speed_of_sound_causal_builtin() # call parent class method - self.eos_fam = lalsim.CreateSimNeutronStarFamily(self.eos) - mmass = lalsim.SimNeutronStarMaximumMass(self.eos_fam) / lal.MSUN_SI - self.mMaxMsun = mmass + self._set_lalsim_family() else: # this test requires these quantities to be built *first* - self.eos_fam = lalsim.CreateSimNeutronStarFamily(self.eos) - mmass = lalsim.SimNeutronStarMaximumMass(self.eos_fam) / lal.MSUN_SI - self.mMaxMsun = mmass + self._set_lalsim_family() valid = self.test_speed_of_sound_causal() # call parent class method if not valid: raise Exception(" EOS : spectral sound speed violates speed of light ") elif not(no_eos_fam): # must create these if not performing the test - self.eos_fam = lalsim.CreateSimNeutronStarFamily(self.eos) - mmass = lalsim.SimNeutronStarMaximumMass(self.eos_fam) / lal.MSUN_SI - self.mMaxMsun = mmass + self._set_lalsim_family() else: self.eos_fam=None self.mMaxMsun = None @@ -681,8 +755,7 @@ def __init__(self,name=None,spec_params=None,verbose=False,use_lal_spec_eos=True cwd = os.getcwd() self.eos=lalsim.SimNeutronStarEOSFromFile(cwd+"/"+name+"_geom.dat") if not(no_eos_fam): - self.eos_fam = lalsim.CreateSimNeutronStarFamily(self.eos) - self.mMaxMsun = lalsim.SimNeutronStarMaximumMass(self.eos_fam) / lal.MSUN_SI + self._set_lalsim_family() else: self.eos_fam=None self.mMaxMsun = None @@ -1044,27 +1117,61 @@ def int_func(x_prime): ### # Les-like -def make_mr_lambda_lal(eos,n_bins=100): +def make_mr_lambda_lal(eos, n_bins=100, branch_id=None): ''' Construct mass-radius curve from EOS Based on modern code resources (https://git.ligo.org/publications/gw170817/bns-eos/blob/master/scripts/eos-params.py) which access low-level structures + + ``branch_id`` is optional for released/single-branch LALSimulation. It is + required for a multibranch family so an overlapping twin-star interval is + never collapsed silently. ''' - fam=lalsim.CreateSimNeutronStarFamily(eos) - max_m = lalsim.SimNeutronStarMaximumMass(fam)/lal.MSUN_SI - min_m = lalsim.SimNeutronStarFamMinimumMass(fam)/lal.MSUN_SI + family = create_family(eos) + if family.number_of_branches > 1 and branch_id is None: + raise ValueError( + "multibranch LAL family requires branch_id; use " + "make_mr_lambda_lal_branches() to construct every stable branch" + ) + max_m = family.maximum_mass(branch_id=branch_id)/lal.MSUN_SI + min_m = family.minimum_mass(branch_id=branch_id)/lal.MSUN_SI mgrid = np.linspace(min_m,max_m, n_bins) mrL_dat = np.zeros((n_bins,3)) mrL_dat[:,0] = mgrid for indx in np.arange(n_bins): mass_now = mgrid[indx] - r = lalsim.SimNeutronStarRadius(mass_now*lal.MSUN_SI,fam)/1000. + r = family.radius(mass_now*lal.MSUN_SI, branch_id=branch_id)/1000. mrL_dat[indx,1] = r - k = lalsim.SimNeutronStarLoveNumberK2(mass_now*lal.MSUN_SI,fam) + k = family.love_number_k2(mass_now*lal.MSUN_SI, branch_id=branch_id) c = mass_now * lal.MRSUN_SI / (r*1000.) mrL_dat[indx,2] = (2. / 3.) * k / c**5. return mrL_dat + +def make_mr_lambda_lal_branches(eos, n_bins=100): + """Return ``{branch_id: [M, R, Lambda]}`` for every stable LAL branch.""" + family = create_family(eos) + return { + branch_id: _make_mr_lambda_for_family(family, n_bins, branch_id) + for branch_id in range(family.number_of_branches) + } + + +def _make_mr_lambda_for_family(family, n_bins, branch_id): + min_m = family.minimum_mass(branch_id=branch_id) / lal.MSUN_SI + max_m = family.maximum_mass(branch_id=branch_id) / lal.MSUN_SI + mgrid = np.linspace(min_m, max_m, n_bins) + result = np.zeros((n_bins, 3)) + result[:, 0] = mgrid + for indx, mass_now in enumerate(mgrid): + mass_si = mass_now * lal.MSUN_SI + radius_m = family.radius(mass_si, branch_id=branch_id) + k2 = family.love_number_k2(mass_si, branch_id=branch_id) + compactness = mass_now * lal.MRSUN_SI / radius_m + result[indx, 1] = radius_m / 1000.0 + result[indx, 2] = (2.0 / 3.0) * k2 / compactness**5 + return result + # Rizzo def make_mr_lambda(eos,use_lal=False): """ @@ -1074,7 +1181,8 @@ def make_mr_lambda(eos,use_lal=False): if use_lal: make_mr_lambda_lal(eos) - fam=lalsim.CreateSimNeutronStarFamily(eos) + family = create_family(eos) + fam = family.family r_cut = 40 # Some EOS we consider for PE purposes will have very large radius! @@ -1108,8 +1216,9 @@ def make_mr_lambda(eos,use_lal=False): # - very frustrating...this data is embedded in the C code fac_max=1.6 r_fin=20. - m_ref = lalsim.SimNeutronStarMaximumMass(fam)/lal.MSUN_SI - r_ref = lalsim.SimNeutronStarRadius(lalsim.SimNeutronStarMaximumMass(fam), fam)/(10**3) + m_ref_si = family.maximum_mass() + m_ref = m_ref_si/lal.MSUN_SI + r_ref = family.radius(m_ref_si)/(10**3) answer=None while r_fin > r_ref or r_fin < 7: #print "Trying min:" diff --git a/MonteCarloMarginalizeCode/Code/RIFT/physics/lalsim_eos_compat.py b/MonteCarloMarginalizeCode/Code/RIFT/physics/lalsim_eos_compat.py new file mode 100644 index 000000000..39d859e3a --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/physics/lalsim_eos_compat.py @@ -0,0 +1,183 @@ +"""Compatibility helpers for the evolving LALSimulation neutron-star API. + +Released LALSimulation versions expose a single-branch +``LALSimNeutronStarFamily`` interface. The reviewed TOV development adds +multipart equations of state, multiple stable branches, and branch-indexed +interpolators. This module keeps the version checks in one place and gives +EOSManager a scalar interface which fails explicitly when a mass has twin-star +solutions. +""" + + +class AmbiguousFamilyBranchError(ValueError): + """Raised when a mass belongs to more than one stable family branch.""" + + +class LALSimNeutronStarFamilyAdapter: + """Version-neutral access to a LALSimulation neutron-star family. + + Parameters + ---------- + eos: + A SWIG ``LALSimNeutronStarEOS`` object. + minimal: + On the reviewed API, request the fast family containing only mass, + radius, and k2. Released APIs do not have this argument and ignore it. + lalsim_module: + Dependency-injection hook used by the interface contract tests. + """ + + _MODERN_REQUIRED = ( + "SimNeutronStarFamNumberOfBranches", + "SimNeutronStarFamMinMassPerBranch", + "SimNeutronStarFamMaxMassPerBranch", + "SimNeutronStarFamRadiusOfMassPerBranch", + "SimNeutronStarFamLoveNumberK2OfMassPerBranch", + ) + + def __init__(self, eos, minimal=True, lalsim_module=None): + if lalsim_module is None: + import lalsimulation as lalsim_module + self.lalsim = lalsim_module + self.eos = eos + self.is_multibranch_api = all( + hasattr(self.lalsim, name) for name in self._MODERN_REQUIRED + ) + if self.is_multibranch_api: + # The reviewed API requires ``min_fam``: 1 selects the PE-oriented + # M/R/k2 solver, while 0 also constructs baryonic mass, k3, and k4. + self.family = self.lalsim.CreateSimNeutronStarFamily( + eos, int(bool(minimal)) + ) + else: + self.family = self.lalsim.CreateSimNeutronStarFamily(eos) + + @classmethod + def from_family(cls, family, lalsim_module=None): + """Wrap an already-created family (mainly for tests and transition code).""" + if lalsim_module is None: + import lalsimulation as lalsim_module + obj = cls.__new__(cls) + obj.lalsim = lalsim_module + obj.eos = None + obj.family = family + obj.is_multibranch_api = all( + hasattr(obj.lalsim, name) for name in cls._MODERN_REQUIRED + ) + return obj + + @property + def number_of_branches(self): + if not self.is_multibranch_api: + return 1 + return int(self.lalsim.SimNeutronStarFamNumberOfBranches(self.family)) + + def minimum_mass(self, branch_id=None): + if self.is_multibranch_api: + if branch_id is None: + fn = getattr(self.lalsim, "SimNeutronStarFamMinMass", None) + if fn is not None: + return fn(self.family) + return min(self.minimum_mass(k) for k in range(self.number_of_branches)) + self._validate_branch_id(branch_id) + return self.lalsim.SimNeutronStarFamMinMassPerBranch( + self.family, int(branch_id) + ) + self._validate_legacy_branch_id(branch_id) + return self.lalsim.SimNeutronStarFamMinimumMass(self.family) + + def maximum_mass(self, branch_id=None): + if self.is_multibranch_api: + if branch_id is None: + fn = getattr(self.lalsim, "SimNeutronStarFamMaxMass", None) + if fn is not None: + return fn(self.family) + return max(self.maximum_mass(k) for k in range(self.number_of_branches)) + self._validate_branch_id(branch_id) + return self.lalsim.SimNeutronStarFamMaxMassPerBranch( + self.family, int(branch_id) + ) + self._validate_legacy_branch_id(branch_id) + return self.lalsim.SimNeutronStarMaximumMass(self.family) + + def branches_for_mass(self, mass_si): + """Return every stable branch whose closed mass interval contains mass_si.""" + if not self.is_multibranch_api: + return [0] if self.minimum_mass() <= mass_si <= self.maximum_mass() else [] + return [ + branch_id + for branch_id in range(self.number_of_branches) + if self.minimum_mass(branch_id) <= mass_si <= self.maximum_mass(branch_id) + ] + + def resolve_branch(self, mass_si, branch_id=None): + candidates = self.branches_for_mass(mass_si) + if branch_id is not None: + branch_id = int(branch_id) + self._validate_branch_id(branch_id) + if branch_id not in candidates: + raise ValueError( + "mass {!r} kg is outside stable branch {} (available branches: {})".format( + mass_si, branch_id, candidates + ) + ) + return branch_id + if len(candidates) == 1: + return candidates[0] + if not candidates: + raise ValueError( + "mass {!r} kg is outside every stable neutron-star branch".format( + mass_si + ) + ) + raise AmbiguousFamilyBranchError( + "mass {!r} kg has twin-star solutions on branches {}; pass branch_id explicitly".format( + mass_si, candidates + ) + ) + + def radius(self, mass_si, branch_id=None): + resolved = self.resolve_branch(mass_si, branch_id) + if self.is_multibranch_api: + return self.lalsim.SimNeutronStarFamRadiusOfMassPerBranch( + mass_si, self.family, resolved + ) + return self.lalsim.SimNeutronStarRadius(mass_si, self.family) + + def love_number_k2(self, mass_si, branch_id=None): + resolved = self.resolve_branch(mass_si, branch_id) + if self.is_multibranch_api: + return self.lalsim.SimNeutronStarFamLoveNumberK2OfMassPerBranch( + mass_si, self.family, resolved + ) + return self.lalsim.SimNeutronStarLoveNumberK2(mass_si, self.family) + + def central_pressure(self, mass_si, branch_id=None): + resolved = self.resolve_branch(mass_si, branch_id) + modern = getattr( + self.lalsim, "SimNeutronStarFamCentralPressureOfMassPerBranch", None + ) + if self.is_multibranch_api and modern is not None: + return modern(mass_si, self.family, resolved) + return self.lalsim.SimNeutronStarCentralPressure(mass_si, self.family) + + def _validate_branch_id(self, branch_id): + branch_id = int(branch_id) + if branch_id < 0 or branch_id >= self.number_of_branches: + raise ValueError( + "branch_id {} outside [0, {})".format( + branch_id, self.number_of_branches + ) + ) + + @staticmethod + def _validate_legacy_branch_id(branch_id): + if branch_id not in (None, 0): + raise ValueError("released LALSimulation family exposes only branch 0") + + +def create_family(eos, minimal=True, lalsim_module=None): + """Return a :class:`LALSimNeutronStarFamilyAdapter` for ``eos``.""" + return LALSimNeutronStarFamilyAdapter( + eos, minimal=minimal, lalsim_module=lalsim_module + ) diff --git a/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py new file mode 100644 index 000000000..0593a993a --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py @@ -0,0 +1,186 @@ +import json + +import numpy as np +import pytest + +from RIFT.physics.lalsim_eos_compat import ( + AmbiguousFamilyBranchError, + LALSimNeutronStarFamilyAdapter, +) + + +class LegacyLALSimulation: + def __init__(self): + self.create_calls = [] + self.file_calls = [] + + def SimNeutronStarEOSFromFile(self, fname): + self.file_calls.append((fname, 0)) + return "clean-eos" + + def SimNeutronStarEOSFromFileChoiceDirtyPT(self, fname, dirty): + self.file_calls.append((fname, dirty)) + return "dirty-eos" + + def CreateSimNeutronStarFamily(self, eos): + self.create_calls.append((eos,)) + return "legacy-family" + + def SimNeutronStarFamMinimumMass(self, family): + return 1.0 + + def SimNeutronStarMaximumMass(self, family): + return 3.0 + + def SimNeutronStarRadius(self, mass, family): + return 10.0 + mass + + def SimNeutronStarLoveNumberK2(self, mass, family): + return 0.1 * mass + + def SimNeutronStarCentralPressure(self, mass, family): + return 100.0 * mass + + +class MultibranchLALSimulation: + bounds = ((1.0, 2.0), (1.5, 3.0)) + + def __init__(self): + self.create_calls = [] + self.file_calls = [] + + def SimNeutronStarEOSFromFile(self, fname): + self.file_calls.append((fname, 0)) + return "clean-eos" + + def SimNeutronStarEOSFromFileChoiceDirtyPT(self, fname, dirty): + self.file_calls.append((fname, dirty)) + return "dirty-eos" + + def CreateSimNeutronStarFamily(self, eos, min_fam): + self.create_calls.append((eos, min_fam)) + return "multibranch-family" + + def SimNeutronStarFamNumberOfBranches(self, family): + return len(self.bounds) + + def SimNeutronStarFamMinMassPerBranch(self, family, branch_id): + return self.bounds[branch_id][0] + + def SimNeutronStarFamMaxMassPerBranch(self, family, branch_id): + return self.bounds[branch_id][1] + + def SimNeutronStarFamMinMass(self, family): + return min(x[0] for x in self.bounds) + + def SimNeutronStarFamMaxMass(self, family): + return max(x[1] for x in self.bounds) + + def SimNeutronStarFamRadiusOfMassPerBranch(self, mass, family, branch_id): + return 10.0 * branch_id + mass + + def SimNeutronStarFamLoveNumberK2OfMassPerBranch( + self, mass, family, branch_id + ): + return branch_id + 0.1 * mass + + def SimNeutronStarFamCentralPressureOfMassPerBranch( + self, mass, family, branch_id + ): + return 100.0 * branch_id + mass + + +def test_released_lalsimulation_uses_one_argument_family_api(): + lalsim = LegacyLALSimulation() + family = LALSimNeutronStarFamilyAdapter( + "eos", minimal=True, lalsim_module=lalsim + ) + + assert lalsim.create_calls == [("eos",)] + assert family.number_of_branches == 1 + assert family.branches_for_mass(2.0) == [0] + assert family.radius(2.0) == 12.0 + assert family.love_number_k2(2.0) == pytest.approx(0.2) + assert family.central_pressure(2.0) == 200.0 + with pytest.raises(ValueError, match="branch_id 1 outside"): + family.radius(2.0, branch_id=1) + + +def test_reviewed_lalsimulation_uses_minimal_multibranch_api(): + lalsim = MultibranchLALSimulation() + family = LALSimNeutronStarFamilyAdapter( + "eos", minimal=True, lalsim_module=lalsim + ) + + assert lalsim.create_calls == [("eos", 1)] + assert family.number_of_branches == 2 + assert family.minimum_mass() == 1.0 + assert family.maximum_mass() == 3.0 + assert family.branches_for_mass(1.25) == [0] + assert family.branches_for_mass(1.75) == [0, 1] + assert family.radius(1.75, branch_id=1) == 11.75 + assert family.love_number_k2(1.75, branch_id=1) == pytest.approx(1.175) + assert family.central_pressure(1.75, branch_id=1) == 101.75 + + +def test_twin_star_mass_requires_an_explicit_branch(): + family = LALSimNeutronStarFamilyAdapter( + "eos", lalsim_module=MultibranchLALSimulation() + ) + + with pytest.raises(AmbiguousFamilyBranchError, match=r"branches \[0, 1\]"): + family.radius(1.75) + with pytest.raises(ValueError, match="outside stable branch 0"): + family.radius(2.5, branch_id=0) + with pytest.raises(ValueError, match="outside every stable"): + family.radius(4.0) + + +def test_eosmanager_file_loader_routes_reviewed_phase_transition_api(monkeypatch): + from RIFT.physics import EOSManager + + fake_lalsim = MultibranchLALSimulation() + monkeypatch.setattr(EOSManager, "lalsim", fake_lalsim) + eos = EOSManager.EOSLALSimulationFromFile( + "new-format.dat", dirty_phase_transitions=True + ) + + assert fake_lalsim.file_calls == [("new-format.dat", 1)] + assert fake_lalsim.create_calls == [("dirty-eos", 1)] + assert eos.eos == "dirty-eos" + assert eos._get_lalsim_family_adapter().number_of_branches == 2 + + +def test_eosmanager_smoke_with_installed_released_lalsimulation(): + from RIFT.physics import EOSManager + + eos = EOSManager.EOSLALSimulation("SLy") + assert eos.branches_for_m(1.4) == [0] + assert eos.radius_from_m(1.4) > 0.0 + assert np.isfinite(eos.lambda_from_m(1.4)) + + +def test_nmb_sequence_dispatch_and_accessors_remain_compatible(tmp_path): + h5py = pytest.importorskip("h5py") + from RIFT.physics import EOSManager + + path = tmp_path / "nmb-sequence.h5" + fields = ["M", "R", "Lambda", "stable"] + sequence = np.array( + [[[1.0, 12.0, 500.0, 1.0], + [1.4, 11.5, 300.0, 1.0], + [2.0, 10.0, 50.0, 1.0]]] + ) + with h5py.File(path, "w") as stream: + stream.attrs["representation"] = "tabular_hc/1" + stream.attrs["schema_version"] = "nmbackend.nss/1" + stream.attrs["fields"] = json.dumps(fields) + stream.create_dataset("sequence", data=sequence) + + eos_sequence = EOSManager.EOSSequenceFromFile( + fname=str(path), load_ns=True, no_sort=True + ) + assert isinstance(eos_sequence, EOSManager.EOSSequenceNMB) + assert eos_sequence.m_max_of_indx(0) == pytest.approx(2.0) + assert eos_sequence.R_of_m_indx(1.4, 0) == pytest.approx(11.5) + assert eos_sequence.lambda_of_m_indx(1.4, 0) == pytest.approx(300.0) From d27d0a4684537e46861848b9b4559d6cda4f8209 Mon Sep 17 00:00:00 2001 From: "R. O'Shaughnessy" Date: Wed, 3 Jun 2026 01:14:10 +0000 Subject: [PATCH 02/16] EOSManager: add EOSSequenceNMB reader for NuclearMatter-Backend NSSequence format (drop-in, branch-aware, reuses EOSSequenceLandry accessors) --- .../Code/RIFT/physics/EOSManager.py | 98 +++++++++++++++++++ 1 file changed, 98 insertions(+) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py b/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py index d131b5be2..13503e099 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py @@ -1605,6 +1605,104 @@ def extract_one_eos_object(self, indx=None,name_eos=None,fail_if=None,**kwargs): my_eos = EOSFromTabularData(name=name_to_use, eos_data=dat_copy,**kwargs) # tabular data inputs need to be cgs and in correct units return my_eos +class EOSSequenceNMB(EOSSequenceLandry): + """Drop-in reader for the NuclearMatter-Backend ``NSSequence`` HDF5 format. + + The NSSequence file stores every quantity as a function of central + pseudo-enthalpy h_c (monotone along the sequence), with an explicit ``stable`` + flag, in a single ``(n_eos, n_pts, n_fields)`` dataset (see + docs/rift-sequence-audit.md in NuclearMatter-Backend). This subclass reads that + file, extracts the **stable rising branch** (M increasing up to M_max) for each + EOS into the same in-memory ``eos_ns_tov`` dict of {M,R,Lambda} structured arrays + that EOSSequenceLandry uses -- so all inherited accessors + (``lambda_of_m_indx``, ``R_of_m_indx``, ``m_max_of_indx``, ``lookup_closest``, + ``oned_order_values``) work unchanged and are branch-safe by construction. + + Only ``load_ns`` is honoured (TOV sequence); the optional microphysical EOS + tables are not read here (use the legacy emitter / EOSSequenceLandry for those). + """ + + @staticmethod + def _stable_rising(M, R, Lam, stable): + ok = np.isfinite(M) & (M > 0) + M, R, Lam, st = M[ok], R[ok], Lam[ok], stable[ok] > 0.5 + if M.size < 2: + return M, R, Lam + imax = int(np.argmax(np.where(st, M, -np.inf))) + M, R, Lam = M[:imax + 1], R[:imax + 1], Lam[:imax + 1] + o = np.argsort(M) + return M[o], R[o], Lam[o] + + def __init__(self, name=None, fname=None, load_eos=False, load_ns=True, + oned_order_name=None, oned_order_mass=None, no_sort=True, + verbose=False, eos_tables_units=None): + import json + import h5py + self.name = name + self.fname = fname + self.eos_ids = None + self.eos_names = None + self.eos_tables = None + self.eos_tables_units = None + self.eos_ns_tov = None + self.oned_order_name = None + self.oned_order_mass = oned_order_mass + self.oned_order_values = None + self.oned_order_indx_original = None + self.oned_order_indx_sorted = None + self.oned_order_sorted = False + self.verbose = verbose + + with h5py.File(self.fname, 'r') as f: + rep = str(f.attrs.get("representation", "tabular_hc/1")) + if not rep.startswith("tabular"): + raise NotImplementedError( + "EOSSequenceNMB: representation {!r} not supported " + "(reserved for future compressed/functional representations)".format(rep)) + fields = json.loads(f.attrs["fields"]) + col = {k: j for j, k in enumerate(fields)} + seq = f["sequence"][:] # (n_eos, n_pts, n_fields) + + n_eos = seq.shape[0] + self.eos_names = np.array(["eos_{}".format(k) for k in range(n_eos)], dtype=str) + self.eos_ids = list(range(n_eos)) + self.eos_ns_tov = {} + for k in range(n_eos): + s = seq[k] + M, R, Lam = self._stable_rising(s[:, col["M"]], s[:, col["R"]], + s[:, col["Lambda"]], s[:, col["stable"]]) + rec = np.zeros(M.size, dtype=[("M", "f8"), ("R", "f8"), ("Lambda", "f8")]) + rec["M"], rec["R"], rec["Lambda"] = M, R, Lam + self.eos_ns_tov["eos_{}".format(k)] = rec + + # Build the 1-D ordering statistic exactly as EOSSequenceLandry does. + create_order = False + if oned_order_name in ('R', 'r'): + create_order, self.oned_order_name = True, 'R' + if oned_order_name in ('Lambda', 'lambda'): + create_order, self.oned_order_name = True, 'Lambda' + if not self.oned_order_mass: + create_order = False + if create_order: + self.oned_order_indx_original = np.arange(len(self.eos_names)) + vals = np.zeros(len(self.eos_names)) + for indx in range(len(self.eos_names)): + if self.oned_order_name == 'Lambda': + vals[indx] = self.lambda_of_m_indx(self.oned_order_mass, indx) + else: + vals[indx] = self.R_of_m_indx(self.oned_order_mass, indx) + self.oned_order_indx_sorted = np.argsort(vals) + if no_sort: + self.oned_order_values = vals + else: + self.eos_names = self.eos_names[self.oned_order_indx_sorted] + self.oned_order_values = vals[self.oned_order_indx_sorted] + self.oned_order_indx_original = self.oned_order_indx_original[self.oned_order_indx_sorted] + self.oned_order_indx_sorted = np.arange(len(self.eos_names)) + self.oned_order_sorted = True + return None + + #### #### General lalsimulation interfacing #### From 23d3ddea822c6c70e4f13879c7f9c95136fed60f Mon Sep 17 00:00:00 2001 From: "R. O'Shaughnessy" Date: Wed, 3 Jun 2026 01:28:52 +0000 Subject: [PATCH 03/16] CIP: auto-detect NSSequence vs EOSSequenceLandry via EOSSequenceFromFile factory (transparent tabular-EOS drop-in for util_RIFT_pseudo_pipe.py --internal-tabular-eos-file) --- .../Code/RIFT/physics/EOSManager.py | 24 +++++++++++++++++++ ...ctIntrinsicPosterior_GenericCoordinates.py | 3 ++- 2 files changed, 26 insertions(+), 1 deletion(-) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py b/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py index 13503e099..f89c654a5 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py @@ -1703,6 +1703,30 @@ def __init__(self, name=None, fname=None, load_eos=False, load_ns=True, return None +def EOSSequenceFromFile(fname=None, **kwargs): + """Open an EOS sequence file, auto-detecting the format. + + Returns an ``EOSSequenceNMB`` for NuclearMatter-Backend ``NSSequence`` files + (identified by the ``representation`` / ``schema_version`` HDF5 attribute) and an + ``EOSSequenceLandry`` otherwise. Both expose the identical consumer API + (``oned_order_values``, ``lambda_of_m_indx``, ``R_of_m_indx``, ``m_max_of_indx``, + ``lookup_closest``), so callers can pass either format transparently. + """ + import h5py + is_nmb = False + try: + with h5py.File(fname, 'r') as f: + a = f.attrs + rep = str(a.get("representation", "")) + schema = str(a.get("schema_version", "")) + is_nmb = rep.startswith("tabular") or schema.startswith("nmbackend") + except Exception: + is_nmb = False + if is_nmb: + return EOSSequenceNMB(fname=fname, **kwargs) + return EOSSequenceLandry(fname=fname, **kwargs) + + #### #### General lalsimulation interfacing #### diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py index 941294d49..b199d8653 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py @@ -1924,7 +1924,8 @@ def fit_gp_sparse(x): if mc_ref > 1e10: mc_ref = mc_ref/lal.MSUN_SI m_ref = mc_ref*np.power(2, 1./5.) # assume equal mass - my_eos_sequence = EOSManager.EOSSequenceLandry(fname=opts.tabular_eos_file, load_ns=True, oned_order_name='Lambda', oned_order_mass=m_ref, no_sort = True) + # auto-detect EOSSequenceLandry vs NuclearMatter-Backend NSSequence format + my_eos_sequence = EOSManager.EOSSequenceFromFile(fname=opts.tabular_eos_file, load_ns=True, oned_order_name='Lambda', oned_order_mass=m_ref, no_sort = True) # Define prior, NOT NORMALIZED prior_map['ordering'] =lambda x: np.ones(x.shape) From a86d448a4e3abc10a9d4afe724d7080e35de9ae0 Mon Sep 17 00:00:00 2001 From: "R. O'Shaughnessy" Date: Wed, 3 Jun 2026 02:06:08 +0000 Subject: [PATCH 04/16] EOSManager: EOSSequencePCA reader for compressed pca_hc/1 files; EOSSequenceFromFile now dispatches pca/tabular/legacy --- .../Code/RIFT/physics/EOSManager.py | 100 +++++++++++++++--- 1 file changed, 88 insertions(+), 12 deletions(-) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py b/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py index f89c654a5..8e22e9060 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py @@ -1703,26 +1703,102 @@ def __init__(self, name=None, fname=None, load_eos=False, load_ns=True, return None +class EOSSequencePCA(EOSSequenceNMB): + """Reader for the compressed NuclearMatter-Backend ``pca_hc/1`` representation. + + The file stores a per-channel PCA decomposition of the M(u), R(u), logLambda(u) + curves (mean + basis ``components`` + per-EOS ``coeffs``). We reconstruct each + EOS's curves, take the stable rising branch, and populate the same in-memory + ``eos_ns_tov`` dict EOSSequenceLandry/EOSSequenceNMB use -- so every inherited + accessor works unchanged. Self-contained (numpy only); no nmbackend dependency. + """ + + def __init__(self, name=None, fname=None, load_eos=False, load_ns=True, + oned_order_name=None, oned_order_mass=None, no_sort=True, + verbose=False, eos_tables_units=None): + import json + import h5py + self.name = name; self.fname = fname + self.eos_ids = None; self.eos_names = None + self.eos_tables = None; self.eos_tables_units = None; self.eos_ns_tov = None + self.oned_order_name = None; self.oned_order_mass = oned_order_mass + self.oned_order_values = None + self.oned_order_indx_original = None; self.oned_order_indx_sorted = None + self.oned_order_sorted = False; self.verbose = verbose + + with h5py.File(self.fname, 'r') as f: + channels = json.loads(f.attrs["channels"]) + mean = f["mean"][:] # (3, n_pts) + comps = f["components"][:] # (3, n_comp, n_pts) + coeffs = f["coeffs"][:] # (n_eos, 3, n_comp) + iM, iR, iL = (channels.index("M"), channels.index("R"), + channels.index("logLambda")) + n_eos = coeffs.shape[0] + self.eos_names = np.array(["eos_{}".format(k) for k in range(n_eos)], dtype=str) + self.eos_ids = list(range(n_eos)) + self.eos_ns_tov = {} + for k in range(n_eos): + rec_curves = mean + np.einsum("ck,ckp->cp", coeffs[k], comps) + M, R, Lam = rec_curves[iM], rec_curves[iR], np.exp(rec_curves[iL]) + stable = np.concatenate([[True], np.diff(M) > 0]) + Mb, Rb, Lb = self._stable_rising(M, R, Lam, stable.astype(float)) + rec = np.zeros(Mb.size, dtype=[("M", "f8"), ("R", "f8"), ("Lambda", "f8")]) + rec["M"], rec["R"], rec["Lambda"] = Mb, Rb, Lb + self.eos_ns_tov["eos_{}".format(k)] = rec + self._build_ordering(oned_order_name, no_sort) + return None + + def _build_ordering(self, oned_order_name, no_sort): + create_order = False + if oned_order_name in ('R', 'r'): + create_order, self.oned_order_name = True, 'R' + if oned_order_name in ('Lambda', 'lambda'): + create_order, self.oned_order_name = True, 'Lambda' + if not self.oned_order_mass: + create_order = False + if not create_order: + return + self.oned_order_indx_original = np.arange(len(self.eos_names)) + vals = np.zeros(len(self.eos_names)) + for indx in range(len(self.eos_names)): + vals[indx] = (self.lambda_of_m_indx(self.oned_order_mass, indx) + if self.oned_order_name == 'Lambda' + else self.R_of_m_indx(self.oned_order_mass, indx)) + self.oned_order_indx_sorted = np.argsort(vals) + if no_sort: + self.oned_order_values = vals + else: + self.eos_names = self.eos_names[self.oned_order_indx_sorted] + self.oned_order_values = vals[self.oned_order_indx_sorted] + self.oned_order_indx_original = self.oned_order_indx_original[self.oned_order_indx_sorted] + self.oned_order_indx_sorted = np.arange(len(self.eos_names)) + self.oned_order_sorted = True + + def EOSSequenceFromFile(fname=None, **kwargs): """Open an EOS sequence file, auto-detecting the format. - Returns an ``EOSSequenceNMB`` for NuclearMatter-Backend ``NSSequence`` files - (identified by the ``representation`` / ``schema_version`` HDF5 attribute) and an - ``EOSSequenceLandry`` otherwise. Both expose the identical consumer API - (``oned_order_values``, ``lambda_of_m_indx``, ``R_of_m_indx``, ``m_max_of_indx``, - ``lookup_closest``), so callers can pass either format transparently. + Dispatches on the HDF5 ``representation`` / ``schema_version`` attribute: + + * ``pca_hc`` (``nmbackend.pca``) -> :class:`EOSSequencePCA` (compressed); + * ``tabular`` (``nmbackend.nss``) -> :class:`EOSSequenceNMB` (tabular); + * anything else -> :class:`EOSSequenceLandry` (legacy/LCEHL). + + All expose the identical consumer API (``oned_order_values``, + ``lambda_of_m_indx``, ``R_of_m_indx``, ``m_max_of_indx``, ``lookup_closest``), so + callers can pass any of the three file types transparently. """ import h5py - is_nmb = False + rep = schema = "" try: with h5py.File(fname, 'r') as f: - a = f.attrs - rep = str(a.get("representation", "")) - schema = str(a.get("schema_version", "")) - is_nmb = rep.startswith("tabular") or schema.startswith("nmbackend") + rep = str(f.attrs.get("representation", "")) + schema = str(f.attrs.get("schema_version", "")) except Exception: - is_nmb = False - if is_nmb: + rep = schema = "" + if rep.startswith("pca") or schema.startswith("nmbackend.pca"): + return EOSSequencePCA(fname=fname, **kwargs) + if rep.startswith("tabular") or schema.startswith("nmbackend"): return EOSSequenceNMB(fname=fname, **kwargs) return EOSSequenceLandry(fname=fname, **kwargs) From 05edff1f49fb2c3c2fd07880949d0ef5406a76a6 Mon Sep 17 00:00:00 2001 From: "R. O'Shaughnessy" Date: Wed, 3 Jun 2026 14:06:12 +0000 Subject: [PATCH 05/16] CIP: fixed single-EOS-realization mode for sequence files (--using-eos nmbseq:: via EOSSequenceSingleIndex) -- enables exact per-EOS evidence for tabular/pca draws --- .../Code/RIFT/physics/EOSManager.py | 37 +++++++++++++++++++ ...ctIntrinsicPosterior_GenericCoordinates.py | 7 ++++ 2 files changed, 44 insertions(+) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py b/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py index 8e22e9060..3342147d9 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py @@ -1775,6 +1775,43 @@ def _build_ordering(self, oned_order_name, no_sort): self.oned_order_sorted = True +class EOSSequenceSingleIndex: + """A SINGLE EOS realization drawn from a sequence file, exposed with the + fixed-EOS (``--using-eos``) interface: ``lambda_from_m(m_Msun)``. + + This enables the exact per-EOS-evidence pattern (one full CIP evidence per + realization, MARG-style) for tabular/compressed sequence files -- the + reference computation against which the ordering-statistic (tabular + hyperpipeline) approximation is validated. + """ + + def __init__(self, fname=None, index=0, name=None): + self.name = name or "nmbseq_{}_{}".format(fname, index) + self.fname = fname + self.index = int(index) + self._seq = EOSSequenceFromFile(fname=fname, load_ns=True, no_sort=True) + if not (0 <= self.index < len(self._seq.eos_names)): + raise ValueError("EOS index {} out of range (n={})".format( + self.index, len(self._seq.eos_names))) + self.mMaxMsun = float(self._seq.m_max_of_indx(self.index)) + + def lambda_from_m(self, m): + # unit auto-detection as in EOSConcrete.lambda_from_m + m_Msun = m / lal.MSUN_SI if m > 1e15 else m + if m_Msun > 0.999 * self.mMaxMsun: + return 1e-8 + val = self._seq.lambda_of_m_indx(m_Msun, self.index) + return float(val) if np.isfinite(val) else 1e-8 + + def lambda_from_m_vector(self, m): + if not isinstance(m, np.ndarray): + return self.lambda_from_m(m) + return np.array([self.lambda_from_m(x) for x in m]) + + def R_from_m(self, m_Msun): + return self._seq.R_of_m_indx(m_Msun, self.index) + + def EOSSequenceFromFile(fname=None, **kwargs): """Open an EOS sequence file, auto-detecting the format. diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py index b199d8653..2f8dfc795 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py @@ -500,6 +500,13 @@ def extract_combination_from_LI(samples_LI, p): spec_params['gamma4']=spec_param_array[3] eos_base = EOSManager.EOSLindblomSpectralSoundSpeedVersusPressure(name=eos_name,spec_params=spec_params,use_lal_spec_eos=not opts.no_use_lal_eos) my_eos = eos_base + elif eos_name.startswith('nmbseq:'): + # fixed single EOS realization from a sequence file (tabular or pca): + # --using-eos nmbseq:: + # The exact per-EOS-evidence ("painful") mode for sequence draws. + _, seq_fname, seq_indx = eos_name.split(':') + my_eos = EOSManager.EOSSequenceSingleIndex(fname=seq_fname, + index=int(seq_indx)) elif 'lal_' in eos_name: eos_name = eos_name.replace('lal_','') my_eos = EOSManager.EOSLALSimulation(name=eos_name) From 606c102f7bb3a2fcf09089c2b3318e49772c6d42 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Fri, 28 Aug 2026 17:29:54 -0400 Subject: [PATCH 06/16] EOSManager: preserve downstream branch contracts --- CHANGES.rst | 6 +- .../Code/RIFT/physics/EOSManager.py | 113 ++++++++++++++++-- ...ctIntrinsicPosterior_GenericCoordinates.py | 24 +++- .../Code/test/test_lalsim_eos_compat.py | 78 ++++++++++++ docs/eos-interface-contract.md | 69 +++++++++++ 5 files changed, 278 insertions(+), 12 deletions(-) create mode 100644 docs/eos-interface-contract.md diff --git a/CHANGES.rst b/CHANGES.rst index afd21236e..e0df8f6c7 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -2,7 +2,11 @@ -------------------------------- - EOSManager supports released and reviewed multipart/multibranch LALSimulation family interfaces while retaining the existing scalar and NuclearMatter-Backend - sequence contracts. See ``docs/eos-interface-contract.md``. + sequence contracts. Fixed-EOS CIP can select a LAL family with + ``--using-eos-branch``. Native NMB/PCA files retain the ``nmbseq:`` primary-branch + contract, with disconnected stable runs split before interpolation. The CIP + piecewise-polytrope constructor keyword is corrected for Kedia-style workflows. + See ``docs/eos-interface-contract.md``. 0.0.17.13 --------- diff --git a/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py b/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py index 3342147d9..c111ae0e3 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py @@ -101,6 +101,19 @@ def branches_for_m(self, m): m = m * lal.MSUN_SI return self._get_lalsim_family_adapter().branches_for_mass(m) + def for_branch(self, branch_id): + """Return a legacy-scalar view restricted to one stable LAL branch. + + Existing RIFT consumers call ``lambda_from_m(m)`` without a branch + keyword. The view preserves that API while making the branch choice + explicit at construction time. + """ + if getattr(self, "_lalsim_family_adapter", None) is None: + raise TypeError( + "branch selection requires a LALSimulation-backed EOS family" + ) + return EOSBranchView(self, branch_id) + def radius_from_m(self, m, branch_id=None): """Return radius in metres; require ``branch_id`` for twin stars.""" if m < 10**15: @@ -260,6 +273,49 @@ def test_speed_of_sound_causal(self, test_only_under_mmax=True,fast_test=True): print(h[indx], vs_internal[indx]) return not np.any(vs_internal>1.1) # allow buffer, so we have some threshold + +class EOSBranchView: + """A selected stable-family branch with the historical scalar EOS surface. + + This wrapper is deliberately small: attributes not related to the stellar + family are delegated to the source EOS, while all mass-radius-tidal queries + are pinned to ``branch_id``. It is the compatibility bridge for callers + such as CIP and hyperpipe that cannot pass a branch keyword on every lookup. + """ + + def __init__(self, source, branch_id): + self.source = source + self.branch_id = int(branch_id) + adapter = source._get_lalsim_family_adapter() + adapter._validate_branch_id(self.branch_id) + self.name = "{}[branch={}]".format(source.name, self.branch_id) + self.eos = source.eos + self.eos_fam = source.eos_fam + self.mMaxMsun = ( + adapter.maximum_mass(self.branch_id) / lal.MSUN_SI + ) + + def branches_for_m(self, m): + available = self.source.branches_for_m(m) + return [self.branch_id] if self.branch_id in available else [] + + def radius_from_m(self, m): + return self.source.radius_from_m(m, branch_id=self.branch_id) + + def lambda_from_m(self, m): + m_msun = m / lal.MSUN_SI if m > 10**15 else m + if m_msun > 0.999 * self.mMaxMsun: + return 1e-8 + return self.source.lambda_from_m(m, branch_id=self.branch_id) + + def lambda_from_m_vector(self, m): + if not isinstance(m, np.ndarray): + return self.lambda_from_m(m) + return np.array([self.lambda_from_m(m_here) for m_here in m]) + + def __getattr__(self, name): + return getattr(self.source, name) + ### ### SERVICE 1: lalsimutils structure ### @@ -1623,15 +1679,44 @@ class EOSSequenceNMB(EOSSequenceLandry): """ @staticmethod - def _stable_rising(M, R, Lam, stable): - ok = np.isfinite(M) & (M > 0) - M, R, Lam, st = M[ok], R[ok], Lam[ok], stable[ok] > 0.5 - if M.size < 2: - return M, R, Lam - imax = int(np.argmax(np.where(st, M, -np.inf))) - M, R, Lam = M[:imax + 1], R[:imax + 1], Lam[:imax + 1] - o = np.argsort(M) - return M[o], R[o], Lam[o] + def _stable_branches(M, R, Lam, stable): + """Split an h_c-ordered sequence into stable, mass-rising branches. + + ``NSSequence`` stores central-enthalpy order. Filtering only at the + global maximum silently mixed disconnected stable branches with the + unstable interval between them. Runs are therefore split whenever the + stored stability flag is false or mass ceases to rise. + """ + M = np.asarray(M); R = np.asarray(R); Lam = np.asarray(Lam) + valid = (np.isfinite(M) & np.isfinite(R) & np.isfinite(Lam) + & (M > 0) & (np.asarray(stable) > 0.5)) + branches = [] + start = None + for k in range(len(M)): + continues = (valid[k] and start is not None + and k > 0 and valid[k - 1] and M[k] > M[k - 1]) + if valid[k] and start is None: + start = k + elif valid[k] and not continues: + if k - start >= 2: + branches.append((M[start:k], R[start:k], Lam[start:k])) + start = k + elif not valid[k] and start is not None: + if k - start >= 2: + branches.append((M[start:k], R[start:k], Lam[start:k])) + start = None + if start is not None and len(M) - start >= 2: + branches.append((M[start:], R[start:], Lam[start:])) + return branches + + @classmethod + def _stable_rising(cls, M, R, Lam, stable): + """Return the primary stable branch used by the NMB v1 RIFT contract.""" + branches = cls._stable_branches(M, R, Lam, stable) + if branches: + return branches[0] + empty = np.array([], dtype=float) + return empty, empty.copy(), empty.copy() def __init__(self, name=None, fname=None, load_eos=False, load_ns=True, oned_order_name=None, oned_order_mass=None, no_sort=True, @@ -1667,8 +1752,14 @@ def __init__(self, name=None, fname=None, load_eos=False, load_ns=True, self.eos_names = np.array(["eos_{}".format(k) for k in range(n_eos)], dtype=str) self.eos_ids = list(range(n_eos)) self.eos_ns_tov = {} + self.stable_branch_counts = np.zeros(n_eos, dtype=int) for k in range(n_eos): s = seq[k] + branches = self._stable_branches( + s[:, col["M"]], s[:, col["R"]], + s[:, col["Lambda"]], s[:, col["stable"]] + ) + self.stable_branch_counts[k] = len(branches) M, R, Lam = self._stable_rising(s[:, col["M"]], s[:, col["R"]], s[:, col["Lambda"]], s[:, col["stable"]]) rec = np.zeros(M.size, dtype=[("M", "f8"), ("R", "f8"), ("Lambda", "f8")]) @@ -1737,10 +1828,14 @@ def __init__(self, name=None, fname=None, load_eos=False, load_ns=True, self.eos_names = np.array(["eos_{}".format(k) for k in range(n_eos)], dtype=str) self.eos_ids = list(range(n_eos)) self.eos_ns_tov = {} + self.stable_branch_counts = np.zeros(n_eos, dtype=int) for k in range(n_eos): rec_curves = mean + np.einsum("ck,ckp->cp", coeffs[k], comps) M, R, Lam = rec_curves[iM], rec_curves[iR], np.exp(rec_curves[iL]) stable = np.concatenate([[True], np.diff(M) > 0]) + self.stable_branch_counts[k] = len( + self._stable_branches(M, R, Lam, stable.astype(float)) + ) Mb, Rb, Lb = self._stable_rising(M, R, Lam, stable.astype(float)) rec = np.zeros(Mb.size, dtype=[("M", "f8"), ("R", "f8"), ("Lambda", "f8")]) rec["M"], rec["R"], rec["Lambda"] = Mb, Rb, Lb diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py index 2f8dfc795..c9a15d21c 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py @@ -335,6 +335,7 @@ def extract_combination_from_LI(samples_LI, p): parser.add_argument("--tabular-eos-file-format",type=str,default=None,help="Format of tabular file of EOS to use. The default prior will be UNIFORM in this table!") parser.add_argument("--tabular-eos-order-statistic",type=str,default=None,help="Order statistic to use. Options will include R1p4, LambdaTildeQ1, and ...}") parser.add_argument("--using-eos", type=str, default=None, help="Name of EOS. Fit parameter list should physically use lambda1, lambda2 information (but need not). If starts with 'file:', uses a filename with EOS parameters ") +parser.add_argument("--using-eos-branch", type=int, default=None, help="Select one stable LALSimulation family branch while preserving the fixed-EOS lambda_from_m(m) interface. Required for an explicitly chosen twin-star branch; not applicable to the primary-branch nmbseq v1 contract.") parser.add_argument("--using-eos-index", type=int, default=None, help="Index of EOS parameters in file.") parser.add_argument("--no-use-lal-eos",action='store_true',help="Do not use LAL EOS interface. Used for spectral EOS. Do not use this.") parser.add_argument("--no-matter1", action='store_true', help="Set the lambda parameters to zero (BBH) but return them") @@ -416,6 +417,16 @@ def extract_combination_from_LI(samples_LI, p): print(" warning: input EOS index, but not using it; presumably you are doing a model-free test ") if not(opts.input_eos_index) and (opts.tabular_eos_file): raise Exception(" Fail: must process EOS input to be able to use it ") +# ensure using_eos_index valid for eos file length +if opts.using_eos and opts.using_eos.startswith('file:') and not(opts.using_eos_index is None): + fname = opts.using_eos.replace('file:', '') + try: + dat = np.loadtxt(fname)[opts.using_eos_index] + except Exception as e: + print(" Fail: EOS index out of range:\n ",e) + sys.exit(0) +if opts.using_eos_branch is not None and opts.using_eos is None: + raise ValueError("--using-eos-branch also requires --using-eos") my_eos=None #option to be used if gridded values not calculated assuming EOS @@ -457,7 +468,7 @@ def extract_combination_from_LI(samples_LI, p): spec_params['gamma1'] = spec_param_array[1] spec_params['gamma2'] = spec_param_array[2] spec_params['gamma3'] = spec_param_array[3] - eos_base = EOSManager.EOSPiecewisePolytrope(name=eos_name,params_dict=spec_params) + eos_base = EOSManager.EOSPiecewisePolytrope(name=eos_name,param_dict=spec_params) my_eos = eos_base else: raise Exception("Unknown method for parametric EOS data file {} : {} ".format(eos_name,opts.eos_param)) @@ -522,6 +533,16 @@ def extract_combination_from_LI(samples_LI, p): else: my_eos = EOSManager.EOSFromDataFile(name=eos_name,fname =EOSManager.dirEOSTablesBase+"/" + eos_name+".dat") + if opts.using_eos_branch is not None: + if not hasattr(my_eos, "for_branch"): + raise ValueError( + "--using-eos-branch requires a LALSimulation-backed EOS. " + "The nmbseq v1 interface intentionally exposes its primary " + "stable branch; multi-branch NMB inference requires the " + "central-enthalpy sequence path." + ) + my_eos = my_eos.for_branch(opts.using_eos_branch) + with open('args.txt','w') as fp: import sys @@ -3626,4 +3647,3 @@ def parse_corr_params(my_str): sys.exit(0) - diff --git a/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py index 0593a993a..e6f89418d 100644 --- a/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py +++ b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py @@ -1,5 +1,6 @@ import json +import lal import numpy as np import pytest @@ -90,6 +91,13 @@ def SimNeutronStarFamCentralPressureOfMassPerBranch( return 100.0 * branch_id + mass +class StellarMassMultibranchLALSimulation(MultibranchLALSimulation): + bounds = tuple( + (lower * lal.MSUN_SI, upper * lal.MSUN_SI) + for lower, upper in MultibranchLALSimulation.bounds + ) + + def test_released_lalsimulation_uses_one_argument_family_api(): lalsim = LegacyLALSimulation() family = LALSimNeutronStarFamilyAdapter( @@ -160,6 +168,44 @@ def test_eosmanager_smoke_with_installed_released_lalsimulation(): assert np.isfinite(eos.lambda_from_m(1.4)) +def test_kedia_parametric_eos_scalar_interfaces_remain_compatible(): + from RIFT.physics import EOSManager + + spectral = EOSManager.EOSLindblomSpectral( + name="contract-spectral", + spec_params=dict(gamma1=1.0, gamma2=1.0, gamma3=0.0, gamma4=0.0), + use_lal_spec_eos=True, + ) + piecewise = EOSManager.EOSPiecewisePolytrope( + name="contract-piecewise", + param_dict=dict( + logP1=34.269, gamma1=2.830, gamma2=3.445, gamma3=3.348 + ), + ) + + assert np.isfinite(spectral.lambda_from_m(1.4)) + assert np.isfinite(piecewise.lambda_from_m(1.4)) + + +def test_selected_branch_view_preserves_legacy_scalar_consumer_api(monkeypatch): + from RIFT.physics import EOSManager + + fake_lalsim = StellarMassMultibranchLALSimulation() + monkeypatch.setattr(EOSManager, "lalsim", fake_lalsim) + eos = EOSManager.EOSLALSimulationFromFile("twin-star.dat") + + primary = eos.for_branch(0) + secondary = eos.for_branch(1) + assert primary.mMaxMsun == pytest.approx(2.0) + assert secondary.mMaxMsun == pytest.approx(3.0) + assert secondary.branches_for_m(1.75) == [1] + assert secondary.radius_from_m(1.75) == pytest.approx( + 10.0 + 1.75 * lal.MSUN_SI + ) + assert np.isfinite(secondary.lambda_from_m(1.75)) + assert primary.lambda_from_m(2.5) == pytest.approx(1e-8) + + def test_nmb_sequence_dispatch_and_accessors_remain_compatible(tmp_path): h5py = pytest.importorskip("h5py") from RIFT.physics import EOSManager @@ -184,3 +230,35 @@ def test_nmb_sequence_dispatch_and_accessors_remain_compatible(tmp_path): assert eos_sequence.m_max_of_indx(0) == pytest.approx(2.0) assert eos_sequence.R_of_m_indx(1.4, 0) == pytest.approx(11.5) assert eos_sequence.lambda_of_m_indx(1.4, 0) == pytest.approx(300.0) + + +def test_nmb_primary_branch_contract_does_not_mix_disconnected_branches(tmp_path): + h5py = pytest.importorskip("h5py") + from RIFT.physics import EOSManager + + path = tmp_path / "nmb-twin-sequence.h5" + fields = ["hc", "M", "R", "Lambda", "stable"] + sequence = np.array( + [[[0.1, 1.0, 13.0, 600.0, 1.0], + [0.2, 2.0, 11.0, 100.0, 1.0], + [0.3, 1.8, 10.8, 80.0, 0.0], + [0.4, 1.6, 10.0, 60.0, 1.0], + [0.5, 2.1, 9.0, 20.0, 1.0]]] + ) + with h5py.File(path, "w") as stream: + stream.attrs["representation"] = "tabular_hc/1" + stream.attrs["schema_version"] = "nmbackend.nss/1" + stream.attrs["fields"] = json.dumps(fields) + stream.create_dataset("sequence", data=sequence) + + eos_sequence = EOSManager.EOSSequenceFromFile( + fname=str(path), load_ns=True, no_sort=True + ) + assert eos_sequence.stable_branch_counts[0] == 2 + assert eos_sequence.m_max_of_indx(0) == pytest.approx(2.0) + expected_primary_radius = np.exp( + np.interp(1.8, [1.0, 2.0], np.log([13.0, 11.0])) + ) + assert eos_sequence.R_of_m_indx(1.8, 0) == pytest.approx( + expected_primary_radius + ) diff --git a/docs/eos-interface-contract.md b/docs/eos-interface-contract.md new file mode 100644 index 000000000..a8815c40c --- /dev/null +++ b/docs/eos-interface-contract.md @@ -0,0 +1,69 @@ +# EOS interface contract across RIFT, LALSimulation, and NuclearMatter-Backend + +RIFT keeps the historical fixed-EOS consumer surface: + +```python +eos.lambda_from_m(m) +eos.lambda_from_m_vector(masses) +eos.mMaxMsun +``` + +Mass arguments may be in solar masses or SI kg, as before. + +## LALSimulation families + +Released LALSimulation has one stable family and needs no user change. The +reviewed multipart TOV API can return several stable branches. RIFT exposes +`branches_for_m(m)` and accepts `branch_id` on direct `radius_from_m` and +`lambda_from_m` calls. An ambiguous twin-star mass raises instead of silently +choosing a solution. + +Legacy scalar consumers select a branch once with `eos.for_branch(branch_id)`. +The fixed-EOS CIP driver exposes the same operation as: + +```text +--using-eos lal_ --using-eos-branch +``` + +For pseudo-pipe workflows, forward the flag with +`--manual-extra-cip-args`. Hydra hyperpipe configurations can put it in the +post driver's `extra-args` when that driver is the fixed-EOS CIP executable. + +## NuclearMatter-Backend sequences + +The `nmbackend.nss/1` (`tabular_hc/1`) and `nmbackend.pca/1` (`pca_hc/1`) +producer tags, field names, and `EOSSequenceFromFile` dispatch remain unchanged. +The nmb-papers exact-evidence command remains: + +```text +--using-eos nmbseq:: +``` + +That v1 consumer contract intentionally projects each EOS onto its **primary +stable mass-rising branch**. RIFT now splits central-enthalpy-ordered data at +unstable or decreasing-mass intervals before doing M-to-Lambda interpolation; +it never mixes disconnected branches. `stable_branch_counts[index]` records +when a native sequence contains more than one stable run. + +`--using-eos-branch` is deliberately rejected for `nmbseq:`. Full NMB +multi-branch inference is not a scalar `M -> Lambda` problem: the required path +is a central-enthalpy likelihood/integration coordinate (one per star), using +the branch-explicit native sequence rather than the legacy Landry projection. +Until that inference path lands, generate or consume the primary branch for +paper-production parity and treat disconnected branches as a separate model. + +## Downstream compatibility + +- Existing Kedia-style spectral and piecewise-polytrope EOS inference keeps the + same constructor and scalar lookup APIs. Ordinary single-branch models need + no new option. +- Existing nmb-papers and hyperpipe fixed-sequence runs keep the same HDF5 tags + and `nmbseq:` syntax. +- Twin-star analyses using reviewed LALSimulation must select a branch for + legacy scalar workflows, or move to the central-enthalpy inference path when + marginalization over branch identity is scientifically required. + +The drift-sentinel registry should eventually declare an EOS contract group +with LALSuite and NuclearMatter-Backend as producers and RIFT/nmb-papers as +consumers. The current registry only covers RIFT/hyperpipe operational archive +and queue boundaries, so it cannot detect EOS schema or callable drift yet. From 13bc0b0c62ba53d38ed171119d0a33b54ebc2c7f Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Fri, 28 Aug 2026 17:34:37 -0400 Subject: [PATCH 07/16] docs: qualify hyperpipe branch forwarding --- docs/eos-interface-contract.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/docs/eos-interface-contract.md b/docs/eos-interface-contract.md index a8815c40c..21bd57621 100644 --- a/docs/eos-interface-contract.md +++ b/docs/eos-interface-contract.md @@ -26,8 +26,9 @@ The fixed-EOS CIP driver exposes the same operation as: ``` For pseudo-pipe workflows, forward the flag with -`--manual-extra-cip-args`. Hydra hyperpipe configurations can put it in the -post driver's `extra-args` when that driver is the fixed-EOS CIP executable. +`--manual-extra-cip-args`. On O4d, Hydra hyperpipe configurations can put it +in the post driver's `extra-args` when that driver is the fixed-EOS CIP +executable. ## NuclearMatter-Backend sequences From 44b64ee29eadbdd5f33fdd5ff2bab9bf4ca4ef44 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Fri, 28 Aug 2026 17:35:21 -0400 Subject: [PATCH 08/16] CIP: validate four-parameter EOS rows --- CHANGES.rst | 3 ++- .../util_ConstructIntrinsicPosterior_GenericCoordinates.py | 5 ++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index e0df8f6c7..fa9dd865e 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -5,7 +5,8 @@ sequence contracts. Fixed-EOS CIP can select a LAL family with ``--using-eos-branch``. Native NMB/PCA files retain the ``nmbseq:`` primary-branch contract, with disconnected stable runs split before interpolation. The CIP - piecewise-polytrope constructor keyword is corrected for Kedia-style workflows. + piecewise-polytrope constructor keyword and four-parameter row guards are corrected + for Kedia-style workflows. See ``docs/eos-interface-contract.md``. 0.0.17.13 diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py index c9a15d21c..7cb98d2cd 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py @@ -453,7 +453,7 @@ def extract_combination_from_LI(samples_LI, p): spec_params['gamma4']=spec_param_array[3] eos_base = EOSManager.EOSLindblomSpectral(name=eos_name,spec_params=spec_params,use_lal_spec_eos=not opts.no_use_lal_eos) my_eos=eos_base - elif opts.eos_param == 'cs_spectral' and len(spec_param_array >=4): + elif opts.eos_param == 'cs_spectral' and len(spec_param_array) >= 4: spec_params ={} spec_params['gamma1']=spec_param_array[0] spec_params['gamma2']=spec_param_array[1] @@ -462,7 +462,7 @@ def extract_combination_from_LI(samples_LI, p): spec_params['gamma4']=spec_param_array[3] eos_base = EOSManager.EOSLindblomSpectralSoundSpeedVersusPressure(name=eos_name,spec_params=spec_params,use_lal_spec_eos=not opts.no_use_lal_eos) my_eos = eos_base - elif opts.eos_param == 'PP' and len(spec_param_array >=4): + elif opts.eos_param == 'PP' and len(spec_param_array) >= 4: spec_params ={} spec_params['logP1'] = spec_param_array[0] spec_params['gamma1'] = spec_param_array[1] @@ -3646,4 +3646,3 @@ def parse_corr_params(my_str): print(" Failed to generate corner for ", extra_plot_coord_names[indx]) sys.exit(0) - From a2d10ac3dac652c63344c3fa2d262e9d07c5ce48 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Fri, 28 Aug 2026 17:38:28 -0400 Subject: [PATCH 09/16] CIP: expose reviewed lalsim EOS files --- CHANGES.rst | 5 +++-- ...nstructIntrinsicPosterior_GenericCoordinates.py | 14 ++++++++++++++ .../Code/test/test_lalsim_eos_compat.py | 6 ++++++ docs/eos-interface-contract.md | 7 +++++++ 4 files changed, 30 insertions(+), 2 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index fa9dd865e..6953d1327 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -2,8 +2,9 @@ -------------------------------- - EOSManager supports released and reviewed multipart/multibranch LALSimulation family interfaces while retaining the existing scalar and NuclearMatter-Backend - sequence contracts. Fixed-EOS CIP can select a LAL family with - ``--using-eos-branch``. Native NMB/PCA files retain the ``nmbseq:`` primary-branch + sequence contracts. Fixed-EOS CIP can load a reviewed table with + ``lalsim_file:`` and select a LAL family with ``--using-eos-branch``. Native + NMB/PCA files retain the ``nmbseq:`` primary-branch contract, with disconnected stable runs split before interpolation. The CIP piecewise-polytrope constructor keyword and four-parameter row guards are corrected for Kedia-style workflows. diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py index 7cb98d2cd..4a463d9dd 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py @@ -336,6 +336,8 @@ def extract_combination_from_LI(samples_LI, p): parser.add_argument("--tabular-eos-order-statistic",type=str,default=None,help="Order statistic to use. Options will include R1p4, LambdaTildeQ1, and ...}") parser.add_argument("--using-eos", type=str, default=None, help="Name of EOS. Fit parameter list should physically use lambda1, lambda2 information (but need not). If starts with 'file:', uses a filename with EOS parameters ") parser.add_argument("--using-eos-branch", type=int, default=None, help="Select one stable LALSimulation family branch while preserving the fixed-EOS lambda_from_m(m) interface. Required for an explicitly chosen twin-star branch; not applicable to the primary-branch nmbseq v1 contract.") +parser.add_argument("--using-eos-dirty-phase-transitions", action='store_true', help="With --using-eos lalsim_file:, request the reviewed LALSimulation correction for numerically imperfect pressure plateaus.") +parser.add_argument("--using-eos-extended-family", action='store_true', help="With --using-eos lalsim_file:, build the reviewed extended family instead of the PE-oriented minimal M/R/k2 family.") parser.add_argument("--using-eos-index", type=int, default=None, help="Index of EOS parameters in file.") parser.add_argument("--no-use-lal-eos",action='store_true',help="Do not use LAL EOS interface. Used for spectral EOS. Do not use this.") parser.add_argument("--no-matter1", action='store_true', help="Set the lambda parameters to zero (BBH) but return them") @@ -427,6 +429,12 @@ def extract_combination_from_LI(samples_LI, p): sys.exit(0) if opts.using_eos_branch is not None and opts.using_eos is None: raise ValueError("--using-eos-branch also requires --using-eos") +if (opts.using_eos_dirty_phase_transitions or opts.using_eos_extended_family) and ( + opts.using_eos is None or not opts.using_eos.startswith('lalsim_file:')): + raise ValueError( + "--using-eos-dirty-phase-transitions and --using-eos-extended-family " + "require --using-eos lalsim_file:" + ) my_eos=None #option to be used if gridded values not calculated assuming EOS @@ -518,6 +526,12 @@ def extract_combination_from_LI(samples_LI, p): _, seq_fname, seq_indx = eos_name.split(':') my_eos = EOSManager.EOSSequenceSingleIndex(fname=seq_fname, index=int(seq_indx)) + elif eos_name.startswith('lalsim_file:'): + my_eos = EOSManager.EOSLALSimulationFromFile( + fname=eos_name.split(':', 1)[1], + dirty_phase_transitions=opts.using_eos_dirty_phase_transitions, + minimal_family=not opts.using_eos_extended_family, + ) elif 'lal_' in eos_name: eos_name = eos_name.replace('lal_','') my_eos = EOSManager.EOSLALSimulation(name=eos_name) diff --git a/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py index e6f89418d..aeb40f33b 100644 --- a/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py +++ b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py @@ -158,6 +158,12 @@ def test_eosmanager_file_loader_routes_reviewed_phase_transition_api(monkeypatch assert eos.eos == "dirty-eos" assert eos._get_lalsim_family_adapter().number_of_branches == 2 + extended = EOSManager.EOSLALSimulationFromFile( + "extended-format.dat", minimal_family=False + ) + assert fake_lalsim.file_calls[-1] == ("extended-format.dat", 0) + assert fake_lalsim.create_calls[-1] == ("clean-eos", 0) + def test_eosmanager_smoke_with_installed_released_lalsimulation(): from RIFT.physics import EOSManager diff --git a/docs/eos-interface-contract.md b/docs/eos-interface-contract.md index 21bd57621..f686324f5 100644 --- a/docs/eos-interface-contract.md +++ b/docs/eos-interface-contract.md @@ -25,6 +25,13 @@ The fixed-EOS CIP driver exposes the same operation as: --using-eos lal_ --using-eos-branch ``` +Reviewed two- or nine-column tables use: + +```text +--using-eos lalsim_file: [--using-eos-dirty-phase-transitions] + [--using-eos-extended-family] [--using-eos-branch ] +``` + For pseudo-pipe workflows, forward the flag with `--manual-extra-cip-args`. On O4d, Hydra hyperpipe configurations can put it in the post driver's `extra-args` when that driver is the fixed-EOS CIP From b010ff6a5a726f9ddd95712fed0a540b3ee702d9 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Tue, 1 Sep 2026 18:15:18 -0400 Subject: [PATCH 10/16] EOS: fail closed on ambiguous branch helpers --- .../Code/RIFT/physics/EOSManager.py | 35 +++++++++-- .../Code/RIFT/physics/lalsim_eos_compat.py | 47 +++++++++++--- ...ctIntrinsicPosterior_GenericCoordinates.py | 6 +- .../Code/test/test_lalsim_eos_compat.py | 61 ++++++++++++++++++- 4 files changed, 134 insertions(+), 15 deletions(-) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py b/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py index c111ae0e3..9cb16c886 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py @@ -149,7 +149,7 @@ def lambda_from_m_vector(self, m, branch_id=None): ]) return out - def estimate_baryon_mass_from_mg(self,m): + def estimate_baryon_mass_from_mg(self, m, branch_id=None): r""" Estimate m_b = m_g + m_g^2/(R_{1.4}/km) based on https://arxiv.org/pdf/1905.03784.pdf Eq. (6) Note baryon mass can be computed exactly with a TOV solution integral (e.g., Eq. 6.21 of Haensel's book) @@ -157,7 +157,7 @@ def estimate_baryon_mass_from_mg(self,m): but lalsuite doesn't provide access to this low-level info !! This function is only for use when LALEOS is created. Use RePrimAnd's baryon_mass_from_mg preferably for most other purposes!! """ - r1p4 = self.radius_from_m(1.4, branch_id=None) / 1e3 + r1p4 = self.radius_from_m(1.4, branch_id=branch_id) / 1e3 return m + (1./r1p4)*m**2 def pressure_density_on_grid_alternate(self,logrho_grid,enforce_causal=False): @@ -228,7 +228,8 @@ def test_speed_of_sound_causal_builtin(self): else: return True - def test_speed_of_sound_causal(self, test_only_under_mmax=True,fast_test=True): + def test_speed_of_sound_causal( + self, test_only_under_mmax=True, fast_test=True, branch_id=None): """ Test if EOS satisfies speed of sound. Relies on low-level lalsimulation interpolation routines to get v(h) and as such is not very reliable @@ -243,12 +244,13 @@ def test_speed_of_sound_causal(self, test_only_under_mmax=True,fast_test=True): eos = self.eos fam = self.eos_fam # Largest NS provides largest attained central pressure - m_max_SI = self.mMaxMsun*lal.MSUN_SI + family = self._get_lalsim_family_adapter() + m_max_SI = family.maximum_mass(branch_id) if branch_id is not None else self.mMaxMsun*lal.MSUN_SI if not test_only_under_mmax: hmax = lalsim.SimNeutronStarEOSMaxPseudoEnthalpy(eos) else: try: - pmax = self._get_lalsim_family_adapter().central_pressure(m_max_SI) + pmax = family.central_pressure(m_max_SI, branch_id=branch_id) hmax = lalsim.SimNeutronStarEOSPseudoEnthalpyOfPressure(pmax,eos) except: # gatch gsl interpolation errors for example @@ -313,6 +315,29 @@ def lambda_from_m_vector(self, m): return self.lambda_from_m(m) return np.array([self.lambda_from_m(m_here) for m_here in m]) + def estimate_baryon_mass_from_mg(self, m): + """Estimate baryonic mass using this branch's 1.4-Msun radius.""" + try: + return self.source.estimate_baryon_mass_from_mg( + m, branch_id=self.branch_id + ) + except ValueError as exc: + raise ValueError( + "cannot estimate baryonic mass on branch {} because that " + "branch does not contain the 1.4-Msun reference star".format( + self.branch_id + ) + ) from exc + + def test_speed_of_sound_causal( + self, test_only_under_mmax=True, fast_test=True): + """Test causality only up to this branch's maximum-mass star.""" + return self.source.test_speed_of_sound_causal( + test_only_under_mmax=test_only_under_mmax, + fast_test=fast_test, + branch_id=self.branch_id, + ) + def __getattr__(self, name): return getattr(self.source, name) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/physics/lalsim_eos_compat.py b/MonteCarloMarginalizeCode/Code/RIFT/physics/lalsim_eos_compat.py index 39d859e3a..3d89c692b 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/physics/lalsim_eos_compat.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/physics/lalsim_eos_compat.py @@ -13,6 +13,24 @@ class AmbiguousFamilyBranchError(ValueError): """Raised when a mass belongs to more than one stable family branch.""" +def validate_fixed_eos_branch_request(branch_id, eos_spec, using_eos_for_prior=False): + """Fail closed when a scalar fixed-EOS branch request cannot be honored. + + EOS hyperprior plugins construct an EOS after the fixed-EOS setup path has + run. Until that plugin contract carries branch identity explicitly, a + ``branch_id`` here would otherwise be accepted and silently ignored. + """ + if branch_id is None: + return + if using_eos_for_prior: + raise ValueError( + "--using-eos-branch is not supported with --using-eos-for-prior; " + "the EOS hyperprior plugin contract does not carry branch identity" + ) + if eos_spec is None: + raise ValueError("--using-eos-branch also requires --using-eos") + + class LALSimNeutronStarFamilyAdapter: """Version-neutral access to a LALSimulation neutron-star family. @@ -35,14 +53,26 @@ class LALSimNeutronStarFamilyAdapter: "SimNeutronStarFamLoveNumberK2OfMassPerBranch", ) + @classmethod + def _uses_multibranch_api(cls, lalsim_module): + present = [hasattr(lalsim_module, name) for name in cls._MODERN_REQUIRED] + if any(present) and not all(present): + missing = [ + name for name, available in zip(cls._MODERN_REQUIRED, present) + if not available + ] + raise RuntimeError( + "partial reviewed LALSimulation family API; missing symbols: {}" + .format(", ".join(missing)) + ) + return all(present) + def __init__(self, eos, minimal=True, lalsim_module=None): if lalsim_module is None: import lalsimulation as lalsim_module self.lalsim = lalsim_module self.eos = eos - self.is_multibranch_api = all( - hasattr(self.lalsim, name) for name in self._MODERN_REQUIRED - ) + self.is_multibranch_api = self._uses_multibranch_api(self.lalsim) if self.is_multibranch_api: # The reviewed API requires ``min_fam``: 1 selects the PE-oriented # M/R/k2 solver, while 0 also constructs baryonic mass, k3, and k4. @@ -61,9 +91,7 @@ def from_family(cls, family, lalsim_module=None): obj.lalsim = lalsim_module obj.eos = None obj.family = family - obj.is_multibranch_api = all( - hasattr(obj.lalsim, name) for name in cls._MODERN_REQUIRED - ) + obj.is_multibranch_api = cls._uses_multibranch_api(obj.lalsim) return obj @property @@ -157,7 +185,12 @@ def central_pressure(self, mass_si, branch_id=None): modern = getattr( self.lalsim, "SimNeutronStarFamCentralPressureOfMassPerBranch", None ) - if self.is_multibranch_api and modern is not None: + if self.is_multibranch_api: + if modern is None: + raise NotImplementedError( + "reviewed LALSimulation family API lacks branch-specific " + "central pressure" + ) return modern(mass_si, self.family, resolved) return self.lalsim.SimNeutronStarCentralPressure(mass_si, self.family) diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py index 4a463d9dd..663a11e2c 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py @@ -427,8 +427,10 @@ def extract_combination_from_LI(samples_LI, p): except Exception as e: print(" Fail: EOS index out of range:\n ",e) sys.exit(0) -if opts.using_eos_branch is not None and opts.using_eos is None: - raise ValueError("--using-eos-branch also requires --using-eos") +from RIFT.physics.lalsim_eos_compat import validate_fixed_eos_branch_request +validate_fixed_eos_branch_request( + opts.using_eos_branch, opts.using_eos, opts.using_eos_for_prior +) if (opts.using_eos_dirty_phase_transitions or opts.using_eos_extended_family) and ( opts.using_eos is None or not opts.using_eos.startswith('lalsim_file:')): raise ValueError( diff --git a/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py index aeb40f33b..634fc7c25 100644 --- a/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py +++ b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py @@ -7,6 +7,7 @@ from RIFT.physics.lalsim_eos_compat import ( AmbiguousFamilyBranchError, LALSimNeutronStarFamilyAdapter, + validate_fixed_eos_branch_request, ) @@ -94,9 +95,15 @@ def SimNeutronStarFamCentralPressureOfMassPerBranch( class StellarMassMultibranchLALSimulation(MultibranchLALSimulation): bounds = tuple( (lower * lal.MSUN_SI, upper * lal.MSUN_SI) - for lower, upper in MultibranchLALSimulation.bounds + for lower, upper in ((1.0, 2.0), (1.3, 3.0)) ) + def SimNeutronStarEOSPseudoEnthalpyOfPressure(self, pressure, eos): + return pressure + + def SimNeutronStarEOSSpeedOfSoundGeometerized(self, enthalpy, eos): + return 0.5 + def test_released_lalsimulation_uses_one_argument_family_api(): lalsim = LegacyLALSimulation() @@ -131,6 +138,16 @@ def test_reviewed_lalsimulation_uses_minimal_multibranch_api(): assert family.central_pressure(1.75, branch_id=1) == 101.75 +def test_partial_reviewed_api_fails_diagnostically(monkeypatch): + lalsim = MultibranchLALSimulation() + monkeypatch.delattr( + MultibranchLALSimulation, + "SimNeutronStarFamLoveNumberK2OfMassPerBranch", + ) + with pytest.raises(RuntimeError, match="partial reviewed LALSimulation"): + LALSimNeutronStarFamilyAdapter("eos", lalsim_module=lalsim) + + def test_twin_star_mass_requires_an_explicit_branch(): family = LALSimNeutronStarFamilyAdapter( "eos", lalsim_module=MultibranchLALSimulation() @@ -212,6 +229,48 @@ def test_selected_branch_view_preserves_legacy_scalar_consumer_api(monkeypatch): assert primary.lambda_from_m(2.5) == pytest.approx(1e-8) +def test_selected_branch_view_preserves_branch_sensitive_helpers(monkeypatch): + from RIFT.physics import EOSManager + + fake_lalsim = StellarMassMultibranchLALSimulation() + monkeypatch.setattr(EOSManager, "lalsim", fake_lalsim) + secondary = EOSManager.EOSLALSimulationFromFile("twin-star.dat").for_branch(1) + + expected_radius_km = (10.0 + 1.4 * lal.MSUN_SI) / 1e3 + assert secondary.estimate_baryon_mass_from_mg(1.4) == pytest.approx( + 1.4 + 1.4**2 / expected_radius_km + ) + assert secondary.test_speed_of_sound_causal() + + +def test_selected_branch_helpers_fail_closed_when_branch_data_are_missing(monkeypatch): + from RIFT.physics import EOSManager + + no_reference_star = StellarMassMultibranchLALSimulation() + no_reference_star.bounds = tuple( + (lower * lal.MSUN_SI, upper * lal.MSUN_SI) + for lower, upper in MultibranchLALSimulation.bounds + ) + monkeypatch.setattr(EOSManager, "lalsim", no_reference_star) + secondary = EOSManager.EOSLALSimulationFromFile("twin-star.dat").for_branch(1) + with pytest.raises(ValueError, match="does not contain the 1.4-Msun"): + secondary.estimate_baryon_mass_from_mg(1.6) + + monkeypatch.delattr( + MultibranchLALSimulation, + "SimNeutronStarFamCentralPressureOfMassPerBranch", + ) + assert secondary.test_speed_of_sound_causal() is False + + +def test_eos_hyperprior_rejects_fixed_branch_request(): + with pytest.raises(ValueError, match="not supported with --using-eos-for-prior"): + validate_fixed_eos_branch_request(1, "file:eos-draws.dat", True) + + assert validate_fixed_eos_branch_request(1, "lalsim_file:eos.dat") is None + assert validate_fixed_eos_branch_request(None, "file:eos-draws.dat", True) is None + + def test_nmb_sequence_dispatch_and_accessors_remain_compatible(tmp_path): h5py = pytest.importorskip("h5py") from RIFT.physics import EOSManager From 5a8f5de4b6790c7dc06fa5aa0e4348edcb8bca54 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Tue, 1 Sep 2026 18:15:26 -0400 Subject: [PATCH 11/16] tests: add opt-in reviewed lalsim build gate --- .../test_lalsim_eos_reviewed_integration.py | 105 ++++++++++++++++++ docs/eos-interface-contract.md | 32 ++++++ 2 files changed, 137 insertions(+) create mode 100644 MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_reviewed_integration.py diff --git a/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_reviewed_integration.py b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_reviewed_integration.py new file mode 100644 index 000000000..d8069f7b6 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_reviewed_integration.py @@ -0,0 +1,105 @@ +"""Opt-in gate against an actual reviewed LALSimulation build. + +This is intentionally separate from the fake-backed compatibility tests. It +skips ordinary CI unless RIFT_REVIEWED_LALSIM_MANIFEST names a build-generated +manifest and fails closed once the gate is enabled. +""" + +import hashlib +import json +import os +from pathlib import Path +import re + +import numpy as np +import pytest + + +MANIFEST_ENV = "RIFT_REVIEWED_LALSIM_MANIFEST" +REQUIRED_SYMBOLS = ( + "SimNeutronStarEOSFromFileChoiceDirtyPT", + "SimNeutronStarFamNumberOfBranches", + "SimNeutronStarFamMinMassPerBranch", + "SimNeutronStarFamMaxMassPerBranch", + "SimNeutronStarFamRadiusOfMassPerBranch", + "SimNeutronStarFamLoveNumberK2OfMassPerBranch", + "SimNeutronStarFamCentralPressureOfMassPerBranch", +) + + +def _sha256(path): + digest = hashlib.sha256() + with path.open("rb") as stream: + for block in iter(lambda: stream.read(1024 * 1024), b""): + digest.update(block) + return digest.hexdigest() + + +def _load_manifest(): + manifest_name = os.environ.get(MANIFEST_ENV) + if not manifest_name: + pytest.skip( + "real reviewed-LALSimulation gate disabled; set {}".format( + MANIFEST_ENV + ) + ) + manifest_path = Path(manifest_name).resolve() + with manifest_path.open() as stream: + manifest = json.load(stream) + ref = manifest.get("lalsuite_ref", "") + assert re.fullmatch(r"[0-9a-f]{40}", ref), ( + "lalsuite_ref must be the exact 40-character commit built for this job" + ) + return manifest_path, manifest + + +def test_actual_reviewed_lalsimulation_tables(record_property): + import lalsimulation as lalsim + from RIFT.physics import EOSManager + from RIFT.physics.lalsim_eos_compat import AmbiguousFamilyBranchError + + manifest_path, manifest = _load_manifest() + record_property("lalsuite_ref", manifest["lalsuite_ref"]) + record_property( + "lalsimulation_version", + getattr(lalsim, "LALSIMULATION_VERSION", "unknown"), + ) + missing = [name for name in REQUIRED_SYMBOLS if not hasattr(lalsim, name)] + assert not missing, "reviewed LALSimulation symbols missing: {}".format(missing) + + fixtures = manifest.get("fixtures", {}) + assert set(fixtures) == {"two_column", "nine_column", "twin_star"} + loaded = {} + expected_columns = {"two_column": 2, "nine_column": 9, "twin_star": None} + for name in ("two_column", "nine_column", "twin_star"): + fixture = fixtures[name] + path = (manifest_path.parent / fixture["path"]).resolve() + assert path.is_file(), "missing {} fixture: {}".format(name, path) + assert _sha256(path) == fixture["sha256"] + data = np.loadtxt(str(path)) + columns = 1 if data.ndim == 1 else data.shape[1] + if expected_columns[name] is not None: + assert columns == expected_columns[name] + loaded[name] = EOSManager.EOSLALSimulationFromFile( + str(path), + dirty_phase_transitions=bool( + fixture.get("dirty_phase_transitions", False) + ), + ) + assert loaded[name]._get_lalsim_family_adapter().number_of_branches >= 1 + + family = loaded["twin_star"]._get_lalsim_family_adapter() + assert family.number_of_branches >= 2 + overlaps = [] + for left in range(family.number_of_branches): + for right in range(left + 1, family.number_of_branches): + lower = max(family.minimum_mass(left), family.minimum_mass(right)) + upper = min(family.maximum_mass(left), family.maximum_mass(right)) + if lower < upper: + overlaps.append((left, right, 0.5 * (lower + upper))) + assert overlaps, "twin_star fixture has no overlapping stable mass branches" + left, right, mass = overlaps[0] + with pytest.raises(AmbiguousFamilyBranchError): + family.radius(mass) + assert family.radius(mass, branch_id=left) > 0 + assert family.radius(mass, branch_id=right) > 0 diff --git a/docs/eos-interface-contract.md b/docs/eos-interface-contract.md index f686324f5..11263fd70 100644 --- a/docs/eos-interface-contract.md +++ b/docs/eos-interface-contract.md @@ -71,6 +71,38 @@ paper-production parity and treat disconnected branches as a separate model. legacy scalar workflows, or move to the central-enthalpy inference path when marginalization over branch identity is scientifically required. +`--using-eos-branch` is rejected with `--using-eos-for-prior`. The current +EOS-hyperprior plugin protocol returns an EOS but does not return branch +identity, so accepting this combination would silently analyze the wrong +branch. A future multibranch hyperprior must extend that plugin contract before +this restriction can be relaxed. + +## Reviewed-LALSimulation integration gate + +The fake-backed compatibility tests check RIFT's dispatch logic, but do not +certify the reviewed SWIG interface. To run the real-build gate, build the +exact reviewed LALSuite commit, activate that Python environment, and set +`RIFT_REVIEWED_LALSIM_MANIFEST` to a JSON file with this shape: + +```json +{ + "lalsuite_ref": "0123456789abcdef0123456789abcdef01234567", + "fixtures": { + "two_column": {"path": "two-column.dat", "sha256": "..."}, + "nine_column": {"path": "nine-column.dat", "sha256": "..."}, + "twin_star": {"path": "twin-star.dat", "sha256": "...", "dirty_phase_transitions": true} + } +} +``` + +Run +`pytest MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_reviewed_integration.py`. +Paths are relative to the manifest. The ref must be the full 40-character +commit actually built by the job, and every fixture hash is mandatory. Once +the manifest enables the gate, missing modern symbols, malformed provenance, +missing fixtures, wrong column counts, or absence of overlapping twin-star +branches are failures. Ordinary CI skips this private-build gate explicitly. + The drift-sentinel registry should eventually declare an EOS contract group with LALSuite and NuclearMatter-Backend as producers and RIFT/nmb-papers as consumers. The current registry only covers RIFT/hyperpipe operational archive From f8fa943ff9874762b60ffdecd082481566f8c9f4 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Tue, 1 Sep 2026 18:18:45 -0400 Subject: [PATCH 12/16] tests: bind lalsim gate to build provenance --- .../test_lalsim_eos_reviewed_integration.py | 56 +++++++++++++++++-- docs/eos-interface-contract.md | 13 +++-- 2 files changed, 60 insertions(+), 9 deletions(-) diff --git a/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_reviewed_integration.py b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_reviewed_integration.py index d8069f7b6..e510d36b7 100644 --- a/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_reviewed_integration.py +++ b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_reviewed_integration.py @@ -17,6 +17,7 @@ MANIFEST_ENV = "RIFT_REVIEWED_LALSIM_MANIFEST" REQUIRED_SYMBOLS = ( + "SimulationVCSInfo", "SimNeutronStarEOSFromFileChoiceDirtyPT", "SimNeutronStarFamNumberOfBranches", "SimNeutronStarFamMinMassPerBranch", @@ -59,14 +60,25 @@ def test_actual_reviewed_lalsimulation_tables(record_property): from RIFT.physics.lalsim_eos_compat import AmbiguousFamilyBranchError manifest_path, manifest = _load_manifest() + missing = [name for name in REQUIRED_SYMBOLS if not hasattr(lalsim, name)] + assert not missing, "reviewed LALSimulation symbols missing: {}".format(missing) + vcs_info = lalsim.SimulationVCSInfo + assert vcs_info.vcsId == manifest["lalsuite_ref"], ( + "manifest ref {} does not match imported LALSimulation build {}".format( + manifest["lalsuite_ref"], vcs_info.vcsId + ) + ) + assert vcs_info.vcsClean == "CLEAN", ( + "reviewed LALSimulation build has uncommitted source modifications: {}" + .format(vcs_info.vcsStatus) + ) record_property("lalsuite_ref", manifest["lalsuite_ref"]) record_property( "lalsimulation_version", getattr(lalsim, "LALSIMULATION_VERSION", "unknown"), ) - missing = [name for name in REQUIRED_SYMBOLS if not hasattr(lalsim, name)] - assert not missing, "reviewed LALSimulation symbols missing: {}".format(missing) - + record_property("lalsimulation_vcs_status", vcs_info.vcsStatus) + record_property("lalsimulation_vcs_tag", vcs_info.vcsTag) fixtures = manifest.get("fixtures", {}) assert set(fixtures) == {"two_column", "nine_column", "twin_star"} loaded = {} @@ -88,6 +100,19 @@ def test_actual_reviewed_lalsimulation_tables(record_property): ) assert loaded[name]._get_lalsim_family_adapter().number_of_branches >= 1 + # The reviewed contract changes both the table loader and CreateFamily's + # second argument. Exercise clean/dirty readers and minimal/extended family + # construction on the real nine-column fixture rather than on a fake. + nine_path = (manifest_path.parent / fixtures["nine_column"]["path"]).resolve() + nine_dirty = EOSManager.EOSLALSimulationFromFile( + str(nine_path), dirty_phase_transitions=True + ) + nine_extended = EOSManager.EOSLALSimulationFromFile( + str(nine_path), minimal_family=False + ) + assert nine_dirty._get_lalsim_family_adapter().number_of_branches >= 1 + assert nine_extended._get_lalsim_family_adapter().number_of_branches >= 1 + family = loaded["twin_star"]._get_lalsim_family_adapter() assert family.number_of_branches >= 2 overlaps = [] @@ -101,5 +126,26 @@ def test_actual_reviewed_lalsimulation_tables(record_property): left, right, mass = overlaps[0] with pytest.raises(AmbiguousFamilyBranchError): family.radius(mass) - assert family.radius(mass, branch_id=left) > 0 - assert family.radius(mass, branch_id=right) > 0 + with pytest.raises(ValueError, match="branch_id .* outside"): + family.radius(mass, branch_id=family.number_of_branches) + outside_left = np.nextafter(family.maximum_mass(left), np.inf) + with pytest.raises(ValueError, match="outside stable branch"): + family.radius(outside_left, branch_id=left) + radii = [family.radius(mass, branch_id=branch) for branch in (left, right)] + love = [ + family.love_number_k2(mass, branch_id=branch) + for branch in (left, right) + ] + pressure = [ + family.central_pressure(mass, branch_id=branch) + for branch in (left, right) + ] + tidal_lambda = [ + loaded["twin_star"].lambda_from_m(mass, branch_id=branch) + for branch in (left, right) + ] + assert all(value > 0 for value in radii + love + pressure + tidal_lambda) + assert not np.isclose(radii[0], radii[1], rtol=1e-10, atol=0) + assert not np.isclose(love[0], love[1], rtol=1e-10, atol=0) + assert not np.isclose(pressure[0], pressure[1], rtol=1e-10, atol=0) + assert not np.isclose(tidal_lambda[0], tidal_lambda[1], rtol=1e-10, atol=0) diff --git a/docs/eos-interface-contract.md b/docs/eos-interface-contract.md index 11263fd70..168aa0d1a 100644 --- a/docs/eos-interface-contract.md +++ b/docs/eos-interface-contract.md @@ -98,10 +98,15 @@ exact reviewed LALSuite commit, activate that Python environment, and set Run `pytest MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_reviewed_integration.py`. Paths are relative to the manifest. The ref must be the full 40-character -commit actually built by the job, and every fixture hash is mandatory. Once -the manifest enables the gate, missing modern symbols, malformed provenance, -missing fixtures, wrong column counts, or absence of overlapping twin-star -branches are failures. Ordinary CI skips this private-build gate explicitly. +commit actually built by the job: the gate requires it to equal +`lalsimulation.SimulationVCSInfo.vcsId` and requires a clean VCS build. Every +fixture hash is mandatory. Once the manifest enables the gate, missing modern +symbols, malformed or mismatched build provenance, missing fixtures, wrong +column counts, or absence of distinct overlapping twin-star solutions are +failures. The gate exercises clean and phase-transition-correcting readers, +minimal and extended family construction, and branch-indexed radius, Love +number, central pressure, and tidal deformability. Ordinary CI skips this +private-build gate explicitly. The drift-sentinel registry should eventually declare an EOS contract group with LALSuite and NuclearMatter-Backend as producers and RIFT/nmb-papers as From 56603c87c41f8f459e6e8803e30e9ae50d4253a4 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Tue, 1 Sep 2026 20:58:51 -0400 Subject: [PATCH 13/16] EOSManager: use reviewed phase-transition SWIG API --- .../Code/RIFT/physics/EOSManager.py | 137 ++++++++++++------ .../Code/RIFT/physics/lalsim_eos_compat.py | 106 ++++++++------ ...ctIntrinsicPosterior_GenericCoordinates.py | 2 +- .../Code/test/test_lalsim_eos_compat.py | 136 +++++++++++------ .../test_lalsim_eos_reviewed_integration.py | 44 ++++-- docs/eos-interface-contract.md | 15 +- 6 files changed, 297 insertions(+), 143 deletions(-) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py b/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py index 9cb16c886..7a2c3583e 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py @@ -77,10 +77,14 @@ def __init__(self,name=None): self.eos_fam = None return None - def _set_lalsim_family(self, minimal=True): - """Create and cache a released-or-multibranch LAL family.""" - self._lalsim_family_adapter = create_family( - self.eos, minimal=minimal, lalsim_module=lalsim + def _set_lalsim_family(self, minimal=True, reviewed_multibranch=False, + log_pressure_min=None): + """Create and cache a released-or-multibranch LAL family.""" + self._lalsim_reviewed_multibranch = bool(reviewed_multibranch) + self._lalsim_family_adapter = create_family( + self.eos, minimal=minimal, lalsim_module=lalsim, + reviewed_multibranch=self._lalsim_reviewed_multibranch, + log_pressure_min=log_pressure_min, ) self.eos_fam = self._lalsim_family_adapter.family self.mMaxMsun = self._lalsim_family_adapter.maximum_mass() / lal.MSUN_SI @@ -89,8 +93,11 @@ def _set_lalsim_family(self, minimal=True): def _get_lalsim_family_adapter(self): adapter = getattr(self, "_lalsim_family_adapter", None) if adapter is None or adapter.family is not self.eos_fam: - adapter = LALSimNeutronStarFamilyAdapter.from_family( - self.eos_fam, lalsim_module=lalsim + adapter = LALSimNeutronStarFamilyAdapter.from_family( + self.eos_fam, lalsim_module=lalsim, + reviewed_multibranch=getattr( + self, "_lalsim_reviewed_multibranch", False + ), ) self._lalsim_family_adapter = adapter return adapter @@ -244,21 +251,38 @@ def test_speed_of_sound_causal( eos = self.eos fam = self.eos_fam # Largest NS provides largest attained central pressure - family = self._get_lalsim_family_adapter() - m_max_SI = family.maximum_mass(branch_id) if branch_id is not None else self.mMaxMsun*lal.MSUN_SI - if not test_only_under_mmax: - hmax = lalsim.SimNeutronStarEOSMaxPseudoEnthalpy(eos) + family = self._get_lalsim_family_adapter() + m_max_SI = family.maximum_mass(branch_id) if branch_id is not None else self.mMaxMsun*lal.MSUN_SI + if not test_only_under_mmax: + if getattr(self, "_lalsim_reviewed_multibranch", False): + hmax = ( + lalsim.SimNeutronStarEOSMultiPartsMaxPseudoEnthalpy(eos) + ) + else: + hmax = lalsim.SimNeutronStarEOSMaxPseudoEnthalpy(eos) else: try: pmax = family.central_pressure(m_max_SI, branch_id=branch_id) - hmax = lalsim.SimNeutronStarEOSPseudoEnthalpyOfPressure(pmax,eos) + if getattr(self, "_lalsim_reviewed_multibranch", False): + hmax = lalsim.SimNeutronStarEOSMultiPartsPseudoEnthalpyOfPressure( + pmax, eos + ) + else: + hmax = lalsim.SimNeutronStarEOSPseudoEnthalpyOfPressure(pmax,eos) except: # gatch gsl interpolation errors for example return False if fast_test: # https://git.ligo.org/lscsoft/lalsuite/blob/lalinference_o2/lalinference/src/LALInference.c#L2513 try: - vsmax = lalsim.SimNeutronStarEOSSpeedOfSoundGeometerized(hmax, eos) + if getattr(self, "_lalsim_reviewed_multibranch", False): + vsmax = ( + lalsim.SimNeutronStarEOSMultiPartsSpeedOfSoundOfPseudoEnthalpy( + hmax, eos + ) / lal.C_SI + ) + else: + vsmax = lalsim.SimNeutronStarEOSSpeedOfSoundGeometerized(hmax, eos) return vsmax <1.1 except: # catch gsl interpolation errors for example @@ -270,7 +294,16 @@ def test_speed_of_sound_causal( # h = np.linspace(0.0001,lalsim.SimNeutronStarEOSMinAcausalPseudoEnthalpy(eos),npts_internal) vs_internal = np.zeros(npts_internal) for indx in np.arange(npts_internal): - vs_internal[indx] = lalsim.SimNeutronStarEOSSpeedOfSoundGeometerized(h[indx],eos) + if getattr(self, "_lalsim_reviewed_multibranch", False): + vs_internal[indx] = ( + lalsim.SimNeutronStarEOSMultiPartsSpeedOfSoundOfPseudoEnthalpy( + h[indx], eos + ) / lal.C_SI + ) + else: + vs_internal[indx] = lalsim.SimNeutronStarEOSSpeedOfSoundGeometerized( + h[indx], eos + ) if rosDebug: print(h[indx], vs_internal[indx]) return not np.any(vs_internal>1.1) # allow buffer, so we have some threshold @@ -357,32 +390,43 @@ def __init__(self,name): class EOSLALSimulationFromFile(EOSConcrete): """Load a released two-column or reviewed nine-column LAL EOS table. - The reviewed LALSimulation reader detects clean phase transitions in both - formats and preserves all thermodynamic columns in the new format. Set - ``dirty_phase_transitions`` to request its opt-in correction of numerically - imperfect pressure plateaus. + When the reviewed interface is installed, this class uses + ``SimNeutronStarEOSFromFilePhaseTransition`` and the matching multipart + family constructor for both formats. That reader always enables its dirty + phase-transition handling; ``dirty_phase_transitions`` remains accepted as + a backward-compatible request but is not a clean/dirty toggle. """ - def __init__(self, fname, name=None, dirty_phase_transitions=False, - skip_family=False, minimal_family=True): + def __init__(self, fname, name=None, dirty_phase_transitions=False, + skip_family=False, minimal_family=True, + family_log_pressure_min=None): self.name = name or os.path.basename(fname) self.fname = fname self.eos = None self.eos_fam = None - dirty_reader = getattr( - lalsim, "SimNeutronStarEOSFromFileChoiceDirtyPT", None - ) - if dirty_phase_transitions: - if dirty_reader is None: - raise NotImplementedError( - "dirty phase-transition correction requires the reviewed " - "LALSimulation multipart EOS interface" - ) - self.eos = dirty_reader(fname, 1) - else: - self.eos = lalsim.SimNeutronStarEOSFromFile(fname) - if not skip_family: - self._set_lalsim_family(minimal=minimal_family) + phase_transition_reader = getattr( + lalsim, "SimNeutronStarEOSFromFilePhaseTransition", None + ) + self._lalsim_reviewed_multibranch = phase_transition_reader is not None + if self._lalsim_reviewed_multibranch: + self.eos = phase_transition_reader(fname) + elif ( + dirty_phase_transitions + or not minimal_family + or family_log_pressure_min is not None + ): + raise NotImplementedError( + "phase-transition, extended-family, and pressure-floor options " + "require the reviewed LALSimulation multipart EOS interface" + ) + else: + self.eos = lalsim.SimNeutronStarEOSFromFile(fname) + if not skip_family: + self._set_lalsim_family( + minimal=minimal_family, + reviewed_multibranch=self._lalsim_reviewed_multibranch, + log_pressure_min=family_log_pressure_min, + ) else: self.mMaxMsun = None @@ -540,7 +584,7 @@ def eos_ls(self): eos_fname = "./" +eos_name + "_geom.dat" # assume write acces np.savetxt(eos_fname, np.transpose((press, edens)), delimiter='\t') eos = lalsim.SimNeutronStarEOSFromFile(eos_fname) - family_adapter = create_family(eos) + family_adapter = create_family(eos, lalsim_module=lalsim) fam = family_adapter.family else: @@ -1198,16 +1242,21 @@ def int_func(x_prime): ### # Les-like -def make_mr_lambda_lal(eos, n_bins=100, branch_id=None): +def make_mr_lambda_lal(eos, n_bins=100, branch_id=None, + reviewed_multibranch=False): ''' Construct mass-radius curve from EOS Based on modern code resources (https://git.ligo.org/publications/gw170817/bns-eos/blob/master/scripts/eos-params.py) which access low-level structures - ``branch_id`` is optional for released/single-branch LALSimulation. It is - required for a multibranch family so an overlapping twin-star interval is - never collapsed silently. + ``branch_id`` is optional for released/single-branch LALSimulation. It is + required for a multibranch family so an overlapping twin-star interval is + never collapsed silently. Set ``reviewed_multibranch`` only for an EOS + returned by ``SimNeutronStarEOSFromFilePhaseTransition``. ''' - family = create_family(eos) + family = create_family( + eos, lalsim_module=lalsim, + reviewed_multibranch=reviewed_multibranch + ) if family.number_of_branches > 1 and branch_id is None: raise ValueError( "multibranch LAL family requires branch_id; use " @@ -1229,9 +1278,11 @@ def make_mr_lambda_lal(eos, n_bins=100, branch_id=None): return mrL_dat -def make_mr_lambda_lal_branches(eos, n_bins=100): - """Return ``{branch_id: [M, R, Lambda]}`` for every stable LAL branch.""" - family = create_family(eos) +def make_mr_lambda_lal_branches(eos, n_bins=100): + """Return ``{branch_id: [M, R, Lambda]}`` for every stable LAL branch.""" + family = create_family( + eos, lalsim_module=lalsim, reviewed_multibranch=True + ) return { branch_id: _make_mr_lambda_for_family(family, n_bins, branch_id) for branch_id in range(family.number_of_branches) @@ -1262,7 +1313,7 @@ def make_mr_lambda(eos,use_lal=False): if use_lal: make_mr_lambda_lal(eos) - family = create_family(eos) + family = create_family(eos, lalsim_module=lalsim) fam = family.family r_cut = 40 # Some EOS we consider for PE purposes will have very large radius! diff --git a/MonteCarloMarginalizeCode/Code/RIFT/physics/lalsim_eos_compat.py b/MonteCarloMarginalizeCode/Code/RIFT/physics/lalsim_eos_compat.py index 3d89c692b..40b32388b 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/physics/lalsim_eos_compat.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/physics/lalsim_eos_compat.py @@ -41,49 +41,76 @@ class LALSimNeutronStarFamilyAdapter: minimal: On the reviewed API, request the fast family containing only mass, radius, and k2. Released APIs do not have this argument and ignore it. + reviewed_multibranch: + True only when ``eos`` came from + ``SimNeutronStarEOSFromFilePhaseTransition``. Reviewed and legacy SWIG + objects coexist in the same module, so module-level symbol detection + cannot safely choose the family constructor. + log_pressure_min: + Optional lower log-central-pressure bound for the reviewed + ``CreateSimNeutronStarFamilyPTWithPcmin`` constructor. lalsim_module: Dependency-injection hook used by the interface contract tests. """ - _MODERN_REQUIRED = ( + _REVIEWED_REQUIRED = ( + "CreateSimNeutronStarFamilyPT", "SimNeutronStarFamNumberOfBranches", - "SimNeutronStarFamMinMassPerBranch", - "SimNeutronStarFamMaxMassPerBranch", - "SimNeutronStarFamRadiusOfMassPerBranch", - "SimNeutronStarFamLoveNumberK2OfMassPerBranch", + "SimNeutronStarFamBranchMinMass", + "SimNeutronStarFamBranchMaxMass", + "SimNeutronStarFamBranchRadius", + "SimNeutronStarFamBranchLoveNumberK2", + "SimNeutronStarFamBranchCentralPressure", ) @classmethod - def _uses_multibranch_api(cls, lalsim_module): - present = [hasattr(lalsim_module, name) for name in cls._MODERN_REQUIRED] - if any(present) and not all(present): + def _require_reviewed_api(cls, lalsim_module): + """Validate the public ``stable_dev-TOV-headers`` SWIG surface.""" + present = [hasattr(lalsim_module, name) for name in cls._REVIEWED_REQUIRED] + if not all(present): missing = [ - name for name, available in zip(cls._MODERN_REQUIRED, present) + name for name, available in zip(cls._REVIEWED_REQUIRED, present) if not available ] raise RuntimeError( - "partial reviewed LALSimulation family API; missing symbols: {}" + "reviewed LALSimulation phase-transition API is incomplete; " + "missing symbols: {}" .format(", ".join(missing)) ) - return all(present) + return True - def __init__(self, eos, minimal=True, lalsim_module=None): + def __init__(self, eos, minimal=True, lalsim_module=None, + reviewed_multibranch=False, log_pressure_min=None): if lalsim_module is None: import lalsimulation as lalsim_module self.lalsim = lalsim_module self.eos = eos - self.is_multibranch_api = self._uses_multibranch_api(self.lalsim) + self.is_multibranch_api = bool(reviewed_multibranch) if self.is_multibranch_api: + self._require_reviewed_api(self.lalsim) # The reviewed API requires ``min_fam``: 1 selects the PE-oriented # M/R/k2 solver, while 0 also constructs baryonic mass, k3, and k4. - self.family = self.lalsim.CreateSimNeutronStarFamily( - eos, int(bool(minimal)) - ) + if log_pressure_min is None: + self.family = self.lalsim.CreateSimNeutronStarFamilyPT( + eos, int(bool(minimal)) + ) + else: + constructor = getattr( + self.lalsim, "CreateSimNeutronStarFamilyPTWithPcmin", None + ) + if constructor is None: + raise NotImplementedError( + "reviewed LALSimulation build does not expose " + "CreateSimNeutronStarFamilyPTWithPcmin" + ) + self.family = constructor( + eos, int(bool(minimal)), float(log_pressure_min) + ) else: self.family = self.lalsim.CreateSimNeutronStarFamily(eos) @classmethod - def from_family(cls, family, lalsim_module=None): + def from_family(cls, family, lalsim_module=None, reviewed_multibranch=False): """Wrap an already-created family (mainly for tests and transition code).""" if lalsim_module is None: import lalsimulation as lalsim_module @@ -91,7 +118,9 @@ def from_family(cls, family, lalsim_module=None): obj.lalsim = lalsim_module obj.eos = None obj.family = family - obj.is_multibranch_api = cls._uses_multibranch_api(obj.lalsim) + obj.is_multibranch_api = bool(reviewed_multibranch) + if obj.is_multibranch_api: + cls._require_reviewed_api(obj.lalsim) return obj @property @@ -103,13 +132,10 @@ def number_of_branches(self): def minimum_mass(self, branch_id=None): if self.is_multibranch_api: if branch_id is None: - fn = getattr(self.lalsim, "SimNeutronStarFamMinMass", None) - if fn is not None: - return fn(self.family) return min(self.minimum_mass(k) for k in range(self.number_of_branches)) self._validate_branch_id(branch_id) - return self.lalsim.SimNeutronStarFamMinMassPerBranch( - self.family, int(branch_id) + return self.lalsim.SimNeutronStarFamBranchMinMass( + int(branch_id), self.family ) self._validate_legacy_branch_id(branch_id) return self.lalsim.SimNeutronStarFamMinimumMass(self.family) @@ -117,13 +143,10 @@ def minimum_mass(self, branch_id=None): def maximum_mass(self, branch_id=None): if self.is_multibranch_api: if branch_id is None: - fn = getattr(self.lalsim, "SimNeutronStarFamMaxMass", None) - if fn is not None: - return fn(self.family) return max(self.maximum_mass(k) for k in range(self.number_of_branches)) self._validate_branch_id(branch_id) - return self.lalsim.SimNeutronStarFamMaxMassPerBranch( - self.family, int(branch_id) + return self.lalsim.SimNeutronStarFamBranchMaxMass( + int(branch_id), self.family ) self._validate_legacy_branch_id(branch_id) return self.lalsim.SimNeutronStarMaximumMass(self.family) @@ -167,31 +190,25 @@ def resolve_branch(self, mass_si, branch_id=None): def radius(self, mass_si, branch_id=None): resolved = self.resolve_branch(mass_si, branch_id) if self.is_multibranch_api: - return self.lalsim.SimNeutronStarFamRadiusOfMassPerBranch( - mass_si, self.family, resolved + return self.lalsim.SimNeutronStarFamBranchRadius( + mass_si, resolved, self.family ) return self.lalsim.SimNeutronStarRadius(mass_si, self.family) def love_number_k2(self, mass_si, branch_id=None): resolved = self.resolve_branch(mass_si, branch_id) if self.is_multibranch_api: - return self.lalsim.SimNeutronStarFamLoveNumberK2OfMassPerBranch( - mass_si, self.family, resolved + return self.lalsim.SimNeutronStarFamBranchLoveNumberK2( + mass_si, resolved, self.family ) return self.lalsim.SimNeutronStarLoveNumberK2(mass_si, self.family) def central_pressure(self, mass_si, branch_id=None): resolved = self.resolve_branch(mass_si, branch_id) - modern = getattr( - self.lalsim, "SimNeutronStarFamCentralPressureOfMassPerBranch", None - ) if self.is_multibranch_api: - if modern is None: - raise NotImplementedError( - "reviewed LALSimulation family API lacks branch-specific " - "central pressure" - ) - return modern(mass_si, self.family, resolved) + return self.lalsim.SimNeutronStarFamBranchCentralPressure( + mass_si, resolved, self.family + ) return self.lalsim.SimNeutronStarCentralPressure(mass_si, self.family) def _validate_branch_id(self, branch_id): @@ -209,8 +226,11 @@ def _validate_legacy_branch_id(branch_id): raise ValueError("released LALSimulation family exposes only branch 0") -def create_family(eos, minimal=True, lalsim_module=None): +def create_family(eos, minimal=True, lalsim_module=None, + reviewed_multibranch=False, log_pressure_min=None): """Return a :class:`LALSimNeutronStarFamilyAdapter` for ``eos``.""" return LALSimNeutronStarFamilyAdapter( - eos, minimal=minimal, lalsim_module=lalsim_module + eos, minimal=minimal, lalsim_module=lalsim_module, + reviewed_multibranch=reviewed_multibranch, + log_pressure_min=log_pressure_min, ) diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py index 663a11e2c..d61b40e48 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py @@ -336,7 +336,7 @@ def extract_combination_from_LI(samples_LI, p): parser.add_argument("--tabular-eos-order-statistic",type=str,default=None,help="Order statistic to use. Options will include R1p4, LambdaTildeQ1, and ...}") parser.add_argument("--using-eos", type=str, default=None, help="Name of EOS. Fit parameter list should physically use lambda1, lambda2 information (but need not). If starts with 'file:', uses a filename with EOS parameters ") parser.add_argument("--using-eos-branch", type=int, default=None, help="Select one stable LALSimulation family branch while preserving the fixed-EOS lambda_from_m(m) interface. Required for an explicitly chosen twin-star branch; not applicable to the primary-branch nmbseq v1 contract.") -parser.add_argument("--using-eos-dirty-phase-transitions", action='store_true', help="With --using-eos lalsim_file:, request the reviewed LALSimulation correction for numerically imperfect pressure plateaus.") +parser.add_argument("--using-eos-dirty-phase-transitions", action='store_true', help="Compatibility flag for --using-eos lalsim_file:. The reviewed PhaseTransition reader always enables its dirty-phase-transition handling.") parser.add_argument("--using-eos-extended-family", action='store_true', help="With --using-eos lalsim_file:, build the reviewed extended family instead of the PE-oriented minimal M/R/k2 family.") parser.add_argument("--using-eos-index", type=int, default=None, help="Index of EOS parameters in file.") parser.add_argument("--no-use-lal-eos",action='store_true',help="Do not use LAL EOS interface. Used for spectral EOS. Do not use this.") diff --git a/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py index 634fc7c25..0770f3422 100644 --- a/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py +++ b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py @@ -20,10 +20,6 @@ def SimNeutronStarEOSFromFile(self, fname): self.file_calls.append((fname, 0)) return "clean-eos" - def SimNeutronStarEOSFromFileChoiceDirtyPT(self, fname, dirty): - self.file_calls.append((fname, dirty)) - return "dirty-eos" - def CreateSimNeutronStarFamily(self, eos): self.create_calls.append((eos,)) return "legacy-family" @@ -51,44 +47,34 @@ def __init__(self): self.create_calls = [] self.file_calls = [] - def SimNeutronStarEOSFromFile(self, fname): - self.file_calls.append((fname, 0)) - return "clean-eos" + def SimNeutronStarEOSFromFilePhaseTransition(self, fname): + self.file_calls.append((fname,)) + return "multipart-eos" - def SimNeutronStarEOSFromFileChoiceDirtyPT(self, fname, dirty): - self.file_calls.append((fname, dirty)) - return "dirty-eos" - - def CreateSimNeutronStarFamily(self, eos, min_fam): + def CreateSimNeutronStarFamilyPT(self, eos, min_fam): self.create_calls.append((eos, min_fam)) return "multibranch-family" + def CreateSimNeutronStarFamilyPTWithPcmin(self, eos, min_fam, log_pc_min): + self.create_calls.append((eos, min_fam, log_pc_min)) + return "multibranch-family-pcmin" + def SimNeutronStarFamNumberOfBranches(self, family): return len(self.bounds) - def SimNeutronStarFamMinMassPerBranch(self, family, branch_id): + def SimNeutronStarFamBranchMinMass(self, branch_id, family): return self.bounds[branch_id][0] - def SimNeutronStarFamMaxMassPerBranch(self, family, branch_id): + def SimNeutronStarFamBranchMaxMass(self, branch_id, family): return self.bounds[branch_id][1] - def SimNeutronStarFamMinMass(self, family): - return min(x[0] for x in self.bounds) - - def SimNeutronStarFamMaxMass(self, family): - return max(x[1] for x in self.bounds) - - def SimNeutronStarFamRadiusOfMassPerBranch(self, mass, family, branch_id): + def SimNeutronStarFamBranchRadius(self, mass, branch_id, family): return 10.0 * branch_id + mass - def SimNeutronStarFamLoveNumberK2OfMassPerBranch( - self, mass, family, branch_id - ): + def SimNeutronStarFamBranchLoveNumberK2(self, mass, branch_id, family): return branch_id + 0.1 * mass - def SimNeutronStarFamCentralPressureOfMassPerBranch( - self, mass, family, branch_id - ): + def SimNeutronStarFamBranchCentralPressure(self, mass, branch_id, family): return 100.0 * branch_id + mass @@ -98,11 +84,39 @@ class StellarMassMultibranchLALSimulation(MultibranchLALSimulation): for lower, upper in ((1.0, 2.0), (1.3, 3.0)) ) - def SimNeutronStarEOSPseudoEnthalpyOfPressure(self, pressure, eos): + def SimNeutronStarEOSMultiPartsPseudoEnthalpyOfPressure(self, pressure, eos): return pressure - def SimNeutronStarEOSSpeedOfSoundGeometerized(self, enthalpy, eos): - return 0.5 + def SimNeutronStarEOSMultiPartsMaxPseudoEnthalpy(self, eos): + return 1.0 + + def SimNeutronStarEOSMultiPartsSpeedOfSoundOfPseudoEnthalpy( + self, enthalpy, eos + ): + return 0.5 * lal.C_SI + + +class CoexistingLALSimulation(MultibranchLALSimulation): + """Reviewed builds retain the released EOS/family object family too.""" + + def CreateSimNeutronStarFamily(self, eos): + self.create_calls.append(("legacy", eos)) + return "legacy-family" + + def SimNeutronStarFamMinimumMass(self, family): + return 1.0 + + def SimNeutronStarMaximumMass(self, family): + return 3.0 + + def SimNeutronStarRadius(self, mass, family): + return 10.0 + mass + + def SimNeutronStarLoveNumberK2(self, mass, family): + return 0.1 * mass + + def SimNeutronStarCentralPressure(self, mass, family): + return 100.0 * mass def test_released_lalsimulation_uses_one_argument_family_api(): @@ -121,10 +135,23 @@ def test_released_lalsimulation_uses_one_argument_family_api(): family.radius(2.0, branch_id=1) +def test_eosmanager_file_loader_keeps_released_api_fallback(monkeypatch): + from RIFT.physics import EOSManager + + fake_lalsim = LegacyLALSimulation() + monkeypatch.setattr(EOSManager, "lalsim", fake_lalsim) + eos = EOSManager.EOSLALSimulationFromFile("released-format.dat") + + assert fake_lalsim.file_calls == [("released-format.dat", 0)] + assert fake_lalsim.create_calls == [("clean-eos",)] + assert eos._get_lalsim_family_adapter().number_of_branches == 1 + + def test_reviewed_lalsimulation_uses_minimal_multibranch_api(): lalsim = MultibranchLALSimulation() family = LALSimNeutronStarFamilyAdapter( - "eos", minimal=True, lalsim_module=lalsim + "eos", minimal=True, lalsim_module=lalsim, + reviewed_multibranch=True, ) assert lalsim.create_calls == [("eos", 1)] @@ -138,19 +165,33 @@ def test_reviewed_lalsimulation_uses_minimal_multibranch_api(): assert family.central_pressure(1.75, branch_id=1) == 101.75 +def test_reviewed_module_keeps_legacy_family_dispatch_explicit(): + lalsim = CoexistingLALSimulation() + family = LALSimNeutronStarFamilyAdapter( + "legacy-eos", lalsim_module=lalsim, reviewed_multibranch=False + ) + + assert lalsim.create_calls == [("legacy", "legacy-eos")] + assert family.number_of_branches == 1 + assert family.radius(2.0) == 12.0 + + def test_partial_reviewed_api_fails_diagnostically(monkeypatch): lalsim = MultibranchLALSimulation() monkeypatch.delattr( MultibranchLALSimulation, - "SimNeutronStarFamLoveNumberK2OfMassPerBranch", + "SimNeutronStarFamBranchLoveNumberK2", ) - with pytest.raises(RuntimeError, match="partial reviewed LALSimulation"): - LALSimNeutronStarFamilyAdapter("eos", lalsim_module=lalsim) + with pytest.raises(RuntimeError, match="phase-transition API is incomplete"): + LALSimNeutronStarFamilyAdapter( + "eos", lalsim_module=lalsim, reviewed_multibranch=True + ) def test_twin_star_mass_requires_an_explicit_branch(): family = LALSimNeutronStarFamilyAdapter( - "eos", lalsim_module=MultibranchLALSimulation() + "eos", lalsim_module=MultibranchLALSimulation(), + reviewed_multibranch=True, ) with pytest.raises(AmbiguousFamilyBranchError, match=r"branches \[0, 1\]"): @@ -170,16 +211,22 @@ def test_eosmanager_file_loader_routes_reviewed_phase_transition_api(monkeypatch "new-format.dat", dirty_phase_transitions=True ) - assert fake_lalsim.file_calls == [("new-format.dat", 1)] - assert fake_lalsim.create_calls == [("dirty-eos", 1)] - assert eos.eos == "dirty-eos" + assert fake_lalsim.file_calls == [("new-format.dat",)] + assert fake_lalsim.create_calls == [("multipart-eos", 1)] + assert eos.eos == "multipart-eos" assert eos._get_lalsim_family_adapter().number_of_branches == 2 extended = EOSManager.EOSLALSimulationFromFile( "extended-format.dat", minimal_family=False ) - assert fake_lalsim.file_calls[-1] == ("extended-format.dat", 0) - assert fake_lalsim.create_calls[-1] == ("clean-eos", 0) + assert fake_lalsim.file_calls[-1] == ("extended-format.dat",) + assert fake_lalsim.create_calls[-1] == ("multipart-eos", 0) + + pressure_floor = EOSManager.EOSLALSimulationFromFile( + "pressure-floor.dat", family_log_pressure_min=12.5 + ) + assert pressure_floor.eos_fam == "multibranch-family-pcmin" + assert fake_lalsim.create_calls[-1] == ("multipart-eos", 1, 12.5) def test_eosmanager_smoke_with_installed_released_lalsimulation(): @@ -228,6 +275,10 @@ def test_selected_branch_view_preserves_legacy_scalar_consumer_api(monkeypatch): assert np.isfinite(secondary.lambda_from_m(1.75)) assert primary.lambda_from_m(2.5) == pytest.approx(1e-8) + curves = EOSManager.make_mr_lambda_lal_branches(eos.eos, n_bins=3) + assert set(curves) == {0, 1} + assert curves[0].shape == (3, 3) + def test_selected_branch_view_preserves_branch_sensitive_helpers(monkeypatch): from RIFT.physics import EOSManager @@ -241,6 +292,7 @@ def test_selected_branch_view_preserves_branch_sensitive_helpers(monkeypatch): 1.4 + 1.4**2 / expected_radius_km ) assert secondary.test_speed_of_sound_causal() + assert secondary.test_speed_of_sound_causal(test_only_under_mmax=False) def test_selected_branch_helpers_fail_closed_when_branch_data_are_missing(monkeypatch): @@ -257,8 +309,8 @@ def test_selected_branch_helpers_fail_closed_when_branch_data_are_missing(monkey secondary.estimate_baryon_mass_from_mg(1.6) monkeypatch.delattr( - MultibranchLALSimulation, - "SimNeutronStarFamCentralPressureOfMassPerBranch", + StellarMassMultibranchLALSimulation, + "SimNeutronStarEOSMultiPartsSpeedOfSoundOfPseudoEnthalpy", ) assert secondary.test_speed_of_sound_causal() is False diff --git a/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_reviewed_integration.py b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_reviewed_integration.py index e510d36b7..722101001 100644 --- a/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_reviewed_integration.py +++ b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_reviewed_integration.py @@ -11,6 +11,7 @@ from pathlib import Path import re +import lal import numpy as np import pytest @@ -18,13 +19,17 @@ MANIFEST_ENV = "RIFT_REVIEWED_LALSIM_MANIFEST" REQUIRED_SYMBOLS = ( "SimulationVCSInfo", - "SimNeutronStarEOSFromFileChoiceDirtyPT", + "SimNeutronStarEOSFromFilePhaseTransition", + "CreateSimNeutronStarFamilyPT", "SimNeutronStarFamNumberOfBranches", - "SimNeutronStarFamMinMassPerBranch", - "SimNeutronStarFamMaxMassPerBranch", - "SimNeutronStarFamRadiusOfMassPerBranch", - "SimNeutronStarFamLoveNumberK2OfMassPerBranch", - "SimNeutronStarFamCentralPressureOfMassPerBranch", + "SimNeutronStarFamBranchMinMass", + "SimNeutronStarFamBranchMaxMass", + "SimNeutronStarFamBranchRadius", + "SimNeutronStarFamBranchLoveNumberK2", + "SimNeutronStarFamBranchCentralPressure", + "SimNeutronStarEOSMultiPartsMaxPseudoEnthalpy", + "SimNeutronStarEOSMultiPartsPseudoEnthalpyOfPressure", + "SimNeutronStarEOSMultiPartsSpeedOfSoundOfPseudoEnthalpy", ) @@ -100,17 +105,18 @@ def test_actual_reviewed_lalsimulation_tables(record_property): ) assert loaded[name]._get_lalsim_family_adapter().number_of_branches >= 1 - # The reviewed contract changes both the table loader and CreateFamily's - # second argument. Exercise clean/dirty readers and minimal/extended family - # construction on the real nine-column fixture rather than on a fake. + # The reviewed PT reader always enables its phase-transition handling; the + # historical dirty_phase_transitions flag is therefore a compatibility + # alias, not a clean/dirty toggle. Exercise both accepted call forms and + # minimal/extended family construction on a real nine-column fixture. nine_path = (manifest_path.parent / fixtures["nine_column"]["path"]).resolve() - nine_dirty = EOSManager.EOSLALSimulationFromFile( + nine_compat_flag = EOSManager.EOSLALSimulationFromFile( str(nine_path), dirty_phase_transitions=True ) nine_extended = EOSManager.EOSLALSimulationFromFile( str(nine_path), minimal_family=False ) - assert nine_dirty._get_lalsim_family_adapter().number_of_branches >= 1 + assert nine_compat_flag._get_lalsim_family_adapter().number_of_branches >= 1 assert nine_extended._get_lalsim_family_adapter().number_of_branches >= 1 family = loaded["twin_star"]._get_lalsim_family_adapter() @@ -140,10 +146,26 @@ def test_actual_reviewed_lalsimulation_tables(record_property): family.central_pressure(mass, branch_id=branch) for branch in (left, right) ] + enthalpy = lalsim.SimNeutronStarEOSMultiPartsPseudoEnthalpyOfPressure( + pressure[0], loaded["twin_star"].eos + ) + sound_speed_si = ( + lalsim.SimNeutronStarEOSMultiPartsSpeedOfSoundOfPseudoEnthalpy( + enthalpy, loaded["twin_star"].eos + ) + ) + assert np.isfinite(enthalpy) + assert np.isfinite(sound_speed_si) and sound_speed_si > 0 + assert sound_speed_si / lal.C_SI < 1.1 + max_enthalpy = lalsim.SimNeutronStarEOSMultiPartsMaxPseudoEnthalpy( + loaded["twin_star"].eos + ) + assert np.isfinite(max_enthalpy) and max_enthalpy >= enthalpy tidal_lambda = [ loaded["twin_star"].lambda_from_m(mass, branch_id=branch) for branch in (left, right) ] + assert loaded["twin_star"].for_branch(left).test_speed_of_sound_causal() assert all(value > 0 for value in radii + love + pressure + tidal_lambda) assert not np.isclose(radii[0], radii[1], rtol=1e-10, atol=0) assert not np.isclose(love[0], love[1], rtol=1e-10, atol=0) diff --git a/docs/eos-interface-contract.md b/docs/eos-interface-contract.md index 168aa0d1a..50cbb1f6b 100644 --- a/docs/eos-interface-contract.md +++ b/docs/eos-interface-contract.md @@ -32,6 +32,15 @@ Reviewed two- or nine-column tables use: [--using-eos-extended-family] [--using-eos-branch ] ``` +On a reviewed build, RIFT dispatches these files through +`SimNeutronStarEOSFromFilePhaseTransition` and constructs their families with +`CreateSimNeutronStarFamilyPT`. Legacy named, spectral, piecewise-polytrope, +and ordinary file EOS objects continue to use the released one-argument family +constructor: both SWIG object families coexist and are selected explicitly. +The reviewed PT reader always enables its dirty-phase-transition handling, so +`--using-eos-dirty-phase-transitions` remains accepted for command-line +compatibility but does not toggle a second reader mode. + For pseudo-pipe workflows, forward the flag with `--manual-extra-cip-args`. On O4d, Hydra hyperpipe configurations can put it in the post driver's `extra-args` when that driver is the fixed-EOS CIP @@ -103,9 +112,9 @@ commit actually built by the job: the gate requires it to equal fixture hash is mandatory. Once the manifest enables the gate, missing modern symbols, malformed or mismatched build provenance, missing fixtures, wrong column counts, or absence of distinct overlapping twin-star solutions are -failures. The gate exercises clean and phase-transition-correcting readers, -minimal and extended family construction, and branch-indexed radius, Love -number, central pressure, and tidal deformability. Ordinary CI skips this +failures. The gate exercises the always-PT multipart reader through both legacy +flag forms, minimal and extended family construction, and branch-indexed radius, +Love number, central pressure, and tidal deformability. Ordinary CI skips this private-build gate explicitly. The drift-sentinel registry should eventually declare an EOS contract group From f238c8b7fd9da61de59fea0dd8148a9b80675616 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Tue, 1 Sep 2026 21:09:28 -0400 Subject: [PATCH 14/16] tests: harden reviewed lalsim O4c gate --- .gitattributes | 1 + .../Code/RIFT/physics/EOSManager.py | 76 +++-- .../Code/RIFT/physics/lalsim_eos_compat.py | 12 +- ...ctIntrinsicPosterior_GenericCoordinates.py | 1 + .../test/run_lalsim_eos_reviewed_fixture.py | 134 ++++++++ .../Code/test/test_lalsim_eos_compat.py | 82 ++++- .../test_lalsim_eos_reviewed_integration.py | 315 ++++++++++++------ docs/eos-interface-contract.md | 65 +++- 8 files changed, 538 insertions(+), 148 deletions(-) create mode 100644 .gitattributes create mode 100644 MonteCarloMarginalizeCode/Code/test/run_lalsim_eos_reviewed_fixture.py diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 000000000..4f73b9467 --- /dev/null +++ b/.gitattributes @@ -0,0 +1 @@ +MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py whitespace=cr-at-eol diff --git a/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py b/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py index 7a2c3583e..ba1442a5a 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py @@ -255,9 +255,14 @@ def test_speed_of_sound_causal( m_max_SI = family.maximum_mass(branch_id) if branch_id is not None else self.mMaxMsun*lal.MSUN_SI if not test_only_under_mmax: if getattr(self, "_lalsim_reviewed_multibranch", False): - hmax = ( - lalsim.SimNeutronStarEOSMultiPartsMaxPseudoEnthalpy(eos) + max_enthalpy = getattr( + lalsim, + "SimNeutronStarEOSMultiPartsMaxPseudoEnthalpy", + None, ) + if max_enthalpy is None: + return False + hmax = max_enthalpy(eos) else: hmax = lalsim.SimNeutronStarEOSMaxPseudoEnthalpy(eos) else: @@ -387,19 +392,20 @@ def __init__(self,name): return None -class EOSLALSimulationFromFile(EOSConcrete): - """Load a released two-column or reviewed nine-column LAL EOS table. - - When the reviewed interface is installed, this class uses - ``SimNeutronStarEOSFromFilePhaseTransition`` and the matching multipart - family constructor for both formats. That reader always enables its dirty - phase-transition handling; ``dirty_phase_transitions`` remains accepted as - a backward-compatible request but is not a clean/dirty toggle. - """ - +class EOSLALSimulationFromFile(EOSConcrete): + """Load a released two-column or reviewed nine-column LAL EOS table. + + The default preserves the released reader and one-argument family API even + when reviewed symbols coexist in the module. ``phase_transition_aware``, + ``dirty_phase_transitions``, an extended family, or a pressure floor opts + into the reviewed multipart reader. ``family_log_pressure_min`` is the + finite natural logarithm ``ln(Pc / Pa)``. + """ + def __init__(self, fname, name=None, dirty_phase_transitions=False, skip_family=False, minimal_family=True, - family_log_pressure_min=None): + family_log_pressure_min=None, + phase_transition_aware=False): self.name = name or os.path.basename(fname) self.fname = fname self.eos = None @@ -407,18 +413,25 @@ def __init__(self, fname, name=None, dirty_phase_transitions=False, phase_transition_reader = getattr( lalsim, "SimNeutronStarEOSFromFilePhaseTransition", None ) - self._lalsim_reviewed_multibranch = phase_transition_reader is not None - if self._lalsim_reviewed_multibranch: - self.eos = phase_transition_reader(fname) - elif ( - dirty_phase_transitions + if family_log_pressure_min is not None: + family_log_pressure_min = float(family_log_pressure_min) + if not np.isfinite(family_log_pressure_min): + raise ValueError( + "family_log_pressure_min must be finite ln(Pc / Pa)" + ) + self._lalsim_reviewed_multibranch = bool( + phase_transition_aware + or dirty_phase_transitions or not minimal_family or family_log_pressure_min is not None - ): - raise NotImplementedError( - "phase-transition, extended-family, and pressure-floor options " - "require the reviewed LALSimulation multipart EOS interface" - ) + ) + if self._lalsim_reviewed_multibranch: + if phase_transition_reader is None: + raise NotImplementedError( + "the requested multipart EOS family requires the reviewed " + "LALSimulation phase-transition interface" + ) + self.eos = phase_transition_reader(fname) else: self.eos = lalsim.SimNeutronStarEOSFromFile(fname) if not skip_family: @@ -1243,7 +1256,7 @@ def int_func(x_prime): # Les-like def make_mr_lambda_lal(eos, n_bins=100, branch_id=None, - reviewed_multibranch=False): + reviewed_multibranch=False, family_adapter=None): ''' Construct mass-radius curve from EOS Based on modern code resources (https://git.ligo.org/publications/gw170817/bns-eos/blob/master/scripts/eos-params.py) which access low-level structures @@ -1253,9 +1266,9 @@ def make_mr_lambda_lal(eos, n_bins=100, branch_id=None, never collapsed silently. Set ``reviewed_multibranch`` only for an EOS returned by ``SimNeutronStarEOSFromFilePhaseTransition``. ''' - family = create_family( + family = family_adapter or create_family( eos, lalsim_module=lalsim, - reviewed_multibranch=reviewed_multibranch + reviewed_multibranch=reviewed_multibranch, ) if family.number_of_branches > 1 and branch_id is None: raise ValueError( @@ -1278,10 +1291,13 @@ def make_mr_lambda_lal(eos, n_bins=100, branch_id=None, return mrL_dat -def make_mr_lambda_lal_branches(eos, n_bins=100): - """Return ``{branch_id: [M, R, Lambda]}`` for every stable LAL branch.""" - family = create_family( - eos, lalsim_module=lalsim, reviewed_multibranch=True +def make_mr_lambda_lal_branches(eos, n_bins=100, + reviewed_multibranch=False, + family_adapter=None): + """Return branch curves; multipart construction is explicit and opt-in.""" + family = family_adapter or create_family( + eos, lalsim_module=lalsim, + reviewed_multibranch=reviewed_multibranch, ) return { branch_id: _make_mr_lambda_for_family(family, n_bins, branch_id) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/physics/lalsim_eos_compat.py b/MonteCarloMarginalizeCode/Code/RIFT/physics/lalsim_eos_compat.py index 40b32388b..45c446aba 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/physics/lalsim_eos_compat.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/physics/lalsim_eos_compat.py @@ -8,6 +8,8 @@ solutions. """ +import math + class AmbiguousFamilyBranchError(ValueError): """Raised when a mass belongs to more than one stable family branch.""" @@ -47,7 +49,8 @@ class LALSimNeutronStarFamilyAdapter: objects coexist in the same module, so module-level symbol detection cannot safely choose the family constructor. log_pressure_min: - Optional lower log-central-pressure bound for the reviewed + Optional lower natural-log central-pressure bound, ``ln(Pc / Pa)``, for + the reviewed ``CreateSimNeutronStarFamilyPTWithPcmin`` constructor. lalsim_module: Dependency-injection hook used by the interface contract tests. @@ -95,6 +98,11 @@ def __init__(self, eos, minimal=True, lalsim_module=None, eos, int(bool(minimal)) ) else: + log_pressure_min = float(log_pressure_min) + if not math.isfinite(log_pressure_min): + raise ValueError( + "log_pressure_min must be finite ln(Pc / Pa)" + ) constructor = getattr( self.lalsim, "CreateSimNeutronStarFamilyPTWithPcmin", None ) @@ -104,7 +112,7 @@ def __init__(self, eos, minimal=True, lalsim_module=None, "CreateSimNeutronStarFamilyPTWithPcmin" ) self.family = constructor( - eos, int(bool(minimal)), float(log_pressure_min) + eos, int(bool(minimal)), log_pressure_min ) else: self.family = self.lalsim.CreateSimNeutronStarFamily(eos) diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py index d61b40e48..5f56dde0e 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py @@ -533,6 +533,7 @@ def extract_combination_from_LI(samples_LI, p): fname=eos_name.split(':', 1)[1], dirty_phase_transitions=opts.using_eos_dirty_phase_transitions, minimal_family=not opts.using_eos_extended_family, + phase_transition_aware=opts.using_eos_branch is not None, ) elif 'lal_' in eos_name: eos_name = eos_name.replace('lal_','') diff --git a/MonteCarloMarginalizeCode/Code/test/run_lalsim_eos_reviewed_fixture.py b/MonteCarloMarginalizeCode/Code/test/run_lalsim_eos_reviewed_fixture.py new file mode 100644 index 000000000..8a65c555e --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/run_lalsim_eos_reviewed_fixture.py @@ -0,0 +1,134 @@ +"""Hard-timeout subprocess target for externally supplied EOS tables.""" + +import argparse +import importlib.util +import json +from pathlib import Path +import sys +import types + +import numpy as np + + +parser = argparse.ArgumentParser() +parser.add_argument("--fixture", required=True) +parser.add_argument("--columns", type=int, choices=(2, 4, 9), required=True) +parser.add_argument("--arrays", action="store_true") +parser.add_argument("--twin", action="store_true") +parser.add_argument("--extended", action="store_true") +parser.add_argument("--eosmanager", action="store_true") +parser.add_argument("--status", required=True) +args = parser.parse_args() + + +def fail(message): + Path(args.status).write_text(json.dumps({"error": str(message)})) + raise SystemExit(2) + + +def record_uncaught(exc_type, exc, traceback): + Path(args.status).write_text( + json.dumps({"error": "{}: {}".format(exc_type.__name__, exc)}) + ) + + +sys.excepthook = record_uncaught + +data = np.loadtxt(args.fixture, ndmin=2) +columns = data.shape[1] +if columns != args.columns: + fail("column mismatch before native loader") + +if args.eosmanager: + try: + import lalframe # noqa: F401 + except ImportError: + lalframe_stub = types.ModuleType("lalframe") + lalframe_stub.__path__ = [] + lalframe_stub.frread = types.ModuleType("lalframe.frread") + sys.modules["lalframe"] = lalframe_stub + sys.modules["lalframe.frread"] = lalframe_stub.frread + from RIFT.physics import EOSManager + import lal + eos = EOSManager.EOSLALSimulationFromFile( + args.fixture, + phase_transition_aware=True, + minimal_family=not args.extended, + ) + family = eos._get_lalsim_family_adapter() +else: + import lal + import lalsimulation as lalsim + source = ( + Path(__file__).resolve().parents[1] + / "RIFT" / "physics" / "lalsim_eos_compat.py" + ) + spec = importlib.util.spec_from_file_location("reviewed_adapter", str(source)) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + if args.arrays: + if args.columns != 4: + fail("array conversion requires four-column wiki input") + conversion = 1.602176634e32 * lal.G_SI / lal.C_SI**4 + energy = lal.CreateREAL8Vector(len(data)) + pressure = lal.CreateREAL8Vector(len(data)) + energy.data[:] = data[:, 2] * conversion + pressure.data[:] = data[:, 3] * conversion + multipart = lalsim.SimNeutronStarEOSFromArraysPhaseTransition( + energy, pressure + ) + else: + multipart = lalsim.SimNeutronStarEOSFromFilePhaseTransition( + args.fixture + ) + family = module.LALSimNeutronStarFamilyAdapter( + multipart, + minimal=not args.extended, + reviewed_multibranch=True, + lalsim_module=lalsim, + ) + +result = {"branches": family.number_of_branches} +if result["branches"] < 1: + fail("no stable branches") +for branch in range(family.number_of_branches): + lower = family.minimum_mass(branch) + upper = family.maximum_mass(branch) + mass = 0.5 * (lower + upper) + if family.radius(mass, branch_id=branch) <= 0: + fail("nonpositive radius") +if args.eosmanager: + mass_si = 0.5 * (family.minimum_mass(0) + family.maximum_mass(0)) + mass_msun = mass_si / lal.MSUN_SI + if 0 not in eos.branches_for_m(mass_msun): + fail("EOSManager branch lookup lost primary branch") + if not np.isfinite(eos.lambda_from_m(mass_msun, branch_id=0)): + fail("EOSManager returned nonfinite tidal deformability") + selected = eos.for_branch(0) + result["causal_under_branch_max"] = bool( + selected.test_speed_of_sound_causal() + ) + result["causal_full_table"] = bool( + eos.test_speed_of_sound_causal(test_only_under_mmax=False) + ) +if args.twin: + overlaps = [] + for left in range(family.number_of_branches): + for right in range(left + 1, family.number_of_branches): + lower = max(family.minimum_mass(left), family.minimum_mass(right)) + upper = min(family.maximum_mass(left), family.maximum_mass(right)) + if lower < upper: + mass = 0.5 * (lower + upper) + radii = [family.radius(mass, branch_id=x) for x in (left, right)] + love = [ + family.love_number_k2(mass, branch_id=x) + for x in (left, right) + ] + if np.isclose(radii[0], radii[1]) or np.isclose( + love[0], love[1] + ): + fail("branch_id ignored") + overlaps.append((left, right)) + if not overlaps: + fail("no overlapping twin branches") +Path(args.status).write_text(json.dumps(result, sort_keys=True)) diff --git a/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py index 0770f3422..41c822c18 100644 --- a/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py +++ b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py @@ -99,6 +99,10 @@ def SimNeutronStarEOSMultiPartsSpeedOfSoundOfPseudoEnthalpy( class CoexistingLALSimulation(MultibranchLALSimulation): """Reviewed builds retain the released EOS/family object family too.""" + def SimNeutronStarEOSFromFile(self, fname): + self.file_calls.append(("legacy", fname)) + return "legacy-file-eos" + def CreateSimNeutronStarFamily(self, eos): self.create_calls.append(("legacy", eos)) return "legacy-family" @@ -164,6 +168,14 @@ def test_reviewed_lalsimulation_uses_minimal_multibranch_api(): assert family.love_number_k2(1.75, branch_id=1) == pytest.approx(1.175) assert family.central_pressure(1.75, branch_id=1) == 101.75 + with pytest.raises(ValueError, match=r"finite ln\(Pc / Pa\)"): + LALSimNeutronStarFamilyAdapter( + "eos", + lalsim_module=lalsim, + reviewed_multibranch=True, + log_pressure_min=np.inf, + ) + def test_reviewed_module_keeps_legacy_family_dispatch_explicit(): lalsim = CoexistingLALSimulation() @@ -228,6 +240,26 @@ def test_eosmanager_file_loader_routes_reviewed_phase_transition_api(monkeypatch assert pressure_floor.eos_fam == "multibranch-family-pcmin" assert fake_lalsim.create_calls[-1] == ("multipart-eos", 1, 12.5) + for invalid in (np.inf, -np.inf, np.nan): + with pytest.raises(ValueError, match=r"finite ln\(Pc / Pa\)"): + EOSManager.EOSLALSimulationFromFile( + "pressure-floor.dat", family_log_pressure_min=invalid + ) + + +def test_eosmanager_file_loader_preserves_legacy_default(monkeypatch): + from RIFT.physics import EOSManager + + fake_lalsim = CoexistingLALSimulation() + monkeypatch.setattr(EOSManager, "lalsim", fake_lalsim) + eos = EOSManager.EOSLALSimulationFromFile("ordinary-two-column.dat") + + assert fake_lalsim.file_calls == [ + ("legacy", "ordinary-two-column.dat") + ] + assert fake_lalsim.create_calls == [("legacy", "legacy-file-eos")] + assert eos._get_lalsim_family_adapter().number_of_branches == 1 + def test_eosmanager_smoke_with_installed_released_lalsimulation(): from RIFT.physics import EOSManager @@ -262,7 +294,9 @@ def test_selected_branch_view_preserves_legacy_scalar_consumer_api(monkeypatch): fake_lalsim = StellarMassMultibranchLALSimulation() monkeypatch.setattr(EOSManager, "lalsim", fake_lalsim) - eos = EOSManager.EOSLALSimulationFromFile("twin-star.dat") + eos = EOSManager.EOSLALSimulationFromFile( + "twin-star.dat", phase_transition_aware=True + ) primary = eos.for_branch(0) secondary = eos.for_branch(1) @@ -275,17 +309,38 @@ def test_selected_branch_view_preserves_legacy_scalar_consumer_api(monkeypatch): assert np.isfinite(secondary.lambda_from_m(1.75)) assert primary.lambda_from_m(2.5) == pytest.approx(1e-8) - curves = EOSManager.make_mr_lambda_lal_branches(eos.eos, n_bins=3) + curves = EOSManager.make_mr_lambda_lal_branches( + eos.eos, n_bins=3, family_adapter=eos._get_lalsim_family_adapter() + ) assert set(curves) == {0, 1} assert curves[0].shape == (3, 3) + explicit = EOSManager.make_mr_lambda_lal( + eos.eos, n_bins=3, branch_id=1, reviewed_multibranch=True + ) + assert explicit.shape == (3, 3) + + +def test_mr_lambda_branches_helper_keeps_legacy_default(monkeypatch): + from RIFT.physics import EOSManager + + fake_lalsim = CoexistingLALSimulation() + monkeypatch.setattr(EOSManager, "lalsim", fake_lalsim) + curves = EOSManager.make_mr_lambda_lal_branches("legacy-eos", n_bins=3) + + assert set(curves) == {0} + assert curves[0].shape == (3, 3) + assert fake_lalsim.create_calls == [("legacy", "legacy-eos")] + def test_selected_branch_view_preserves_branch_sensitive_helpers(monkeypatch): from RIFT.physics import EOSManager fake_lalsim = StellarMassMultibranchLALSimulation() monkeypatch.setattr(EOSManager, "lalsim", fake_lalsim) - secondary = EOSManager.EOSLALSimulationFromFile("twin-star.dat").for_branch(1) + secondary = EOSManager.EOSLALSimulationFromFile( + "twin-star.dat", phase_transition_aware=True + ).for_branch(1) expected_radius_km = (10.0 + 1.4 * lal.MSUN_SI) / 1e3 assert secondary.estimate_baryon_mass_from_mg(1.4) == pytest.approx( @@ -304,7 +359,9 @@ def test_selected_branch_helpers_fail_closed_when_branch_data_are_missing(monkey for lower, upper in MultibranchLALSimulation.bounds ) monkeypatch.setattr(EOSManager, "lalsim", no_reference_star) - secondary = EOSManager.EOSLALSimulationFromFile("twin-star.dat").for_branch(1) + secondary = EOSManager.EOSLALSimulationFromFile( + "twin-star.dat", phase_transition_aware=True + ).for_branch(1) with pytest.raises(ValueError, match="does not contain the 1.4-Msun"): secondary.estimate_baryon_mass_from_mg(1.6) @@ -315,6 +372,23 @@ def test_selected_branch_helpers_fail_closed_when_branch_data_are_missing(monkey assert secondary.test_speed_of_sound_causal() is False +def test_multipart_full_table_causality_requires_max_enthalpy(monkeypatch): + from RIFT.physics import EOSManager + + fake_lalsim = StellarMassMultibranchLALSimulation() + monkeypatch.setattr(EOSManager, "lalsim", fake_lalsim) + secondary = EOSManager.EOSLALSimulationFromFile( + "twin-star.dat", phase_transition_aware=True + ).for_branch(1) + monkeypatch.delattr( + StellarMassMultibranchLALSimulation, + "SimNeutronStarEOSMultiPartsMaxPseudoEnthalpy", + ) + assert secondary.test_speed_of_sound_causal( + test_only_under_mmax=False + ) is False + + def test_eos_hyperprior_rejects_fixed_branch_request(): with pytest.raises(ValueError, match="not supported with --using-eos-for-prior"): validate_fixed_eos_branch_request(1, "file:eos-draws.dat", True) diff --git a/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_reviewed_integration.py b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_reviewed_integration.py index 722101001..9a0554752 100644 --- a/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_reviewed_integration.py +++ b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_reviewed_integration.py @@ -6,12 +6,15 @@ """ import hashlib +import importlib.util import json import os from pathlib import Path import re +import subprocess +import sys +import tempfile -import lal import numpy as np import pytest @@ -19,6 +22,8 @@ MANIFEST_ENV = "RIFT_REVIEWED_LALSIM_MANIFEST" REQUIRED_SYMBOLS = ( "SimulationVCSInfo", + "SimNeutronStarEOSMultiPartsByName", + "SimNeutronStarEOSFromArraysPhaseTransition", "SimNeutronStarEOSFromFilePhaseTransition", "CreateSimNeutronStarFamilyPT", "SimNeutronStarFamNumberOfBranches", @@ -59,115 +64,231 @@ def _load_manifest(): return manifest_path, manifest -def test_actual_reviewed_lalsimulation_tables(record_property): - import lalsimulation as lalsim - from RIFT.physics import EOSManager - from RIFT.physics.lalsim_eos_compat import AmbiguousFamilyBranchError +def _load_adapter_module(): + """Load the pure adapter without importing RIFT's heavyweight package.""" + source = ( + Path(__file__).resolve().parents[1] + / "RIFT" / "physics" / "lalsim_eos_compat.py" + ) + spec = importlib.util.spec_from_file_location( + "rift_lalsim_eos_compat_gate", str(source) + ) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module - manifest_path, manifest = _load_manifest() + +def _validate_reviewed_build(lalsim, manifest, record_property): missing = [name for name in REQUIRED_SYMBOLS if not hasattr(lalsim, name)] assert not missing, "reviewed LALSimulation symbols missing: {}".format(missing) vcs_info = lalsim.SimulationVCSInfo - assert vcs_info.vcsId == manifest["lalsuite_ref"], ( - "manifest ref {} does not match imported LALSimulation build {}".format( - manifest["lalsuite_ref"], vcs_info.vcsId + assert vcs_info.vcsId == manifest["lalsuite_ref"] + assert vcs_info.vcsClean == "CLEAN" + record_property("lalsuite_ref", vcs_info.vcsId) + record_property("lalsimulation_vcs_status", vcs_info.vcsStatus) + record_property("lalsimulation_vcs_tag", vcs_info.vcsTag) + + +def _run_fixture_subprocess(command): + """Run native parsing with bounded time, output, and returned status.""" + with tempfile.TemporaryDirectory(prefix="rift-reviewed-eos-") as tmpdir: + status = Path(tmpdir) / "status.json" + result = subprocess.run( + command + ["--status", str(status)], + stdout=subprocess.DEVNULL, + stderr=subprocess.DEVNULL, + timeout=60, + check=False, + ) + if status.exists(): + with status.open("rb") as stream: + detail = stream.read(65536).decode("utf-8", errors="replace") + else: + detail = "" + assert result.returncode == 0, detail or ( + "fixture subprocess failed with return code {}".format( + result.returncode ) ) - assert vcs_info.vcsClean == "CLEAN", ( - "reviewed LALSimulation build has uncommitted source modifications: {}" - .format(vcs_info.vcsStatus) - ) - record_property("lalsuite_ref", manifest["lalsuite_ref"]) - record_property( - "lalsimulation_version", - getattr(lalsim, "LALSIMULATION_VERSION", "unknown"), + + +def _run_expected_upstream_crash(command): + result = subprocess.run( + command, + stdout=subprocess.DEVNULL, + stderr=subprocess.DEVNULL, + timeout=60, + check=False, ) - record_property("lalsimulation_vcs_status", vcs_info.vcsStatus) - record_property("lalsimulation_vcs_tag", vcs_info.vcsTag) - fixtures = manifest.get("fixtures", {}) - assert set(fixtures) == {"two_column", "nine_column", "twin_star"} - loaded = {} - expected_columns = {"two_column": 2, "nine_column": 9, "twin_star": None} - for name in ("two_column", "nine_column", "twin_star"): - fixture = fixtures[name] - path = (manifest_path.parent / fixture["path"]).resolve() - assert path.is_file(), "missing {} fixture: {}".format(name, path) - assert _sha256(path) == fixture["sha256"] - data = np.loadtxt(str(path)) - columns = 1 if data.ndim == 1 else data.shape[1] - if expected_columns[name] is not None: - assert columns == expected_columns[name] - loaded[name] = EOSManager.EOSLALSimulationFromFile( - str(path), - dirty_phase_transitions=bool( - fixture.get("dirty_phase_transitions", False) - ), + if result.returncode == 0: + pytest.fail( + "known two-transition native crash is fixed; promote this fixture " + "to a required passing gate" ) - assert loaded[name]._get_lalsim_family_adapter().number_of_branches >= 1 + assert result.returncode in (-11, 139), ( + "expected SIGSEGV/139, got {}".format(result.returncode) + ) + pytest.xfail("upstream reviewed LALSimulation two-transition SIGSEGV") - # The reviewed PT reader always enables its phase-transition handling; the - # historical dirty_phase_transitions flag is therefore a compatibility - # alias, not a clean/dirty toggle. Exercise both accepted call forms and - # minimal/extended family construction on a real nine-column fixture. - nine_path = (manifest_path.parent / fixtures["nine_column"]["path"]).resolve() - nine_compat_flag = EOSManager.EOSLALSimulationFromFile( - str(nine_path), dirty_phase_transitions=True + +def test_fixture_subprocess_status_detail_is_bounded(): + code = ( + "from pathlib import Path; import sys; " + "Path(sys.argv[-1]).write_text('x' * 1000000); raise SystemExit(3)" ) - nine_extended = EOSManager.EOSLALSimulationFromFile( - str(nine_path), minimal_family=False + with pytest.raises(AssertionError) as caught: + _run_fixture_subprocess([sys.executable, "-c", code]) + assert len(str(caught.value)) <= 66000 + + +def test_actual_reviewed_lalsimulation_builtin(record_property): + """Safe acceptance using LALSuite's trusted built-in SLY table.""" + import gc + import lal + import lalsimulation as lalsim + + _, manifest = _load_manifest() + _validate_reviewed_build(lalsim, manifest, record_property) + Adapter = _load_adapter_module().LALSimNeutronStarFamilyAdapter + + multipart_eos = lalsim.SimNeutronStarEOSMultiPartsByName("SLY") + minimal = Adapter( + multipart_eos, + minimal=True, + reviewed_multibranch=True, + lalsim_module=lalsim, ) - assert nine_compat_flag._get_lalsim_family_adapter().number_of_branches >= 1 - assert nine_extended._get_lalsim_family_adapter().number_of_branches >= 1 - - family = loaded["twin_star"]._get_lalsim_family_adapter() - assert family.number_of_branches >= 2 - overlaps = [] - for left in range(family.number_of_branches): - for right in range(left + 1, family.number_of_branches): - lower = max(family.minimum_mass(left), family.minimum_mass(right)) - upper = min(family.maximum_mass(left), family.maximum_mass(right)) - if lower < upper: - overlaps.append((left, right, 0.5 * (lower + upper))) - assert overlaps, "twin_star fixture has no overlapping stable mass branches" - left, right, mass = overlaps[0] - with pytest.raises(AmbiguousFamilyBranchError): - family.radius(mass) - with pytest.raises(ValueError, match="branch_id .* outside"): - family.radius(mass, branch_id=family.number_of_branches) - outside_left = np.nextafter(family.maximum_mass(left), np.inf) - with pytest.raises(ValueError, match="outside stable branch"): - family.radius(outside_left, branch_id=left) - radii = [family.radius(mass, branch_id=branch) for branch in (left, right)] - love = [ - family.love_number_k2(mass, branch_id=branch) - for branch in (left, right) - ] - pressure = [ - family.central_pressure(mass, branch_id=branch) - for branch in (left, right) - ] + assert minimal.number_of_branches >= 1 + mass = 0.5 * (minimal.minimum_mass(0) + minimal.maximum_mass(0)) + radius = minimal.radius(mass, branch_id=0) + love = minimal.love_number_k2(mass, branch_id=0) + pressure = minimal.central_pressure(mass, branch_id=0) enthalpy = lalsim.SimNeutronStarEOSMultiPartsPseudoEnthalpyOfPressure( - pressure[0], loaded["twin_star"].eos + pressure, multipart_eos ) - sound_speed_si = ( - lalsim.SimNeutronStarEOSMultiPartsSpeedOfSoundOfPseudoEnthalpy( - enthalpy, loaded["twin_star"].eos + sound_si = lalsim.SimNeutronStarEOSMultiPartsSpeedOfSoundOfPseudoEnthalpy( + enthalpy, multipart_eos + ) + assert all( + np.isfinite(value) and value > 0 + for value in ( + radius, + love, + pressure, + enthalpy, + sound_si, + lalsim.SimNeutronStarEOSMultiPartsMaxPseudoEnthalpy(multipart_eos), ) ) - assert np.isfinite(enthalpy) - assert np.isfinite(sound_speed_si) and sound_speed_si > 0 - assert sound_speed_si / lal.C_SI < 1.1 - max_enthalpy = lalsim.SimNeutronStarEOSMultiPartsMaxPseudoEnthalpy( - loaded["twin_star"].eos + assert sound_si > 1.0 + assert sound_si / lal.C_SI < 1.1 + + legacy_eos = lalsim.SimNeutronStarEOSByName("SLY") + legacy = Adapter( + legacy_eos, reviewed_multibranch=False, lalsim_module=lalsim ) - assert np.isfinite(max_enthalpy) and max_enthalpy >= enthalpy - tidal_lambda = [ - loaded["twin_star"].lambda_from_m(mass, branch_id=branch) - for branch in (left, right) - ] - assert loaded["twin_star"].for_branch(left).test_speed_of_sound_causal() - assert all(value > 0 for value in radii + love + pressure + tidal_lambda) - assert not np.isclose(radii[0], radii[1], rtol=1e-10, atol=0) - assert not np.isclose(love[0], love[1], rtol=1e-10, atol=0) - assert not np.isclose(pressure[0], pressure[1], rtol=1e-10, atol=0) - assert not np.isclose(tidal_lambda[0], tidal_lambda[1], rtol=1e-10, atol=0) + legacy_mass = 0.5 * (legacy.minimum_mass() + legacy.maximum_mass()) + assert legacy.radius(legacy_mass, branch_id=0) > 0 + + extended = Adapter( + multipart_eos, + minimal=False, + reviewed_multibranch=True, + lalsim_module=lalsim, + ) + for name in ( + "SimNeutronStarFamBranchBaryonicMass", + "SimNeutronStarFamBranchLoveNumberK3", + "SimNeutronStarFamBranchLoveNumberK4", + ): + fn = getattr(lalsim, name, None) + if fn is None: + continue + assert np.isfinite(fn(mass, 0, extended.family)) + with pytest.raises(Exception): + fn(mass, 0, minimal.family) + + for _ in range(8): + eos_here = lalsim.SimNeutronStarEOSMultiPartsByName("SLY") + family_here = Adapter( + eos_here, reviewed_multibranch=True, lalsim_module=lalsim + ) + assert family_here.number_of_branches >= 1 + del family_here, eos_here + gc.collect() + + +def test_actual_reviewed_lalsimulation_tables(record_property): + import lalsimulation as lalsim + + manifest_path, manifest = _load_manifest() + _validate_reviewed_build(lalsim, manifest, record_property) + fixtures = manifest.get("fixtures") + if not fixtures: + pytest.skip("external reviewed EOS fixtures not supplied in manifest") + assert {"two_column", "nine_column"}.issubset(fixtures) + runner = Path(__file__).with_name("run_lalsim_eos_reviewed_fixture.py") + fixture_specs = [("two_column", 2), ("nine_column", 9)] + if "twin_star" in fixtures: + fixture_specs.append( + ("twin_star", int(fixtures["twin_star"]["columns"])) + ) + for name, expected_columns in fixture_specs: + assert expected_columns in (2, 9), ( + "file-loader fixtures must have 2 or 9 columns; four-column wiki " + "arrays require a separately provenance-recorded transform" + ) + fixture = fixtures[name] + path = (manifest_path.parent / fixture["path"]).resolve() + assert path.is_file(), "missing {} fixture: {}".format(name, path) + assert _sha256(path) == fixture["sha256"] + data = np.loadtxt(str(path), ndmin=2) + columns = data.shape[1] + assert columns == expected_columns + command = [ + sys.executable, + str(runner), + "--fixture", + str(path), + "--columns", + str(expected_columns), + ] + if name == "twin_star": + command.append("--twin") + _run_fixture_subprocess(command) + + nine_path = (manifest_path.parent / fixtures["nine_column"]["path"]).resolve() + for extra in ("--extended", "--eosmanager"): + _run_fixture_subprocess( + [ + sys.executable, + str(runner), + "--fixture", + str(nine_path), + "--columns", + "9", + extra, + ] + ) + + crash = manifest.get("known_upstream_crash") + if crash: + crash_path = (manifest_path.parent / crash["path"]).resolve() + assert _sha256(crash_path) == crash["sha256"] + crash_data = np.loadtxt(str(crash_path), ndmin=2) + assert crash_data.shape[1] == 4 + expected_codes = crash.get("expected_returncodes", [-11, 139]) + assert expected_codes == [-11, 139] + _run_expected_upstream_crash( + [ + sys.executable, + str(runner), + "--fixture", + str(crash_path), + "--columns", + "4", + "--arrays", + "--status", + os.devnull, + ] + ) diff --git a/docs/eos-interface-contract.md b/docs/eos-interface-contract.md index 50cbb1f6b..e413e8742 100644 --- a/docs/eos-interface-contract.md +++ b/docs/eos-interface-contract.md @@ -32,15 +32,21 @@ Reviewed two- or nine-column tables use: [--using-eos-extended-family] [--using-eos-branch ] ``` -On a reviewed build, RIFT dispatches these files through -`SimNeutronStarEOSFromFilePhaseTransition` and constructs their families with -`CreateSimNeutronStarFamilyPT`. Legacy named, spectral, piecewise-polytrope, -and ordinary file EOS objects continue to use the released one-argument family -constructor: both SWIG object families coexist and are selected explicitly. +By default, `lalsim_file:` preserves released `SimNeutronStarEOSFromFile` +behavior even on a reviewed build. Requesting a branch, the compatibility +dirty-phase-transition option, or the extended family opts into +`SimNeutronStarEOSFromFilePhaseTransition` and `CreateSimNeutronStarFamilyPT`. +Legacy named, spectral, piecewise-polytrope, and ordinary file EOS objects +continue to use the released one-argument family constructor: both SWIG object +families coexist and are selected explicitly. The reviewed PT reader always enables its dirty-phase-transition handling, so `--using-eos-dirty-phase-transitions` remains accepted for command-line compatibility but does not toggle a second reader mode. +The optional Python API `family_log_pressure_min` selects +`CreateSimNeutronStarFamilyPTWithPcmin`; its value is the finite natural +logarithm of central pressure in pascals, `ln(Pc / Pa)`. + For pseudo-pipe workflows, forward the flag with `--manual-extra-cip-args`. On O4d, Hydra hyperpipe configurations can put it in the post driver's `extra-args` when that driver is the fixed-EOS CIP @@ -89,33 +95,62 @@ this restriction can be relaxed. ## Reviewed-LALSimulation integration gate The fake-backed compatibility tests check RIFT's dispatch logic, but do not -certify the reviewed SWIG interface. To run the real-build gate, build the -exact reviewed LALSuite commit, activate that Python environment, and set +certify the reviewed SWIG interface. Fixture-free acceptance uses trusted +`SimNeutronStarEOSMultiPartsByName("SLY")`, checks coexistence with the legacy +named-EOS API, extended-family fields when exported, and repeated family +destruction. To run it, build the exact reviewed LALSuite commit, activate that +Python environment, and set `RIFT_REVIEWED_LALSIM_MANIFEST` to a JSON file with this shape: ```json { - "lalsuite_ref": "0123456789abcdef0123456789abcdef01234567", + "lalsuite_ref": "974c0ef468b76e8298e67fd8baf71ed259cc5fee", "fixtures": { "two_column": {"path": "two-column.dat", "sha256": "..."}, "nine_column": {"path": "nine-column.dat", "sha256": "..."}, - "twin_star": {"path": "twin-star.dat", "sha256": "...", "dirty_phase_transitions": true} + "twin_star": { + "path": "twin_star.dat", + "sha256": "1bc0fa5f92788cfb10cea3a7f04d5e6587a342043c94f1d6e7fbef2f7294ac16", + "columns": 2, + "raw_source_sha256": "3905405889b7da968e07a7ec97d37ca92c710a15af9407e014d6b30869b9608e", + "transform": "pressure=raw col4*1.602176634e32*G/c^4; energy=raw col3*1.602176634e32*G/c^4" + } + }, + "known_upstream_crash": { + "path": "upstream_2pt_raw.dat", + "sha256": "e8439b356e0815cfd710c3e4f0ccb7330a9ee7e42050ae6c5963c9d5890f8db7", + "expected_returncodes": [-11, 139] } } ``` Run `pytest MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_reviewed_integration.py`. -Paths are relative to the manifest. The ref must be the full 40-character +The `fixtures` object is optional, and `twin_star` is optional within it. +Without fixtures, builtin acceptance runs and external-table coverage is an +explicit skip. File-loader fixtures must be exact two- or nine-column numeric +tables; four-column wiki arrays require a separate provenance-recorded +transformation. Each table is loaded in a subprocess with a 60-second timeout, +stdout and stderr sent to `DEVNULL`, and only a status file capped at 64 KiB +returned to pytest. Paths are relative to the manifest. The ref must be the +full 40-character commit actually built by the job: the gate requires it to equal `lalsimulation.SimulationVCSInfo.vcsId` and requires a clean VCS build. Every -fixture hash is mandatory. Once the manifest enables the gate, missing modern +supplied fixture hash is mandatory. Once the manifest enables the gate, missing modern symbols, malformed or mismatched build provenance, missing fixtures, wrong column counts, or absence of distinct overlapping twin-star solutions are -failures. The gate exercises the always-PT multipart reader through both legacy -flag forms, minimal and extended family construction, and branch-indexed radius, -Love number, central pressure, and tidal deformability. Ordinary CI skips this -private-build gate explicitly. +failures. Ordinary CI skips this private-build gate explicitly. + +The recorded one-transition twin table is a two-column file derived using +`pressure = raw column 4 * 1.602176634e32 * G / c^4` and +`energy = raw column 3 * 1.602176634e32 * G / c^4`; both its derived and raw +SHA-256 values are retained above. The optional `known_upstream_crash` record +uses `SimNeutronStarEOSFromArraysPhaseTransition` with the same geometrized +factor. At reviewed commit `974c0ef`, the authoritative two-transition RC139 +array reaches `CreateSimNeutronStarFamilyPT` and exits with SIGSEGV (`-11`, or +shell `139`). The gate records that as an expected upstream `xfail`; if it +starts passing, the test fails with instructions to promote it to a required +passing fixture rather than silently accepting changed behavior. The drift-sentinel registry should eventually declare an EOS contract group with LALSuite and NuclearMatter-Backend as producers and RIFT/nmb-papers as From f7868535432dc0f0b12fb97aded5cb7b4805cb72 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Tue, 1 Sep 2026 21:11:19 -0400 Subject: [PATCH 15/16] tests: isolate known lalsim crash evidence --- .../test_lalsim_eos_reviewed_integration.py | 50 +++++++++++-------- 1 file changed, 30 insertions(+), 20 deletions(-) diff --git a/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_reviewed_integration.py b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_reviewed_integration.py index 9a0554752..efbc686a6 100644 --- a/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_reviewed_integration.py +++ b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_reviewed_integration.py @@ -271,24 +271,34 @@ def test_actual_reviewed_lalsimulation_tables(record_property): ] ) + + +def test_known_upstream_two_transition_crash(record_property): + """Keep the upstream SIGSEGV visible without masking passing fixtures.""" + import lalsimulation as lalsim + + manifest_path, manifest = _load_manifest() + _validate_reviewed_build(lalsim, manifest, record_property) crash = manifest.get("known_upstream_crash") - if crash: - crash_path = (manifest_path.parent / crash["path"]).resolve() - assert _sha256(crash_path) == crash["sha256"] - crash_data = np.loadtxt(str(crash_path), ndmin=2) - assert crash_data.shape[1] == 4 - expected_codes = crash.get("expected_returncodes", [-11, 139]) - assert expected_codes == [-11, 139] - _run_expected_upstream_crash( - [ - sys.executable, - str(runner), - "--fixture", - str(crash_path), - "--columns", - "4", - "--arrays", - "--status", - os.devnull, - ] - ) + if not crash: + pytest.skip("known upstream two-transition fixture not supplied") + crash_path = (manifest_path.parent / crash["path"]).resolve() + assert _sha256(crash_path) == crash["sha256"] + crash_data = np.loadtxt(str(crash_path), ndmin=2) + assert crash_data.shape[1] == 4 + expected_codes = crash.get("expected_returncodes", [-11, 139]) + assert expected_codes == [-11, 139] + runner = Path(__file__).with_name("run_lalsim_eos_reviewed_fixture.py") + _run_expected_upstream_crash( + [ + sys.executable, + str(runner), + "--fixture", + str(crash_path), + "--columns", + "4", + "--arrays", + "--status", + os.devnull, + ] + ) From 4f01edef7040675a3ba689b65a8687ecef7c3091 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Tue, 1 Sep 2026 21:13:20 -0400 Subject: [PATCH 16/16] EOS: align multipart helper keywords --- .../Code/RIFT/physics/EOSManager.py | 57 +++++++++++++++---- .../Code/test/test_lalsim_eos_compat.py | 33 ++++++++++- docs/eos-interface-contract.md | 7 +++ 3 files changed, 83 insertions(+), 14 deletions(-) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py b/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py index ba1442a5a..8192d0e08 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py @@ -13,8 +13,9 @@ rosDebug=False import numpy as np -import os -import sys +import os +import sys +import warnings import lal import lalsimulation as lalsim from scipy.integrate import quad @@ -1254,21 +1255,47 @@ def int_func(x_prime): ### Utilities ### -# Les-like +# Les-like +def _resolve_multipart_helper_keyword(multipart, reviewed_multibranch): + """Resolve the canonical helper keyword and its deprecated O4c alias.""" + if reviewed_multibranch is not None: + warnings.warn( + "reviewed_multibranch is deprecated; use multipart", + DeprecationWarning, + stacklevel=3, + ) + if ( + multipart is not None + and bool(multipart) != bool(reviewed_multibranch) + ): + raise ValueError( + "conflicting multipart and reviewed_multibranch values" + ) + multipart = reviewed_multibranch + return False if multipart is None else bool(multipart) + + def make_mr_lambda_lal(eos, n_bins=100, branch_id=None, - reviewed_multibranch=False, family_adapter=None): + multipart=None, family_adapter=None, + reviewed_multibranch=None): ''' Construct mass-radius curve from EOS Based on modern code resources (https://git.ligo.org/publications/gw170817/bns-eos/blob/master/scripts/eos-params.py) which access low-level structures ``branch_id`` is optional for released/single-branch LALSimulation. It is required for a multibranch family so an overlapping twin-star interval is - never collapsed silently. Set ``reviewed_multibranch`` only for an EOS - returned by ``SimNeutronStarEOSFromFilePhaseTransition``. - ''' + never collapsed silently. Set ``multipart=True`` only for an EOS returned + by ``SimNeutronStarEOSFromFilePhaseTransition``. + + ``reviewed_multibranch`` is a deprecated O4c keyword alias for + ``multipart``. Supplying both with different values is an error. + ''' + multipart = _resolve_multipart_helper_keyword( + multipart, reviewed_multibranch + ) family = family_adapter or create_family( eos, lalsim_module=lalsim, - reviewed_multibranch=reviewed_multibranch, + reviewed_multibranch=multipart, ) if family.number_of_branches > 1 and branch_id is None: raise ValueError( @@ -1292,12 +1319,18 @@ def make_mr_lambda_lal(eos, n_bins=100, branch_id=None, def make_mr_lambda_lal_branches(eos, n_bins=100, - reviewed_multibranch=False, - family_adapter=None): - """Return branch curves; multipart construction is explicit and opt-in.""" + multipart=None, family_adapter=None, + reviewed_multibranch=None): + """Return branch curves with canonical ``multipart`` dispatch. + + ``reviewed_multibranch`` remains as a deprecated O4c keyword alias. + """ + multipart = _resolve_multipart_helper_keyword( + multipart, reviewed_multibranch + ) family = family_adapter or create_family( eos, lalsim_module=lalsim, - reviewed_multibranch=reviewed_multibranch, + reviewed_multibranch=multipart, ) return { branch_id: _make_mr_lambda_for_family(family, n_bins, branch_id) diff --git a/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py index 41c822c18..e949b6560 100644 --- a/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py +++ b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py @@ -310,16 +310,45 @@ def test_selected_branch_view_preserves_legacy_scalar_consumer_api(monkeypatch): assert primary.lambda_from_m(2.5) == pytest.approx(1e-8) curves = EOSManager.make_mr_lambda_lal_branches( - eos.eos, n_bins=3, family_adapter=eos._get_lalsim_family_adapter() + eos.eos, n_bins=3, multipart=True ) assert set(curves) == {0, 1} assert curves[0].shape == (3, 3) explicit = EOSManager.make_mr_lambda_lal( - eos.eos, n_bins=3, branch_id=1, reviewed_multibranch=True + eos.eos, n_bins=3, branch_id=1, multipart=True ) assert explicit.shape == (3, 3) + with pytest.warns(DeprecationWarning, match="use multipart"): + alias = EOSManager.make_mr_lambda_lal( + eos.eos, + n_bins=3, + branch_id=1, + reviewed_multibranch=True, + ) + assert alias.shape == (3, 3) + + with pytest.warns(DeprecationWarning, match="use multipart"): + alias_branches = EOSManager.make_mr_lambda_lal_branches( + eos.eos, n_bins=3, reviewed_multibranch=True + ) + assert set(alias_branches) == {0, 1} + + for helper, kwargs in ( + (EOSManager.make_mr_lambda_lal, {"branch_id": 1}), + (EOSManager.make_mr_lambda_lal_branches, {}), + ): + with pytest.warns(DeprecationWarning, match="use multipart"): + with pytest.raises(ValueError, match="conflicting multipart"): + helper( + eos.eos, + n_bins=3, + multipart=True, + reviewed_multibranch=False, + **kwargs + ) + def test_mr_lambda_branches_helper_keeps_legacy_default(monkeypatch): from RIFT.physics import EOSManager diff --git a/docs/eos-interface-contract.md b/docs/eos-interface-contract.md index e413e8742..47ac0b854 100644 --- a/docs/eos-interface-contract.md +++ b/docs/eos-interface-contract.md @@ -47,6 +47,13 @@ The optional Python API `family_log_pressure_min` selects `CreateSimNeutronStarFamilyPTWithPcmin`; its value is the finite natural logarithm of central pressure in pascals, `ln(Pc / Pa)`. +The mass-radius-Lambda helpers use the same canonical keyword as O4d: +`make_mr_lambda_lal(..., multipart=True)` and +`make_mr_lambda_lal_branches(..., multipart=True)`. The older O4c keyword +`reviewed_multibranch` remains a deprecated alias. Supplying both keywords with +different values raises an error rather than selecting an object family +ambiguously. + For pseudo-pipe workflows, forward the flag with `--manual-extra-cip-args`. On O4d, Hydra hyperpipe configurations can put it in the post driver's `extra-args` when that driver is the fixed-EOS CIP