Moving em_progenitors.py and ns_sequence data from tmpltbank to new neutron_stars package

bin/pycbc_create_injections

@@ -206,7 +206,7 @@ dists = distributions.read_distributions_from_config(cp, opts.dist_section)
 
 # construct class that will draw the samples
 randomsampler = JointDistribution(variable_params, *dists,
-                               **{"constraints" : constraints})
+                               **{"constraints": constraints})
 
 if opts.ninjections:
     draw_size = opts.ninjections

pycbc/conversions.py

@@ -23,7 +23,7 @@
 # =============================================================================
 #
 """
-This modules provides a library of functions that calculate waveform parameters
+This module provides a library of functions that calculate waveform parameters
 from other parameters. All exposed functions in this module's namespace return
 one parameter given a set of inputs.
 """
@@ -33,7 +33,10 @@
 import lal
 from pycbc.detector import Detector
 import pycbc.cosmology
-from .coordinates import spherical_to_cartesian as _spherical_to_cartesian
+from .coordinates import (
+    spherical_to_cartesian as _spherical_to_cartesian,
+    cartesian_to_spherical as _cartesian_to_spherical)
+from pycbc import neutron_stars as ns
 
 pykerr = pycbc.libutils.import_optional('pykerr')
 lalsim = pycbc.libutils.import_optional('lalsimulation')
@@ -417,6 +420,92 @@ def lambda_from_mass_tov_file(mass, tov_file, distance=0.):
     lambdav = numpy.interp(mass_src, mass_from_file, lambda_from_file)
     return lambdav
 
+
+def remnant_mass_from_mass1_mass2_spherical_spin_eos(
+        mass1, mass2, spin1a=0.0, spin1pol=0.0, eos='2H'):
+    """
+    Function that determines the remnant disk mass of an NS-BH system
+    using the fit to numerical-relativity results discussed in
+    Foucart, Hinderer & Nissanke, PRD 98, 081501(R) (2018).
+    The BH spin may be misaligned with the orbital angular momentum.
+    In such cases the ISSO is approximated following the approach of
+    Stone, Loeb & Berger, PRD 87, 084053 (2013), which was originally
+    devised for a previous NS-BH remnant mass fit of
+    Foucart, PRD 86, 124007 (2012).
+    Note: NS spin is assumed to be 0!
+
+    Parameters
+    -----------
+    mass1 : float
+        The mass of the black hole, in solar masses.
+    mass2 : float
+        The mass of the neutron star, in solar masses.
+    spin1a : float, optional
+        The dimensionless magnitude of the spin of mass1. Default = 0.
+    spin1pol : float, optional
+        The tilt angle of the spin of mass1. Default = 0 (aligned w L).
+    eos : str, optional
+        Name of the equation of state being adopted. Default is '2H'.
+
+    Returns
+    ----------
+    remnant_mass: float
+        The remnant mass in solar masses
+    """
+    mass1, mass2, spin1a, spin1pol, input_is_array = ensurearray(
+        mass1, mass2, spin1a, spin1pol)
+    # mass1 must be greater than mass2
+    try:
+        if any(mass2 > mass1) and input_is_array:
+            raise ValueError(f'Require mass1 >= mass2')
+    except TypeError:
+        if mass2 > mass1 and not input_is_array:
+            raise ValueError(f'Require mass1 >= mass2. {mass1} < {mass2}')
+    ns_compactness, ns_b_mass = ns.initialize_eos(mass2, eos)
+    eta = eta_from_mass1_mass2(mass1, mass2)
+    remnant_mass = ns.foucart18(
+        eta, ns_compactness, ns_b_mass, spin1a, spin1pol)
+    return formatreturn(remnant_mass, input_is_array)
+
+
+def remnant_mass_from_mass1_mass2_cartesian_spin_eos(
+        mass1, mass2, spin1x=0.0, spin1y=0.0, spin1z=0.0, eos='2H'):
+    """
+    Function that determines the remnant disk mass of an NS-BH system
+    using the fit to numerical-relativity results discussed in
+    Foucart, Hinderer & Nissanke, PRD 98, 081501(R) (2018).
+    The BH spin may be misaligned with the orbital angular momentum.
+    In such cases the ISSO is approximated following the approach of
+    Stone, Loeb & Berger, PRD 87, 084053 (2013), which was originally
+    devised for a previous NS-BH remnant mass fit of
+    Foucart, PRD 86, 124007 (2012).
+    Note: NS spin is assumed to be 0!
+
+    Parameters
+    -----------
+    mass1 : float
+        The mass of the black hole, in solar masses.
+    mass2 : float
+        The mass of the neutron star, in solar masses.
+    spin1x : float, optional
+        The dimensionless x-component of the spin of mass1. Default = 0.
+    spin1y : float, optional
+        The dimensionless y-component of the spin of mass1. Default = 0.
+    spin1z : float, optional
+        The dimensionless z-component of the spin of mass1. Default = 0.
+    eos: str, optional
+        Name of the equation of state being adopted. Default is '2H'.
+
+    Returns
+    ----------
+    remnant_mass: float
+        The remnant mass in solar masses
+    """
+    spin1a, _, spin1pol = _cartesian_to_spherical(spin1x, spin1y, spin1z)
+    return remnant_mass_from_mass1_mass2_spherical_spin_eos(
+        mass1, mass2, spin1a, spin1pol, eos=eos)
+
+
 #
 # =============================================================================
 #
