EOSManager: support reviewed LALSimulation families (O4c)

.gitattributes

@@ -0,0 +1 @@
+MonteCarloMarginalizeCode/Code/RIFT/physics/EOSManager.py whitespace=cr-at-eol

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:<path>`` 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

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,
+    )

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:<path>. 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:<path>, 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:<path>"
+    )
 
 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:<sequence_file.h5>:<index>
+        # 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)
-
-