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

Total running time of the script: ( 0 minutes 0.000 seconds)

Gallery generated by Sphinx-Gallery