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)
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
Cell In[2], line 6
2 catalog = 'sxs'
3 sim_id = 180
4
5 if catalog == 'sxs':
----> 6 ebbh = Waveform_SXS(
7 path='./',
8 download=True,
9 ID=sim_id,
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/PyART/catalogs/sxs.py:173, in Waveform_SXS.__init__(self, path, ID, order, level, cut_N, cut_U, ellmax, load, download, downloads, load_m0, nu_rescale, src, ignore_deprecation, basename)
171 break
172 elif lvn == ref_lv_min:
--> 173 raise RuntimeError(
174 f"No data for ref-levels:[{ref_lv_min:d},{ref_lv_max:d}] found"
175 )
177 else:
178 raise RuntimeError(f"Invalid input for level: {self.level}")
RuntimeError: No data for ref-levels:[1,7] found
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")
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)
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