Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions docs/source/utilities.rst
Original file line number Diff line number Diff line change
Expand Up @@ -74,8 +74,8 @@ This function takes in quite a few parameters, due to its overlapping usage with
chif_guess=0.5,
t0=t0,
tf=tf,
force_Omega0_optimization=False, #True if we want to perform a least squares fit for Omega0.
include_Omega0_as_parameter=True, #Use Omega0 as a third parameter in the optimization algorithm, along with the final mass and spin.
force_Omega0_optimization=False, #True if we want to perform a least squares fit for Omega_0.
include_Omega0_as_parameter=True, #Use Omega_0 as a third parameter in the optimization algorithm, along with the final mass and spin.
start_with_wide_search=False, #If True we start with the scipy differential_evolution algorithm and finish with the scipy minimize algorithm
t_shift_range = np.arange(-10,10,0.1),#Time shift array to be used in the mismatch calculation
)
Expand Down
116 changes: 58 additions & 58 deletions gwBOB/BOB_terms.py

Large diffs are not rendered by default.

140 changes: 70 additions & 70 deletions gwBOB/BOB_terms_jax.py

Large diffs are not rendered by default.

20 changes: 10 additions & 10 deletions gwBOB/BOB_utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -627,9 +627,9 @@ def get_correct_Phi_and_Omega(self):
ValueError: If ``what_should_BOB_create`` is not one of the
supported flavors.
'''
#Even in the cases of strain_using_news, we still want to use the news frequency in all of the Omega0 optimizations because the analytical news frequency term
#Even in the cases of strain_using_news, we still want to use the news frequency in all of the Omega_0 optimizations because the analytical news frequency term
#is built assuming the BOB amplitude best describes the news. While in principle, the accuracy could be improved for strain_using_news (and all X_using_Y cases)
#by optimizing Omega0 against the NR strain frequency, this would be unphysical.
#by optimizing Omega_0 against the NR strain frequency, this would be unphysical.
if('psi4' in self._wf_config.what_to_create):
if(self.minf_t0 is True):
Phi,Omega = BOB_terms.BOB_psi4_phase(self)
Expand Down Expand Up @@ -829,7 +829,7 @@ def fit_Omega0(self):
ValueError: If ``minf_t0`` is False (use ``fit_t0`` for finite t0).
'''
if(self.minf_t0 is False):
raise ValueError("You are setup for a finite t0 right now. Omega0 fitting is only defined for t0 = infinity.")
raise ValueError("You are setup for a finite t0 right now. Omega_0 fitting is only defined for t0 = infinity.")
if(self._wf_config.end_after_tpeak<self.end_fit_after_tpeak):
logger.warning("end_after_tpeak is less than end_fit_after_tpeak. Setting end_fit_after_tpeak to end_after_tpeak")
self.end_fit_after_tpeak = self._wf_config.end_after_tpeak
Expand Down Expand Up @@ -921,7 +921,7 @@ def construct_BOB_finite_t0(self,N):
'''
#Perform parameter sanity checks
if(self.optimize_Omega0):
raise ValueError("Cannot optimize Omega0 for finite t0 values.")
raise ValueError("Cannot optimize Omega_0 for finite t0 values.")

if(self._fit_config.optimize_t0_and_Omega0):
self.fit_t0_and_Omega0()
Expand Down Expand Up @@ -1364,11 +1364,11 @@ def initialize_with_sxs_data(self,sxs_id,l=2,m=2,download=True,resample_dt = 0.1

self.chif = np.linalg.norm(self.chif)
self.chif_with_sign = self.chif*sign
self.Omega_ISCO = np.abs(gen_utils.get_Omega_isco(self.chif_with_sign,self.mf))
self.Omega_ISCO = np.abs(gen_utils.get_Omega_isco(self.mf, self.chif_with_sign))
self.Omega_0 = self.Omega_ISCO
self.l = l
self.m = m
w_r,tau = gen_utils.get_qnm(self.chif,self.mf,self.l,np.abs(self.m),n=0,sign=sign)
w_r,tau = gen_utils.get_qnm(self.mf, self.chif, self.l, np.abs(self.m), n=0, sign=sign)
self.w_r = np.abs(w_r)
self.tau = np.abs(tau)
self.Omega_QNM = self.w_r/np.abs(self.m)
Expand Down Expand Up @@ -1524,12 +1524,12 @@ def initialize_with_cce_data(self,cce_id,l=2,m=2,perform_superrest_transformatio
sign = np.sign(self.chif[2])
self.chif = np.linalg.norm(self.chif)
self.chif_with_sign = self.chif*sign
self.Omega_ISCO = np.abs(gen_utils.get_Omega_isco(self.chif_with_sign,self.mf))
self.Omega_ISCO = np.abs(gen_utils.get_Omega_isco(self.mf, self.chif_with_sign))
self.Omega_0 = self.Omega_ISCO
self.l = l
self.m = m

w_r,tau = gen_utils.get_qnm(self.chif,self.mf,self.l,np.abs(self.m),n=0,sign=sign)
w_r,tau = gen_utils.get_qnm(self.mf, self.chif, self.l, np.abs(self.m), n=0, sign=sign)
self.w_r = np.abs(w_r)
self.tau = np.abs(tau)
self.Omega_QNM = self.w_r/np.abs(self.m)
Expand Down Expand Up @@ -1665,12 +1665,12 @@ def initialize_with_NR_mode(self,t_NR,y_psi4,y_strain,mf,chif,l=2,m=2,w_r = -1,t
self.l = l
self.m = m

self.Omega_ISCO = np.abs(gen_utils.get_Omega_isco(self.chif_with_sign,self.mf))
self.Omega_ISCO = np.abs(gen_utils.get_Omega_isco(self.mf, self.chif_with_sign))
self.Omega_0 = self.Omega_ISCO

if(w_r<0 or tau<0):
logger.info("Calculating Kerr QNM parameters from provided Mf and chif")
w_r,tau = gen_utils.get_qnm(self.chif,self.mf,self.l,np.abs(self.m),n=0,sign=sign)
w_r,tau = gen_utils.get_qnm(self.mf, self.chif, self.l, np.abs(self.m), n=0, sign=sign)
self.w_r = np.abs(w_r)
self.tau = np.abs(tau)
else:
Expand Down
88 changes: 42 additions & 46 deletions gwBOB/gen_utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -149,16 +149,16 @@ def get_frequency(ts):
if(freq.y[find_nearest_index(freq.t,tp)]<0):
freq.y = -freq.y
return kuibit_ts(ts.t,freq.y)
def get_r_isco(chi_with_sign,M):
def get_r_isco(mf, chif_with_sign):
'''
Compute the prograde ISCO radius for a Kerr black hole. For negative chi_with_sign values, we return the retrograde value
Compute the prograde ISCO radius for a Kerr black hole. For negative chif_with_sign values, we return the retrograde value

