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