Example usage of BHPTNRSur1dq1e4 surrogate model.¶

Notebook author: Tousif Islam (tislam@umassd.edu)

Contact: Tousif Islam (tislam@umassd.edu), Scott Field (sfield@umassd.edu), Vijay Varma (vijay.varma392@gmail.com)

In [1]:

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

import gwsurrogate as gws

import numpy as np import matplotlib.pyplot as plt %matplotlib inline import gwsurrogate as gws

cannot import LAL
setting __package__ to gwsurrogate.new so relative imports work
__name__ = gwsurrogate.new.spline_evaluation
__package__= gwsurrogate.new
setting __package__ to gwsurrogate.new so relative imports work
setting __package__ to gwsurrogate.new so relative imports work

/data/tislam/miniconda3/envs/gws/lib/python3.7/site-packages/gwtools/rotations.py:63: UserWarning: Could not import GWFrames, needed for rotations module
  _warnings.warn("Could not import GWFrames, needed for rotations module")
/data/tislam/miniconda3/envs/gws/lib/python3.7/site-packages/gwtools/__init__.py:11: UserWarning: Could not import rotations, decompositions, or fitfuncs. These are not needed by GWSurrogate.
  _warnings.warn("Could not import rotations, decompositions, or fitfuncs. These are not needed by GWSurrogate.")

Download surrogate data, this only needs to be done once¶

In [2]:

# This can take a few minutes
gws.catalog.pull('BHPTNRSur1dq1e4')

# This can take a few minutes gws.catalog.pull('BHPTNRSur1dq1e4')

BHPTNRSur1dq1e4.h5 file exits, moving to /data/tislam/git/gwsurrogate/gwsurrogate/surrogate_downloads/backup/2022Oct27_15h:23m:36s_BHPTNRSur1dq1e4.h5.

Out[2]:

'/data/tislam/git/gwsurrogate/gwsurrogate/surrogate_downloadsBHPTNRSur1dq1e4.h5'

Load the surrogate, this only needs to be done once at the start of a script¶

In [3]:

import sys
path_to_gws = '/data/tislam/git/gwsurrogate/gwsurrogate/'
sys.path.append(path_to_gws)
path_to_surrogate = path_to_gws+'surrogate_downloads/BHPTNRSur1dq1e4.h5'
model = gws.EvaluateSurrogate(path_to_surrogate, ell_m=[(2,2),(2,1),(3,1),(3,2),(3,3),(4,2),(4,3),(4,4),(5,3),(5,4),(5,5)]) 

import sys path_to_gws = '/data/tislam/git/gwsurrogate/gwsurrogate/' sys.path.append(path_to_gws) path_to_surrogate = path_to_gws+'surrogate_downloads/BHPTNRSur1dq1e4.h5' model = gws.EvaluateSurrogate(path_to_surrogate, ell_m=[(2,2),(2,1),(3,1),(3,2),(3,3),(4,2),(4,3),(4,4),(5,3),(5,4),(5,5)])

loading surrogate mode... l2_m2

>>> Found surrogate ID from file name: BHPTNRSur1dq1e4

>>> Warning: Guessing quadrature weights to be identical with 0.200000
Cannot load greedy points...OK
Special case: using spline for parametric model at each EI node
num_fits_amp = 7
num_fits_phase = 7
setting norm fitparams to None...

/data/tislam/git/gwsurrogate/gwsurrogate/surrogateIO.py:559: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  self.fitparams_amp = np.array(fitparams_amp)
/data/tislam/git/gwsurrogate/gwsurrogate/surrogateIO.py:560: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  self.fitparams_phase = np.array(fitparams_phase)

loading surrogate mode... l2_m1

>>> Found surrogate ID from file name: BHPTNRSur1dq1e4

>>> Warning: Guessing quadrature weights to be identical with 0.200000
Cannot load greedy points...OK
Special case: using spline for parametric model at each EI node
num_fits_re = 13
num_fits_im = 13
setting norm fitparams to None...

/data/tislam/git/gwsurrogate/gwsurrogate/surrogateIO.py:585: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  self.fitparams_re = np.array(fitparams_re)
/data/tislam/git/gwsurrogate/gwsurrogate/surrogateIO.py:586: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  self.fitparams_im = np.array(fitparams_im)

