# Copyright (C) 2023, 2024, 2025 Cecilio García Quirós
"""
Generator classes for the gwtemplate waveform interface.
This is the supported external interface for phenomxpy.
"""
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
from dataclasses import dataclass
[docs]
@dataclass
class MassTime:
"""Time samples in mass units, with seconds available via total_mass."""
value: object # numpy.ndarray or cupy.ndarray (times in units of M)
total_mass: float # total mass in solar masses
@property
def unit(self):
return u.Msun
@property
def seconds(self):
"""Same times in seconds (no astropy wrapper)."""
return MasstoSecond(self.value, self.total_mass)
[docs]
def to(self, unit):
unit = u.Unit(unit)
if unit == u.Msun:
return self.value * u.Msun
if unit == u.s:
return self.seconds * u.s
raise ValueError(f"Unsupported unit: {unit}")
def __len__(self):
return len(self.value)
def __getitem__(self, key):
return MassTime(self.value[key], self.total_mass)
[docs]
class PyIMRPhenomT(CompactBinaryCoalescenceGenerator):
"""
Generator classes for the gwtemplate waveform interface.
"""
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_time_domain_polarizations`, 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.
The epoch is the time of the peak of the 22 amplitude.
Returns
-------
MassTime
Start time of the currently initialized waveform model. `.seconds` attribute returns the start time in seconds. `.value` attribute returns the start time in mass units.
"""
if not hasattr(self, "model"):
raise RuntimeError("'epoch' is not available. Call initialize_model(**parameters) first.")
return MassTime(self.model.pWF.tmin, self.model.pWF.total_mass)
@property
def duration(self):
"""
Duration of the currently initialized waveform model.
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, although it can be evaluated beyond that time.
Returns
-------
MassTime
Duration of the currently initialized waveform model. `.seconds` attribute returns the duration in seconds. `.value` attribute returns the duration in mass units.
"""
if not hasattr(self, "model"):
raise RuntimeError("'duration' is not available. Call initialize_model(**parameters) first.")
if hasattr(self.model.pWF, "length"):
return MassTime(self.model.pWF.length * self.model.pWF.delta_t, self.model.pWF.total_mass)
return MassTime(self.model.pWF.tEnd - self.model.pWF.tmin, self.model.pWF.total_mass)
[docs]
def initialize_model(self, **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.
"""
parameters_checked = self.parameter_check(**parameters)
# Convert from astropy dictionary to phenomxpy dictionary
self.waveform_dict = convert_params_from_astropy(parameters_checked)
self.model = self._class(**self.waveform_dict)
return True
def _generate_time_domain_polarizations(self, times=None, deltaT=None):
"""
Compute time domain polarizations.
Parameters
----------
times: array-like or MassTime, optional
Mutually exclusive with ``deltaT``.
If omitted, ``deltaT`` must be supplied to build a uniform grid.
deltaT: astropy Quantity, optional
Uniform sampling interval. Mutually exclusive with ``times``.
If it matches the ``deltaT`` used at initialization, the model's
internal time grid is reused (avoiding a repeat ODE solve for
precessing models such as ``PyIMRPhenomTPHM``).
Returns
-------
tuple of (gwpy.timeseries.TimeSeries, gwpy.timeseries.TimeSeries)
``(hp, hc)`` as gwpy ``TimeSeries`` with unit ``strain``.
Sample times are in seconds.
"""
if times is None and deltaT is None:
raise ValueError("Either 'times' or 'deltaT' must be provided.")
if times is not None and deltaT is not None:
raise ValueError("Only one of 'times' or 'deltaT' must be provided.")
if times is not None and bool(self.model.pWF.condition) is True:
raise ValueError("Conditioning is not supported for custom time arrays")
if times is None:
if hasattr(self.model.pWF, "delta_t_sec") and self.model.pWF.delta_t_sec == deltaT.value:
# If deltaT is the same as the one used to initialize the model, we can use the time array stored in the model,
# this would make PyIMRPhenomTPHM faster since it won't need to solve again the ODE.
times = None
else:
# If deltaT is different than the one used to initialize the model, we need to generate a new time array
times = self.uniform_time_array(deltaT)
# Compute polarizations (evaluation in time array)
hp, hc, times = self.model.compute_polarizations(times=times)
# Return gwpy time 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_frequency_domain_polarizations(self):
"""
Compute Fourier domain polarizations.
Evaluates conditioned time-domain polarizations internally and
Fourier-transforms them. The frequency grid follows the model
initialization (``deltaF``, ``f_min``, ``f_max``); a custom frequency
array is not yet supported.
Returns
-------
tuple of (gwpy.frequencyseries.FrequencySeries, gwpy.frequencyseries.FrequencySeries)
``(hp, hc)`` as gwpy ``FrequencySeries`` with unit ``strain * s``. Sample frequencies are in Hz.
"""
# 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
# Return gwpy frequency 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]
def generate_time_domain_polarizations(self, **kwargs):
return self._generate_time_domain_polarizations(**kwargs)
[docs]
def generate_frequency_domain_polarizations(self):
return self._generate_frequency_domain_polarizations()
[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)