# Copyright (C) 2023, 2024, 2025 Cecilio García Quirós
"""
Generator classes for a future gwsignal plugin interface.
Upstream LALSuite does not yet register external Python generators; a modified
LALSuite build is required to use this module through ``lalsimulation.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, MasstoSecond, SecondtoMass, to_numpy
from phenomxpy.common import convert_params_from_astropy
import numpy as np
[docs]
class PyIMRPhenomT(CompactBinaryCoalescenceGenerator):
"""
Generator class for a future gwsignal plugin interface.
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_astropy(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)