Note
Click here 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.
8 import os
9
10 import matplotlib.pyplot as plt
11 import numpy as np
12
13 from pyfstat import (
14 AllSkyInjectionParametersGenerator,
15 InjectionParametersGenerator,
16 Writer,
17 )
18
19 label = "PyFstat_example_InjectionParametersGenerator"
20 outdir = os.path.join("PyFstat_example_data", label)
21
22 # Properties of the GW data
23 gw_data = {
24 "sqrtSX": 1e-23,
25 "tstart": 1000000000,
26 "duration": 86400,
27 "detectors": "H1,L1",
28 "Band": 1,
29 "Tsft": 1800,
30 }
31
32 print("Drawing random signal parameters...")
33
34 # Draw random signal phase parameters.
35 # The AllSkyInjectionParametersGenerator covers [Alpha,Delta] priors automatically.
36 # The rest can be a mix of nontrivial prior distributions and fixed values.
37 phase_params_generator = AllSkyInjectionParametersGenerator(
38 priors={
39 "F0": {"uniform": {"low": 29.0, "high": 31.0}},
40 "F1": -1e-10,
41 "F2": 0,
42 },
43 seed=23,
44 )
45 phase_parameters = phase_params_generator.draw()
46 phase_parameters["tref"] = gw_data["tstart"]
47
48 # Draw random signal amplitude parameters.
49 # Here we use the plain InjectionParametersGenerator class.
50 amplitude_params_generator = InjectionParametersGenerator(
51 priors={
52 "h0": {"normal": {"loc": 1e-24, "scale": 1e-24}},
53 "cosi": {"uniform": {"low": 0.0, "high": 1.0}},
54 "phi": {"uniform": {"low": 0.0, "high": 2 * np.pi}},
55 "psi": {"uniform": {"low": 0.0, "high": np.pi}},
56 },
57 seed=42,
58 )
59 amplitude_parameters = amplitude_params_generator.draw()
60
61 # Now we can pass the parameter dictionaries to the Writer class and make SFTs.
62 data = Writer(
63 label=label,
64 outdir=outdir,
65 **gw_data,
66 **phase_parameters,
67 **amplitude_parameters,
68 )
69 data.make_data()
70
71 # Now we draw many phase parameters and check the sky distribution
72 Ndraws = 10000
73 phase_parameters = [phase_params_generator.draw() for n in range(Ndraws)]
74 Alphas = np.array([p["Alpha"] for p in phase_parameters])
75 Deltas = np.array([p["Delta"] for p in phase_parameters])
76 plotfile = os.path.join(outdir, label + "_allsky.png")
77 print(f"Plotting sky distribution of {Ndraws} points to file: {plotfile}")
78 plt.subplot(111, projection="aitoff")
79 plt.plot(Alphas - np.pi, Deltas, ".", markersize=1)
80 plt.savefig(plotfile, dpi=300)
81 plt.close()
82 plotfile = os.path.join(outdir, label + "_alpha_hist.png")
83 print(f"Plotting Alpha distribution of {Ndraws} points to file: {plotfile}")
84 plt.hist(Alphas, 50)
85 plt.xlabel("Alpha")
86 plt.ylabel("draws")
87 plt.savefig(plotfile, dpi=100)
88 plt.close()
89 plotfile = os.path.join(outdir, label + "_delta_hist.png")
90 print(f"Plotting Delta distribution of {Ndraws} points to file: {plotfile}")
91 plt.hist(Deltas, 50)
92 plt.xlabel("Delta")
93 plt.ylabel("draws")
94 plt.savefig(plotfile, dpi=100)
95 plt.close()
96 plotfile = os.path.join(outdir, label + "_sindelta_hist.png")
97 print(f"Plotting sin(Delta) distribution of {Ndraws} points to file: {plotfile}")
98 plt.hist(np.sin(Deltas), 50)
99 plt.xlabel("sin(Delta)")
100 plt.ylabel("draws")
101 plt.savefig(plotfile, dpi=100)
102 plt.close()
Total running time of the script: ( 0 minutes 0.000 seconds)