Time of coalescence incorrect for time-domain waveforms #503

Closed bilby-bot closed 4 weeks ago

In GitLab by @git.ligo:serguei.ossokine on Dec 16, 2019, 12:16

This is an issue I encountered when running with parallel bilby but I think the underlying cause is in Bilby itself. If one uses a time-domain model for analysis, the coalescence time (i.e. geocentric end time) can be off by almost 100 ms. I have tested this by injecting an SEOBNRv4 signal and recovering with SEOBNRv4 and SEOBNRv4_ROM. While the ROM result is close to the true value, the time-domain is way off. With @git.ligo:roberto.cotesta we have looked at the difference in handling the time-domain vs frequency-domain waveforms. One thing that we think is missing is that XLALSimInspiralFD does not take into account the fact that the waveform does not start at 0, but rather at a negative time such that the peak of the waveform occurs at 0. Thus, we think there is an additional rotation in the Fourier domain that should be applied. I tried a preliminary version of the code and obtained the following results for the SEOBNRv4 injection:

t_geo

As one can see, the geocentric end time posterior is now much closer to the injected value than before. However, it's still not centred on the true value. Thus I think there might be something else missing.

Could the core developers take a look?

In GitLab by @git.ligo:gregory.ashton on Dec 16, 2019, 21:53

@git.ligo:serguei.ossokine is this related to #443? I.e. is it the same issue, or different?

I'm currently struggling to understand if this is a waveforms interface issue or a bilby issue. I personally can't see a "bug" right now in bilby with how we use lalsimulation. I'm going to email the waveform and PE groups to try and get some perspective.

In GitLab by @git.ligo:serguei.ossokine on Dec 16, 2019, 22:23

Hi Greg, apologies for not being clearer. I think this is more of a waveform interface issue. This issue is related but not the same as #443 . Both issue stem from using XLALSimInspiralFD. In the present case, I think this is what causes the difference in the time of coalescence for time-domain waveforms. Now lalinference calls XLALChooseTDWaveform directly and then does all the conditioning, padding, FFTing and time-shifting internally. I think it would be good to take a good look at what they do, but unfortunately I didn't have enough time to parse it.

In GitLab by @git.ligo:rory-smith on Dec 16, 2019, 22:55

For what it's worth, my 2¢ are that all waveforms (FD and TD) called via XLALSimInspiralFD should produce similar outputs (meaning they have a common reference time). Given that the FD waveforms do have a common reference time, what is the rationale for not having one with TD waveforms? My concern with considering this an interface-issue is that it falls on, e.g., PE-code developers to figure out what convention particular TD waveform families have, which can easily lead to errors creeping in. Because of this, I think we should consider modifying XLALSimInspiralFD so that all of the TD waveforms are aligned in time in a way that's consistent with the FD waveforms.

In GitLab by @git.ligo:jolien-creighton on Dec 16, 2019, 23:02

I'm pretty sure that XLALSimInspiralFD does not assume that waveforms start at 0. In fact, it is the convention that 0 is roughly the coalescence / peak time. Is the problem actually within the ROM? Also, please check to make sure that you use the .epoch of the returned timeseries as this can be shifted in order to compensate with phase rotation in the frequency domain that is needed to keep the waveform from wrapping in time.

In GitLab by @git.ligo:rory-smith on Dec 16, 2019, 23:15

Interesting, thanks @git.ligo:jolien-creighton ! It would be good if someone more familiar with the bilby likelihood could confirm if .epoch is used to shift the signals, because there's potentially a v. easy fix for the TD-waveform likelihoods here.

In GitLab by @git.ligo:gregory.ashton on Dec 17, 2019, 00:14

@git.ligo:rory-smith the wrapping of all lalsim interfaces is done in the utils: https://git.ligo.org/lscsoft/bilby/blob/master/bilby/gw/utils.py#L793

I don't see anywhere that epoch is used

In GitLab by @git.ligo:rory-smith on Dec 17, 2019, 00:44

Thanks @git.ligo:gregory.ashton I would expect that epoch would be used in the likelihood itself, as it's supposed to be used to shift the waveform to whatever epoch happens to be in the time series (I recall something similar is done in LI). Could you point me to the relevant line(s) in the likelihood function where the time shifting is done?

