Upload New File

17:23 INFO    : Waveform generator initiated with

  frequency_domain_source_model: bilby.gw.source.lal_binary_black_hole

  time_domain_source_model: None

  parameter_conversion: bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters

17:23 INFO    : No prior given, using default BBH priors in /home/weichangfeng/.local/lib/python3.8/site-packages/bilby/gw/prior_files/precessing_spins_bbh.prior.

17:23 INFO    : Running for label 'BH_encounter_0922_rdm10000_d5000', output will be saved to 'outdir_BH_encounter'

17:23 INFO    : Using lal version 7.1.2

17:23 INFO    : Using lal git version Branch: None;Tag: lalsuite-v6.82;Id: cf792129c2473f42ce6c6ee21d8234254cefd337;;Builder: Unknown User <>;Repository status: UNCLEAN: Modified working tree

17:23 INFO    : Using lalsimulation version 2.5.1

17:23 INFO    : Using lalsimulation git version Branch: None;Tag: lalsuite-v6.82;Id: cf792129c2473f42ce6c6ee21d8234254cefd337;;Builder: Unknown User <>;Repository status: UNCLEAN: Modified working tree

17:23 INFO    : Search parameters:

17:23 INFO    :   mass_1 = Uniform(minimum=30, maximum=160, name='mass_1', latex_label='$m_1$', unit='$M_{\\odot}$', boundary=None)

17:23 INFO    :   mass_2 = Uniform(minimum=30, maximum=160, name='mass_2', latex_label='$m_2$', unit='$M_{\\odot}$', boundary=None)

17:23 INFO    :   luminosity_distance = Uniform(minimum=1000, maximum=3000, name='luminosity_distance', latex_label='$d_L$', unit='$Mpc_{\\odot}$', boundary=None)

17:23 INFO    :   dec = Cosine(minimum=-1.5707963267948966, maximum=1.5707963267948966, name='dec', latex_label='$\\mathrm{DEC}$', unit=None, boundary=None)

17:23 INFO    :   ra = Uniform(minimum=0, maximum=6.283185307179586, name='ra', latex_label='$\\mathrm{RA}$', unit=None, boundary='periodic')

17:23 INFO    :   theta_jn = Sine(minimum=0, maximum=3.141592653589793, name='theta_jn', latex_label='$\\theta_{JN}$', unit=None, boundary=None)

17:23 INFO    :   psi = Uniform(minimum=0.0, maximum=3.141592653589793, name='psi', latex_label='$\\psi$', unit=None, boundary='periodic')

17:23 INFO    :   phase = Uniform(minimum=0, maximum=6.283185307179586, name='phase', latex_label='$\\phi$', unit=None, boundary='periodic')

17:23 INFO    :   geocent_time = Uniform(minimum=1126259644.65, maximum=1126259644.85, name='geocent_time', latex_label='$t_c$', unit='$s$', boundary=None)

17:23 INFO    :   a_1 = 0.0

17:23 INFO    :   a_2 = 0.0

17:23 INFO    :   tilt_1 = 0.0

17:23 INFO    :   tilt_2 = 0.0

17:23 INFO    :   phi_12 = 0.0

17:23 INFO    :   phi_jl = 0.0

17:23 INFO    : Single likelihood evaluation took 2.218e-03 s

17:50 INFO    : Using sampler Dynesty with kwargs {'bound': 'multi', 'sample': 'rwalk', 'verbose': True, 'periodic': None, 'reflective': None, 'check_point_delta_t': 1800, 'nlive': 1000, 'first_update': None, 'walks': 100, 'npdim': None, 'rstate': None, 'queue_size': 1, 'pool': None, 'use_pool': None, 'live_points': None, 'logl_args': None, 'logl_kwargs': None, 'ptform_args': None, 'ptform_kwargs': None, 'enlarge': 1.5, 'bootstrap': None, 'vol_dec': 0.5, 'vol_check': 8.0, 'facc': 0.2, 'slices': 5, 'update_interval': 600, 'print_func': <bound method Dynesty._print_func of <bilby.core.sampler.dynesty.Dynesty object at 0x7f8f2eb57af0>>, 'dlogz': 0.1, 'maxiter': None, 'maxcall': None, 'logl_max': inf, 'add_live': True, 'print_progress': True, 'save_bounds': False, 'n_effective': None, 'maxmcmc': 5000, 'nact': 5}

17:50 INFO    : Checkpoint every check_point_delta_t = 1800s

17:50 INFO    : Using dynesty version 1.0.1

17:50 INFO    : Using the bilby-implemented rwalk sample method with ACT estimated walks

17:50 INFO    : Resume file outdir_BH_encounter/BH_encounter_0922_rdm10000_d5000_resume.pickle does not exist.

17:50 INFO    : Generating initial points from the prior

18:43 INFO    : Written checkpoint file outdir_BH_encounter/BH_encounter_0922_rdm10000_d5000_resume.pickle

18:43 WARNING : singular matrix

18:43 WARNING : Failed to create dynesty state plot at checkpoint

