pycbc_generate_hwinj

[albertoaguilarUIB]

Hi all, I'm running PyCBC=1.18.3 and 2.0.3 and trying to generate an injection and insert it into our original frames. I run this command:

pycbc_generate_hwinj \
--instruments H1 L1 V1 \
--psd-estimation median \
--frame-files H1:GW200129_065458_H1_16kHz_32seconds.gwf L1:GW200129_065458_L1_16kHz_32seconds.gwf V1:GW200129_065458_V1_16kHz_32seconds.gwf \
--approximant IMRPhenomTHM \
--order pseudoFourPN  \
--waveform-low-frequency-cutoff 20. \
--mass1 39.639  \
--mass2 34.055  \
--spin1z 0.208  \
--spin2z 0.242  \
--inclination 0.0  \
--polarization 0.0  \
--ra 5.540212805434207  \
--dec 0.13678653139409228 \
--taper TAPER_START \
--network-snr 25 \
--low-frequency-cutoff 20.0 \
--high-frequency-cutoff 1024.0 \
--sample-rate H1:16384. L1:16384. V1:16384. \
--geocentric-end-time 1264316116 \
--gps-start-time 1264316101 \
--gps-end-time   1264316133 \
--strain-high-pass 20. \
--psd-model H1:aLIGOZeroDetHighPower L1:aLIGOZeroDetHighPower V1:AdvVirgo\
--psd-output H1:'./H1-PSD.txt' L1:'./L1-PSD.txt' V1:'./V1-PSD.txt' \
--channel-name H1:GWOSC-16KHZ_R1_STRAIN L1:GWOSC-16KHZ_R1_STRAIN V1:GWOSC-16KHZ_R1_STRAIN

In the first case (1.18.3) I get the following error:

2022-07-20 20:29:48,388 : Creating XML file
2022-07-20 20:29:48,435 : Generating waveform at 40.000Mpc beginning at 20.000Hz for SNR calculation
2022-07-20 20:29:48,460 : Reading Frames
XLAL Error - XLALFrStreamSeek (LALFrStream.c:564): Invalid time
XLAL Error - XLALFrStreamReadREAL8TimeSeries (LALFrStreamReadTS_source.c:211): Internal function call failed: Invalid time
Traceback (most recent call last):
  File "/mnt/data/home/maite/py3_env/bin/pycbc_generate_hwinj", line 386, in <module>
    dyn_range_fac=DYN_RANGE_FAC)
  File "/mnt/data/home/maite/py3_env/lib/python3.6/site-packages/pycbc/strain/strain.py", line 436, in from_cli_multi_ifos
    inj_filter_rejector_dict[ifo], **kwargs)
  File "/mnt/data/home/maite/py3_env/lib/python3.6/site-packages/pycbc/strain/strain.py", line 425, in from_cli_single_ifo
    inj_filter_rejector=inj_filter_rejector, **kwargs)
  File "/mnt/data/home/maite/py3_env/lib/python3.6/site-packages/pycbc/strain/strain.py", line 227, in from_cli
    sieve=sieve)
  File "/mnt/data/home/maite/py3_env/lib/python3.6/site-packages/pycbc/frame/frame.py", line 255, in read_frame
    return _read_channel(channels, stream, start_time, duration)
  File "/mnt/data/home/maite/py3_env/lib/python3.6/site-packages/pycbc/frame/frame.py", line 82, in _read_channel
    data = read_func(stream, channel, start, duration, 0)
RuntimeError: Internal function call failed: Invalid time

And in GRAVITY-igvn environment (pycbc v2.0.3) I get the following one:

2022-07-20 21:15:43,058 : Creating XML file
Traceback (most recent call last):
  File "/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn/bin/pycbc_generate_hwinj", line 277, in <module>
    proc_id = utils.process.register_to_xmldoc(
  File "/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py39-20220707/lib/python3.9/site-packages/ligo/lw/utils/process.py", line 100, in register_to_xmldoc
    process = proctable.RowType.initialized(program = program, process_id = proctable.get_next_id(), **kwargs)
  File "/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py39-20220707/lib/python3.9/site-packages/ligo/lw/lsctables.py", line 586, in initialized
    cvs_entry_time = lal.UTCToGPS(time.strptime(cvs_entry_time, "%Y-%m-%d %H:%M:%S +0000")) if cvs_entry_time is not None else None,
  File "/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py39-20220707/lib/python3.9/_strptime.py", line 562, in _strptime_time
    tt = _strptime(data_string, format)[0]
  File "/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py39-20220707/lib/python3.9/_strptime.py", line 349, in _strptime
    raise ValueError("time data %r does not match format %r" %
ValueError: time data '' does not match format '%Y-%m-%d %H:%M:%S +0000'

Any help would be appreciated! TIA,
Jose Alberto

[spxiwh]

Hi @albertoaguilarUIB apologies for the slow response!

The issue here is that the command you are running (in 1.18.3) is trying to read in more than 32 seconds of data. The reason for this is that normally the user is doing something like a high-pass filter (as you are doing with --strain-high-pass 20, and so we try to read in an additional 8 seconds of data at the start and the end, which is then corrupted by the high-pass filter, and is stripped away, so you have 32 seconds of usable data after this). However, with only 32 second frame files, this doesn't work.

You could download the longer frame files from GWOSC, or you could add:

--pad-data 0

to the command. This will ensure that only 32s is read in, but the start and end will be a little corrupted by the high-pass filter.

You also have both --psd-estimation and psd-model options. These are mutually exclusive. I tested with:

--psd-estimation median \
--psd-segment-length 2 \
--psd-segment-stride 1 \

and things worked.

There has been some broken XML code in the more recent version (because of changes in the underlying library). The particular error you report was fixed in 2.0.4/2.0.5, but this code still didn't run on master. I'll be putting in a pull request soon to fix this for 2.0.6

[albertoaguilarUIB]

Thanks @spxiwh,

The pad_data could be the key to continue using my 32 seconds frames, I didn't realized that I needed those extra 8 seconds. And regarding the psd options I already saw it and select only the second one, the psd_model, but can maybe change into the first one to extract the psd from the frames.

I'm running my tests now in igwn-py39 environment in CIT with v.2.0.3, generating fake strain (aLIGOZeroDetHighPower PSD for H1, L1 and V1) and making/injecting zero-noise time series with signals using IMRPhenomTHM/TPHM approximants.

I have another question regarding the approximants and other arguments, for the IMRPhenomTHM approximant I'm passing the mode_array argument in '22 21 33' or [[2,2], [2,1], [3,3]] formats but not sure if get_td_waveform() function generates and add all of them (and also maybe the negative ones by symmetry). Also wonder if we can pass some other flags/arguments for the precessing approximants, e.g. IMRPhemonTPHM, like PhenomXPrecVersion or PhenomXPFinalSpinMod precession flags.

I'm just passing the mode_array but in the lalsimulation documentation I read:

For selecting a particular list of modes to be returned or to be employed in the polarisations construction, the user can follow the usual procedure:
*Create a mode array object with lalsimulation.SimInspiralCreateModeArray
*Activate the desired modes with lalsim.SimInspiralModeArrayActivateMode
*Insert the mode array into a LAL dictionary with lalsim.SimInspiralWaveformParamsInsertModeArray
*Pass the LAL ditionary to ChooseTDWaveform or ChooseTDModes. For a user specified mode array, only the implemented modes in the model will be computed.

I attach two plots with the H1 injected signals using different mode_arrays. (I am using THMv2 instead of THM in these examples, but they have the same arguments).

Thanks.

[spxiwh]

Hi @albertoaguilarUIB, it is the responsibility of PyCBC to ensure that if you send a mode_array to our waveform routines it will be correctly sent to lalsimulation. Indeed if you see here:

https://github.com/gwastro/pycbc/blob/master/pycbc/waveform/waveform.py#L130

you can see us doing points 1-3 of what you describe. Here:

https://github.com/gwastro/pycbc/blob/master/pycbc/waveform/waveform.py#L188

the resulting LAL dictionary is send to ChooseTDWaveform (similar for ChooseFDWaveform). In this way lalsimulation has been sent a correct mode array, and what it does with it depends on whoever has coded the waveform.

Regarding other flags. Yes, PyCBC can do this. The full list of flags/options or anything else sent via the "LALDict" is given here:

https://github.com/gwastro/pycbc/blob/master/pycbc/waveform/waveform.py#L89

If what you need is not defined here, you just need to choose an appropriate "parameter name" and add it here, and also in our parameter.py module:

https://github.com/gwastro/pycbc/blob/master/pycbc/waveform/parameters.py

Again, what lalsimulation does with any such parameters is up to lalsimulation and the waveform model in question.

.... Also note that we do define things like chi_p in the parameters library, as such parameters are often useful for parameter estimation. We do not automatically translate this to the component spins that lalsimulation needs, so calling the PyCBC waveform routines with chi_p and not spin[1,2][x,y,z] will result in a non-spinning binary (and hopefully a warning that chi_p is being ignored.

The pycbc_generate_hwinj frontend, may not know about these parameters/options. pycbc_condition_strain is another option for adding injections to frame data (it can be run to read in data (similar to generate_hwinj) without conditioning it at all, and giving it a XML/HDF injection file with the injection to add). Doing this with HDF injection files offers more flexibility if you are wanting to set "non-standard" options.

[albertoaguilarUIB]

Hi Ian,

I've finnaly edited waveform.py and added those sub-version parameters (PhenomXPrecVersion and PhenomXPFinalSpinMod). I will add it also to the parameters.py, even that it's already working fine (if I don't pass them it continues with the default approximant).

Added in waveform.py:
#approximant sub-model
if p['phenomXPrecVersion'] is not None:
lalsimulation.SimInspiralWaveformParamsInsertPhenomXPrecVersion(lal_pars, p['phenomXPrecVersion'])
if p['phenomXPFinalSpinMod'] is not None:
lalsimulation.SimInspiralWaveformParamsInsertPhenomXPFinalSpinMod(lal_pars, p['phenomXPFinalSpinMod'])

I have another question about mode_array, it can translate things like [[2,2], [2,1], [3,3]] or 22, 21, 33, but how can you pass lalsimulation the negative modes? If I leave it empty, eg. for THM, I get the default modes [[2,2], [2,1], [3,3], [4,4], [5,5], [2,-2], [2,-1], [3,-3], [4,-4], [5,-5]]. I'm hand passing it directly to the get_td_waveform() routine.

Thanks for the advise about chi_p...

And I can also try pycbc_condition_strain, seems interesting for me, but I've been focused in pycbc_generate_hwinj (and in my proper work outside this world).

TIA,
Jose Alberto

[spxiwh]

@albertoaguilarUIB I think the negative-m modes are given in the same as positive-m modes, so something like:

from pycbc.waveform import get_td_waveform


for apx in ['IMRPhenomTHM']:
    hp, hc = get_td_waveform(approximant=apx,
                                 mass1=10,
                                 mass2=10,
                                 spin1z=0.9,
                                 delta_t=1.0/16384,
                                 f_lower=40,
                                 mode_array=[[7,-7]])

should be valid (if not scientifically that interesting). I note that this call fails with:

Mode (7,7) is not available by the model.

Implying that the "-7" bit has been ignored. I think that's happening in the model itself. If you see the source code for this model here:

https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/_l_a_l_sim_i_m_r_phenom_t_h_m_8c_source.html#l00728

You can see that this model does not distinguish between positive and negative m modes. So if you have the [2, 2] mode, you must have the [2, -2] mode as well. The code seems to just ignore the sign on the m, when checking the mode array, so giving [2,2] or [2,-2] (or both) would achieve the same thing.