Added draft of the data release.

LICENSE

+Copyright 2022 Weichagnfeng Guo, Daniel Williams,  Ik Siong Heng, Hunter Gabbard, Yeong-Bok Bae, Gungwon Kang, Zong-Hong Zhu
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
\ No newline at end of file

README.md

+# Mimicking Mergers: Mistaking Black Hole Captures as Mergers
+
+This repository contains the code and data release to accompany "Mimicking Mergers: Mistaking Black Hole Captures as Mergers" (https://arxiv.org/abs/2203.06969).
+We provide instructions to reproduce the analyses below.
+
+If you reuse code or data from this publication we would appreciate it if you cite the paper: citation details are provided at the end of this README file.
+The data release is made available under the MIT license, which is detailed in a file called "LICENSE".
+
+We use three different analyses in this work. 
+One uses the VItamin python package, and one the bilby package.
+The final one is an analysis of results from VItamin.
+We detail the process for running each seperately.
+
+## VItamin analyses
+
+The code for the VItamin analyses can be found in the VItamin directory.
+In order to run these analyses you will require VItamin, installed from a development branch.
+Please note that at present VItamin is not compatible with Python versions later than 3.6.
+
+This can be downloaded using `git`:
+
+```
+git clone https://github.com/hagabbar/vitamin_b
+git checkout weichangfeng_vitamin_c
+```
+Then the dependencies should be installed by running
+```
+pip install -r requirements.txt
+```
+in the root of the repository.
+VItamin may then be installed by running
+```
+pip install .
+```
+
+This will allow the use of a BBH model we used in this work which has been trained with a wider prior on component masses.
+The model itself is put in the 'bbh_mimics' home directory.
+The instruction for the analysis and specifying a model can be found here:
+https://hagabbar.github.io/vitamin_c/quickstart.html
+The VItamin analysis outputs the posterior corner plot for the input data.
+
+
+## BH capture data creation
+
+We generate BH capture data using the 'vitamin/BH_capture_data_creation.ipynb' notebook.
+This notebook also includes code for signal whitening and generating noise realisations, in order to produce timeseries which are suitable for analysis with bilby.
+In this notebook, the distributions of sky location are specified and contain 100 samples, the distance injection and the noise realisation is fixed, and the waveform is specified. 
+Therefore, the script generates 100 data files per waveform. 
+Then we perform VItamin analysis on them, and look for the posterior which is visually similar to that of BBH. If there is not, we adjust distance injection for another turn of data creation and make the following analysis, until obtaining the appropriate posterior.
+The process should also be performed for other waveforms with different mass ratios.
+ 
+## Bilby analyses
+
+In this section, we introduce the code realted to Bilby part of our work, which are put in the 'bilby' directory.
+We employ 'bilby/PE_dynesty.py' to apply Bayesian inference on BH capture data with the dynesty sampler.
+Using this script, we first reproduce one BH capture signal; then inject the signal into simulated detector noise.
+We then specify the BBH model, prior and likelihood, and finally perform parameter estimation on it. 
+We can easily switch the spinning model to the non-spinning one, by setting a delta distribution prior of six spins at zero.
+
+## JS divergence analysis
+
+In this section, we introduce the codes for JS divergence analysis, which are put in the 'jsd' directory.
+Currently VItamin only outputs the posterior corner plot, though it calculates samples in the function. 
+For JS divergence analysis, we make it output samples file as well, by modifying VItamin file 'run_vitamin.py' located in 'vitamin_c/vitamin_c/'.
+We add the following codes in line 1491:
+>    for i in range(num_timeseries):
+>       for j in range(len(params['inf_pars'])):
+>            samples_for_one_timeseries = samples[i].transpose()
+>            samples_for_one_parameter=samples_for_one_timeseries[j]
+>            np.savetxt(f'vitamin_sample_timeseries-{i}_parameter-{j}.txt', samples_for_one_parameter, fmt="%.8f")
+Then VItamin can generate both corner plot and sample data file for the input data.
+
+### BH capture data creation and analysis
+We use the script `BH_capture_data_creation.ipynb`, from `jsd` directory, to inject one BH capture signal into 100 noise realisations, and finally generate 100 data file. We do this for all waveforms and obtain their samples data files by VItamin analysis.
+
+### BBH data creation and analysis
+We use the script `BBH_data_creation.ipynb`, from `jsd` directory, to inject one BBH signal into 100 noise realisations, and finally generate 100 data files. 
+In this analysis, we have two kinds of BBH file to produce. 
+First, it's the reference BBH, of which the injection is set in advance, and can be easily generate with 100 noise realisations. 
+The other one is the BBH recovered by the average peak values of posterior in Sec 3.2. 
+We generate four recovered BBH signals corresponding to the four waveforms.
+Each of them is injected into 100 noise realisations.
+We then obtain their samples data files by VItamin analysis.
+
+### JS divergence calculation 
+
+We use the script `JSD_calculation.py` to load all the samples data files, and calculate their JS divergences.
+The distributions of JS divergence are plotted as histogram.
+We also compute JS divergence at 90% confidence level, and some useful percentage in this script. 
+
+
+# Citation
+
+We would appreciate it if you cite this work if you make use of any of the code or data in this release.
+For your convenience below is a bibtex citation.
+
+```
+@ARTICLE{Guo+2022Mimick,
+       author = {{Guo}, Weichangfeng and {Williams}, Daniel and {Heng}, Ik Siong and {Gabbard}, Hunter and {Bae}, Yeong-Bok and {Kang}, Gungwon and {Zhu}, Zong-Hong},
+        title = "{Mimicking Mergers: Mistaking Black Hole Captures as Mergers}",
+      journal = {arXiv e-prints},
+     keywords = {General Relativity and Quantum Cosmology, Astrophysics - High Energy Astrophysical Phenomena},
+         year = 2022,
+        month = mar,
+          eid = {arXiv:2203.06969},
+        pages = {arXiv:2203.06969},
+archivePrefix = {arXiv},
+       eprint = {2203.06969},
+ primaryClass = {gr-qc},
+       adsurl = {https://ui.adsabs.harvard.edu/abs/2022arXiv220306969G},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+```

bilby/parameter_estimation_dynesty.py

+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()