Coverage for pesummary/gw/file/calibration.py: 57.6%
59 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-05-02 08:42 +0000
1# Licensed under an MIT style license -- see LICENSE.md
3import numpy as np
4from scipy.interpolate import interp1d
5from pesummary import conf
6from pesummary.utils.utils import check_file_exists_and_rename
7from pesummary.utils.dict import Dict
9__author__ = ["Charlie Hoy <charlie.hoy@ligo.org>"]
12def _spline_angle_xform(delta_psi):
13 """Returns the angle in degrees corresponding to the spline
14 calibration parameters delta_psi. Code taken from lalinference.bayespputils
16 Parameters
17 ----------
18 delta_psi: array
19 calibration phase uncertainty
20 """
21 rotation = (2.0 + 1.0j * delta_psi) / (2.0 - 1.0j * delta_psi)
22 return 180.0 / np.pi * np.arctan2(np.imag(rotation), np.real(rotation))
25def _interpolate_spline_model(
26 frequencies, data, interpolated_frequencies, nfreqs=100, xform=None,
27 level=0.9, pbar=None
28):
29 """Interpolate calibration posterior estimates for a spline model in log
30 space. Code based upon same function in lalinference.bayespputils
32 Parameters
33 ----------
34 frequencies: array
35 The spline control points
36 data: ndarray
37 Array of posterior samples at each spline control point
38 interpolated_frequencies: array
39 Array of frequencies you wish to evaluate the interpolant for
40 nfreqs: int, optional
41 Number of points to evaluate the interpolates spline. Default 100
42 xform: func, optional
43 Function to transform the spline
44 """
45 interpolated_data = np.zeros((np.asarray(data).shape[0], nfreqs))
46 for num, samp in enumerate(data):
47 interp = interp1d(
48 frequencies, samp, kind="cubic", fill_value=0., bounds_error=False
49 )(interpolated_frequencies)
50 if xform is not None:
51 interp = xform(interp)
52 interpolated_data[num] = interp
53 if pbar is not None:
54 pbar.update(1)
56 mean = np.mean(interpolated_data, axis=0)
57 lower = np.quantile(interpolated_data, (1 - level) / 2., axis=0)
58 upper = np.quantile(interpolated_data, (1 + level) / 2., axis=0)
59 return mean, lower, upper
62def interpolate_calibration_posterior_from_samples(
63 log_frequencies, amplitudes, phases, nfreqs=100, level=0.9, **kwargs
64):
65 """Interpolate calibration posterior estimates for a spline model in log
66 space and return the amplitude and phase uncertainties. Code based upon same
67 function in lalinference.bayespputils
69 Parameters
70 ----------
71 log_frequencies: array
72 The spline control points.
73 amplitudes: ndarray
74 Array of amplitude posterior samples at each of the spline control
75 points
76 phases: ndarray
77 Array of phase posterior samples at each of the spline control points
78 nfreqs: int, optional
79 Number of points to evaluate the interpolates spline. Default 100
80 **kwargs: dict
81 All kwargs passed to _interpolate_spline_model
82 """
83 frequencies = np.exp(log_frequencies)
84 interpolated_frequencies = np.logspace(
85 np.min(log_frequencies), np.max(log_frequencies), nfreqs, base=np.e
86 )
87 amp_mean, amp_lower, amp_upper = (
88 1 + np.array(
89 _interpolate_spline_model(
90 log_frequencies, np.column_stack(amplitudes),
91 np.log(interpolated_frequencies), nfreqs=nfreqs, xform=None,
92 level=level, **kwargs
93 )
94 )
95 )
96 phase_mean, phase_lower, phase_upper = np.array(
97 _interpolate_spline_model(
98 log_frequencies, np.column_stack(phases),
99 np.log(interpolated_frequencies), nfreqs=nfreqs,
100 xform=_spline_angle_xform, level=level, **kwargs
101 )
102 ) * (np.pi / 180)
103 return np.column_stack(
104 [
105 interpolated_frequencies, amp_mean, phase_mean, amp_lower,
106 phase_lower, amp_upper, phase_upper
107 ]
108 )
111class CalibrationDict(Dict):
112 """Class to handle a dictionary of calibration data
114 Parameters
115 ----------
116 detectors: list
117 list of detectors
118 data: nd list
119 list of calibration samples for each detector. Each of the columns
120 should represent Frequency, Median Mag, Phase (Rad), -1 Sigma Mag,
121 -1 Sigma Phase, +1 Sigma Mag, +1 Sigma Phase
123 Attributes
124 ----------
125 detectors: list
126 list of detectors stored in the dictionary
128 Methods
129 -------
130 plot:
131 Generate a plot based on the calibration samples stored
132 """
133 def __init__(self, *args):
134 _columns = [
135 "frequencies", "magnitude", "phase", "magnitude_lower",
136 "phase_lower", "magnitude_upper", "phase_upper"
137 ]
138 super(CalibrationDict, self).__init__(
139 *args, value_class=Calibration, value_columns=_columns
140 )
142 @property
143 def detectors(self):
144 return list(self.keys())
147class Calibration(np.ndarray):
148 """Class to handle Calibration data
149 """
150 def __new__(cls, input_array):
151 obj = np.asarray(input_array).view(cls)
152 if obj.shape[1] != 7:
153 raise ValueError(
154 "Invalid input data. See the docs for instructions"
155 )
156 return obj
158 @classmethod
159 def read(cls, path_to_file, IFO=None, **kwargs):
160 """Read in a file and initialize the Calibration class
162 Parameters
163 ----------
164 path_to_file: str
165 the path to the file you wish to load
166 IFO: str, optional
167 name of the IFO which relates to the input file
168 **kwargs: dict
169 all kwargs are passed to the np.genfromtxt method
170 """
171 try:
172 f = np.genfromtxt(path_to_file, **kwargs)
173 return cls(f)
174 except Exception:
175 raise
177 @classmethod
178 def from_spline_posterior_samples(
179 cls, log_frequencies, amplitudes, phases, **kwargs
180 ):
181 """Interpolate calibration posterior estimates for a spline model in log
182 space and initialize the Calibration class
184 Parameters
185 ----------
186 """
187 samples = interpolate_calibration_posterior_from_samples(
188 log_frequencies, amplitudes, phases, level=0.68, nfreqs=300,
189 **kwargs
190 )
191 return cls(samples)
193 def save_to_file(self, file_name, comments="#", delimiter=conf.delimiter):
194 """Save the calibration data to file
196 Parameters
197 ----------
198 file_name: str
199 name of the file name that you wish to use
200 comments: str, optional
201 String that will be prepended to the header and footer strings, to
202 mark them as comments. Default is '#'.
203 delimiter: str, optional
204 String or character separating columns.
205 """
206 check_file_exists_and_rename(file_name)
207 header = [
208 "Frequency", "Median Mag", "Phase (Rad)", "-1 Sigma Mag",
209 "-1 Sigma Phase", "+1 Sigma Mag", "+1 Sigma Phase"
210 ]
211 np.savetxt(
212 file_name, self, delimiter=delimiter, comments=comments,
213 header=delimiter.join(header)
214 )
216 def __array_finalize__(self, obj):
217 if obj is None:
218 return