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/CHANGES.rst b/CHANGES.rst index c824160a7..6953d1327 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,3 +1,15 @@ +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. 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. + 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..8192d0e08 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py @@ -13,14 +13,20 @@ 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 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 +78,66 @@ 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, 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 + 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, + reviewed_multibranch=getattr( + self, "_lalsim_reviewed_multibranch", False + ), + ) + 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 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: + 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,7 +145,19 @@ def lambda_from_m(self, m): return dimensionless_lam - def estimate_baryon_mass_from_mg(self,m): + 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, 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) @@ -97,7 +165,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=branch_id) / 1e3 return m + (1./r1p4)*m**2 def pressure_density_on_grid_alternate(self,logrho_grid,enforce_causal=False): @@ -168,7 +236,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 @@ -183,20 +252,43 @@ 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 - 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): + max_enthalpy = getattr( + lalsim, + "SimNeutronStarEOSMultiPartsMaxPseudoEnthalpy", + None, + ) + if max_enthalpy is None: + return False + hmax = max_enthalpy(eos) + else: + hmax = lalsim.SimNeutronStarEOSMaxPseudoEnthalpy(eos) else: try: - pmax = lalsim.SimNeutronStarCentralPressure(m_max_SI,fam) - hmax = lalsim.SimNeutronStarEOSPseudoEnthalpyOfPressure(pmax,eos) + pmax = family.central_pressure(m_max_SI, branch_id=branch_id) + 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 @@ -208,11 +300,86 @@ def test_speed_of_sound_causal(self, test_only_under_mmax=True,fast_test=True): # 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 + +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 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) + ### ### SERVICE 1: lalsimutils structure ### @@ -222,11 +389,62 @@ 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 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, + phase_transition_aware=False): + self.name = name or os.path.basename(fname) + self.fname = fname + self.eos = None + self.eos_fam = None + phase_transition_reader = getattr( + lalsim, "SimNeutronStarEOSFromFilePhaseTransition", None + ) + 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 + ) + 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: + self._set_lalsim_family( + minimal=minimal_family, + reviewed_multibranch=self._lalsim_reviewed_multibranch, + log_pressure_min=family_log_pressure_min, + ) + else: + self.mMaxMsun = None + + ### ### SERVICE 2: EOSFromFile @@ -325,7 +543,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 +598,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, lalsim_module=lalsim) + 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 +706,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 +741,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 +894,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 @@ -1043,28 +1255,104 @@ def int_func(x_prime): ### Utilities ### -# Les-like -def make_mr_lambda_lal(eos,n_bins=100): +# 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, + 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 - ''' - fam=lalsim.CreateSimNeutronStarFamily(eos) - max_m = lalsim.SimNeutronStarMaximumMass(fam)/lal.MSUN_SI - min_m = lalsim.SimNeutronStarFamMinimumMass(fam)/lal.MSUN_SI + + ``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 ``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=multipart, + ) + 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, + 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=multipart, + ) + 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 +1362,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, lalsim_module=lalsim) + fam = family.family r_cut = 40 # Some EOS we consider for PE purposes will have very large radius! @@ -1108,8 +1397,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:" @@ -1496,6 +1786,280 @@ 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_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, + 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 = {} + 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")]) + 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 + + +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 = {} + 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 + 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 + + +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. + + 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 + rep = schema = "" + try: + with h5py.File(fname, 'r') as f: + rep = str(f.attrs.get("representation", "")) + schema = str(f.attrs.get("schema_version", "")) + except Exception: + 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) + + #### #### General lalsimulation interfacing #### 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..45c446aba --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/physics/lalsim_eos_compat.py @@ -0,0 +1,244 @@ +"""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. +""" + +import math + + +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. + + 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. + 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 natural-log central-pressure bound, ``ln(Pc / Pa)``, for + the reviewed + ``CreateSimNeutronStarFamilyPTWithPcmin`` constructor. + lalsim_module: + Dependency-injection hook used by the interface contract tests. + """ + + _REVIEWED_REQUIRED = ( + "CreateSimNeutronStarFamilyPT", + "SimNeutronStarFamNumberOfBranches", + "SimNeutronStarFamBranchMinMass", + "SimNeutronStarFamBranchMaxMass", + "SimNeutronStarFamBranchRadius", + "SimNeutronStarFamBranchLoveNumberK2", + "SimNeutronStarFamBranchCentralPressure", + ) + + @classmethod + 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._REVIEWED_REQUIRED, present) + if not available + ] + raise RuntimeError( + "reviewed LALSimulation phase-transition API is incomplete; " + "missing symbols: {}" + .format(", ".join(missing)) + ) + return True + + 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 = 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. + if log_pressure_min is None: + self.family = self.lalsim.CreateSimNeutronStarFamilyPT( + 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 + ) + if constructor is None: + raise NotImplementedError( + "reviewed LALSimulation build does not expose " + "CreateSimNeutronStarFamilyPTWithPcmin" + ) + self.family = constructor( + eos, int(bool(minimal)), log_pressure_min + ) + else: + self.family = self.lalsim.CreateSimNeutronStarFamily(eos) + + @classmethod + 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 + obj = cls.__new__(cls) + obj.lalsim = lalsim_module + obj.eos = None + obj.family = family + obj.is_multibranch_api = bool(reviewed_multibranch) + if obj.is_multibranch_api: + cls._require_reviewed_api(obj.lalsim) + 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: + return min(self.minimum_mass(k) for k in range(self.number_of_branches)) + self._validate_branch_id(branch_id) + return self.lalsim.SimNeutronStarFamBranchMinMass( + int(branch_id), self.family + ) + 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: + return max(self.maximum_mass(k) for k in range(self.number_of_branches)) + self._validate_branch_id(branch_id) + return self.lalsim.SimNeutronStarFamBranchMaxMass( + int(branch_id), self.family + ) + 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.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.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) + if self.is_multibranch_api: + return self.lalsim.SimNeutronStarFamBranchCentralPressure( + mass_si, resolved, self.family + ) + 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, + reviewed_multibranch=False, log_pressure_min=None): + """Return a :class:`LALSimNeutronStarFamilyAdapter` for ``eos``.""" + return LALSimNeutronStarFamilyAdapter( + 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 941294d49..5f56dde0e 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py @@ -335,6 +335,9 @@ 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-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.") parser.add_argument("--no-matter1", action='store_true', help="Set the lambda parameters to zero (BBH) but return them") @@ -416,6 +419,24 @@ 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) +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( + "--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 @@ -442,7 +463,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] @@ -451,13 +472,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 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] 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)) @@ -500,6 +521,20 @@ 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 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, + phase_transition_aware=opts.using_eos_branch is not None, + ) elif 'lal_' in eos_name: eos_name = eos_name.replace('lal_','') my_eos = EOSManager.EOSLALSimulation(name=eos_name) @@ -515,6 +550,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 @@ -1924,7 +1969,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) @@ -3617,5 +3663,3 @@ def parse_corr_params(my_str): print(" Failed to generate corner for ", extra_plot_coord_names[indx]) sys.exit(0) - - 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 new file mode 100644 index 000000000..e949b6560 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_compat.py @@ -0,0 +1,484 @@ +import json + +import lal +import numpy as np +import pytest + +from RIFT.physics.lalsim_eos_compat import ( + AmbiguousFamilyBranchError, + LALSimNeutronStarFamilyAdapter, + validate_fixed_eos_branch_request, +) + + +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 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 SimNeutronStarEOSFromFilePhaseTransition(self, fname): + self.file_calls.append((fname,)) + return "multipart-eos" + + 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 SimNeutronStarFamBranchMinMass(self, branch_id, family): + return self.bounds[branch_id][0] + + def SimNeutronStarFamBranchMaxMass(self, branch_id, family): + return self.bounds[branch_id][1] + + def SimNeutronStarFamBranchRadius(self, mass, branch_id, family): + return 10.0 * branch_id + mass + + def SimNeutronStarFamBranchLoveNumberK2(self, mass, branch_id, family): + return branch_id + 0.1 * mass + + def SimNeutronStarFamBranchCentralPressure(self, mass, branch_id, family): + return 100.0 * branch_id + mass + + +class StellarMassMultibranchLALSimulation(MultibranchLALSimulation): + bounds = tuple( + (lower * lal.MSUN_SI, upper * lal.MSUN_SI) + for lower, upper in ((1.0, 2.0), (1.3, 3.0)) + ) + + def SimNeutronStarEOSMultiPartsPseudoEnthalpyOfPressure(self, pressure, eos): + return pressure + + 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 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" + + 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(): + 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_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, + reviewed_multibranch=True, + ) + + 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 + + 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() + 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, + "SimNeutronStarFamBranchLoveNumberK2", + ) + 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(), + reviewed_multibranch=True, + ) + + 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",)] + 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",) + 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) + + 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 + + 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_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", phase_transition_aware=True + ) + + 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) + + curves = EOSManager.make_mr_lambda_lal_branches( + 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, 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 + + 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", 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( + 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): + 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", 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) + + monkeypatch.delattr( + StellarMassMultibranchLALSimulation, + "SimNeutronStarEOSMultiPartsSpeedOfSoundOfPseudoEnthalpy", + ) + 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) + + 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 + + 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) + + +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/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..efbc686a6 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_lalsim_eos_reviewed_integration.py @@ -0,0 +1,304 @@ +"""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 importlib.util +import json +import os +from pathlib import Path +import re +import subprocess +import sys +import tempfile + +import numpy as np +import pytest + + +MANIFEST_ENV = "RIFT_REVIEWED_LALSIM_MANIFEST" +REQUIRED_SYMBOLS = ( + "SimulationVCSInfo", + "SimNeutronStarEOSMultiPartsByName", + "SimNeutronStarEOSFromArraysPhaseTransition", + "SimNeutronStarEOSFromFilePhaseTransition", + "CreateSimNeutronStarFamilyPT", + "SimNeutronStarFamNumberOfBranches", + "SimNeutronStarFamBranchMinMass", + "SimNeutronStarFamBranchMaxMass", + "SimNeutronStarFamBranchRadius", + "SimNeutronStarFamBranchLoveNumberK2", + "SimNeutronStarFamBranchCentralPressure", + "SimNeutronStarEOSMultiPartsMaxPseudoEnthalpy", + "SimNeutronStarEOSMultiPartsPseudoEnthalpyOfPressure", + "SimNeutronStarEOSMultiPartsSpeedOfSoundOfPseudoEnthalpy", +) + + +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 _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 + + +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"] + 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 + ) + ) + + +def _run_expected_upstream_crash(command): + result = subprocess.run( + command, + stdout=subprocess.DEVNULL, + stderr=subprocess.DEVNULL, + timeout=60, + check=False, + ) + if result.returncode == 0: + pytest.fail( + "known two-transition native crash is fixed; promote this fixture " + "to a required passing gate" + ) + assert result.returncode in (-11, 139), ( + "expected SIGSEGV/139, got {}".format(result.returncode) + ) + pytest.xfail("upstream reviewed LALSimulation two-transition SIGSEGV") + + +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)" + ) + 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 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, multipart_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 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 + ) + 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, + ] + ) + + + +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 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, + ] + ) diff --git a/docs/eos-interface-contract.md b/docs/eos-interface-contract.md new file mode 100644 index 000000000..47ac0b854 --- /dev/null +++ b/docs/eos-interface-contract.md @@ -0,0 +1,165 @@ +# 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 +``` + +Reviewed two- or nine-column tables use: + +```text +--using-eos lalsim_file: [--using-eos-dirty-phase-transitions] + [--using-eos-extended-family] [--using-eos-branch ] +``` + +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)`. + +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 +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. + +`--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. 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": "974c0ef468b76e8298e67fd8baf71ed259cc5fee", + "fixtures": { + "two_column": {"path": "two-column.dat", "sha256": "..."}, + "nine_column": {"path": "nine-column.dat", "sha256": "..."}, + "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`. +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 +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. 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 +consumers. The current registry only covers RIFT/hyperpipe operational archive +and queue boundaries, so it cannot detect EOS schema or callable drift yet.