18:43 WARNING : singular matrix

18:43 WARNING : Failed to create dynesty unit state plot at checkpoint

18:43 WARNING : singular matrix

18:43 WARNING : Failed to create dynesty run plot at checkpoint

21:24 INFO    : Written checkpoint file outdir_BH_encounter/BH_encounter_0922_rdm10000_d5000_resume.pickle

22:49 INFO    : Written checkpoint file outdir_BH_encounter/BH_encounter_0922_rdm10000_d5000_resume.pickle

23:20 INFO    : Written checkpoint file outdir_BH_encounter/BH_encounter_0922_rdm10000_d5000_resume.pickle

23:20 WARNING : singular matrix

23:20 WARNING : Failed to create dynesty state plot at checkpoint

23:20 WARNING : singular matrix

23:20 WARNING : Failed to create dynesty unit state plot at checkpoint

23:50 INFO    : Written checkpoint file outdir_BH_encounter/BH_encounter_0922_rdm10000_d5000_resume.pickle

00:21 INFO    : Written checkpoint file outdir_BH_encounter/BH_encounter_0922_rdm10000_d5000_resume.pickle

00:51 INFO    : Written checkpoint file outdir_BH_encounter/BH_encounter_0922_rdm10000_d5000_resume.pickle

01:21 INFO    : Written checkpoint file outdir_BH_encounter/BH_encounter_0922_rdm10000_d5000_resume.pickle

01:51 INFO    : Written checkpoint file outdir_BH_encounter/BH_encounter_0922_rdm10000_d5000_resume.pickle

02:22 INFO    : Written checkpoint file outdir_BH_encounter/BH_encounter_0922_rdm10000_d5000_resume.pickle

02:22 INFO    : Writing 315 current samples to outdir_BH_encounter/BH_encounter_0922_rdm10000_d5000_samples.dat

02:43 INFO    : Written checkpoint file outdir_BH_encounter/BH_encounter_0922_rdm10000_d5000_resume.pickle

02:43 INFO    : Writing 5288 current samples to outdir_BH_encounter/BH_encounter_0922_rdm10000_d5000_samples.dat

02:44 INFO    : Sampling time: 8:52:56.480534

02:44 INFO    : Summary of results:

nsamples: 30302

ln_noise_evidence: -188.255

ln_evidence: -62.977 +/-  0.220

ln_bayes_factor: 125.278 +/-  0.220

import minke.sources

import lalsimulation

import lal

import numpy as np

import minke.distribution

import matplotlib.pyplot as plt

import bilby

outdir = 'outdir_BH_encounter'

label = 'BH_encounter_0907'

bilby.core.utils.setup_logger(outdir=outdir, label=label)

np.random.seed(10000)

nr_path = '../h_psi4/{}'

waveforms_merge = {

    1: "h_m1_L0.9_l2m2_r300.dat",

    2: "h_m2_L0.87_l2m2_r300.dat",

    4: "h_m4_L0.5_l2m2_r300.dat",

    8: "h_m8_L0.35_l2m2_r280.dat",

    16: "h_m16_L0.18_l2m2_r300.dat"

}

nr_waveform = waveforms_merge[1]

duration = 1

sampling_frequency = 256

Nt=duration*sampling_frequency

t0=0.22

ref_geocent_time = 1126259642.5

start_time = ref_geocent_time-duration/2.0

epoch = str(start_time)

num_samples = 1

total_mass = 150

locations=(np.array([3.69]),np.array([-0.94]),np.array([1.54]))

distances=np.array([4000.0])

def generate_for_detector(source, ifos, sampling_frequency, epoch, distance, total_mass, ra, dec, psi):

    """

    Generate an injection for a given waveform.

    Parameters

    ----------

    source : ``minke.Source`` object

       The source which should generate the waveform.

    ifos : list

       A list of interferometer initialisms, e.g. ['L1', 'H1']

    sampling_frequency : int

       The sampling_frequency in hertz for the generated signal.

    epoch : str

       The epoch (start time) for the injection.

       Note that this should be provided as a string to prevent overflows.

    distance : float

       The distance for the injection, in megaparsecs.

    total_mass : float

       The total mass for the injected signal.

    ra : float

        The right ascension of the source, in radians.

    dec : float

        The declination of the source, in radians.

    psi : float

        The polarisation angle of the source, in radians.

    Notes

    -----

    This is very much a draft implementation, and there are a number of 

    things which should be tidied-up before this is moved into Minke,

    including bundling total mass with the source.

    """

    nr_waveform = source.datafile

    waveform = {}

    waveform['signal'] = {}

    waveform['times'] = {}

    for ifo in ifos:

        det = lalsimulation.DetectorPrefixToLALDetector(ifo)

        hp, hx = source._generate(half=True, epoch=epoch, rate=sampling_frequency)[:2]

        h_tot = lalsimulation.SimDetectorStrainREAL8TimeSeries(hp, hx, ra, dec, psi, det)

        waveform['signal'][ifo] = h_tot.data.data.tolist()

        waveform['times'][ifo] = np.linspace(0, len(h_tot.data.data)*h_tot.deltaT, len(h_tot.data.data)).tolist()

    return waveform

