Example usage of NRSur7dq4v2 surrogate model.¶
In [1]:
Copied!
import numpy as np
import matplotlib.pyplot as P
%matplotlib inline
import gwsurrogate
import numpy as np
import matplotlib.pyplot as P
%matplotlib inline
import gwsurrogate
Download surrogate data, this only needs to be done once¶
In [2]:
Copied!
# This can take a few minutes
gwsurrogate.catalog.pull('NRSur7dq4v2')
# This can take a few minutes
gwsurrogate.catalog.pull('NRSur7dq4v2')
Out[2]:
'/home/aravishankar_umassd_edu/sxs/gwsurrogate/gwsurrogate/surrogate_downloads/NRSur7dq4v2.h5?rlkey=mci9eq1y71p3jwjsxrlgcmbps&st=amueac9c&dl=1'
Load the surrogate, this only needs to be done once at the start of a script¶
In [3]:
Copied!
sur = gwsurrogate.LoadSurrogate('NRSur7dq4v2')
sur = gwsurrogate.LoadSurrogate('NRSur7dq4v2')
Loaded NRSur7dq4v2 model
Read the documentation¶
In [4]:
Copied!
help(sur)
help(sur)
Help on NRSur7dq4v2 in module gwsurrogate.surrogate object:
class NRSur7dq4v2(SurrogateEvaluator)
| NRSur7dq4v2(h5filename, basis_tol_opts=None)
|
| A class for the NRSur7dq4v2 surrogate model, a domain-decomposed modification
| of the surrogate presented in Varma et al. 2019, arxiv1905.09300.
|
| Evaluates gravitational waveforms generated by precessing binary black hole
| systems with generic mass ratios and spins.
|
| This model includes the following spin-weighted spherical harmonic modes:
| 2<=ell<=5, -ell<=m<=ell.
|
| The parameter space of validity is:
| q \in [1, 6], and |chi1|,|chi2| \in [-1, 1], with generic directions.
| where q is the mass ratio and chi1/chi2 are the spin vectors of the
| heavier/lighter BH, respectively.
|
| The surrogate has been trained in the range
| q \in [1, 4] and |chi1|/|chi2| \in [-0.8, 0.8], but produces reasonable
| waveforms in the above range and has been tested against existing
| NR waveforms in that range.
|
| See the __call__ method on how to evaluate waveforms.
|
| Method resolution order:
| NRSur7dq4v2
| SurrogateEvaluator
| builtins.object
|
| Methods defined here:
|
| __init__(self, h5filename, basis_tol_opts=None)
| name: Name of the surrogate
| domain_type: 'Time' or 'Frequency'
| keywords: keywords for this model. For allowed keys see
| self._check_keywords_and_set_defaults.default_keywords.
| If keywords['Precessing'] = False, will automatically
| determine the m<0 modes from the m>0 modes.
| soft_param_lims: Parameter bounds beyond which a warning is raised.
| hard_param_lims: Parameter bounds beyond which an error is raised.
| Should be in format [qMax, chimax]
| Setting soft_param_lims/hard_param_lims to None will
| skip that particular check.
|
| coorbital_basis_sizes(self, ellMax=5)
| Returns the number of basis functions used in the surrogate fits for
| each datapiece. This is useful when wanting to reduce dimensionality
| of the surrogate fits by truncating the number of basis functions
| during evaluation.
|
| ----------------------------------------------------------------------
| Methods inherited from SurrogateEvaluator:
|
| __call__(self, q, chiA0, chiB0, M=None, dist_mpc=None, f_low=None, f_ref=None, dt=None, df=None, times=None, freqs=None, mode_list=None, ellMax=None, inclination=None, phi_ref=0, precessing_opts=None, tidal_opts=None, par_dict=None, units='dimensionless', skip_param_checks=False, taper_end_duration=None)
| INPUT
| =====
| q : Mass ratio, mA/mB >= 1.
| chiA0: Dimensionless spin vector of the heavier black hole at
| reference epoch.
| chiB0: Dimensionless spin vector of the lighter black hole at
| reference epoch.
|
| This follows the same convention as LAL, where the spin
| components are defined as:
| \chi_z = \chi \cdot \hat{L}, where L is the orbital angular
| momentum vector at the epoch.
| \chi_x = \chi \cdot \hat{n}, where n = body2 -> body1 is the
| separation vector at the epoch. body1 is the heavier body.
| \chi_y = \chi \cdot \hat{L \cross n}.
| These spin components are frame-independent as they are
| defined using vector inner products. This is equivalent to
| specifying the spins in the coorbital frame used in the
| surrogate papers.
|
| M, dist_mpc: Either specify both M and dist_mpc or neither.
| M : Total mass (solar masses). Default: None.
| dist_mpc : Distance to binary system (MegaParsecs). Default: None.
|
| f_low : Instantaneous initial frequency of the (2, 2) mode. In
| practice, this is estimated to be twice the initial orbital
| frequency in the coprecessing frame. Note: the coprecessing
| frame is the minimal rotation frame of arXiv:1110.2965.
|
| f_low should be in cycles/M if units = 'dimensionless',
| should be in Hertz if units = 'mks'.
| If 0, the entire waveform is returned.
| Default: None, must be specified by user.
|
| NOTE: For some models like NRSur7dq4, f_low=0 is recommended.
| The role of f_low is only to truncate the lower frequencies
| before returning the waveform. Since this model is already
| very short, this truncation is not required. On the other hand,
| f_ref is used to set the reference epoch, and can be freely
| specified.
|
| WARNING: Using f_low=0 with a small dt (like 0.1M) can lead to
| very expensive evaluation for hybridized surrogates like
| NRHybSur3dq8.
|
| f_ref: Frequency used to set the reference epoch at which the
| reference frame is defined and the spins are specified.
| See below for definition of the reference frame.
| Should be in cycles/M if units = 'dimensionless', should be
| in Hertz if units = 'mks'.
| Default: If f_ref is not given, we set f_ref = f_low. If
| f_low is 0, this corresponds to the initial index.
|
| For time domain models, f_ref is used to determine a t_ref,
| such that the orbital frequency in the coprecessing frame
| equals f_ref/2 at t=t_ref.
|
| dt, df : Time/Frequency step size, specify at most one of dt/df,
| depending on whether the surrogate is a time/frequency domain
| surrogate.
| Default: None. If None, the internal domain of the surrogate is
| used, which can be nonuniformly sampled.
| dt (df) Should be in M (cycles/M) if units = 'dimensionless',
| should be in seconds (Hertz) if units = 'mks'. Do not specify
| times/freqs if using dt/df.
|
|
| times, freqs:
| Array of time/frequency samples at which to evaluate the
| waveform, depending on whether the surrogate is a
| time/frequency domain surrogate. time (freqs) should be in
| M (cycles/M) if units = 'dimensionless', should be in
| seconds (Hertz) if units = 'mks'. Do not specify dt/df if
| using times/freqs. Default None.
|
| ellMax: Maximum ell index for modes to include. All available m
| indicies for each ell will be included automatically. The
| m<0 modes will automatically be included for nonprecessing
| models.
| Default: None, in which case all available ells will be
| included.
|
| mode_list : A list of (ell, m) modes tuples to be included. Valid only
| for nonprecessing models.
|
| Example: mode_list = [(2,2),(2,1)].
| Default: None, in which case all available modes are included.
|
| At most one of ellMax and mode_list can be specified.
|
| Note: mode_list is allowed only for nonprecessing models; for
| precessing models use ellMax. For precessing systems, all m
| indices of a given ell index mix with each other, so there is
| no clear hierarchy. To get the individual modes just don't
| specify inclination and a dictionary of modes will be returned.
|
| Note: When the inclination is set, the m<0 modes are
| automatically included. For example, passing mode_list = [(2,2)]
| will include the (2,2) and (2,-2) modes in the computation of
| the strain.
|
| Note: When the inclination is None, the m<0 modes are
| automatically generated.
|
| inclination : Inclination angle between the orbital angular momentum
| direction at the reference epoch and the line-of-sight to the
| observer. If inclination is None, the mode data is returned
| as a dictionary.
| Default: None.
|
| phi_ref : The azimuthal angle on the sky of the source frame following
| the LAL convention.
| Default: 0.
|
| If inclination/phi_ref are specified, the complex strain (h =
| hplus -i hcross) evaluated at (inclination, pi/2 - phi_ref) on
| the sky of the reference frame is returned. This follows the
| same convention as LAL. See below for definition of the
| reference frame.
|
| precessing_opts:
| A dictionary containing optional parameters for a precessing
| surrogate model. Default: None.
| Allowed keys are:
| init_orbphase: The orbital phase in the coprecessing frame
| at the reference epoch.
| Default: 0, in which case the coorbital frame and
| coprecessing frame are the same.
| init_quat: The unit quaternion (length 4 vector) giving the
| rotation from the coprecessing frame to the inertial frame
| at the reference epoch.
| Default: None, in which case the coprecessing frame is the
| same as the inertial frame.
| return_dynamics:
| Return the frame dynamics and spin evolution along with
| the waveform. Default: False.
| Example: precessing_opts = {
| 'init_orbphase': 0,
| 'init_quat': [1,0,0,0],
| 'return_dynamics': True
| }
|
| tidal_opts:
| A dictionary containing optional parameters for a tidal
| surrogate model. Default: None.
| Allowed keys are:
| Lambda1: The tidal deformability parameter for the heavier
| object.
| Lambda2: The tidal deformability parameter for the lighter
| object.
| Example: tidal_opts = {'Lambda1': 200, 'Lambda2': 300}
|
|
| par_dict: A dictionary containing any additional parameters needed for a
| particular surrogate model. Default: None.
|
| units: 'dimensionless' or 'mks'. Default: 'dimensionless'.
| If 'dimensionless': Any of f_low, f_ref, dt, df, times and
| freqs, if specified, must be in dimensionless units. That
| is, dt/times should be in units of M, while f_ref, f_low
| and df/freqs should be in units of cycles/M.
| M and dist_mpc must be None. The waveform and domain are
| returned as dimensionless quantities as well.
| If 'mks': Any of f_low, f_ref, dt, df, times and freqs, if
| specified, must be in MKS units. That is, dt/times should
| be in seconds, while f_ref, f_low and df/freqs should be
| in Hz. M and dist_mpc must be specified. The waveform and
| domain are returned in MKS units as well.
|
|
| skip_param_checks :
| Skip sanity checks for inputs. Use this if you want to
| extrapolate outside allowed range. Default: False.
|
| taper_end_duration:
| Taper the last TAPER_END_DURATION (M) of a time-domain waveform
| in units of M. For exmple, passing 40 will taper the last 40M.
| When set to None, no taper is applied
| Default: None.
|
| RETURNS
| =====
|
| domain, h, dynamics
|
|
| domain : Array of time/frequency samples corresponding to h and
| dynamics, depending on whether the surrogate is a
| time/frequency domain model. This is the same as times/freqs
| if times/freqs are given as an inputs.
| For time domain models the time is set to 0 at the peak of
| the waveform. The time (frequency) values are in M (cycles/M)
| if units = 'dimensionless', they are in seconds (Hertz) if
| units = 'mks'
|
| h : The waveform.
| If inclination is specified, the complex strain (h = hplus
| -i hcross) evaluated at (inclination, pi/2 - phi_ref) on
| the sky of the reference frame is returned. This follows
| the LAL convention, see below for details. This includes
| all modes given in the ellMax/mode_list argument. For
| nonprecessing systems the m<0 modes are automatically
| deduced from the m>0 modes. To see if a model is precessing
| check self.keywords.
|
| Else, h is a dictionary of available modes with (l, m)
| tuples as keys. For example, h22 = h[(2,2)].
|
| If M and dist_mpc are given, the physical waveform
| at that distance is returned. Else, it is returned in
| code units: r*h/M extrapolated to future null-infinity.
|
| dynamics: A dict containing the frame dynamics and spin evolution. This
| is None for nonprecessing models. This is also None if
| return_dynamics in precessing_opts is False (Default).
|
| The dynamics include (L=len(domain)):
|
| q_copr = dynamics['q_copr']
| The quaternion representing the coprecessing frame with
| shape (4, L)
| orbphase = dynamics['orbphase']
| The orbital phase in the coprecessing frame with length L.
| chiA = dynamics['chiA']
| The inertial frame chiA with shape (L, 3)
| chiB = dynamics['chiB']
| The inertial frame chiB with shape (L, 3)
|
|
| IMPORTANT NOTES:
| ===============
|
| The reference frame (or inertial frame) is defined as follows:
| The +ve z-axis is along the orbital angular momentum at the reference
| epoch. The separation vector from the lighter BH to the heavier BH at
| the reference epoch is along the +ve x-axis. The y-axis completes the
| right-handed triad. The reference epoch is set using f_ref.
|
| Now, if inclination/phi_ref are given, the waveform is evaluated at
| (inclination, pi/2 - phi_ref) in the reference frame. This agrees with
| the LAL convention. See LIGO DCC document T1800226 for the LAL frame
| diagram.
|
| ----------------------------------------------------------------------
| Data descriptors inherited from SurrogateEvaluator:
|
| __dict__
| dictionary for instance variables (if defined)
|
| __weakref__
| list of weak references to the object (if defined)
Evaluate the waveform¶
Evaluate waveform modes in dimensionless units (default)¶
In [5]:
Copied!
q = 4
chiA = [-0.2, 0.4, 0.1]
chiB = [-0.5, 0.2, -0.4]
dt = 0.1 # timestep size, Units of M
f_low = 0 # initial frequency, f_low=0 returns the full surrogate
t, h, dyn = sur(q, chiA, chiB, dt=dt, f_low=f_low) # dyn stands for dynamics, do dyn.keys() to see contents
q = 4
chiA = [-0.2, 0.4, 0.1]
chiB = [-0.5, 0.2, -0.4]
dt = 0.1 # timestep size, Units of M
f_low = 0 # initial frequency, f_low=0 returns the full surrogate
t, h, dyn = sur(q, chiA, chiB, dt=dt, f_low=f_low) # dyn stands for dynamics, do dyn.keys() to see contents
In [6]:
Copied!
# Let's see all available modes
print( sorted(h.keys()) )
# Let's see all available modes
print( sorted(h.keys()) )
[(2, -2), (2, -1), (2, 0), (2, 1), (2, 2), (3, -3), (3, -2), (3, -1), (3, 0), (3, 1), (3, 2), (3, 3), (4, -4), (4, -3), (4, -2), (4, -1), (4, 0), (4, 1), (4, 2), (4, 3), (4, 4), (5, -5), (5, -4), (5, -3), (5, -2), (5, -1), (5, 0), (5, 1), (5, 2), (5, 3), (5, 4), (5, 5)]
In [7]:
Copied!
P.plot(t, h[(2,2)].real, label='l2m2 real')
P.plot(t, h[(2,1)].real, label='l2m1 real')
P.plot(t, h[(3,3)].real, label='l3m3 real')
P.ylabel(r'Re[$h_{lm}$]', fontsize=18)
P.xlabel(r't [M]', fontsize=18)
P.legend()
P.plot(t, h[(2,2)].real, label='l2m2 real')
P.plot(t, h[(2,1)].real, label='l2m1 real')
P.plot(t, h[(3,3)].real, label='l3m3 real')
P.ylabel(r'Re[$h_{lm}$]', fontsize=18)
P.xlabel(r't [M]', fontsize=18)
P.legend()
Out[7]:
<matplotlib.legend.Legend at 0x741ceb0ceaa0>
Evaluate waveform modes in physical units¶
In [8]:
Copied!
q = 4
chiA = [-0.2, 0.4, 0.1]
chiB = [-0.5, 0.2, -0.4]
f_ref = 20 # Reference frequecny in Hz. The spins are assumed to specified at this frequency
f_low = 0 # initial frequency, f_low=0 returns the full surrogate
M = 70 # Total masss in solar masses
dist_mpc = 100 # distance in megaparsecs
dt = 1./4096 # step size in seconds
ellMax = 4 # Highest ell index for modes to use
# dyn stands for dynamics, do dyn.keys() to see contents
t, h, dyn = sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, ellMax=ellMax, M=M, dist_mpc=dist_mpc, units='mks')
P.plot(t, h[(2,2)].real, label='l2m2 real')
P.plot(t, h[(2,1)].real, label='l2m1 real')
P.ylabel(r'Re[$h_{lm}$]', fontsize=18)
P.xlabel(r't [s]', fontsize=18)
P.legend()
q = 4
chiA = [-0.2, 0.4, 0.1]
chiB = [-0.5, 0.2, -0.4]
f_ref = 20 # Reference frequecny in Hz. The spins are assumed to specified at this frequency
f_low = 0 # initial frequency, f_low=0 returns the full surrogate
M = 70 # Total masss in solar masses
dist_mpc = 100 # distance in megaparsecs
dt = 1./4096 # step size in seconds
ellMax = 4 # Highest ell index for modes to use
# dyn stands for dynamics, do dyn.keys() to see contents
t, h, dyn = sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, ellMax=ellMax, M=M, dist_mpc=dist_mpc, units='mks')
P.plot(t, h[(2,2)].real, label='l2m2 real')
P.plot(t, h[(2,1)].real, label='l2m1 real')
P.ylabel(r'Re[$h_{lm}$]', fontsize=18)
P.xlabel(r't [s]', fontsize=18)
P.legend()
Out[8]:
<matplotlib.legend.Legend at 0x741cdbc1e260>
Evaluate waveform at a point on the sky¶
In [9]:
Copied!
q = 4
chiA = [-0.2, 0.4, 0.1]
chiB = [-0.5, 0.2, -0.4]
f_ref = 20 # Reference frequecny in Hz. The spins are assumed to specified at this frequency
f_low = 0 # initial frequency, f_low=0 returns the full surrogate
M = 70 # Total masss in solar masses
dist_mpc = 100 # distance in megaparsecs
dt = 1./4096 # step size in seconds
ellMax = 4 # Highest ell index for modes to use
inclination = np.pi/4
phi_ref = np.pi/5
# Will only include modes ell<=ellMax
# Returns h = h_+ -i h_x at (inclination, phi_ref) in the sky of the source frame
# dyn stands for dynamics, do dyn.keys() to see contents
t, h, dyn = sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, ellMax=ellMax, M=M, dist_mpc=dist_mpc,
inclination=inclination, phi_ref=phi_ref, units='mks')
P.plot(t, h.real)
P.ylabel(r'$h_{+}$ $(\iota, \phi_{ref})$', fontsize=18)
P.xlabel(r't [s]', fontsize=18)
q = 4
chiA = [-0.2, 0.4, 0.1]
chiB = [-0.5, 0.2, -0.4]
f_ref = 20 # Reference frequecny in Hz. The spins are assumed to specified at this frequency
f_low = 0 # initial frequency, f_low=0 returns the full surrogate
M = 70 # Total masss in solar masses
dist_mpc = 100 # distance in megaparsecs
dt = 1./4096 # step size in seconds
ellMax = 4 # Highest ell index for modes to use
inclination = np.pi/4
phi_ref = np.pi/5
# Will only include modes ell<=ellMax
# Returns h = h_+ -i h_x at (inclination, phi_ref) in the sky of the source frame
# dyn stands for dynamics, do dyn.keys() to see contents
t, h, dyn = sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, ellMax=ellMax, M=M, dist_mpc=dist_mpc,
inclination=inclination, phi_ref=phi_ref, units='mks')
P.plot(t, h.real)
P.ylabel(r'$h_{+}$ $(\iota, \phi_{ref})$', fontsize=18)
P.xlabel(r't [s]', fontsize=18)
Out[9]:
Text(0.5, 0, 't [s]')
Dynamics Surrogate¶
In [10]:
Copied!
q = 4
chiA = [-0.2, 0.4, 0.1] # unless f_ref is given the spins are assumed to be given at the start of the waveform
chiB = [-0.5, 0.2, -0.4]
dt = 0.1 # step size, Units of M
f_low = 0 # initial frequency, f_low=0 returns the full surrogate
t, h, dyn = sur(q, chiA, chiB, dt=dt, f_low=f_low, precessing_opts={'return_dynamics': True}) # dyn stands for dynamics, do dyn.keys() to see contents
q = 4
chiA = [-0.2, 0.4, 0.1] # unless f_ref is given the spins are assumed to be given at the start of the waveform
chiB = [-0.5, 0.2, -0.4]
dt = 0.1 # step size, Units of M
f_low = 0 # initial frequency, f_low=0 returns the full surrogate
t, h, dyn = sur(q, chiA, chiB, dt=dt, f_low=f_low, precessing_opts={'return_dynamics': True}) # dyn stands for dynamics, do dyn.keys() to see contents
In [11]:
Copied!
# See all available dynamics data
dyn.keys()
# See all available dynamics data
dyn.keys()
Out[11]:
dict_keys(['chiA', 'chiB', 'chiA_copr', 'chiB_copr', 'q_copr', 'orbphase'])
In [12]:
Copied!
P.figure(1)
P.plot(t, dyn['chiA'][:,0], label=r'$\chi_{Ax}$')
P.plot(t, dyn['chiA'][:,1], label=r'$\chi_{Ay}$')
P.plot(t, dyn['chiA'][:,2], label=r'$\chi_{Az}$')
P.ylabel(r'$\chi_A$', fontsize=18)
P.xlabel(r't', fontsize=18)
P.title('Spin of heavier BH')
P.legend(fontsize=14)
P.figure(2)
P.plot(t, dyn['chiB'][:,0], label=r'$\chi_{Bx}$')
P.plot(t, dyn['chiB'][:,1], label=r'$\chi_{By}$')
P.plot(t, dyn['chiB'][:,2], label=r'$\chi_{Bz}$')
P.ylabel(r'$\chi_B$', fontsize=18)
P.xlabel(r't', fontsize=18)
P.title('Spin of lighter BH')
P.legend(fontsize=14)
P.figure(3)
P.plot(t, dyn['orbphase'][:,])
P.ylabel(r'$\phi_{\mathrm{orb}}$', fontsize=18)
P.xlabel(r't', fontsize=18)
P.title('Orbital phase')
P.figure(4)
P.plot(t, dyn['q_copr'][0,:], label=r'$\hat{Q}_0$')
P.plot(t, dyn['q_copr'][1,:], label=r'$\hat{Q}_1$')
P.plot(t, dyn['q_copr'][2,:], label=r'$\hat{Q}_2$')
P.plot(t, dyn['q_copr'][3,:], label=r'$\hat{Q}_3$')
P.ylabel(r'$\hat{Q}$', fontsize=18)
P.xlabel(r't', fontsize=18)
P.title('Coprecessing frame quaternions')
P.legend(fontsize=14)
P.figure(1)
P.plot(t, dyn['chiA'][:,0], label=r'$\chi_{Ax}$')
P.plot(t, dyn['chiA'][:,1], label=r'$\chi_{Ay}$')
P.plot(t, dyn['chiA'][:,2], label=r'$\chi_{Az}$')
P.ylabel(r'$\chi_A$', fontsize=18)
P.xlabel(r't', fontsize=18)
P.title('Spin of heavier BH')
P.legend(fontsize=14)
P.figure(2)
P.plot(t, dyn['chiB'][:,0], label=r'$\chi_{Bx}$')
P.plot(t, dyn['chiB'][:,1], label=r'$\chi_{By}$')
P.plot(t, dyn['chiB'][:,2], label=r'$\chi_{Bz}$')
P.ylabel(r'$\chi_B$', fontsize=18)
P.xlabel(r't', fontsize=18)
P.title('Spin of lighter BH')
P.legend(fontsize=14)
P.figure(3)
P.plot(t, dyn['orbphase'][:,])
P.ylabel(r'$\phi_{\mathrm{orb}}$', fontsize=18)
P.xlabel(r't', fontsize=18)
P.title('Orbital phase')
P.figure(4)
P.plot(t, dyn['q_copr'][0,:], label=r'$\hat{Q}_0$')
P.plot(t, dyn['q_copr'][1,:], label=r'$\hat{Q}_1$')
P.plot(t, dyn['q_copr'][2,:], label=r'$\hat{Q}_2$')
P.plot(t, dyn['q_copr'][3,:], label=r'$\hat{Q}_3$')
P.ylabel(r'$\hat{Q}$', fontsize=18)
P.xlabel(r't', fontsize=18)
P.title('Coprecessing frame quaternions')
P.legend(fontsize=14)
Out[12]:
<matplotlib.legend.Legend at 0x741cdabbe770>
Comparison with NRSur7dq4¶
Since NRSur7dq4v2 can be used identically to NRSur7dq4 on the modes they model in common, we compare them here
In [ ]:
Copied!
# This can take a few minutes -- only needed if you don't have NRSur7dq4 already
gwsurrogate.catalog.pull('NRSur7dq4')
# This can take a few minutes -- only needed if you don't have NRSur7dq4 already
gwsurrogate.catalog.pull('NRSur7dq4')
In [13]:
Copied!
surv1 = gwsurrogate.LoadSurrogate('NRSur7dq4')
q = 4
chiA = [-0.2, 0.4, 0.1]
chiB = [-0.5, 0.2, -0.4]
dt = 0.1
f_low = 0
t, hv1, dynv1 = surv1(q, chiA, chiB, dt=dt, f_low=f_low, precessing_opts={'return_dynamics': True})
surv1 = gwsurrogate.LoadSurrogate('NRSur7dq4')
q = 4
chiA = [-0.2, 0.4, 0.1]
chiB = [-0.5, 0.2, -0.4]
dt = 0.1
f_low = 0
t, hv1, dynv1 = surv1(q, chiA, chiB, dt=dt, f_low=f_low, precessing_opts={'return_dynamics': True})
Loaded NRSur7dq4 model
In [14]:
Copied!
fig0, axes0 = P.subplots(2,2, figsize=(14,10))
# Comparing real parts of the (2,2) mode from both surrogates
axes0[0,0].plot(t, h[(2,2)].real, label='NRSur7dq4v2')
axes0[0,0].plot(t, hv1[(2,2)].real, '--' , label='NRSur7dq4')
axes0[0,0].set_ylabel(r'Re[$h_{22}$]', fontsize=18)
axes0[0,0].set_xlabel(r't [M]', fontsize=18)
axes0[0,0].set_title('Surrogate Waveforms')
axes0[0,0].legend(fontsize=14)
# Plotting absolute difference of real parts
axes0[0,1].plot(t, h[(2,2)].real-hv1[(2,2)].real)
axes0[0,1].set_ylabel(r'Difference', fontsize=18)
axes0[0,1].set_xlabel(r't [M]', fontsize=18)
axes0[0,1].set_title('Difference between Surrogates Waveforms')
# Comparing imaginary parts of the (2,2) mode from both surrogates
axes0[1,0].plot(t, h[(2,2)].imag, label='NRSur7dq4v2')
axes0[1,0].plot(t, hv1[(2,2)].imag, '--' , label='NRSur7dq4')
axes0[1,0].set_ylabel(r'Im[$h_{22}$]', fontsize=18)
axes0[1,0].set_xlabel(r't [M]', fontsize=18)
axes0[1,0].set_title('Surrogate Waveforms')
axes0[1,0].legend(fontsize=14)
# Plotting absolute difference of imag parts
axes0[1,1].plot(t, h[(2,2)].imag-hv1[(2,2)].imag)
axes0[1,1].set_ylabel(r'Difference', fontsize=18)
axes0[1,1].set_xlabel(r't [M]', fontsize=18)
axes0[1,1].set_title('Difference between Surrogates Waveforms')
fig0.suptitle('Comparison of $h_{22}$ from NRSur7dq4v2 and NRSur7dq4', fontsize=16)
P.tight_layout()
P.show()
fig0, axes0 = P.subplots(2,2, figsize=(14,10))
# Comparing real parts of the (2,2) mode from both surrogates
axes0[0,0].plot(t, h[(2,2)].real, label='NRSur7dq4v2')
axes0[0,0].plot(t, hv1[(2,2)].real, '--' , label='NRSur7dq4')
axes0[0,0].set_ylabel(r'Re[$h_{22}$]', fontsize=18)
axes0[0,0].set_xlabel(r't [M]', fontsize=18)
axes0[0,0].set_title('Surrogate Waveforms')
axes0[0,0].legend(fontsize=14)
# Plotting absolute difference of real parts
axes0[0,1].plot(t, h[(2,2)].real-hv1[(2,2)].real)
axes0[0,1].set_ylabel(r'Difference', fontsize=18)
axes0[0,1].set_xlabel(r't [M]', fontsize=18)
axes0[0,1].set_title('Difference between Surrogates Waveforms')
# Comparing imaginary parts of the (2,2) mode from both surrogates
axes0[1,0].plot(t, h[(2,2)].imag, label='NRSur7dq4v2')
axes0[1,0].plot(t, hv1[(2,2)].imag, '--' , label='NRSur7dq4')
axes0[1,0].set_ylabel(r'Im[$h_{22}$]', fontsize=18)
axes0[1,0].set_xlabel(r't [M]', fontsize=18)
axes0[1,0].set_title('Surrogate Waveforms')
axes0[1,0].legend(fontsize=14)
# Plotting absolute difference of imag parts
axes0[1,1].plot(t, h[(2,2)].imag-hv1[(2,2)].imag)
axes0[1,1].set_ylabel(r'Difference', fontsize=18)
axes0[1,1].set_xlabel(r't [M]', fontsize=18)
axes0[1,1].set_title('Difference between Surrogates Waveforms')
fig0.suptitle('Comparison of $h_{22}$ from NRSur7dq4v2 and NRSur7dq4', fontsize=16)
P.tight_layout()
P.show()
Note the larger differences in the merger and ringdown, which results from the independent ringdown modeling of NRSur7dq4v2.
Now let us check for the same in the dynamics data.
In [15]:
Copied!
fig1, axes1 = P.subplots(3,2, figsize=(14,12))
# Comparing x-components of the spin of the heavier BH from both surrogates
axes1[0,0].plot(t, dyn['chiA'][:,0], label='NRSur7dq4v2')
axes1[0,0].plot(t, dynv1['chiA'][:,0], label='NRSur7dq4')
axes1[0,0].set_ylabel(r'$\chi_{Ax}$', fontsize=18)
axes1[0,0].set_xlabel(r't [M]', fontsize=18)
axes1[0,0].set_title(r'$x$-component of spin of heavier BH')
axes1[0,0].legend(fontsize=14)
# Plotting absolute difference of x-components
axes1[0,1].plot(t, dyn['chiA'][:,0]-dynv1['chiA'][:,0])
axes1[0,1].set_ylabel(r'Difference', fontsize=18)
axes1[0,1].set_xlabel(r't [M]', fontsize=18)
axes1[0,1].set_title(r'Difference between $x$-components')
# Comparing orbital phases from both surrogates
axes1[1,0].plot(t, dyn['orbphase'][:,], label='NRSur7dq4v2')
axes1[1,0].plot(t, dynv1['orbphase'][:,], label='NRSur7dq4')
axes1[1,0].set_ylabel(r'$\phi_{\mathrm{orb}}$', fontsize=18)
axes1[1,0].set_xlabel(r't [M]', fontsize=18)
axes1[1,0].set_title(r'Orbital phase')
axes1[1,0].legend(fontsize=14)
# Plotting absolute difference of orbital phases
axes1[1,1].plot(t, dyn['orbphase'][:,]-dynv1['orbphase'][:,])
axes1[1,1].set_ylabel(r'Difference', fontsize=18)
axes1[1,1].set_xlabel(r't [M]', fontsize=18)
axes1[1,1].set_title(r'Difference between orbital phases')
# Comparing first component of the coprecessing frame quaternion from both surrogates
axes1[2,0].plot(t, dyn['q_copr'][0,:], label='NRSur7dq4v2')
axes1[2,0].plot(t, dynv1['q_copr'][0,:], label='NRSur7dq4')
axes1[2,0].set_ylabel(r'$\hat{Q}_0$', fontsize=18)
axes1[2,0].set_xlabel(r't [M]', fontsize=18)
axes1[2,0].set_title(r'Coprecessing frame quaternion component $\hat{Q}_0$')
axes1[2,0].legend(fontsize=14)
# Plotting absolute difference of first component of the coprecessing frame quaternion
axes1[2,1].plot(t, dyn['q_copr'][0,:]-dynv1['q_copr'][0,:])
axes1[2,1].set_ylabel(r'Difference', fontsize=18)
axes1[2,1].set_xlabel(r't [M]', fontsize=18)
axes1[2,1].set_title(r'Difference between quaternion components')
fig1.suptitle('Comparison of Dynamics from NRSur7dq4v2 and NRSur7dq4', fontsize=16)
P.tight_layout()
P.show()
fig1, axes1 = P.subplots(3,2, figsize=(14,12))
# Comparing x-components of the spin of the heavier BH from both surrogates
axes1[0,0].plot(t, dyn['chiA'][:,0], label='NRSur7dq4v2')
axes1[0,0].plot(t, dynv1['chiA'][:,0], label='NRSur7dq4')
axes1[0,0].set_ylabel(r'$\chi_{Ax}$', fontsize=18)
axes1[0,0].set_xlabel(r't [M]', fontsize=18)
axes1[0,0].set_title(r'$x$-component of spin of heavier BH')
axes1[0,0].legend(fontsize=14)
# Plotting absolute difference of x-components
axes1[0,1].plot(t, dyn['chiA'][:,0]-dynv1['chiA'][:,0])
axes1[0,1].set_ylabel(r'Difference', fontsize=18)
axes1[0,1].set_xlabel(r't [M]', fontsize=18)
axes1[0,1].set_title(r'Difference between $x$-components')
# Comparing orbital phases from both surrogates
axes1[1,0].plot(t, dyn['orbphase'][:,], label='NRSur7dq4v2')
axes1[1,0].plot(t, dynv1['orbphase'][:,], label='NRSur7dq4')
axes1[1,0].set_ylabel(r'$\phi_{\mathrm{orb}}$', fontsize=18)
axes1[1,0].set_xlabel(r't [M]', fontsize=18)
axes1[1,0].set_title(r'Orbital phase')
axes1[1,0].legend(fontsize=14)
# Plotting absolute difference of orbital phases
axes1[1,1].plot(t, dyn['orbphase'][:,]-dynv1['orbphase'][:,])
axes1[1,1].set_ylabel(r'Difference', fontsize=18)
axes1[1,1].set_xlabel(r't [M]', fontsize=18)
axes1[1,1].set_title(r'Difference between orbital phases')
# Comparing first component of the coprecessing frame quaternion from both surrogates
axes1[2,0].plot(t, dyn['q_copr'][0,:], label='NRSur7dq4v2')
axes1[2,0].plot(t, dynv1['q_copr'][0,:], label='NRSur7dq4')
axes1[2,0].set_ylabel(r'$\hat{Q}_0$', fontsize=18)
axes1[2,0].set_xlabel(r't [M]', fontsize=18)
axes1[2,0].set_title(r'Coprecessing frame quaternion component $\hat{Q}_0$')
axes1[2,0].legend(fontsize=14)
# Plotting absolute difference of first component of the coprecessing frame quaternion
axes1[2,1].plot(t, dyn['q_copr'][0,:]-dynv1['q_copr'][0,:])
axes1[2,1].set_ylabel(r'Difference', fontsize=18)
axes1[2,1].set_xlabel(r't [M]', fontsize=18)
axes1[2,1].set_title(r'Difference between quaternion components')
fig1.suptitle('Comparison of Dynamics from NRSur7dq4v2 and NRSur7dq4', fontsize=16)
P.tight_layout()
P.show()
Again, we see a large difference in the ringdown. This is due to the tapering of the coprecessing frame in the ringdown of NRSur7dq4v2.
Basis-size reduction¶
Basis size reduction can be performed for specific data pieces by loading a model preset with gwsurrogate.LoadSurrogate(). We first describe the basis-size dictionary format, then show advanced users how to register a custom model preset, and finally demonstrate the packaged "Fast" preset.
Basis-size specification format¶
A basis-size preset is a dictionary containing <datapiece-name>: <basis-size>, pairs. For example, custom_basis_sizes = {'2_2_Re+_sd_0': 16, '3_0_imag_sd_1': 3}
Every data piece has an inspiral subdomain size indicated by sd_0 and a ringdown subdomain basis size indicated by sd_1. Every mode with $m=0$ has two kinds of datapieces: imag (imaginary part) and real (real part). Every mode with $m\neq0$ has four kinds of datapieces: Im+, Im-, Re+, Re-, where $\pm$ indicates the sign of the linear combination of $h_{\ell, m}$ and $h_{\ell,-m}$ modes used.
Further details regarding these datapieces and basis truncation can be found in arXiv:2609.07873.
Basis size inquiry¶
In order to reduce basis sizes for chosen datapieces, one must first know the basis size used by NRSur7dq4v2 for these datapieces, as these serve as the upper limit for the basis sizes that can be requested. This can be done using the surrogate's coorbital_basis_sizes() method, which accepts an optional argument ellMax which can be used to specify the maximum $\ell$ value up to which basis sizes must be recovered.
In [16]:
Copied!
max_basis_sizes = sur.coorbital_basis_sizes(ellMax=3)
print(r"Max basis sizes for datapieces up to $\ell=3$: ")
for dp_name, size in max_basis_sizes.items():
print(f"{dp_name}: {size}")
max_basis_sizes = sur.coorbital_basis_sizes(ellMax=3)
print(r"Max basis sizes for datapieces up to $\ell=3$: ")
for dp_name, size in max_basis_sizes.items():
print(f"{dp_name}: {size}")
Max basis sizes for datapieces up to $\ell=3$: 2_0_real_sd_0: 26 2_0_imag_sd_0: 33 2_0_real_sd_1: 12 2_0_imag_sd_1: 20 2_1_Re+_sd_0: 36 2_1_Re+_sd_1: 25 2_1_Re-_sd_0: 24 2_1_Re-_sd_1: 36 2_1_Im+_sd_0: 44 2_1_Im+_sd_1: 28 2_1_Im-_sd_0: 21 2_1_Im-_sd_1: 17 2_2_Re+_sd_0: 64 2_2_Re+_sd_1: 6 2_2_Re-_sd_0: 24 2_2_Re-_sd_1: 4 2_2_Im+_sd_0: 5 2_2_Im+_sd_1: 4 2_2_Im-_sd_0: 27 2_2_Im-_sd_1: 14 3_0_real_sd_0: 37 3_0_imag_sd_0: 27 3_0_real_sd_1: 13 3_0_imag_sd_1: 24 3_1_Re+_sd_0: 24 3_1_Re+_sd_1: 42 3_1_Re-_sd_0: 35 3_1_Re-_sd_1: 30 3_1_Im+_sd_0: 22 3_1_Im+_sd_1: 40 3_1_Im-_sd_0: 38 3_1_Im-_sd_1: 18 3_2_Re+_sd_0: 27 3_2_Re+_sd_1: 10 3_2_Re-_sd_0: 30 3_2_Re-_sd_1: 9 3_2_Im+_sd_0: 33 3_2_Im+_sd_1: 7 3_2_Im-_sd_0: 10 3_2_Im-_sd_1: 22 3_3_Re+_sd_0: 22 3_3_Re+_sd_1: 8 3_3_Re-_sd_0: 7 3_3_Re-_sd_1: 10 3_3_Im+_sd_0: 22 3_3_Im+_sd_1: 8 3_3_Im-_sd_0: 20 3_3_Im-_sd_1: 7
Custom basis-size reduction example¶
Now that we see that 2_2_Re+_sd_0 has a basis size of 64, let us try to restrict it to various degrees to see the effect on the $(2,2)$ mode. Custom presets are an advanced, provisional interface: adding an entry to MODEL_PRESETS affects the current Python process only, and a distinctive name should be used to avoid replacing a packaged preset.
Reducing basis sizes introduces additional approximation error whose magnitude is not known in advance. Before using a custom preset for scientific analysis, compare it against the full NRSur7dq4v2 model across the relevant parameter region and waveform modes. Waveform mismatches are recommended for assessing this additional truncation error. This comparison does not independently measure the full model's error relative to numerical relativity. The residual comparison below is illustrative and tests only one parameter-space point.
In [17]:
Copied!
from gwsurrogate.new._model_presets import MODEL_PRESETS
dp_name = '2_2_Re+_sd_0'
basis_sizes = [32,16,8,4,2]
custom_preset_name = 'TutorialCustomBasisSizes'
# We use the same binary black hole parameters as before
q = 4
chiA = [-0.2, 0.4, 0.1]
chiB = [-0.5, 0.2, -0.4]
dt = 0.1
f_low = 0
t, h, _ = sur(q, chiA, chiB, dt=dt, f_low=f_low)
h22re_waveforms = {} # We will store the real part of the (2,2) mode waveforms in this dictionary
h22re_waveforms[64] = np.real(h[(2,2)]) # Append the full surrogate waveform with basis size 64 to the dictionary
for size in basis_sizes:
del sur # Delete the surrogate to free up memory before creating a new one with smaller basis size
custom_basis_sizes = {dp_name: size}
MODEL_PRESETS['NRSur7dq4v2'][custom_preset_name] = custom_basis_sizes
sur = gwsurrogate.LoadSurrogate('NRSur7dq4v2', model_preset=custom_preset_name)
_, h, _ = sur(q, chiA, chiB, dt=dt, f_low=f_low)
h22re_waveforms[size] = np.real(h[(2,2)])
from gwsurrogate.new._model_presets import MODEL_PRESETS
dp_name = '2_2_Re+_sd_0'
basis_sizes = [32,16,8,4,2]
custom_preset_name = 'TutorialCustomBasisSizes'
# We use the same binary black hole parameters as before
q = 4
chiA = [-0.2, 0.4, 0.1]
chiB = [-0.5, 0.2, -0.4]
dt = 0.1
f_low = 0
t, h, _ = sur(q, chiA, chiB, dt=dt, f_low=f_low)
h22re_waveforms = {} # We will store the real part of the (2,2) mode waveforms in this dictionary
h22re_waveforms[64] = np.real(h[(2,2)]) # Append the full surrogate waveform with basis size 64 to the dictionary
for size in basis_sizes:
del sur # Delete the surrogate to free up memory before creating a new one with smaller basis size
custom_basis_sizes = {dp_name: size}
MODEL_PRESETS['NRSur7dq4v2'][custom_preset_name] = custom_basis_sizes
sur = gwsurrogate.LoadSurrogate('NRSur7dq4v2', model_preset=custom_preset_name)
_, h, _ = sur(q, chiA, chiB, dt=dt, f_low=f_low)
h22re_waveforms[size] = np.real(h[(2,2)])
Using model preset: TutorialCustomBasisSizes Loaded NRSur7dq4v2 model Using model preset: TutorialCustomBasisSizes Loaded NRSur7dq4v2 model Using model preset: TutorialCustomBasisSizes Loaded NRSur7dq4v2 model Using model preset: TutorialCustomBasisSizes Loaded NRSur7dq4v2 model Using model preset: TutorialCustomBasisSizes Loaded NRSur7dq4v2 model
In [18]:
Copied!
fig_basred, ax_basred = P.subplots(figsize=(10,6))
ax_basred.plot(t, h22re_waveforms[64], label='Basis size: 64 (Full Basis v2)')
for size in basis_sizes:
ax_basred.plot(t, h22re_waveforms[size], label=f'Basis size: {size}')
ax_basred.set_xlabel(r't [M]', fontsize=18)
ax_basred.set_ylabel(r'Re[$h_{22}$]', fontsize=18)
ax_basred.set_title('Effect of Basis Size Reduction on $h_{22}$', fontsize=16)
ax_basred.legend(fontsize=14)
P.show()
fig_basred, ax_basred = P.subplots(figsize=(10,6))
ax_basred.plot(t, h22re_waveforms[64], label='Basis size: 64 (Full Basis v2)')
for size in basis_sizes:
ax_basred.plot(t, h22re_waveforms[size], label=f'Basis size: {size}')
ax_basred.set_xlabel(r't [M]', fontsize=18)
ax_basred.set_ylabel(r'Re[$h_{22}$]', fontsize=18)
ax_basred.set_title('Effect of Basis Size Reduction on $h_{22}$', fontsize=16)
ax_basred.legend(fontsize=14)
P.show()
The waveforms are visibly very similar. In order to tell the difference, let us look at the residuals when subtracted from the full NRSur7dq4v2 waveform
In [19]:
Copied!
fig_basred, ax_basred = P.subplots(figsize=(10,6))
for size in basis_sizes:
ax_basred.plot(t, h22re_waveforms[64]-h22re_waveforms[size], label=f'Basis size: {size}', zorder =size) # Plotting the residuals with different zorder to ensure smaller basis sizes are plotted on top
ax_basred.set_xlabel(r't [M]', fontsize=18)
ax_basred.set_ylabel(r'Residual in Re[$h_{22}$]', fontsize=18)
ax_basred.set_title('Effect of Basis Size Reduction on $h_{22}$', fontsize=16)
ax_basred.legend(fontsize=14, loc='upper left')
P.show()
fig_basred, ax_basred = P.subplots(figsize=(10,6))
for size in basis_sizes:
ax_basred.plot(t, h22re_waveforms[64]-h22re_waveforms[size], label=f'Basis size: {size}', zorder =size) # Plotting the residuals with different zorder to ensure smaller basis sizes are plotted on top
ax_basred.set_xlabel(r't [M]', fontsize=18)
ax_basred.set_ylabel(r'Residual in Re[$h_{22}$]', fontsize=18)
ax_basred.set_title('Effect of Basis Size Reduction on $h_{22}$', fontsize=16)
ax_basred.legend(fontsize=14, loc='upper left')
P.show()
Now we see the effects of basis size reduction. It is clear that the smaller the basis size chosen, the largest the overall deviation from the full NRSur7dq4v2 waveform. Note that since we only changed the basis sizes of the sd_0 datapiece, the waveforms all agree exactly in the ringdown past $t=10M$.
Model presets¶
The package currently contains one model preset referred to as "Fast". This preset reduces basis sizes to trade a small amount of accuracy for faster evaluation, as assessed in arXiv:2609.07873. To use it, pass "Fast" as model_preset when loading NRSur7dq4v2.
In [20]:
Copied!
del sur # Delete the surrogate previously loaded so that we can reload it with a different basis size
# We use the same binary black hole parameters as before
q = 4
chiA = [-0.2, 0.4, 0.1]
chiB = [-0.5, 0.2, -0.4]
dt = 0.1
f_low = 0
# We first load the surrogate with the full basis as a reference
sur = gwsurrogate.LoadSurrogate('NRSur7dq4v2')
t, h_full, _ = sur(q, chiA, chiB, dt=dt, f_low=f_low)
del sur # Delete the surrogate to free up memory before creating a new one with smaller basis size
sur = gwsurrogate.LoadSurrogate('NRSur7dq4v2', model_preset="Fast")
t, h_fast, _ = sur(q, chiA, chiB, dt=dt, f_low=f_low)
del sur # Delete the surrogate previously loaded so that we can reload it with a different basis size
# We use the same binary black hole parameters as before
q = 4
chiA = [-0.2, 0.4, 0.1]
chiB = [-0.5, 0.2, -0.4]
dt = 0.1
f_low = 0
# We first load the surrogate with the full basis as a reference
sur = gwsurrogate.LoadSurrogate('NRSur7dq4v2')
t, h_full, _ = sur(q, chiA, chiB, dt=dt, f_low=f_low)
del sur # Delete the surrogate to free up memory before creating a new one with smaller basis size
sur = gwsurrogate.LoadSurrogate('NRSur7dq4v2', model_preset="Fast")
t, h_fast, _ = sur(q, chiA, chiB, dt=dt, f_low=f_low)
Loaded NRSur7dq4v2 model Using model preset: Fast Loaded NRSur7dq4v2 model
In [21]:
Copied!
fig_fast, ax_fast = P.subplots(4, 2, figsize=(15,12), tight_layout=True)
modes_to_compare = [(2,2), (3,2), (4,4), (5,-2)] # List of modes to compare
for i, (ell, m) in enumerate(modes_to_compare):
# Use the left column to compare the waveforms from the full basis and the fast basis for different modes
ax_fast[i,0].plot(t, h_full[(ell,m)].real, label='Full Basis v2')
ax_fast[i,0].plot(t, h_fast[(ell,m)].real, label='Fast Basis')
ax_fast[i,0].set_xlabel(r't [M]', fontsize=18)
ax_fast[i,0].set_ylabel(r'Re[$h_{%d,%d}$]' % (ell,m), fontsize=18)
ax_fast[i,0].set_title('Effect of Fast Basis on $h_{%d,%d}$' % (ell,m), fontsize=16)
ax_fast[i,0].legend(fontsize=14)
# Use the right column to plot the residuals
ax_fast[i,1].plot(t, h_full[(ell,m)].real - h_fast[(ell,m)].real)
ax_fast[i,1].set_xlabel(r't [M]', fontsize=18)
ax_fast[i,1].set_ylabel(r'Re[$\Delta h_{%d,%d}$]' % (ell,m), fontsize=18)
ax_fast[i,1].set_title('Residuals for $h_{%d,%d}$' % (ell,m), fontsize=16)
P.show()
fig_fast, ax_fast = P.subplots(4, 2, figsize=(15,12), tight_layout=True)
modes_to_compare = [(2,2), (3,2), (4,4), (5,-2)] # List of modes to compare
for i, (ell, m) in enumerate(modes_to_compare):
# Use the left column to compare the waveforms from the full basis and the fast basis for different modes
ax_fast[i,0].plot(t, h_full[(ell,m)].real, label='Full Basis v2')
ax_fast[i,0].plot(t, h_fast[(ell,m)].real, label='Fast Basis')
ax_fast[i,0].set_xlabel(r't [M]', fontsize=18)
ax_fast[i,0].set_ylabel(r'Re[$h_{%d,%d}$]' % (ell,m), fontsize=18)
ax_fast[i,0].set_title('Effect of Fast Basis on $h_{%d,%d}$' % (ell,m), fontsize=16)
ax_fast[i,0].legend(fontsize=14)
# Use the right column to plot the residuals
ax_fast[i,1].plot(t, h_full[(ell,m)].real - h_fast[(ell,m)].real)
ax_fast[i,1].set_xlabel(r't [M]', fontsize=18)
ax_fast[i,1].set_ylabel(r'Re[$\Delta h_{%d,%d}$]' % (ell,m), fontsize=18)
ax_fast[i,1].set_title('Residuals for $h_{%d,%d}$' % (ell,m), fontsize=16)
P.show()
