From b014802b12225ff6809b6bf7c64283e4f8be26db Mon Sep 17 00:00:00 2001 From: Ravi Kumar Date: Thu, 19 Dec 2024 16:46:02 +0000 Subject: [PATCH 1/8] bug fix for xy spins --- ml4gw/waveforms/cbc/phenom_p.py | 40 +++++++++++++++++---------------- 1 file changed, 21 insertions(+), 19 deletions(-) diff --git a/ml4gw/waveforms/cbc/phenom_p.py b/ml4gw/waveforms/cbc/phenom_p.py index df91ef2f..46f30527 100644 --- a/ml4gw/waveforms/cbc/phenom_p.py +++ b/ml4gw/waveforms/cbc/phenom_p.py @@ -1,14 +1,19 @@ +""" + Based on the JAX implementation of IMRPhenomPv2 from + https://github.com/tedwards2412/ripple/blob/main/src/ripplegw/waveforms/IMRPhenomPv2.py +""" + from typing import Dict, Optional, Tuple import torch from jaxtyping import Float from torch import Tensor -from ml4gw.constants import MPC_SEC, MTSUN_SI, PI -from ml4gw.types import BatchTensor, FrequencySeries1d -from ml4gw.waveforms.conversion import rotate_y, rotate_z +from constants import MPC_SEC, MTSUN_SI, PI +from types import BatchTensor, FrequencySeries1d +from conversion import rotate_y, rotate_z -from .phenom_d import IMRPhenomD +from phenom_d import IMRPhenomD class IMRPhenomPv2(IMRPhenomD): @@ -83,7 +88,6 @@ def forward( s1x, s2x = s2x, s1x s1y, s2y = s2y, s1y s1z, s2z = s2z, s1z - ( chi1_l, chi2_l, @@ -320,11 +324,9 @@ def PhenomPOneFrequency( phic: Orbital phase at peak of the underlying non precessing model M: Total mass (Solar masses) """ - M_s = M * MTSUN_SI Mf = torch.outer(M_s, fs) fRD, _ = self.phP_get_fRD_fdamp(m1, m2, chi1, chi2, chip) - phase, _ = self.phenom_d_phase(Mf, m1, m2, eta, eta2, chi1, chi2, xi) phase = (phase.mT - (phic + PI / 4.0)).mT Amp = self.phenom_d_amp( @@ -336,10 +338,12 @@ def PhenomPOneFrequency( # phase -= 2. * phic; # line 1316 ??? hPhenom = Amp * (torch.exp(-1j * phase)) - fRDs = torch.outer( - fRD, torch.linspace(0.5, 1.5, 101, device=fRD.device) - ) - delta_fRds = torch.median(torch.diff(fRDs, axis=1), axis=1)[0] + # calculating derivative of phase with frequency following + # https://git.ligo.org/lscsoft/lalsuite/-/blame/master/lalsimulation/lib/LALSimIMRPhenomP.c?page=2#L1057 + n_fixed = 1000 + x = torch.linspace(0.8, 1.2, n_fixed, device=fRD.device) + fRDs = torch.outer(fRD, x) + delta_fRds = (1.2 * fRD - 0.8 * fRD) / (n_fixed - 1) MfRDs = torch.zeros_like(fRDs) for i in range(fRD.shape[0]): MfRDs[i, :] = torch.outer(M_s, fRDs[i, :])[i, :] @@ -350,12 +354,10 @@ def PhenomPOneFrequency( diffRDphase = (diff[:, 1:] + diff[:, :-1]) / ( 2 * delta_fRds.unsqueeze(1) ) - diffRDphase = -diffRDphase[:, 50] - # MfRD = torch.outer(M_s, fRD) - # Dphase = torch.diag( - # -self.phenom_d_phase( - # MfRD, m1, m2, eta, eta2, chi1, chi2, xi)[1] * M_s - # ).view(-1, 1) + # interpolate at x = 1, as thats the same as f = fRD + diffRDphase = -self.interpolate( + torch.tensor([1]), x[1:-1], diffRDphase + ) return hPhenom, diffRDphase # Utility functions @@ -752,9 +754,9 @@ def FinalSpin_inplane( eta2 = eta * eta # m1 > m2, the convention used in phenomD # (not the convention of internal phenomP) - mass_ratio = m1 / m2 + q_factor = m1 / M af_parallel = self.FinalSpin0815(eta, eta2, chi1_l, chi2_l) - Sperp = chip * mass_ratio * mass_ratio + Sperp = chip * q_factor * q_factor af = torch.copysign( torch.ones_like(af_parallel), af_parallel ) * torch.sqrt(Sperp * Sperp + af_parallel * af_parallel) From 4405fb02a59c088babca861f94242ccb29b46759 Mon Sep 17 00:00:00 2001 From: Ravi Kumar Date: Thu, 19 Dec 2024 16:51:10 +0000 Subject: [PATCH 2/8] run pre-commit --- ml4gw/waveforms/cbc/phenom_p.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/ml4gw/waveforms/cbc/phenom_p.py b/ml4gw/waveforms/cbc/phenom_p.py index 46f30527..8d3fc81d 100644 --- a/ml4gw/waveforms/cbc/phenom_p.py +++ b/ml4gw/waveforms/cbc/phenom_p.py @@ -1,19 +1,17 @@ """ - Based on the JAX implementation of IMRPhenomPv2 from + Based on the JAX implementation of IMRPhenomPv2 from https://github.com/tedwards2412/ripple/blob/main/src/ripplegw/waveforms/IMRPhenomPv2.py """ +from types import BatchTensor, FrequencySeries1d from typing import Dict, Optional, Tuple import torch -from jaxtyping import Float -from torch import Tensor - from constants import MPC_SEC, MTSUN_SI, PI -from types import BatchTensor, FrequencySeries1d from conversion import rotate_y, rotate_z - +from jaxtyping import Float from phenom_d import IMRPhenomD +from torch import Tensor class IMRPhenomPv2(IMRPhenomD): @@ -339,7 +337,7 @@ def PhenomPOneFrequency( hPhenom = Amp * (torch.exp(-1j * phase)) # calculating derivative of phase with frequency following - # https://git.ligo.org/lscsoft/lalsuite/-/blame/master/lalsimulation/lib/LALSimIMRPhenomP.c?page=2#L1057 + # https://git.ligo.org/lscsoft/lalsuite/-/blame/master/lalsimulation/lib/LALSimIMRPhenomP.c?page=2#L1057 # noqa: E501 n_fixed = 1000 x = torch.linspace(0.8, 1.2, n_fixed, device=fRD.device) fRDs = torch.outer(fRD, x) From c8e449b52fb1ac245017e6503578ac57345fb32e Mon Sep 17 00:00:00 2001 From: Ravi Kumar Date: Sun, 29 Dec 2024 00:10:12 +0000 Subject: [PATCH 3/8] pass fRD and fDM through to phenom_d phase and amp functions --- ml4gw/waveforms/cbc/phenom_d.py | 57 +++++++++++++++++++++++++-------- ml4gw/waveforms/cbc/phenom_p.py | 31 +++++++++++++++--- 2 files changed, 69 insertions(+), 19 deletions(-) diff --git a/ml4gw/waveforms/cbc/phenom_d.py b/ml4gw/waveforms/cbc/phenom_d.py index 9a883825..c37388d9 100644 --- a/ml4gw/waveforms/cbc/phenom_d.py +++ b/ml4gw/waveforms/cbc/phenom_d.py @@ -157,21 +157,27 @@ def phenom_d_amp( chi22, xi, distance, + fRD=None, # used for passing ringdown frequency from phenom_p + fDM=None, # used for passing damping frequency from phenom_p ): ins_amp, ins_Damp = self.phenom_d_inspiral_amp( Mf, eta, eta2, Seta, xi, chi1, chi2, chi12, chi22 ) int_amp, int_Damp = self.phenom_d_int_amp( - Mf, eta, eta2, Seta, chi1, chi2, chi12, chi22, xi + Mf, eta, eta2, Seta, chi1, chi2, chi12, chi22, xi, fRD, fDM ) mrd_amp, mrd_Damp = self.phenom_d_mrd_amp( - Mf, eta, eta2, chi1, chi2, xi + Mf, eta, eta2, chi1, chi2, xi, fRD, fDM ) gamma2 = self.gamma2_fun(eta, eta2, xi) gamma3 = self.gamma3_fun(eta, eta2, xi) - fRD, fDM = self.fring_fdamp(eta, eta2, chi1, chi2) + + # merger ringdown + if (fRD is None) or (fDM is None): + fRD, fDM = self.fring_fdamp(eta, eta2, chi1, chi2) Mf_peak = self.fmaxCalc(fRD, fDM, gamma2, gamma3) + # Geometric peak and joining frequencies Mf_peak = (torch.ones_like(Mf).mT * Mf_peak).mT Mf_join_ins = 0.014 * torch.ones_like(Mf) @@ -201,10 +207,22 @@ def phenom_d_amp( return amp, Damp def phenom_d_int_amp( - self, Mf, eta, eta2, Seta, chi1, chi2, chi12, chi22, xi + self, + Mf, + eta, + eta2, + Seta, + chi1, + chi2, + chi12, + chi22, + xi, + fRD=None, # used for passing ringdown frequency from phenom_p + fDM=None, # used for passing damping frequency from phenom_p ): # merger ringdown - fRD, fDM = self.fring_fdamp(eta, eta2, chi1, chi2) + if (fRD is None) or (fDM is None): + fRD, fDM = self.fring_fdamp(eta, eta2, chi1, chi2) # Geometric frequency definition from PhenomD header file AMP_fJoin_INS = 0.014 @@ -220,7 +238,9 @@ def phenom_d_int_amp( v1, d1 = self.phenom_d_inspiral_amp( Mf1, eta, eta2, Seta, xi, chi1, chi2, chi12, chi22 ) - v3, d2 = self.phenom_d_mrd_amp(Mf3, eta, eta2, chi1, chi2, xi) + v3, d2 = self.phenom_d_mrd_amp( + Mf3, eta, eta2, chi1, chi2, xi, fRD, fDM + ) v2 = ( torch.ones_like(Mf).mT * self.AmpIntColFitCoeff(eta, eta2, xi) ).mT @@ -239,9 +259,12 @@ def phenom_d_int_amp( ) return amp, Damp - def phenom_d_mrd_amp(self, Mf, eta, eta2, chi1, chi2, xi): + def phenom_d_mrd_amp( + self, Mf, eta, eta2, chi1, chi2, xi, fRD=None, fDM=None + ): # merger ringdown - fRD, fDM = self.fring_fdamp(eta, eta2, chi1, chi2) + if (fRD is None) or (fDM is None): + fRD, fDM = self.fring_fdamp(eta, eta2, chi1, chi2) gamma1 = self.gamma1_fun(eta, eta2, xi) gamma2 = self.gamma2_fun(eta, eta2, xi) @@ -384,17 +407,20 @@ def phenom_d_inspiral_amp( return amp, Damp - def phenom_d_phase(self, Mf, mass_1, mass_2, eta, eta2, chi1, chi2, xi): + def phenom_d_phase( + self, Mf, mass_1, mass_2, eta, eta2, chi1, chi2, xi, fRD=None, fDM=None + ): ins_phase, ins_Dphase = self.phenom_d_inspiral_phase( Mf, mass_1, mass_2, eta, eta2, xi, chi1, chi2 ) int_phase, int_Dphase = self.phenom_d_int_phase(Mf, eta, eta2, xi) mrd_phase, mrd_Dphase = self.phenom_d_mrd_phase( - Mf, eta, eta2, chi1, chi2, xi + Mf, eta, eta2, chi1, chi2, xi, fRD, fDM ) # merger ringdown - fRD, fDM = self.fring_fdamp(eta, eta2, chi1, chi2) + if (fRD is None) or (fDM is None): + fRD, fDM = self.fring_fdamp(eta, eta2, chi1, chi2) # definitions in Eq. (35) of arXiv:1508.07253 # PHI_fJoin_INS in header LALSimIMRPhenomD.h # C1 continuity at intermediate region i.e. f_1 @@ -415,7 +441,7 @@ def phenom_d_phase(self, Mf, mass_1, mass_2, eta, eta2, chi1, chi2, xi): fRDJoin, eta, eta2, xi ) mrd_phase_rd, mrd_Dphase_rd = self.phenom_d_mrd_phase( - fRDJoin, eta, eta2, chi1, chi2, xi + fRDJoin, eta, eta2, chi1, chi2, xi, fRD, fDM ) PhiIntTempVal = (int_phase_rd.mT / eta).mT + C1Int + C2Int * fRDJoin # C2MRD = int_Dphase_rd - mrd_Dphase_rd @@ -454,7 +480,9 @@ def phenom_d_phase(self, Mf, mass_1, mass_2, eta, eta2, chi1, chi2, xi): return phasing, Dphasing - def phenom_d_mrd_phase(self, Mf, eta, eta2, chi1, chi2, xi): + def phenom_d_mrd_phase( + self, Mf, eta, eta2, chi1, chi2, xi, fRD=None, fDM=None + ): alpha1 = self.alpha1Fit(eta, eta2, xi) alpha2 = self.alpha2Fit(eta, eta2, xi) alpha3 = self.alpha3Fit(eta, eta2, xi) @@ -462,7 +490,8 @@ def phenom_d_mrd_phase(self, Mf, eta, eta2, chi1, chi2, xi): alpha5 = self.alpha5Fit(eta, eta2, xi) # merger ringdown - fRD, fDM = self.fring_fdamp(eta, eta2, chi1, chi2) + if (fRD is None) or (fDM is None): + fRD, fDM = self.fring_fdamp(eta, eta2, chi1, chi2) f_minus_alpha5_fRD = (Mf.t() - alpha5 * fRD).t() # Leading 1/eta is not multiplied at this stage diff --git a/ml4gw/waveforms/cbc/phenom_p.py b/ml4gw/waveforms/cbc/phenom_p.py index 8d3fc81d..5c65548f 100644 --- a/ml4gw/waveforms/cbc/phenom_p.py +++ b/ml4gw/waveforms/cbc/phenom_p.py @@ -324,16 +324,37 @@ def PhenomPOneFrequency( """ M_s = M * MTSUN_SI Mf = torch.outer(M_s, fs) - fRD, _ = self.phP_get_fRD_fdamp(m1, m2, chi1, chi2, chip) - phase, _ = self.phenom_d_phase(Mf, m1, m2, eta, eta2, chi1, chi2, xi) + fRD, fDM = self.phP_get_fRD_fdamp(m1, m2, chi1, chi2, chip) + # pass M_s * ringdown and M_s * damping frequency to PhenomD functions + MfRD, MfDM = M_s * fRD, M_s * fDM + + phase, _ = self.phenom_d_phase( + Mf, m1, m2, eta, eta2, chi1, chi2, xi, MfRD, MfDM + ) phase = (phase.mT - (phic + PI / 4.0)).mT + # why are they subtracting 2*phic? + # https://git.ligo.org/lscsoft/lalsuite/-/blob/master/lalsimulation/lib/LALSimIMRPhenomP.c#L1316 + Amp = self.phenom_d_amp( - Mf, m1, m2, eta, eta2, Seta, chi1, chi2, chi12, chi22, xi, distance + Mf, + m1, + m2, + eta, + eta2, + Seta, + chi1, + chi2, + chi12, + chi22, + xi, + distance, + MfRD, + MfDM, )[0] Amp0 = self.get_Amp0(Mf, eta) dist_s = distance * MPC_SEC Amp = ((Amp0 * Amp).mT * (M_s**2.0) / dist_s).mT - # phase -= 2. * phic; # line 1316 ??? + hPhenom = Amp * (torch.exp(-1j * phase)) # calculating derivative of phase with frequency following @@ -346,7 +367,7 @@ def PhenomPOneFrequency( for i in range(fRD.shape[0]): MfRDs[i, :] = torch.outer(M_s, fRDs[i, :])[i, :] RD_phase = self.phenom_d_phase( - MfRDs, m1, m2, eta, eta2, chi1, chi2, xi + MfRDs, m1, m2, eta, eta2, chi1, chi2, xi, MfRD, MfDM )[0] diff = torch.diff(RD_phase, axis=1) diffRDphase = (diff[:, 1:] + diff[:, :-1]) / ( From cd0e8ea1207cc34c81fa6eb6846ef4cb872a75ca Mon Sep 17 00:00:00 2001 From: Ravi Kumar Date: Sun, 29 Dec 2024 02:17:08 +0000 Subject: [PATCH 4/8] switch to relative imports --- ml4gw/waveforms/cbc/phenom_p.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/ml4gw/waveforms/cbc/phenom_p.py b/ml4gw/waveforms/cbc/phenom_p.py index 5c65548f..6a14a32b 100644 --- a/ml4gw/waveforms/cbc/phenom_p.py +++ b/ml4gw/waveforms/cbc/phenom_p.py @@ -3,16 +3,17 @@ https://github.com/tedwards2412/ripple/blob/main/src/ripplegw/waveforms/IMRPhenomPv2.py """ -from types import BatchTensor, FrequencySeries1d from typing import Dict, Optional, Tuple import torch -from constants import MPC_SEC, MTSUN_SI, PI -from conversion import rotate_y, rotate_z from jaxtyping import Float -from phenom_d import IMRPhenomD from torch import Tensor +from ...constants import MPC_SEC, MTSUN_SI, PI +from ...types import BatchTensor, FrequencySeries1d +from ..conversion import rotate_y, rotate_z +from .phenom_d import IMRPhenomD + class IMRPhenomPv2(IMRPhenomD): def __init__(self): From 70b72b6d97da838cf9b860dd57ddd566622d863f Mon Sep 17 00:00:00 2001 From: Ravi Kumar Date: Sun, 29 Dec 2024 02:18:25 +0000 Subject: [PATCH 5/8] use chirp mass to mass ratio conversion function --- ml4gw/waveforms/cbc/phenom_p.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/ml4gw/waveforms/cbc/phenom_p.py b/ml4gw/waveforms/cbc/phenom_p.py index 6a14a32b..39c4d0cb 100644 --- a/ml4gw/waveforms/cbc/phenom_p.py +++ b/ml4gw/waveforms/cbc/phenom_p.py @@ -11,7 +11,11 @@ from ...constants import MPC_SEC, MTSUN_SI, PI from ...types import BatchTensor, FrequencySeries1d -from ..conversion import rotate_y, rotate_z +from ..conversion import ( + chirp_mass_and_mass_ratio_to_components, + rotate_y, + rotate_z, +) from .phenom_d import IMRPhenomD @@ -79,8 +83,9 @@ def forward( if tc is None: tc = torch.zeros_like(chirp_mass) - m2 = chirp_mass * (1.0 + mass_ratio) ** 0.2 / mass_ratio**0.6 - m1 = m2 * mass_ratio + m1, m2 = chirp_mass_and_mass_ratio_to_components( + chirp_mass, mass_ratio + ) # # flip m1 m2. For some reason LAL uses this convention for PhenomPv2 m1, m2 = m2, m1 From fb7661c9c88400ad0bab1b73563b7d2674821a96 Mon Sep 17 00:00:00 2001 From: Ravi Kumar Date: Sun, 29 Dec 2024 03:26:12 +0000 Subject: [PATCH 6/8] create tests for phenom_p --- tests/waveforms/cbc/test_cbc_waveforms.py | 105 +++++++++++++--------- 1 file changed, 65 insertions(+), 40 deletions(-) diff --git a/tests/waveforms/cbc/test_cbc_waveforms.py b/tests/waveforms/cbc/test_cbc_waveforms.py index b2542e96..ccaa060d 100644 --- a/tests/waveforms/cbc/test_cbc_waveforms.py +++ b/tests/waveforms/cbc/test_cbc_waveforms.py @@ -7,7 +7,10 @@ from torch.distributions import Uniform import ml4gw.waveforms as waveforms -from ml4gw.waveforms.conversion import chirp_mass_and_mass_ratio_to_components +from ml4gw.waveforms.conversion import ( + bilby_spins_to_lalsim, + chirp_mass_and_mass_ratio_to_components, +) @pytest.fixture(params=[256, 1024, 2048]) @@ -299,7 +302,6 @@ def test_phenom_d( ) hc_ml4gw = hc_ml4gw[0] - hp_ml4gw = hp_ml4gw[0] hp_lal_data = hp_lal.data.data[lal_mask] @@ -322,41 +324,64 @@ def test_phenom_d( ) -def test_phenom_p(): - chirp_mass = torch.tensor([15.0, 30.0]) - mass_ratio = torch.tensor([0.99, 0.5]) - chi1z = torch.tensor([0.0, 0.5]) - chi2z = torch.tensor([-0.1, 0.1]) - distance = torch.tensor([100.0, 1000.0]) - sample_rate = 2048 - - mass_2 = chirp_mass * (1 + mass_ratio) ** 0.2 / mass_ratio**0.6 - mass_1 = mass_2 * mass_ratio +def test_phenom_p( + chirp_mass, + mass_ratio, + distance, + phase, + sample_rate, + f_ref, + theta_jn, + phi_jl, + tilt_1, + tilt_2, + phi_12, + a_1, + a_2, +): + mass_1, mass_2 = chirp_mass_and_mass_ratio_to_components( + chirp_mass, mass_ratio + ) + ( + inclination, + chi1x, + chi1y, + chi1z, + chi2x, + chi2y, + chi2z, + ) = bilby_spins_to_lalsim( + theta_jn, + phi_jl, + tilt_1, + tilt_2, + phi_12, + a_1, + a_2, + mass_1, + mass_2, + f_ref, + phase, + ) - f_ref = 20.0 - phic = 0.0 tc = 0.0 - inclination = 0.0 - for i in range(chirp_mass.shape[0]): - m1, m2 = mass_1[i], mass_2[i] - mr = mass_ratio[i] - if m2 > m1: - m1, m2 = m2, m1 - mr = 1 / mr + # compare each waveform with lalsimulation + for i in range(len(chirp_mass)): + # construct lalinference params params = dict( - m1=m1.item() * lal.MSUN_SI, - m2=m2.item() * lal.MSUN_SI, - S1x=0, - S1y=0, + m1=mass_1[i].item() * lal.MSUN_SI, + m2=mass_2[i].item() * lal.MSUN_SI, + S1x=chi1x[i].item(), + S1y=chi1y[i].item(), S1z=chi1z[i].item(), - S2x=0, - S2y=0, + S2x=chi2x[i].item(), + S2y=chi2y[i].item(), S2z=chi2z[i].item(), distance=(distance[i].item() * u.Mpc).to("m").value, - inclination=inclination, - phiRef=phic, + inclination=inclination[i].item(), + phiRef=phase[i].item(), longAscNodes=0.0, eccentricity=0.0, meanPerAno=0.0, @@ -388,16 +413,16 @@ def test_phenom_p(): hc_ml4gw, hp_ml4gw = waveforms.IMRPhenomPv2()( torch_freqs, chirp_mass[i][None], - torch.tensor([mr]), - torch.tensor([0.0]), - torch.tensor([0.0]), + mass_ratio[i][None], + chi1x[i][None], + chi1y[i][None], chi1z[i][None], - torch.tensor([0.0]), - torch.tensor([0.0]), + chi2x[i][None], + chi2y[i][None], chi2z[i][None], distance[i][None], - torch.tensor([phic]), - torch.tensor([inclination]), + phase[i][None], + inclination[i][None], f_ref, torch.tensor([tc]), ) @@ -409,14 +434,14 @@ def test_phenom_p(): hc_lal_data = hc_lal.data.data[lal_mask] assert np.allclose( - 1e21 * hp_lal_data.real, 1e21 * hp_ml4gw.real.numpy(), atol=2e-3 + 1e21 * hp_lal_data.real, 1e21 * hp_ml4gw.real.numpy(), atol=1e-3 ) assert np.allclose( - 1e21 * hp_lal_data.imag, 1e21 * hp_ml4gw.imag.numpy(), atol=2e-3 + 1e21 * hp_lal_data.imag, 1e21 * hp_ml4gw.imag.numpy(), atol=1e-3 ) assert np.allclose( - 1e21 * hc_lal_data.real, 1e21 * hc_ml4gw.real.numpy(), atol=2e-3 + 1e21 * hc_lal_data.real, 1e21 * hc_ml4gw.real.numpy(), atol=1e-3 ) assert np.allclose( - 1e21 * hc_lal_data.imag, 1e21 * hc_ml4gw.imag.numpy(), atol=2e-3 + 1e21 * hc_lal_data.imag, 1e21 * hc_ml4gw.imag.numpy(), atol=1e-3 ) From d99253eb62cda03996a9f6a409552e3013995a83 Mon Sep 17 00:00:00 2001 From: Ravi Kumar Date: Sun, 29 Dec 2024 13:41:07 +0000 Subject: [PATCH 7/8] lower tolerance for phenom_p --- tests/waveforms/cbc/test_cbc_waveforms.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/waveforms/cbc/test_cbc_waveforms.py b/tests/waveforms/cbc/test_cbc_waveforms.py index ccaa060d..0b0c68d6 100644 --- a/tests/waveforms/cbc/test_cbc_waveforms.py +++ b/tests/waveforms/cbc/test_cbc_waveforms.py @@ -434,14 +434,14 @@ def test_phenom_p( hc_lal_data = hc_lal.data.data[lal_mask] assert np.allclose( - 1e21 * hp_lal_data.real, 1e21 * hp_ml4gw.real.numpy(), atol=1e-3 + 1e21 * hp_lal_data.real, 1e21 * hp_ml4gw.real.numpy(), atol=1e-2 ) assert np.allclose( - 1e21 * hp_lal_data.imag, 1e21 * hp_ml4gw.imag.numpy(), atol=1e-3 + 1e21 * hp_lal_data.imag, 1e21 * hp_ml4gw.imag.numpy(), atol=1e-2 ) assert np.allclose( - 1e21 * hc_lal_data.real, 1e21 * hc_ml4gw.real.numpy(), atol=1e-3 + 1e21 * hc_lal_data.real, 1e21 * hc_ml4gw.real.numpy(), atol=1e-2 ) assert np.allclose( - 1e21 * hc_lal_data.imag, 1e21 * hc_ml4gw.imag.numpy(), atol=1e-3 + 1e21 * hc_lal_data.imag, 1e21 * hc_ml4gw.imag.numpy(), atol=1e-2 ) From 971475be85f58e1b1ef87c5d2b18dabcf02d76a7 Mon Sep 17 00:00:00 2001 From: Ravi Kumar Date: Sun, 29 Dec 2024 15:25:41 +0000 Subject: [PATCH 8/8] comments and additional checks for fRD and fDM --- ml4gw/waveforms/cbc/phenom_d.py | 41 +++++++++++++++++++++++++++------ 1 file changed, 34 insertions(+), 7 deletions(-) diff --git a/ml4gw/waveforms/cbc/phenom_d.py b/ml4gw/waveforms/cbc/phenom_d.py index c37388d9..8b4911af 100644 --- a/ml4gw/waveforms/cbc/phenom_d.py +++ b/ml4gw/waveforms/cbc/phenom_d.py @@ -174,10 +174,14 @@ def phenom_d_amp( gamma3 = self.gamma3_fun(eta, eta2, xi) # merger ringdown - if (fRD is None) or (fDM is None): + if (fRD is None) != (fDM is None): + raise ValueError( + "Both fRD and fDM must either be provided or both be None" + ) + if (fRD is None) and (fDM is None): fRD, fDM = self.fring_fdamp(eta, eta2, chi1, chi2) - Mf_peak = self.fmaxCalc(fRD, fDM, gamma2, gamma3) + Mf_peak = self.fmaxCalc(fRD, fDM, gamma2, gamma3) # Geometric peak and joining frequencies Mf_peak = (torch.ones_like(Mf).mT * Mf_peak).mT Mf_join_ins = 0.014 * torch.ones_like(Mf) @@ -221,8 +225,13 @@ def phenom_d_int_amp( fDM=None, # used for passing damping frequency from phenom_p ): # merger ringdown - if (fRD is None) or (fDM is None): + if (fRD is None) != (fDM is None): + raise ValueError( + "Both fRD and fDM must either be provided or both be None" + ) + if (fRD is None) and (fDM is None): fRD, fDM = self.fring_fdamp(eta, eta2, chi1, chi2) + # Geometric frequency definition from PhenomD header file AMP_fJoin_INS = 0.014 @@ -259,11 +268,17 @@ def phenom_d_int_amp( ) return amp, Damp + # fRD and fDM are to be passed for generating phenom_p waveforms + # and remain None for phenom_d def phenom_d_mrd_amp( self, Mf, eta, eta2, chi1, chi2, xi, fRD=None, fDM=None ): # merger ringdown - if (fRD is None) or (fDM is None): + if (fRD is None) != (fDM is None): + raise ValueError( + "Both fRD and fDM must either be provided or both be None" + ) + if (fRD is None) and (fDM is None): fRD, fDM = self.fring_fdamp(eta, eta2, chi1, chi2) gamma1 = self.gamma1_fun(eta, eta2, xi) @@ -407,6 +422,8 @@ def phenom_d_inspiral_amp( return amp, Damp + # fRD and fDM are to be passed for generating phenom_p waveforms + # and remain None for phenom_d def phenom_d_phase( self, Mf, mass_1, mass_2, eta, eta2, chi1, chi2, xi, fRD=None, fDM=None ): @@ -419,7 +436,11 @@ def phenom_d_phase( ) # merger ringdown - if (fRD is None) or (fDM is None): + if (fRD is None) != (fDM is None): + raise ValueError( + "Both fRD and fDM must either be provided or both be None" + ) + if (fRD is None) and (fDM is None): fRD, fDM = self.fring_fdamp(eta, eta2, chi1, chi2) # definitions in Eq. (35) of arXiv:1508.07253 # PHI_fJoin_INS in header LALSimIMRPhenomD.h @@ -480,6 +501,8 @@ def phenom_d_phase( return phasing, Dphasing + # fRD and fDM are to be passed for generating phenom_p waveforms + # and remain None for phenom_d def phenom_d_mrd_phase( self, Mf, eta, eta2, chi1, chi2, xi, fRD=None, fDM=None ): @@ -490,10 +513,14 @@ def phenom_d_mrd_phase( alpha5 = self.alpha5Fit(eta, eta2, xi) # merger ringdown - if (fRD is None) or (fDM is None): + if (fRD is None) != (fDM is None): + raise ValueError( + "Both fRD and fDM must either be provided or both be None" + ) + if (fRD is None) and (fDM is None): fRD, fDM = self.fring_fdamp(eta, eta2, chi1, chi2) - f_minus_alpha5_fRD = (Mf.t() - alpha5 * fRD).t() + f_minus_alpha5_fRD = (Mf.t() - alpha5 * fRD).t() # Leading 1/eta is not multiplied at this stage mrd_phasing = (Mf.t() * alpha1).t() mrd_phasing -= (1 / Mf.t() * alpha2).t()