Parameters
----------
chi_with_sign : float
Dimensionless spin of the remnant. Negative value indicates a spin pointed opposite the direction of the initial anular momentum of the binary.
M : float
mf : float
Mass of the remnant (in geometric units).
chif_with_sign : float
Dimensionless spin of the remnant. Negative value indicates a spin pointed opposite the direction of the initial anular momentum of the binary.

Returns
-------
Expand All @@ -167,29 +167,29 @@ def get_r_isco(chi_with_sign,M):
'''
#Bardeen Press Teukolskly eq 2.21
#defined for prograde orbits
if(chi_with_sign>=0):
if(chif_with_sign>=0):
sign = 1.0
else:
sign = -1.0
chi = np.abs(chi_with_sign)
a = chi*M
a_M = a/M
chi = np.abs(chif_with_sign)
a = chi*mf
a_M = a/mf

z1 = 1 + (((1-a_M**2)**(1./3.)) * ((1+a_M)**(1./3.) + (1-a_M)**(1./3.))) #good
z2 = (3*(a_M**2) + z1**2)**0.5 #good
r_isco = M * (3 + z2 - (sign)*((3-z1)*(3+z1+2*z2))**0.5) #good
r_isco = mf * (3 + z2 - (sign)*((3-z1)*(3+z1+2*z2))**0.5) #good
return r_isco
def get_Omega_isco(chi_with_sign,M):
def get_Omega_isco(mf, chif_with_sign):
'''
Compute the orbital angular velocity at the ISCO for a Kerr black hole.
For negative chi_with_sign values, we return the retrograde value
For negative chif_with_sign values, we return the retrograde value

