1 """ Evaluates Numerical Relativity fits such as remnant BH mass, spin, etc, for
9 from lal
import MSUN_SI
11 from .NRSur7dq4Remnant
import NRSur7dq4Remnant
12 from .NRSur3dq8Remnant
import NRSur3dq8Remnant
13 from .
import quaternion_utils
17 """ Saves attributes of a particular fit.
20 def __init__(self, fit_class, desc, refs, model_keywords):
21 """ Initializes a FitAttributes object. See usage examples below.
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.
32 allowed_model_keywords = [
'nonspinning',
'aligned_spin',
'precessing', \
35 for key
in model_keywords:
36 if key
not in allowed_model_keywords:
37 raise Exception(
"Invalid model_keyword: %s"%key)
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'],
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'],
72 Function to truncate fit values when they go beyond the following physical
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.
79 - val_dict: A dictionary of fits with one or more of the above keys. See
82 - behavior (str) can be one of the following:
83 * "TRUNCATE": Rescale output magnitude to maximum limit, with an
85 * "TRUNCATESILENT": Rescale output magnitude to maximum limit, without
87 * "KEEP": Keep output magnitude, even if unphysical, but
88 still give the info message
89 * "IGNORE": Keep output magnitude, even if unphysical, without
91 * "ERROR": Abort with an error if output magnitude is
96 if behavior ==
"IGNORE":
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) +
"].")
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:
116 if fit_type
in [
"FinalSpin",
"RecoilKick"]:
117 message +=
"%s magnitude=%.7e"%(fit_type, magnitude)
119 message +=
"Dimensionless %s=%.7e"%(fit_type, magnitude)
120 message +=
" is over max limit=%.2e."%(limit)
122 if behavior ==
"ERROR":
123 raise RuntimeError(message +
" Adapt " \
124 "extra_params_dict['physical_limit_violation_behavior'] to "\
125 "continue without an error.")
127 elif behavior
in [
"TRUNCATE",
"TRUNCATESILENT"]:
128 message +=
" Setting to max limit."
129 val_dict[fit_type] *= limit/magnitude
132 if behavior
in [
"KEEP",
"TRUNCATE"]:
133 warnings.warn(message)
139 """ Does some sanity checks on extra_params_dict.
140 If any of the default_keywords are not specified, sets them to the
149 'eccentricity':
None,
150 'mean_anomaly':
None,
151 'unlimited_extrapolation':
False,
152 'physical_limit_violation_behavior':
"ERROR",
155 if extra_params_dict
is None:
156 extra_params_dict = {}
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()) +
"].")
165 for key
in default_keywords:
166 if key
not in extra_params_dict.keys():
167 extra_params_dict[key] = default_keywords[key]
169 return extra_params_dict
174 def eval_nrfit(m1, m2, chiA_vec, chiB_vec, model_name, fit_types_list, f_ref=-1,
175 extra_params_dict=None):
177 Evaluates Numerical Relativity fits for a given model.
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
185 mass of the object 2 in kg. \n
187 dimensionless spin (3-vector) of object 1 at the reference epoch. \n
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
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
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
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)
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"]
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.
227 Any additional parameters required for a specific model.
229 Allowed keys for extra_params_dict are:
231 Tidal deformability of object 1. Default: None.
233 Tidal deformability of object 2. Default: None.
235 Eccentricity at the reference epoch. Default: None.
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"):
248 Abort with an error message.
250 Rescale output magnitude to maximum limit, with a warning.
252 Rescale output magnitude to maximum limit, without
255 Keep output magnitude, but with a warning.
257 Keep output magnitude, without warning.
260 - Returns: return_dict. \n
261 Dictionary of values corresponding to the keys in fit_types_list. \n
262 Example: mf = return_dict["FinalMass"]
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()) +
"].")
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))
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))
280 if m1 <= 0
or m2 <= 0:
281 raise ValueError(
"Got nonpositive mass: m1=%.3e, m2=%.3e"%(m1,m2))
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.")
288 if np.linalg.norm(chiA_vec) > 1:
289 raise ValueError(
"Invalid spin magnitude |chiA_vec|=%.3f."%( \
290 np.linalg.norm(chiA_vec)))
292 if np.linalg.norm(chiB_vec) > 1:
293 raise ValueError(
"Invalid spin magnitude |chiB_vec|=%.3f."%( \
294 np.linalg.norm(chiB_vec)))
296 if not type(fit_types_list) == list:
297 raise TypeError(
"fit_types_list should be a list.")
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.")
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 "
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")
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")
326 swapped_labels =
False
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
344 val_dict = fits_collection[model_name].
fit_class(m1, m2, chiA_vec, chiB_vec,
345 f_ref, fit_types_list, extra_params_dict)
349 extra_params_dict[
'physical_limit_violation_behavior'])
352 if 'FinalMass' in val_dict.keys():
353 val_dict[
'FinalMass'] *= (m1 + m2)
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)
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.
def __init__(self, fit_class, desc, refs, model_keywords)
Initializes a FitAttributes object.
def truncate_output_to_physical_limits(val_dict, behavior)
Function to truncate fit values when they go beyond the following physical limits:
def check_extra_params_and_set_defaults(extra_params_dict)
Does some sanity checks on extra_params_dict.
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.