Utilities needed for aligned-spin NR-hybrid surrogate models. More...

Prototypes
int	ReadHDF5DoubleDataset (REAL8 param, LALH5File sub, const char *name)
	Reads a REAL8 value from a H5 file/group. More...

int	ReadHDF5IntDataset (int param, LALH5File sub, const char *name)
	Reads an INT8 value from a H5 file/group. More...

static int	NRHybSur_LoadDataPiece (DataPiece *data_piece, LALH5File file, const char *sub_grp_name)
	Loads a single waveform data piece from a H5 file. More...

static int	NRHybSur_LoadSingleModeData (ModeDataPieces *mode_data_pieces, LALH5File file, const int mode_idx, const gsl_matrix_long *mode_list)
	Loads all data pieces of a single waveform mode. More...

int	NRHybSur_Init (NRHybSurData NR_hybsur_data, LALH5File file)
	Initialize a NRHybSurData structure from an open H5 file. More...

static REAL8	kernel (const gsl_vector x1, const gsl_vector x2, const GPRHyperParams hyperparams, gsl_vector dummy_worker)
	The GPR Kernel evaluated at a single pair of points. More...

static REAL8	gp_predict (const gsl_vector xst, const GPRHyperParams hyperparams, const gsl_matrix x_train, gsl_vector dummy_worker)
	Evaluate a GPR fit. More...

REAL8	NRHybSur_eval_fit (const NRHybSurFitData fit_data, const gsl_vector fit_params, const gsl_matrix x_train, gsl_vector dummy_worker)
	Evaluate a NRHybSur fit. More...

static int	NRHybSur_eval_data_piece (gsl_vector *result, const gsl_vector fit_params, const DataPiece data_piece, const gsl_matrix x_train, gsl_vector *dummy_worker)
	Evaluate a single NRHybSur waveform data piece. More...

static gsl_vector *	spline_array_interp (const gsl_vector xout, const gsl_vector x, const gsl_vector *y)
	Do cubic spline interpolation using a gsl_interp_cspline. More...

static REAL8	compute_omega_22 (const gsl_vector times, const gsl_vector phi_22, const int idx)
	Computes omega_22 by forward differences. More...

static int	search_omega_22 (int omega_idx, const gsl_vector times, const gsl_vector *phi_22, const REAL8 omegaM_22_val)
	Finds closest index such that omegaM_22[index] = omegaM_22_val. More...

static int	NRHybSur_eval_phase_22_sparse (gsl_vector *phi_22_sparse, const REAL8 eta, const gsl_vector fit_params, gsl_vector dummy_dp, const gsl_matrix x_train, gsl_vector dummy_worker, const NRHybSurData NR_hybsur_data)
	Evaluate the phase of the (2,2) mode on the sprase surrogate domain. More...

static int	NRHybSur_upsample_phi_22 (gsl_vector phi_22_dense, gsl_vector times_dense, const REAL8 deltaTOverM, const REAL8 omegaM_22_min, const gsl_vector phi_22_sparse, const gsl_vector domain)
	Upsamples sparse \( \phi_{22} \) such that time step size is deltaTOverM, and initial frequency is roughly omegaM_22_min. More...

int	NRHybSur_eval_phase_22 (gsl_vector phi_22, gsl_vector output_times, const REAL8 eta, const gsl_vector fit_params, const REAL8 omegaM_22_min, const REAL8 deltaTOverM, const REAL8 phiRef, const REAL8 omegaM_22_ref, gsl_vector dummy_dp, const gsl_matrix x_train, gsl_vector dummy_worker, const NRHybSurData *NR_hybsur_data)
	Evaluate the phase of the (2,2) mode. More...

int	NRHybSur_eval_mode_data_pieces (EvaluatedDataPieces *this_mode_eval_dp, const UINT4 ell, const UINT4 m, const ModeDataPieces data_pieces, const gsl_vector output_times, const gsl_vector fit_params, gsl_vector dummy_dp, const gsl_matrix x_train, gsl_vector dummy_worker, const NRHybSurData NR_hybsur_data)
	Evaluate waveform data pieces of a single mode. More...