Parameters
----------
chi_with_sign : float
Dimensionless spin of the remnant. Negative value indicates a spin pointed opposite the direction of the initial anular momentum of the binary.
M : float
mf : float
Mass of the remnant (in geometric units).
chif_with_sign : float
Dimensionless spin of the remnant. Negative value indicates a spin pointed opposite the direction of the initial anular momentum of the binary.

Returns
-------
Expand All @@ -198,26 +198,26 @@ def get_Omega_isco(chi_with_sign,M):
'''
#Bardeen Press Teukolskly eq 2.16
#defined for prograde orbits
if(chi_with_sign>=0):
if(chif_with_sign>=0):
sign = 1.0
else:
sign = -1.0
chi = np.abs(chi_with_sign)
r_isco = get_r_isco(chi_with_sign,M)
a = chi*M
Omega = sign*np.sqrt(M)/(r_isco**1.5 + sign*a*np.sqrt(M)) # = dphi/dt
chi = np.abs(chif_with_sign)
r_isco = get_r_isco(mf, chif_with_sign)
a = chi*mf
Omega = sign*np.sqrt(mf)/(r_isco**1.5 + sign*a*np.sqrt(mf)) # = dphi/dt
return Omega

def get_qnm(chif,Mf,l,m,n=0,sign=1):
def get_qnm(mf, chif, l, m, n=0, sign=1):
'''
Get the fundamental quasinormal mode frequency components for a Kerr black hole.

Parameters
----------
mf : float
Final black-hole mass (geometric units).
chif : float
Dimensionless final spin magnitude (chi).
Mf : float
Final black-hole mass (geometric units).
l : int
Spherical-harmonic index ``l``.
m : int
Expand All @@ -232,15 +232,14 @@ def get_qnm(chif,Mf,l,m,n=0,sign=1):
tuple[float, float]
``(w_r, tau)`` where ``w_r`` is the real angular frequency and ``tau`` is the damping time.
'''
#omega_qnm, all_C, ells = qnmfits.read_qnms.qnm_from_tuple((l,m,n,1),chif,M=M)
if(sign==-1):
grav_lmn = qnm.modes_cache(s=-2,l=l,m=-m,n=n)
omega_qnm, A, C = grav_lmn(a=chif)#qnm package uses M = 1 so a = chi here
omega_qnm = -np.conj(omega_qnm)
else:
grav_lmn = qnm.modes_cache(s=-2,l=l,m=m,n=n)
omega_qnm, A, C = grav_lmn(a=chif)#qnm package uses M = 1 so a = chi here
omega_qnm /= Mf #rescale to remnant black hole mass
omega_qnm /= mf #rescale to remnant black hole mass
w_r = np.abs(omega_qnm.real)
imag_qnm = np.abs(omega_qnm.imag)
tau = 1./imag_qnm
Expand Down Expand Up @@ -454,17 +453,17 @@ def estimate_parameters(BOB,
tf : float, optional
End time after peak for mismatch window, by default 75.
force_Omega0_optimization : bool, optional
If True, enforce Omega0 optimization during construction, by default False.
If True, enforce Omega_0 optimization during construction, by default False.
NR_data : kuibit_ts, optional
Reference data to compare against. If None, inferred from ``BOB`` selection.
make_current_naturally : bool, optional
If True, construct current quadrupole via natural method, by default False.
make_mass_naturally : bool, optional
If True, construct mass quadrupole via natural method, by default False.
include_Omega0_as_parameter : bool, optional
Include Omega0 as an optimization parameter, by default False.
Include Omega_0 as an optimization parameter, by default False.
include_2Omega0_as_parameters : bool, optional
Include both lm and lmm Omega0 as parameters (quadrupole builds), by default False.
Include both lm and lmm Omega_0 as parameters (quadrupole builds), by default False.
perform_phase_alignment_first : bool, optional
If True, perform phase alignment before quadrupole combination, by default False.
start_with_wide_search : bool, optional
Expand Down Expand Up @@ -545,7 +544,7 @@ def create_guess(x):
if(include_Omega0_as_parameter):
#keep this for ordinary (l,m) &(l,-m) modes
BOB.Omega_0 = lm_Omega0_guess
w_r,tau = get_qnm(BOB.chif,BOB.mf,BOB.l,np.abs(BOB.m),sign=np.sign(BOB.chif_with_sign))
w_r,tau = get_qnm(BOB.mf, BOB.chif, BOB.l, np.abs(BOB.m), sign=np.sign(BOB.chif_with_sign))
BOB.Omega_QNM = w_r/np.abs(BOB.m)
BOB.Phi_0 = 0
BOB.tau = tau
Expand Down Expand Up @@ -922,13 +921,13 @@ def load_lower_lev_SXS(sim):
else:
raise ValueError("only one Level found")
return sim_lower
def Omega_0_fit_psi4(Mf,chif_with_sign):
def Omega_0_fit_psi4(mf, chif_with_sign):
'''
Omega0 for psi4 using the fit from Kankani and McWilliams (2025)
Omega_0 for psi4 using the fit from Kankani and McWilliams (2025)

