Using waveform interfaces

Examples for the supported external interfaces:

• phenomxpy.plugins.gwtemplate_plugin generator classes

• GeneratePolarizations with interface="gwtemplate"

%matplotlib inline

import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np

from phenomxpy.common import GeneratePolarizations
from phenomxpy.plugins.gwtemplate_plugin import PyIMRPhenomTPHM

params = {
    "mass1": 30.0 * u.solMass,
    "mass2": 10.0 * u.solMass,
    "spin1x": 0.0 * u.dimensionless_unscaled,
    "spin1y": 0.0 * u.dimensionless_unscaled,
    "spin1z": 0.5 * u.dimensionless_unscaled,
    "spin2x": 0.0 * u.dimensionless_unscaled,
    "spin2y": 0.0 * u.dimensionless_unscaled,
    "spin2z": -0.3 * u.dimensionless_unscaled,
    "deltaT": 1.0 / 512.0 * u.s,
    "f22_start": 20.0 * u.Hz,
    "f22_ref": 20.0 * u.Hz,
    "phi_ref": 0.0 * u.rad,
    "distance": 100.0 * u.Mpc,
    "inclination": 0.5 * u.rad,
}

gwtemplate plugin

gw_params = params.copy()

gen = PyIMRPhenomTPHM()
gen.initialize_model(**gw_params)

# Uniform time array
hp, hc = gen.generate_time_domain_polarizations(deltaT=gw_params["deltaT"])

# Generate polarization with custom time array
# pick random points from hp.times and sort them in ascending order
uniform_times = hp.times.value
custom_times = np.sort(np.random.choice(uniform_times[500:2000], size=1000))
hp_custom, hc_custom = gen.generate_time_domain_polarizations(times=custom_times)
# Check the length of the polarizations is the same as the input custom_times
assert len(hp_custom) == len(hc_custom) == len(custom_times)

common_idx = np.searchsorted(uniform_times, custom_times)
assert np.allclose(uniform_times[common_idx], custom_times)

hc_uniform_common = hc.value[common_idx]
hc_diff = np.abs(hc_uniform_common - hc_custom.value)

fig, axs = plt.subplots(1, 2, figsize=(12, 4))
axs[0].plot(uniform_times, hp.value, label="uniform")
axs[0].plot(custom_times, hp_custom.value, ".", label="custom")
axs[0].set_xlabel("time [s]")
axs[0].set_ylabel("h+")
axs[0].legend()

axs[1].plot(custom_times, hc_diff)
axs[1].set_xlabel("time [s]")
axs[1].set_ylabel("Absolute difference")

plt.tight_layout()
plt.show()

Chunked time-domain generation

Use generate_time_domain_polarizations_chunked to evaluate the waveform in memory-friendly sections. Pass deltaT for a uniform grid, or times for a custom (non-uniform) grid.

# Uniform grid: chunk by points_per_chunk
hp_chunked, hc_chunked = gen.generate_time_domain_polarizations_chunked(
    points_per_chunk=256,
    deltaT=gw_params["deltaT"],
)
assert len(hp_chunked) == len(hp)
assert np.allclose(hp.value, hp_chunked.value)

# Custom grid: pass times explicitly
window = uniform_times[500:2000]
select_in_window = np.sort(np.random.choice(len(window), size=min(200, len(window)), replace=False))
custom_times = window[select_in_window]

hp_custom, hc_custom = gen.generate_time_domain_polarizations(times=custom_times)
hp_custom_chunked, hc_custom_chunked = gen.generate_time_domain_polarizations_chunked(
    points_per_chunk=64,
    times=custom_times,
)
assert len(hp_custom_chunked) == len(custom_times)
assert np.allclose(hp_custom.value, hp_custom_chunked.value)

hp_uniform_err = np.abs(hp - hp_chunked)
hp_custom_err = np.abs(hp_custom - hp_custom_chunked)

fig, axs = plt.subplots(2, 2, figsize=(12, 7))

axs[0, 0].plot(hp, label="reference")
axs[0, 0].plot(hp_chunked, "--", label="chunked")
axs[0, 0].set_xlabel("time [s]")
axs[0, 0].set_ylabel("h+")
axs[0, 0].set_title("Uniform grid")
axs[0, 0].legend()

axs[0, 1].plot(hp_uniform_err)
axs[0, 1].set_xlabel("time [s]")
axs[0, 1].set_ylabel("Absolute error")
axs[0, 1].set_title("Uniform grid: absolute error")

