Walkthrough notebook

A guided tour of waveform generation, parameter derivatives, hybridisation and hybridised derivatives, run against the released model. It ships with the package as examples/usage/00_walkthrough.ipynb.

NRHJSur3dq8: waveforms, derivatives and hybridisation¶

A guided tour of the released surrogate through its public API. Everything here runs against the installed package and the released model file. There is no PYTHONPATH, no training file and no repository checkout involved, so this notebook doubles as a check that the release is self contained.

Four things are demonstrated, in order:

Waveform generation, in geometric and in physical units
Parameter derivatives, differentiated analytically by the model
Hybridisation, extending the inspiral below the stored content
Hybridised derivatives, which need a particular route

A fifth section shows the model's own uncertainty band, because a surrogate that cannot say how confident it is has only told you half of the answer.

The topic by topic scripts in this directory (01_load_and_modes.py through 10_finding_the_artifacts.py) go deeper on each lever.

In [1]:

%matplotlib inline

import warnings

# LAL prints a redirection notice on import under IPython. It is noise here.
warnings.filterwarnings("ignore", "Wswiglal-redir-stdio")

import matplotlib.pyplot as plt
import numpy as np

from nrhjsurrogate import NRHJSurAA

plt.rcParams.update({"figure.dpi": 110, "font.size": 9,
                     "axes.grid": True, "grid.alpha": 0.3})

model = NRHJSurAA.load_h5("NRHJSur3dq8_AA_v3")

print("modes served :", sorted(model.mode_names))

modes served : [(2, 1), (2, 2), (3, 2), (3, 3), (4, 3), (4, 4), (5, 5)]

1. Waveform generation¶

The geometric methods are the model's native language: time in units of the total mass $M$, and modes as $r\,h_{\ell m}/M$. Nothing is scaled, so this is the cheapest way to ask for content.

A binary is named by three numbers, the mass ratio and the two aligned spin components.

In [2]:

mass_ratio, spin1z, spin2z = 4.0, 0.3, -0.2

times, modes = model.get_td_waveform_modes_geometric(
    mass_ratio, spin1z, spin2z)

print("samples        :", times.size)
print("span (M)       : %.1f to %.1f" % (times[0], times[-1]))
print("h[(2,2)] shape :", modes[(2, 2)].shape)

samples        : 15775
span (M)       : -3819.1 to 124.4
h[(2,2)] shape : (15775,)

In [3]:

figure, axes = plt.subplots(2, 1, figsize=(7.0, 4.4), sharex=True)

axes[0].plot(times, modes[(2, 2)].real, lw=0.7, color="C0")
axes[0].plot(times, np.abs(modes[(2, 2)]), lw=1.0, color="C3")
axes[0].set_ylabel(r"$r\,h_{22}/M$")

for mode in [(2, 1), (3, 3), (4, 4)]:
    axes[1].semilogy(times, np.abs(modes[mode]), lw=0.8,
                     label=r"$(\ell,m)=(%d,%d)$" % mode)
axes[1].semilogy(times, np.abs(modes[(2, 2)]), lw=0.8, color="k",
                 label=r"$(\ell,m)=(2,2)$")
axes[1].set_ylim(1e-6, 1.0)
axes[1].set_ylabel(r"$|r\,h_{\ell m}/M|$")
axes[1].set_xlabel(r"$t/M$, peak of the $(2,2)$ amplitude at $t=0$")
axes[1].legend(ncol=4, loc="lower left", framealpha=1.0)

figure.tight_layout()
plt.show()

No description has been provided for this image

Physical units and the polarizations¶

get_td_waveform mirrors the LAL spelling: component masses in solar masses, a sample spacing and a starting frequency in Hz, a distance in Mpc and an inclination. It returns the two polarizations a detector would see, and like LAL it returns those alone, with the time array implied by delta_t. The geometric sibling get_td_waveform_geometric hands back the time array too, because there the grid is the model's own rather than one you asked for.

The two polarizations below look amplitude modulated, and that is geometry rather than precession. This model is non-eccentric and aligned-spin, so there is no precession in it at all. At an inclination between face-on and edge-on the polarizations are elliptical: $h_+$ carries a factor $(1+\cos^2\iota)/2$ and $h_\times$ a factor $\cos\iota$, so their combined envelope $|h_+ - i h_\times|$ breathes twice per orbit. It is strongest edge-on, where $h_\times$ vanishes entirely and the envelope drops to zero twice a cycle, and it disappears face-on, where the polarization is circular and the envelope is flat. Measured on this model, the peak-to-peak swing of the envelope about its trend is 0.006 face-on, 0.45 at 60 degrees and 1.78 edge-on, and a $(2,2)$-only waveform shows the same thing, which is the check that it is not higher-mode structure either.

In [4]:

delta_t = 1.0 / 4096.0

# Like LAL, this returns the polarizations alone: the time array is implied
# by delta_t, and the peak of the (2,2) amplitude sits at t = 0.
h_plus, h_cross = model.get_td_waveform(
    mass1=40.0, mass2=10.0, spin1z=0.3, spin2z=-0.2,
    delta_t=delta_t, f_min=30.0,
    inclination=np.pi / 3.0, distance=400.0)

seconds = np.arange(h_plus.size) * delta_t
seconds -= seconds[int(np.argmax(np.abs(h_plus - 1j * h_cross)))]

print("duration : %.3f s over %d samples" % (seconds[-1] - seconds[0],
                                             seconds.size))

figure, axis = plt.subplots(figsize=(7.0, 2.6))
axis.plot(seconds, h_plus, lw=0.7, label=r"$h_+$")
axis.plot(seconds, h_cross, lw=0.7, label=r"$h_\times$")
axis.set_xlabel("time (s)")
axis.set_ylabel("strain")
axis.legend(loc="upper left", framealpha=1.0)
figure.tight_layout()
plt.show()

duration : 0.627 s over 2568 samples

2. Parameter derivatives¶

The model differentiates its own stored content analytically. These are not finite differences of the waveform: the gradient kernel produces the waveform and its derivative in one fused pass, which is why asking for derivatives returns the waveform too.

derivatives[(l, m)] has shape (3, N), the three rows being $\partial/\partial q$, $\partial/\partial \chi_{1z}$ and $\partial/\partial \chi_{2z}$. The same numbers are available by parameter name through modes_grad_named, which returns a namedtuple per mode, so named[2, 2].spin1z is that row without remembering its index.

In [5]:

times, modes, derivatives = model.get_td_waveform_modes_geometric(
    mass_ratio, spin1z, spin2z, with_derivatives=True)

print("h[(2,2)]  shape :", modes[(2, 2)].shape)
print("dh[(2,2)] shape :", derivatives[(2, 2)].shape)

# The same derivatives by parameter name: a namedtuple per mode.
_, _, named = model.modes_grad_named(mass_ratio, spin1z, spin2z)
print("named fields    :", named[2, 2]._fields)
print("named is rows   :", np.array_equal(named[2, 2].spin1z,
                                          derivatives[(2, 2)][1]))

figure, axes = plt.subplots(3, 1, figsize=(7.0, 5.0), sharex=True)
labels = [r"$\partial h_{22}/\partial q$",
          r"$\partial h_{22}/\partial \chi_{1z}$",
          r"$\partial h_{22}/\partial \chi_{2z}$"]
rows = (named[2, 2].mass_ratio, named[2, 2].spin1z, named[2, 2].spin2z)
for axis, label, row in zip(axes, labels, rows):
    axis.plot(times, row.real, lw=0.7)
    axis.set_ylabel(label)
axes[-1].set_xlabel(r"$t/M$")
figure.tight_layout()
plt.show()

h[(2,2)]  shape : (15775,)
dh[(2,2)] shape : (3, 15775)

named fields    : ('mass_ratio', 'spin1z', 'spin2z')
named is rows   : True

After the peak the derivatives come from the merger-ringdown arm. In the figure above that stretch looks ragged, but that is the drawing: fifteen thousand samples across four thousand $M$ alias the fast merger cycles into what reads as noise. On its own axis the post-peak derivative is a damped oscillation under a smooth envelope, shown below with its magnitude. Note where each envelope starts: at the peak the amplitude is pinned, so the spin derivatives begin at zero there and grow into the ringdown, while the mass-ratio derivative does not.

In [6]:

after_peak = times > 0.0

figure, axes = plt.subplots(3, 1, figsize=(7.0, 5.0), sharex=True)
for axis, label, row in zip(axes, labels, rows):
    axis.plot(times[after_peak], row[after_peak].real, lw=0.8)
    axis.plot(times[after_peak], np.abs(row[after_peak]), lw=0.8,
              color="C3")
    axis.set_ylabel(label)
axes[0].legend(["real part", "magnitude"], loc="upper right",
               framealpha=1.0)
axes[-1].set_xlabel(r"$t/M$, after the peak only")
figure.tight_layout()
plt.show()

Checking the derivative against a finite difference¶

This is a check you can repeat, not a claim. Both waveforms are taken on the same fixed interior grid, because the inspiral's own domain moves with the mass ratio and comparing across two different grids would measure that instead.

In [7]:

check_grid = np.linspace(-3000.0, -50.0, 4096)
_, _, here = model.get_td_waveform_modes_geometric(
    mass_ratio, spin1z, spin2z, times=check_grid, with_derivatives=True)

step = 1e-4
_, plus, _ = model.get_td_waveform_modes_geometric(
    mass_ratio + step, spin1z, spin2z, times=check_grid,
    with_derivatives=True)
_, minus, _ = model.get_td_waveform_modes_geometric(
    mass_ratio - step, spin1z, spin2z, times=check_grid,
    with_derivatives=True)

finite_difference = (plus[(2, 2)] - minus[(2, 2)]) / (2.0 * step)
analytic = here[(2, 2)][0]
relative = np.abs(analytic - finite_difference) / np.abs(analytic).max()

print("median disagreement : %.3e" % np.median(relative))
print("worst  disagreement : %.3e" % relative.max())

median disagreement : 1.777e-05
worst  disagreement : 1.171e-04

3. Hybridisation¶

The stored artifact carries a finite stretch of inspiral. Ask to start below what it carries and the model extends the inspiral with a post-Newtonian description joined onto the surrogate content. Ask for a start already in reach and the waveform is trimmed there, with the overlapping content bit for bit unchanged.

The hybridization argument is a policy with three named states, so that a caller can declare an intent rather than hope:

state	meaning
`"extend_if_below_reach"`	the default: trim when in reach, extend when below
`"never"`	refuse to extend, so the content stays numerical-relativity based
`"require"`	refuse unless an extension actually happens

In [8]:

native_times, native_modes = model.get_td_waveform_modes_geometric(
    mass_ratio, spin1z, spin2z)
print("native span (M) : %.1f to %.1f" % (native_times[0], native_times[-1]))

for start_frequency in (0.004, 0.002, 0.0015):
    extended_times, _ = model.get_td_waveform_modes_geometric(
        mass_ratio, spin1z, spin2z, f_low=start_frequency)
    print("M f = %.4f     : %7d samples, starts at %9.1f M"
          % (start_frequency, extended_times.size, extended_times[0]))

# In reach, f_low only trims: the overlap is untouched.
trimmed_times, trimmed_modes = model.get_td_waveform_modes_geometric(
    mass_ratio, spin1z, spin2z, f_low=0.01)
overlap = np.isin(native_times, trimmed_times)
print("in-reach trim is exact :",
      np.array_equal(trimmed_modes[(2, 2)], native_modes[(2, 2)][overlap]))

native span (M) : -3819.1 to 124.4

M f = 0.0040     :   55881 samples, starts at  -13845.5 M

M f = 0.0020     :  367642 samples, starts at  -91785.9 M

M f = 0.0015     :  796583 samples, starts at -199021.0 M
in-reach trim is exact : True

In [9]:

# Each policy refuses the case it was named to refuse.
for frequency, policy in [(0.001, "never"), (0.02, "require")]:
    try:
        model.get_td_waveform_modes_geometric(
            mass_ratio, spin1z, spin2z, f_low=frequency,
            hybridization=policy)
    except ValueError as error:
        print("%-8s at M f = %.3f : %s" % (policy, frequency, error))

never    at M f = 0.001 : f_low=0.001 is below this point's native start (M f = 0.006311) and hybridization='never' guarantees NR-only content; use hybridization='extend_if_below_reach' or 'require' to engage the PN extension arm

require  at M f = 0.020 : hybridization='require' but f_low=0.02 is inside the native span at (mass_ratio, spin1z, spin2z) = (4, 0.3, -0.2); the PN extension arm would not engage there; use hybridization='extend_if_below_reach' (which trims the native-span reconstruction to f_low)

In [10]:

extended_times, extended_modes = model.get_td_waveform_modes_geometric(
    mass_ratio, spin1z, spin2z, f_low=0.0015)

# The extension is 50x longer than the stored content, so a linear time
# axis would squash the stored part into the axis line. Plotting against
# time BEFORE the peak on a log axis shows both, with the merger on the
# right where it belongs.
figure, axis = plt.subplots(figsize=(7.0, 2.9))
axis.plot(-extended_times, np.abs(extended_modes[(2, 2)]), lw=0.8,
          color="C1", label="extended to $Mf = 0.0015$")
axis.plot(-native_times, np.abs(native_modes[(2, 2)]), lw=1.4, color="C0",
          label="stored content")
axis.axvline(-native_times[0], color="k", ls=":", lw=0.8)
axis.annotate("edge of the stored content", xy=(-native_times[0], 6.4e-2),
              xytext=(700.0, 2.6e-2), fontsize=8,
              arrowprops=dict(arrowstyle="->", lw=0.6))
axis.set_xscale("log")
axis.set_yscale("log")
axis.invert_xaxis()
axis.set_xlim(3e5, 3.0)
axis.set_ylim(1.5e-2, 4e-1)
axis.set_xlabel(r"time before the peak, $-t/M$ (merger to the right)")
axis.set_ylabel(r"$|r\,h_{22}/M|$")
axis.legend(loc="upper left", framealpha=1.0)
figure.tight_layout()
plt.show()

4. Hybridised derivatives¶

Derivatives and hybridisation do not compose on the default route. Ask for both and the model refuses, and names the entry point that does the job rather than quietly handing back a narrower answer:

In [11]:

try:
    model.get_td_waveform_modes_geometric(
        mass_ratio, spin1z, spin2z, f_low=0.004, with_derivatives=True)
except NotImplementedError as error:
    print("default route :", error)

default route : f_low gradients on timing_mode='aa': use hybrid_gradients.hybrid_modes_grad

hybrid_modes_grad is that entry point. Two properties of it matter before you use the result:

it is inspiral only. There is no merger or ringdown in the returned arrays, by construction and not by accident.
it works on the raw hybrid clock, running up to the amplitude peak at $t \approx 0$, rather than the peak-aligned clock the plain geometric methods use.

In [12]:

from nrhjsurrogate.driver.hybrid_gradients import hybrid_modes_grad

# The call warns about the caveat discussed below. Captured and printed as
# text so it reads as part of the story rather than as a stack trace.
with warnings.catch_warnings(record=True) as raised:
    warnings.simplefilter("always")
    hybrid_times, hybrid_modes, hybrid_derivatives = hybrid_modes_grad(
        model, (mass_ratio, spin1z, spin2z), f_low=0.005)

for warning in raised:
    print("WARNED :", str(warning.message).split(". ")[0], "...")
print()

print("samples         :", hybrid_times.size)
print("span (M)        : %.1f to %.1f  (inspiral only)"
      % (hybrid_times[0], hybrid_times[-1]))
print("dh[(2,2)] shape :", hybrid_derivatives[(2, 2)].shape)

figure, axes = plt.subplots(2, 1, figsize=(7.0, 4.0), sharex=True)
axes[0].plot(hybrid_times, hybrid_modes[(2, 2)].real, lw=0.5)
axes[0].set_ylabel(r"$r\,h_{22}/M$")
axes[1].plot(hybrid_times, hybrid_derivatives[(2, 2)][0].real, lw=0.5,
             color="C3")
axes[1].set_ylabel(r"$\partial h_{22}/\partial q$")
axes[1].set_xlabel(r"$t/M$, raw hybrid clock")
figure.tight_layout()
plt.show()

WARNED : PN_RATIO_POWERLAW (NRHJ_PN_RATIO_POWERLAW) shape the extension the VALUES path serves but are not implemented in the gradient chain (hybrid_gradients.extension_transport and its stage-C mode assembly), so dh is the derivative of a slightly different waveform than h -- measured pre-unification at 1.7e-3..5.2e-3 of peak, worst on (4,4); the mode-assembly share remains ...

samples         : 29699
span (M)        : -7424.6 to -0.3  (inspiral only)
dh[(2,2)] shape : (3, 29699)

Two cautions worth reading before a Fisher matrix¶

The extension's derivative is not exactly the derivative of the extension. Some shaping applied to the extended values is not implemented in the gradient chain, so dh is the derivative of a slightly different waveform than h. The model measures the gap at roughly 0.2 to 0.5 percent of peak, worst on the (4,4) mode, and warns at the call site rather than letting you discover it downstream. Using h and dh together is inconsistent by about that much, systematically, and there is no second gradient arm to fall back to. The warning is diagnostic and can be turned into an error with NRHJ_GRADIENT_FLAG_GUARD=raise.

The "aa-int-pn" route returns a different clock with derivatives than without. Values alone come back peak aligned; the same call with with_derivatives=True comes back on the generator's own clock, offset by the full span. The waveform content is identical, only the labelling of time differs, so realigning is a subtraction. Do it before combining the result with anything taken on the peak-aligned clock.

In [13]:

values_times, _ = model.get_td_waveform_modes_geometric(
    mass_ratio, spin1z, spin2z, timing_mode="aa-int-pn")
derivative_times, _, _ = model.get_td_waveform_modes_geometric(
    mass_ratio, spin1z, spin2z, timing_mode="aa-int-pn",
    with_derivatives=True)

print("values only : %9.1f to %9.1f M" % (values_times[0], values_times[-1]))
print("derivatives : %9.1f to %9.1f M" % (derivative_times[0],
                                          derivative_times[-1]))

realigned = derivative_times + values_times[0]
print("realigned   : %9.1f to %9.1f M" % (realigned[0], realigned[-1]))
print("matches the values clock exactly :",
      np.array_equal(realigned, values_times))

values only :   -3818.0 to       0.0 M
derivatives :       0.0 to    3818.0 M
realigned   :   -3818.0 to       0.0 M
matches the values clock exactly : True

5. The model's own uncertainty¶

get_waveform_error_estimate returns, for every served mode, a 1-sigma amplitude uncertainty in the mode's own units and a phase uncertainty in radians.

Read the caveat with the number. This is the surrogate's own predictive uncertainty, a statement from its fits. It is not a measurement of how far the model sits from numerical relativity, and outside the training region it is an extrapolation that should be read as a lower bound. The merger-ringdown part is the model's weakest, and an envelope that is correctly computed there is a different statement from an arm that is accurate there.

Read the two returned arrays as one number, not two. The pair is built as (envelope, envelope / |h|), so the "phase" array is exactly the relative error of the mode. That is the usual convention of charging the whole error budget to one channel or the other: a small complex perturbation of size $|\delta h|$ is either an amplitude error of $|\delta h|$ or, equivalently, a phase error of $|\delta h| / |h|$ radians. It is one envelope wearing two hats, and adding an amplitude term and a phase term into the same likelihood would count it twice.

Across most of the inspiral that band is far thinner than a plotted line, so the relative form is the readable one. The band itself is drawn further down, near the peak and deliberately exaggerated, so its shape is visible at all.

In [14]:

samples, error_estimate = model.get_waveform_error_estimate(
    mass_ratio, spin1z, spin2z)
delta_amplitude, delta_phase = error_estimate[(2, 2)]

_, band_modes = model.get_td_waveform_modes_geometric(
    mass_ratio, spin1z, spin2z, times=samples)
amplitude = np.abs(band_modes[(2, 2)])

# Across the whole span the absolute band is far thinner than the plotted
# line, which is itself the headline: plot the FRACTIONAL amplitude
# uncertainty instead, so the number is visible and directly readable.
# The two returned arrays are ONE quantity in two dresses, exactly:
print("delta_phase == delta_amplitude / amplitude :",
      np.array_equal(delta_phase, delta_amplitude / amplitude))

figure, axis = plt.subplots(figsize=(7.0, 2.7))
axis.semilogy(samples, delta_phase, lw=0.8, color="C0")
axis.axvline(0.0, color="k", ls=":", lw=0.8)
axis.set_ylabel("relative $1\\sigma$\n(= radians of phase)")
axis.set_xlabel(r"$t/M$")
figure.tight_layout()
plt.show()

inspiral = samples < -100.0
print("relative 1-sigma, inspiral median : %.2e"
      % np.median(delta_phase[inspiral]))
print("relative 1-sigma, at the peak     : %.2e"
      % delta_phase[np.argmin(np.abs(samples))])

# The band IS there, it is simply narrow. Near the peak, on a linear axis:
near_peak = (samples > -300.0)
figure, axis = plt.subplots(figsize=(7.0, 2.4))
axis.plot(samples[near_peak], amplitude[near_peak], lw=0.9, color="C0",
          label="amplitude")
axis.fill_between(samples[near_peak],
                  (amplitude - 50.0 * delta_amplitude)[near_peak],
                  (amplitude + 50.0 * delta_amplitude)[near_peak],
                  alpha=0.35, color="C0", lw=0,
                  label=r"$1\sigma$, exaggerated $\times 50$")
axis.set_xlabel(r"$t/M$")
axis.set_ylabel(r"$|r\,h_{22}/M|$")
axis.legend(loc="upper left", framealpha=1.0)
figure.tight_layout()
plt.show()

delta_phase == delta_amplitude / amplitude : True

relative 1-sigma, inspiral median : 4.73e-04
relative 1-sigma, at the peak     : 6.14e-04

6. Comparing against another model¶

A surrogate's own uncertainty says what its fits think of themselves. An independent check is a different model. NRHybSur3dq8 covers the same non-eccentric aligned-spin corner of parameter space and is built from the same numerical-relativity catalogue by a different method, so it is the natural comparison.

This section needs lalsimulation and the NRHybSur3dq8 data file on LAL_DATA_PATH. If either is missing the section reports that and moves on.

Compare modes, not polarizations. Two bookkeeping conventions differ between the two models and neither is a disagreement about physics:

a constant orbital phase offset, because each model puts zero orbital phase in a different place. It enters mode $m$ as $e^{-im\Delta\phi}$, so removing one constant per waveform removes it everywhere.
a per-mode factor $e^{i(m-2)\pi/2}$ between the two mode conventions. It is unity on $(2,2)$, a quarter turn on the odd-$m$ modes and a sign flip on $(4,4)$.

The second does not reduce to a single phase, so it does not cancel out of a sum over modes. Comparing mode by mode, with one constant phase removed per mode, sidesteps both cleanly.

In [15]:

try:
    import lal
    import lalsimulation as lalsim
    comparison_available = True
except ImportError as error:
    comparison_available = False
    print("skipping: lalsimulation is not installed (%s)" % error)


def aligned_overlap(times_a, mode_a, times_b, mode_b, sample_spacing):
    '''Overlap of two modes, maximised over a time shift and one phase.

    Each is put on its own clock with the amplitude peak at zero, then both
    are interpolated onto a shared grid. The maximisation is the standard
    one: the cross-correlation gives the best shift, and its complex
    argument at that shift gives the phase.
    '''
    times_a = times_a - times_a[int(np.argmax(np.abs(mode_a)))]
    times_b = times_b - times_b[int(np.argmax(np.abs(mode_b)))]
    start = max(times_a[0], times_b[0]) + 0.02
    stop = min(times_a[-1], times_b[-1]) - 0.005
    grid = np.arange(start, stop, sample_spacing)
    a = (np.interp(grid, times_a, mode_a.real)
         + 1j * np.interp(grid, times_a, mode_a.imag))
    b = (np.interp(grid, times_b, mode_b.real)
         + 1j * np.interp(grid, times_b, mode_b.imag))
    correlation = np.fft.ifft(np.fft.fft(a) * np.conj(np.fft.fft(b)))
    best = int(np.argmax(np.abs(correlation)))
    norm = np.sqrt((np.abs(a) ** 2).sum() * (np.abs(b) ** 2).sum())
    return np.abs(correlation[best]) / norm, np.angle(correlation[best])

In [16]:

if comparison_available:
    mass1, mass2, spin1, spin2 = 40.0, 10.0, 0.3, -0.2
    delta_t, f_min, distance = 1.0 / 4096.0, 30.0, 400.0

    approximant = lalsim.SimInspiralGetApproximantFromString("NRHybSur3dq8")
    try:
        node = lalsim.SimInspiralChooseTDModes(
            0.0, delta_t, mass1 * lal.MSUN_SI, mass2 * lal.MSUN_SI,
            0.0, 0.0, spin1, 0.0, 0.0, spin2,
            f_min, f_min, distance * 1e6 * lal.PC_SI, None, 5, approximant)
    except Exception as error:
        node = None
        print("skipping: NRHybSur3dq8 could not be generated (%s)" % error)
        print("the data file must be on LAL_DATA_PATH")

    reference_modes = {}
    while node is not None:
        reference_times = (np.arange(node.mode.data.length) * delta_t
                           + float(node.mode.epoch))
        reference_modes[(node.l, node.m)] = np.array(node.mode.data.data)
        node = node.next

if comparison_available and reference_modes:
    our_times, our_modes = model.get_td_waveform_modes(
        mass1=mass1, mass2=mass2, spin1z=spin1, spin2z=spin2,
        delta_t=delta_t, f_min=f_min, distance=distance)

    print("mode      overlap     mismatch")
    mismatches = {}
    for mode in sorted(our_modes):
        if mode not in reference_modes:
            continue
        overlap, _ = aligned_overlap(our_times, our_modes[mode],
                                     reference_times, reference_modes[mode],
                                     delta_t)
        mismatches[mode] = 1.0 - overlap
        print("  %-7s %.6f    %.2e" % (str(mode), overlap, 1.0 - overlap))

mode      overlap     mismatch
  (2, 1)  0.999978    2.22e-05
  (2, 2)  0.999969    3.10e-05
  (3, 2)  0.999897    1.03e-04
  (3, 3)  0.999915    8.47e-05
  (4, 3)  0.999699    3.01e-04
  (4, 4)  0.999567    4.33e-04
  (5, 5)  0.999565    4.35e-04

In [17]:

if comparison_available and mismatches:
    figure, axis = plt.subplots(figsize=(7.0, 2.7))
    labels = [r"$(%d,%d)$" % m for m in mismatches]
    axis.bar(range(len(mismatches)), list(mismatches.values()),
             color="C0", width=0.6)
    axis.set_xticks(range(len(mismatches)))
    axis.set_xticklabels(labels)
    axis.set_yscale("log")
    axis.set_ylabel("mismatch against\nNRHybSur3dq8")
    axis.grid(axis="x", visible=False)
    figure.tight_layout()
    plt.show()

    print("worst mode : %s at %.2e"
          % (max(mismatches, key=mismatches.get), max(mismatches.values())))
    print("%d + %d Msun, spins (%.1f, %.1f), from %.0f Hz"
          % (mass1, mass2, spin1, spin2, f_min))

worst mode : (5, 5) at 4.35e-04
40 + 10 Msun, spins (0.3, -0.2), from 30 Hz

Every served mode agrees with the other model to a few parts in $10^4$ or better, and the $(2,2)$ to a few parts in $10^5$. Two models built from the same simulations by different methods will not agree perfectly, and the gap above is a bound on the pair rather than an error attributable to either one: saying which is closer to the truth needs numerical relativity as the arbiter, not the other model.

If you compare polarizations instead¶

The polarizations this model returns are in LAL's convention: the line of sight sits at azimuth $\pi/2 - \phi_{\rm ref}$ in the spin-weighted harmonic and the modes are projected in LAL's source frame, exactly as SimInspiralChooseTDWaveform(inclination, phiRef) does. So get_td_waveform(inclination=i, phi_ref=p) and LAL's call with inclination=i, phiRef=p describe the same observer, and comparing the two needs only the time alignment, with no per-mode phase to remove.

Two things are worth knowing. First, a polarization comparison with the absolute phase held fixed does not reach the per-mode figure above, which freed one phase per mode: the section below measures it at these parameters, and the residual is the two codes' anchoring of the reference phase plus the $m=0$ modes LAL includes, not the convention. The convention itself is exact: LAL's own modes through this model's projection reproduce LAL's own polarizations to $10^{-15}$. Second, the earlier convention is still there for anyone who needs bit-for-bit continuity with it: source_frame="native" together with azimuth=0.0 reproduces the pre-convention output exactly.

If you would rather not think about conventions at all, do what the model's own validation campaigns do: take modes from both models and project them yourself through one function, so whatever convention it uses cancels between the two.

The polarizations themselves, against NRHybSur3dq8¶

The per-mode figure above freed one phase per mode. This is the stricter test, at the parameters of the polarization plot in section 1: both models are asked for the same observer, the inclination and reference phase passed to each call verbatim, and the two strains compared with only a time shift between them and no phase freed at all. That is what a detector would see, and it is the test a convention mismatch cannot pass. The same comparison with one global phase freed is shown alongside, because the difference between the two numbers is exactly the disagreement in where the reference phase was anchored.

In [18]:

def aligned_polarization_overlap(times_a, strain_a, times_b, strain_b,
                                 sample_spacing):
    '''Overlap of two complex strains h_+ - i h_x, maximised over a time
    shift with the phase held fixed, and again with one global phase freed.

    Returns (overlap with the phase fixed, overlap with the phase freed,
    the freed phase in radians, the time shift for each, and the two clocks
    with the amplitude peak at zero) so the strains can be drawn aligned.'''
    times_a = times_a - times_a[int(np.argmax(np.abs(strain_a)))]
    times_b = times_b - times_b[int(np.argmax(np.abs(strain_b)))]
    start = max(times_a[0], times_b[0]) + 0.02
    stop = min(times_a[-1], times_b[-1]) - 0.005
    grid = np.arange(start, stop, sample_spacing)
    a = (np.interp(grid, times_a, strain_a.real)
         + 1j * np.interp(grid, times_a, strain_a.imag))
    b = (np.interp(grid, times_b, strain_b.real)
         + 1j * np.interp(grid, times_b, strain_b.imag))
    correlation = np.fft.ifft(np.fft.fft(a) * np.conj(np.fft.fft(b)))
    norm = np.sqrt((np.abs(a) ** 2).sum() * (np.abs(b) ** 2).sum())

    def as_seconds(index):
        return (index - grid.size if index > grid.size // 2 else index) \
            * sample_spacing

    fixed = int(np.argmax(correlation.real))
    freed = int(np.argmax(np.abs(correlation)))
    return (correlation.real[fixed] / norm,
            np.abs(correlation[freed]) / norm, np.angle(correlation[freed]),
            as_seconds(fixed), as_seconds(freed), times_a, times_b)


def compare_polarizations(f_ref):
    '''Both models at the observer of the section 1 plot, reference
    frequency f_ref. LAL is asked for exactly the same observer.'''
    inclination, phi_ref = np.pi / 3.0, 0.0
    lal_plus, lal_cross = lalsim.SimInspiralChooseTDWaveform(
        mass1 * lal.MSUN_SI, mass2 * lal.MSUN_SI,
        0.0, 0.0, spin1, 0.0, 0.0, spin2,
        distance * 1e6 * lal.PC_SI, inclination, phi_ref,
        0.0, 0.0, 0.0, delta_t, f_min, f_ref, None, approximant)
    lal_times = (np.arange(lal_plus.data.length) * delta_t
                 + float(lal_plus.epoch))
    lal_strain = (np.array(lal_plus.data.data)
                  - 1j * np.array(lal_cross.data.data))
    our_plus, our_cross = model.get_td_waveform(
        mass1=mass1, mass2=mass2, spin1z=spin1, spin2z=spin2,
        delta_t=delta_t, f_min=f_min, f_ref=f_ref,
        inclination=inclination, phi_ref=phi_ref, distance=distance)
    our_times = np.arange(our_plus.size) * delta_t
    our_strain = our_plus - 1j * our_cross
    result = aligned_polarization_overlap(our_times, our_strain,
                                          lal_times, lal_strain, delta_t)
    return dict(our_strain=our_strain, lal_strain=lal_strain,
                overlap_fixed=result[0], overlap_freed=result[1],
                freed_phase=result[2], shift_fixed=result[3],
                shift_freed=result[4], our_clock=result[5],
                lal_clock=result[6])


if comparison_available and reference_modes:
    # f_ref = f_min is what the section 1 plot used (the default); 40 Hz is
    # a reference frequency inside the band.
    comparisons = {f_ref: compare_polarizations(f_ref)
                   for f_ref in (f_min, 40.0)}
    print("polarization mismatch against NRHybSur3dq8, inclination pi/3, "
          "phi_ref 0")
    print("  f_ref    phase held fixed   one phase freed   freed phase")
    for f_ref, c in comparisons.items():
        print("  %3.0f Hz   %.2e           %.2e          %+.3f rad"
              % (f_ref, 1.0 - c["overlap_fixed"], 1.0 - c["overlap_freed"],
                 c["freed_phase"]))
    print("%d + %d Msun, spins (%.1f, %.1f), f_min %.0f Hz, 1/delta_t %.0f Hz"
          % (mass1, mass2, spin1, spin2, f_min, 1.0 / delta_t))

polarization mismatch against NRHybSur3dq8, inclination pi/3, phi_ref 0
  f_ref    phase held fixed   one phase freed   freed phase
   30 Hz   1.81e-02           2.41e-03          +0.278 rad
   40 Hz   2.51e-03           3.73e-04          +0.065 rad
40 + 10 Msun, spins (0.3, -0.2), f_min 30 Hz, 1/delta_t 4096 Hz

In [19]:

if comparison_available and reference_modes:
    drawn = comparisons[f_min]
    clock = drawn["our_clock"]
    window = (clock > -0.12) & (clock < 0.02)
    our_plus = drawn["our_strain"].real
    lal_time_aligned = np.interp(
        clock, drawn["lal_clock"] + drawn["shift_fixed"],
        drawn["lal_strain"].real)
    lal_phase_aligned = np.interp(
        clock, drawn["lal_clock"] + drawn["shift_freed"],
        (drawn["lal_strain"] * np.exp(1j * drawn["freed_phase"])).real)

    figure, axes = plt.subplots(3, 1, figsize=(7.0, 6.2), sharex=True)
    axes[0].plot(clock[window], our_plus[window], lw=1.0, color="C0",
                 label="NRHJSur3dq8, this model")
    axes[0].plot(clock[window], lal_time_aligned[window], lw=0.8,
                 color="k", ls="--", label="NRHybSur3dq8, through LAL")
    axes[0].set_ylabel(r"$h_+$")
    axes[0].legend(loc="upper left", framealpha=1.0)

    axes[1].plot(clock[window], (our_plus - lal_time_aligned)[window],
                 lw=0.8, color="C3")
    axes[1].set_ylabel(r"$h_+$ difference," + "\n" + "phase held fixed")

    axes[2].plot(clock[window], (our_plus - lal_phase_aligned)[window],
                 lw=0.8, color="C3")
    axes[2].set_ylabel(r"$h_+$ difference," + "\n" + "one phase freed")
    axes[2].set_xlabel(r"time (s), peak of $|h_+ - i h_\times|$ at 0")
    for axis in axes[1:]:
        axis.set_ylim(axes[0].get_ylim())
    figure.tight_layout()
    plt.show()

With the phase held fixed the two strains at the section 1 parameters sit $1.8 \times 10^{-2}$ apart, and freeing one global phase takes that to $2.4 \times 10^{-3}$. The middle panel shows what the number means: the inspiral agrees, and the merger cycles are displaced by a fraction of a cycle. The freed phase, 0.28 rad, is the disagreement between the two codes about where the reference phase is anchored when the reference frequency is the starting frequency, which is the default in both. With the reference frequency moved inside the band, to 40 Hz, the same comparison gives $2.5 \times 10^{-3}$ with the phase held fixed and the freed phase drops to 0.07 rad, which is one sample of orbital phase at this sampling rate. The mode-level agreement is the same at every reference frequency, so the whole polarization-level gap is that one orbital phase. Set f_ref inside the band when the absolute phase matters, and compare modes when it does not.

Where to go next¶

01_load_and_modes.py, 02_physical_units.py for the loading and unit levers
03_hybridisation.py, 05_parameter_derivatives.py for the full argument set
06_ada_model.py for the AdA family, which shares the merger-ringdown arm
07_fisher_matrix.py for using the derivatives in a Fisher forecast
08_model_error.py for the uncertainty in more depth
09_backends.py for the compiled evaluator, 10_finding_the_artifacts.py for where the model files live

nrhjsur-verify runs the shipped reference checks against your installation.