Parameters
----------
Mf : float
mf : float
Remnant mass.
chif_with_sign : float
Remnant spin. (negative values indicate a final spin pointing opposite to the initial orbital angular momentum)
Expand All @@ -937,19 +936,18 @@ def Omega_0_fit_psi4(Mf,chif_with_sign):
-------
float
Omega_0 fit.

'''
A = 1.42968337
B = 0.08424419
C = -1.22848524
return A*Mf + B*chif_with_sign + C
def Omega_0_fit_news(Mf,chif_with_sign):
return A*mf + B*chif_with_sign + C
def Omega_0_fit_news(mf, chif_with_sign):
'''
Omega0 for news using the fit from Kankani and McWilliams (2025)
Omega_0 for news using the fit from Kankani and McWilliams (2025)

Parameters
----------
Mf : float
mf : float
Remnant mass.
chif_with_sign : float
Remnant spin. (negative values indicate a final spin pointing opposite to the initial orbital angular momentum)
Expand All @@ -958,19 +956,18 @@ def Omega_0_fit_news(Mf,chif_with_sign):
-------
float
Omega_0 fit.

'''
A = 0.33568227
B = 0.03450997
C = -0.18763176
return A*Mf + B*chif_with_sign + C
def Omega_0_fit_strain(Mf,chif_with_sign):
return A*mf + B*chif_with_sign + C
def Omega_0_fit_strain(mf, chif_with_sign):
'''
Omega0 for strain using the fit from Kankani and McWilliams (2025)
Omega_0 for strain using the fit from Kankani and McWilliams (2025)

Parameters
----------
Mf : float
mf : float
Remnant mass.
chif_with_sign : float
Remnant spin. (negative values indicate a final spin pointing opposite to the initial orbital angular momentum)
Expand All @@ -979,12 +976,11 @@ def Omega_0_fit_strain(Mf,chif_with_sign):
-------
float
Omega_0 fit.

