Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Feature - calibrated argument #53

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Feature - calibrated argument #53

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
57fb9a9
added calibration parameter to eval surrogate call and corresponding …
ssingh-sgr Sep 6, 2024
de9a8ee
moved the surrogate file
ssingh-sgr Sep 6, 2024
14db080
removing the duplicate surrogate.py file
ssingh-sgr Sep 10, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
48 changes: 23 additions & 25 deletions gwsurrogate/surrogate.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -140,7 +140,7 @@ def __init__(self, path, deg=3, subdir='', closeQ=True):
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
def __call__(self, q, M=None, dist=None, phi_ref=None,\
f_low=None, times=None, units='dimensionless',\
singlemode_call=True):
calibrated=True, singlemode_call=True):
"""Return single mode surrogate evaluation for...

Input
Expand All @@ -155,7 +155,7 @@ def __call__(self, q, M=None, dist=None, phi_ref=None,\
times --- array of times at which surrogate is to be evaluated
units --- units (mks or dimensionless) of input array of time samples
singlemode_call --- Never set this by hand, will be False if this routine is called by the multimode evaluator

calibrated --- If set to false, the calibration coefficients alpha and beta are set to 1
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

calibration is only applicable for the models EMRISur1dq1e4 and BHPTNRSur1dq1e4. Can you please add a short remark about this in the docstring?


More information
================
Expand Down Expand Up @@ -183,7 +183,7 @@ def __call__(self, q, M=None, dist=None, phi_ref=None,\

if (M is not None) and (dist is not None) and (times is not None) and (units != 'mks'):
raise ValueError('passing values of M, dist, and times suggest mks units should be used!')

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

At this part of the code there's a bunch of checks to make sure the user passes something reasonable.

I think we should check that if the model is not one of EMRISur1dq1e4 or BHPTNRSur1dq1e4 and calibration=False an error is raised. Because this would suggest the user is trying to override the calibration default for a model that doesn't support this parameter in the first place.

It might also be safer to assign calibration=None if the model is not one of EMRISur1dq1e4 or BHPTNRSur1dq1e4. This is just another safeguard that the code doesn't somehow depend on the value of calibration when it shouldn't

if self.surrogate_mode_type == 'coorb_waveform_basis' and singlemode_call:
msg = 'directly calling a coorb_waveform_basis surrogate will return a waveform in the co-orbital frame.\n'
msg += 'Please use EvaluateSurrogate to evaluate co-orbital frame surrogates.\n'
Expand All @@ -196,18 +196,21 @@ def __call__(self, q, M=None, dist=None, phi_ref=None,\
# at a given value of parameter space. These calibration parameters will scale the
# time and amplitude. Currently, these models use calibration:
# BHPTNRSur1dq1e4, EMRISur1dq1e4
if(self.surrogateID == 'BHPTNRSur1dq1e4'):
if not calibrated:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

perhaps nest this entire if-statement under if(self.surrogateID in ['EMRISur1dq1e4','BHPTNRSur1dq1e4']):. So it could be something like

if(self.surrogateID in ['EMRISur1dq1e4','BHPTNRSur1dq1e4']):
  # CODE YOU ALREADY HAVE
else: 
   alpha_nr_calibration, beta_nr_calibration = None, None

Which is again to be extra careful about not defining calibration parameter values for models that don't use them

alpha_nr_calibration, beta_nr_calibration = 1,1
elif(self.surrogateID == 'BHPTNRSur1dq1e4'):
alpha_nr_calibration, beta_nr_calibration = self.compute_BHPT_calibration_params(q)

if(self.surrogateID == 'EMRISur1dq1e4'):
elif(self.surrogateID == 'EMRISur1dq1e4'):
# see Eq 4 of https://arxiv.org/abs/1910.10473
# if alpha = 1 we have the "raw" waveform computed from
# point-particle black hole perturbation theory
# alpha rescales all modes in the same way in time and amplitude
nu = q/(1.+q)**2.
alpha_emri = 1.0-1.352854*nu-1.223006*nu*nu+8.601968*nu*nu*nu-46.74562*nu*nu*nu*nu
alpha_nr_calibration = 1.0-1.352854*nu-1.223006*nu*nu+8.601968*nu*nu*nu-46.74562*nu*nu*nu*nu
beta_nr_calibration = alpha_nr_calibration
else:
alpha_emri = None
alpha_nr_calibration = None
beta_nr_calibration = None


### if (M,distance) provided, a physical mode in mks units is returned ###
Expand All @@ -220,9 +223,7 @@ def __call__(self, q, M=None, dist=None, phi_ref=None,\
t_scale = 1.0

# any model-specific amplitude scalings should go here
if(self.surrogateID == 'EMRISur1dq1e4'):
amp0 = alpha_emri * amp0
if(self.surrogateID == 'BHPTNRSur1dq1e4'):
if(self.surrogateID in ['EMRISur1dq1e4','BHPTNRSur1dq1e4']):
amp0 = alpha_nr_calibration * amp0


Expand All @@ -231,10 +232,8 @@ def __call__(self, q, M=None, dist=None, phi_ref=None,\
t = times
else:
t = self.time() # shallow copy of self.times
if self.surrogateID == 'EMRISur1dq1e4':
t = t * alpha_emri # this will deep copy, preserving data in self.times
elif self.surrogateID == 'BHPTNRSur1dq1e4':
t = t * beta_nr_calibration
if self.surrogateID in ['EMRISur1dq1e4','BHPTNRSur1dq1e4']:
t = t * beta_nr_calibration # this will deep copy, preserving data in self.times

### if input times are dimensionless, convert to MKS if a physical surrogate is requested ###
if units == 'dimensionless':
Expand All @@ -245,9 +244,7 @@ def __call__(self, q, M=None, dist=None, phi_ref=None,\
times = times / t_scale

# input time grid needs to be mapped back to "raw EMRI" time grid before evaluation
if times is not None and self.surrogateID == 'EMRISur1dq1e4':
times = times / alpha_emri
if times is not None and self.surrogateID == 'BHPTNRSur1dq1e4':
if times is not None and self.surrogateID in ['EMRISur1dq1e4','BHPTNRSur1dq1e4']:
times = times / beta_nr_calibration

# convert from input to internal surrogate parameter values, and check within training region #
Expand Down Expand Up @@ -1059,7 +1056,7 @@ def __init__(self, path, deg=3, ell_m=None, excluded='DEFAULT', use_orbital_plan
def __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):
ell=None, m=None, mode_sum=True,fake_neg_modes=True, calibrated=True):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same comments as above about the docstring discussing calibrated

"""Return surrogate evaluation for...

