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