'''
A = 0.01663248
B = 0.01798275
C = 0.07882578
return A*Mf + B*chif_with_sign + C
return A*mf + B*chif_with_sign + C



Expand Down
2 changes: 1 addition & 1 deletion tests/unit/test_BOB_orchestrator.py
Original file line number Diff line number Diff line change
Expand Up @@ -187,7 +187,7 @@ def test_construct_BOB_finite_t0_rejects_optimize_Omega0(self):
bob = BOB_utils.BOB()
bob.optimize_Omega0 = True
# Force the finite-t0 branch by direct call (skipping the setter)
with pytest.raises(ValueError, match="Cannot optimize Omega0 for finite t0"):
with pytest.raises(ValueError, match="Cannot optimize Omega_0 for finite t0"):
bob.construct_BOB_finite_t0(N=1)

def test_setting_set_initial_time_before_what_to_create_raises(self):
Expand Down
38 changes: 19 additions & 19 deletions tests/unit/test_gen_utils_math.py
Original file line number Diff line number Diff line change
Expand Up @@ -16,34 +16,34 @@


def test_get_r_isco_values():
"""ISCO radius for representative (chi, M) pairs.
"""ISCO radius for representative (mf, chif) pairs.

Reference values: chi=0 is Schwarzschild (r_isco = 6M); other values
Reference values: chif=0 is Schwarzschild (r_isco = 6M); other values
were computed once and treated as the trusted reference.
"""
chi_vals = np.array([0.0, 0.5, 0.9])
M_vals = np.array([1.0, 2.0, 5.0])
mf_vals = np.array([1.0, 2.0, 5.0])
chif_vals = np.array([0.0, 0.5, 0.9])
expected = [
6.0, # (chi=0, M=1) -> Schwarzschild ISCO = 6M
8.466005059061652, # (chi=0.5, M=2)
11.604415208809435, # (chi=0.9, M=5)
6.0, # (mf=1, chif=0) -> Schwarzschild ISCO = 6M
8.466005059061652, # (mf=2, chif=0.5)
11.604415208809435, # (mf=5, chif=0.9)
]
for chi, M, exp in zip(chi_vals, M_vals, expected):
result = gen_utils.get_r_isco(chi, M)
for mf, chif, exp in zip(mf_vals, chif_vals, expected):
result = gen_utils.get_r_isco(mf, chif)
assert np.isclose(result, exp, rtol=1e-8)


def test_get_Omega_isco_values():
"""ISCO orbital frequency for representative (chi, M) pairs."""
chi_vals = np.array([0.0, 0.5, 0.9])
M_vals = np.array([1.0, 2.0, 5.0])
"""ISCO orbital frequency for representative (mf, chif) pairs."""
mf_vals = np.array([1.0, 2.0, 5.0])
chif_vals = np.array([0.0, 0.5, 0.9])
expected = [
0.06804138174397717,
0.05429417949013838,
0.0450883417670616,
]
for chi, M, exp in zip(chi_vals, M_vals, expected):
result = gen_utils.get_Omega_isco(chi, M)
for mf, chif, exp in zip(mf_vals, chif_vals, expected):
result = gen_utils.get_Omega_isco(mf, chif)
assert np.isclose(result, exp, rtol=1e-8)


Expand All @@ -52,8 +52,8 @@ def test_get_qnm():

Spans ell=2,3, n=0,1, retrograde and prograde branches.
"""
chi_vals = np.array([0.0, 0.0, 0.0, 0.5, 0.5, 0.5])
M_vals = np.array([1.0, 1.0, 1.0, 1.0, 1.0, 2.0])
mf_vals = np.array([1.0, 1.0, 1.0, 1.0, 1.0, 2.0])
chif_vals = np.array([0.0, 0.0, 0.0, 0.5, 0.5, 0.5])
l_vals = np.array([2, 3, 2, 2, 2, 2])
m_vals = np.array([2, 2, 2, 2, 2, 2])
n_vals = np.array([0, 0, 1, 0, 0, 0])
Expand All @@ -68,10 +68,10 @@ def test_get_qnm():
11.231973996651769, 11.676945396785948, 23.353890793571896,
])

for chi, M, l, m, n, sgn, exp_w, exp_tau in zip(
chi_vals, M_vals, l_vals, m_vals, n_vals, sign_vals,
for mf, chif, l, m, n, sgn, exp_w, exp_tau in zip(
mf_vals, chif_vals, l_vals, m_vals, n_vals, sign_vals,
expected_w_r_vals, expected_tau_vals,
):
result_w, result_tau = gen_utils.get_qnm(chi, M, l, m, n=n, sign=sgn)
result_w, result_tau = gen_utils.get_qnm(mf, chif, l, m, n=n, sign=sgn)
assert np.isclose(result_w, exp_w, rtol=1e-8)
assert np.isclose(result_tau, exp_tau, rtol=1e-8)
Loading