void	NRHybSur_DestroyEvaluatedDataPieces (gsl_vector phi_22, EvaluatedDataPieces *evaluated_mode_dps, const UINT4 num_modes_incl)
	Destroy phi_22 and an EvaluatedDataPieces structure. More...

int	NRHybSur_sanity_check_sample_rate (REAL8 deltaT, REAL8 m1, REAL8 m2, REAL8 chi1z, REAL8 chi2z, UINT4 max_ell)
	Sanity check (warning only, not error) that the sample rate is high enough to capture the ringdown frequencies, by ensuring Nyquist frequency is greater than the QNM frequency of the (max_ell,max_ell,0) mode, where max_ell is the maximum ell index among the included modes. More...

void	NRHybSur_set_default_modes (LALValue ModeArray, const NRHybSurData NR_hybsur_data)
	Activates all modes of an NRHybSur model. More...

int	NRHybSur_check_mode_array (UINT4 num_modes_incl, UINT4 max_ell, LALValue ModeArray, const NRHybSurData NR_hybsur_data)
	Sanity checks on ModeArray. More...

Detailed Description

Utilities needed for aligned-spin NR-hybrid surrogate models.

Author: Vijay Varma

Called from: LALSimIMRNRHybSur3dq8.h

Definition in file LALSimNRHybSurUtilities.c.

Go to the source code of this file.

Function Documentation

◆ ReadHDF5DoubleDataset()

int ReadHDF5DoubleDataset	(	REAL8 *	param,
		LALH5File *	sub,
		const char *	name
	)

Reads a REAL8 value from a H5 file/group.

Parameters

param	Output: REAL8 value from H5 file/group.
sub	H5 file or group.
name	Name of REAL8 dataset within file/group.

Definition at line 68 of file LALSimNRHybSurUtilities.c.

◆ ReadHDF5IntDataset()

int ReadHDF5IntDataset	(	int *	param,
		LALH5File *	sub,
		const char *	name
	)

Reads an INT8 value from a H5 file/group.

Parameters

param	Output: int value from H5 file/group.
sub	H5 file or group.
name	Name of int dataset within file/group.

Definition at line 92 of file LALSimNRHybSurUtilities.c.

◆ NRHybSur_LoadDataPiece()

static int NRHybSur_LoadDataPiece	(	DataPiece **	data_piece,
		LALH5File *	file,
		const char *	sub_grp_name
	)

static

Loads a single waveform data piece from a H5 file.

Parameters

data_piece	Output: Waveform data piece. *data_piece should be NULL. Space will be allocated.
file	Opened HDF5 file.
sub_grp_name	H5 group name.

Definition at line 122 of file LALSimNRHybSurUtilities.c.

◆ NRHybSur_LoadSingleModeData()

static int NRHybSur_LoadSingleModeData	(	ModeDataPieces **	mode_data_pieces,
		LALH5File *	file,
		const int	mode_idx,
		const gsl_matrix_long *	mode_list
	)

static

Loads all data pieces of a single waveform mode.

Parameters

mode_data_pieces	Output: Waveform data pieces of a given mode. Space will be allocated to **mode_data_pieces.
file	Opened HDF5 file.
mode_idx	Index corresponding to the mode.
mode_list	List of all modes.

Definition at line 249 of file LALSimNRHybSurUtilities.c.

◆ NRHybSur_Init()

int NRHybSur_Init	(	NRHybSurData *	NR_hybsur_data,
		LALH5File *	file
	)

Initialize a NRHybSurData structure from an open H5 file.

This will typically only be called once. For example from NRHybSur3dq8_Init_LALDATA.

The HDF5 file should have the following structure:

Top level:
'GPR_X_train': h5 dataset, shape=(N,dim) where N=size of training dataset.
'domain': h5 dataset, size = number of time samples, can be sparse.
'TaylorT3_t_ref': h5 dataset, INT8, used in TaylorT3 piece.
'mode_list': h5 dataset, INT8, shape=(M,2), where M=num. of modes. Each element is (ell,m) of that mode. Only m>=0 modes are modeled. First mode should always be (2,2).
'phaseAlignIdx: h5 dataset, INT8. Needed in TaylorT3 piece.
'l2_m0': h5 group, surrogate data needed for this mode.
'l2_m1':
'l2_m2':
- Subgroups for each mode:
  - 'l2_m2':
    'amp': h5 group, surrogate data needed for amplitude data piece.
    'phase' : h5 group, surrogate data needed for phase data piece.
  - 'l2_m1' and all other modes:
    're': h5 group, data needed for real part of coorbital frame waveform data piece.
    'im': h5 group, data needed for imag part of coorbital frame waveform data piece.
  - Subgroups for each data piece:
    'ei_basis': h5 dataset, shape=(n_nodes, L), where L=len(domain).
    'n_nodes': h5 dataset, INT8. Num. of empirical time nodes for this data piece.
    'node_num_0': h5 group, fit data for this node.
    'node_num_1':
    - Subgroups for each node:
      'alpha': h5 dataset, size=N, where N=size of training dataset.
      'constant_value': h5 dataset, REAL8, constant value in kernel.
      'data_mean': h5 dataset, REAL8. Mean of fit data.
      'data_std': h5 dataset, REAL8. Standard deviation of fit data.
      'length_scale': h5 dataset, size=dom, where dim=dimensions of the model. Length scale parameters in the kernel.
      'lin_coef': h5 dataset, size=dim. Coefs of of linear fit.
      'lin_intercept': h5 dataset, REAL8. Intercept of linear fit.
      'L': h5 dataset, shape=(N,N), where N=size of training dataset. Not used right now, needed to evaluate error estimate.
      'noise_level': h5 dataset, REAL8. Noise parameter, not needed for evaluating fit.
      'y_train_mean': h5 dataset, REAL8. Mean value before doing GPR fit, usually zero. We already remove the mean and store it in data_mean.

Parameters

NR_hybsur_data	Output: Struct to save surrogate data.
file	Opened HDF5 file.

Definition at line 403 of file LALSimNRHybSurUtilities.c.

◆ kernel()

static REAL8 kernel	(	const gsl_vector *	x1,
		const gsl_vector *	x2,
		const GPRHyperParams *	hyperparams,
		gsl_vector *	dummy_worker
	)

static

The GPR Kernel evaluated at a single pair of points.

We follow sklearn closely. We use the kernel given in Eq.(S3) of arxiv:1809.09125 :

\[ K(x, x') = \sigma_k^2 exp(-1/2 \sum_i^D (x^{i} - x'^{i})^2/\sigma_i^2) \]

where D is the dimension of the model.

Note that Eq.(S3) also includes the WhiteKernel, which has a noise parameter, but we don't need that here since we only need to evaluate \( K_{x* x} \) when evaluating the fits (See the mean value in Eq.(S2) of same paper). The other term we need is alpha = \( K_{x x}^{-1} {\bf f}\), which involves the WhiteKernel, but is precomputed offline. alpha is a vector of size N, where N is the number of cases in the training data set.

Parameters

x1	Parameter space point 1.
x2	Parameter space point 2.
hyperparams	GPR hyperparameters.
dummy_worker	Dummy worker array for computations.

Definition at line 520 of file LALSimNRHybSurUtilities.c.

◆ gp_predict()

static REAL8 gp_predict	(	const gsl_vector *	xst,
		const GPRHyperParams *	hyperparams,
		const gsl_matrix *	x_train,
		gsl_vector *	dummy_worker
	)

static

Evaluate a GPR fit.

See Eq.(S2) of arxiv:1809.09125.

Parameters

xst	Point \( x_* \) where fit will be evaluated.
hyperparams	Hyperparams for the GPR kernel.
x_train	Training set points.
dummy_worker	Dummy worker array for computations.

Definition at line 548 of file LALSimNRHybSurUtilities.c.

◆ NRHybSur_eval_fit()

REAL8 NRHybSur_eval_fit	(	const NRHybSurFitData *	fit_data,
		const gsl_vector *	fit_params,
		const gsl_matrix *	x_train,
		gsl_vector *	dummy_worker
	)

Evaluate a NRHybSur fit.

Parameters

