Delete

#!/usr/bin/env python

"""

Tutorial to demonstrate running parameter estimation on a reduced parameter

space for an injected signal.

This example estimates the masses using a uniform prior in both component masses

and distance using a uniform in comoving volume prior on luminosity distance

between luminosity distances of 100Mpc and 5Gpc, the cosmology is Planck15.

"""

import numpy as np

import bilby

# Set the duration and sampling frequency of the data segment that we're

# going to inject the signal into

duration = 1.

sampling_frequency = 256.

ref_geocent_time = 1126259642.5

# Specify the output directory and the name of the simulation.

outdir = 'outdir_BBH'

label = 'bbh_in_bilby_and_vitamin_0703'

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

# Set up a random seed for result reproducibility.  This is optional!

np.random.seed(88170235)

# We are going to inject a binary black hole waveform.  We first establish a

# dictionary of parameters that includes all of the different waveform

# parameters, including masses of the two black holes (mass_1, mass_2),

# spins of both black holes (a, tilt, phi), etc.

injection_parameters = dict(

        mass_1=81.9, mass_2=70.91, a_1=0., a_2=0., tilt_1=0., tilt_2=0.,

        phi_12=0., phi_jl=0., luminosity_distance=1600, theta_jn=1.51, psi=1.54,

        phase=0., geocent_time=1126259642.5+0.22, ra=3.69, dec=-0.94) 

# Fixed arguments passed into the source model

waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',

                          reference_frequency=20., minimum_frequency=20.)

# Create the waveform_generator using a LAL BinaryBlackHole source function

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)

# Set up interferometers.  In this case we'll use two interferometers

# (LIGO-Hanford (H1), LIGO-Livingston (L1). These default to their design

# sensitivity

ifos = bilby.gw.detector.InterferometerList(['H1', 'L1', 'V1'])

ifos.set_strain_data_from_power_spectral_densities(

    sampling_frequency=sampling_frequency, duration=duration,

    start_time=ref_geocent_time - 0.5)

ifos.inject_signal(waveform_generator=waveform_generator,

                   parameters=injection_parameters)

# Set up a PriorDict, which inherits from dict.

# By default we will sample all terms in the signal models.  However, this will

# take a long time for the calculation, so for this example we will set almost

# all of the priors to be equall to their injected values.  This implies the

# prior is a delta function at the true, injected value.  In reality, the

# sampler implementation is smart enough to not sample any parameter that has

# a delta-function prior.

# The above list does *not* include mass_1, mass_2, theta_jn and luminosity

# distance, which means those are the parameters that will be included in the

# sampler.  If we do nothing, then the default priors get used.

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$')

# Initialise the likelihood by passing in the interferometer data (ifos) and

# the waveform generator

likelihood = bilby.gw.GravitationalWaveTransient(

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

# Run sampler.  In this case we're going to use the `dynesty` sampler

#result = bilby.run_sampler(

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

#    nlive=2000, sample='rwalk',walks =100,nact=50,check_point_delta_t=1800,check_point_plot=True)  

result = bilby.run_sampler(

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

    npoints=1000)

# Make a corner plot.

result.plot_corner()

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_d4000'

label = 'BH_encounter_20210722'

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

np.random.seed(88170235)

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)

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

    fixed_waveforms = []

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

    for i in range(num_samples):

        fixed_waveform = {}

        fixed_waveform['signal'] = {}

        for ifo in ifos_list:

            len_wfm = len(waveforms[i]['signal'][ifo])

            peak_index = waveforms[i]['signal'][ifo].index(max(waveforms[i]['signal'][ifo]))

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

            if supposed_peak_index > peak_index:

                for j in range(supposed_peak_index-peak_index):

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

                if len(fixed_waveform['signal'][ifo]) >= 256:

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

                else: 

                    for k in range(256-len(fixed_waveform['signal'][ifo])):

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

            else:

                fixed_waveform['signal'][ifo] = fixed_waveform['signal'][ifo][peak_index-supposed_peak_index:len_wfm]

                if len(fixed_waveform['signal'][ifo]) >= 256:

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

                else: 

                    for j in range(256-len(fixed_waveform['signal'][ifo])):

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

        fixed_waveforms.append(fixed_waveform)

    return fixed_waveforms            

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

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.plot_corner()

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_d6000'

label = 'BH_encounter_20210722'

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

np.random.seed(88170235)

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([6000.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)

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

    fixed_waveforms = []

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

    for i in range(num_samples):

        fixed_waveform = {}

        fixed_waveform['signal'] = {}

        for ifo in ifos_list:

            len_wfm = len(waveforms[i]['signal'][ifo])

            peak_index = waveforms[i]['signal'][ifo].index(max(waveforms[i]['signal'][ifo]))

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

            if supposed_peak_index > peak_index:

                for j in range(supposed_peak_index-peak_index):

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

                if len(fixed_waveform['signal'][ifo]) >= 256:

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

                else: 

                    for k in range(256-len(fixed_waveform['signal'][ifo])):

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

            else:

                fixed_waveform['signal'][ifo] = fixed_waveform['signal'][ifo][peak_index-supposed_peak_index:len_wfm]

                if len(fixed_waveform['signal'][ifo]) >= 256:

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

                else: 

                    for j in range(256-len(fixed_waveform['signal'][ifo])):

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

        fixed_waveforms.append(fixed_waveform)

    return fixed_waveforms            

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

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.plot_corner()