Optimizing Initial Conditions for EOB Waveforms¶
This tutorial demonstrates how to optimize the initial conditions of Effective One Body (EOB) waveforms to best match Numerical Relativity (NR) simulations.
The optimizer searches for the best initial energy (E₀) and angular momentum (pₚₕ₀) that minimize the mismatch between EOB and NR waveforms.
Setup¶
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import os
import numpy as np
import matplotlib.pyplot as plt
from PyART.analysis.opt_ic import Optimizer
from PyART.catalogs.sxs import Waveform_SXS
from PyART.catalogs.rit import Waveform_RIT
from PyART.analysis.match import Matcher
/opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/PyART/analysis/match.py:15: UserWarning: Wswiglal-redir-stdio:
SWIGLAL standard output/error redirection is enabled in IPython.
This may lead to performance penalties. To disable locally, use:
with lal.no_swig_redirect_standard_output_error():
...
To disable globally, use:
lal.swig_redirect_standard_output_error(False)
Note however that this will likely lead to error messages from
LAL functions being either misdirected or lost when called from
Jupyter notebooks.
To suppress this warning, use:
import warnings
warnings.filterwarnings("ignore", "Wswiglal-redir-stdio")
import lal
import lal
WARNING: TEOBResumS not installed.
Load NR Waveform¶
First, we load a numerical relativity waveform that we want to match with EOB:
# Load an SXS waveform
catalog = 'sxs'
sim_id = 180
if catalog == 'sxs':
ebbh = Waveform_SXS(
path='./',
download=True,
ID=sim_id,
order="Extrapolated_N3.dir",
ellmax=7,
nu_rescale=True
)
# Remove junk radiation from the beginning
ebbh.cut(200)
else:
# For RIT or other catalogs
ebbh = Waveform_RIT(
path='./local_data/rit/',
download=True,
ID=sim_id,
nu_rescale=True
)
# Display metadata
print('Waveform Metadata:')
print('=' * 50)
for k, v in ebbh.metadata.items():
print(f'{k:15s} : {v}')
print('=' * 50)
Waveform Metadata:
==================================================
name : SXS:BBH:0180
ref_time : 250.0
m1 : 0.499999985387
m2 : 0.499999985116
M : 0.999999970503
q : 1.000000000542
nu : 0.24999999999999994
S1 : [-6.95604332e-13 3.90134701e-13 -9.17626593e-10]
S2 : [ 5.69749564e-13 -5.89768887e-13 -5.39321192e-10]
chi1x : -2.78241749107e-12
chi1y : 1.5605388945e-12
chi1z : -3.67050658542e-09
chi2x : 2.27899839037e-12
chi2y : -2.35907568866e-12
chi2z : -2.1572848954e-09
LambdaAl2 : 0.0
LambdaBl2 : 0.0
pos1 : [-9.23165599e+00 6.45733324e-01 3.68318004e-10]
pos2 : [ 9.23165601e+00 -6.45734196e-01 5.40638974e-10]
r0 : 18.50842452509741
e0 : 5.11e-05
f0v : [-4.00190520e-14 -5.11164262e-14 3.90474262e-03]
f0 : 0.003904742624312793
E0 : 0.9937350479750683
P0 : [ 5.000e-15 3.539e-13 -3.778e-13]
J0 : [-1.56989012e-06 -3.87578732e-07 1.18461067e+00]
Jz0 : 1.184610674783749
E0byM : 0.9937350772872718
pph0 : 4.738442984502489
Mf : 0.951614826833
afv : [-2.14539981e-13 -8.96386037e-12 6.86429827e-01]
af : 0.686429826547
scat_angle : None
flags : ['nonspinning', 'equal-mass', 'quasi-circular']
==================================================
Configure Optimizer Settings¶
The optimizer requires several settings:
Mismatch settings: Parameters for computing mismatches
Minimizer: Algorithm and parameters for optimization
Bounds: Search ranges for initial conditions
Iteration settings: How to adaptively adjust bounds
# Mismatch computation settings
mm_settings = {
'cut_second_waveform': True,
'initial_frequency_mm': 10,
'M': 100, # Total mass in solar masses
'final_frequency_mm': 1024,
'taper_alpha': 0.50,
'taper_start': 0.10,
}
# Add catalog-specific alignment settings
if catalog == 'rit':
mm_settings['pre_align_shift'] = 100.
elif catalog == 'sxs':
mm_settings['pre_align_shift'] = 0.
# Optimizer settings
kind_ic = 'E0pph0' # Optimize E0 (energy) and pph0 (angular momentum)
# Bounds for the optimization
bounds = {'E0byM': [None, None], 'pph0': [None, None]}
# Adaptive bounds iteration
bounds_iter = {
'eps_initial': {'E0byM': 5e-3, 'pph0': 1e-2},
'eps_factors': {'E0byM': 4, 'pph0': 2},
'bad_mm': 1e-2,
'max_iter': 3
}
# Minimizer configuration (using dual annealing)
minimizer = {
'kind': 'dual_annealing',
'opt_maxfun': 100, # Maximum function evaluations
'opt_max_iter': 1, # Maximum optimization iterations
'opt_seed': 190521
}
print("Optimizer configured with:")
print(f" Kind: {kind_ic}")
print(f" Minimizer: {minimizer['kind']}")
print(f" Target mass: {mm_settings['M']} Msun")
print(f" Good mismatch threshold: 5e-3")
Optimizer configured with:
Kind: E0pph0
Minimizer: dual_annealing
Target mass: 100 Msun
Good mismatch threshold: 5e-3
Run the Optimizer¶
Now we create the optimizer and let it find the best initial conditions:
# Create and run optimizer
opt = Optimizer(
ebbh,
kind_ic=kind_ic,
mm_settings=mm_settings,
use_nqc=False,
minimizer=minimizer,
opt_good_mm=5e-3,
opt_bounds=bounds,
bounds_iter=bounds_iter,
debug=False,
json_file=None,
overwrite=False
)
print("\nOptimization Results:")
print('=' * 50)
print(f"Optimized E0/M : {opt.opt_Waveform.metadata.get('E0byM', 'N/A')}")
print(f"Optimized pph0 : {opt.opt_Waveform.metadata.get('pph0', 'N/A')}")
print(f"Final mismatch : {opt.opt_mismatch:.5e}")
print('=' * 50)
WARNING:root:Error occurred in EOB wave generation:
name 'EOB' is not defined
WARNING:root:Error occurred in EOB wave generation:
name 'EOB' is not defined
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[4], line 2
1 # Create and run optimizer
----> 2 opt = Optimizer(
3 ebbh,
4 kind_ic=kind_ic,
5 mm_settings=mm_settings,
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/PyART/analysis/opt_ic.py:270, in Optimizer.__init__(self, ref_Waveform, kind_ic, vrs, map_function, use_nqc, r0_eob, model_opts, opt_max_iter, opt_good_mm, opt_bounds, bounds_iter, minimizer, use_matcher_cache, json_file, overwrite, json_save_dyn, mm_settings, objective_settings, verbose, debug)
264 logging.info(
265 f"{dashes}\nOptimization iteration #{j:d}\n{dashes}"
266 )
267 if (
268 i == 1 and j == 1 and opt_data is None
269 ): # if first iter of both loops
--> 270 opt_data = self.optimize_mismatch(use_ref_guess=True)
271 else:
272 opt_data_new = self.optimize_mismatch(use_ref_guess=False)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/PyART/analysis/opt_ic.py:969, in Optimizer.optimize_mismatch(self, use_ref_guess, verbose)
966 f = lambda x: self.__func_to_minimize(x, kys, verbose=verbose, cache=cache)
968 t0_annealing = time.perf_counter()
--> 969 opts, mm_opt = self.minimize(f, x0, bounds_array, kys)
971 if verbose:
972 self._log_progress_line(
973 " >> mismatch - iter : {:.3e} - {:3d}".format(
974 mm_opt, self.annealing_counter
975 )
976 )
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/PyART/analysis/opt_ic.py:1094, in Optimizer.__minimize_annealing_(self, f, x0, bounds_array, kys)
1091 seed = self.minimizer.get("opt_seed", 190521)
1092 x0 = x0
-> 1094 opt_result = optimize.dual_annealing(
1095 f,
1096 maxfun=maxiter,
1097 seed=seed,
1098 x0=x0,
1099 bounds=bounds_array,
1100 )
1102 opt_pars = opt_result["x"]
1103 opts = {kys[i]: opt_pars[i] for i in range(len(kys))}
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/scipy/_lib/_util.py:352, in _transition_to_rng.<locals>.decorator.<locals>.wrapper(*args, **kwargs)
345 message = (
346 "The NumPy global RNG was seeded by calling "
347 f"`np.random.seed`. Beginning in {end_version}, this "
348 "function will no longer use the global RNG."
349 ) + cmn_msg
350 warnings.warn(message, FutureWarning, stacklevel=2)
--> 352 return fun(*args, **kwargs)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/scipy/optimize/_dual_annealing.py:674, in dual_annealing(func, bounds, args, maxiter, minimizer_kwargs, initial_temp, restart_temp_ratio, visit, accept, maxfun, rng, no_local_search, callback, x0)
672 # Initialization of the energy state
673 energy_state = EnergyState(lower, upper, callback)
--> 674 energy_state.reset(func_wrapper, rng_gen, x0)
675 # Minimum value of annealing temperature reached to perform
676 # re-annealing
677 temperature_restart = initial_temp * restart_temp_ratio
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/scipy/optimize/_dual_annealing.py:173, in EnergyState.reset(self, func_wrapper, rng_gen, x0)
171 reinit_counter = 0
172 while init_error:
--> 173 self.current_energy = func_wrapper.fun(self.current_location)
174 if self.current_energy is None:
175 raise ValueError('Objective function is returning None')
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/scipy/optimize/_dual_annealing.py:381, in ObjectiveFunWrapper.fun(self, x)
379 def fun(self, x):
380 self.nfev += 1
--> 381 return self.func(x, *self.args)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/PyART/analysis/opt_ic.py:966, in Optimizer.optimize_mismatch.<locals>.<lambda>(x)
964 x0 = [vs0[ky] for ky in kys]
965 bounds_array = np.array([[bounds[ky][0], bounds[ky][1]] for ky in kys])
--> 966 f = lambda x: self.__func_to_minimize(x, kys, verbose=verbose, cache=cache)
968 t0_annealing = time.perf_counter()
969 opts, mm_opt = self.minimize(f, x0, bounds_array, kys)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/PyART/analysis/opt_ic.py:880, in Optimizer.__func_to_minimize(self, x, kys, verbose, cache)
878 chi2 = ref_meta["chi2z"]
879 rvec = np.linspace(2, 20, num=200)
--> 880 Vmin = PotentialMinimum(rvec, pph0, q, chi1, chi2)
881 dV = Vmin - vs["E0byM"]
882 else:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/PyART/models/teob.py:450, in PotentialMinimum(rvec, pph, q, chi1, chi2)
430 def PotentialMinimum(rvec, pph, q, chi1, chi2):
431 """
432 Compute the minimum of the EOB radial potential for given parameters.
433 Parameters
(...) 448 The minimum value of the EOB radial potential.
449 """
--> 450 V = RadialPotential(rvec, pph, q, chi1, chi2)
451 peaks, _ = find_peaks(-V, height=-1)
452 if len(peaks) > 0:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/PyART/models/teob.py:393, in RadialPotential(r, pph, q, chi1, chi2)
372 def RadialPotential(r, pph, q, chi1, chi2):
373 """
374 Compute the EOB radial potential for given parameters.
375 Parameters
(...) 391 The EOB radial potential values.
392 """
--> 393 return np.array([SpinHamiltonian(ri, pph, q, chi1, chi2) for ri in r])
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/PyART/models/teob.py:393, in <listcomp>(.0)
372 def RadialPotential(r, pph, q, chi1, chi2):
373 """
374 Compute the EOB radial potential for given parameters.
375 Parameters
(...) 391 The EOB radial potential values.
392 """
--> 393 return np.array([SpinHamiltonian(ri, pph, q, chi1, chi2) for ri in r])
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/PyART/models/teob.py:366, in SpinHamiltonian(r, pph, q, chi1, chi2, prstar)
342 def SpinHamiltonian(r, pph, q, chi1, chi2, prstar=0.0):
343 """
344 Compute the EOB Hamiltonian for given parameters.
345
(...) 364 The EOB Hamiltonian value.
365 """
--> 366 hatH = EOB.eob_ham_s_py(r, q, pph, prstar, chi1, chi2)
367 nu = q / (1 + q) ** 2
368 E0 = nu * hatH[0]
NameError: name 'EOB' is not defined
Visualize: Mismatch vs Total Mass¶
Let’s see how the optimized waveform performs across different total masses:
# Compute mismatch for a range of masses
masses = np.linspace(20, 200, num=19)
mm = np.zeros_like(masses)
for i, M in enumerate(masses):
mm_settings['M'] = M
matcher = Matcher(ebbh, opt.opt_Waveform, settings=mm_settings)
mm[i] = matcher.mismatch
if i % 5 == 0:
print(f'M = {M:6.1f} Msun, mismatch = {mm[i]:.3e}')
# Plot
plt.figure(figsize=(10, 6))
plt.plot(masses, mm, linewidth=2, marker='o', markersize=6)
plt.yscale('log')
plt.xlabel(r'Total Mass $M$ [$M_\odot$]', fontsize=16)
plt.ylabel(r'Mismatch $\bar{\mathcal{F}}$', fontsize=16)
plt.title('Optimized EOB vs NR Mismatch', fontsize=18)
plt.ylim(1e-4, 1e-1)
plt.grid(True, alpha=0.3, which='both')
plt.axhline(y=5e-3, color='r', linestyle='--', alpha=0.5, label='Good mismatch threshold')
plt.axhline(y=1e-2, color='orange', linestyle='--', alpha=0.5, label='Acceptable threshold')
plt.legend(fontsize=12)
plt.tight_layout()
plt.show()
print(f'\nMinimum mismatch: {mm.min():.3e} at M = {masses[mm.argmin()]:.1f} Msun')
Compare Waveforms¶
Let’s visually compare the NR and optimized EOB waveforms:
# Get merger times
nr_mrg, _, _, _ = ebbh.find_max()
eob_mrg, _, _, _ = opt.opt_Waveform.find_max()
plt.figure(figsize=(14, 5))
# Real part
plt.subplot(1, 2, 1)
plt.plot(ebbh.u - nr_mrg, ebbh.hlm[(2,2)]['real'],
label='NR', linewidth=2, alpha=0.8)
plt.plot(opt.opt_Waveform.u - eob_mrg, opt.opt_Waveform.hlm[(2,2)]['real'],
label='Optimized EOB', linewidth=2, alpha=0.8)
plt.xlabel('Time (M)', fontsize=14)
plt.ylabel(r'Re[$h_{22}$]', fontsize=14)
plt.title('Waveform Comparison', fontsize=16)
plt.legend(fontsize=12)
plt.grid(True, alpha=0.3)
plt.xlim([-500, 100])
# Amplitude
plt.subplot(1, 2, 2)
plt.plot(ebbh.u - nr_mrg, ebbh.hlm[(2,2)]['A'],
label='NR', linewidth=2, alpha=0.8)
plt.plot(opt.opt_Waveform.u - eob_mrg, opt.opt_Waveform.hlm[(2,2)]['A'],
label='Optimized EOB', linewidth=2, alpha=0.8)
plt.xlabel('Time (M)', fontsize=14)
plt.ylabel(r'$|h_{22}|$', fontsize=14)
plt.title('Amplitude Comparison', fontsize=16)
plt.legend(fontsize=12)
plt.grid(True, alpha=0.3)
plt.xlim([-500, 100])
plt.tight_layout()
plt.show()
Summary¶
This tutorial demonstrated:
Loading NR waveforms from various catalogs (SXS, RIT)
Configuring the optimizer with appropriate settings
Running the optimization to find best-fit initial conditions
Evaluating the optimized waveform across different masses
Visualizing the agreement between NR and optimized EOB
Key Points¶
The optimizer uses dual annealing by default for robust global optimization
Initial conditions (E₀, pₚₕ₀) significantly affect EOB waveform quality
Mismatches typically vary with total mass
Good mismatches (< 5×10⁻³) indicate excellent EOB-NR agreement
Next Steps¶
Try different catalogs (RIT, CoRe, ICCUB)
Experiment with different initial condition parameterizations (e0f0 vs E0pph0)
Use NQC (Non-Quasi-Circular) corrections for eccentric binaries
Optimize for multiple modes simultaneously