fit_data	Data for fit.
fit_params	Parameter space point to evaluate the fit at. size=D, the dimension of the model.
x_train	Training set points.
dummy_worker	Dummy worker array for computations.

Definition at line 581 of file LALSimNRHybSurUtilities.c.

◆ NRHybSur_eval_data_piece()

static int NRHybSur_eval_data_piece	(	gsl_vector **	result,
		const gsl_vector *	fit_params,
		const DataPiece *	data_piece,
		const gsl_matrix *	x_train,
		gsl_vector *	dummy_worker
	)

static

Evaluate a single NRHybSur waveform data piece.

Parameters

result	Output: Should have already been assigned space.
fit_params	Parameter space point to evaluate the fit at. size=D, the dimension of the model.
data_piece	The waveform data piece to evaluate
x_train	Training set points.
dummy_worker	Dummy worker array for computations.

Definition at line 611 of file LALSimNRHybSurUtilities.c.

◆ spline_array_interp()

static gsl_vector* spline_array_interp	(	const gsl_vector *	xout,
		const gsl_vector *	x,
		const gsl_vector *	y
	)

static

Do cubic spline interpolation using a gsl_interp_cspline.

Results differ slightly from scipy.interpolate.InterpolatedUnivariateSpline due to different boundary conditions.

gwsurrogate by default will use gsl_interp_cspline, and results should agree exactly. However, gwsurrogate also has the option to use scipy interpolation, in which case the waveform values at the start and end can be slightly different. gwsurrogate will only use the scipy interpolation if for some reasong the gsl build fails.

Parameters

xout	The vector of points onto which we want to interpolate.
x	The x values of the data to interpolate.
y	The y values of the data to interpolate.

Definition at line 650 of file LALSimNRHybSurUtilities.c.

◆ compute_omega_22()

static REAL8 compute_omega_22	(	const gsl_vector *	times,
		const gsl_vector *	phi_22,
		const int	idx
	)

static

Computes omega_22 by forward differences.

Parameters

times	Time array.
phi_22	Phase of (2,2) mode.
idx	Index to compute omega_22 at.

Definition at line 674 of file LALSimNRHybSurUtilities.c.

◆ search_omega_22()

static int search_omega_22	(	int *	omega_idx,
		const gsl_vector *	times,
		const gsl_vector *	phi_22,
		const REAL8	omegaM_22_val
	)

static

Finds closest index such that omegaM_22[index] = omegaM_22_val.

Parameters

omega_idx	Output: closest index such that omegaM_22[index] = omegaM_22_val.
times	Time array.
phi_22	Phase of (2,2) mode.
omegaM_22_val	Desired frequency of (2,2) mode.

Definition at line 694 of file LALSimNRHybSurUtilities.c.

◆ NRHybSur_eval_phase_22_sparse()

static int NRHybSur_eval_phase_22_sparse	(	gsl_vector **	phi_22_sparse,
		const REAL8	eta,
		const gsl_vector *	fit_params,
		gsl_vector *	dummy_dp,
		const gsl_matrix *	x_train,
		gsl_vector *	dummy_worker,
		const NRHybSurData *	NR_hybsur_data
	)

static

Evaluate the phase of the (2,2) mode on the sprase surrogate domain.

The surrogate actually models the phase residual \( \phi^{res}_{22} \) defined in Eq.(44) of arxiv:1812.07865. Here we first evaluate that and then add the 0PN TaylorT3 phase to get the (2,2) mode phase.

Parameters

phi_22_sparse	Output: (2,2) mode phase on sparse domain.
eta	Symmetric mass ratio.
fit_params	Parameter space point to evaluate the fit at. size=D, the dimension of the model.
dummy_dp	Dummy vector to store phase evaluation.
x_train	Training set points.
dummy_worker	Dummy worker array for computations.
NR_hybsur_data	Loaded surrogate data.

Definition at line 740 of file LALSimNRHybSurUtilities.c.

◆ NRHybSur_upsample_phi_22()

static int NRHybSur_upsample_phi_22	(	gsl_vector **	phi_22_dense,
		gsl_vector **	times_dense,
		const REAL8	deltaTOverM,
		const REAL8	omegaM_22_min,
		const gsl_vector *	phi_22_sparse,
		const gsl_vector *	domain
	)

