Error while generating noise_from_psd

[lpathak97]

I am trying to generate the time-domain data and using advanced LIGO and virgo psds. Here is a code which can be used to reproduce the error:

def generate_sim_data(ifos, f_lower, sampling_frequency, cosmology, duration_offset=1, seed=10, **inj_params):

    hp, hc = get_td_waveform(approximant=inj_params['approximant'],
                         mass1=(1 + inj_params['z'])*inj_params['mass1'],
                         mass2=(1 + inj_params['z'])*inj_params['mass2'],
                         spin1z=inj_params['spin1z'],
                         spin2z=inj_params['spin2z'],
                         inclination=inj_params['inclination'],
                         coa_phase=inj_params['coa_phase'],
                         distance=cosmology.luminosity_distance(inj_params['z']).value,
                         delta_t=1.0/sampling_frequency,
                         f_lower=f_lower)

    det = {}
    for ifo in ifos:
        det[ifo] = Detector(ifo)

    end_time = inj_params['tc'] 
    declination = inj_params['dec']
    right_ascension = inj_params['ra']
    polarization = inj_params['polarization']
    hp.start_time += end_time
    hc.start_time += end_time

    signal = {}
    for ifo in ifos:
        
        signal[ifo] = det[ifo].project_wave(hp, hc, right_ascension, declination, polarization)

    delta_t = 1.0 / sampling_frequency

    noise = {}
    psd = {}
    data = {}

    mf_snr = {}
    network_snr = 0

    print(duration_offset)

    for ifo in ifos:

        signal[ifo] = det[ifo].project_wave(hp, hc, right_ascension, declination, polarization)
        duration = signal[ifo].duration + duration_offset
        tsamples = int(duration / delta_t)

        delta_f = 1.0 / duration
        flen = int(sampling_frequency / delta_f) + 1

        if(ifo=='L1' or ifo == 'H1'):

            psd[ifo] = pycbc.psd.from_txt('AplusDesign.txt', flen, delta_f, f_lower, is_asd_file=True)

        else:

            psd[ifo] = pycbc.psd.from_txt('avirgo_O5low_NEW.txt', flen, delta_f, f_lower, is_asd_file=True)

        noise[ifo] = pycbc.noise.noise_from_psd(tsamples, delta_t, psd[ifo], seed=seed)
        noise[ifo] = TimeSeries(noise[ifo], epoch=signal[ifo]._epoch)
        signal[ifo].resize(len(noise[ifo]))
        data[ifo] = signal[ifo] + noise[ifo]

        hp_temp, hc_temp = get_td_waveform(approximant=inj_params['approximant'],
                            mass1=(1 + inj_params['z'])*inj_params['mass1'],
                            mass2=(1 + inj_params['z'])*inj_params['mass2'],
                            spin1z=inj_params['spin1z'],
                            spin2z=inj_params['spin2z'],
                            delta_t=1.0/sampling_frequency,
                            f_lower=f_lower)
        
        hp_temp.resize(len(data[ifo]))
        mf_snr[ifo] = abs(pycbc.filter.matched_filter(hp_temp, data[ifo], psd=psd[ifo], low_frequency_cutoff=f_lower))
        network_snr += abs(mf_snr[ifo]).max()**2
        
    network_snr = np.sqrt(network_snr)

    return signal, noise, data, psd

H0 = 67.8
Om0 = 0.308

cosmology = FlatLambdaCDM(H0=H0 * u.km / u.s / u.Mpc, Om0=Om0)

inj_params = {'approximant' : 'SpinTaylorT4', 'mass1' : 1.4, 'mass2' : 1.4, 'spin1z' : 0.0, 'spin2z' : 0.0, 'inclination' : np.pi/4, \
                'ra': np.pi/4, 'dec': np.pi/4, 'polarization': 0, 'coa_phase': 0, 'tc': 1e7, 'z': 0.01}

ifos = ['L1', 'H1', 'V1']
f_lower = 10
sampling_frequency = 4096
duration_offset = 0
signal, noise, data, psd = generate_sim_data(ifos, f_lower, sampling_frequency, cosmology, duration_offset=duration_offset, seed=10, **inj_params)

When I run it for f_lower = 10 Hz, I get the following error:

XLAL Error - XLALSimNoise (LALSimNoise.c:162): Invalid argument
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[19], line 13
     11 sampling_frequency = 4096
     12 duration_offset = 0
---> 13 signal, noise, data, psd = generate_sim_data(ifos, f_lower, sampling_frequency, cosmology, duration_offset=duration_offset, seed=10, **inj_params)

Cell In[6], line 58
     54 else:
     56     psd[ifo] = pycbc.psd.from_txt('avirgo_O5low_NEW.txt', flen, delta_f, f_lower, is_asd_file=True)
---> 58 noise[ifo] = pycbc.noise.noise_from_psd(tsamples, delta_t, psd[ifo], seed=seed)
     59 noise[ifo] = TimeSeries(noise[ifo], epoch=signal[ifo]._epoch)
     60 signal[ifo].resize(len(noise[ifo]))

File ~/miniconda3/envs/gwsim_env/lib/python3.10/site-packages/pycbc/noise/gaussian.py:121, in noise_from_psd(length, delta_t, psd, seed)
    118 segment = TimeSeries(zeros(N), delta_t=delta_t).lal()
    119 length_generated = 0
--> 121 lalsimulation.SimNoise(segment, 0, psd, randomness)
    122 while (length_generated < length):
    123     if (length_generated + stride) < length:

RuntimeError: Invalid argument

It happens for f_lower = 10.1 Hz as well. Can anyone please have a look at the code and let me know if it is expected or point to a possible bug in the above code?

Thanks

[rahuldhurkunde]

@lpathak97 The way you are importingpycbc.psd might be conflicting with your psd dictionary which is then passed as an invalid argument to the noise_from_psd. Below is a minimally working example:

import pycbc 
from pycbc import waveform, detector, noise
import pycbc.psd as PSD
import matplotlib.pyplot as plt

sampling_frequency = 4096
duration = 1000.0
delta_t = 1.0 / sampling_frequency
delta_f = 1.0 / duration
f_lower = 10
flen = int(sampling_frequency / delta_f) + 1
tsamples = int(duration / delta_t)

psd = PSD.analytical.aLIGOLateLowSensitivityP1200087(flen, delta_f, f_lower)
noise = pycbc.noise.noise_from_psd(tsamples, delta_t, psd, seed=1234)

plt.plot(noise)