@@ -1537,5 +1626,7 @@ def nltides_gw_phase_diff_isco(f_low, f0, amplitude, n, m1, m2):
            'optimal_dec_from_detector', 'optimal_ra_from_detector',
            'chi_eff_from_spherical', 'chi_p_from_spherical',
            'nltides_gw_phase_diff_isco', 'spin_from_pulsar_freq',
-           'freqlmn_from_other_lmn', 'taulmn_from_other_lmn'
+           'freqlmn_from_other_lmn', 'taulmn_from_other_lmn',
+           'remnant_mass_from_mass1_mass2_spherical_spin_eos',
+           'remnant_mass_from_mass1_mass2_cartesian_spin_eos'
           ]

pycbc/distributions/constraints.py

@@ -16,6 +16,7 @@
 This modules provides classes for evaluating multi-dimensional constraints.
 """
 
+import re
 import scipy.spatial
 import numpy
 import h5py
@@ -37,7 +38,15 @@ def __init__(self, constraint_arg, static_args=None, transforms=None,
                 static_args.items(), key=lambda x: len(x[0]), reverse=True))
             )
         for arg, val in static_args.items():
-            constraint_arg = constraint_arg.replace(arg, str(val))
+            swp = f"'{val}'" if isinstance(val, str) else str(val)
+            # Substitute static arg name for value if it appears in the
+            # constraint_arg string at the beginning of a word and is not
+            # followed by an underscore or equals sign.
+            # This ensures that static_args that are also kwargs in function calls are
+            # handled correctly, i.e., the kwarg is not touched while its value is replaced
+            # with the static_arg value.
+            constraint_arg = re.sub(
+                r'\b{}(?!\_|\=)'.format(arg), swp, constraint_arg)
         self.constraint_arg = constraint_arg
         self.transforms = transforms
         for kwarg in kwargs.keys():

pycbc/io/record.py

@@ -1802,7 +1802,8 @@ class WaveformArray(_FieldArrayWithDefaults):
         parameters.spin_px, parameters.spin_py, parameters.spin_pz,
         parameters.spin_sx, parameters.spin_sy, parameters.spin_sz,
         parameters.spin1_a, parameters.spin1_azimuthal, parameters.spin1_polar,
-        parameters.spin2_a, parameters.spin2_azimuthal, parameters.spin2_polar]
+        parameters.spin2_a, parameters.spin2_azimuthal, parameters.spin2_polar,
+        parameters.remnant_mass]
 
     @property
     def primary_mass(self):
@@ -1912,5 +1913,12 @@ def spin2_polar(self):
         return coordinates.cartesian_to_spherical_polar(
                                      self.spin2x, self.spin2y, self.spin2z)
 
+    @property
+    def remnant_mass(self):
+        """Returns the remnant mass for an NS-BH binary."""
+        return conversions.remnant_mass_from_mass1_mass2_cartesian_spin_eos(
+                                     self.mass1, self.mass2, self.spin1x,
+                                     self.spin1y, self.spin1z)
+
 
 __all__ = ['FieldArray', 'WaveformArray']

pycbc/neutron_stars/__init__.py

@@ -0,0 +1,9 @@
+import os.path
+# Setup the directory with the NS equilibrium sequence(s)
+NS_DATA_DIRECTORY = os.path.join(
+    os.path.dirname(__file__), 'ns_data')
+NS_SEQUENCES = [
+    f.replace('equil_', '').replace('.dat', '')
+    for f in os.listdir(NS_DATA_DIRECTORY) if f.endswith('.dat')]
+from pycbc.neutron_stars.eos_utils import *
+from pycbc.neutron_stars.pg_isso_solver import *