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 set_up_logger,
18 )
19
20 label = "PyFstat_example_InjectionParametersGenerator"
21 outdir = os.path.join("PyFstat_example_data", label)
22 logger = set_up_logger(label=label, outdir=outdir)
23
24 # Properties of the GW data
25 gw_data = {
26 "sqrtSX": 1e-23,
27 "tstart": 1000000000,
28 "duration": 86400,
29 "detectors": "H1,L1",
30 "Band": 1,
31 "Tsft": 1800,
32 }
33
34 logger.info("Drawing random signal parameters...")
35
36 # Draw random signal phase parameters.
37 # The AllSkyInjectionParametersGenerator covers [Alpha,Delta] priors automatically.
38 # The rest can be a mix of nontrivial prior distributions and fixed values.
39 phase_params_generator = AllSkyInjectionParametersGenerator(
40 priors={
41 "F0": {"uniform": {"low": 29.0, "high": 31.0}},
42 "F1": -1e-10,
43 "F2": 0,
44 },
45 seed=23,
46 )
47 phase_parameters = phase_params_generator.draw()
48 phase_parameters["tref"] = gw_data["tstart"]
49
50 # Draw random signal amplitude parameters.
51 # Here we use the plain InjectionParametersGenerator class.
52 amplitude_params_generator = InjectionParametersGenerator(
53 priors={
54 "h0": {"normal": {"loc": 1e-24, "scale": 1e-24}},
55 "cosi": {"uniform": {"low": 0.0, "high": 1.0}},
56 "phi": {"uniform": {"low": 0.0, "high": 2 * np.pi}},
57 "psi": {"uniform": {"low": 0.0, "high": np.pi}},
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() for n in range(Ndraws)]
76 Alphas = np.array([p["Alpha"] for p in phase_parameters])
77 Deltas = np.array([p["Delta"] for p in phase_parameters])
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()
Total running time of the script: ( 0 minutes 0.000 seconds)