Full PE scripts¶
1. LISA¶
full_pe_MBHB_LISA.py¶
1import bilby
2import lal
3import numpy as np
4import taichi as ti
5# local_rank_id = int(os.environ['MPI_LOCALRANKID'])
6# gpus_list = os.environ['GPU_DEVICE_ORDINAL'].split(',')
7# selected_gpu = gpus_list[local_rank_id%len(gpus_list)]
8# os.environ['CUDA_VISIBLE_DEVICES'] = selected_gpu
9# ti.init(arch=ti.gpu, default_fp=ti.f64, cpu_max_num_threads=1, offline_cache=False)
10ti.init(arch=ti.cpu, default_fp=ti.f64, cpu_max_num_threads=1, offline_cache=False)
11
12from pespace.detector.antenna import InterferometerAntenna, FDResponseModelMarset2018
13from pespace.detector.tdi import TDIChannelData, FDMichelsonConstantEqualArm
14from pespace.detector.orbit import available_orbit_models
15from pespace.inference.interface_bilby import LikelihoodBilbyInterface
16from tiwave.waveforms import IMRPhenomXAS
17
18########################################################################################
19# set injection
20tdi_gen = "2.0"
21tdi_chan = ("A", "E", "T")
22
23dt = 10
24f_min = 1e-4
25f_max = 0.5*(1/dt)
26f_ref = f_min
27t_start = 0.0
28
29num_tsamples = 2**np.ceil(np.log2(7*lal.DAYJUL_SI/dt))
30duration = num_tsamples * dt
31before_tc = 0.8 * duration
32after_tc = 0.2 * duration
33tc = t_start + before_tc
34
35params = dict(
36 total_mass=3e6,
37 mass_ratio=0.6,
38 chi_1=0.75,
39 chi_2=0.62,
40 luminosity_distance=56000.0,
41 inclination=0.4,
42 reference_phase=1.3,
43 ecliptic_longitude=1.375,
44 ecliptic_latitude=-1.2108,
45 polarization=2.659,
46 coalescence_time=tc,
47)
48
49tdi_data = TDIChannelData()
50tdi_data.set_fd_data_from_zero(
51 channels=tdi_chan,
52 duration=duration,
53 delta_time=dt,
54 start_time=t_start,
55 minimum_frequency=f_min,
56 maximum_frequency=f_max,
57 )
58tdi_data.set_fd_noise_power_density_from_model("LISA_SciRDv1", tdi_generation=tdi_gen)
59noise = tdi_data.get_fd_noise_realization()
60tdi_data.add_into_fd_data(noise)
61
62orbit_model = available_orbit_models['LISA_analytic']
63response_model = FDResponseModelMarset2018()
64tdi_combination = FDMichelsonConstantEqualArm(generation=tdi_gen, orthogonal=True)
65lisa = InterferometerAntenna(
66 name="lisa",
67 tdi_data=tdi_data,
68 orbit_model=orbit_model,
69 response_model=response_model,
70 tdi_combination=tdi_combination,
71)
72
73wf_xas = IMRPhenomXAS(tdi_data.frequency_samples, f_ref, parameter_conversion=bilby.gw.conversion.generate_component_masses)
74wf_xas.update_waveform(params)
75lisa.inject_signal(
76 wf_xas.waveform_container,
77 params["ecliptic_longitude"],
78 params["ecliptic_latitude"],
79 params["polarization"],
80 params["coalescence_time"],
81)
82
83########################################################################################
84# set sampling
85likelihood = LikelihoodBilbyInterface(
86 waveform=wf_xas,
87 detector=lisa,
88 channels=('A', 'E', 'T')
89 )
90
91priors = {}
92priors['chi_1'] = bilby.core.prior.Uniform(name='chi_1', minimum=-0.99, maximum=0.99)
93priors['chi_2'] = bilby.core.prior.Uniform(name='chi_2', minimum=-0.99, maximum=0.99)
94priors['ecliptic_longitude'] = bilby.core.prior.Uniform(name='ecliptic_longitude', minimum=0, maximum=2 * lal.PI, boundary='periodic')
95priors['ecliptic_latitude'] = bilby.core.prior.Cosine(name='ecliptic_latitude')
96priors['inclination'] = bilby.core.prior.Sine(name='inclination')
97priors['polarization'] = bilby.core.prior.Uniform(name='polarization', minimum=0, maximum=lal.PI, boundary='periodic')
98priors['reference_phase'] = bilby.core.prior.Uniform(name='reference_phase', minimum=0, maximum=2 * lal.PI, boundary='periodic')
99priors['coalescence_time'] = bilby.core.prior.Uniform(name='coalescence_time', minimum=tc-100, maximum=tc+100)
100priors['chirp_mass'] = bilby.gw.prior.UniformInComponentsChirpMass(name='chirp_mass', minimum=5e5, maximum=2e6)
101priors['mass_ratio'] = bilby.gw.prior.UniformInComponentsMassRatio(name='mass_ratio', minimum=0.05, maximum=0.99)
102priors['luminosity_distance'] = bilby.core.prior.LogUniform(name='luminosity_distance', minimum=1e4, maximum=5e6)
103priors = bilby.core.prior.PriorDict(dictionary=priors)
104
105label = "inj_MBHB_LISA"
106outdir = f"outdir_{label}"
107result = bilby.run_sampler(
108 likelihood=likelihood,
109 priors=priors,
110 outdir=outdir,
111 label=label,
112 sampler='pymultinest',
113 npoints=2048,
114 )
2. LISA-Taiji-TianQin network and waveform with higher modes¶
full_pe_MBHB_network_HM.py¶
1import bilby
2import lal
3import numpy as np
4import taichi as ti
5# local_rank_id = int(os.environ['MPI_LOCALRANKID'])
6# gpus_list = os.environ['GPU_DEVICE_ORDINAL'].split(',')
7# selected_gpu = gpus_list[local_rank_id%len(gpus_list)]
8# os.environ['CUDA_VISIBLE_DEVICES'] = selected_gpu
9# ti.init(arch=ti.gpu, default_fp=ti.f64, cpu_max_num_threads=1, offline_cache=False)
10ti.init(arch=ti.cpu, default_fp=ti.f64, cpu_max_num_threads=1, offline_cache=False)
11
12from pespace.detector.antenna import InterferometerAntenna, FDResponseModelMarset2018
13from pespace.detector.tdi import TDIChannelData, FDMichelsonConstantEqualArm
14from pespace.detector.orbit import available_orbit_models
15from pespace.inference.interface_bilby import LikelihoodBilbyInterface
16from tiwave.waveforms import IMRPhenomXHM
17
18########################################################################################
19# common settings
20# set injection
21tdi_gen = "2.0"
22tdi_chan = ("A", "E", "T")
23
24dt = 10
25f_min = 1e-4
26f_max = 0.5*(1/dt)
27f_ref = f_min
28t_start = 0.0
29
30num_tsamples = 2**np.ceil(np.log2(7*lal.DAYJUL_SI/dt))
31duration = num_tsamples * dt
32before_tc = 0.8 * duration
33after_tc = 0.2 * duration
34tc = t_start + before_tc
35
36params = dict(
37 total_mass=3e6,
38 mass_ratio=0.6,
39 chi_1=0.75,
40 chi_2=0.62,
41 luminosity_distance=56000.0,
42 inclination=0.4,
43 reference_phase=1.3,
44 ecliptic_longitude=1.375,
45 ecliptic_latitude=-1.2108,
46 polarization=2.659,
47 coalescence_time=tc,
48)
49
50# detector settings
51tdi_data_lisa = TDIChannelData()
52tdi_data_lisa.set_fd_data_from_zero(
53 channels=tdi_chan,
54 duration=duration,
55 delta_time=dt,
56 start_time=t_start,
57 minimum_frequency=f_min,
58 maximum_frequency=f_max,
59 )
60tdi_data_lisa.set_fd_noise_power_density_from_model("LISA_SciRDv1", tdi_generation=tdi_gen)
61noise_lisa = tdi_data_lisa.get_fd_noise_realization()
62tdi_data_lisa.add_into_fd_data(noise_lisa)
63lisa = InterferometerAntenna(
64 name="lisa",
65 tdi_data=tdi_data_lisa,
66 orbit_model=available_orbit_models['LISA_analytic'],
67 response_model=FDResponseModelMarset2018(),
68 tdi_combination=FDMichelsonConstantEqualArm(generation=tdi_gen, orthogonal=True),
69)
70
71tdi_data_taiji = TDIChannelData()
72tdi_data_taiji.set_fd_data_from_zero(
73 channels=tdi_chan,
74 duration=duration,
75 delta_time=dt,
76 start_time=t_start,
77 minimum_frequency=f_min,
78 maximum_frequency=f_max,
79 )
80tdi_data_taiji.set_fd_noise_power_density_from_model("Taiji_TDC", tdi_generation=tdi_gen)
81noise_taiji = tdi_data_taiji.get_fd_noise_realization()
82tdi_data_taiji.add_into_fd_data(noise_taiji)
83taiji = InterferometerAntenna(
84 name="taiji",
85 tdi_data=tdi_data_taiji,
86 orbit_model=available_orbit_models['Taiji_analytic'],
87 response_model=FDResponseModelMarset2018(),
88 tdi_combination=FDMichelsonConstantEqualArm(generation=tdi_gen, orthogonal=True),
89)
90
91tdi_data_tianqin = TDIChannelData()
92tdi_data_tianqin.set_fd_data_from_zero(
93 channels=tdi_chan,
94 duration=duration,
95 delta_time=dt,
96 start_time=t_start,
97 minimum_frequency=f_min,
98 maximum_frequency=f_max,
99 )
100tdi_data_tianqin.set_fd_noise_power_density_from_model("Tianqin_GWSpace", tdi_generation=tdi_gen)
101noise_tianqin = tdi_data_tianqin.get_fd_noise_realization()
102tdi_data_tianqin.add_into_fd_data(noise_tianqin)
103tianqin = InterferometerAntenna(
104 name="tianqin",
105 tdi_data=tdi_data_tianqin,
106 orbit_model=available_orbit_models['Tianqin_analytic'],
107 response_model=FDResponseModelMarset2018(),
108 tdi_combination=FDMichelsonConstantEqualArm(generation=tdi_gen, orthogonal=True),
109)
110
111common_freqs = tdi_data_lisa.data_info.frequency_samples_array
112wf_xhm = IMRPhenomXHM(common_freqs, f_ref, parameter_conversion=bilby.gw.conversion.generate_component_masses)
113wf_xhm.update_waveform(params)
114lisa.inject_signal(
115 wf_xhm.waveform_container,
116 params["ecliptic_longitude"],
117 params["ecliptic_latitude"],
118 params["polarization"],
119 params["coalescence_time"],
120)
121taiji.inject_signal(
122 wf_xhm.waveform_container,
123 params["ecliptic_longitude"],
124 params["ecliptic_latitude"],
125 params["polarization"],
126 params["coalescence_time"],
127)
128tianqin.inject_signal(
129 wf_xhm.waveform_container,
130 params["ecliptic_longitude"],
131 params["ecliptic_latitude"],
132 params["polarization"],
133 params["coalescence_time"],
134)
135
136########################################################################################
137# sampling settings
138likelihood = LikelihoodBilbyInterface(
139 waveform=wf_xhm,
140 detector=(lisa, taiji, tianqin,),
141 channels=('A', 'E', 'T')
142 )
143
144priors = {}
145priors['chi_1'] = bilby.core.prior.Uniform(name='chi_1', minimum=-0.99, maximum=0.99)
146priors['chi_2'] = bilby.core.prior.Uniform(name='chi_2', minimum=-0.99, maximum=0.99)
147priors['ecliptic_longitude'] = bilby.core.prior.Uniform(name='ecliptic_longitude', minimum=0, maximum=2 * lal.PI, boundary='periodic')
148priors['ecliptic_latitude'] = bilby.core.prior.Cosine(name='ecliptic_latitude')
149priors['inclination'] = bilby.core.prior.Sine(name='inclination')
150priors['polarization'] = bilby.core.prior.Uniform(name='polarization', minimum=0, maximum=lal.PI, boundary='periodic')
151priors['reference_phase'] = bilby.core.prior.Uniform(name='reference_phase', minimum=0, maximum=2 * lal.PI, boundary='periodic')
152priors['coalescence_time'] = bilby.core.prior.Uniform(name='coalescence_time', minimum=tc-100, maximum=tc+100)
153priors['chirp_mass'] = bilby.gw.prior.UniformInComponentsChirpMass(name='chirp_mass', minimum=5e5, maximum=2e6)
154priors['mass_ratio'] = bilby.gw.prior.UniformInComponentsMassRatio(name='mass_ratio', minimum=0.05, maximum=0.99)
155priors['luminosity_distance'] = bilby.core.prior.LogUniform(name='luminosity_distance', minimum=1e4, maximum=5e6)
156priors = bilby.core.prior.PriorDict(dictionary=priors)
157
158label = "inj_MBHB_network_HM"
159outdir = f"outdir_{label}"
160result = bilby.run_sampler(
161 likelihood=likelihood,
162 priors=priors,
163 outdir=outdir,
164 label=label,
165 sampler='pymultinest',
166 npoints=2048,
167 )