Note
Go to the end to download the full example code.
Short transient MCMC search
MCMC search for a Short transient CW signal.
8 import os
9
10 import numpy as np
11 import PyFstat_example_make_data_for_short_transient_search as data
12
13 import pyfstat
14 from pyfstat.utils import get_predict_fstat_parameters_from_dict
15
16 if __name__ == "__main__":
17 label = "PyFstatExampleShortTransientMCMCSearch"
18 logger = pyfstat.set_up_logger(label=label, outdir=data.outdir)
19
20 if not os.path.isdir(data.outdir) or not np.any(
21 [f.endswith(".sft") for f in os.listdir(data.outdir)]
22 ):
23 raise RuntimeError(
24 "Please first run PyFstat_example_make_data_for_short_transient_search.py !"
25 )
26
27 inj = {
28 "tref": data.tstart,
29 "F0": data.F0,
30 "F1": data.F1,
31 "F2": data.F2,
32 "Alpha": data.Alpha,
33 "Delta": data.Delta,
34 "transient_tstart": data.transient_tstart,
35 "transient_duration": data.transient_duration,
36 }
37
38 DeltaF0 = 1e-2
39 DeltaF1 = 1e-9
40
41 theta_prior = {
42 "F0": {
43 "type": "unif",
44 "lower": inj["F0"] - DeltaF0 / 2.0,
45 "upper": inj["F0"] + DeltaF0 / 2.0,
46 },
47 "F1": {
48 "type": "unif",
49 "lower": inj["F1"] - DeltaF1 / 2.0,
50 "upper": inj["F1"] + DeltaF1 / 2.0,
51 },
52 "F2": inj["F2"],
53 "Alpha": inj["Alpha"],
54 "Delta": inj["Delta"],
55 "transient_tstart": {
56 "type": "unif",
57 "lower": data.tstart,
58 "upper": data.tstart + data.duration - 2 * data.Tsft,
59 },
60 "transient_duration": {
61 "type": "unif",
62 "lower": 2 * data.Tsft,
63 "upper": data.duration - 2 * data.Tsft,
64 },
65 }
66
67 ntemps = 2
68 log10beta_min = -1
69 nwalkers = 100
70 nsteps = [200, 200]
71
72 BSGL = False
73 transientWindowType = "rect"
74
75 mcmc = pyfstat.MCMCTransientSearch(
76 label=label + ("BSGL" if BSGL else ""),
77 outdir=data.outdir,
78 sftfilepattern=os.path.join(data.outdir, f"*{data.label}*sft"),
79 theta_prior=theta_prior,
80 tref=inj["tref"],
81 nsteps=nsteps,
82 nwalkers=nwalkers,
83 ntemps=ntemps,
84 log10beta_min=log10beta_min,
85 transientWindowType=transientWindowType,
86 BSGL=BSGL,
87 )
88 mcmc.run(walker_plot_args={"plot_det_stat": True, "injection_parameters": inj})
89 mcmc.print_summary()
90 mcmc.plot_corner(add_prior=True, truths=inj)
91 mcmc.plot_prior_posterior(injection_parameters=inj)
92
93 # plot cumulative 2F, first building a dict as required for PredictFStat
94 d, maxtwoF = mcmc.get_max_twoF()
95 for key, val in mcmc.theta_prior.items():
96 if key not in d:
97 d[key] = val
98 d["h0"] = data.h0
99 d["cosi"] = data.cosi
100 d["psi"] = data.psi
101 PFS_input = get_predict_fstat_parameters_from_dict(
102 d, transientWindowType=transientWindowType
103 )
104 mcmc.plot_cumulative_max(PFS_input=PFS_input)