Note
Go to the end to download the full example code.
Directed grid search: Linear spindown
Search for CW signal including one spindown parameter using a parameter space grid (i.e. no MCMC).
9 import os
10
11 import numpy as np
12
13 import pyfstat
14
15 label = "PyFstatExampleGridSearchF0F1"
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 sqrtSX = 1e-23
21 tstart = 1000000000
22 duration = 10 * 86400
23 tend = tstart + duration
24 tref = 0.5 * (tstart + tend)
25 IFOs = "H1"
26
27 # parameters for injected signals
28 depth = 20
29 inj = {
30 "tref": tref,
31 "F0": 30.0,
32 "F1": -1e-10,
33 "F2": 0,
34 "Alpha": 1.0,
35 "Delta": 1.5,
36 "h0": sqrtSX / depth,
37 "cosi": 0.0,
38 }
39
40 data = pyfstat.Writer(
41 label=label,
42 outdir=outdir,
43 tstart=tstart,
44 duration=duration,
45 sqrtSX=sqrtSX,
46 detectors=IFOs,
47 **inj,
48 )
49 data.make_data()
50
51 m = 0.01
52 dF0 = np.sqrt(12 * m) / (np.pi * duration)
53 dF1 = np.sqrt(180 * m) / (np.pi * duration**2)
54 dF2 = 1e-17
55 N = 100
56 DeltaF0 = N * dF0
57 DeltaF1 = N * dF1
58 F0s = [inj["F0"] - DeltaF0 / 2.0, inj["F0"] + DeltaF0 / 2.0, dF0]
59 F1s = [inj["F1"] - DeltaF1 / 2.0, inj["F1"] + DeltaF1 / 2.0, dF1]
60 F2s = [inj["F2"]]
61 Alphas = [inj["Alpha"]]
62 Deltas = [inj["Delta"]]
63 search = pyfstat.GridSearch(
64 label=label,
65 outdir=outdir,
66 sftfilepattern=data.sftfilepath,
67 F0s=F0s,
68 F1s=F1s,
69 F2s=F2s,
70 Alphas=Alphas,
71 Deltas=Deltas,
72 tref=tref,
73 minStartTime=tstart,
74 maxStartTime=tend,
75 )
76 search.run()
77
78 # report details of the maximum point
79 max_dict = search.get_max_twoF()
80 logger.info(
81 "max2F={:.4f} from GridSearch, offsets from injection: {:s}.".format(
82 max_dict["twoF"],
83 ", ".join(
84 [
85 "{:.4e} in {:s}".format(max_dict[key] - inj[key], key)
86 for key in max_dict.keys()
87 if not key == "twoF"
88 ]
89 ),
90 )
91 )
92 search.generate_loudest()
93
94 logger.info("Plotting 2F(F0)...")
95 search.plot_1D(xkey="F0", xlabel="freq [Hz]", ylabel="$2\\mathcal{F}$")
96 logger.info("Plotting 2F(F1)...")
97 search.plot_1D(xkey="F1")
98 logger.info("Plotting 2F(F0,F1)...")
99 search.plot_2D(xkey="F0", ykey="F1", colorbar=True)
100
101 logger.info("Making gridcorner plot...")
102 F0_vals = np.unique(search.data["F0"]) - inj["F0"]
103 F1_vals = np.unique(search.data["F1"]) - inj["F1"]
104 twoF = search.data["twoF"].reshape((len(F0_vals), len(F1_vals)))
105 xyz = [F0_vals, F1_vals]
106 labels = [
107 "$f - f_0$",
108 "$\\dot{f} - \\dot{f}_0$",
109 "$\\widetilde{2\\mathcal{F}}$",
110 ]
111 fig, axes = pyfstat.gridcorner(
112 twoF, xyz, projection="log_mean", labels=labels, whspace=0.1, factor=1.8
113 )
114 fig.savefig(os.path.join(outdir, label + "_projection_matrix.png"))