Moving em_progenitors.py and ns_sequence data from tmpltbank to new n…

…eutron_stars package (#3937)

* Remnant mass function moved from em_progenitors to conversions. Relevant NS functions moved to neutron_stars. Calls to these functions fixed to reflect new locations

* no change, pushing to resolve issue with checkout

* importing pycbc.tmpltbank rather than pycbc in attempt to resolve import issues (currently unresolved)

* moved neutron stars code and eos data to new package and created __init__ file

* Renamed em_progenitors as ns_functions and moved to neutron_stars package. Edited how said functions are imported elsewhere where necessary. Added neutron_stars package path to setup file

* First suite of changes suggested in PR. Removes some cases of double imports and old code comments from during development

* Removed the logging and sys calls. Changed the exception to a ValueError with clear error messages.

Co-authored-by: Andrew R. Williamson <[email protected]>

* retrigger checks

* EM-bright work

* Constraint string substitutions for static args

* Speed up remnant mass calculation (1st attempt)

* Adding derivatives in ns_functions for root solver

* Add to comment

* Missing bracket

* Fix bugs

* Enable spherical or cartesian spins for m_rem calc

* Reorganising neutron_stars module

* Some fixes to isso solver

* Adding 2D interpolation for ISSO radii

* Refined bivariate interpolator for ISSO calc

* Change in setup.py, needed for tox

* Fix docstrings

* Fixing doc strings

* Remove old pycbc_dark_vs_bright_injections code

* Added newline at end of file to conform to codeclimate check

* removed trailing whitespace

* removed trailing whitespaces

* Removed imports outside of toplevel

* removed trailing whitespaces

* Empty-Commit, retrigger checks

* correctly importing pickle and os.path at top of pg_isso_solver.py

* Creating 3 variables in which previously repeated lines of calculation are performed (in need of new names, potenatially better solution available)

* fixed issue with unnecessary indents

* changed variables to functions after negecting to realise they relied on called variables. Still in need of more informative names, potentially room for more elegant solution

* removed trailing whitespaces

* Fixing error in trig functions (missing return)

* removed functions aimed at avoiding repeated code. Deemed unnecessary

* rectified mislabeled parameter chi_lims to bounds as required by function

* introducing latest version of remnant mass function in coord_utils get_randm_mass function. Minor rewrite to earlier portion sof function to acommodate this change

* added short description of concat_grid function

* Moved ns_sequence_file initialisation above it's first use to avoid error

* fixed issue with methods used to call variables in remnant mass function

* introduced eos as a variable in the get_random_mass function in order to work with updated remnant_mass funtion

* adding new lines between function definitions

* re-introduced eta_nsbh in get_random_mass as it is required later in the function

* moved module level imports to top of init file

* moved module level imports back to bottom of init file as moving to top seemed to break code

* removed trailing whitespace

* Initialised input_is_array variable and added functionality to check if input variable is an array. If it is, input_is_array is set to True

* resolved issue with undefined variable erroneously added in code which checks if an input is an array

* addrd r prefix to doctring in an attempt to resolve issues related to use of backslashes

* Replace code which checks if input is an array with more reliable method

Co-authored-by: Tito Dal Canton <[email protected]>

* Updated docstrings to more accurately represent the forms of inputs and outputs in a number of functions

* Fixed typo in docstring

Co-authored-by: Tito Dal Canton <[email protected]>

* Reformatted code to agree with codeclimate

Co-authored-by: Tito Dal Canton <[email protected]>

* Removed hardcoded output path for generate_isso_bivariate_interp()

* Fix wrong import name

* Updated comments in constraint function

Co-authored-by: Francesco Pannarale <[email protected]>

* Remove ISSO interpolant, will bring back a better version later

* Forgot a few things

* Update pycbc/tmpltbank/coord_utils.py

Co-authored-by: Tito Dal Canton <[email protected]>

Co-authored-by: Samuel Higginbotham <[email protected]>
Co-authored-by: Andrew R. Williamson <[email protected]>
Co-authored-by: Andrew Williamson <[email protected]>
Co-authored-by: Tito Dal Canton <[email protected]>
Co-authored-by: Francesco Pannarale <[email protected]>
Co-authored-by: Tito Dal Canton <[email protected]>

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 *