Note
Go to the end to download the full example code.
Follow up example
Multi-stage MCMC follow up of a CW signal produced by an isolated source using a ladder of coherent times.
9 import os
10
11 import matplotlib.pyplot as plt
12 import numpy as np
13
14 import pyfstat
15
16 label = "PyFstatExampleSemiCoherentDirectedFollowup"
17 outdir = os.path.join("PyFstat_example_data", label)
18 logger = pyfstat.set_up_logger(label=label, outdir=outdir)
19
20 # Properties of the GW data
21 data_parameters = {
22 "sqrtSX": 1e-23,
23 "tstart": 1000000000,
24 "duration": 100 * 86400,
25 "detectors": "H1",
26 }
27 tend = data_parameters["tstart"] + data_parameters["duration"]
28 mid_time = 0.5 * (data_parameters["tstart"] + tend)
29
30 # Properties of the signal
31 depth = 40
32 signal_parameters = {
33 "F0": 30.0,
34 "F1": -1e-10,
35 "F2": 0,
36 "Alpha": np.radians(83.6292),
37 "Delta": np.radians(22.0144),
38 "tref": mid_time,
39 "h0": data_parameters["sqrtSX"] / depth,
40 "cosi": 1.0,
41 }
42
43 data = pyfstat.Writer(
44 label=label, outdir=outdir, **data_parameters, **signal_parameters
45 )
46 data.make_data()
47
48 # The predicted twoF (expectation over noise realizations) can be accessed by
49 twoF = data.predict_fstat()
50 logger.info("Predicted twoF value: {}\n".format(twoF))
51
52 # Search
53 VF0 = VF1 = 1e5
54 DeltaF0 = np.sqrt(VF0) * np.sqrt(3) / (np.pi * data_parameters["duration"])
55 DeltaF1 = np.sqrt(VF1) * np.sqrt(180) / (np.pi * data_parameters["duration"] ** 2)
56 theta_prior = {
57 "F0": {
58 "type": "unif",
59 "lower": signal_parameters["F0"] - DeltaF0 / 2.0,
60 "upper": signal_parameters["F0"] + DeltaF0 / 2,
61 },
62 "F1": {
63 "type": "unif",
64 "lower": signal_parameters["F1"] - DeltaF1 / 2.0,
65 "upper": signal_parameters["F1"] + DeltaF1 / 2,
66 },
67 }
68 for key in "F2", "Alpha", "Delta":
69 theta_prior[key] = signal_parameters[key]
70
71
72 ntemps = 3
73 log10beta_min = -0.5
74 nwalkers = 100
75 nsteps = [100, 100]
76
77 mcmc = pyfstat.MCMCFollowUpSearch(
78 label=label,
79 outdir=outdir,
80 sftfilepattern=data.sftfilepath,
81 theta_prior=theta_prior,
82 tref=mid_time,
83 minStartTime=data_parameters["tstart"],
84 maxStartTime=tend,
85 nwalkers=nwalkers,
86 nsteps=nsteps,
87 ntemps=ntemps,
88 log10beta_min=log10beta_min,
89 )
90
91 NstarMax = 1000
92 Nsegs0 = 100
93 walkers_fig, walkers_axes = plt.subplots(nrows=2, figsize=(3.4, 3.5))
94 mcmc.run(
95 NstarMax=NstarMax,
96 Nsegs0=Nsegs0,
97 plot_walkers=True,
98 walker_plot_args={
99 "labelpad": 0.01,
100 "plot_det_stat": False,
101 "fig": walkers_fig,
102 "axes": walkers_axes,
103 "injection_parameters": signal_parameters,
104 },
105 )
106 walkers_fig.savefig(os.path.join(outdir, label + "_walkers.png"))
107 plt.close(walkers_fig)
108
109 mcmc.print_summary()
110 mcmc.plot_corner(add_prior=True, truths=signal_parameters)
111 mcmc.plot_prior_posterior(injection_parameters=signal_parameters)