BasicTutorial.py 3.36 KB
Newer Older
plasky's avatar
bluh    
plasky committed
1
import numpy as np
plasky's avatar
plasky committed
2
import pylab as plt
3

4
import dynesty.plotting as dyplot
5
import corner
6
import peyote
plasky's avatar
bluh    
plasky committed
7

8
peyote.utils.setup_logger()
9

plasky's avatar
plasky committed
10
time_duration = 1.
11
12
sampling_frequency = 4096.

13
simulation_parameters = dict(
14
15
    mass_1=36.,
    mass_2=29.,
16
17
18
19
20
21
    spin11=0,
    spin12=0,
    spin13=0,
    spin21=0,
    spin22=0,
    spin23=0,
Gregory Ashton's avatar
Gregory Ashton committed
22
    luminosity_distance=100.,
23
    iota=0.4,
Gregory Ashton's avatar
Gregory Ashton committed
24
25
26
27
28
29
    phase=1.3,
    waveform_approximant='IMRPhenomPv2',
    reference_frequency=50.,
    ra=1.375,
    dec=-1.2108,
    geocent_time=1126259642.413,
30
    psi=2.659
31
)
plasky's avatar
bluh    
plasky committed
32

moritz's avatar
moritz committed
33
# Create the waveform_generator using a LAL BinaryBlackHole source function
34
waveform_generator = peyote.waveform_generator.WaveformGenerator(
35
36
37
    sampling_frequency=sampling_frequency,
    time_duration=time_duration,
    source_model=peyote.source.lal_binary_black_hole)
38
waveform_generator.parameters = simulation_parameters
39
hf_signal = waveform_generator.frequency_domain_strain()
Gregory Ashton's avatar
Gregory Ashton committed
40

41
# Simulate the data in H1
42
43
H1 = peyote.detector.H1
H1_hf_noise, frequencies = H1.power_spectral_density.get_noise_realisation(
44
    sampling_frequency, time_duration)
45
46
H1.set_data(sampling_frequency, time_duration,
            frequency_domain_strain=H1_hf_noise)
47
H1.inject_signal(waveform_generator)
48
49
50
51
52
53
H1.set_spectral_densities()

# Simulate the data in L1
L1 = peyote.detector.L1
L1_hf_noise, frequencies = L1.power_spectral_density.get_noise_realisation(
    sampling_frequency, time_duration)
54
55
L1.set_data(sampling_frequency, time_duration,
            frequency_domain_strain=L1_hf_noise)
56
L1.inject_signal(waveform_generator)
57
58
59
60
L1.set_spectral_densities()

IFOs = [H1, L1]

61
# Plot the noise and signal
Gregory Ashton's avatar
Gregory Ashton committed
62
63
64
65
66
67
68
69
70
fig, ax = plt.subplots()
plt.loglog(frequencies, np.abs(H1_hf_noise), lw=1.5, label='H1 noise+signal')
plt.loglog(frequencies, np.abs(L1_hf_noise), lw=1.5, label='L1 noise+signal')
plt.loglog(frequencies, np.abs(hf_signal['plus']), lw=0.8, label='signal')
plt.xlim(10, 1000)
plt.legend()
plt.xlabel(r'frequency')
plt.ylabel(r'strain')
fig.savefig('data')
71

72
likelihood = peyote.likelihood.Likelihood(IFOs, waveform_generator)
Gregory Ashton's avatar
Gregory Ashton committed
73

74
# New way way of doing it, still not perfect
75
76
77
78
79
simulation_parameters = peyote.parameter.Parameter.parse_floats_to_parameters(simulation_parameters)
simulation_parameters['mass_1'].prior = peyote.prior.Uniform(lower=35, upper=37)
simulation_parameters['mass_1'].is_fixed = False
simulation_parameters['luminosity_distance'].prior = peyote.prior.Uniform(lower=30, upper=200)
simulation_parameters['luminosity_distance'].is_fixed = False
80
#waveform_generator.set_values(simulation_parameters)
81
result = peyote.sampler.run_sampler(likelihood, simulation_parameters, sampler='nestle', verbose=True)
82
83
truths = [simulation_parameters[x].value for x in result.search_parameter_keys]

84
85
86
87
88
89
90
91
# Old way of doing it, still works
# prior = simulation_parameters.copy()
# prior['mass_1'] = peyote.parameter.Parameter(
#     'mass_1', prior=peyote.prior.Uniform(lower=35, upper=37),
#     latex_label='$m_1$')
# prior['luminosity_distance'] = peyote.parameter.Parameter(
#     'luminosity_distance', prior=peyote.prior.Uniform(lower=30, upper=200),
#     latex_label='$d_L$')
92
93
# result = peyote.run_sampler(likelihood, prior, sampler='nestle', verbose=True)
# truths = [simulation_parameters[x] for x in result.search_parameter_keys]
Gregory Ashton's avatar
Gregory Ashton committed
94

Gregory Ashton's avatar
Gregory Ashton committed
95

96
fig = corner.corner(result.samples, truths=truths, labels=result.search_parameter_keys)
Gregory Ashton's avatar
Gregory Ashton committed
97
98
99
100
fig.savefig('corner')

fig, axes = dyplot.traceplot(result['sampler_output'])
fig.savefig('trace')