LALSimulation  5.4.0.1-8083555
eval_fits.py
Go to the documentation of this file.
1 """ Evaluates Numerical Relativity fits such as remnant BH mass, spin, etc, for
2  various models.
3 
4 Vijay Varma, 2019.
5 """
6 import numpy as np
7 import warnings
8 
9 from lal import MSUN_SI
10 
11 from .NRSur7dq4Remnant import NRSur7dq4Remnant
12 from .NRSur3dq8Remnant import NRSur3dq8Remnant
13 from . import quaternion_utils
14 
15 #=============================================================================
16 class FitAttributes(object):
17  """ Saves attributes of a particular fit.
18  """
19  #-------------------------------------------------------------------------
20  def __init__(self, fit_class, desc, refs, model_keywords):
21  """ Initializes a FitAttributes object. See usage examples below.
22 
23  - fit_class: Class to load for a given model. type: should be a
24  daughter class of nrfits.NRFits.
25  - desc: A brief description of the model, type: str.
26  - refs: Paper reference, type: str.
27  - model_keywords: List of strings, with one or more of
28  allowed_model_keywords.
29  """
30 
31  # Types of model keywords allowed (add more if required)
32  allowed_model_keywords = ['nonspinning', 'aligned_spin', 'precessing', \
33  'eccentric', 'tidal']
34 
35  for key in model_keywords:
36  if key not in allowed_model_keywords:
37  raise Exception("Invalid model_keyword: %s"%key)
38 
39  self.fit_classfit_class = fit_class
40  self.descdesc = desc
41  self.refsrefs = refs
42  self.model_keywordsmodel_keywords = model_keywords
43 
44 
45 #=============================================================================
46 
47 # Collection of all implemented fits
48 fits_collection = {}
49 
50 fits_collection['NRSur3dq8Remnant'] = FitAttributes( \
51  fit_class = NRSur3dq8Remnant(),
52  desc = 'Fits for remnant mass, spin and kick velocity for nonprecessing'
53  ' BBH systems, trained on mass ratios q <= 8 and dimensionless spin '
54  'magnitudes <= 0.81. This model was called surfinBH3dq8 in the paper.',
55  refs = 'arxiv:1809.09125',
56  model_keywords = ['aligned_spin'],
57  )
58 
59 fits_collection['NRSur7dq4Remnant'] = FitAttributes( \
60  fit_class = NRSur7dq4Remnant(),
61  desc = 'Fits for remnant mass, spin and kick velocity for generically'
62  ' precessing BBH systems, trained on mass ratios q <= 4.01 and '
63  'dimensionless spin magnitudes <= 0.81.',
64  refs = 'arxiv:1905.09300',
65  model_keywords = ['precessing'],
66  )
67 
68 
69 #=============================================================================
70 def truncate_output_to_physical_limits(val_dict, behavior):
71  """
72  Function to truncate fit values when they go beyond the following physical
73  limits:
74  - FinalMass: 1, FinalMass is expected as a fraction of total mass here.
75  - FinalSpin: magnitude = 1, the Kerr limit for dimensionless spins.
76  - RecoilKick: magnitude = 1, as this is in units of c.
77 
78  Inputs: \n
79  - val_dict: A dictionary of fits with one or more of the above keys. See
80  eval_nrfit.
81 
82  - behavior (str) can be one of the following:
83  * "TRUNCATE": Rescale output magnitude to maximum limit, with an
84  info message.
85  * "TRUNCATESILENT": Rescale output magnitude to maximum limit, without
86  info message
87  * "KEEP": Keep output magnitude, even if unphysical, but
88  still give the info message
89  * "IGNORE": Keep output magnitude, even if unphysical, without
90  info message
91  * "ERROR": Abort with an error if output magnitude is
92  unphysical.
93  """
94 
95  # Don't need to do anything
96  if behavior == "IGNORE":
97  return val_dict
98 
99  allowed_behaviors = ["ERROR","IGNORE","KEEP","TRUNCATE","TRUNCATESILENT"]
100  if behavior not in allowed_behaviors:
101  raise ValueError("Invalid physical_limit_violation_behavior=%s"%behavior
102  + ". Should be one of ["+ ", ".join(allowed_behaviors) + "].")
103 
104  physical_limits = {
105  "FinalMass": 1, # Remnant can't be heavier than the total mass
106  "FinalSpin": 1, # Kerr limit
107  "RecoilKick": 1, # Speed of light limit
108  }
109 
110  for fit_type in val_dict.keys():
111  limit = physical_limits[fit_type]
112  magnitude = np.linalg.norm(val_dict[fit_type])
113  if magnitude > limit:
114  # Error/Warning message
115  message = ""
116  if fit_type in ["FinalSpin", "RecoilKick"]:
117  message += "%s magnitude=%.7e"%(fit_type, magnitude)
118  else:
119  message += "Dimensionless %s=%.7e"%(fit_type, magnitude)
120  message += " is over max limit=%.2e."%(limit)
121 
122  if behavior == "ERROR": # raise error
123  raise RuntimeError(message + " Adapt " \
124  "extra_params_dict['physical_limit_violation_behavior'] to "\
125  "continue without an error.")
126 
127  elif behavior in ["TRUNCATE", "TRUNCATESILENT"]: # rescale magnitude
128  message += " Setting to max limit."
129  val_dict[fit_type] *= limit/magnitude
130 
131  # Only these cases get a warning
132  if behavior in ["KEEP", "TRUNCATE"]:
133  warnings.warn(message)
134 
135  return val_dict
136 
137 #=============================================================================
138 def check_extra_params_and_set_defaults(extra_params_dict):
139  """ Does some sanity checks on extra_params_dict.
140  If any of the default_keywords are not specified, sets them to the
141  default values.
142  """
143 
144  # NOTE: To add a new key to extra_params_dict, set a default value here
145  # and update the documentation of eval_nrfit
146  default_keywords = {
147  'Lambda1': None,
148  'Lambda2': None,
149  'eccentricity': None,
150  'mean_anomaly': None,
151  'unlimited_extrapolation': False,
152  'physical_limit_violation_behavior': "ERROR",
153  }
154 
155  if extra_params_dict is None:
156  extra_params_dict = {}
157 
158  # Sanity checks
159  for key in extra_params_dict.keys():
160  if key not in default_keywords.keys():
161  raise ValueError('Invalid key %s in extra_params_dict. '%(key)
162  + "Should be one of ["+ ", ".join(default_keywords.keys()) + "].")
163 
164  # set to default if keyword is not specified
165  for key in default_keywords:
166  if key not in extra_params_dict.keys():
167  extra_params_dict[key] = default_keywords[key]
168 
169  return extra_params_dict
170 
171 
172 
173 #=============================================================================
174 def eval_nrfit(m1, m2, chiA_vec, chiB_vec, model_name, fit_types_list, f_ref=-1,
175  extra_params_dict=None):
176  """
177  Evaluates Numerical Relativity fits for a given model.
178 
179  - m1:
180  mass of object 1 in kg. \n
181  This needs to be in kg for models that allow f_ref != -1. When
182  f_ref = -1, if other units are used for m1/m2, the remnant mass will be
183  returned in the same units. \n
184  - m2:
185  mass of the object 2 in kg. \n
186  - chiA_vec:
187  dimensionless spin (3-vector) of object 1 at the reference epoch. \n
188  - chiB_vec:
189  dimensionless spin (3-vector) of object 2 at the reference epoch. \n
190  This follows the same convention as the waveform interface (see
191  https://dcc.ligo.org/T1800226/public), where the spin components are
192  defined as:
193  * \f$\chi_z = \chi \cdot \hat{L}\f$, where \f$L\f$ is the orbital
194  angular momentum vector at the reference epoch.
195  * \f$\chi_x = \chi \cdot \hat{n}\f$, where \f$n =\f$ body2 -> body1
196  is the separation vector pointing from body2 to body1 at the
197  reference epoch.
198  * \f$\chi_y = \chi \cdot (\hat{L} \times \hat{n})\f$. \n
199  These spin components are frame-independent as they are defined using
200  vector inner products. One can also think of the above spins as being
201  defined in the following reference frame: The positive z-axis is along
202  the orbital angular momentum at the reference epoch. The separation
203  vector from the object 2 to object 1 at the reference epoch is along
204  the positive x-axis. The y-axis completes the right-handed triad. The
205  returned vectors such as final spin and kick velocity are also defined
206  in the same frame.
207  - model_name:
208  model used for fit. \n
209  See lalsimulation.nrfits.eval_fits.fits_collection.keys() for all
210  available fits. Details about a particular model, including its
211  validity, can be found by doing the following: \n
212  >> from lalsimulation.nrfits.eval_fits import fits_collection \n
213  >> help(fits_collection[model_name].fit_class)
214  - fit_types_list:
215  List of fits to evaluate. \n
216  Allowed values for elements are
217  "FinalMass", "FinalSpin", "RecoilKick", and "PeakLuminosity". \n
218  Example: fit_types_list = ["FinalMass", "FinalSpin"]
219  - f_ref:
220  reference frequency (in Hz) used to set the reference epoch. \n
221  The reference epoch is set such that the orbital frequency is equal to
222  f_ref/2. If f_ref = -1, the reference epoch is taken to be the
223  time/frequency at which the fits were constructed. This can be
224  different for different models, see the documentation for the specific
225  model. Default: f_ref = -1.
226  - extra_params_dict:
227  Any additional parameters required for a specific model.
228  Default: None. \n
229  Allowed keys for extra_params_dict are:
230  * Lambda1:
231  Tidal deformability of object 1. Default: None.
232  * Lambda2:
233  Tidal deformability of object 2. Default: None.
234  * eccentricity:
235  Eccentricity at the reference epoch. Default: None.
236  * mean_anomaly:
237  Mean anomaly at the reference epoch. Default: None.
238  * unlimited_extrapolation:
239  Some models will raise an error if evaluated well outside
240  their calibration range. If unlimited_extrapolation = True,
241  then these errors are not raised and the model will still
242  produce an output. USE AT YOUR OWN RISK!! Default: False.
243  * physical_limit_violation_behavior:
244  What to do if the fit values violate physical limits
245  FinalMass > total mass, |FinalSpin| > 1, or |RecoilKick| > 1. \n
246  Allowed options are (Default: "ERROR"):
247  - "ERROR":
248  Abort with an error message.
249  - "TRUNCATE":
250  Rescale output magnitude to maximum limit, with a warning.
251  - "TRUNCATESILENT":
252  Rescale output magnitude to maximum limit, without
253  warning.
254  - "KEEP":
255  Keep output magnitude, but with a warning.
256  - "IGNORE":
257  Keep output magnitude, without warning.
258 
259 
260  - Returns: return_dict. \n
261  Dictionary of values corresponding to the keys in fit_types_list. \n
262  Example: mf = return_dict["FinalMass"]
263  """
264 
265  ### Sanity checks
266  if model_name not in fits_collection.keys():
267  raise ValueError("Invalid model_name=%s. "%model_name \
268  + "Should be one of ["+ ", ".join(fits_collection.keys()) + "].")
269 
270  if m1 < 0.09 * MSUN_SI and f_ref != -1:
271  warnings.warn("Small value of m1 = %e (kg) = %e (Msun) requested. "
272  "When f_ref != -1, component masses must be in kgs, perhaps you "
273  "are using different units?"%(m1, m1/MSUN_SI))
274 
275  if m2 < 0.09 * MSUN_SI and f_ref != -1:
276  warnings.warn("Small value of m2 = %e (kg) = %e (Msun) requested. "
277  "When f_ref != -1, component masses must be in kgs, perhaps you "
278  "are using different units?"%(m2, m2/MSUN_SI))
279 
280  if m1 <= 0 or m2 <= 0:
281  raise ValueError("Got nonpositive mass: m1=%.3e, m2=%.3e"%(m1,m2))
282 
283  chiA_vec = np.atleast_1d(chiA_vec)
284  chiB_vec = np.atleast_1d(chiB_vec)
285  if len(chiA_vec) != 3 or len(chiB_vec) != 3:
286  raise TypeError("Expected input spins to be 3-vectors.")
287 
288  if np.linalg.norm(chiA_vec) > 1:
289  raise ValueError("Invalid spin magnitude |chiA_vec|=%.3f."%( \
290  np.linalg.norm(chiA_vec)))
291 
292  if np.linalg.norm(chiB_vec) > 1:
293  raise ValueError("Invalid spin magnitude |chiB_vec|=%.3f."%( \
294  np.linalg.norm(chiB_vec)))
295 
296  if not type(fit_types_list) == list:
297  raise TypeError("fit_types_list should be a list.")
298 
299 
300  # do sanity checks on extra_params_dict and set default values if
301  # required.
302  extra_params_dict = check_extra_params_and_set_defaults(extra_params_dict)
303 
304  # Some further sanity checks to makes sure extra_params_dict is
305  # compatible with the given model
306  if (extra_params_dict['Lambda1'] is not None) \
307  or (extra_params_dict['Lambda2'] is not None):
308  if 'tidal' not in fits_collection[model_name].model_keywords:
309  raise ValueError("This model does not allow Lambda1/Lambda2.")
310 
311  if (extra_params_dict['eccentricity'] is not None) \
312  or (extra_params_dict['mean_anomaly'] is not None):
313  if 'eccentric' not in fits_collection[model_name].model_keywords:
314  raise ValueError("This model does not allow eccentricity or "
315  "mean_anomaly.")
316 
317  if 'aligned_spin' in fits_collection[model_name].model_keywords:
318  if np.linalg.norm(chiA_vec[:2]) > 0 or np.linalg.norm(chiB_vec[:2]) > 0:
319  raise ValueError("This model only allows nonprecessing spins")
320 
321  if 'nonspinning' in fits_collection[model_name].model_keywords:
322  if np.linalg.norm(chiA_vec) > 0 or np.linalg.norm(chiB_vec) > 0:
323  raise ValueError("This model only allows zero spins")
324 
325 
326  swapped_labels = False
327  # If m1 < m2, we switch the labels of the two objects, and then rotate the
328  # spins by pi about the z-direction. This amounts to a rigid rotation of
329  # the full system about the z-axis by pi. After computing the remnant
330  # properties, the final spin and kick vectors are rotated by pi to undo
331  # this switch.
332  if m1 < m2:
333  temp = m1
334  m1 = m2
335  m2 = temp
336  temp = chiA_vec
337  chiA_vec = chiB_vec
338  chiB_vec = temp
339  chiA_vec = quaternion_utils.rotate_in_plane(chiA_vec, np.pi)
340  chiB_vec = quaternion_utils.rotate_in_plane(chiB_vec, np.pi)
341  swapped_labels = True
342 
343  # Compute NR fit quantities
344  val_dict = fits_collection[model_name].fit_class(m1, m2, chiA_vec, chiB_vec,
345  f_ref, fit_types_list, extra_params_dict)
346 
347  # If the output breaks physical limits, truncate it if requested.
348  val_dict = truncate_output_to_physical_limits(val_dict, \
349  extra_params_dict['physical_limit_violation_behavior'])
350 
351  # So far FinalMass was dimless, now rescale to same units as total mass
352  if 'FinalMass' in val_dict.keys():
353  val_dict['FinalMass'] *= (m1 + m2)
354 
355  # Rotate remnant vectors by pi if the component lables were swapped
356  if swapped_labels:
357  if 'FinalSpin' in val_dict.keys():
358  val_dict['FinalSpin'] = quaternion_utils.rotate_in_plane( \
359  val_dict['FinalSpin'], np.pi)
360  if 'RecoilKick' in val_dict.keys():
361  val_dict['RecoilKick'] = quaternion_utils.rotate_in_plane( \
362  val_dict['RecoilKick'], np.pi)
363 
364  return val_dict
Class for NRSur3dq8Remnant model for the remnant mass, spin and kick velocity for nonprecessing BBH s...
Class for NRSur7dq4Remnant model for the mass, spin and kick of the final BH of generically precessin...
Saves attributes of a particular fit.
Definition: eval_fits.py:18
def __init__(self, fit_class, desc, refs, model_keywords)
Initializes a FitAttributes object.
Definition: eval_fits.py:29
def truncate_output_to_physical_limits(val_dict, behavior)
Function to truncate fit values when they go beyond the following physical limits:
Definition: eval_fits.py:93
def check_extra_params_and_set_defaults(extra_params_dict)
Does some sanity checks on extra_params_dict.
Definition: eval_fits.py:142
def eval_nrfit(m1, m2, chiA_vec, chiB_vec, model_name, fit_types_list, f_ref=-1, extra_params_dict=None)
Evaluates Numerical Relativity fits for a given model.
Definition: eval_fits.py:263