axs[1, 0].plot(hp_custom, label="reference")
axs[1, 0].plot(hp_custom_chunked, "--", label="chunked")
axs[1, 0].set_xlabel("time [s]")
axs[1, 0].set_ylabel("h+")
axs[1, 0].set_title("Non-uniform grid")
axs[1, 0].legend()

axs[1, 1].plot(hp_custom_err)
axs[1, 1].set_xlabel("time [s]")
axs[1, 1].set_ylabel("Absolute error")
axs[1, 1].set_title("Non-uniform grid: absolute error")

plt.tight_layout()
plt.show()

GeneratePolarizations

gp_params = params.copy()
gp_params["interface"] = "gwtemplate"

# Uniform time array
hp, hc = GeneratePolarizations(parameters=gp_params, approximant="PyIMRPhenomTPHM", domain="TD")

# Custom time array passed in the parameters dictionary
uniform_times = hp.times.value
custom_times = np.sort(np.random.choice(uniform_times[500:2000], size=1000))

custom_params = gp_params.copy()
custom_params["times"] = custom_times
hp_custom, hc_custom = GeneratePolarizations(parameters=custom_params, approximant="PyIMRPhenomTPHM", domain="TD")
assert len(hp_custom) == len(hc_custom) == len(custom_times)

common_idx = np.searchsorted(uniform_times, custom_times)
assert np.allclose(uniform_times[common_idx], custom_times)

hc_uniform_common = hc.value[common_idx]
hc_diff = np.abs(hc_uniform_common - hc_custom.value)

fig, axs = plt.subplots(1, 2, figsize=(12, 4))
axs[0].plot(uniform_times, hp.value, label="uniform")
axs[0].plot(custom_times, hp_custom.value, ".", label="custom")
axs[0].set_xlabel("time [s]")
axs[0].set_ylabel("h+")
axs[0].legend()

axs[1].plot(custom_times, hc_diff)
axs[1].set_xlabel("time [s]")
axs[1].set_ylabel("Absolute difference")

plt.tight_layout()
plt.show()

Chunked time-domain generation

Pass chunk_options to GeneratePolarizations to evaluate the waveform in memory-friendly sections. Use the default uniform grid via deltaT in the parameters, or pass a custom times array for a non-uniform grid.

# Uniform grid: chunk by points_per_chunk
hp_chunked, hc_chunked = GeneratePolarizations(
    parameters=gp_params,
    approximant="PyIMRPhenomTPHM",
    domain="TD",
    chunk_options={"points_per_chunk": 256},
)
assert len(hp_chunked) == len(hp)
assert np.allclose(hp.value, hp_chunked.value)

# Custom grid: pass times in the parameters dictionary
window = uniform_times[500:2000]
select_in_window = np.sort(np.random.choice(len(window), size=min(200, len(window)), replace=False))
custom_times = window[select_in_window]

custom_params = gp_params.copy()
custom_params["times"] = custom_times

hp_custom, hc_custom = GeneratePolarizations(parameters=custom_params, approximant="PyIMRPhenomTPHM", domain="TD")
hp_custom_chunked, hc_custom_chunked = GeneratePolarizations(
    parameters=custom_params,
    approximant="PyIMRPhenomTPHM",
    domain="TD",
    chunk_options={"points_per_chunk": 64},
)
assert len(hp_custom_chunked) == len(custom_times)
assert np.allclose(hp_custom.value, hp_custom_chunked.value)

hp_uniform_err = np.abs(hp - hp_chunked)
hp_custom_err = np.abs(hp_custom - hp_custom_chunked)

fig, axs = plt.subplots(2, 2, figsize=(12, 7))

axs[0, 0].plot(hp, label="reference")
axs[0, 0].plot(hp_chunked, "--", label="chunked")
axs[0, 0].set_xlabel("time [s]")
axs[0, 0].set_ylabel("h+")
axs[0, 0].set_title("Uniform grid")
axs[0, 0].legend()

axs[0, 1].plot(hp_uniform_err)
axs[0, 1].set_xlabel("time [s]")
axs[0, 1].set_ylabel("Absolute error")
axs[0, 1].set_title("Uniform grid: absolute error")

axs[1, 0].plot(hp_custom, label="reference")
axs[1, 0].plot(hp_custom_chunked, "--", label="chunked")
axs[1, 0].set_xlabel("time [s]")
axs[1, 0].set_ylabel("h+")
axs[1, 0].set_title("Non-uniform grid")
axs[1, 0].legend()

axs[1, 1].plot(hp_custom_err)
axs[1, 1].set_xlabel("time [s]")
axs[1, 1].set_ylabel("Absolute error")
axs[1, 1].set_title("Non-uniform grid: absolute error")

plt.tight_layout()
plt.show()