waveforms = []

ifos_list = ['H1','L1','V1']

for distance, location in zip(distances, np.vstack(locations).T):

    hyper_object = minke.sources.Hyperbolic(

    datafile = nr_path.format(nr_waveform),

    total_mass = total_mass,

    distance = distance,

    extraction=300

    )

    hyper_object.datafile = nr_waveform

    waveform = generate_for_detector(hyper_object, ifos_list, sampling_frequency, epoch, distance, total_mass, *location)

    waveforms.append(waveform)

peak_index = waveforms[0]['signal']['H1'].index(max(waveforms[0]['signal']['H1']))

supposed_peak_index = int((0.5+0.25)/duration*256)

delta_peak = supposed_peak_index - peak_index

delta_start_time = delta_peak/256

new_start_time = start_time + delta_start_time

new_epoch = str(new_start_time)

new_waveforms = []

for distance, location in zip(distances, np.vstack(locations).T):

    hyper_object = minke.sources.Hyperbolic(

    datafile = nr_path.format(nr_waveform),

    total_mass = total_mass,

    distance = distance,

    extraction=300

    )

    hyper_object.datafile = nr_waveform

    new_waveform = generate_for_detector(hyper_object, ifos_list, sampling_frequency, new_epoch, distance, total_mass, *location)

    new_waveforms.append(waveform)

new_peak_index = new_waveforms[0]['signal']['H1'].index(max(new_waveforms[0]['signal']['H1']))

def fix_waveform(ifos,num_samples,waveforms,delta_peak):

    fixed_waveforms = []

    for i in range(num_samples):

        fixed_waveform = {}

        fixed_waveform['signal'] = {}

        for ifo in ifos_list:

            fixed_waveform['signal'][ifo] = waveforms[i]['signal'][ifo]   

            for j in range(delta_peak):

                fixed_waveform['signal'][ifo] = [0] + fixed_waveform['signal'][ifo]           

            fixed_waveform['signal'][ifo] = fixed_waveform['signal'][ifo][0:256]   

        fixed_waveforms.append(fixed_waveform)

    return fixed_waveforms    

fixed_waveforms = fix_waveform(ifos_list,num_samples,waveforms,delta_peak)

waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',

                          reference_frequency=20., minimum_frequency=20.)

waveform_generator = bilby.gw.WaveformGenerator(

    duration=duration, sampling_frequency=sampling_frequency,

    frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,

    parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters,

    waveform_arguments=waveform_arguments,

    start_time=ref_geocent_time - 0.5)

ifos = bilby.gw.detector.InterferometerList(ifos_list)

ifos.set_strain_data_from_power_spectral_densities(

            sampling_frequency=sampling_frequency, duration=duration,

            start_time=ref_geocent_time-duration/2.0)

priors = bilby.gw.prior.BBHPriorDict()

priors.pop('chirp_mass')

priors.pop('mass_ratio')

priors['mass_1'] = bilby.core.prior.Uniform(name='mass_1', minimum=30, maximum=160, unit='$M_{\\odot}$')

priors['mass_2'] = bilby.core.prior.Uniform(name='mass_2', minimum=30, maximum=160, unit='$M_{\\odot}$')

priors['luminosity_distance'] = bilby.core.prior.Uniform(name='luminosity_distance', minimum=1000, maximum=3000, unit='$Mpc_{\\odot}$')

priors['psi'] = bilby.core.prior.Uniform(name='psi', minimum=0.0, maximum=np.pi, boundary='periodic')

for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl']:

    priors[key] = 0.0

priors['geocent_time'] = bilby.core.prior.Uniform(

    minimum=ref_geocent_time + 0.15,

    maximum=ref_geocent_time + 0.35,

    name='geocent_time', latex_label='$t_c$', unit='$s$')

def fft_signal(signal):

    signal_fd = np.fft.rfft(signal)/Nt

    return signal_fd

signal_fd_array=np.zeros((3,129))

k = 0

for ifo in ifos_list:

    signal_fd_array[k,:] = fft_signal(fixed_waveforms[0]['signal'][ifo]) 

    k += 1

k = 0

for det in ifos_list:

    ifos[k].set_strain_data_from_frequency_domain_strain(signal_fd_array[k],

                                                    sampling_frequency=sampling_frequency,

                                                    duration=duration,

                                                    start_time=start_time)

    k += 1

likelihood = bilby.gw.GravitationalWaveTransient(

    interferometers=ifos, waveform_generator=waveform_generator,priors=priors)

#result = bilby.run_sampler(

#    likelihood=likelihood, priors=priors, label=label,outdir=outdir,sampler='dynesty', 

#    npoints=1000)

result = bilby.run_sampler(

    likelihood=likelihood, priors=priors, label=label,outdir=outdir,sampler='dynesty', 

    nlive=1000, sample='rwalk',walks =100,nact=5,check_point_delta_t=1800,check_point_plot=True)  

result.plot_corner()