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