Note
Go to the end to download the full example code.
Glitch robust MCMC search
MCMC search employing a signal hypothesis allowing for a glitch to be present in the data. The setup corresponds to a targeted search, and the simulated signal contains a single glitch.
10 import os
11 import time
12
13 import numpy as np
14 from PyFstat_example_make_data_for_search_on_1_glitch import (
15 F0,
16 F1,
17 F2,
18 Alpha,
19 Delta,
20 delta_F0,
21 dtglitch,
22 duration,
23 outdir,
24 tref,
25 tstart,
26 )
27
28 import pyfstat
29
30 label = "PyFstatExampleGlitchRobustDirectedMCMCSearchOn1Glitch"
31 logger = pyfstat.set_up_logger(label=label, outdir=outdir)
32
33 Nstar = 1000
34 F0_width = np.sqrt(Nstar) * np.sqrt(12) / (np.pi * duration)
35 F1_width = np.sqrt(Nstar) * np.sqrt(180) / (np.pi * duration**2)
36
37 theta_prior = {
38 "F0": {"type": "unif", "lower": F0 - F0_width / 2.0, "upper": F0 + F0_width / 2.0},
39 "F1": {"type": "unif", "lower": F1 - F1_width / 2.0, "upper": F1 + F1_width / 2.0},
40 "F2": F2,
41 "delta_F0": {"type": "unif", "lower": 0, "upper": 1e-5},
42 "delta_F1": 0,
43 "tglitch": {
44 "type": "unif",
45 "lower": tstart + 0.1 * duration,
46 "upper": tstart + 0.9 * duration,
47 },
48 "Alpha": Alpha,
49 "Delta": Delta,
50 }
51
52 ntemps = 3
53 log10beta_min = -0.5
54 nwalkers = 100
55 nsteps = [250, 250]
56
57 mcmc = pyfstat.MCMCGlitchSearch(
58 label=label,
59 outdir=outdir,
60 sftfilepattern=os.path.join(outdir, "*1glitch*sft"),
61 theta_prior=theta_prior,
62 tref=tref,
63 minStartTime=tstart,
64 maxStartTime=tstart + duration,
65 nsteps=nsteps,
66 nwalkers=nwalkers,
67 ntemps=ntemps,
68 log10beta_min=log10beta_min,
69 nglitch=1,
70 )
71 mcmc.transform_dictionary["F0"] = dict(
72 subtractor=F0, multiplier=1e6, symbol="$f-f_\\mathrm{s}$"
73 )
74 mcmc.unit_dictionary["F0"] = "$\\mu$Hz"
75 mcmc.transform_dictionary["F1"] = dict(
76 subtractor=F1, multiplier=1e12, symbol="$\\dot{f}-\\dot{f}_\\mathrm{s}$"
77 )
78 mcmc.unit_dictionary["F1"] = "$p$Hz/s"
79 mcmc.transform_dictionary["delta_F0"] = dict(
80 multiplier=1e6, subtractor=delta_F0, symbol="$\\delta f-\\delta f_\\mathrm{s}$"
81 )
82 mcmc.unit_dictionary["delta_F0"] = "$\\mu$Hz/s"
83 mcmc.transform_dictionary["tglitch"]["subtractor"] = tstart + dtglitch
84 mcmc.transform_dictionary["tglitch"][
85 "label"
86 ] = "$t^\\mathrm{g}-t^\\mathrm{g}_\\mathrm{s}$\n[d]"
87
88 t1 = time.time()
89 mcmc.run(save_loudest=False) # uses CFSv2 which doesn't support glitch parameters
90 dT = time.time() - t1
91 mcmc.print_summary()
92
93 logger.info("Making corner plot...")
94 mcmc.plot_corner(
95 label_offset=0.25,
96 truths={"F0": F0, "F1": F1, "delta_F0": delta_F0, "tglitch": tstart + dtglitch},
97 quantiles=(0.16, 0.5, 0.84),
98 hist_kwargs=dict(lw=1.5, zorder=-1),
99 truth_color="C3",
100 )
101
102 mcmc.plot_cumulative_max(savefig=True)
103
104 logger.info(f"Prior widths = {F0_width}, {F1_width}")
105 logger.info(f"Actual run time = {dT} s")