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.14/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)
###########################################
### Running Optimizer ###
###########################################
Reference waveform : SXS:BBH:0180
(q, chi1z, chi2z) : (1.00, -0.00, -0.00)
binary type : nonspinning, equal-mass, quasi-circular
Variables for ICs : ['E0byM', 'pph0']
*********************************************
Search bounds (eps) iteration #1
*********************************************
---------------------------------------------
Optimization iteration #1
---------------------------------------------
Original mismatch : 1.000e+00
Optimization interval :
E0byM : [9.888e-01,9.987e-01]
pph0 : [4.691e+00,4.786e+00]
Initial guess :
E0byM : 0.993735077287272
pph0 : 4.738442984502489
---------------------------------------------------------------------------
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,
6 use_nqc=False,
7 minimizer=minimizer,
8 opt_good_mm=5e-3,
9 opt_bounds=bounds,
10 bounds_iter=bounds_iter,
11 debug=False,
12 json_file=None,
13 overwrite=False
14 )
16 print("\nOptimization Results:")
17 print('=' * 50)
File /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/PyART/analysis/opt_ic.py:218, 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, verbose, debug)
214 print(f"{dashes}\nOptimization iteration #{j:d}\n{dashes}")
215 if (
216 i == 1 and j == 1 and opt_data is None
217 ): # if first iter of both loops
--> 218 opt_data = self.optimize_mismatch(use_ref_guess=True)
219 else:
220 opt_data_new = self.optimize_mismatch(use_ref_guess=False)
File /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/PyART/analysis/opt_ic.py:788, in Optimizer.optimize_mismatch(self, use_ref_guess, verbose)
785 f = lambda x: self.__func_to_minimize(x, kys, verbose=verbose, cache=cache)
787 t0_annealing = time.perf_counter()
--> 788 opts, mm_opt = self.minimize(f, x0, bounds_array, kys)
790 if verbose:
791 print(
792 " >> mismatch - iter : {:.3e} - {:3d}".format(
793 mm_opt, self.annealing_counter
794 ),
795 end="\r",
796 )
File /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/PyART/analysis/opt_ic.py:911, in Optimizer.__minimize_annealing_(self, f, x0, bounds_array, kys)
908 seed = self.minimizer.get("opt_seed", 190521)
909 x0 = x0
--> 911 opt_result = optimize.dual_annealing(
912 f,
913 maxfun=maxiter,
914 seed=seed,
915 x0=x0,
916 bounds=bounds_array,
917 )
919 opt_pars = opt_result["x"]
920 opts = {kys[i]: opt_pars[i] for i in range(len(kys))}
File /opt/hostedtoolcache/Python/3.11.14/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.14/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.14/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.14/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.14/x64/lib/python3.11/site-packages/PyART/analysis/opt_ic.py:785, in Optimizer.optimize_mismatch.<locals>.<lambda>(x)
783 x0 = [vs0[ky] for ky in kys]
784 bounds_array = np.array([[bounds[ky][0], bounds[ky][1]] for ky in kys])
--> 785 f = lambda x: self.__func_to_minimize(x, kys, verbose=verbose, cache=cache)
787 t0_annealing = time.perf_counter()
788 opts, mm_opt = self.minimize(f, x0, bounds_array, kys)
File /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/PyART/analysis/opt_ic.py:704, in Optimizer.__func_to_minimize(self, x, kys, verbose, cache)
702 chi2 = ref_meta["chi2z"]
703 rvec = np.linspace(2, 20, num=200)
--> 704 Vmin = PotentialMinimum(rvec, pph0, q, chi1, chi2)
705 dV = Vmin - vs["E0byM"]
706 else:
File /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/PyART/models/teob.py:440, in PotentialMinimum(rvec, pph, q, chi1, chi2)
420 def PotentialMinimum(rvec, pph, q, chi1, chi2):
421 """
422 Compute the minimum of the EOB radial potential for given parameters.
423 Parameters
(...) 438 The minimum value of the EOB radial potential.
439 """
--> 440 V = RadialPotential(rvec, pph, q, chi1, chi2)
441 peaks, _ = find_peaks(-V, height=-1)
442 if len(peaks) > 0:
File /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/PyART/models/teob.py:383, in RadialPotential(r, pph, q, chi1, chi2)
362 def RadialPotential(r, pph, q, chi1, chi2):
363 """
364 Compute the EOB radial potential for given parameters.
365 Parameters
(...) 381 The EOB radial potential values.
382 """
--> 383 return np.array([SpinHamiltonian(ri, pph, q, chi1, chi2) for ri in r])
File /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/PyART/models/teob.py:383, in <listcomp>(.0)
362 def RadialPotential(r, pph, q, chi1, chi2):
363 """
364 Compute the EOB radial potential for given parameters.
365 Parameters
(...) 381 The EOB radial potential values.
382 """
--> 383 return np.array([SpinHamiltonian(ri, pph, q, chi1, chi2) for ri in r])
File /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/PyART/models/teob.py:356, in SpinHamiltonian(r, pph, q, chi1, chi2, prstar)
332 def SpinHamiltonian(r, pph, q, chi1, chi2, prstar=0.0):
333 """
334 Compute the EOB Hamiltonian for given parameters.
335
(...) 354 The EOB Hamiltonian value.
355 """
--> 356 hatH = EOB.eob_ham_s_py(r, q, pph, prstar, chi1, chi2)
357 nu = q / (1 + q) ** 2
358 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