Upload New File

# Signal whitening.

def whiten_signal(signal_fd,ifo):

    #signal_fd = np.fft.rfft(signal_td)/sampling_frequency

    whitened_signal_fd = signal_fd/ifo.amplitude_spectral_density_array

    # The following signal whitening will contain 100 sets of white noise since we don't fix the random seed here.

    whitened_signal_td = np.sqrt(2.0*Nt)*np.fft.irfft(whitened_signal_fd) + 1.0*np.random.randn(256)

    return whitened_signal_td

        injection_parameters['geocent_time'],injection_parameters['ra'],injection_parameters['dec'],

        injection_parameters['psi']]

# Record the injections.

# Record the injections in the data file.

x_data_sum = np.zeros((num_samples,7))

for i in range(num_samples):

    x_data_sum[i] = x_data_list

import bilby

get_ipython().run_line_magic('matplotlib', 'inline')

# Specify the waveform we use.

mass_ratio= 1

# In this study we generate one signal that combines with 100 noise realisations.

num_samples = 100

# The BH capture waveform sets we use in the data production.

nr_path = './h_psi4/{}'

waveforms_merge = {

    1: "h_m1_L0.9_l2m2_r300.dat",

    8: "h_m8_L0.35_l2m2_r280.dat",

    16: "h_m16_L0.18_l2m2_r300.dat"

}

# The extraction is corresponding to the waveform.

extraction = {

    1: 300,

    2: 300,

    16: 300

}

# In this study we generate 100 data files, which contain one signal combining 100 noise realisations.	.

num_samples = 100

# Specify the waveform by mass ratio.

mass_ratio= 1

nr_waveform = waveforms_merge[mass_ratio]

# Set the injections.

duration = 1

sampling_frequency = 256

Nt=duration*sampling_frequency

}

# Only one set of injections is used in this script. 

# Only one set of sky location injections is used in this script. 

total_mass = 150

ra = np.random.uniform(0.89,0.89,num_samples)

dec = np.random.uniform(-0.94,-0.94,num_samples)

waveforms = []

# Set the detector network we use.

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

