Figure 2B

Here we describe how to reproduce Figure 2B

Step 1: Run the PyCBC Inference Bayesian analysis code

We use pycbc_inference to make the data for this figure, but it first needs to know about our stretched/squeezed waveform as described in the manuscript.

To do this we add the following function into pycbc/waveform/waveform.py

from scipy import interpolate
MPC_S_SI = lal.PC_SI * 1E6 / lal.C_SI
from pycbc.waveform.utils import td_taper
def get_dopplered_tf2_waveform(dtype=numpy.complex128, return_hc=True,
                               is_sequence=False,
                               **params):
    cparams = copy.deepcopy(params)
    cparams['approximant'] = 'TaylorF2'
    if is_sequence:
        hp, hc = get_fd_waveform_sequence(**cparams)
    else:
        hp, hc = get_fd_waveform(**cparams)
    non_gr_fac_n = params['non_gr_fac_n']
    non_gr_fac_a_at_30 = params['non_gr_fac_a_at_30']
    non_gr_fac_a = non_gr_fac_a_at_30 * 30**non_gr_fac_n
    tsize = int(1.0 / params['delta_t'] /  params['delta_f'])
    fsize = tsize // 2 + 1
    hp.resize(fsize)
    td_wav = hp.to_timeseries()
    td_wav.roll(- int(len(td_wav) * 0.2))
    sts_shifted = numpy.arange(len(td_wav), dtype=float) * td_wav.delta_t + float(td_wav._epoch)
    sts_shifted += int(len(td_wav) * 0.2) * td_wav.delta_t
    sts_shifted_red = sts_shifted[sts_shifted<0]
    mc = conversions.mchirp_from_mass1_mass2(params['mass1'],
                                             params['mass2'])
    f_t_basic = (256. / 5. * (mc * lal.MTSUN_SI)**(5./3.) * numpy.pi**(8./3.)
                     * (-sts_shifted_red))**(-3./8.)
    new_times = (sts_shifted_red - non_gr_fac_a * 40 * MPC_S_SI
                     * f_t_basic**(-non_gr_fac_n))
    new_times = numpy.concatenate([new_times,sts_shifted[sts_shifted >= 0]])
    interfunc = interpolate.interp1d(new_times, td_wav.data,
                                     bounds_error=False,
                                     fill_value='extrapolate')
    new_data = interfunc(sts_shifted)
    td_wav.data[:] = new_data[:]
    td_wav = td_taper(td_wav, -len(td_wav) * td_wav.delta_t, -0.9 * len(td_wav) * td_wav.delta_t, side='left')
    td_wav = td_taper(td_wav, -0.1*len(td_wav) * td_wav.delta_t, 0, side='right')
    td_wav.roll(+ int(len(td_wav) * 0.2))

    hp = td_wav.to_frequencyseries()
    hp = hp.astype(dtype)
    
    if return_hc:
        hc.resize(fsize)
        td_wav = hc.to_timeseries()
        td_wav.roll(- int(len(td_wav) * 0.2))
        interfunc = interpolate.interp1d(new_times, td_wav.data, 
                                         bounds_error=False, 
                                         fill_value='extrapolate')
        new_data = interfunc(sts_shifted)
        td_wav.data[:] = new_data[:]
        td_wav = td_taper(td_wav, -len(td_wav) * td_wav.delta_t, -0.9 * len(td_wav) * td_wav.delta_t, side='left')
        td_wav = td_taper(td_wav, -0.1*len(td_wav) * td_wav.delta_t, 0, side='right')
        td_wav.roll(+ int(len(td_wav) * 0.2))

        hc = td_wav.to_frequencyseries()
        hc = hc.astype(dtype)
    else:
        hc = None

    return hp, hc

and declare it as a waveform by adding:

    if apx == 'TaylorF2':
        apx_new = 'TaylorF2_DOPPLER'
        cpu_fd[apx_new] = get_dopplered_tf2_waveform

at this line

https://github.com/gwastro/pycbc/blob/v2.0.3/pycbc/waveform/waveform.py#L1003

This will correspond to superluminal cgw. We simulate subluminal cgw by changing the line:

    new_times = (sts_shifted_red - non_gr_fac_a * 40 * MPC_S_SI
                     * f_t_basic**(-non_gr_fac_n))

    new_times = (sts_shifted_red + non_gr_fac_a * 40 * MPC_S_SI
                     * f_t_basic**(-non_gr_fac_n))

(doubtless this would need to be more elegant if being run on production LVK/LISA data!).

Then we need a configuration file, ours can be downloaded here:

https://icg-gravwaves.github.io/probing_speed_of_gravity/Figure_2B/inference_emcee.ini

Then we need to run pycbc_inference. We did this using MPI according to:

mpiexec -n 50 pycbc_inference --config-files inference_emcee.ini --output-file results_superluminal.hdf --verbose --use-mpi

This produces results_superluminal.hdf and results_subluminal.hdf, which we upload to the repository here:

https://icg-gravwaves.github.io/probing_speed_of_gravity/Figure_2B/results_superluminal.hdf

https://icg-gravwaves.github.io/probing_speed_of_gravity/Figure_2B/results_subluminal.hdf

Note that this is slow (takes weeks on 50 cores). The waveform generator should be made faster is using this in the future.

Step 2: Make the plot

We then make the plot in a few stages. We first extract the necessary data from the HDF files. This is done by downloaded the edited executable:

https://icg-gravwaves.github.io/probing_speed_of_gravity/Figure_2B/pycbc_inference_plot_posterior

and then running

./pycbc_inference_plot_posterior --input-file results_superluminal.hdf --plot-scatter --output-file scatter_2param_mk2.pdf --z-arg snr --parameters non_gr_fac_n "log10(non_gr_fac_a_at_30)"
mv figure2_data.npz figure2_data_superluminal.npz

./pycbc_inference_plot_posterior --input-file results_subluminal.hdf --plot-scatter --output-file scatter_2param_mk2.pdf --z-arg snr --parameters non_gr_fac_n "log10(non_gr_fac_a_at_30)"

mv figure2_data.npz figure2_data_subluminal.npz

to generate two .npz files. These files are also available here:

https://icg-gravwaves.github.io/probing_speed_of_gravity/Figure_2B/figure2_data_subluminal.npz

https://icg-gravwaves.github.io/probing_speed_of_gravity/Figure_2B/figure2_data_superluminal.npz

Finally the plot is generated from this be running the script

https://icg-gravwaves.github.io/probing_speed_of_gravity/Figure_2B/make_plot.py

Step 3: Also make the video

The make_plot.py command also makes the still frames used to make the animation. We combined this into an mp4 using ffmpeg:

ffmpeg -r 1 -f image2 -i figure2b_video_panel_%d.png -vcodec libx264 -crf 25  -pix_fmt yuv420p -vf pad="width=ceil(iw/2)*2:height=ceil(ih/2)*2" -start_number 0 figure2b.mp4