Cannot obtain the same PE result (pycbc_inference) with my customed waveform on GW150914 compared with using IMRPhenomXAS

[chendawe]

Hi!

I used get_fd_waveform() to generate hp and hc with the parameter settings example:

events = dict(
mass1 = 30.33,
mass2 = 31.30,
spin1z = 0.06424,
spin2z = 0.08528,
distance = 443.41,
inclination = 221,
coa_phase = 100000,
delta_t = -0.02247,
long_asc_nodes=0,
mean_per_ano=0,
f_lower = 20,
delta_f = 1/16,
f_final = 512,
f_ref=20
)

The gr case is:

hp_gr1, hc_gr1 = wf.get_fd_waveform(approximant='IMRPhenomXAS',
mass1=events['mass2'],
mass2=events['mass1'],
spin1z=events['spin2z'],
spin2z=events['spin1z'],
delta_f=events["delta_f"],
distance=events['distance'],
# mode_array = [[2,-2]],
inclination=events['inclination'],
f_final=events['f_final'],
f_lower=events["f_lower"],
f_ref=events["f_ref"],
coa_phase=events["coa_phase"],
)

My modified IMRPhenomXAS is:

hp_gr_my, hc_gr_my = wf.get_fd_waveform(**{
'approximant': 'dchi_modified_IMRPhenomXAS',
'mass1': events['mass1'],
'mass2': events['mass2'],
'spin1z': events['spin1z'],
'spin2z': events['spin2z'],
'delta_f': events["delta_f"],
'distance': events['distance'],
# 'mode_array': [[2, -2]],
'inclination': events['inclination'],
'f_final': events['f_final'],
'f_lower': events["f_lower"],
'f_ref': events["f_ref"],
'coa_phase': events["coa_phase"],
})

I did not include

long_asc_nodes

and

mean_per_ano

for now.

The results of using various parameters combinations show that the approximants 'IMRPhenomXAS' and 'dchi_modified_IMRPhenomXAS' (based on WF4py) generate exactly the same GR waveform results up to some negelegible numerical errors (1 - inner product between the results of two cases are ~1e-8 and ~1e-10; the numerical differences of phase, amplitude and the real and imag parts of hp and hc are also negeligible).

Besides, I set the epoch of my 'dchi_modified_IMRPhenomXAS' the same as the default setting:

epoch=_lal.LIGOTimeGPS(-16.0)

However, the PE results of using 'dchi_modified_IMRPhenomXAS' on GW150914 still did not fit well with the result of using 'IMRPhenomXAS':

IMRPhenomXAS

dchi_modified_IMRPhenomXAS:

---

My questions are:

1. When hp and hc are generated, should there be more parameters to be passed into FrequencySeries()?
2. After hp and hs are generated, and projected to the detector, should this be also included in the plugin?
3. Are polarization, long_asc_nodes and mean_per_ano or else related to the problems?

These are the parameters properly dealt with in my generator:

approximant
f_lower
f_ref
mass1, mass2 (srcmchirp, q)
spin1z
spin2z
distance (comoving_volume)
inclination
coa_phase

These are the parameters not dealt with my generator, but apperaing in the pycbc_inference .config:

trigger_time
delta_tc
ra
dec

[spxiwh]

Do the two waveforms agree to 1E-8 precision when compared with PyCBC?

https://pycbc.org/pycbc/latest/html/waveform.html#calculating-the-match-between-waveforms

[ahnitz]

@chendawe For a more detailed conversation, I might suggest joining the slack. Send me an email, if you don't know how.

When hp and hc are generated, should there be more parameters to be passed into FrequencySeries()?

The plugin should return hp/hc as FreuencySeries, where delta_f is the one requested by the user and the epoch is set so that 0 is at the reference time (usually merger by convention). This point should also be aligned in the waveform so that it is at the first sample if you were to FFT it to the time domain. That's pretty much it. The rest of the arguments are for your use and by convention (e..g masses are always mass1/mass2).

After hp and hs are generated, and projected to the detector, should this be also included in the plugin?

Projection happens internally. Unless you are modifying GR in a way that changes the propagation effects. You should just return hp/hc. Nothing for you to do.

Are polarization, long_asc_nodes and mean_per_ano or else related to the problems?

No. Not unless you are doing something that creates non-circular orbits and need to describe the additional degrees of freedom.

@chendawe If you want a simple test case, I'd first have a very basic plugin that does nothing but create a waveform from a different approximant and directly returns it. That should give you the same results either way and might be a good starting point for you to track down the issue.

e.g. somethign like

def my_new_waveform(**params):
       params.pop('approximant')
       return get_fd_waveform(approximant='IMRPhenomXAS', **params)

That should essentially just a rename of IMRPhenomXAS then and should give you the same results back (assuming all settings are really the same). From there you can see where the difference is introduced by graduate change to your plugin.

[chendawe]

@spxiwh

Yes!

[chendawe]

@ahnitz
Very grateful for your detailed comment!
And sorry for that slack is not available in my country.

But I will try to use the renamed IMRPhenomXAS for further debugs. Great thanks!

I will report my debug results later!

[ahnitz]

@chendawe OK.

From past experience, I've encountered a number of ways to get this kind of issue. I'll just list them off in case you haven't check all of them.

(1) the sampler didn't converge correctly (this one is a bit tricky to debug so usually I like to eliminate a problem with the likelihood / waveform model first). If this is reliably inconsistent, it is probably not the issue (e.g. rerun with different seed, rerun with more live points / walkers).

(2) inconsistent parameters for the PE (e.g. accidentally truncating the waveform or having it wrap around the dataset (recall Fourier domain has circular continuity at boundaries) because the data analyzed was shorter than the actual template. Using low / high pass choices for the data and the model which are inconsistent, etc. Even accidentally changing a prior / prior range without noticing. Or perhaps waveform generation just failed completely for some parameters and you left on the option to ignore waveform failures and didn't notice.

(3) waveform model not actually returning something identical / consistent. This could be many things. I've had a case where someone didn't realize what the plus / cross polarizations were so treated them as the same thing (obviously this does not work). It can be subtle things that might not show up in match though such as breaking the alignment conventions (the alignment of the merger at 0/ start of the vector I discussed above). A very common one here is a mislabeled parameter. Our code is very forgiving (in this case to a fault) on parameter names and allows you to use whatever you want. The downside is that if you pass something that is actually a typo, we don't know that it is not what you intended and happily pass it along. This can sometimes lead to people thinking they are using a parameter when they are not.

[chendawe]

@ahnitz

So kind of you for your illuminating suggestions!

(1) I've tried to rerun many times or PE over other events (based on the configs in 4-ogc), it seems not like a problem from convergence. I think as the waveform results are generated the same, then on the same config should two models achieve the same results including the successful convergence.

(2) The length and wrap might be the problem, cause in my impression, PyCBC will complete the waveform by a 2^n > f_final rule.
I did not implement this, but I tried to set f_final = 2048 in the config. The problem did not disappear.

Real parts of hp(s) with f_final = 1234Hz:

(3)
I believe my function did return the hp and hc correctly, for I've checked the phase, amplitude and hp, hc numerical errors.
And I derived hp and hc from the phase and amplitude by the formulas (20) in IMRPhenomXAS which introduces iota and coa_phase. These two parameters are deviated to ensure hp and hc match well with the original IMRPhenomXAS results too.

The name of parameters are all PyCBC conventions are all conserved

However about

the alignment of the merger at 0/ start of the vector I discussed above

, this might be the problem cause I do not intentionally deal with anything on this especially if the alignment is always not a trivial procedure.
But I am sure that my hp and hc are the same as the ones generated by get_fd_waveform("IMRPhenomXAS"). And I tried iFFT once which also gave me the same TD waveforms but with different time start points:

This is why I tried to set epoch=_lal.LIGOTimeGPS(-16.0, 0) mannualy (default _lal.LIGOTimeGPS(0, 0)). This align the iFFT-ed IMRPhenomXAS and my IMRPhenomXAS TD results well, though does not solve the PE problem.

Anyway I will still keep my eye on these possibilities mentioned in the list when tracing back!

Oh and:

[model]
name = marginalized_polarization
low-frequency-cutoff = 20
check-for-valid-times =
shift-psd-times-to-valid =
ignore-failed-waveforms =

I will check the possible ignorings too!

[chendawe]

@spxiwh

I tried my waveform on the grid:

num = 5
mass1_s = np.linspace(3, 100, num)
mass2_s = np.linspace(3, 100, num)
spin1z_s = np.linspace(-0.99, 0.99, num)
spin2z_s = np.linspace(-0.99, 0.99, num)

and the match distribution is like:

Would this be crucial?
If so, does this come from some false codings underneath my realisation of IMRPhenomXAS, or the numerical difference by using Python (WF4Py) and C language (LALSuite)?

[chendawe]

@ahnitz

The result of using renamed IMRPhenomXAS is no problem;

I filtered those waveforms generated by my model with mismatches over 1e-4 by replacing them with PyCBC implementations, the result work well too;
Then I think I found the cause that the choices of different versions of phase and amplitude changes the mismatches distribution greatly:

InsPhaseVersion=114, IntPhaseVersion=105, IntAmpVersion=105

InsPhaseVersion=115, IntPhaseVersion=105, IntAmpVersion=104

InsPhaseVersion=115, IntPhaseVersion=105, IntAmpVersion=105

---

Then I choose InsPhaseVersion=115, IntPhaseVersion=105, IntAmpVersion=105 to PE GW150914 and GW170729, which works pretty well (GW170729's lominosity distance median is 300 Mpc less deviated from official release 2400 from https://peviewer.igwn.org/, seeming statistically acceptable?)
However for GW170104, GW170608, PE results deviated quite a lot:
GW170104:

GW170608:

These might be traced back to my too broad srcmchirp and q settings being [10, 100] and [1, 10] not following 4-ogc, that seems even causing the luminosity distance out of LIGO O2 range. Any way I'm exploring this possibility, or another possibility that these events are in the mass and spin regime where my now version settings are not appropriate, resulting in relatively high mismatches in the regime?

---

Finally, my question is:
Does PyCBC automatically choose different phase and amplitude versions for the best performances in different mass and spin regimes?