pylal.dq.noisebudget

1 # Copyright (C) 2011 Duncan Macleod 2 # 3 # This program is free software; you can redistribute it and/or modify it 4 # under the terms of the GNU General Public License as published by the 5 # Free Software Foundation; either version 3 of the License, or (at your 6 # option) any later version. 7 # 8 # This program is distributed in the hope that it will be useful, but 9 # WITHOUT ANY WARRANTY; without even the implied warranty of 10 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 11 # Public License for more details. 12 # 13 # You should have received a copy of the GNU General Public License along 14 # with this program; if not, write to the Free Software Foundation, Inc., 15 # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 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 # Preamble 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 # set metadata 38 __author__ = "Duncan M. Macleod <duncan.macleod@ligo.org>" 39 __version__ = git_version.id 40 __date__ = git_version.date 41 42 # ============================================================================= 43 # NoiseBudget class 44 # ============================================================================= 45

46 -class NoiseTerm(object):

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 # set types 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 # set defaults 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

85 - def __init__(self, **kwargs):

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 # define methods from seriesutils 103 # 104

105 - def fromNDS(self, *args, **kwargs):

106 self.timeseries = seriesutils.fromNDS(*args, **kwargs) 107 self.data = self.timeseries.data.data 108 self.deltaT = self.timeseries.deltaT 109 self.epoch = self.timeseries.epoch 110 self.f0 = self.timeseries.f0

111 fromNDS.__doc__ = seriesutils.fromNDS.__doc__ 112

113 - def fromcache(self, *args, **kwargs):

114 self.timeseries = seriesutils.fromlalcache(*args, **kwargs) 115 self.data = self.timeseries.data.data 116 self.deltaT = self.timeseries.deltaT 117 self.epoch = self.timeseries.epoch 118 self.f0 = self.timeseries.f0

119 fromcache.__doc__ = seriesutils.fromlalcache.__doc__ 120

121 - def fromframefile(self, *args, **kwargs):

122 self.timeseries = seriesutils.fromframefile(*args, **kwargs) 123 self.data = self.timeseries.data.data 124 self.deltaT = self.timeseries.deltaT 125 self.epoch = self.timeseries.epoch 126 self.f0 = self.timeseries.f0

127 fromframefile.__doc__ = seriesutils.fromframefile.__doc__ 128

129 - def compute_average_spectrum(self, *args, **kwargs):

130 self.frequencyseries =\ 131 seriesutils.compute_average_spectrum(*args, **kwargs) 132 self.spectrum = self.frequencyseries.data.data 133 self.epoch = self.frequencyseries.epoch 134 self.f0 = self.frequencyseries.f0 135 self.deltaF = self.frequencyseries.deltaF

136 compute_average_spectrum.__doc__=\ 137 seriesutils.compute_average_spectrum.__doc__ 138 139 # 140 # transfer functions and calibration 141 # 142

143 - def apply_calibration(self, func):

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

158 - def apply_spectrum_calibration(self, func):

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

178 - def loadreference(self, referencefile, epoch=seriesutils.lal.LIGOTimeGPS(),\ 179 sampleUnits=seriesutils.lal.lalStrainUnit, fcol=0,acol=1):

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 # restrict to frequency band 213 # 214

215 - def apply_frequency_band(self, fmin, fmax):

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 # plotter 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 # labels 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 # extract params 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 # build limits 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 # set plot object 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 # finalize plot 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 # save 288 plot.savefig(outfile, bbox_inches=bbox_inches,\ 289 bbox_extra_artists=plot.ax.texts)

290 291 # ============================================================================= 292 # NoiseBudget class 293 # ============================================================================= 294

295 -class NoiseBudget(list):

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 # init

313 - def __init__(self, *args, **kwargs):

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 # init 324 list.__init__([]) 325 326 # set arguments 327 for arg in args: 328 self.append(NoiseTerm(arg)) 329 330 # set keyword arguments 331 for arg,val in kwargs.items(): 332 setattr(self, arg, self.typemap[arg](val))

333 334 335 # append

336 - def append(self, *args, **kwargs):

337 list.append(self, *args, **kwargs) 338 self.numTerms += 1

339

340 - def sort(self, *args, **kwargs):

341 kwargs.setdefault("key", lambda t: t.name) 342 list.sort(self, *args, **kwargs)

343 344 # 345 # methods 346 # 347 348 # set target

349 - def set_target(self, noiseterm):

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

357 - def get_target(self):

358 """ 359 Returns the target of this NoiseBudget. Returns NoneType if no target 360 has been set. 361 """ 362 return self.target

363 364 # calculate sum of NoiseTerms

365 - def compute_noise_sum(self, name="Noise sum"):

366 """ 367 Calculate the quadrature sum of terms in the NoiseBudget. All terms whose 368 sum attribute evaluates to True will be included. 369 """ 370 # get metadata from first sum term 371 sum_terms = [t for t in self if t.sum] 372 if len(sum_terms): 373 t = sum_terms[0] 374 epoch = t.epoch 375 deltaF = t.deltaF 376 f0 = t.f0 377 else: 378 raise RuntimeError("No NoiseTerms were set to be included in the "+\ 379 "budget sum.") 380 self.noise_sum = NoiseTerm(name=name) 381 data = (numpy.asarray([t.spectrum for t in\ 382 sum_terms])**2).sum(axis=0)**(1/2) 383 self.noise_sum.frequencyseries =\ 384 seriesutils.fromarray(data, str(name), epoch=epoch,\ 385 deltaT=deltaF, f0=f0, frequencyseries=True) 386 self.noise_sum.fromarray

387 388 # compute difference between noise target and noise_sum

389 - def compute_deficit(self, func=lambda (s,t): abs(1-(t/s)**2)**(1/2)):

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 # compute difference between noise target and noise_sum

421 - def compute_ratio_deficit(self):

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 # plot

431 - def plot(self, outfile, **params):

432 433 """ 434 Plot this NoiseBudget 435 """ 436 437 # extract params 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 # build limits 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 # labels 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 # copy params for reference lines 461 refparams = copy.deepcopy(params) 462 refparams.setdefault('linestyle', '--') 463 refparams['linewidth'] = 0.3 464 465 # set plot object 466 plot = plotutils.DataPlot(xlabel, ylabel, title, subtitle) 467 reference_plotted = False 468 469 # plot target and noise sum 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 # plot reference 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 # plot terms in the budget 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 # finalize plot 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 # add text box 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 # save 521 plot.savefig(outfile, bbox_inches=bbox_inches,\ 522 bbox_extra_artists=plot.ax.texts)

523

524 - def plot_deficit(self, outfile, **params):

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

531 - def plot_ratio_deficit(self, outfile, **params):

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

Source Code for Module pylal.dq.noisebudget