loading surrogate mode... l3_m1

>>> Found surrogate ID from file name: BHPTNRSur1dq1e4

>>> Warning: Guessing quadrature weights to be identical with 0.200000
Cannot load greedy points...OK
Special case: using spline for parametric model at each EI node
num_fits_re = 13
num_fits_im = 13
setting norm fitparams to None...
loading surrogate mode... l3_m2

>>> Found surrogate ID from file name: BHPTNRSur1dq1e4

>>> Warning: Guessing quadrature weights to be identical with 0.200000
Cannot load greedy points...OK
Special case: using spline for parametric model at each EI node
num_fits_re = 13
num_fits_im = 13
setting norm fitparams to None...
loading surrogate mode... l3_m3

>>> Found surrogate ID from file name: BHPTNRSur1dq1e4

>>> Warning: Guessing quadrature weights to be identical with 0.200000
Cannot load greedy points...OK
Special case: using spline for parametric model at each EI node
num_fits_re = 13
num_fits_im = 13
setting norm fitparams to None...
loading surrogate mode... l4_m2

>>> Found surrogate ID from file name: BHPTNRSur1dq1e4

>>> Warning: Guessing quadrature weights to be identical with 0.200000
Cannot load greedy points...OK
Special case: using spline for parametric model at each EI node
num_fits_re = 13
num_fits_im = 13
setting norm fitparams to None...
loading surrogate mode... l4_m3

>>> Found surrogate ID from file name: BHPTNRSur1dq1e4

>>> Warning: Guessing quadrature weights to be identical with 0.200000
Cannot load greedy points...OK
Special case: using spline for parametric model at each EI node
num_fits_re = 13
num_fits_im = 13
setting norm fitparams to None...
loading surrogate mode... l4_m4

>>> Found surrogate ID from file name: BHPTNRSur1dq1e4

>>> Warning: Guessing quadrature weights to be identical with 0.200000
Cannot load greedy points...OK
Special case: using spline for parametric model at each EI node
num_fits_re = 13
num_fits_im = 13
setting norm fitparams to None...
loading surrogate mode... l5_m3

>>> Found surrogate ID from file name: BHPTNRSur1dq1e4

>>> Warning: Guessing quadrature weights to be identical with 0.200000
Cannot load greedy points...OK
Special case: using spline for parametric model at each EI node
num_fits_re = 13
num_fits_im = 13
setting norm fitparams to None...
loading surrogate mode... l5_m4

>>> Found surrogate ID from file name: BHPTNRSur1dq1e4

>>> Warning: Guessing quadrature weights to be identical with 0.200000
Cannot load greedy points...OK
Special case: using spline for parametric model at each EI node
num_fits_re = 13
num_fits_im = 13
setting norm fitparams to None...
loading surrogate mode... l5_m5

>>> Found surrogate ID from file name: BHPTNRSur1dq1e4

>>> Warning: Guessing quadrature weights to be identical with 0.200000
Cannot load greedy points...OK
Special case: using spline for parametric model at each EI node
num_fits_re = 13
num_fits_im = 13
setting norm fitparams to None...
Surrogate interval [[0.39794001]
 [4.        ]]
Surrogate time grid [-30500.         -30499.8        -30499.6        ...    114.40000011
    114.60000011    114.80000011]
Surrogate parameterization map from q to log10(q)
  
  Surrogates with this parameterization expect its user intput 
  to be the mass ratio q. 
  
  The surrogate will map q to the internal surrogate's 
  parameterization which is log10(q)
  
  The surrogates training interval is quoted in log10(q).
  

Read the documentation¶

In [4]:

help(model)

help(model)

Help on EvaluateSurrogate in module gwsurrogate.surrogate object:

class EvaluateSurrogate(builtins.object)
 |  EvaluateSurrogate(path, deg=3, ell_m=None, excluded='DEFAULT', use_orbital_plane_symmetry=True)
 |  
 |  Evaluate multi-mode surrogates
 |  
 |  Methods defined here:
 |  
 |  __call__(self, q, M=None, dist=None, theta=None, phi=None, z_rot=None, f_low=None, times=None, units='dimensionless', ell=None, m=None, mode_sum=True, fake_neg_modes=True)
 |      Return surrogate evaluation for...
 |      
 |      INPUT
 |      =====
 |      q              --- binary parameter values EXCLUDING total mass M.
 |                         In 1D, mass ratio (dimensionless) must be supplied.
 |                         In nD, the surrogate's internal parameterization is assumed.
 |      M              --- total mass (solar masses)
 |      dist           --- distance to binary system (megaparsecs)
 |      theta/phi      --- evaluate hp and hc modes at this location on sphere
 |      z_rot          --- physical rotation about angular momentum (z-)axis (radians)
 |      flow           --- instantaneous initial frequency, will check if flow_surrogate < flow mode-by-mode
 |      times          --- array of times at which surrogate is to be evaluated
 |      units          --- units ('mks' or 'dimensionless') of input array of time samples
 |      ell            --- list or array of N ell modes to evaluate for (if none, all modes are returned)
 |      m              --- for each ell, supply a matching m value
 |      mode_sum       --- if true, all modes are summed, if false all modes are returned in an array
 |      fake_neg_modes --- if true, include m<0 modes deduced from m>0 mode. all m in [ell,m] input should be non-negative
 |      
 |      NOTE: if only requesting one mode, this should be ell=[2],m=[2]
 |      
 |       Note about Angles
 |       =================
 |       For circular orbits, the binary's orbital angular momentum is taken to
 |       be the z-axis. Theta and phi is location on the sphere relative to this
 |       coordinate system.
 |  
 |  __init__(self, path, deg=3, ell_m=None, excluded='DEFAULT', use_orbital_plane_symmetry=True)
 |      Loads a surrogate.
 |      
 |      path: the path to the surrogate
 |      deg: the degree of the splines representing the basis (default 3, cubic).
 |          Unless there is good reason to use deg !=3 one should not change this.
 |          Some surrogates (e.g. 4d2s) are validated with this in mind.
 |      ell_m: A list of (ell, m) modes to load, for example [(2,2),(3,3)].
 |          None (default) loads all modes.
 |      excluded: A list of (ell, m) modes to skip loading.
 |          The default ('DEFAULT') excludes any modes with an 'EXCLUDED' dataset.
 |          Use [] or None to load these modes as well.
 |      use_orbital_plane_symmetry: If set to true (i) CreateManyEvaluateSingleModeSurrogates
 |          will explicitly check that m<0 do not exist in the data file and (ii) m<0 modes
 |          are inferred from m>0 modes. If set to false no symmetry is assumed -- typical
 |          of precessing models. When False, fake_neg_modes must be false.
 |  
 |  add_l2m2_mode_if_not_in_modelist(self, list_of_modes)
 |      Adds the (2,2) mode if it is not in ell, m list. This is 
 |      required for some models which evaluates coorbital frame waveform
 |      for the higher modes (2,2) mode info is required here.
 |      
 |      This routine also ensures (2,2) is the front of the list
 |  
 |  all_model_modes(self, minus_m=False)
 |      from single mode keys deduce all available model modes.
 |      If minus_m=True, include (ell,-m) whenever (ell,m) is available .
 |  
 |  coorbital_to_inertial(self, coorb_re, coorb_im, m, orbital_phase)
 |      Takes the real and imaginary part of the waveform and
 |      combine them to obtain the coorbital frame waveform;
 |      then transform the wf into the inertial frame
 |  
 |  evaluate_on_sphere(self, ell, m, theta, phi, hp_mode, hc_mode)
 |      evaluate on the sphere
 |  
 |  evaluate_single_mode(self, q, M, dist, f_low, times, units, ell, m)
 |      light wrapper around single mode evaluator
 |  
 |  evaluate_single_mode_by_symmetry(self, q, M, dist, f_low, times, units, ell, m)
 |      evaluate m<0 mode from m>0 mode and relationship between these
 |  
 |  generate_mode_eval_list(self, ell=None, m=None, minus_m=False)
 |      generate list of (ell,m) modes to evaluate for.
 |      
 |      1) ell=m=None: use all available model modes
 |      2) ell=NUM, m=None: all modes up to ell_max = NUM. unmodelled modes set to zero
 |      3) list of [ell], [m] pairs: only use modes (ell,m). unmodelled modes set to zero
 |         ex: ell=[3,2] and m=[2,2] generates a (3,2) and (2,2) mode.
 |      
 |      These three options produce a list of (ell,m) modes.
 |      
 |      Set minus_m=True to generate m<0 modes from m>0 modes.
 |  
 |  h_sphere_builder(self, q, M=None, dist=None, ell=None, m=None)
 |      Returns a function for evaluations of h(t,theta,phi;q,M,d) which include
 |      all available modes.
 |      
 |      This new function h(t,theta,phi;q,M,d)
 |      can be evaluated for rotations about z-axis and at any set of
 |      points on the sphere. modes_to_evalute are also returned
 |  
 |  match_surrogate(self, t_ref, h_ref, q, M=None, dist=None, theta=None, t_ref_units='dimensionless', ell=None, m=None, fake_neg_modes=True, t_low_adj=0.0125, t_up_adj=0.0125, speed='slow')
 |      match discrete complex polarization (t_ref,h_ref) to surrogate waveform for
 |      given input values. Inputs have same meaning as those passed to __call__
 |      
 |      Minimization (i.e. match) over time shifts and z-axis rotations
 |  
 |  single_mode(self, mode)
 |      Returns a single-mode object for mode=(ell,m).
 |      This object stores information for the (ell,m)-mode surrogate
 |  
 |  sort_mode_list(self, mode_list)
 |      sort modes as (2,-2), (2,-1), ..., (2,2), (3,-3),(3,-2)...
 |  
 |  ----------------------------------------------------------------------
 |  Data descriptors defined here:
 |  
 |  __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]:

# Evaluate and plot the 2,2 mode. 
# By default, the modes are evaluated on the sphere, and negative modes are generated from 
# known relationships. So we need to set both options to false to get only the (2,2) mode.
modes, times, hp, hc = model(q=20, ell=[2], m=[2], mode_sum=False, fake_neg_modes=False)
print('You have evaluated the (%i,%i) mode'%(modes[0][0],modes[0][1]))

plt.plot(times, hp, label='Re[h]')
plt.plot(times,abs(hp+1j*hc),'r', label='Amp[h]')
plt.title('The (%i,%i) mode'%(modes[0][0],modes[0][1]))
plt.ylabel('$Re[h]$', fontsize=18)
plt.xlabel('t/M ')
plt.legend(fontsize=12)
plt.show()

# Evaluate and plot the 2,2 mode. # By default, the modes are evaluated on the sphere, and negative modes are generated from # known relationships. So we need to set both options to false to get only the (2,2) mode. modes, times, hp, hc = model(q=20, ell=[2], m=[2], mode_sum=False, fake_neg_modes=False) print('You have evaluated the (%i,%i) mode'%(modes[0][0],modes[0][1])) plt.plot(times, hp, label='Re[h]') plt.plot(times,abs(hp+1j*hc),'r', label='Amp[h]') plt.title('The (%i,%i) mode'%(modes[0][0],modes[0][1])) plt.ylabel('$Re[h]$', fontsize=18) plt.xlabel('t/M ') plt.legend(fontsize=12) plt.show()

You have evaluated the (2,2) mode

Evaluate waveform modes in physical units¶

In [6]:

mass=100 # Total masss in solar masses
dist_mpc=500 # distance in megaparsec

modes, times, hp, hc = model(q=20, M=mass, dist=dist_mpc, ell=[2], m=[2], mode_sum=False, fake_neg_modes=False)
print('You have evaluated the (%i,%i) mode in physical unit'%(modes[0][0],modes[0][1]))

plt.plot(times, hp, label='Re[h]')
plt.plot(times,abs(hp+1j*hc),'r', label='Amp[h]')
plt.ylabel('$h_{+}$', fontsize=18)
plt.xlabel('t [s]', fontsize=18)
plt.legend(fontsize=12)
plt.show()