# Generate the mock signal.

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

    hyper_object = minke.sources.Hyperbolic(

# The start time was specified when generating mock signal, while the merger time t0 and signal length could not be. We do the following process to make the signal have a start time of 1126259642.5s in GPS time, and a t0 of 0.22s. 

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

supposed_peak_index = int((0.5+0.22)/duration*sampling_frequency)

delta_peak = supposed_peak_index - peak_index

delta_start_time = delta_peak/Nt

new_waveforms = []

# By reproducing mock signal with a new start time, we make the signal have a supposed t0=0.22s, while maintaining a correct start time.

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

    hyper_object = minke.sources.Hyperbolic(

            sampling_frequency=sampling_frequency, duration=duration,

            start_time=ref_geocent_time-duration/2.0)

# Signal whitening.

def whiten_signal(signal,ifo):

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

    whitened_signal_fd = signal_fd/ifo.amplitude_spectral_density_array

    # The following signal whitening will contain 100 sets of white noise since we don't fix the random seed here.

    whitened_signal_td = np.sqrt(2.0*Nt)*np.fft.irfft(whitened_signal_fd) + 1.0*np.random.randn(int(Nt))

    return whitened_signal_td

# Produce 100 sets of detector network(H1,L1,V1), for the following signal whitening.

np.random.seed(10000)

whitened_signal_td_array=np.zeros((num_samples,3,int(Nt)))

    name_list.append(name)

    data_set_list.append(f1[f'{name}'])

# Calculate m1,m2 injection.

m2=total_mass/(mass_ratio+1)

m1=total_mass-m2

# Record the injections.

# Record the injections in the data file.

x_data_sum = np.zeros((num_samples,7))

for i in range(num_samples):

    x_data_sum[i] = [m1,m2,distance_set[mass_ratio],0.22,locations[2][i],locations[0][i],locations[1][i]]

import bilby

# # Load the posterior samples from VItamin.

# # Load the posterior samples from VItamin, which can be divied into three types.

num_samples=100

# Posterior of reference BBH

m1_BBH=np.zeros((num_samples,10000))

m2_BBH=np.zeros((num_samples,10000))

t0_BBH=np.zeros((num_samples,10000))

# Posterior of BH capture,

m1_1=np.zeros((num_samples,10000))

m2_1=np.zeros((num_samples,10000))

t0_1=np.zeros((num_samples,10000))

m2_16=np.zeros((num_samples,10000))

t0_16=np.zeros((num_samples,10000))

# Posterior of the recovered signal of BH capture.

m1_recover_1=np.zeros((num_samples,10000))

m2_recover_1=np.zeros((num_samples,10000))

t0_recover_1=np.zeros((num_samples,10000))

    t0_recover_16[i] = np.loadtxt(f'/mnt/d/vitamin/encounter-like_BBH/m16/results/samples/vitamin_sample_timeseries-{i}_parameter-3.txt')  

# Create JSD arrays for the following calculation.

jsd_m1_BBH=np.zeros((int(num_samples/2)))

jsd_m2_BBH=np.zeros((int(num_samples/2)))

jsd_t0_BBH=np.zeros((int(num_samples/2)))

    return js_array

# Calculate the JSDs.

# JSD between reference BBH with different noise. 

for i in range(int(num_samples/2)):

    jsd_t0_BBH[i] = calculate_js(t0_BBH[2*i],t0_BBH[2*i+1]).mean()

    jsd_m1_BBH[i] = calculate_js(m1_BBH[2*i],m1_BBH[2*i+1]).mean()

    jsd_m2_BBH[i] = calculate_js(m2_BBH[2*i],m2_BBH[2*i+1]).mean()

# JSD between reference BBH and BH capture.

for i in range(num_samples):

    jsd_m1_ref_1[i] = calculate_js(m1_BBH[i],m1_1[i]).mean()

    jsd_m2_ref_1[i] = calculate_js(m2_BBH[i],m2_1[i]).mean()

    jsd_m2_ref_16[i] = calculate_js(m2_BBH[i],m2_16[i]).mean()

    jsd_t0_ref_16[i] = calculate_js(t0_BBH[i],t0_16[i]).mean()

# JSD between BH capture and its recovered BBH.

for i in range(num_samples):

    jsd_m1_recover_1[i] = calculate_js(m1_1[i],m1_recover_1[i]).mean()

    jsd_m2_recover_1[i] = calculate_js(m2_1[i],m2_recover_1[i]).mean()

    jsd_t0_recover_16[i] = calculate_js(t0_16[i],t0_recover_16[i]).mean()

# # Histogram weights

# # Histogram weights that make the heights be the fraction of samples.

weights = np.ones_like(jsd_m1_ref_1)/float(len(jsd_m1_1))

weights_BBH = np.ones_like(jsd_m1_BBH)/float(len(jsd_m1_BBH))

label = 'BH_encounter_0907'

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

np.random.seed(10000)

# Specify the waveform we use.

# In the Bilby analysis we only use one data.

num_samples = 1

# The BH capture waveform sets we use in the data production..

nr_path = '../h_psi4/{}'

waveforms_merge = {

    1: "h_m1_L0.9_l2m2_r300.dat",

    8: "h_m8_L0.35_l2m2_r280.dat",

    16: "h_m16_L0.18_l2m2_r300.dat"

}

nr_waveform = waveforms_merge[1]

# The extraction is corresponding to the waveform.

extraction = {

    1: 300,

    2: 300,

    4: 300,

    8: 280,

    16: 300

}

# Specify the waveform by mass ratio.

mass_ratio= 1

nr_waveform = waveforms_merge[mass_ratio]

# Set the injections.

duration = 1

sampling_frequency = 256

Nt=duration*sampling_frequency

start_time = ref_geocent_time-duration/2.0

epoch = str(start_time)

# Set injection parameters.

num_samples = 1

total_mass = 150

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

distances=np.array([4000.0])

distances=np.array([5000.0])

# Generate BH capture signal.

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

waveforms = []

# Set the detector network we use.

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

# Generate the mock signal.

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

    hyper_object = minke.sources.Hyperbolic(

# The start time was specified when generating mock signal, while the merger time t0 and signal length could not be. We do the following process to make the signal have a start time of 1126259642.5s in GPS time, and a t0 of 0.22s. 

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

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_waveforms = []

# By reproducing mock signal with a new start time, we make the signal have a supposed t0=0.22s, while maintaining a correct start time.

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

    hyper_object = minke.sources.Hyperbolic(

# Specify the model that we use to perform parameter estimation

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,

    start_time=ref_geocent_time - 0.5)

# Set the detector network and the noise.

np.random.seed(10000)

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)

# Set the priors.

# Set the priors. For the non-spinning BBH analysis, we set prior of six spins as zero.

priors = bilby.gw.prior.BBHPriorDict()

priors.pop('chirp_mass')

priors.pop('mass_ratio')

likelihood = bilby.gw.GravitationalWaveTransient(

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

# We use dynesty as the sampler.

# We use dynesty as the sampler and perform parameter estimation.

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)