INPUT
Expand All @@ -1078,6 +1075,7 @@ def __call__(self, q, M=None, dist=None, theta=None,phi=None,
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
calibrated --- if false, set the calibration coefficients alpha and beta to 1

NOTE: if only requesting one mode, this should be ell=[2],m=[2]

Expand Down Expand Up @@ -1148,9 +1146,9 @@ def __call__(self, q, M=None, dist=None, theta=None,phi=None,

# if model is BHPTNRSur1dq1e4 and mode not 22, hp/hc are in the coorbital frame
if is_modeled:
t_mode, hp_mode, hc_mode = self.evaluate_single_mode(q,M,dist,f_low,times,units,ell,m)
t_mode, hp_mode, hc_mode = self.evaluate_single_mode(q,M,dist,f_low,times,units,ell,m,calibrated)
else: # then we must have neg_modeled=True and fake_neg_modes=True
t_mode, hp_mode, hc_mode = self.evaluate_single_mode_by_symmetry(q,M,dist,f_low,times,units,ell,m)
t_mode, hp_mode, hc_mode = self.evaluate_single_mode_by_symmetry(q,M,dist,f_low,times,units,ell,m,calibrated)

# BHPTNRSur1dq1e4 models the real and imag parts of the higher order modes
# in the coorbital frame. We have to make sure we apply a frame transformation
Expand Down Expand Up @@ -1240,23 +1238,23 @@ def evaluate_on_sphere(self,ell,m,theta,phi,hp_mode,hc_mode):
return hp_mode, hc_mode

#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
def evaluate_single_mode(self,q, M, dist, f_low, times, units,ell,m):
def evaluate_single_mode(self,q, M, dist, f_low, times, units,ell,m, calibrated):
""" light wrapper around single mode evaluator"""

t_mode, hp_mode, hc_mode = self.single_mode_dict[(ell,m)](q, M, dist, None, f_low, times,
units,singlemode_call=False)
units, calibrated, singlemode_call=False)

return t_mode, hp_mode, hc_mode


#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# TODO: this routine will evaluate m>0 to get the m<0 case!!! It would be much better
# to simply use m>0 modes that were already evaluated (could be ~2x faster)
def evaluate_single_mode_by_symmetry(self,q, M, dist, f_low, times, units,ell,m):
def evaluate_single_mode_by_symmetry(self,q, M, dist, f_low, times, units,ell,m, calibrated):
""" evaluate m<0 mode from m>0 mode and relationship between these"""

if m<0:
t_mode, hp_mode, hc_mode = self.evaluate_single_mode(q, M, dist, f_low, times, units,ell,-m)
t_mode, hp_mode, hc_mode = self.evaluate_single_mode(q, M, dist, f_low, times, units,ell,-m, calibrated)
hp_mode, hc_mode = self._generate_minus_m_mode(hp_mode,hc_mode,ell,-m)
else:
raise ValueError('m must be negative.')
Expand Down
Loading