static

Upsamples sparse \( \phi_{22} \) such that time step size is deltaTOverM, and initial frequency is roughly omegaM_22_min.

The initial frequency will be refined later on.

Parameters

phi_22_dense	Output: Densely sampled (2,2) mode phase.
times_dense	Output: Dense time array.
deltaTOverM	Time step in M.
omegaM_22_min	Start frequency of (2,2) mode in rad/M.
phi_22_sparse	(2,2) mode phase on sparse domain.
domain	Sparse surrogate domain.

Definition at line 794 of file LALSimNRHybSurUtilities.c.

◆ NRHybSur_eval_phase_22()

int NRHybSur_eval_phase_22	(	gsl_vector **	phi_22,
		gsl_vector **	output_times,
		const REAL8	eta,
		const gsl_vector *	fit_params,
		const REAL8	omegaM_22_min,
		const REAL8	deltaTOverM,
		const REAL8	phiRef,
		const REAL8	omegaM_22_ref,
		gsl_vector *	dummy_dp,
		const gsl_matrix *	x_train,
		gsl_vector *	dummy_worker,
		const NRHybSurData *	NR_hybsur_data
	)

Evaluate the phase of the (2,2) mode.

The surrogate actually models the phase residual \( \phi^{res}_{22} \) defined in Eq.(44) of arxiv:1812.07865. Here we first evaluate that and then add the 0PN TaylorT3 phase to get the (2,2) mode phase.

Sets the orbital phase to phiRef at the reference frequency omegaM_22_ref. The orbital phase is obtained as phi_22/2, so this leaves a pi ambiguity. But the surrogate data is already aligned such that the heavier BH is on the +ve x-axis at t=-1000M. See Sec.VI.A.4 of arxiv:1812.07865, the resolves the pi ambiguity. This means that the after the realignment, the orbital phase at reference frequency omegaM_22_ref is phiRef, or the heavier BH is at azimuthal angle = phiRef from the +ve x-axis.

Only uses data at (2,2) mode frequencies >= omegaM_22_min. This determines the start time. The start time, along with the step size deltaTOverM, is used to determine the output_times. Uses cubic spline interpolation to interpolate from the surrogate's time array to output_times.

Parameters

phi_22	Output: (2,2) mode phase.
output_times	Output: Time array.
eta	Symmetric mass ratio.
fit_params	Parameter space point to evaluate the fit at. size=D, the dimension of the model.
omegaM_22_min	Start frequency of (2,2) mode in rad/M.
deltaTOverM	Time step in M.
phiRef	Orbital phase at reference frequency.
omegaM_22_ref	Reference freq of (2,2) mode in rad/M.
dummy_dp	Dummy vector to store phase evaluation.
x_train	Training set points.
dummy_worker	Dummy worker array for computations.
NR_hybsur_data	Loaded surrogate data.

Definition at line 887 of file LALSimNRHybSurUtilities.c.

◆ NRHybSur_eval_mode_data_pieces()

int NRHybSur_eval_mode_data_pieces	(	EvaluatedDataPieces **	this_mode_eval_dp,
		const UINT4	ell,
		const UINT4	m,
		const ModeDataPieces *	data_pieces,
		const gsl_vector *	output_times,
		const gsl_vector *	fit_params,
		gsl_vector *	dummy_dp,
		const gsl_matrix *	x_train,
		gsl_vector *	dummy_worker,
		const NRHybSurData *	NR_hybsur_data
	)

Evaluate waveform data pieces of a single mode.

For (2,2) mode we model the amplitude and phase, but the phase is evaluated using NRHybSur_eval_phase_22, since it is required for all modes to transform from the coorbital frame to the inertial frame. For all other modes we evaluate the real and imaginary parts of the coorbital frame strain, defined in Eq.(39) of arxiv:1812.07865.

Only the required data pieces for a given mode will be evaluated.

Parameters

