search

bilby-dev / bilby

A unified framework for stochastic sampling packages and gravitational-wave inference in Python. Note that we are currently transitioning from git.ligo.org/lscsoft/bilby, please bear with us!
https://bilby-dev.github.io/bilby/
MIT License
55 stars 63 forks source link

[Question] About Using My Own Waveforms in Bilby #821

Open CalNaTNS opened 2 hours ago

CalNaTNS commented 2 hours ago

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 commented 2 hours ago

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 commented 20 minutes ago

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