# Copyright (C) 2023, 2024, 2025 Cecilio García Quirós
"""
Define generators to interface with gwsignal.
"""
try:
from gwtemplate.waveform import CompactBinaryCoalescenceGenerator
import astropy.units as u
from gwpy.timeseries import TimeSeries
from gwpy.frequencyseries import FrequencySeries
except:
raise ImportWarning("Cannot use gwtemplate interface for phenomxpy.")
import phenomxpy
from phenomxpy.utils import parse_options_from_approximant_name, convert_params, MasstoSecond, SecondtoMass
import numpy as np
default_dict = {
"mass1": 1.0 * u.solMass,
"mass2": 1.0 * u.solMass,
"spin1x": 0.0 * u.dimensionless_unscaled,
"spin1y": 0.0 * u.dimensionless_unscaled,
"spin1z": 0.0 * u.dimensionless_unscaled,
"spin2x": 0.0 * u.dimensionless_unscaled,
"spin2y": 0.0 * u.dimensionless_unscaled,
"spin2z": 0.0 * u.dimensionless_unscaled,
"deltaT": 1.0 / 64.0 * u.s,
"f22_start": 20.0 * u.Hz,
"f22_ref": 20.0 * u.Hz,
"phi_ref": 0.0 * u.rad,
"distance": 1.0 * u.Mpc,
"inclination": 0.0 * u.rad,
"eccentricity": 0.0 * u.dimensionless_unscaled,
"longAscNodes": 0.0 * u.rad,
"meanPerAno": 0.0 * u.rad,
}
[docs]
def strip_units(waveform_dict):
"""
Remove units from astropy dictionary.
"""
new_dc = {}
for key in waveform_dict.keys():
new_dc[key] = waveform_dict[key].value if isinstance(waveform_dict[key], u.Quantity) else waveform_dict[key]
return new_dc
[docs]
def convert_params_from_gwsignal(input_params):
"""
Convert gwsignal (astropy dictionary) to phenomxpy.
Parameters are transformed to the units in the default_dict of gwsignal.
Units are removed.
Parameters
----------
input_params: gwsiganl dictionary
parameters with astropy units.
Returns
-------
dict
Parameters in phenomxpy format without units.
"""
# Deep copy to not modify the original dictionary
params = input_params.copy()
# Conversion to units used in phenomxpy. It uses the same units as defined in gwsignal.core.parameter_conventions.default_dict
for key in params:
if key in default_dict and hasattr(params[key], "unit") and params[key].unit is not u.dimensionless_unscaled:
params[key] = params[key].to(default_dict[key].unit)
# Renaming some keys
for old_key, new_key in {
"f22_start": "f_min",
"f22_ref": "f_ref",
"deltaT": "delta_t",
"deltaF": "delta_f",
"ModeArray": "mode_array",
}.items():
if old_key in params:
params[new_key] = params.pop(old_key)
return convert_params(strip_units(params))
[docs]
def to_numpy(x):
"""
Convert to numpy array.
"""
return x.get() if hasattr(x, "get") else x
[docs]
class PyIMRPhenomT(CompactBinaryCoalescenceGenerator):
"""
gwsignal generator wrapper for the IMRPhenomT model.
This is the parent class from which all the other IMRPhenomT* models subclass.
Example usage:
.. code-block:: python
from lalsimulation import gwsignal as gws
gen = PyIMRPhenomT()
hp, hc = gws.core.waveform.GenerateTDWaveform(params, gen)
"""
def __init__(self, **kwargs):
super().__init__(**kwargs)
self._domain = "time"
self._implemented_domain = "time"
self._generation_domain = None
self._approximant = "IMRPhenomT"
# This imports the corresponding Phenom class but doesn't create the instance.
# Instances are created inside the methods `generate_td_waveform`, etc.
self._class = getattr(phenomxpy, self._approximant)
@property
def metadata(self):
return self._class.metadata()
@property
def epoch(self):
"""
Start time of the currently initialized waveform model in seconds.
The epoch is the time of the peak of the 22 amplitude.
"""
if not hasattr(self, "model"):
raise RuntimeError("'epoch' is not available. Call initialize_model(**parameters) first.")
return self.model.epoch * u.s
@property
def duration(self):
"""
Duration of the currently initialized waveform model in seconds.
The approximate end of the waveform is 500M after the peak of the 22 amplitude.
If the model is initialized with a deltaT, the duration includes the last sample.
If the model is initialized without deltaT, the duration still reports 500M after merger but it can be evaluated beyond that time.
"""
if not hasattr(self, "model"):
raise RuntimeError("'duration' is not available. Call initialize_model(**parameters) first.")
if hasattr(self.model.pWF, "length"):
return self.model.pWF.length * MasstoSecond(self.model.pWF.delta_t, self.model.pWF.total_mass) * u.s
return (MasstoSecond(self.model.pWF.tEnd, self.model.pWF.total_mass) - self.model.epoch) * u.s
[docs]
def initialize_model(self, times=None, **parameters):
"""
Initialize the waveform model with the given parameters.
Parameters
----------
times: xp.ndarray
Time array in seconds. Only used here to check for consistent input parameters.
**parameters: dict
Parameters for the waveform model.
Returns
-------
bool
True if the model is initialized successfully.
"""
if times is not None and parameters:
raise ValueError("Cannot pass both 'times' and 'parameters' at the same time.")
if parameters:
parameters = self.parameter_check(**parameters)
# Convert from gwsignal dictionary to phenomxpy dictionary
self.waveform_dict = convert_params_from_gwsignal(parameters)
self.model = self._class(**self.waveform_dict)
return True
elif times is not None:
if hasattr(self, "model") is False:
raise RuntimeError("Model not initialized. Call initialize_model(**parameters) first.")
elif hasattr(self, "model") is False:
raise RuntimeError("Model not initialized. Call initialize_model(**parameters) first.")
return True
def _generate_td_polarizations(self, **parameters):
"""
Compute time domain polarizations.
"""
# Intialize the model.
self.initialize_model(**parameters)
# Compute polarizations (evaluation in time array)
hp, hc, times = self.model.compute_polarizations(times=parameters.get("times", None))
# Return gwpy series
times = to_numpy(times)
return (
TimeSeries(to_numpy(hp), times=times, name="hp", unit=u.Unit("strain")),
TimeSeries(to_numpy(hc), times=times, name="hc", unit=u.Unit("strain")),
)
def _generate_fd_polarizations(self, **parameters):
"""
Compute Fourier domain polarizations.
"""
# Intialize the model.
self.initialize_model(**parameters)
# Compute polarizations in Fourier domain. Computes the conditioned TD polarizations and Fourier transform them
hp, hc, frequencies = self.model.compute_fd_polarizations() # FIXME currently it doesn't accept custom frequency array
# Retrun gwpy series
frequencies = to_numpy(frequencies)
return (
FrequencySeries(
to_numpy(hp), frequencies=frequencies, name="hp", unit=u.Unit("strain") * u.s, epoch=self.model.epoch
),
FrequencySeries(
to_numpy(hc), frequencies=frequencies, name="hc", unit=u.Unit("strain") * u.s, epoch=self.model.epoch
),
)
[docs]
class PyIMRPhenomTHM(PyIMRPhenomT):
"""
gwsignal generator wrapper for the IMRPhenomTHM model.
"""
def __init__(self, **kwargs):
super().__init__(**kwargs)
self._approximant = "IMRPhenomTHM"
self._class = getattr(phenomxpy, self._approximant)
[docs]
class PyIMRPhenomTP(PyIMRPhenomT):
"""
gwsignal generator wrapper for the IMRPhenomTP model.
"""
def __init__(self, **kwargs):
super().__init__(**kwargs)
self._approximant = "IMRPhenomTP"
self._class = getattr(phenomxpy, self._approximant)
[docs]
class PyIMRPhenomTPHM(PyIMRPhenomT):
"""
gwsignal generator wrapper for the IMRPhenomTPHM model.
"""
def __init__(self, **kwargs):
super().__init__(**kwargs)
self._approximant = "IMRPhenomTPHM"
self._class = getattr(phenomxpy, self._approximant)
[docs]
class PyIMRPhenomTFamily(PyIMRPhenomT):
"""
gwsignal generator wrapper for the IMRPhenomTFamily.
This class permits to use all the models in the IMRPhenomT family with the same generator.
One needs to specify the `approximant` name in the input parameters.
Example usage:
.. code-block:: python
gen = PyIMRPhenomTFamily()
params["approximant"] = "IMRPhenomTHM"
hp, hc = gws.core.waveform.GenerateTDWaveform(params, gen)
params["approximant"] = "IMRPhenomTPHM"
hp, hc = gws.core.waveform.GenerateTDWaveform(params, gen)
"""
def __init__(self, **kwargs):
super().__init__(**kwargs)
self._approximant = "IMRPhenomTFamily"
@property
def metadata(self):
return {
"type": "aligned_spin, precessing",
"f_ref_spin": True,
"modes": True,
"polarizations": True,
"implemented_domain": "time",
"approximant": "IMRPhenomTFamily",
"implementation": "",
"conditioning_routines": "",
}
def _evaluate_class(self, **parameters):
# Update parameters with options read from approximant name
# Returns approximant name without options, i.e. IMRPhenomT, IMRPhenomTHM, ...
approx_wo_options = parse_options_from_approximant_name(parameters)
# Choose appropiate phenomxpy class from approximant name
try:
self._class = getattr(phenomxpy, approx_wo_options)
except KeyError as e:
raise KeyError(f"Using PyIMRPhenomTFamily needs to provide `approximant` key in `parameters`.") from e
self.waveform_dict = convert_params_from_gwsignal(parameters)
return self._class(**self.waveform_dict)