Note
Go to the end to download the full example code.
Randomly sampling parameter space points
Application of dedicated classes to sample software injection parameters according to the specified parameter space priors.
9 import os
10
11 import matplotlib.pyplot as plt
12 import numpy as np
13
14 from pyfstat import (
15 AllSkyInjectionParametersGenerator,
16 InjectionParametersGenerator,
17 Writer,
18 isotropic_amplitude_distribution,
19 set_up_logger,
20 )
21
22 label = "PyFstatExampleInjectionParametersGenerator"
23 outdir = os.path.join("PyFstat_example_data", label)
24 logger = set_up_logger(label=label, outdir=outdir)
25
26 # Properties of the GW data
27 gw_data = {
28 "sqrtSX": 1e-23,
29 "tstart": 1000000000,
30 "duration": 86400,
31 "detectors": "H1,L1",
32 "Band": 1,
33 "Tsft": 1800,
34 }
35
36 logger.info("Drawing random signal parameters...")
37
38 # Draw random signal phase parameters.
39 # The AllSkyInjectionParametersGenerator covers [Alpha,Delta] priors automatically.
40 # The rest can be a mix of nontrivial prior distributions and fixed values.
41 phase_params_generator = AllSkyInjectionParametersGenerator(
42 priors={
43 "F0": {"stats.uniform": {"loc": 29.0, "scale": 2.0}},
44 "F1": -1e-10,
45 "F2": 0,
46 },
47 seed=23,
48 )
49 phase_parameters = phase_params_generator.draw()
50 phase_parameters["tref"] = gw_data["tstart"]
51
52 # Draw random signal amplitude parameters.
53 # Here we use the plain InjectionParametersGenerator class.
54 amplitude_params_generator = InjectionParametersGenerator(
55 priors={
56 "h0": {"stats.norm": {"loc": 1e-24, "scale": 1e-26}},
57 **isotropic_amplitude_distribution,
58 },
59 seed=42,
60 )
61 amplitude_parameters = amplitude_params_generator.draw()
62
63 # Now we can pass the parameter dictionaries to the Writer class and make SFTs.
64 data = Writer(
65 label=label,
66 outdir=outdir,
67 **gw_data,
68 **phase_parameters,
69 **amplitude_parameters,
70 )
71 data.make_data()
72
73 # Now we draw many phase parameters and check the sky distribution
74 Ndraws = 10000
75 phase_parameters = phase_params_generator.draw_many(size=Ndraws)
76 Alphas = phase_parameters["Alpha"]
77 Deltas = phase_parameters["Delta"]
78 plotfile = os.path.join(outdir, label + "_allsky.png")
79 logger.info(f"Plotting sky distribution of {Ndraws} points to file: {plotfile}")
80 plt.subplot(111, projection="aitoff")
81 plt.plot(Alphas - np.pi, Deltas, ".", markersize=1)
82 plt.savefig(plotfile, dpi=300)
83 plt.close()
84 plotfile = os.path.join(outdir, label + "_alpha_hist.png")
85 logger.info(f"Plotting Alpha distribution of {Ndraws} points to file: {plotfile}")
86 plt.hist(Alphas, 50)
87 plt.xlabel("Alpha")
88 plt.ylabel("draws")
89 plt.savefig(plotfile, dpi=100)
90 plt.close()
91 plotfile = os.path.join(outdir, label + "_delta_hist.png")
92 logger.info(f"Plotting Delta distribution of {Ndraws} points to file: {plotfile}")
93 plt.hist(Deltas, 50)
94 plt.xlabel("Delta")
95 plt.ylabel("draws")
96 plt.savefig(plotfile, dpi=100)
97 plt.close()
98 plotfile = os.path.join(outdir, label + "_sindelta_hist.png")
99 logger.info(f"Plotting sin(Delta) distribution of {Ndraws} points to file: {plotfile}")
100 plt.hist(np.sin(Deltas), 50)
101 plt.xlabel("sin(Delta)")
102 plt.ylabel("draws")
103 plt.savefig(plotfile, dpi=100)
104 plt.close()