In GitLab by @git.ligo:gregory.ashton on Dec 17, 2019, 01:03

I believe the time-shifting happens here: https://git.ligo.org/lscsoft/bilby/blob/master/bilby/gw/detector/interferometer.py#L316

In GitLab by @git.ligo:rory-smith on Dec 17, 2019, 01:37

Thanks, @git.ligo:gregory.ashton I think the issue is on exactly line 316. I'd bet that

dt = dt_geocent + time_shift

Just needs to be changed to something like

dt = dt_geocent + time_shift ± output_of_XLALSimInspiralFD.epoch

The ± sign being there because I can never remember which way the shift is applied. Do you know a convenient way of adding output_of_XLALSimInspiralFD.epoch to the interferometer class?

In GitLab by @git.ligo:rory-smith on Dec 17, 2019, 05:30

I wrote a simple script below that computes the epoch for a waveform computed using SimInspiralFD and SimInspiralChooseFDWaveform. I found that the the epoch is indeed different between the two waveform interfaces. Moreover, the definition of epoch is consistent for all FD waveforms generated using SimInspiralChooseFDWaveform (but not SimInspiralFD). For reference, the standard definition in SimInspiralChooseFDWaveform is epoch=-T where T is the data/signal duration. In light of this, there's a time offset that needs to be applied to all waveforms from SimInspiralFD which is just dt = T - abs(waveform.epoch) By applying this offset, the definition of epoch should be the same between SimInspiralFD and SimInspiralChooseFDWaveform

import lal
import lalsimulation as lalsim

longAscNodes = 0
eccentricity = 0 # Assume circular orbit
meanPerAno = 0

WFdict = lal.CreateDict()

lalsim.SimInspiralWaveformParamsInsertFrameAxis(WFdict, 2)
lalsim.SimInspiralWaveformParamsInsertPNSpinOrder(WFdict, -1)
lalsim.SimInspiralWaveformParamsInsertPNTidalOrder(WFdict, -1)
lalsim.SimInspiralWaveformParamsInsertPNPhaseOrder(WFdict, -1)
lalsim.SimInspiralWaveformParamsInsertPNAmplitudeOrder(WFdict, -1)
lalsim.SimInspiralWaveformParamsInsertTidalLambda1(WFdict, 0.0)
lalsim.SimInspiralWaveformParamsInsertTidalLambda2(WFdict, 0.0)

s1 = [0.0, 0, 0.0]
s2 = [0.0, 0.0, 0.0]
dist = 500 # Distance in Mpc
iota = np.pi*0.4  
phi_c = 0. # Overall phase
T = 4.
deltaF = 1./T 
f_high = 4096. 
Fs = 2*f_high
deltaT = 1./Fs
f_ref = 40. # Reference frequency at which spins are calculated
f_low = 20.

approx = lalsim.IMRPhenomD # Approximation scheme in time domain

m1=30
m2=30

hplus_INSPIRAL_FD, hcross_INSPIRAL_FD = lalsim.SimInspiralFD(m1 * lal.MSUN_SI, m2 * lal.MSUN_SI,\
                                                          s1[0], s1[1], s1[2],\
                                                          s2[0], s2[1], s2[2],\
                                                          dist * 1e6 * lal.PC_SI, iota, phi_c,\
                                                          longAscNodes, eccentricity, meanPerAno,
                                                          deltaF, f_low, f_high, f_ref,\
                                                          WFdict, approx)

hplus_CHOOSE_FD, hcross_CHOOSE_FD = lalsim.SimInspiralChooseFDWaveform(m1 * lal.MSUN_SI, m2 * lal.MSUN_SI,\
                                                          0, 0., 0,\
                                                          0, 0, 0,\
                                                          dist * 1e6 * lal.PC_SI, iota, phi_c,\
                                                          longAscNodes, eccentricity, meanPerAno,\
                                                          deltaF, f_low, f_high, f_ref, WFdict, approx)
print(hplus_CHOOSE_FD.epoch)
print(hplus_INSPIRAL_FD.epoch)

In GitLab by @git.ligo:harald.pfeiffer on Dec 17, 2019, 22:12

