1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 """
18 Class definitions and helper functions for building a GW
19 interferometer or system noise budget.
20 Bugs/change requests: https://bugs.ligo.org/redmine/projects/detchar/issues/new
21 """
22
23
24
25
26
27 from __future__ import division
28
29 import numpy
30 import re
31 import copy
32
33 from pylal import seriesutils
34 from pylal import plotutils
35 from pylal import git_version
36
37
38 __author__ = "Duncan M. Macleod <duncan.macleod@ligo.org>"
39 __version__ = git_version.id
40 __date__ = git_version.date
41
42
43
44
45
47 """
48 Object representing one term in a noise budget/projection. It holds a few
49 parameters and has methods to read data and generate ASD spectrum borrowed
50 from pylal.seriesutils, and simple plotting from pylal.plotutils
51 """
52
53 typemap = {\
54 'channel': str,\
55 'frame_type': str,\
56 'name': str,\
57 'sum': bool,\
58 'data': numpy.array,\
59 'deltaT': float,\
60 'spectrum': numpy.array,\
61 'deltaF': float,\
62 'f0': float,\
63 'epoch': seriesutils.lal.LIGOTimeGPS,\
64 'ref_spectrum': numpy.array,\
65 'ref_deltaF': float,\
66 'ref_f0': float,\
67 'ref_epoch': seriesutils.lal.LIGOTimeGPS\
68 }
69
70 name = None
71 channel = None
72 frame_type = None
73 sum = False
74 data = None
75 deltaT = 0
76 spectrum = None
77 deltaF = 0
78 f0 = 0
79 epoch = seriesutils.lal.LIGOTimeGPS()
80 ref_spectrum = None
81 ref_deltaF = 0
82 ref_f0 = 0
83 ref_epoch = seriesutils.lal.LIGOTimeGPS()
84
86 """"
87 Initialise a NoiseTerm.
88
89 Keyword arguments:
90
91 name : str
92 description of this NoiseTerm
93 channel : str
94 data channel from which this noise term is estimated
95 frame_type : str
96 LDR frame type for files containing data for this NoiseTerm
97 """
98 for arg,val in kwargs.items():
99 setattr(self, arg, self.typemap[arg](val))
100
101
102
103
104
105 - def fromNDS(self, *args, **kwargs):
111 fromNDS.__doc__ = seriesutils.fromNDS.__doc__
112
119 fromcache.__doc__ = seriesutils.fromlalcache.__doc__
120
127 fromframefile.__doc__ = seriesutils.fromframefile.__doc__
128
136 compute_average_spectrum.__doc__=\
137 seriesutils.compute_average_spectrum.__doc__
138
139
140
141
142
144 """
145 Apply a calibration function to the timeseries.
146
147 Arguments:
148
149 func : callable
150 any callable function that accepts numpy.array as input
151 """
152 if not hasattr(self, "data"):
153 raise RuntimeError("No data has been loaded for this channel. "+\
154 "Please use one of the data getting methods f"+\
155 "or this NoiseTerm before applying calibration.")
156 self.data = func(self.data)
157
159 """
160 Apply the transfer function to the spectrum.
161
162 Arguments:
163
164 func : [ callable | numpy.array ]
165 any callable function that accepts numpy.array as input or
166 array by which to multiply spectrum.
167 """
168 if not hasattr(self, "spectrum"):
169 raise RuntimeError("No spectrum has been loaded for this channel."+\
170 " Please use one of the spectrum getting "+\
171 "methods for this NoiseTerm before applying "+\
172 "calibration.")
173 if hasattr(func, "shape"):
174 self.spectrum *= func
175 else:
176 self.spectrum = func(self.spectrum)
177
180 """
181 Read the reference spectrum from a two-column file of frequency and ASD.
182
183 Arguments:
184
185 referencefile : str
186 path to file with space-delimited columns of frequency and ASD
187
188 Keyword arguments:
189
190 epoch : lal.LIGOTimeGPS
191 GPS epoch of this reference
192 sampleUnits : lal.LALUnit
193 amplitude unit for reference spectrum
194 fcol : int
195 number of column holding frequency array, default 0
196 acol : int
197 number of column holding amplitude array, default 1
198 """
199 f, data = numpy.loadtxt(referencefile, dtype=float, unpack=True,\
200 usecols=[fcol, acol])
201 deltaF = f[1]-f[0]
202 f0 = f[0]
203 self.ref_frequencyseries = seriesutils.fromarray(data, str(self.name),
204 epoch=epoch,\
205 deltaT=deltaF, f0=f0,\
206 frequencyseries=True)
207 self.ref_spectrum = data
208 self.ref_f0 = f0
209 self.ref_deltaF = deltaF
210
211
212
213
214
216 """
217 Restrict the spectrum for this NoiseTerm to the given semiopen
218 [fmin, fmax) interval.
219
220 Arguments:
221
222 fmin : float
223 minimum frequency for this NoiseTerm
224 fmax : float
225 maximum frequency for this NoiseTerm
226 """
227 f = numpy.arange(len(self.spectrum))*self.deltaF + self.f0
228 numpy.putmask(self.spectrum, ((f<flim[0])|(f>=flim[1])), 0)
229
230
231
232
233
234 - def plot(self, outfile, **params):
235 """
236 Plot the spectrum of this NoiseTerm.
237
238 Arguments:
239
240 outfile : str
241 path to output file for this plot
242 """
243
244 xlabel = params.pop("xlabel", "Frequency (Hz)")
245 ylabel = params.pop("ylabel", "Strain amplitude spectral density "
246 "$/\sqrt{\mbox{Hz}}$")
247 title = params.pop("title", "%s noise curve"\
248 % plotutils.display_name(self.name))
249 subtitle = params.pop("subtitle", "")
250
251
252 bbox_inches = params.pop("bbox_inches", "tight")
253 hidden_colorbar = params.pop("hidden_colorbar", False)
254 logx = params.pop("logx", True)
255 logy = params.pop("logy", True)
256
257
258 minx = self.deltaF
259 maxx = self.spectrum.size*self.deltaF + self.f0
260 xlim = params.pop("xlim", [minx, maxx])
261 ylim = params.pop("ylim", None)
262
263
264 plot = plotutils.DataPlot(xlabel, ylabel, title, subtitle)
265
266 f = numpy.arange(self.spectrum.size) * self.deltaF +\
267 self.f0
268 plot.add_content(f, self.spectrum,\
269 label=plotutils.display_name(self.name),\
270 color=self.color, **params)
271
272 if self.ref_spectrum is not None:
273 params["linestyle"] = "--"
274 params['linewidth'] = 0.3
275 f = numpy.arange(self.ref_spectrum.size) *\
276 self.ref_frequencyseries.deltaF + self.frequencyseries.f0
277 plot.add_content(f, self.ref_spectrum,\
278 label=plotutils.display_name(self.name),\
279 color=self.color, **params)
280
281
282 plot.finalize(logx=logx, logy=logy, loc='upper right',\
283 hidden_colorbar=hidden_colorbar)
284 if xlim: plot.ax.set_xlim(xlim)
285 if ylim: plot.ax.set_ylim(ylim)
286
287
288 plot.savefig(outfile, bbox_inches=bbox_inches,\
289 bbox_extra_artists=plot.ax.texts)
290
291
292
293
294
296 """
297 Object representing a list of NoiseTerms comprising a noise budget
298 estimation.
299 """
300 typemap = {\
301 'name': str,\
302 'numTerms': int,\
303 'target': NoiseTerm,\
304 'noise_sum': NoiseTerm\
305 }
306
307 name = None
308 numTerms = 0
309 target = None
310 noise_sum = None
311
312
314 """
315 Initialise this NoiseBudget. Arguments should be NoiseTerms to add to the
316 NoiseBudget.
317
318 Keyword arguments:
319
320 name : str
321 name for this NoiseBudget
322 """
323
324 list.__init__([])
325
326
327 for arg in args:
328 self.append(NoiseTerm(arg))
329
330
331 for arg,val in kwargs.items():
332 setattr(self, arg, self.typemap[arg](val))
333
334
335
336 - def append(self, *args, **kwargs):
339
340 - def sort(self, *args, **kwargs):
341 kwargs.setdefault("key", lambda t: t.name)
342 list.sort(self, *args, **kwargs)
343
344
345
346
347
348
350 """
351 Set the target of this NoiseBudget to be the noiseterm. Argument should be
352 a NoiseTerm object.
353 """
354 if not isinstance(noiseterm, NoiseTerm):
355 raise TypeError("Target must be a NoiseTerm.")
356 self.target = noiseterm
358 """
359 Returns the target of this NoiseBudget. Returns NoneType if no target
360 has been set.
361 """
362 return self.target
363
364
387
388
390 """
391 Calculate the deficit of thise NoiseBudget as the normalised quadrature
392 difference ratio the sum of the NoiseTerms and the target. Defaults to:
393
394 deficit = $\sqrt{1 - \left(\frac{target}{sum of terms}\right)^2}
395
396 Keyword arguments:
397
398 func : callable
399 any callable function that accepts noise sum and target arrays
400 """
401 if self.target is None:
402 raise RuntimeError("NoiseBudget target has not been set. "+\
403 "Run NoiseBudget.set_target(mytarget).")
404 if self.noise_sum is None:
405 raise RuntimeError("NoiseBudget sum has not be computed. "+\
406 "Run NoiseBudget.compute_noise_sum().")
407 deficit = NoiseTerm('%d budget deficit' % self.name)
408 deficit.epoch = self.target.epoch
409 deficit.deltaF = self.target.deltaF
410 deficit.f0 = self.target.f0
411 data = func(self.noise_sum.spectrum, self.target.spectrum)
412 deficit.frequencyseries = seriesutils.fromarray(data,\
413 name=deficit.name,\
414 epoch=deficit.epoch,\
415 deltaT=deficit.deltaF,\
416 f0=deficit.f0,\
417 frequencyseries=True)
418 return deficit
419
420
422 """
423 Calculate the deficit of thise NoiseBudget as the ratio of the sum of the
424 NoiseTerms and the target.
425
426 ratio_deficit = $\frac{target}{sum of terms}
427 """
428 return self.compute_deficit(func=lambda (s,t): t/s)
429
430
431 - def plot(self, outfile, **params):
432
433 """
434 Plot this NoiseBudget
435 """
436
437
438 bbox_inches = params.pop("bbox_inches", "tight")
439 hidden_colorbar = params.pop("hidden_colorbar", False)
440 logx = params.pop("logx", True)
441 logy = params.pop("logy", True)
442
443 plot_component_ref = params.pop("plot_component_ref", False)
444 if plot_component_ref == 'true':
445 plot_component_ref = True
446
447
448 minx = min(t.deltaF for t in self)
449 maxx = max(t.spectrum.size*t.deltaF + t.f0 for t in self)
450 xlim = params.pop("xlim", [minx, maxx])
451 ylim = params.pop("ylim", None)
452
453
454 xlabel = params.pop("xlabel", "Frequency (Hz)")
455 ylabel = params.pop("ylabel", "Strain amplitude spectral density "
456 "$/\sqrt{\mbox{Hz}}$")
457 title = params.pop("title", plotutils.display_name(self.name))
458 subtitle = params.pop("subtitle", "")
459
460
461 refparams = copy.deepcopy(params)
462 refparams.setdefault('linestyle', '--')
463 refparams['linewidth'] = 0.3
464
465
466 plot = plotutils.DataPlot(xlabel, ylabel, title, subtitle)
467 reference_plotted = False
468
469
470 for term in ["target", "noise_sum"]:
471 if hasattr(self, term) and getattr(self, term) is not None:
472 term = getattr(self, term)
473 f = numpy.arange(term.data.size) *term.frequencyseries.deltaF +\
474 term.frequencyseries.f0
475
476 plot.add_content(f, term.data,\
477 label=plotutils.display_name(term.name),\
478 linecolor=term.color, **params)
479
480
481 if term.ref_spectrum is not None:
482 f = numpy.arange(term.ref_spectrum.size) *\
483 self.ref_frequencyseries.deltaF +\
484 self.ref_frequencyeries.f0
485 plot.add_content(f, term.ref_spectrum,\
486 label=plotutils.display_name(self.name),\
487 color=self.color, **refparams)
488 reference_plotted = True
489
490
491 for term in self:
492 f = numpy.arange(term.spectrum.size) * term.deltaF +\
493 term.f0
494 plot.add_content(f, term.spectrum,\
495 label=plotutils.display_name(term.name),\
496 color=term.color, **params)
497 if plot_component_ref and term.ref_spectrum is not None:
498 f = numpy.arange(term.ref_spectrum.size) *\
499 term.ref_frequencyseries.deltaF +\
500 term.ref_frequencyseries.f0
501 plot.add_content(f, term.ref_spectrum,\
502 label=None, color=term.color, **refparams)
503 reference_plotted = True
504
505
506 plot.finalize(logx=logx, logy=logy, loc='upper right',\
507 hidden_colorbar=hidden_colorbar)
508 if xlim: plot.ax.set_xlim(xlim)
509 if ylim: plot.ax.set_ylim(ylim)
510
511
512 if reference_plotted:
513 plot.ax.text(0.005, 0.99, 'Dashed line indicates reference',\
514 horizontalalignment='left', verticalalignment='top',\
515 transform=plot.ax.transAxes,\
516 backgroundcolor='white',\
517 bbox=dict(facecolor='white', alpha=0.6,\
518 edgecolor='black'))
519
520
521 plot.savefig(outfile, bbox_inches=bbox_inches,\
522 bbox_extra_artists=plot.ax.texts)
523
525 """
526 Plot the difference between the noise budget sum and its target.
527 """
528 deficit = self.compute_deficit()
529 deficit.plot(outfile, **params)
530
532 """
533 Plot the difference between the noise budget sum and its target.
534 """
535 deficit = compute_deficit_ratio()
536 deficit.plot(outfile, **params)
537