this_mode_eval_dp	Output: evaluated waveform data pieces of a single mode.
ell	\(\ell\) index of mode.
m	m index of mode.
data_pieces	Surrogate data pieces of this mode.
output_times	Time vector to evaluate at.
fit_params	Parameter space point to evaluate the fit at. size=D, the dimension of the model.
dummy_dp	Dummy vector to store phase evaluation.
x_train	Training set points.
dummy_worker	Dummy worker array for computations.
NR_hybsur_data	Loaded surrogate data.

Definition at line 978 of file LALSimNRHybSurUtilities.c.

◆ NRHybSur_DestroyEvaluatedDataPieces()

void NRHybSur_DestroyEvaluatedDataPieces	(	gsl_vector *	phi_22,
		EvaluatedDataPieces **	evaluated_mode_dps,
		const UINT4	num_modes_incl
	)

Destroy phi_22 and an EvaluatedDataPieces structure.

Free all associated memory.

Parameters

phi_22	\(\phi_{22}\) data piece.
evaluated_mode_dps	All other data pieces.
num_modes_incl	Number of models included.

Definition at line 1056 of file LALSimNRHybSurUtilities.c.

◆ NRHybSur_sanity_check_sample_rate()

int NRHybSur_sanity_check_sample_rate	(	REAL8	deltaT,
		REAL8	m1,
		REAL8	m2,
		REAL8	chi1z,
		REAL8	chi2z,
		UINT4	max_ell
	)

Sanity check (warning only, not error) that the sample rate is high enough to capture the ringdown frequencies, by ensuring Nyquist frequency is greater than the QNM frequency of the (max_ell,max_ell,0) mode, where max_ell is the maximum ell index among the included modes.

Parameters

deltaT	Sampling interval (s).
m1	Mass of Bh1 (kg).
m2	Mass of Bh2 (kg).
chi1z	Dimensionless spin of Bh1.
chi2z	Dimensionless spin of Bh2.
max_ell	Max ell index included.

Definition at line 1101 of file LALSimNRHybSurUtilities.c.

◆ NRHybSur_set_default_modes()

void NRHybSur_set_default_modes	(	LALValue *	ModeArray,
		const NRHybSurData *	NR_hybsur_data
	)

Activates all modes of an NRHybSur model.

For NRHybSur3dq8 that is \( \ell \leq 4, m \geq 0 \), and (5,5), but not (4,1) or (4,0).

Parameters

ModeArray	Output: Container for modes.
NR_hybsur_data	Loaded surrogate data.

Definition at line 1151 of file LALSimNRHybSurUtilities.c.

◆ NRHybSur_check_mode_array()

int NRHybSur_check_mode_array	(	UINT4 *	num_modes_incl,
		UINT4 *	max_ell,
		LALValue *	ModeArray,
		const NRHybSurData *	NR_hybsur_data
	)

Sanity checks on ModeArray.

Will raise an error if an unavailable mode is requested. Note that we only look for m>=0 modes in ModeArray, and will ignore m<0 modes even if present. The m<0 modes automatically get added when evaluting the waveform.

Parameters

num_modes_incl	Output: Number of modes to include.
max_ell	Output: Max ell index included.
ModeArray	Container for modes.
NR_hybsur_data	Loaded surrogate data.

Definition at line 1175 of file LALSimNRHybSurUtilities.c.

Prototypes

Detailed Description

Function Documentation

◆ ReadHDF5DoubleDataset()

◆ ReadHDF5IntDataset()

◆ NRHybSur_LoadDataPiece()

◆ NRHybSur_LoadSingleModeData()

◆ NRHybSur_Init()

◆ kernel()

◆ gp_predict()

◆ NRHybSur_eval_fit()

◆ NRHybSur_eval_data_piece()

◆ spline_array_interp()

◆ compute_omega_22()

◆ search_omega_22()

◆ NRHybSur_eval_phase_22_sparse()

◆ NRHybSur_upsample_phi_22()

◆ NRHybSur_eval_phase_22()

◆ NRHybSur_eval_mode_data_pieces()

◆ NRHybSur_DestroyEvaluatedDataPieces()

◆ NRHybSur_sanity_check_sample_rate()

◆ NRHybSur_set_default_modes()

◆ NRHybSur_check_mode_array()