For info: There are some time-domain waveforms where "t=0" does not at all correspond to merger time. Specifically, the surrogate models NRsur7dq2, NRsur7dq4 set t=0 at the start of the waveforms, some 30-40 GW cycles before merger.

In GitLab by @git.ligo:jolien-creighton on Dec 17, 2019, 22:17

Incidentally, the correct way to do this IMO is to use the LALSimulation routines XLALSimDetectorStrainREAL8TimeSeries to generate the strain on a detector and XLALSimAddInjectionREAL8TimeSeries to add the strain to the data (I guess for PE purposes, a routine that subtracted the strain from the data would be in order...). These are designed to work with the other routines in LALSimulation. The first one additionally takes into account Doppler shifting due to the Earth's rotation (which affects long duration signals) as well as the non-long-wavelength-limit effects on the response pattern at high frequencies (not too important for CBCs in aLIGO but will become more important especially for IFOs with longer arms). The second one deals with the subsample interpolation.

See also the stand-alone programs lalsim-inspiral, lalsim-detector-strain, lalsim-inject, and specifically the example given in the latter.

In GitLab by @git.ligo:rory-smith on Dec 18, 2019, 09:09

Hi @git.ligo:serguei.ossokine @git.ligo:gregory.ashton I implemented my patch and performed a run on with SEOBNRv4 on GW150914. The tc recovery using my patch are consistent with the IMRPhenomPv2 results. In contrast, the master branch of the code shows a significant bias the tc recovery:

The full set of comparisons plots are here on the summarypage : https://ldas-jobs.ligo.caltech.edu/~rory.smith/pbilby_review/150914/pb_li_comparison_with_SEOBNRv4_2/html/Comparison_geocent_time.html

In GitLab by @git.ligo:gregory.ashton on Jan 14, 2020, 23:38

@git.ligo:rory-smith is the point that pbilby_seobnrv4_fix is closer to the IMR results? Is it not concerning that it still seems to be shifted by a small fraction?

In GitLab by @git.ligo:rory-smith on Jan 15, 2020, 07:27

Not really. If you look at a typical event, the spread on tc from different waveform families is quite broad, e.g., https://ldas-jobs.ligo.caltech.edu/~mpitkin/projects/O3/S190521g/offline/html/Comparison_geocent_time.html

In GitLab by @git.ligo:serguei.ossokine on Jan 15, 2020, 10:31

@git.ligo:rory-smith @git.ligo:gregory.ashton I think the most convincing thing would be to run SEOBNRv4 with lalinference on this event. Then we can know for sure. @git.ligo:rory-smith could you set up such a run?

In GitLab by @git.ligo:rory-smith on Jan 20, 2020, 11:29

@git.ligo:serguei.ossokine : @git.ligo:roberto.cotesta has a nice demonstration using one of the NR surrogates which shows that the LI and pbilby tc posteriors now agree.

In GitLab by @git.ligo:roberto.cotesta on Feb 3, 2020, 08:48

Hi all,

I've used the NRHybSurr3dq8 both with LI MCMC and with pbilby and I get a very consistent tc. This is with Rory's patch. combined_1d_posterior_geocent_time

I think that we can close this issue if this patch is merged to master.

The PE run is here. Unfortunately the priors for the other parameters are a bit different so we cannot really compare the two runs for the other parameters.

In GitLab by @git.ligo:roberto.cotesta on Feb 21, 2020, 23:18

@git.ligo:gregory.ashton can we close this issue and merge the commits? It is important for the exceptional events and we have already checked that is working properly

In GitLab by @git.ligo:gregory.ashton on Feb 25, 2020, 23:52

@git.ligo:roberto.cotesta this was fixed in !736 and has been merged to master. Closing.

In GitLab by @git.ligo:gregory.ashton on Feb 25, 2020, 23:52

closed

In GitLab by @git.ligo:roberto.cotesta on Feb 25, 2020, 23:54

@git.ligo:gregory.ashton thanks!

In GitLab by @git.ligo:roberto.cotesta on Feb 27, 2020, 20:26

mentioned in commit b3cb81f4ee060a372a62aa4b01f749ea2864934b

bilby-dev / bilby

Time of coalescence incorrect for time-domain waveforms #503