mass=100 # Total masss in solar masses dist_mpc=500 # distance in megaparsec modes, times, hp, hc = model(q=20, M=mass, dist=dist_mpc, ell=[2], m=[2], mode_sum=False, fake_neg_modes=False) print('You have evaluated the (%i,%i) mode in physical unit'%(modes[0][0],modes[0][1])) plt.plot(times, hp, label='Re[h]') plt.plot(times,abs(hp+1j*hc),'r', label='Amp[h]') plt.ylabel('$h_{+}$', fontsize=18) plt.xlabel('t [s]', fontsize=18) plt.legend(fontsize=12) plt.show()

You have evaluated the (2,2) mode in physical unit

Evaluate waveform at a point on the sky¶

In [7]:

mass=100 # Total masss in solar masses
dist_mpc=500 # distance in megaparsec
inclination = np.pi/4 
phi_ref = np.pi/5

modes, times, hp, hc = model(q=20, M=mass, dist=dist_mpc, theta=inclination, phi=phi_ref, ell=[2], m=[2], mode_sum=False, fake_neg_modes=False)
print('You have evaluated the (%i,%i) mode in physical unit'%(modes[0][0],modes[0][1]))

plt.plot(times, hp, label='Re[h]')
plt.plot(times,abs(hp+1j*hc),'r', label='Amp[h]')
plt.title('The (%i,%i) mode'%(modes[0][0],modes[0][1]))
plt.xlabel('t [s]', fontsize=18)
plt.ylabel('$h_{+}$ $(\iota, \phi_{ref})$', fontsize=18)
plt.legend(fontsize=12)
plt.show()

mass=100 # Total masss in solar masses dist_mpc=500 # distance in megaparsec inclination = np.pi/4 phi_ref = np.pi/5 modes, times, hp, hc = model(q=20, M=mass, dist=dist_mpc, theta=inclination, phi=phi_ref, ell=[2], m=[2], mode_sum=False, fake_neg_modes=False) print('You have evaluated the (%i,%i) mode in physical unit'%(modes[0][0],modes[0][1])) plt.plot(times, hp, label='Re[h]') plt.plot(times,abs(hp+1j*hc),'r', label='Amp[h]') plt.title('The (%i,%i) mode'%(modes[0][0],modes[0][1])) plt.xlabel('t [s]', fontsize=18) plt.ylabel('$h_{+}$ $(\iota, \phi_{ref})$', fontsize=18) plt.legend(fontsize=12) plt.show()

You have evaluated the (2,2) mode in physical unit

In [8]:

mass=100 # Total masss in solar masses
dist_mpc=500 # distance in megaparsec
inclination = np.pi/4 
phi_ref = np.pi/5

# when mode_sum=True, (2,2) mode must be included in ell, m
times, hp, hc = model(q=20, M=mass, dist=dist_mpc, theta=inclination, phi=phi_ref, ell=[2,3], m=[2,3], mode_sum=True, fake_neg_modes=False)
print('You have evaluated the (%i,%i) mode in physical unit'%(modes[0][0],modes[0][1]))

plt.plot(times, hp, label='Re[h]')
plt.plot(times,abs(hp+1j*hc),'r', label='Amp[h]')
plt.title('The (%i,%i) mode'%(modes[0][0],modes[0][1]))
plt.xlabel('t [s]', fontsize=18)
plt.ylabel('$h_{+}$ $(\iota, \phi_{ref})$', fontsize=18)
plt.legend(fontsize=12)
plt.show()

mass=100 # Total masss in solar masses dist_mpc=500 # distance in megaparsec inclination = np.pi/4 phi_ref = np.pi/5 # when mode_sum=True, (2,2) mode must be included in ell, m times, hp, hc = model(q=20, M=mass, dist=dist_mpc, theta=inclination, phi=phi_ref, ell=[2,3], m=[2,3], mode_sum=True, fake_neg_modes=False) print('You have evaluated the (%i,%i) mode in physical unit'%(modes[0][0],modes[0][1])) plt.plot(times, hp, label='Re[h]') plt.plot(times,abs(hp+1j*hc),'r', label='Amp[h]') plt.title('The (%i,%i) mode'%(modes[0][0],modes[0][1])) plt.xlabel('t [s]', fontsize=18) plt.ylabel('$h_{+}$ $(\iota, \phi_{ref})$', fontsize=18) plt.legend(fontsize=12) plt.show()

You have evaluated the (2,2) mode in physical unit

In [ ]: