[Question] About Using My Own Waveforms in Bilby

[CalNaTNS]

Hi, I am using some waveforms not included in LALsuite (as well as LALsuite-extra) to analyze several GW events. However, I encountered some troubles when trying to use them in bilby as I am following the example create_own_time_source_model since it doesn’t make use of the parameters conversion function and other things related to real waveforms like SEOBNRv4PHM.
Here's a part of my codes:

# Add the geocent time prior
priors["geocent_time"] = bilby.core.prior.Uniform(
    trigger_time - 0.1, trigger_time + 0.1, name="geocent_time"
)

def SEOBNRE_waveform(time,mass_1,mass_2,spin_1x,spin_1y,spin_1z,spin_2x,spin_2y,spin_2z,eccentricity,luminosity_distance,iota,phase,):
    params = (mass_1,mass_2,spin_1x,spin_1y,spin_1z,spin_2x,spin_2y,spin_2z,eccentricity,luminosity_distance,0,iota,phase,0)
    fMin = minimum_frequency
    srate=sampling_frequency
    Mf_ref=0.002
    waveform, dynamics = calculate_waveform_ep(params, fMin, Mf_ref = Mf_ref, srate=srate, is_coframe=False)
    hpc = waveform.hpc
    plus = hpc.real
    cross = -hpc.imag
    tpeak = waveform.hpc.time[waveform.h22.argpeak]
    plus = TimeSeries(plus,epoch=-tpeak)
    cross = TimeSeries(cross,epoch=-tpeak)
    return {'plus':plus,'cross':cross}


# In this step we define a `waveform_generator`. This is the object which
# creates the frequency-domain strain. In this instance, we are using the
# `lal_binary_black_hole model` source model. We also pass other parameters:
# the waveform approximant and reference frequency and a parameter conversion
# which allows us to sample in chirp mass and ratio rather than component mass
waveform_generator = bilby.gw.WaveformGenerator(
    sampling_frequency=sampling_frequency,
    duration = duration,
    time_domain_source_model=SEOBNRE_waveform,
    parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters,
    waveform_arguments={
        "waveform_approximant": "SEOBNRE",
        "reference_frequency": 11,
    },
)

And here are some errors:

(base) [tanghonglue@swarm01 SEOBNRE]$ vim slurm-8126090.out 

  File "/project/tanghonglue/Injection/SEOBNRE/GW150914.py", line 124, in <module>
    result = bilby.run_sampler(
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/core/sampler/__init__.py", line 264, in run_sampler
    sampler = sampler_class(
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/core/sampler/dynesty.py", line 221, in __init__
    super(Dynesty, self).__init__(
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/core/sampler/base_sampler.py", line 262, in __init__
    self._verify_parameters()
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/core/sampler/base_sampler.py", line 431, in _verify_parameters
    self.log_likelihood(theta)
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/core/sampler/base_sampler.py", line 889, in log_likelihood
    return Sampler.log_likelihood(self, theta)
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/core/sampler/base_sampler.py", line 551, in log_likelihood
    return self.likelihood.log_likelihood()
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/gw/likelihood/base.py", line 853, in log_likelihood
    return self.log_likelihood_ratio() + self.noise_log_likelihood()
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/gw/likelihood/base.py", line 408, in log_likelihood_ratio
    self.waveform_generator.frequency_domain_strain(self.parameters)
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/gw/waveform_generator.py", line 131, in frequency_domain_strain
    return self._calculate_strain(model=self.frequency_domain_source_model,
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/gw/waveform_generator.py", line 171, in _calculate_strain
    self.parameters = parameters
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/gw/waveform_generator.py", line 237, in parameters
00:20 bilby INFO    : Analysis priors:
    new_parameters.pop(key)
KeyError: 'spin_1y'
00:20 bilby INFO    : a_1=Uniform(minimum=0, maximum=0.99, name='a_1', latex_label='$a_1$', unit=None, boundary=None)

and

Traceback (most recent call last):
  File "/project/tanghonglue/Injection/SEOBNRE/GW150914.py", line 124, in <module>
    result = bilby.run_sampler(
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/core/sampler/__init__.py", line 264, in run_sampler
    sampler = sampler_class(
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/core/sampler/dynesty.py", line 221, in __init__
    super(Dynesty, self).__init__(
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/core/sampler/base_sampler.py", line 262, in __init__
    self._verify_parameters()
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/core/sampler/base_sampler.py", line 431, in _verify_parameters
    self.log_likelihood(theta)
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/core/sampler/base_sampler.py", line 889, in log_likelihood
    return Sampler.log_likelihood(self, theta)
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/core/sampler/base_sampler.py", line 551, in log_likelihood
    return self.likelihood.log_likelihood()
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/gw/likelihood/base.py", line 853, in log_likelihood
    return self.log_likelihood_ratio() + self.noise_log_likelihood()
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/gw/likelihood/base.py", line 408, in log_likelihood_ratio
    self.waveform_generator.frequency_domain_strain(self.parameters)
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/gw/waveform_generator.py", line 131, in frequency_domain_strain
    return self._calculate_strain(model=self.frequency_domain_source_model,
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/gw/waveform_generator.py", line 171, in _calculate_strain
    self.parameters = parameters
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/gw/waveform_generator.py", line 237, in parameters
    new_parameters.pop(key)
KeyError: 'iota'
00:20 bilby INFO    : Analysis priors:
00:20 bilby INFO    : a_1=Uniform(minimum=0, maximum=0.99, name='a_1', latex_label='$a_1$', unit=None, boundary=None)
00:20 bilby INFO    : a_2=Uniform(minimum=0, maximum=0.99, name='a_2', latex_label='$a_2$', unit=None, boundary=None)
00:20 bilby INFO    : chirp_mass=bilby.gw.prior.UniformInComponentsChirpMass(minimum=50, maximum=150, name='chirp_mass', latex_label='$\\mathcal{M}$', unit='$M_{\\odot}$', boundary=None)
00:20 bilby INFO    : eccentricity=LogUniform(minimum=0.0001, maximum=0.4, name='eccentricity', latex_label='$e$', unit=None, boundary=None)

Thanks!

[CalNaTNS]

the function calculate_waveform_ep make use of several parameters I will listed down here and it returns evaluated time-domain waveforms includes plus and cross strain.
Parameters used:
mass_1, mass_2, which are component mass of BH in solar mass,
Spin_1x，Spin_1y,..., Spin_2z which are dimentionless spin vector chi of BH
e0, the initial eccentricity at given reference frequency Mf_ref
dL, the luminosity distance in Mpc
zeta_rad, the relativistic anomaly zeta in r = p / (1 + e cos(zeta) (Actually it is not used and I can't find an equivalent parameter in LALInspiral)
iota_rad and beta_rad which are the inclination angle and the phiRef.

[CalNaTNS]

I also tried to convert these parameters in the function SEOBNRE_waveform in order to use the same prior file as the default BBH precession-spin prior. However, it doesn't work and I got some errors that I can't undestand.
Here's my codes:

from bilby.gw.conversion import chirp_mass_and_mass_ratio_to_component_masses, chirp_mass_and_mass_ratio_to_total_mass,bilby_to_lalsimulation_spins

def mass_conversion(parameters):

   # This function is to convert bwtween parameters 

    converted_parameters = parameters.copy()

    converted_parameters['mass_1'],converted_parameters['mass_2'] = chirp_mass_and_mass_ratio_to_component_masses(chirp_mass=converted_parameters['chirp_mass'],mass_ratio=converted_parameters['mass_ratio'])

    converted_parameters['total_mass']=chirp_mass_and_mass_ratio_to_total_mass(chirp_mass=converted_parameters['chirp_mass'],mass_ratio=converted_parameters['mass_ratio'])

    return converted_parameters



def SEOBNRE_waveform(time,chirp_mass,mass_ratio,a_1,a_2,tilt_1,tilt_2,phi_12,phi_jl,psi,theta_jn,luminosity_distance,eccentricity,phase):
    mass_1,mass_2 = chirp_mass_and_mass_ratio_to_component_masses(chirp_mass=chirp_mass,mass_ratio=mass_ratio)
    iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z = bilby_to_lalsimulation_spins(theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1, mass_2,reference_frequency,phase)
    params = (mass_1,mass_2,spin_1x,spin_1y,spin_1z,spin_2x,spin_2y,spin_2z,eccentricity,luminosity_distance,0,iota,phase,0)
    fMin = minimum_frequency
    srate=sampling_frequency
    Mf_ref=0.002
    waveform, dynamics = calculate_waveform_ep(params, fMin, Mf_ref = Mf_ref, srate=srate, is_coframe=False)
    hpc = waveform.hpc
    plus = hpc.real
    cross = -hpc.imag
    tpeak = waveform.hpc.time[waveform.h22.argpeak]
    plus = TimeSeries(plus,epoch=-tpeak)
    cross = TimeSeries(cross,epoch=-tpeak)
    return {'plus':plus,'cross':cross}

waveform_generator = bilby.gw.WaveformGenerator(
    sampling_frequency=sampling_frequency,
    duration = duration,
    time_domain_source_model=SEOBNRE_waveform,

)
priors = bilby.gw.prior.BBHPriorDict(filename="GW190521.prior",conversion_function=mass_conversion)

# Add the geocent time prior
priors["geocent_time"] = bilby.core.prior.Uniform(
    trigger_time - 0.1, trigger_time + 0.1, name="geocent_time"
)

and the output file:

14:37 bilby INFO    : Analysis likelihood noise evidence: -7029.469110814864
Traceback (most recent call last):
  File "/project/tanghonglue/Injection/SEOBNRE/GW190521.py", line 139, in <module>
    result = bilby.run_sampler(
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/core/sampler/__init__.py", line 264, in run_sampler
    sampler = sampler_class(
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/core/sampler/dynesty.py", line 221, in __init__
    super(Dynesty, self).__init__(
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/core/sampler/base_sampler.py", line 262, in __init__
    self._verify_parameters()
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/core/sampler/base_sampler.py", line 431, in _verify_parameters
    self.log_likelihood(theta)
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/core/sampler/base_sampler.py", line 889, in log_likelihood
    return Sampler.log_likelihood(self, theta)
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/core/sampler/base_sampler.py", line 551, in log_likelihood
    return self.likelihood.log_likelihood()
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/gw/likelihood/base.py", line 853, in log_likelihood
    return self.log_likelihood_ratio() + self.noise_log_likelihood()
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/gw/likelihood/base.py", line 420, in log_likelihood_ratio
    per_detector_snr = self.calculate_snrs(
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/gw/likelihood/base.py", line 275, in calculate_snrs
    signal = self._compute_full_waveform(
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/gw/likelihood/base.py", line 709, in _compute_full_waveform
    return interferometer.get_detector_response(signal_polarizations, self.parameters)
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/gw/detector/interferometer.py", line 321, in get_detector_response
    signal_ifo = sum(signal.values()) * mask
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/astropy/units/quantity.py", line 1223, in __mul__
    return super().__mul__(other)
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/gwpy/timeseries/core.py", line 786, in __array_ufunc__
    out = super().__array_ufunc__(ufunc, method, *inputs, **kwargs)
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/gwpy/types/array.py", line 410, in __array_ufunc__
    out = super().__array_ufunc__(function, method, *inputs, **kwargs)
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/astropy/units/quantity.py", line 696, in __array_ufunc__
    raise e
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/astropy/units/quantity.py", line 671, in __array_ufunc__
    result = super().__array_ufunc__(function, method, *arrays, **kwargs)
ValueError: operands could not be broadcast together with shapes (749,) (16385,)
Traceback (most recent call last):
  File "/project/tanghonglue/Injection/SEOBNRE/GW190521.py", line 139, in <module>
    result = bilby.run_sampler(
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/core/sampler/__init__.py", line 264, in run_sampler
    sampler = sampler_class(
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/core/sampler/dynesty.py", line 221, in __init__
    super(Dynesty, self).__init__(
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/core/sampler/base_sampler.py", line 262, in __init__
    self._verify_parameters()
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/core/sampler/base_sampler.py", line 431, in _verify_parameters
    self.log_likelihood(theta)
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/core/sampler/base_sampler.py", line 889, in log_likelihood
    return Sampler.log_likelihood(self, theta)
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/core/sampler/base_sampler.py", line 551, in log_likelihood
    return self.likelihood.log_likelihood()
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/gw/likelihood/base.py", line 853, in log_likelihood
    return self.log_likelihood_ratio() + self.noise_log_likelihood()
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/gw/likelihood/base.py", line 420, in log_likelihood_ratio
    per_detector_snr = self.calculate_snrs(
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/gw/likelihood/base.py", line 275, in calculate_snrs
    signal = self._compute_full_waveform(
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/gw/likelihood/base.py", line 709, in _compute_full_waveform
    return interferometer.get_detector_response(signal_polarizations, self.parameters)
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/bilby/gw/detector/interferometer.py", line 321, in get_detector_response
    signal_ifo = sum(signal.values()) * mask
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/astropy/units/quantity.py", line 1223, in __mul__
    return super().__mul__(other)
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/gwpy/timeseries/core.py", line 786, in __array_ufunc__
    out = super().__array_ufunc__(ufunc, method, *inputs, **kwargs)
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/gwpy/types/array.py", line 410, in __array_ufunc__
    out = super().__array_ufunc__(function, method, *inputs, **kwargs)
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/astropy/units/quantity.py", line 696, in __array_ufunc__
    raise e
  File "/home/tanghonglue/project/miniconda3/envs/bilby_mp/lib/python3.9/site-packages/astropy/units/quantity.py", line 671, in __array_ufunc__
    result = super().__array_ufunc__(function, method, *arrays, **kwargs)
ValueError: operands could not be broadcast together with shapes (732,) (16385,)
srun: error: n0373: task 6: Exited with exit code 1
srun: error: n0373: task 13: Exited with exit code 1

[ColmTalbot]

There are quite a few differences between your source model and the example in the example you mentioned.

• the shape of the returned waveform modes must match the shape of the time argument, note that the output array is defined to have the same length as time and then only the relevant slice of the array is updated with the waveform (here)
• the returned objects should be numpy arrays not TimeSeries. I suspect there will be another error down the line due to using the time series object.

Two things that aren't addressed in that example:

• the typical definition for signal generation is that the reference time should happen at the beginning of the array, this is typically done using numpy.roll so that the inspiral is at the end of the template and the ringdown happens at the beginning.
• often waveform need to be tapered at the beginning and end to avoid Gibbs phenomena

[CalNaTNS]

Thanks for your time and help! However, I have some further questions.

1. I noted that in PyCBC, the position where the amplitude of the WHOLE waveform reaches the peak is set to zero using numpy.roll, while several waveforms choose that of the 2,2 mode of the waveform. Although the difference may be slight, I still want to know whether such a difference makes a difference.
2. After generating the waveform, I suppose I need to expand it to the correct length by adding zeros to it. Can I expand the waveform either at the beginning or at the end, since I will roll the array and these zeros are just located in the middle of the array? I think they are the same.
3. Is it necessary to taper at the end of a GW waveform since the amplitude at ringdown is typically small? I noted that in PyCBC only the beginning of a waveform is tapered after a waveform longer than needed is created.(here) If necessary, is it correct that I taper the waveform after the post-trigger-duration (e.g. 2s after the merger)?

By the way, I always met a Bilby UserWarning:

The sampling was stopped short due to maxiter maxcall limit the delta(log(z))criterion is not achieved: posterior may be poorly sampled
        warnings.warn('The sampling was stopped short due to

It occurs much more often at the end of sampling e.g. when dlogz<10, therefore I suppose it is because the sampling converges much more difficultly here. Does such a problem REALLY affect the posterior greatly, or I can ignore it? How can I fix such a problem?

I would greatly appreciate any guidance and insight you offer, and thank you again for your previous help!