pylal.dq.dqPlotUtils

1 #!/usr/bin/env python 2 3 # Copyright (C) 2011 Duncan Macleod 4 # 5 # This program is free software; you can redistribute it and/or modify it 6 # under the terms of the GNU General Public License as published by the 7 # Free Software Foundation; either version 3 of the License, or (at your 8 # option) any later version. 9 # 10 # This program is distributed in the hope that it will be useful, but 11 # WITHOUT ANY WARRANTY; without even the implied warranty of 12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 13 # Public License for more details. 14 # 15 # You should have received a copy of the GNU General Public License along 16 # with this program; if not, write to the Free Software Foundation, Inc., 17 # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 18 19 # ============================================================================= 20 # Preamble 21 # ============================================================================= 22 23 from __future__ import division 24 import math,re,numpy,itertools,copy,matplotlib,sys,warnings 25 26 # test matplotlib backend and reset if needed 27 from os import getenv 28 _display = getenv('DISPLAY','') 29 _backend_warn = """No display detected, moving to 'Agg' backend in matplotlib.""" 30 if not _display and not matplotlib.get_backend().lower() == 'agg': 31 warnings.warn(_backend_warn) 32 matplotlib.use('Agg', warn=False) 33 import pylab 34 35 try: from mpl_toolkits.basemap import Basemap 36 except ImportError,e: warnings.warn(str(e)) 37 try: from mpl_toolkits import axes_grid 38 except ImportError,e: warnings.warn(str(e)) 39 40 from datetime import datetime 41 from glue import segments,git_version 42 from pylal.xlal.datatypes.ligotimegps import LIGOTimeGPS 43 from lalburst import date 44 from pylal import plotutils 45 from pylal.dq.dqTriggerUtils import def_get_time,get_column 46 47 __author__ = "Duncan Macleod <duncan.macleod@astro.cf.ac.uk>" 48 __version__ = "git id %s" % git_version.id 49 __date__ = git_version.date 50 51 """ 52 This module provides plotting routines for use in data quality investigations. All routines are written to work in as general a way as possible with ligolw tables and lsctables compatible columns, and to plot in a similar pythonic way to pylal.plotutils. 53 """

54 55 # ============================================================================= 56 # Set plot parameters aux helper functions 57 # ============================================================================= 58 59 -def set_rcParams():

60 61 """ 62 Customise the figure parameters. 63 """ 64 65 # customise plot appearance 66 pylab.rcParams.update({"text.usetex": True, 67 "text.verticalalignment": "center", 68 "lines.linewidth": 2, 69 "xtick.labelsize": 16, 70 "ytick.labelsize": 16, 71 "axes.titlesize": 22, 72 "axes.labelsize": 16, 73 "axes.linewidth": 1, 74 "grid.linewidth": 1, 75 "legend.fontsize": 16, 76 "legend.loc": "best", 77 "figure.figsize": [12,6], 78 "figure.dpi": 80, 79 "image.origin": 'lower', 80 "axes.grid": True, 81 "axes.axisbelow": False })

82

83 -def set_ticks(ax, calendar_time=False):

84 85 """ 86 Format the x- and y-axis ticks to ensure minor ticks appear when needed 87 and the x-axis is set for spaces of 4 rather than 5. 88 89 Arguments: 90 91 ax : matplotlib.axes.AxesSubplot 92 Axes object to format 93 """ 94 95 # 96 # make sure we get minor ticks if there are no major ticks in the range 97 # 98 99 # xticks 100 ticks = list(ax.get_xticks()) 101 xlim = ax.get_xlim() 102 for i,tick in enumerate(ticks[::-1]): 103 if not xlim[0] <= tick <= xlim[1]: ticks.pop(-1) 104 if len(ticks)<=1: 105 ax.xaxis.set_minor_formatter(pylab.matplotlib.ticker.ScalarFormatter()) 106 107 # yticks 108 ticks = list(ax.get_yticks()) 109 ylim = ax.get_ylim() 110 for i,tick in enumerate(ticks[::-1]): 111 if not ylim[0] <= tick <= ylim[1]: ticks.pop(-1) 112 if len(ticks)<=1: 113 ax.yaxis.set_minor_formatter(pylab.matplotlib.ticker.ScalarFormatter()) 114 115 # set xticks in time format, python2.5 is not new enough for 116 # flexibility, recoding part of AutoDateFormatter to get it 117 if calendar_time: 118 dateLocator = pylab.matplotlib.dates.AutoDateLocator() 119 dateLocator.set_axis(ax.xaxis) 120 dateLocator.refresh() 121 scale = float( dateLocator._get_unit() ) 122 if ( scale == 365.0 ): 123 dateFormatter = pylab.matplotlib.dates.DateFormatter("$%Y$") 124 elif ( scale == 30.0 ): 125 dateFormatter = pylab.matplotlib.dates.DateFormatter("$%y$/$%m$") 126 elif ( (scale == 1.0) or (scale == 7.0) ): 127 dateFormatter = pylab.matplotlib.dates.DateFormatter("$%m$/$%d$") 128 elif ( scale == (1.0/24.0) ): 129 dateFormatter = pylab.matplotlib.dates.DateFormatter("$%d$-$%H$") 130 elif ( scale == (1.0/(24*60)) ): 131 dateFormatter = pylab.matplotlib.dates.DateFormatter("$%H$:$%M$") 132 elif ( scale == (1.0/(24*3600)) ): 133 dateFormatter = pylab.matplotlib.dates.DateFormatter("$%H$:$%M$") 134 135 ax.xaxis.set_major_locator(dateLocator) 136 ax.xaxis.set_major_formatter(dateFormatter) 137 138 else: 139 # set xticks for 4 hours rather than 5 140 xticks = ax.get_xticks() 141 if len(xticks)>1 and xticks[1]-xticks[0]==5: 142 ax.xaxis.set_major_locator(\ 143 pylab.matplotlib.ticker.MultipleLocator(base=2)) 144 145 146 # set tick linewidth 147 for line in ax.get_xticklines() + ax.get_yticklines(): 148 line.set_markersize(10) 149 line.set_markeredgewidth(1)

150

151 -def display_name(columnName):

152 153 """ 154 Format the string columnName (e.g. xml table column) into latex format for 155 an axis label. 156 157 Examples: 158 159 >>> display_name('snr') 160 'SNR' 161 162 >>> display_name('bank_chisq_dof') 163 'Bank $\\chi^2$ DOF' 164 165 Arguments: 166 167 columnName : string 168 string to format 169 """ 170 171 acro = ['snr', 'ra','dof', 'id', 'ms', 'far'] 172 greek = ['alpha', 'beta', 'gamma', 'delta', 'epsilon', 'zeta', 'eta',\ 173 'theta', 'iota', 'kappa', 'lamda', 'mu', 'nu', 'xi', 'omicron',\ 174 'pi', 'rho', 'sigma', 'tau', 'upsilon', 'phi', 'chi', 'psi', 'omega'] 175 unit = ['ns'] 176 sub = ['flow', 'fhigh', 'hrss', 'mtotal', 'mchirp'] 177 178 words = [] 179 for w in re.split('\s', columnName): 180 if w.isupper(): words.append(w) 181 else: words.extend(re.split('_', w)) 182 183 # parse words 184 for i,w in enumerate(words): 185 # get acronym in lower case 186 if w in acro: 187 words[i] = w.upper() 188 # get numerical unit 189 elif w in unit: 190 words[i] = '$(%s)$' % w 191 # get character with subscript text 192 elif w in sub: 193 words[i] = '%s$_{\mbox{\\small %s}}$' % (w[0], w[1:]) 194 # get greek word 195 elif w in greek: 196 words[i] = '$\%s$' % w 197 # get starting with greek word 198 elif re.match('(%s)' % '|'.join(greek), w): 199 if w[-1].isdigit(): 200 words[i] = '$\%s_{%s}$''' % tuple(re.findall(r"[a-zA-Z]+|\d+",w)) 201 elif w.endswith('sq'): 202 words[i] = '$\%s^2$' % w.rstrip('sq') 203 # get everything else 204 else: 205 if w.isupper(): 206 words[i] = w 207 else: 208 words[i] = w.title() 209 # escape underscore 210 words[i] = re.sub('(?<!\\\\)_', '\_', words[i]) 211 212 return ' '.join(words)

213

214 -def gps2datenum(gpstime):

215 216 """ 217 Convert GPS time into floating in standard python time format 218 (days since Jan 01 0000), don't seem to use properly leap seconds 219 """ 220 221 # set time of 0 GPS in datenum units 222 zeroGPS = 722820.0 223 ## correct for leap seconds assuming that all time stamps are within 224 # a range not including a leap 225 # select the first time stamp 226 if isinstance(gpstime,float) or isinstance(gpstime,int): 227 repTime = gpstime 228 else: 229 if len(gpstime)>0: 230 repTime = gpstime[0] 231 else: 232 return gpstime 233 234 zeroGPS = zeroGPS + float(date.XLALLeapSeconds(LIGOTimeGPS(0)) -\ 235 date.XLALLeapSeconds(LIGOTimeGPS(repTime)))/86400 236 # convert to datenum (days) 237 datenum = gpstime/86400 + zeroGPS 238 239 return datenum

240

241 -def calendar_time_unit(duration):

242 243 if duration > 63072000: 244 return "year" 245 elif duration > 5184000: 246 return "year/month" 247 elif duration > 604800: 248 return "month/day" 249 elif duration > 7200: 250 return "day-hour" 251 elif duration > 600: 252 return "hour:minute" 253 else: 254 return "hour:minute:second"

255

256 -def time_unit(duration):

257 258 """ 259 Work out renormalisation for the time axis, makes the label more 260 appropriate. Returns unit (in seconds) and string descriptor 261 262 Example: 263 264 >>> time_unit(100) 265 (1, 'seconds') 266 267 >>> time_unit(604800) 268 (86400, 'days') 269 270 Arguments: 271 272 duration : float 273 plot duration to normalise 274 """ 275 276 # set plot time unit whether it's used or not 277 if (duration) < 1000: 278 unit = 1 279 elif (duration) < 20000: 280 unit = 60 281 elif (duration) >= 20000 and (duration) < 604800: 282 unit = 3600 283 else: 284 unit = 86400 285 timestring = {1:'seconds', 60:'minutes',3600:'hours',86400:'days'} 286 287 return unit, timestring[unit]

288

289 -def getTrigAttribute(trig, col):

290 291 """ 292 Returns the value of trig.col or trig.get_col() for the given string col, 293 and the object trig. If col='time' is given, trig.get_peak() is returned for 294 *Burst* objects, trig.get_end() for *Inspiral* objects and trig.get_start() 295 for *Ringdown* objects. Raises KeyError if cannot execute. 296 297 Arguments: 298 299 trig : [ lsctables.SnglBurst | lscatbles.SnglInspiral | 300 lsctables.SnglRingdown ] 301 xml table entry from which to extract parameter 302 col : string 303 glue.ligolw.table column name to extract 304 """ 305 306 # return time 307 if col=='time': 308 if re.search('Burst', str(trig)): 309 return trig.get_peak() 310 elif re.search('Inspiral', str(trig)): 311 return trig.get_end() 312 elif re.search('Ringdown', str(trig)): 313 return trig.get_start() 314 else: 315 return trig.time 316 317 # return simple column 318 if col in trig.__slots__: 319 return getattr(trig, col) 320 321 # return get_XXX() parameter 322 try: 323 return eval('trig.get_%s()', col) 324 325 # if we get here, panic 326 except: 327 raise KeyError, "Column '%s' not found in %s." % (col, type(trig))

328

329 -def add_colorbar(self, mappable, cax=None, clim=None, log=False, base=10,\ 330 label=None, **kwargs):

331 332 # set axes for colorbar 333 #if not cax: 334 # div = axes_grid.make_axes_locatable(self.ax) 335 # cax = div.new_horizontal(size="2%", pad=0.05, visible=True) 336 #cax = div.append_axes("right", size="2%", pad=0.05) 337 338 cmin, cmax = clim 339 340 # construct colorbar tick formatter, using logic not supported before python 341 # 2.5 342 colorticksize = None 343 if log and numpy.power(base,cmax)-numpy.power(base,cmin) > 4: 344 formatter = pylab.matplotlib.ticker.FuncFormatter(lambda x,pos:\ 345 numpy.power(base,x)>=1 and\ 346 "$%d$" % round(numpy.power(base, x)) or\ 347 numpy.power(base,x)>=0.01 and\ 348 "$%.2f$" % numpy.power(base, x) or\ 349 numpy.power(base,x)<0.01 and "$%f$") 350 elif log and numpy.power(base,cmax)-numpy.power(base,cmin) > 0.4: 351 formatter = pylab.matplotlib.ticker.FuncFormatter(lambda x,pos: "$%.2f$"\ 352 % numpy.power(base, x)) 353 elif log: 354 colorticksize = 'small' 355 formatter = pylab.matplotlib.ticker.FuncFormatter(lambda x,pos: '$%.2e$'\ 356 % numpy.power(base, x)) 357 elif not log and cmax-cmin > 4: 358 formatter = pylab.matplotlib.ticker.FuncFormatter(lambda x,pos:\ 359 x>=1 and "$%d$" % round(x) or\ 360 x>=0.01 and "$%.2f$" % x or\ 361 x==0.0 and "$0$" or\ 362 x and "$%f$") 363 else: 364 formatter = None 365 366 if clim: 367 colorticks = numpy.linspace(cmin, cmax, 4) 368 else: 369 colorticks = None 370 371 self.colorbar = self.ax.figure.colorbar(mappable, cax=cax, format=formatter,\ 372 ticks=colorticks, pad=0.005,\ 373 fraction=0.02, aspect=40, **kwargs) 374 if colorticksize is not None: 375 for l in self.colorbar.ax.yaxis.get_ticklabels(): 376 l.set_size(colorticksize) 377 # colorbar = pylab.colorbar(mappable, cax=cax, format=formatter,\ 378 #ticks=colorticks, **kwargs) 379 self.colorbar.set_label(label) 380 self.colorbar.draw_all()

381

382 #try: 383 # self.colorbar = colorbar 384 #except: 385 # return colorbar 386 387 -def add_hidden_colorbar(self):

388 389 # set axes for colorbar 390 div = axes_grid.make_axes_locatable(self.ax) 391 cax = div.new_horizontal(size="2%", pad=0.05, visible=False)

392

393 # ============================================================================= 394 # Abstract classes for plots 395 # ============================================================================= 396 397 # ============================================================================= 398 # Class for segment plot 399 400 -class PlotSegmentsPlot(plotutils.BasicPlot):

401 402 """ 403 Horizontal bar segment plot. Based originally on PlotSegmentsPlot class in 404 pylal/bin/plotsegments.py 405 """ 406 407 color_code = {'H1':'r', 'H2':'b', 'L1':'g', 'V1':'m', 'G1':'k'} 408

409 - def __init__(self, xlabel="", ylabel="", title="", subtitle="", t0=0, unit=1,\ 410 calendar_time=False):

411 """ 412 Create a fresh plot. Provide t0 to provide a reference time to use as 413 zero. 414 """ 415 plotutils.BasicPlot.__init__(self, xlabel, ylabel, title) 416 self.ax.set_title(title, x=0.5, y=1.025) 417 self.ax.text(0.5, 1.035, subtitle, horizontalalignment='center', 418 transform=self.ax.transAxes, verticalalignment='top') 419 self.segdict = segments.segmentlistdict() 420 self.keys = [] 421 if calendar_time: 422 self._time_transform =\ 423 lambda seg: segments.segment(gps2datenum(float(seg[0])),\ 424 gps2datenum(float(seg[1]))) 425 else: 426 self._time_transform =\ 427 lambda seg: segments.segment(float(seg[0]-t0)/unit,\ 428 float(seg[1]-t0)/unit)

429 430 431

432 - def add_content(self, segdict, keys=None, t0=0, unit=1):

433 if not keys: 434 keys = sorted(segdict.keys()) 435 for key in keys: 436 self.segdict[key] = segdict[key] 437 self.keys.extend(keys)

438

439 - def highlight_segment(self, seg, **plot_args):

440 """ 441 Highlight a particular segment with dashed lines. 442 """ 443 a,b = self._time_transform(seg) 444 plot_args.setdefault('linestyle', '--') 445 plot_args.setdefault('color','r') 446 self.ax.axvline(a, **plot_args) 447 self.ax.axvline(b, **plot_args)

448 449 @plotutils.method_callable_once

450 - def finalize(self, labels_inset=False, hidden_colorbar=False):

451 452 for row,key in enumerate(self.keys): 453 if self.color_code.has_key(key): 454 col = self.color_code[key] 455 else: 456 col = 'b' 457 for seg in self.segdict[key]: 458 a,b = self._time_transform(seg) 459 self.ax.fill([a, b, b, a, a],\ 460 [row-0.4, row-0.4, row+0.4, row+0.4, row-0.4], 'b') 461 if labels_inset: 462 self.ax.text(0.01,(row+1)/(len(self.keys)+1), re.sub('\\+_+','\_',key),\ 463 horizontalalignment='left', verticalalignment='center',\ 464 transform=self.ax.transAxes, backgroundcolor='white',\ 465 bbox=dict(facecolor='white', alpha=0.5, edgecolor='none')) 466 467 ticks = pylab.arange(len(self.keys)) 468 self.ax.set_yticks(ticks) 469 if labels_inset: 470 self.ax.set_yticklabels(ticks, color='white') 471 else: 472 self.ax.set_yticklabels([re.sub(r'\\+_+', '\_', k)\ 473 for k in self.keys], size='small') 474 self.ax.set_ylim(-1, len(self.keys)) 475 476 # add hidden colorbar 477 if hidden_colorbar: 478 self.add_hidden_colorbar()

479 480 PlotSegmentsPlot.add_hidden_colorbar = add_hidden_colorbar

481 482 # ============================================================================= 483 # Class for standard scatter plot 484 485 -class ScatterPlot(plotutils.BasicPlot):

486 487 """ 488 A simple scatter plot, taking x- and y-axis data. 489 """ 490

491 - def __init__(self, xlabel="", ylabel="", title="", subtitle=""):

492 plotutils.BasicPlot.__init__(self, xlabel, ylabel, title) 493 self.ax.set_title(title, x=0.5, y=1.025) 494 self.ax.text(0.5, 1.035, subtitle, horizontalalignment='center', 495 transform=self.ax.transAxes, verticalalignment='top') 496 self.x_data_sets = [] 497 self.y_data_sets = [] 498 self.kwarg_sets = []

499

500 - def add_content(self, x_data, y_data, **kwargs):

501 self.x_data_sets.append(x_data) 502 self.y_data_sets.append(y_data) 503 self.kwarg_sets.append(kwargs)

504 505 @plotutils.method_callable_once

506 - def finalize(self, loc='best', logx=False, logy=False, hidden_colorbar=False):

507 # make plot 508 for x_vals, y_vals, plot_kwargs, c in \ 509 itertools.izip(self.x_data_sets, self.y_data_sets, self.kwarg_sets,\ 510 plotutils.default_colors()): 511 plot_kwargs.setdefault("marker", "o") 512 plot_kwargs.setdefault("s", 20) 513 plot_kwargs.setdefault("linewidth", 1) 514 515 self.ax.scatter(x_vals, y_vals, c=c, **plot_kwargs) 516 517 # add hidden colorbar to make plots the right aspect 518 if hidden_colorbar: 519 self.add_hidden_colorbar() 520 521 # add legend if there are any non-trivial labels 522 self.add_legend_if_labels_exist(loc=loc) 523 524 leg = self.ax.get_legend() 525 # set transparent legend 526 if leg: 527 legfr = leg.get_frame() 528 legfr.set_alpha(0.5) 529 530 # set axes 531 if logx: 532 self.ax.set_xscale('log') 533 if logy: 534 self.ax.set_yscale('log')

535 536 ScatterPlot.add_hidden_colorbar = add_hidden_colorbar

537 538 # ============================================================================= 539 # Class for scatter plot with colorbar 540 541 -class ColorbarScatterPlot(plotutils.BasicPlot):

542 543 """ 544 A scatter plot of x- versus y-data, coloured by z-data. A single colorbar 545 is used, from the union of the ranges of z-data for each content set, 546 unless strict limits are passed as clim=(cmin, cmax) to finalize(). 547 """ 548

549 - def __init__(self, xlabel="", ylabel="", zlabel="", title="", subtitle=""):

550 plotutils.BasicPlot.__init__(self, xlabel, ylabel, title) 551 self.ax.set_title(title, x=0.5, y=1.025) 552 self.ax.text(0.5, 1.035, subtitle, horizontalalignment='center', 553 transform=self.ax.transAxes, verticalalignment='top') 554 555 self.x_data_sets = [] 556 self.y_data_sets = [] 557 self.z_data_sets = [] 558 self.rank_data_sets = [] 559 self.kwarg_sets = [] 560 self.color_label = zlabel

561

562 - def add_content(self, x_data, y_data, z_data, rank_data, **kwargs):

563 self.x_data_sets.append(x_data) 564 self.y_data_sets.append(y_data) 565 self.z_data_sets.append(z_data) 566 self.rank_data_sets.append(rank_data) 567 self.kwarg_sets.append(kwargs)

568 569 @plotutils.method_callable_once

570 - def finalize(self, loc='best', logx=False, logy=False, logz=False, clim=None,\ 571 base=10):

572 573 # set up things we'll need later 574 p = [] 575 576 # set colorbar limits 577 if clim: 578 cmin,cmax = clim 579 else: 580 try: 581 cmin = min([min(z_vals) for z_vals in self.z_data_sets])*0.99 582 cmax = max([max(z_vals) for z_vals in self.z_data_sets])*1.01 583 # catch no triggers 584 except ValueError: 585 cmin = 1 586 cmax = 10 587 588 # reset logs 589 if logz: 590 cmin = numpy.math.log(cmin, base) 591 cmax = numpy.math.log(cmax, base) 592 593 # make plot 594 for x_vals, y_vals, z_vals, rank_vals, plot_kwargs in\ 595 itertools.izip(self.x_data_sets, self.y_data_sets, self.z_data_sets,\ 596 self.rank_data_sets, self.kwarg_sets): 597 598 if logz: z_vals = [numpy.math.log(z, base) for z in z_vals] 599 600 plot_kwargs.setdefault("marker", "o") 601 plot_kwargs.setdefault("vmin", cmin) 602 plot_kwargs.setdefault("vmax", cmax) 603 604 # sort data by z-value 605 zipped = zip(x_vals, y_vals, z_vals, rank_vals) 606 zipped.sort(key=lambda (x,y,z,r): r) 607 x_vals, y_vals, z_vals, rank_vals = map(list, zip(*zipped)) 608 609 p.append(self.ax.scatter(x_vals, y_vals, c=z_vals, **plot_kwargs)) 610 611 if len(p)<1: 612 p.append(self.ax.scatter([1], [1], c=[cmin], vmin=cmin,\ 613 vmax=cmax, visible=False)) 614 615 616 # write colorbar 617 #self.set_colorbar(p[-1], [cmin, cmax], logz, base) 618 self.add_colorbar(p[-1], clim=[cmin, cmax], log=logz,\ 619 label=self.color_label, base=10) 620 621 # set axes 622 if logx: 623 self.ax.set_xscale('log') 624 if logy: 625 self.ax.set_yscale('log') 626 627 self.add_legend_if_labels_exist(loc=loc) 628 629 # set transparent legend 630 leg = self.ax.get_legend() 631 if leg: 632 legfr = leg.get_frame() 633 legfr.set_alpha(0.5)

634

635 - def set_colorbar(self, mappable, clim=None, log=True, base=10):

636 637 cmin, cmax = clim 638 639 # construct colorbar tick formatter, using logic not supported before python 640 # 2.5 641 if log and numpy.power(base,cmax)-numpy.power(base,cmin) > 4: 642 formatter = pylab.matplotlib.ticker.FuncFormatter(lambda x,pos:\ 643 numpy.power(base,x)>=1 and\ 644 "$%d$" % round(numpy.power(base, x)) or\ 645 numpy.power(base,x)>=0.01 and\ 646 "$%.2f$" % numpy.power(base, x) or\ 647 numpy.power(base,x)<0.01 and "$%f$") 648 elif log and numpy.power(base,cmax)-numpy.power(base,cmin) > 0.4: 649 formatter = pylab.matplotlib.ticker.FuncFormatter(lambda x,pos: "$%.2f$"\ 650 % numpy.power(base, x)) 651 elif log: 652 formatter = pylab.matplotlib.ticker.FuncFormatter(lambda x,pos: "$%.4e$"\ 653 % numpy.power(base, x)) 654 else: 655 formatter = pylab.matplotlib.ticker.FuncFormatter(lambda x,pos: "$%.4e$"\ 656 % x) 657 658 if clim: 659 colorticks = numpy.linspace(cmin, cmax, 4) 660 self.colorbar = self.ax.figure.colorbar(mappable, format=formatter,\ 661 ticks=colorticks) 662 else: 663 self.colorbar = self.ax.figure.colorbar(mappable, format=formatter) 664 665 self.colorbar.set_label(self.color_label) 666 self.colorbar.draw_all()

667 668 ColorbarScatterPlot.add_colorbar = add_colorbar

669 670 # ============================================================================= 671 # Extension of ColorbarScatterPlot to plot DetChar-style scatter plot 672 673 -class DetCharScatterPlot(ColorbarScatterPlot):

674 675 """ 676 A 'DetChar' style scatter plot, whereby those triggers under a threshold 677 on the colour column are plotted much smaller than others, allowing line 678 features to be shown easily. 679 """ 680 681 @plotutils.method_callable_once

682 - def finalize(self, loc='best', logx=False, logy=False, logz=True, base=10,\ 683 clim=None, rankthreshold=None):

684 685 p = [] 686 # set colorbar limits 687 if clim: 688 cmin,cmax = clim 689 else: 690 try: 691 cmin = min([min(z_vals) for z_vals in self.z_data_sets])*0.99 692 cmax = max([max(z_vals) for z_vals in self.z_data_sets])*1.01 693 # catch no triggers 694 except ValueError: 695 cmin = 1 696 cmax = 10 697 698 # reset logs 699 if logz: 700 cmin = numpy.math.log(cmin, base) 701 cmax = numpy.math.log(cmax, base) 702 703 # make plot 704 for x_vals, y_vals, z_vals, rank_vals, plot_kwargs in\ 705 itertools.izip(self.x_data_sets, self.y_data_sets, self.z_data_sets,\ 706 self.rank_data_sets, self.kwarg_sets): 707 plot_kwargs.setdefault("vmin", cmin) 708 plot_kwargs.setdefault("vmax", cmax) 709 plot_kwargs.setdefault("s", 15) 710 plot_kwargs.setdefault("marker", "o") 711 plot_kwargs.setdefault("edgecolor", "k") 712 713 if logz: 714 z_vals = [numpy.math.log(z, base) for z in z_vals] 715 716 # sort data by z-value 717 zipped = zip(x_vals, y_vals, z_vals, rank_vals) 718 zipped.sort(key=lambda (x,y,z,r): r) 719 x_vals, y_vals, z_vals, rank_vals = map(list, zip(*zipped)) 720 721 bins = ['low', 'high'] 722 x_bins = {} 723 y_bins = {} 724 z_bins = {} 725 for bin in bins: 726 x_bins[str(bin)] = [] 727 y_bins[str(bin)] = [] 728 z_bins[str(bin)] = [] 729 for i in xrange(len(x_vals)): 730 if rankthreshold and rank_vals[i] < rankthreshold: 731 x_bins[bins[0]].append(float(x_vals[i])) 732 y_bins[bins[0]].append(float(y_vals[i])) 733 z_bins[bins[0]].append(float(z_vals[i])) 734 else: 735 x_bins[bins[1]].append(float(x_vals[i])) 736 y_bins[bins[1]].append(float(y_vals[i])) 737 z_bins[bins[1]].append(float(z_vals[i])) 738 739 # plot bins 740 for i,bin in enumerate(bins): 741 if bin == bins[0]: 742 args = copy.deepcopy(plot_kwargs) 743 args['s']/=4 744 if not (args.has_key('marker') and args['marker']=='x'): 745 args['edgecolor'] = 'none' 746 if len(x_bins[bin])>=1: 747 p.append(self.ax.scatter(x_bins[bin], y_bins[bin], c=z_bins[bin],\ 748 **args)) 749 else: 750 if len(x_bins[bin])>=1: 751 p.append(self.ax.scatter(x_bins[bin], y_bins[bin], c=z_bins[bin],\ 752 **plot_kwargs)) 753 args = plot_kwargs 754 755 if len(p)<1: 756 p.append(self.ax.scatter([1], [1], c=[cmin], vmin=cmin,\ 757 vmax=cmax, visible=False)) 758 759 # write colorbar 760 self.add_colorbar(p[-1], clim=[cmin, cmax], log=logz,\ 761 label=self.color_label, base=10) 762 #self.set_colorbar(p[-1], [cmin, cmax], logz, base) 763 764 # set axes 765 if logx: 766 self.ax.set_xscale('log') 767 if logy: 768 self.ax.set_yscale('log') 769 770 self.add_legend_if_labels_exist(loc=loc)

771

772 # ============================================================================= 773 # Class for line histogram 774 775 -class LineHistogram(ColorbarScatterPlot, plotutils.BasicPlot):

776 777 """ 778 A simple line histogram plot. The values of each histogram bin are plotted 779 using pylab.plot(), with points centred on the x values and height equal 780 to the y values. 781 782 Cumulative, and rate options can be passeed to the finalize() method to 783 format each trace individually. 784 785 """ 786

787 - def __init__(self, xlabel="", ylabel="", title="", subtitle="",\ 788 colorlabel=""):

789 plotutils.BasicPlot.__init__(self, xlabel, ylabel, title) 790 self.ax.set_title(title, x=0.5, y=1.025) 791 self.ax.text(0.5, 1.035, subtitle, horizontalalignment='center', 792 transform=self.ax.transAxes, verticalalignment='top') 793 self.data_sets = [] 794 self.livetimes = [] 795 self.kwarg_sets = [] 796 self.color_label = colorlabel 797 self.color_values = []

798

799 - def add_content(self, data, livetime=1, cval=None, **kwargs):

800 self.data_sets.append(data) 801 self.livetimes.append(livetime) 802 self.kwarg_sets.append(kwargs) 803 self.color_values.append(cval)

804 805 @plotutils.method_callable_once

806 - def finalize(self, loc='best', num_bins=100, cumulative=False, rate=False,\ 807 logx=False, logy=False, fill=False, base=10, xlim=None,\ 808 colorbar=None, clim=None, hidden_colorbar=False):

809 810 # determine binning 811 if xlim: 812 min_stat, max_stat = xlim 813 else: 814 min_stat, max_stat = plotutils.determine_common_bin_limits(self.data_sets) 815 816 if min_stat!=max_stat: 817 max_stat *= 1.001 818 if logx: 819 bins = numpy.logspace(numpy.math.log(min_stat, base),\ 820 numpy.math.log(max_stat, base),\ 821 num_bins+1, endpoint=True) 822 else: 823 bins = numpy.linspace(min_stat, max_stat, num_bins + 1, endpoint=True) 824 else: 825 bins = [min_stat] 826 827 if logy: 828 ymin = 1E-100 829 else: 830 ymin = 0 831 832 # get colors 833 if colorbar: 834 if isinstance(colorbar, str): 835 ctype = colorbar 836 cmap = pylab.matplotlib.colors.LinearSegmentedColormap('clrs',\ 837 pylab.matplotlib.cm.jet._segmentdata) 838 else: 839 ctype, cmap = colorbar 840 assert re.match('(lin|log)', ctype, re.I),\ 841 "colorbar must have type 'linear', or 'log'" 842 try: 843 colors = pylab.matplotlib.colors.makeMappingArray(100000,\ 844 cmap) 845 except IndexError: 846 xind = numpy.linspace(0, 1, 100000) 847 colors = numpy.clip(numpy.array(cmap(xind), dtype=numpy.float), 0, 1) 848 849 if clim: 850 cmin,cmax = clim 851 else: 852 try: 853 cmin = min(self.color_values) 854 cmax = max(self.color_values) 855 except ValueError: 856 cmin = 1 857 cmax = 10 858 859 if re.search('log', ctype, re.I): 860 cmin = numpy.math.log(cmin, base) 861 cmax = numpy.math.log(cmax, base) 862 self.color_values = [numpy.math.log(x, base) for x in self.color_values] 863 color_idx = lambda x: int((x-cmin)/(cmax-cmin)*(len(colors)-1)) 864 self.color_values = [colors[color_idx(x)] for x in self.color_values] 865 866 # add hidden colorbar to make plots the right aspect 867 elif hidden_colorbar: 868 self.add_hidden_colorbar() 869 870 p = [] 871 for i, (data_set, livetime, plot_kwargs, col) in\ 872 enumerate(itertools.izip(self.data_sets, self.livetimes,\ 873 self.kwarg_sets, self.color_values)): 874 875 # 876 # make histogram 877 # 878 879 # get version 880 v = [int(i) for i in numpy.version.version.split('.')] 881 if v[1] < 1: 882 y, x = numpy.histogram(data_set, bins=bins) 883 elif v[1] < 3: 884 y, x = numpy.histogram(data_set, bins=bins, new=True) 885 x = x[:-1] 886 else: 887 y, x = numpy.histogram(data_set, bins=bins) 888 x = x[:-1] 889 890 # get cumulative sum 891 if cumulative: 892 y = y[::-1].cumsum()[::-1] 893 894 # convert to rate 895 if rate: 896 if livetime > 0: 897 y = y/livetime 898 ymin /= livetime 899 else: 900 y = numpy.empty(y.shape) 901 y.fill(numpy.nan) 902 ymin = numpy.nan 903 904 # reset zeros on logscale, tried with numpy, unreliable 905 if logy: 906 y = y.astype(float) 907 numpy.putmask(y, y==0, ymin) 908 909 # set color 910 if colorbar: 911 plot_kwargs['color'] = col 912 913 # plot 914 if fill: 915 plot_kwargs.setdefault('linewidth', 1) 916 plot_kwargs.setdefault('alpha', 0.8) 917 p.append(self.ax.plot(x, y, **plot_kwargs)) 918 self.ax.fill_between(x, 1E-100, y, **plot_kwargs) 919 else: 920 p.append(self.ax.plot(x, y, **plot_kwargs)) 921 922 # set axes 923 if logx: 924 self.ax.set_xscale('log') 925 if logy: 926 self.ax.set_yscale('log') 927 928 # set colorbar 929 if colorbar: 930 self.add_colorbar(self.ax.scatter([1], [1], c=1, cmap=cmap,\ 931 vmin=cmin, vmax=cmax, visible=False),\ 932 clim=[cmin, cmax], label=self.color_label) 933 elif hidden_colorbar: 934 self.add_hidden_colorbar() 935 936 # add legend if there are any non-trivial labels 937 self.add_legend_if_labels_exist(loc=loc) 938 939 # fix legend 940 leg = self.ax.get_legend() 941 if leg: 942 for l in leg.get_lines(): 943 l.set_linewidth(4)

944 945 LineHistogram.add_hidden_colorbar = add_hidden_colorbar

946 947 # ============================================================================= 948 # Extension of VerticalBarHistogram to include log axes 949 950 -class VerticalBarHistogram(plotutils.VerticalBarHistogram):

951

952 - def __init__(self, xlabel="", ylabel="", title="", subtitle=""):

953 plotutils.BasicPlot.__init__(self, xlabel, ylabel, title) 954 self.ax.set_title(title, x=0.5, y=1.025) 955 self.ax.text(0.5, 1.035, subtitle, horizontalalignment='center', 956 transform=self.ax.transAxes, verticalalignment='top') 957 self.data_sets = [] 958 self.kwarg_sets = []

959 960 @plotutils.method_callable_once

961 - def finalize(self, num_bins=20, normed=False, logx=False, logy=False,\ 962 base=10, hidden_colorbar=False, loc='best'):

963 964 # determine binning 965 min_stat, max_stat = plotutils.determine_common_bin_limits(self.data_sets) 966 if logx and min_stat!=0 and max_stat!=0: 967 bins = numpy.logspace(numpy.math.log(min_stat, base),\ 968 numpy.math.log(max_stat, base),\ 969 num_bins+1, endpoint=True) 970 else: 971 bins = numpy.linspace(min_stat, max_stat, num_bins + 1, endpoint=True) 972 973 # determine bar width; gets silly for more than a few data sets 974 if logx: 975 width = [bins[i+1]-bins[i] for i in xrange(len(bins)-1)] 976 width.append(width[-1]) 977 else: 978 width = (1 - 0.1 * len(self.data_sets)) * (bins[1] - bins[0]) 979 980 width = numpy.asarray(width)/2 981 982 # get version 983 v = map(int, numpy.version.version.split('.')) 984 if v[1] >= 1: width = width[:-1] 985 986 # set base of plot in log scale 987 if logy: 988 ymin = (base**-1)*5 989 else: 990 ymin = 0 991 992 # make plot 993 legends = [] 994 plot_list = [] 995 for i, (data_set, plot_kwargs) in \ 996 enumerate(itertools.izip(self.data_sets, self.kwarg_sets)): 997 # set default values 998 plot_kwargs.setdefault("alpha", 0.6) 999 plot_kwargs.setdefault("align", "center") 1000 plot_kwargs.setdefault("width", width) 1001 if logy: 1002 plot_kwargs.setdefault("bottom", ymin) 1003 else: 1004 plot_kwargs.setdefault("bottom", ymin) 1005 1006 # make histogram 1007 if v[1] < 1: 1008 y, x = numpy.histogram(data_set, bins=bins, normed=normed) 1009 elif v[1] < 3: 1010 y, x = numpy.histogram(data_set, bins=bins, new=True, normed=normed) 1011 x = x[:-1] 1012 else: 1013 y, x = numpy.histogram(data_set, bins=bins, normed=normed) 1014 x = x[:-1] 1015 1016 if logy: 1017 y = y-ymin 1018 1019 # stagger bins for pure aesthetics 1020 #x += 0.1 * i * max_stat / num_bins 1021 1022 # plot 1023 self.ax.bar(x, y, **plot_kwargs) 1024 1025 # set axes 1026 if logx: 1027 self.ax.set_xscale('log') 1028 if logy: 1029 self.ax.set_yscale('log') 1030 1031 # set logy minimum 1032 self.ax.set_ybound(lower=ymin) 1033 1034 # add legend if there are any non-trivial labels 1035 self.add_legend_if_labels_exist(loc=loc) 1036 # set transparent legend 1037 leg = self.ax.get_legend() 1038 if leg: leg.get_frame().set_alpha(0.5) 1039 1040 1041 # add hidden colorbar 1042 if hidden_colorbar: 1043 self.add_hidden_colorbar()

1044 1045 VerticalBarHistogram.add_hidden_colorbar = add_hidden_colorbar

1046 1047 # ============================================================================= 1048 # Class for time series plot 1049 1050 -class DataPlot(plotutils.BasicPlot):

1051 1052 """ 1053 Time-series data plot. Just a nice looking line plot. 1054 """ 1055

1056 - def __init__(self, xlabel="", ylabel="", title="", subtitle=""):

1057 plotutils.BasicPlot.__init__(self, xlabel, ylabel, title) 1058 self.ax.set_title(title, x=0.5, y=1.025) 1059 self.ax.text(0.5, 1.035, subtitle, horizontalalignment='center', 1060 transform=self.ax.transAxes, verticalalignment='top') 1061 self.x_data_sets = [] 1062 self.y_data_sets = [] 1063 self.kwarg_sets = []

1064

1065 - def add_content(self, x_data, y_data, **kwargs):

1066 self.x_data_sets.append(numpy.ma.asarray(x_data)) 1067 self.y_data_sets.append(numpy.ma.asarray(y_data)) 1068 self.kwarg_sets.append(kwargs)

1069 1070 @plotutils.method_callable_once

1071 - def finalize(self, loc='best', logx=False, logy=False, hidden_colorbar=False):

1072 # make plot 1073 plots = [] 1074 markersizes = [] 1075 markerscale = 4 1076 1077 for x_vals, y_vals, plot_kwargs, c in \ 1078 itertools.izip(self.x_data_sets, self.y_data_sets,\ 1079 self.kwarg_sets, plotutils.default_colors()): 1080 1081 # magnify the markers on the legend 1082 plot_kwargs.setdefault('markersize', 5) 1083 markersizes.append(plot_kwargs['markersize']) 1084 plot_kwargs['markersize'] = min(20, plot_kwargs['markersize']*markerscale) 1085 # plot 1086 if logx: 1087 x_vals = x_vals.astype(float) 1088 numpy.putmask(x_vals, x_vals==0, 1E-100) 1089 if logy: 1090 y_vals = y_vals.astype(float) 1091 numpy.putmask(y_vals, y_vals==0, 1E-100) 1092 1093 plots.append(self.ax.plot(x_vals, y_vals, **plot_kwargs)) 1094 1095 # add legend if there are any non-trivial labels 1096 self.add_legend_if_labels_exist(loc=loc) 1097 leg = self.ax.get_legend() 1098 # magnify the lines on the legend 1099 if leg: 1100 try: 1101 for l in leg.get_lines(): 1102 if l.get_linewidth(): 1103 l.set_linewidth(4) 1104 except AttributeError: 1105 pass 1106 1107 # reset markersizes on plot 1108 for i,plot in enumerate(plots): 1109 l = plot[0] 1110 l.set_markersize(markersizes[i]) 1111 1112 # set transparent legend 1113 if leg: 1114 legfr = leg.get_frame() 1115 legfr.set_alpha(0.5) 1116 1117 # add hidden colorbar 1118 if hidden_colorbar: self.add_hidden_colorbar() 1119 1120 # set axes 1121 if logx: self.ax.set_xscale('log') 1122 if logy: self.ax.set_yscale('log')

1123 1124 DataPlot.add_hidden_colorbar = add_hidden_colorbar

1125 1126 # ============================================================================= 1127 # Class for color map plot 1128 1129 -class ColorMap(ColorbarScatterPlot):

1130

1131 - def __init__(self, xlabel="", ylabel="", zlabel="", title="", subtitle=""):

1132 plotutils.BasicPlot.__init__(self, xlabel, ylabel, title) 1133 self.color_label = zlabel 1134 self.ax.set_title(title, x=0.5, y=1.025) 1135 self.ax.text(0.5, 1.035, subtitle, horizontalalignment='center', 1136 transform=self.ax.transAxes, verticalalignment='top') 1137 self.data_sets = [] 1138 self.extent_sets = [] 1139 self.kwarg_sets = []

1140

1141 - def add_content(self, data, extent, **kwargs):

1142 self.data_sets.append(numpy.asarray(data)) 1143 self.extent_sets.append(extent) 1144 self.kwarg_sets.append(kwargs)

1145 1146 @plotutils.method_callable_once

1147 - def finalize(self, loc='best', logx=False, logy=False, logz=False, clim=None,\ 1148 origin=pylab.rcParams['image.origin'], base=10):

1149 1150 # set colorbar limits 1151 if clim: 1152 cmin,cmax = clim 1153 else: 1154 cmin = min([z_vals.min() for z_vals in self.data_sets])*0.99 1155 cmax = max([z_vals.max() for z_vals in self.data_sets])*1.01 1156 1157 # reset logs 1158 if logz: 1159 cmin = numpy.math.log(cmin, base) 1160 cmax = numpy.math.log(cmax, base) 1161 1162 p = [] 1163 1164 for i, (data_set, extent, plot_kwargs) in \ 1165 enumerate(itertools.izip(self.data_sets, self.extent_sets,\ 1166 self.kwarg_sets)): 1167 1168 plot_kwargs.setdefault("vmin", cmin) 1169 plot_kwargs.setdefault("vmax", cmax) 1170 plot_kwargs.setdefault("norm", None) 1171 plot_kwargs.setdefault("interpolation", "kaiser") 1172 1173 if logz: 1174 if base==10: 1175 data_set = numpy.log10(data_set) 1176 elif base==numpy.e: 1177 data_set = numpy.log(data_set) 1178 elif base==2: 1179 data_set = numpy.log2(data_set) 1180 else: 1181 raise AttributeError("Can only use base = 2, e, or 10.") 1182 1183 p.append(self.ax.imshow(data_set, extent=extent, origin=origin,\ 1184 **plot_kwargs)) 1185 1186 if len(p)==0: 1187 p.append(self.ax.imshow([[1]], vmin=cmin,\ 1188 vmax=cmax, visible=False)) 1189 1190 # write colorbar 1191 self.add_colorbar(p[-1], clim=[cmin, cmax], log=logz,\ 1192 label=self.color_label, base=10) 1193 1194 # set axes 1195 if logx: self.ax.set_xscale('log') 1196 if logy: self.ax.set_yscale('log')

1197

1198 # ============================================================================= 1199 # Projection for sky position plotting 1200 1201 -class SkyPositionsPlot(ScatterPlot):

1202 1203 """ 1204 Plot of sky positions plotted onto mpl_toolkits basemap projection. 1205 """ 1206 1207 @plotutils.method_callable_once

1208 - def finalize(self, loc='lower right', projection='moll', centre=None,\ 1209 range=[(0, 0), (1, 1)]):

1210 1211 """ 1212 Finalize and draw plot. Can only be called once. 1213 1214 Keyword arguments: 1215 1216 toc : [ str | int ] 1217 compatible value for legend loc, default: 'lower right' 1218 projection : str 1219 valid projection name for mpl_toolkits.basemap.Basemap 1220 centre : list 1221 (ra, dec) in degrees for point to centre on plot, only works for 1222 some projections 1223 range : list 1224 [(xmin, ymin), (xmax, ymax)] pair of tuples for lower left and upper 1225 right corners of plot area. Full areai (entire sphere) is 1226 [(0,0), (1,1)]. Narrower range is zoomed in, wider range zoomed out. 1227 """ 1228 1229 # set up projection 1230 projection = projection.lower() 1231 if not centre: 1232 centre = (0,0) 1233 centre = [(-(centre[0]-180)+180) % 360, centre[1]] 1234 m1 = Basemap(projection=projection, lon_0=centre[0], lat_0=centre[1],\ 1235 resolution=None, ax=self.ax) 1236 1237 # set up range 1238 width = m1.urcrnrx 1239 height = m1.urcrnry 1240 range = [list(range[0]), list(range[1])] 1241 maprange = [[-width/2, -height/2], [width/2, height/2]] 1242 if range[0][0] != None: maprange[0][0] = width/2*(range[0][0]-1) 1243 if range[0][1] != None: maprange[0][1] = height/2*(range[0][1]-1) 1244 if range[1][0] != None: maprange[1][0] = width/2*(range[1][0]) 1245 if range[1][1] != None: maprange[1][1] = height/2*(range[1][1]) 1246 1247 # set projection 1248 m = Basemap(projection=projection, lon_0=centre[0], lat_0=centre[1],\ 1249 llcrnrx=maprange[0][0], llcrnry=maprange[0][1],\ 1250 urcrnrx=maprange[1][0], urcrnry=maprange[1][1],\ 1251 resolution=None, ax=self.ax) 1252 1253 # turn on 'fulldisk' property when using full disk 1254 if maprange == [[-width/2, -height/2], [width/2, height/2]]: 1255 m._fulldisk = True 1256 1257 xrange = (m.llcrnrx, m.urcrnrx) 1258 yrange = (m.llcrnry, m.urcrnry) 1259 1260 # make plot 1261 for x_vals, y_vals, plot_kwargs, c in\ 1262 itertools.izip(self.x_data_sets, self.y_data_sets, self.kwarg_sets,\ 1263 plotutils.default_colors()): 1264 1265 plot_kwargs.setdefault("marker", "o") 1266 plot_kwargs.setdefault("edgecolor", "k") 1267 plot_kwargs.setdefault("facecolor", c) 1268 1269 # project data 1270 convert = lambda x: (x>=180 and x-360) or (x<180 and x) 1271 x_vals = [-convert(x) for x in x_vals] 1272 if projection in ['moll', 'hammer', 'orth']: 1273 convert = lambda x: (x>=0 and x) or (x<0 and x+360) 1274 x_vals, y_vals = m(x_vals, y_vals) 1275 m.scatter(x_vals, y_vals, **plot_kwargs) 1276 1277 # finish projection 1278 m.drawmapboundary() 1279 1280 # set labels 1281 if projection in ['ortho']: 1282 plabels = [0, 0, 0, 0] 1283 mlabels = [0, 0, 0, 0] 1284 else: 1285 plabels = [1, 0, 0, 0] 1286 mlabels = [0, 0, 0, 0] 1287 1288 # draw parallels 1289 parallels = numpy.arange(-90., 120., 30.) 1290 m.drawparallels(parallels, labels=plabels,\ 1291 labelstyle='+/-', latmax=90) 1292 1293 # draw meridians 1294 if projection in ['moll', 'hammer', 'ortho']: 1295 meridians = numpy.arange(0., 360., 45.) 1296 else: 1297 meridians = numpy.arange(-180, 181, 45) 1298 m.drawmeridians(meridians, labels=mlabels,\ 1299 latmax=90, labelstyle='+/-') 1300 1301 # label parallels for certain projections 1302 if projection in ['ortho']: 1303 for lon,lat in zip([0.]*len(parallels), parallels): 1304 #if lonrange[0]<=lon<=lonrange[1] and latrange[0]<=lat<=latrange[1]: 1305 x, y = m(lon, lat) 1306 lon, lat = m1(x, y, inverse=True) 1307 if x<=10**20 and y<=10**20\ 1308 and xrange[0]<x<xrange[1] and yrange[0]<=y<=yrange[1]: 1309 m.ax.text(x, y, r"$%0.0f^\circ$" % lat) 1310 1311 # label meridians for all projections 1312 for lon,lat in zip(meridians, [0.]*len(meridians)): 1313 tlon = (-(lon-180)+180) % 360 1314 x, y = m(lon, lat) 1315 lon, lat = m1(x, y, inverse=True) 1316 1317 if x<=10**20 and y<=10**20\ 1318 and xrange[0]<x<xrange[1] and yrange[0]<=y<=yrange[1]: 1319 m.ax.text(x, y, r"$%0.0f^\circ$" % tlon) 1320 1321 # set legend 1322 self.add_legend_if_labels_exist(loc=loc, scatterpoints=1)

1323

1324 -class ColorbarSkyPositionsPlot(ColorbarScatterPlot):

1325 1326 @plotutils.method_callable_once

1327 - def finalize(self, loc='lower right', projection='moll', centre=None,\ 1328 range=[(0, 0), (1, 1)], logz=False, clim=None,\ 1329 base=10):

1330 1331 """ 1332 Finalize and draw plot. Can only be called once. 1333 1334 Keyword arguments: 1335 1336 toc : [ str | int ] 1337 compatible value for legend loc, default: 'lower right' 1338 projection : str 1339 valid projection name for mpl_toolkits.basemap.Basemap 1340 centre : list 1341 (ra, dec) in degrees for point to centre on plot, only works for 1342 some projections 1343 range : list 1344 [(xmin, ymin), (xmax, ymax)] pair of tuples for lower left and upper 1345 right corners of plot area. Full area (entire sphere) is 1346 [(0,0), (1,1)]. Narrower range is zoomed in, wider range zoomed out. 1347 logz : [ True | False ] 1348 plot z-axis in log scale 1349 clim : 1350 """ 1351 1352 p = [] 1353 1354 # set up projection 1355 projection = projection.lower() 1356 if not centre: 1357 centre = (0,0) 1358 m1 = Basemap(projection=projection, lon_0=centre[0], lat_0=centre[1],\ 1359 resolution=None, ax=self.ax) 1360 1361 # set up range 1362 width = m1.urcrnrx 1363 height = m1.urcrnry 1364 range = [list(range[0]), list(range[1])] 1365 maprange = [[-width/2, -height/2], [width/2, height/2]] 1366 if range[0][0] != None: maprange[0][0] = width/2*(range[0][0]-1) 1367 if range[0][1] != None: maprange[0][1] = height/2*(range[0][1]-1) 1368 if range[1][0] != None: maprange[1][0] = width/2*(range[1][0]) 1369 if range[1][1] != None: maprange[1][1] = height/2*(range[1][1]) 1370 1371 # set projection 1372 m = Basemap(projection=projection, lon_0=centre[0], lat_0=centre[1],\ 1373 llcrnrx=maprange[0][0], llcrnry=maprange[0][1],\ 1374 urcrnrx=maprange[1][0], urcrnry=maprange[1][1],\ 1375 resolution=None, ax=self.ax) 1376 1377 # turn on 'fulldisk' property when using full disk 1378 if maprange == [[-width/2, -height/2], [width/2, height/2]]: 1379 m._fulldisk = True 1380 1381 xrange = (m.llcrnrx, m.urcrnrx) 1382 yrange = (m.llcrnry, m.urcrnry) 1383 1384 # set colorbar limits 1385 if clim: 1386 cmin,cmax = clim 1387 else: 1388 try: 1389 cmin = min([min(z_vals) for z_vals in self.z_data_sets])*0.99 1390 cmax = max([max(z_vals) for z_vals in self.z_data_sets])*1.01 1391 # catch no triggers 1392 except ValueError: 1393 cmin = 1 1394 cmax = 10 1395 1396 # reset logs 1397 if logz: 1398 cmin = numpy.math.log(cmin, base) 1399 cmax = numpy.math.log(cmax, base) 1400 1401 # make plot 1402 for x_vals, y_vals, z_vals, plot_kwargs in\ 1403 itertools.izip(self.x_data_sets, self.y_data_sets, self.z_data_sets,\ 1404 self.kwarg_sets): 1405 1406 if logz: z_vals = [numpy.math.log(z, base) for z in z_vals] 1407 1408 plot_kwargs.setdefault("marker", "o") 1409 plot_kwargs.setdefault("edgecolor", "k") 1410 plot_kwargs.setdefault("vmin", cmin) 1411 plot_kwargs.setdefault("vmax", cmax) 1412 1413 # sort data by z-value 1414 zipped = zip(x_vals, y_vals, z_vals) 1415 zipped.sort(key=lambda (x,y,z): z) 1416 x_vals, y_vals, z_vals = map(list, zip(*zipped)) 1417 1418 # project data 1419 if projection in ['moll', 'hammer', 'orth']: 1420 convert = lambda x: (x>=0 and x) or (x<0 and x+360) 1421 else: 1422 convert = lambda x: (x>=180 and x-360) or (x<180 and x) 1423 x_vals = [convert(x) for x in x_vals] 1424 x_vals, y_vals = m(x_vals, y_vals) 1425 1426 # pull out z value if given 1427 if plot_kwargs.has_key('facecolor'): 1428 self.ax.scatter(x_vals, y_vals, **plot_kwargs) 1429 else: 1430 p.append(self.ax.scatter(x_vals, y_vals, c=z_vals, **plot_kwargs)) 1431 1432 # finish projection 1433 m.drawmapboundary() 1434 1435 # set labels 1436 if projection in ['ortho']: 1437 plabels = [0, 0, 0, 0] 1438 mlabels = [0, 0, 0, 0] 1439 else: 1440 plabels = [1, 0, 0, 0] 1441 mlabels = [0, 0, 0, 1] 1442 1443 # draw parallels 1444 parallels = numpy.arange(-90., 120., 30.) 1445 m.drawparallels(parallels, labels=plabels,\ 1446 labelstyle='+/-', latmax=90) 1447 1448 # draw meridians 1449 if projection in ['moll', 'hammer', 'ortho']: 1450 meridians = numpy.arange(0., 360., 45.) 1451 else: 1452 meridians = numpy.arange(-180, 181, 45) 1453 m.drawmeridians(meridians, labels=mlabels,\ 1454 latmax=90, labelstyle='+/-') 1455 1456 # label parallels for certain projections 1457 if projection in ['ortho']: 1458 for lon,lat in zip([0.]*len(parallels), parallels): 1459 x, y = m(lon, lat) 1460 if x<=10**20 and y<=10**20\ 1461 and xrange[0]<x<xrange[1] and yrange[0]<=y<=yrange[1]: 1462 self.ax.text(x, y, r"$%0.0f^\circ$" % lat) 1463 # label meridians for certain projections 1464 if projection in ['moll', 'hammer', 'ortho']: 1465 for lon,lat in zip(meridians, [0.]*len(meridians)): 1466 tlon = lon 1467 while tlon < 0: tlon += 360 1468 x, y = m(lon, lat) 1469 if x<=10**20 and y<=10**20\ 1470 and xrange[0]<x<xrange[1] and yrange[0]<=y<=yrange[1]: 1471 self.ax.text(x, y, r"$%0.0f^\circ$" % tlon) 1472 1473 # write colorbar 1474 self.set_colorbar(p[-1], [cmin, cmax], logz, base) 1475 1476 # add legend if there are any non-trivial labels 1477 self.add_legend_if_labels_exist(loc=loc, scatterpoints=1)

1478

1479 # ============================================================================= 1480 # Wrappers to translate ligolw table into plot 1481 # ============================================================================= 1482 1483 # ============================================================================= 1484 # Plot time series of data set 1485 1486 -def plot_data_series(data, outfile, x_format='time', zero=None, \ 1487 zeroindicator=False, **kwargs):

1488 1489 """ 1490 Plot the time series / spectrum of a given set (or given sets) of data. 1491 1492 Arguments: 1493 1494 data : list 1495 list of (ChannelName,x_data,y_data) tuples with channel name (or data 1496 source) and time/freq, amplitude arrays for each channel. Channels are 1497 plotted in the order given. 1498 outfile : str 1499 output plot path 1500 1501 Keyword Arguments: 1502 1503 x_format : [ 'time' | 'frequency' ] 1504 type of data for x_axis, allows formatting of axes 1505 zero : [ float | int | LIGOTimeGPS ] 1506 time around which to centre time series plot 1507 zeroindicator : [ False | True ] 1508 indicate zero time with veritcal dashed line, default: False 1509 1510 Unnamed keyword arguments: 1511 1512 logx : [ True | False ] 1513 boolean option to display x-axis in log scale. 1514 logy : [ True | False ] 1515 boolean option to display y-axis in log scale. 1516 xlim : tuple 1517 (xmin, xmax) limits for x-axis 1518 ylim : tuple 1519 (ymin, ymax) limits for y-axis 1520 xlabel : string 1521 label for x-axis 1522 ylabel : string 1523 label for y-axis 1524 title : string 1525 title for plot 1526 subtitle : string 1527 subtitle for plot 1528 1529 All other given arguments will be passed to matplotlib.axes.Axes.plot. 1530 """ 1531 1532 # format times 1533 if not kwargs.has_key('xlim') or not kwargs['xlim']: 1534 start = min([min(d[1]) for d in data]) 1535 end = max([max(d[1]) for d in data]) 1536 kwargs['xlim'] = [start,end] 1537 1538 if not zero: 1539 zero = kwargs['xlim'][0] 1540 1541 # set plot time unit whether it's used or not 1542 if x_format=='time': 1543 unit, timestr = time_unit(kwargs['xlim'][1]-kwargs['xlim'][0]) 1544 1545 # set labels 1546 if x_format=='time': 1547 zero = LIGOTimeGPS('%.3f' % zero) 1548 if zero.nanoseconds==0: 1549 tlabel = datetime(*date.XLALGPSToUTC(LIGOTimeGPS(zero))[:6])\ 1550 .strftime("%B %d %Y, %H:%M:%S %ZUTC") 1551 else: 1552 tlabel = datetime(*date.XLALGPSToUTC(LIGOTimeGPS(zero.seconds))[:6])\ 1553 .strftime("%B %d %Y, %H:%M:%S %ZUTC") 1554 tlabel = tlabel.replace(' UTC', '.%.3s UTC' % zero.nanoseconds) 1555 xlabel = kwargs.pop('xlabel',\ 1556 'Time (%s) since %s (%s)' % (timestr, tlabel, zero)) 1557 title = kwargs.pop('title', 'Time series') 1558 else: 1559 xlabel = kwargs.pop('xlabel', 'Frequency (Hz)') 1560 title = kwargs.pop('title', 'Frequency spectrum') 1561 1562 ylabel = kwargs.pop('ylabel', 'Amplitude') 1563 subtitle = kwargs.pop('subtitle', '') 1564 1565 # customise plot appearance 1566 set_rcParams() 1567 1568 # get limits 1569 xlim = kwargs.pop('xlim', None) 1570 ylim = kwargs.pop('ylim', None) 1571 calendar_time = kwargs.pop('calendar_time', False) 1572 1573 # get axis scales 1574 logx = kwargs.pop('logx', False) 1575 logy = kwargs.pop('logy', False) 1576 1577 # get legend loc 1578 loc = kwargs.pop('loc', 'best') 1579 1580 # get colorbar options 1581 hidden_colorbar = kwargs.pop('hidden_colorbar', False) 1582 1583 # get savefig option 1584 bbox_inches = kwargs.pop('bbox_inches', 'tight') 1585 1586 # generate plot object 1587 plot = DataPlot(xlabel, ylabel, title, subtitle) 1588 1589 # set plot params 1590 style = kwargs.pop('style', '-') 1591 if style in ['-', '--', '-.', ':']: 1592 kwargs.setdefault('linestyle', style) 1593 kwargs.setdefault('linewidth', 2) 1594 kwargs.pop('marker', None) 1595 else: 1596 kwargs.setdefault('marker', style) 1597 kwargs.setdefault('markersize', 5) 1598 kwargs.setdefault('linewidth', 0) 1599 kwargs.pop('linestyle', ' ') 1600 1601 # get uniq data sets that aren't errors 1602 allchannels = [] 1603 channels = [] 1604 for i,(c,_,_) in enumerate(data): 1605 allchannels.append(c) 1606 if not re.search('(min|max)\Z', str(c)): 1607 channels.append((i,c)) 1608 1609 # add data 1610 for i,c in channels: 1611 x_data,y_data = data[i][1:] 1612 if x_format=='time': 1613 if calendar_time: 1614 x_data = gps2datenum(numpy.array(map(float, x_data))) 1615 else: 1616 x_data = (x_data.astype(float)-float(zero))/unit 1617 lab = str(c) 1618 if lab != '_': lab = lab.replace('_', '\_') 1619 plot.add_content(x_data, y_data, label=lab,**kwargs) 1620 1621 # finalize plot 1622 plot.finalize(logx=logx, logy=logy, loc=loc, hidden_colorbar=hidden_colorbar) 1623 1624 # plot errors 1625 l = None 1626 for (i,channel),c in itertools.izip(channels, plotutils.default_colors()): 1627 try: 1628 minidx = allchannels.index('%s_min' % str(channel)) 1629 maxidx = allchannels.index('%s_max' % str(channel)) 1630 except ValueError: 1631 continue 1632 y = [] 1633 for idx in [minidx,maxidx]: 1634 x_data,y_data = data[idx][1:] 1635 y.append(y_data) 1636 if x_format=='time': 1637 x_data = (numpy.array(map(float, x_data))-float(zero))/unit 1638 l = l!=None and l or float(kwargs.pop('linewidth', 1))/4 1639 plot.ax.plot(x_data, y_data, color=c, linewidth=l, **kwargs) 1640 plot.ax.fill_between(x_data, y[1], y[0], color=c, alpha=0.25) 1641 1642 # set axes 1643 plot.ax.autoscale_view(tight=True, scalex=True, scaley=True) 1644 if ylim: 1645 plot.ax.set_ylim(tuple(ylim)) 1646 1647 # FIXME add zero indicator 1648 if x_format=='time': 1649 axis_lims = plot.ax.get_ylim() 1650 if zeroindicator: 1651 plot.ax.plot([0, 0], [axis_lims[0], axis_lims[1]], 'r--', linewidth=2) 1652 plot.ax.set_ylim([ axis_lims[0], axis_lims[1] ]) 1653 1654 # set x axis 1655 if xlim and calendar_time: 1656 plot.ax.set_xlim([gps2datenum(float(xlim[0])),\ 1657 gps2datenum(float(xlim[1]))]) 1658 elif xlim: 1659 plot.ax.set_xlim([ float(xlim[0]-zero)/unit, float(xlim[1]-zero)/unit ]) 1660 else: 1661 # set global axis limits 1662 if xlim: 1663 plot.ax.set_xlim(xlim) 1664 if ylim: 1665 plot.ax.set_ylim(ylim) 1666 1667 1668 set_ticks(plot.ax, calendar_time=calendar_time) 1669 1670 plot.savefig(outfile, bbox_inches=bbox_inches,\ 1671 bbox_extra_artists=plot.ax.texts) 1672 pylab.close(plot.fig)

1673

1674 # ============================================================================= 1675 # Plot a histogram of any column 1676 1677 -def plot_trigger_hist(triggers, outfile, column='snr', num_bins=1000,\ 1678 color_column=None, color_bins=10,\ 1679 seglist=None, flag='unknown', start=None, end=None,\ 1680 livetime=None, etg=None, **kwargs):

1681 1682 """ 1683 Wrapper for dqPlotUtils.LineHistogram to plot a histogram of the value in 1684 any column of the ligolw table triggers. If a glue.segments.segmentlist 1685 seglist is given, the histogram is presented before and after removal of 1686 triggers falling inside any segment in the list. 1687 1688 Arguments: 1689 1690 triggers : glue.ligolw.table.Table 1691 ligolw table containing triggers 1692 outfile : string 1693 string path for output plot 1694 1695 Keyword arguments: 1696 1697 column : string 1698 valid column of triggers table to plot as histrogram 1699 num_bins : int 1700 number of histogram bins to use 1701 seglist : glue.segments.segmentlist 1702 list of segments with which to veto triggers 1703 flag : string 1704 display name of segmentlist, normally the name of the DQ flag 1705 start : [ float | int | LIGOTimeGPS] 1706 GPS start time (exclude triggers and segments before this time) 1707 end : [ float | int | LIGOTimeGPS] 1708 GPS end time (exclude triggers and segments after this time) 1709 livetime : [ float | int | LIGOTimeGPS ] 1710 span of time from which triggers and segments are valid, used to 1711 display histogram counts in terms of rate (Hz) for easy comparisons 1712 etg : string 1713 display name of trigger generator, defaults based on triggers tableName 1714 1715 Unnamed keyword arguments: 1716 1717 cumulative : [ True | False ] 1718 plot cumulative histogram 1719 rate : [ True | False ] 1720 plot rate histogram (normalises with given or calculated livetime) 1721 fill : [ True | False ] 1722 fill below the histogram curves, default colors: 1723 red (vetoed), green (not vetoed). 1724 logx : [ True | False ] 1725 boolean option to display x-axis in log scale. 1726 logy : [ True | False ] 1727 boolean option to display y-axis in log scale. 1728 xlim : tuple 1729 (xmin, xmax) limits for x-axis 1730 ylim : tuple 1731 (ymin, ymax) limits for y-axis 1732 xlabel : string 1733 label for x-axis 1734 ylabel : string 1735 label for y-axis 1736 title : string 1737 title for plot 1738 subtitle : string 1739 subtitle for plot 1740 greyscale : [ True | False ] 1741 use (non-greyscale) colour scheme suitable for greyscale plots 1742 1743 All other given arguments will be passed to matplotlib.axes.Axes.plot and 1744 matplotlib.axes.Axes.fill_between. 1745 """ 1746 1747 get_time = def_get_time(triggers.tableName) 1748 1749 # calculate livetime 1750 if not start or not end: 1751 times = [get_time(t) for t in triggers] 1752 if not start: 1753 start = min(times) 1754 if not end: 1755 end = max(times) 1756 if not livetime: 1757 livetime = end-start 1758 livetime = float(livetime) 1759 1760 # format seglist 1761 if seglist==None: 1762 seglist = segments.segmentlist() 1763 else: 1764 seglist = segments.segmentlist(seglist) 1765 1766 # format colours 1767 logz = kwargs.pop('logz', False) 1768 clim = kwargs.pop('clim', None) 1769 1770 if color_column: 1771 colData = get_column(triggers, color_column).astype(float) 1772 if not clim: 1773 if len(colData)>0: 1774 clim = [colData.min(), colData.max()] 1775 else: 1776 clim = [1,10] 1777 if isinstance(color_bins, int): 1778 num_color_bins = color_bins 1779 if logz: 1780 color_bins = numpy.logspace(numpy.math.log10(clim[0]),\ 1781 numpy.math.log10(clim[1]),\ 1782 num=num_color_bins) 1783 else: 1784 color_bins = numpy.linspace(0, clim[1],\ 1785 num=num_color_bins) 1786 else: 1787 num_color_bins = len(color_bins) 1788 else: 1789 num_color_bins = 1 1790 1791 # get data 1792 tdata = get_column(triggers, 'time') 1793 preData = [] 1794 postData = [] 1795 for i in range(num_color_bins): 1796 if color_column: 1797 preData.append(get_column(triggers, column).astype(float)\ 1798 [colData>=color_bins[i]]) 1799 else: 1800 preData.append(get_column(triggers, column).astype(float)) 1801 postData.append([p for j,p in enumerate(preData[i])\ 1802 if tdata[j] not in seglist]) 1803 1804 # get veto livetime 1805 vetoLivetime = livetime-float(abs(seglist)) 1806 1807 # set some random plot parameters 1808 greyscale = kwargs.pop('greyscale', False) 1809 if seglist: 1810 color = ['r','g'] 1811 label = ['Before vetoes', 'After vetoes'] 1812 else: 1813 color = ['b'] 1814 label = ['_'] 1815 if greyscale: 1816 color = ['k','k'] 1817 linestyle = ['-','--'] 1818 else: 1819 linestyle = ['-','-'] 1820 1821 # fix names for latex 1822 if flag: 1823 flag = flag.replace('_','\_') 1824 else: 1825 flag = 'Unknown' 1826 1827 # customise plot appearance 1828 set_rcParams() 1829 1830 # get limits 1831 xlim = kwargs.pop('xlim', None) 1832 ylim = kwargs.pop('ylim', None) 1833 1834 # get axis scales 1835 logx = kwargs.pop('logx', False) 1836 logy = kwargs.pop('logy', False) 1837 1838 # get colorbar options 1839 hidden_colorbar = kwargs.pop('hidden_colorbar', False) 1840 1841 # get savefig option 1842 bbox_inches = kwargs.pop('bbox_inches', 'tight') 1843 1844 # get fill 1845 fill = kwargs.pop('fill', False) 1846 1847 # get extras 1848 cumulative = kwargs.pop('cumulative', False) 1849 rate = kwargs.pop('rate', False) 1850 1851 # set labels 1852 xlabel = kwargs.pop('xlabel', display_name(column)) 1853 if rate and cumulative: 1854 ylabel = kwargs.pop('ylabel', 'Cumulative rate (Hz)') 1855 elif rate: 1856 ylabel = kwargs.pop('ylabel', 'Rate (Hz)') 1857 elif not rate and cumulative: 1858 ylabel = kwargs.pop('ylabel', 'Cumulative number') 1859 elif not rate and not cumulative: 1860 ylabel = kwargs.pop('ylabel', 'Number') 1861 1862 # get ETG 1863 if not etg: 1864 if re.search('burst', triggers.tableName.lower()): 1865 etg = 'Burst' 1866 elif re.search('inspiral', triggers.tableName.lower()): 1867 etg = 'Inspiral' 1868 elif re.search('ringdown', triggers.tableName.lower()): 1869 etg = 'Ringdown' 1870 else: 1871 etg = 'Unknown' 1872 1873 title = '%s triggers' % (etg.replace('_','\_')) 1874 if seglist: 1875 title += ' and %s segments' % (flag) 1876 title = kwargs.pop('title', title) 1877 if start and end: 1878 subtitle = '%d-%d' % (start, math.ceil(end)) 1879 else: 1880 subtitle = "" 1881 subtitle = kwargs.pop('subtitle', subtitle) 1882 zlabel = kwargs.pop('zlabel',\ 1883 color_column and display_name(color_column) or "") 1884 1885 # generate plot object 1886 plot = LineHistogram(xlabel, ylabel, title, subtitle, zlabel) 1887 1888 # add each data set 1889 if color_column: 1890 for i in range(len(preData)): 1891 plot.add_content(preData[i], livetime=livetime, cval=color_bins[i],\ 1892 linestyle=linestyle[0], label=label[0], **kwargs) 1893 else: 1894 plot.add_content(preData[0], livetime=livetime, color=color[0],\ 1895 linestyle=linestyle[0], label=label[0], **kwargs) 1896 if seglist: 1897 plot.add_content(postData[0], livetime=vetoLivetime, color=color[1],\ 1898 linestyle=linestyle[1], label=label[1], **kwargs) 1899 1900 1901 1902 # finalize plot with histograms 1903 if not num_bins: num_bins=100 1904 if color_column: 1905 colorbar = logz and 'log' or 'linear' 1906 else: 1907 colorbar = None 1908 plot.finalize(num_bins=num_bins, logx=logx, logy=logy, cumulative=cumulative,\ 1909 rate=rate, fill=fill, colorbar=colorbar,\ 1910 hidden_colorbar=hidden_colorbar) 1911 plot.ax.autoscale_view(tight=True, scalex=True, scaley=True) 1912 1913 # set lower y axis limit 1914 if rate: 1915 ymin = 1/livetime 1916 elif logy: 1917 ymin = 0.5 1918 else: 1919 ymin = plot.ax.get_ylim()[0] 1920 plot.ax.set_ybound(lower=ymin) 1921 1922 # set global axis limits 1923 if xlim: 1924 plot.ax.set_xlim(xlim) 1925 if ylim: 1926 plot.ax.set_ylim(ylim) 1927 1928 # set global ticks 1929 set_ticks(plot.ax) 1930 1931 # save figure 1932 plot.savefig(outfile, bbox_inches=bbox_inches,\ 1933 bbox_extra_artists=plot.ax.texts)

1934

1935 # ============================================================================= 1936 # Plot one column against another column coloured by any third column 1937 1938 -def plot_triggers(triggers, outfile, reftriggers=None, xcolumn='time',\ 1939 ycolumn='snr', zcolumn=None, rankcolumn=None, etg=None, start=None, end=None,\ 1940 zero=None, seglist=None, flag=None, **kwargs):

1941 1942 """ 1943 Plots ycolumn against xcolumn for columns in given 1944 Sngl{Burst,Inspiral}Table object triggers, coloured by the zcolumn 1945 highlighting those entries falling inside one of the entries in the 1946 glue.segments.segmentlist object segments, if given. 1947 1948 'time' given as a column name is a special case, since s and ns times are 1949 stored separately in the SnglTable structures. In this case the 1950 trigger.get_xxxx() function is called. 1951 1952 Arguments: 1953 1954 triggers : glue.ligolw.table.Table 1955 ligolw table containing triggers 1956 outfile : string 1957 string path for output plot 1958 1959 Keyword arguments: 1960 1961 xcolumn : string 1962 valid column of triggers table to plot on x-axis 1963 ycolumn : string 1964 valid column of triggers table to plot on y-axis 1965 zcolumn : string 1966 valid column of triggers table to use for colorbar (optional). 1967 rankcolumn : string 1968 valid column of triggers table to use for ranking events (optional). 1969 etg : string 1970 display name of trigger generator, defaults based on triggers tableName 1971 start : [ float | int | LIGOTimeGPS ] 1972 GPS start time of plot 1973 end : [ float | int | LIGOTimeGPS ] 1974 GPS end time of plot 1975 zero : [ float | int | LIGOTimeGPS ] 1976 time around which to centre plot 1977 seglist : glue.segments.segmentlist 1978 list of segments with which to veto triggers 1979 flag : string 1980 display name of segmentlist, normally the name of the DQ flag 1981 1982 Unnamed keyword arguments: 1983 1984 detchar : [ True | False ] 1985 use 'DetChar' style for scatter plot with colorbar, triggers below given 1986 dcthreshold are small with no edges, whilst other triggers are normal 1987 dcthreshold : float 1988 threshold below which scatter points are small with no edges when using 1989 DetChar plotting style 1990 logx : [ True | False ] 1991 boolean option to display x-axis in log scale. 1992 logy : [ True | False ] 1993 boolean option to display y-axis in log scale. 1994 logz : [ True | False ] 1995 boolean option to display z-axis in log scale. 1996 xlim : tuple 1997 (xmin, xmax) limits for x-axis. Triggers outside range are removed. 1998 ylim : tuple 1999 (ymin, ymax) limits for y-axis. Triggers outside range are removed. 2000 zlim : tuple 2001 (zmin, zmax) limits for z-axis. Triggers outside range are removed. 2002 clim : tuple 2003 (cmin, cmax) limits for color scale. Triggers outside range are moved 2004 onto boundary. 2005 xlabel : string 2006 label for x-axis 2007 ylabel : string 2008 label for y-axis 2009 zlabel : string 2010 label for z-axis 2011 title : string 2012 title for plot 2013 subtitle : string 2014 subtitle for plot 2015 greyscale : [ True | False ] 2016 use (non-greyscale) colour scheme suitable for greyscale plots 2017 2018 All other given arguments will be passed to matplotlib.axes.Axes.scatter. 2019 """ 2020 2021 from pylal import plotutils 2022 2023 # test multiple tables 2024 if not len(triggers)==0 and \ 2025 (isinstance(triggers[0], tuple) or isinstance(triggers[0], list)): 2026 assert not zcolumn,\ 2027 "Can only plot single table when using colorbar plot" 2028 tables = [t[1] for t in triggers] 2029 tablelabel = [t[0] for t in triggers] 2030 for i,t in enumerate(tablelabel): 2031 if t!='_': 2032 tablelabel[i] = t.replace('_','\_') 2033 else: 2034 tables = [triggers] 2035 tablelabel = '_' 2036 2037 # get time column 2038 get_time = [] 2039 for t in tables: 2040 get_time.append(def_get_time(t.tableName)) 2041 2042 # set start and end time if needed 2043 if not start or not end: 2044 times = [get_time[i](t) for i in xrange(len(tables)) for t in tables[i]] 2045 if not start and len(times)>=1: 2046 start = int(math.floor(min(times))) 2047 elif not start: 2048 start = 0 2049 if not end and len(times)>=1: 2050 end = int(math.ceil(max(times))) 2051 elif not end: 2052 end = 1 2053 2054 # set zero time 2055 if not zero: 2056 zero = start 2057 2058 # get time params 2059 unit,timestr = time_unit(end-start) 2060 2061 # set up segmentlist 2062 if seglist: 2063 segs = segments.segmentlist(seglist) 2064 if not seglist: 2065 segs = segments.segmentlist() 2066 2067 # get axis limits 2068 xlim = kwargs.pop('xlim', None) 2069 ylim = kwargs.pop('ylim', None) 2070 zlim = kwargs.pop('zlim', None) 2071 clim = kwargs.pop('clim', zlim) 2072 calendar_time = kwargs.pop('calendar_time', None) 2073 2074 # set up columns 2075 columns = list(map(str.lower, [xcolumn, ycolumn])) 2076 if zcolumn: 2077 columns.append(zcolumn.lower()) 2078 if rankcolumn: 2079 columns.append(rankcolumn.lower()) 2080 else: columns.append(zcolumn.lower()) 2081 2082 # set up limits 2083 limits = [xlim, ylim, zlim, None] 2084 for i,col in enumerate(columns): 2085 if re.search('time\Z', col) and not limits[i]: 2086 limits[i] = [start,end] 2087 2088 # get veto info 2089 if seglist: 2090 tdata = get_column(triggers, 'time') 2091 2092 # get all data 2093 vetoData = [] 2094 nvetoData = [] 2095 for j,tab in enumerate(tables): 2096 vetoData.append({}) 2097 nvetoData.append({}) 2098 # get veto info 2099 if seglist: 2100 tdata = get_column(tab, 'time') 2101 # get data 2102 for i,col in enumerate(columns): 2103 nvetoData[j][col] = get_column(tab, col).astype(float) 2104 # apply limits and vetoes 2105 condition = True 2106 for i,col in enumerate(columns): 2107 if limits[i]: 2108 condition = condition & (limits[i][0] <= nvetoData[j][col])\ 2109 & (nvetoData[j][col] <= limits[i][1]) 2110 for col in nvetoData[j].keys(): 2111 nvetoData[j][col] = nvetoData[j][col][condition] 2112 if seglist: 2113 vetoData[j][col] = numpy.asarray([d for i,d in\ 2114 enumerate(nvetoData[j][col]) if\ 2115 tdata[i] in seglist]) 2116 else: 2117 vetoData[j][col] = numpy.array([]) 2118 2119 data = {} 2120 2121 2122 2123 # normalize zcolumn by time-averaged value 2124 whitenedFlag = kwargs.pop('whitened', False) 2125 if zcolumn and whitenedFlag: 2126 # get ref data if provided 2127 refData = {} 2128 if reftriggers: 2129 for i,col in enumerate(columns): 2130 refData[col] = get_column(reftriggers, col).astype(float) 2131 else: 2132 for i,col in enumerate(columns): 2133 refData[col] = nvetoData[0][col] 2134 2135 uniqYvalues = numpy.unique1d(refData[ycolumn]) 2136 # building look back table by hand, is included in unique1d for numpy >= v1.3 2137 for yVal in uniqYvalues: 2138 backTable = numpy.where(yVal == nvetoData[0][ycolumn]) 2139 zMedian = numpy.median(refData[zcolumn][yVal == refData[ycolumn]]) 2140 for iTrig in backTable[0]: 2141 nvetoData[0][zcolumn][iTrig] /= zMedian 2142 2143 # filter zcolumn by provided poles/zeros filter as a function of ycolumn 2144 flatenedFlag = kwargs.pop('filter', False) 2145 if zcolumn and flatenedFlag: 2146 # get filter params 2147 polesList = kwargs.pop('poles', None) 2148 zerosList = kwargs.pop('zeros', None) 2149 amplitude = kwargs.pop('amplitude', 1) 2150 nvetoData[0][zcolumn] *= amplitude 2151 for filtPole in polesList: 2152 nvetoData[0][zcolumn] /= abs(nvetoData[0][ycolumn] - filtPole) 2153 for filtZero in zerosList: 2154 nvetoData[0][zcolumn] *= abs(nvetoData[0][ycolumn] - filtZero) 2155 nvetoData[0][zcolumn].astype(float) 2156 2157 # flaten zcolumn by 1/sqrt of sum given rational fraction mononomes as a 2158 # function of ycolumn 2159 flatenedFlag = kwargs.pop('flaten', False) 2160 if zcolumn and flatenedFlag: 2161 # get filter params 2162 expList = kwargs.pop('exponents', None) 2163 constList = kwargs.pop('constants', None) 2164 filter = numpy.zeros(len(nvetoData[0][zcolumn])) 2165 for iTerm, exponent in enumerate(expList): 2166 filter += pow(constList[iTerm]*numpy.power(nvetoData[0][ycolumn],expList[iTerm]),2) 2167 filter = numpy.sqrt(filter) 2168 nvetoData[0][zcolumn] /= filter 2169 2170 # median/min/max of ycolumn binned by exact xcolumn values 2171 minmaxmedianFlag = kwargs.pop('minmaxmedian', False) 2172 if minmaxmedianFlag: 2173 uniqXvalues = numpy.unique1d(nvetoData[j][xcolumn]) 2174 # building look back table by hand, is included in unique1d for numpy >= v1.3 2175 for xVal in uniqXvalues: 2176 backTable = numpy.where(xVal == nvetoData[j][xcolumn]) 2177 if len(backTable[0]) > 3: 2178 nvetoData[j][ycolumn][backTable[0][0]] =\ 2179 numpy.median(nvetoData[j][ycolumn][xVal == nvetoData[j][xcolumn]]) 2180 nvetoData[j][ycolumn][backTable[0][1]] =\ 2181 numpy.min(nvetoData[j][ycolumn][xVal == nvetoData[j][xcolumn]]) 2182 nvetoData[j][ycolumn][0][backTable[0][2]] =\ 2183 numpy.max(nvetoData[j][ycolumn][xVal == nvetoData[j][xcolumn]]) 2184 for iTrig in backTable[0][3:]: 2185 nvetoData[j][ycolumn][iTrig] = numpy.nan 2186 2187 # down-sample (xaxis) the triggers by plotting only median z-value over averaging window 2188 medianDuration = float(kwargs.pop('medianduration', 0)) 2189 stdDuration = float(kwargs.pop('stdduration', 0)) 2190 if medianDuration or stdDuration: 2191 uniqYvalues = numpy.unique1d(nvetoData[0][ycolumn]) 2192 if medianDuration: 2193 avDuration = medianDuration 2194 elif stdDuration: 2195 avDuration = stdDuration 2196 tedges = numpy.arange(float(start), float(end), avDuration) 2197 tedges = numpy.append(tedges, float(end)) 2198 # array for repacking triggers into time-frequency bins 2199 repackTrig = {} 2200 for yVal in uniqYvalues: 2201 repackTrig[yVal] = {} 2202 for iBin in range(len(tedges)-1): 2203 repackTrig[yVal][iBin] = [] 2204 # building new set of triggers 2205 newTrigs = {} 2206 newTrigs[xcolumn] = [] 2207 newTrigs[ycolumn] = [] 2208 newTrigs[zcolumn] = [] 2209 for iTrig in range(len(nvetoData[0][zcolumn])): 2210 trigVal = nvetoData[0][ycolumn][iTrig] 2211 trigBin = math.floor((nvetoData[0][xcolumn][iTrig]-float(start))/avDuration) 2212 repackTrig[trigVal][trigBin].append(nvetoData[0][zcolumn][iTrig]) 2213 for yVal in uniqYvalues: 2214 for iBin in range(len(tedges)-1): 2215 # keep only if at least a few elements, std doesn't make sens otherwise 2216 if medianDuration and len(repackTrig[yVal][iBin]) > 0: 2217 newTrigs[xcolumn].append( (tedges[iBin]+tedges[iBin+1])/2 ) 2218 newTrigs[ycolumn].append(yVal) 2219 newTrigs[zcolumn].append(numpy.median(numpy.array(repackTrig[yVal][iBin]))) 2220 elif stdDuration and len(repackTrig[yVal][iBin]) > 5: 2221 newTrigs[xcolumn].append( (tedges[iBin]+tedges[iBin+1])/2 ) 2222 newTrigs[ycolumn].append(yVal) 2223 newTrigs[zcolumn].append(numpy.std(numpy.array(repackTrig[yVal][iBin]))/ \ 2224 numpy.median(numpy.array(repackTrig[yVal][iBin]))) 2225 if newTrigs[zcolumn][-1] == 0: 2226 newTrigs[zcolumn][-1] = numpy.nan 2227 for i,col in enumerate(columns): 2228 nvetoData[0][col] = numpy.array(newTrigs[col]) 2229 2230 # get limits 2231 for i,col in enumerate(columns): 2232 if not limits[i]: 2233 limits[i] = [0,0] 2234 for j in xrange(len(tables)): 2235 data[col] = numpy.concatenate((nvetoData[j][col], vetoData[j][col])) 2236 if len(data[col])>=1: 2237 limits[i][0] = min(data[col].min()*0.99, limits[i][0]) 2238 limits[i][1] = max(data[col].max()*1.01, limits[i][1]) 2239 2240 # renormalise time and set time axis label unless given 2241 if re.search('time\Z', col): 2242 renormalise = True 2243 if kwargs.has_key('xlabel'): renormalise = False 2244 if renormalise: 2245 for j in xrange(len(tables)): 2246 if calendar_time: 2247 vetoData[j][col] = gps2datenum(vetoData[j][col]) 2248 nvetoData[j][col] = gps2datenum(nvetoData[j][col]) 2249 else: 2250 vetoData[j][col] = (vetoData[j][col]-float(zero))/unit 2251 nvetoData[j][col] = (nvetoData[j][col]-float(zero))/unit 2252 if calendar_time: 2253 limits[i] = [gps2datenum(float(limits[i][0])),\ 2254 gps2datenum(float(limits[i][1]))] 2255 else: 2256 limits[i] = [float(limits[i][0]-zero)/unit,\ 2257 float(limits[i][1]-zero)/unit] 2258 2259 # set labels 2260 label = {} 2261 for i,col in enumerate(['xcolumn', 'ycolumn', 'zcolumn']): 2262 if i >= len(columns): continue 2263 if re.search('time\Z', columns[i]): 2264 zerostr = datetime(*date.XLALGPSToUTC(LIGOTimeGPS(zero))[:6])\ 2265 .strftime("%B %d %Y, %H:%M:%S %ZUTC") 2266 lab = 'Time (%s) since %s (%s)' % (timestr, zerostr, zero) 2267 label[columns[i]] = kwargs.pop('%slabel' % col[0], lab) 2268 else: 2269 label[columns[i]] = kwargs.pop('%slabel' % col[0],\ 2270 display_name(columns[i])) 2271 2272 # find loudest event 2273 loudest = {} 2274 if len(columns)==4 and\ 2275 len(nvetoData[0][columns[0]])+len(vetoData[0][columns[0]])>=1: 2276 # find loudest vetoed event 2277 vetomax = 0 2278 if len(vetoData[0][columns[3]])>=1: 2279 vetomax = vetoData[0][columns[3]].max() 2280 nvetomax = 0 2281 # find loudest unvetoed event 2282 if len(nvetoData[0][columns[3]])>=1: 2283 nvetomax = nvetoData[0][columns[3]].max() 2284 if vetomax == nvetomax == 0: 2285 pass 2286 # depending on which one is loudest, find loudest overall event 2287 elif vetomax > nvetomax: 2288 index = vetoData[0][columns[3]].argmax() 2289 for col in columns: 2290 loudest[col] = vetoData[0][col][index] 2291 else: 2292 index = nvetoData[0][columns[3]].argmax() 2293 for col in columns: 2294 loudest[col] = nvetoData[0][col][index] 2295 2296 # fix flag for use with latex 2297 if flag: 2298 flag = flag.replace('_', '\_') 2299 else: 2300 flag = 'Unknown' 2301 2302 # get ETG 2303 if not etg: 2304 if re.search('burst', tables[0].tableName.lower()): 2305 etg = 'Burst' 2306 elif re.search('inspiral', tables[0].tableName.lower()): 2307 etg = 'Inspiral' 2308 elif re.search('ringdown', tables[0].tableName.lower()): 2309 etg = 'Ringdown' 2310 else: 2311 etg = 'Unknown' 2312 else: 2313 etg = etg.replace('_', '\_') 2314 2315 # customise plot appearance 2316 set_rcParams() 2317 2318 # set title 2319 title = '%s triggers' % etg 2320 if seglist: 2321 title += ' and %s segments' % (flag) 2322 title = kwargs.pop('title', title) 2323 2324 if len(columns)==4 and loudest: 2325 subtitle = "Loudest event by %s:" % display_name(columns[-1]) 2326 for col in columns: 2327 maxcol = loudest[col] 2328 if re.search('time\Z', col) and renormalise: 2329 maxcol = maxcol*unit+zero 2330 loudstr = "%s=%.2f" % (display_name(col), maxcol) 2331 if not re.search(loudstr, subtitle): 2332 subtitle += ' %s' % loudstr 2333 elif start and end: 2334 subtitle = '%s-%s' % (start, end) 2335 else: 2336 subtitle = '' 2337 subtitle = kwargs.pop('subtitle', subtitle) 2338 2339 # get axis scales 2340 logx = kwargs.pop('logx', False) 2341 logy = kwargs.pop('logy', False) 2342 logz = kwargs.pop('logz', False) 2343 2344 # get colorbar options 2345 hidden_colorbar = kwargs.pop('hidden_colorbar', False) 2346 2347 # get savefig option 2348 bbox_inches = kwargs.pop('bbox_inches', 'tight') 2349 2350 # get detchar plot params 2351 detchar = kwargs.pop('detchar', False) 2352 dcthresh = float(kwargs.pop('dcthreshold', 10)) 2353 2354 # get greyscale param 2355 greyscale = kwargs.pop('greyscale', False) 2356 if greyscale and not kwargs.has_key('cmap'): 2357 kwargs['cmap'] = pylab.matplotlib.colors.LinearSegmentedColormap('clrs',\ 2358 pylab.matplotlib.cm.hot._segmentdata) 2359 2360 # initialise standard scatter plot 2361 if len(columns)==2: 2362 plotutils.default_colors = lambda: itertools.cycle(('b', 'r', 'g', 'c', 'm', 'y', 'k')) 2363 plot = ScatterPlot(label[columns[0]], label[columns[1]], title, subtitle) 2364 for j in xrange(len(tables)): 2365 if len(nvetoData[j][columns[0]])>=1: 2366 plot.add_content(nvetoData[j][columns[0]], nvetoData[j][columns[1]],\ 2367 label=tablelabel[j], **kwargs) 2368 # add veto triggers 2369 if len(vetoData[j][columns[0]])>=1: 2370 plot.add_content(vetoData[j][columns[0]], vetoData[j][columns[1]],\ 2371 label=tablelabel[j], marker='x', color='r') 2372 # finalise 2373 plot.finalize(logx=logx, logy=logy, hidden_colorbar=hidden_colorbar) 2374 # initialise scatter plot with colorbar 2375 elif len(columns)==4: 2376 # initialize color bar plot 2377 if detchar: 2378 plot = DetCharScatterPlot(label[columns[0]], label[columns[1]],\ 2379 label[columns[2]], title, subtitle) 2380 else: 2381 plot = ColorbarScatterPlot(label[columns[0]], label[columns[1]],\ 2382 label[columns[2]], title, subtitle) 2383 2384 # add non veto triggers 2385 if len(nvetoData[0][columns[0]])>=1: 2386 plot.add_content(nvetoData[0][columns[0]], nvetoData[0][columns[1]],\ 2387 nvetoData[0][columns[2]], nvetoData[0][columns[3]],\ 2388 **kwargs) 2389 # add veto triggers 2390 if len(vetoData[0][columns[0]])>=1: 2391 plot.add_content(vetoData[0][columns[0]], vetoData[0][columns[1]],\ 2392 vetoData[0][columns[2]], vetoData[0][columns[2]],\ 2393 marker='x', edgecolor='r',\ 2394 **kwargs) 2395 # finalise 2396 if detchar: 2397 plot.finalize(logx=logx, logy=logy, logz=logz, clim=clim,\ 2398 rankthreshold=dcthresh) 2399 else: 2400 plot.finalize(logx=logx, logy=logy, logz=logz, clim=clim) 2401 # add loudest event to plot 2402 if loudest: 2403 plot.ax.plot([loudest[columns[0]]], [loudest[columns[1]]],\ 2404 marker='*', color='gold', markersize=15) 2405 2406 # set axes 2407 plot.ax.autoscale_view(tight=True, scalex=True, scaley=True) 2408 2409 if limits[0]: 2410 plot.ax.set_xlim(limits[0]) 2411 if limits[1]: 2412 plot.ax.set_ylim(limits[1]) 2413 2414 # reset ticks 2415 set_ticks(plot.ax, calendar_time=calendar_time) 2416 2417 # get both major and minor grid lines 2418 plot.savefig(outfile, bbox_inches=bbox_inches,\ 2419 bbox_extra_artists=plot.ax.texts) 2420 pylab.close(plot.fig)

2421

2422 # ============================================================================= 2423 # Plot a histogram of segment duration 2424 2425 -def plot_segment_hist(segs, outfile, keys=None, num_bins=100, coltype=int,\ 2426 **kwargs):

2427 2428 """ 2429 segments. 2430 Plots a histogram of segment duration for the glue.segments.segmentlist 2431 2432 Arguments: 2433 2434 segs : [ glue.segments.segmentlist | glue.segments.segmentlistdict ] 2435 list of segments with which to veto triggers, use dict for multiple 2436 datasets 2437 outfile : string 2438 string path for output plot 2439 2440 Keyword arguments: 2441 2442 flag : string 2443 display name for segments, normally the name of the DQ flag 2444 logx : [ True | False ] 2445 boolean option to display x-axis in log scale. 2446 logy : [ True | False ] 2447 boolean option to display y-axis in log scale. 2448 """ 2449 2450 # customise plot appearance 2451 set_rcParams() 2452 2453 # get limits 2454 xlim = kwargs.pop('xlim', None) 2455 ylim = kwargs.pop('ylim', None) 2456 2457 # get labels 2458 xlabel = kwargs.pop('xlabel', 'Length of segment (seconds)') 2459 ylabel = kwargs.pop('ylabel', 'Number of segments') 2460 title = kwargs.pop('title', 'Segment Duration Histogram') 2461 subtitle = kwargs.pop('subtitle', "") 2462 2463 # get axis scale 2464 logx = kwargs.pop('logx', False) 2465 logy = kwargs.pop('logy', False) 2466 2467 # get savefig option 2468 bbox_inches = kwargs.pop('bbox_inches', 'tight') 2469 2470 # get colorbar option 2471 hidden_colorbar = kwargs.pop('hidden_colorbar', False) 2472 2473 # format mutltiple segments 2474 if isinstance(segs,list): 2475 if keys and isinstance(keys, 'str'): 2476 segs = segments.segmentlistdict({keys:segs}) 2477 keys = [keys] 2478 elif keys: 2479 segs = segments.segmentlistdict({keys[0]:segs}) 2480 else: 2481 segs = segments.segmentlistdict({'_':segs}) 2482 keys = ['_'] 2483 else: 2484 if not keys: 2485 keys = segs.keys() 2486 for i,flag in enumerate(keys): 2487 keys[i] = display_name(flag) 2488 if keys[i]!=flag: 2489 segs[keys[i]] = segs[flag] 2490 del segs[flag] 2491 2492 if kwargs.has_key('color'): 2493 colors = kwargs.pop('color').split(',') 2494 else: 2495 colors = zip(*zip(range(len(keys)), plotutils.default_colors()))[-1] 2496 2497 # generate plot object 2498 plot = VerticalBarHistogram(xlabel, ylabel, title, subtitle) 2499 2500 # add each segmentlist 2501 for flag,c in zip(keys, colors): 2502 plot.add_content([float(abs(seg)) for seg in segs[flag]],\ 2503 label=flag, color=c, **kwargs) 2504 2505 # finalize plot with histograms 2506 plot.finalize(num_bins=num_bins, logx=logx, logy=logy,\ 2507 hidden_colorbar=hidden_colorbar) 2508 2509 # set limits 2510 plot.ax.autoscale_view(tight=True, scalex=True, scaley=True) 2511 if ylim: 2512 ylim = map(float, ylim) 2513 plot.ax.set_ylim(ylim) 2514 if xlim: 2515 xlim = map(float, xlim) 2516 plot.ax.set_xlim(xlim) 2517 2518 # save figure 2519 plot.savefig(outfile, bbox_inches=bbox_inches,\ 2520 bbox_extra_artists=plot.ax.texts) 2521 pylab.close(plot.fig)

2522

2523 # ============================================================================= 2524 # Plot rate versus time in bins 2525 2526 -def plot_trigger_rate(triggers, outfile, average=600, start=None, end=None,\ 2527 zero=None, bincolumn='peak_frequency', bins=[],\ 2528 etg='Unknown', **kwargs):

2529 2530 """ 2531 Plot rate versus time for the given ligolw table triggers, binned by the 2532 given bincolumn using the bins list. 2533 2534 Arguments: 2535 2536 triggers : glue.ligolw.table 2537 LIGOLW table containing a list of triggers 2538 outfile : string 2539 string path for output plot 2540 2541 Keyword arguments: 2542 2543 average : float 2544 Length (seconds) of rate segment 2545 start : [ float | int | LIGOTimeGPS ] 2546 GPS start time 2547 end : [ float | int | LIGOTimeGPS ] 2548 GPS end time 2549 zero : [ float | int | LIGOTimeGPS ] 2550 GPS time to use for 0 on time axis 2551 bincolumn : string 2552 valid column of the trigger table to use for binning 2553 bins : list 2554 list of tuples defining the rate bins 2555 etg : string 2556 display name of trigger generator 2557 logy : [ True | False ] 2558 boolean option to display y-axis in log scale 2559 ylim : tuple 2560 (ymin, ymax) limits for rate axis 2561 """ 2562 2563 tableName = triggers.tableName.lower() 2564 get_time = def_get_time(tableName) 2565 2566 # set start and end times 2567 if not start and not end: 2568 times = [get_time(t) for t in triggers] 2569 if not start: 2570 start = min(times) 2571 if not end: 2572 end = max(times) 2573 2574 if not zero: 2575 zero = start 2576 2577 # set plot time unit whether it's used or not 2578 unit, timestr = time_unit(end-start) 2579 2580 # set ETG 2581 if not etg: 2582 if re.search('burst', tableName): 2583 etg = 'Burst' 2584 elif re.search('inspiral', tableName): 2585 etg = 'Inspiral' 2586 elif re.search('ringdown', tableName): 2587 etg = 'Ringdown' 2588 else: 2589 etg = 'Unknown' 2590 2591 # get limits 2592 calendar_time = kwargs.pop('calendar_time', False) 2593 if calendar_time: 2594 xlim = kwargs.pop('xlim', [gps2datenum(float(start)),\ 2595 gps2datenum(float(end))]) 2596 else: 2597 xlim = kwargs.pop('xlim', [float(start-zero)/unit, float(end-zero)/unit]) 2598 ylim = kwargs.pop('ylim', None) 2599 2600 # get axis scales 2601 logx = kwargs.pop('logx', False) 2602 logy = kwargs.pop('logy', False) 2603 2604 # get averaging time 2605 average = kwargs.pop('average',average) 2606 2607 # get colorbar options 2608 hidden_colorbar = kwargs.pop('hidden_colorbar', False) 2609 2610 # get savefig option 2611 bbox_inches = kwargs.pop('bbox_inches', 'tight') 2612 2613 # format ybins 2614 if not bins: 2615 bins = [[0,float('inf')]] 2616 ybins = [map(float, bin) for bin in bins] 2617 2618 # bin data 2619 tbins = {} 2620 rate = {} 2621 for bin in ybins: 2622 if calendar_time: 2623 tbins[bin[0]] = list(gps2datenum(numpy.arange(float(start), float(end),\ 2624 average))) 2625 else: 2626 tbins[bin[0]] = list(numpy.arange(0,float(end-start), average)/unit) 2627 rate[bin[0]] = list(numpy.zeros(len(tbins[bin[0]]))) 2628 2629 for trig in triggers: 2630 x = int(float(getTrigAttribute(trig, 'time')-start)//average) 2631 y = getTrigAttribute(trig, bincolumn) 2632 for bin in ybins: 2633 if bin[0] <= y < bin[1]: 2634 if x>=len(rate[bin[0]]): 2635 print "trigger after end time, something is wrong", x, len(rate[bin[0]]) 2636 continue 2637 rate[bin[0]][x] += 1/average 2638 break 2639 2640 # if logscale includes zeros, pylab.scatter will break, so remove zeros 2641 if logy: 2642 for bin in ybins: 2643 removes = 0 2644 numtbins = len(tbins[bin[0]]) 2645 for rbin in xrange(0,numtbins): 2646 if rate[bin[0]][rbin-removes]==0: 2647 rate[bin[0]].pop(rbin-removes) 2648 tbins[bin[0]].pop(rbin-removes) 2649 removes+=1 2650 2651 # set labels 2652 etg = etg.replace('_', '\_') 2653 zero = LIGOTimeGPS('%.3f' % zero) 2654 if zero.nanoseconds==0: 2655 tlabel = datetime(*date.XLALGPSToUTC(LIGOTimeGPS(zero))[:6])\ 2656 .strftime("%B %d %Y, %H:%M:%S %ZUTC") 2657 else: 2658 tlabel = datetime(*date.XLALGPSToUTC(LIGOTimeGPS(zero.seconds))[:6])\ 2659 .strftime("%B %d %Y, %H:%M:%S %ZUTC") 2660 tlabel = tlabel.replace(' UTC', '.%.3s UTC' % zero.nanoseconds) 2661 xlabel = kwargs.pop('xlabel',\ 2662 'Time (%s) since %s (%s)' % (timestr, tlabel, zero)) 2663 ylabel = kwargs.pop('ylabel', 'Rate (Hz)') 2664 title = kwargs.pop('title', '%s triggers binned by %s'\ 2665 % (etg, display_name(bincolumn))) 2666 if start and end: 2667 subtitle = '%s-%s' % (start, end) 2668 else: 2669 subtitle = " " 2670 subtitle = kwargs.pop('subtitle', subtitle) 2671 2672 # customise plot appearance 2673 set_rcParams() 2674 2675 # generate plot object 2676 plot = ScatterPlot(xlabel, ylabel, title, subtitle) 2677 2678 # plot rates 2679 for bin in ybins: 2680 if logy: 2681 if len(rate[bin[0]])>0: 2682 plot.add_content(tbins[bin[0]], rate[bin[0]],\ 2683 label='-'.join(map(str, bin)), **kwargs) 2684 else: 2685 plot.add_content([1],[0.1], label='-'.join(map(str, bin)),\ 2686 visible=False) 2687 else: 2688 plot.add_content(tbins[bin[0]], rate[bin[0]], label='-'.join(map(str, bin)),\ 2689 **kwargs) 2690 2691 # finalise plot 2692 plot.finalize(logx=logx, logy=logy, hidden_colorbar=hidden_colorbar) 2693 2694 # set limits 2695 plot.ax.autoscale_view(tight=True, scalex=True, scaley=True) 2696 plot.ax.set_xlim(xlim) 2697 if ylim: 2698 plot.ax.set_ylim(ylim) 2699 2700 # normalize ticks 2701 set_ticks(plot.ax, calendar_time) 2702 2703 # save 2704 plot.savefig(outfile, bbox_inches=bbox_inches,\ 2705 bbox_extra_artists=plot.ax.texts) 2706 pylab.close(plot.fig)

2707

2708 # ============================================================================= 2709 # Plot RMS versus time in bins 2710 2711 -def plot_trigger_rms(triggers, outfile, average=600, start=None, end=None,\ 2712 zero=None, rmscolumn='snr', bincolumn='peak_frequency',\ 2713 bins=[], etg='Unknown', **kwargs):

2714 2715 """ 2716 Plot RMS versus time for the given ligolw table triggers, binned by the 2717 given bincolumn using the bins list. 2718 2719 Arguments: 2720 2721 triggers : glue.ligolw.table 2722 LIGOLW table containing a list of triggers 2723 outfile : string 2724 string path for output plot 2725 2726 Keyword arguments: 2727 2728 average : float 2729 Length (seconds) of RMS segment 2730 start : [ float | int | LIGOTimeGPS ] 2731 GPS start time 2732 end : [ float | int | LIGOTimeGPS ] 2733 GPS end time 2734 zero : [ float | int | LIGOTimeGPS ] 2735 GPS time to use for 0 on time axis 2736 rmscolumn : string 2737 valid column of the trigger table to RMS over 2738 bincolumn : string 2739 valid column of the trigger table to use for binning 2740 bins : list 2741 list of tuples defining the rate bins 2742 etg : string 2743 display name of trigger generator 2744 logy : [ True | False ] 2745 boolean option to display y-axis in log scale 2746 ylim : tuple 2747 (ymin, ymax) limits for rate axis 2748 """ 2749 2750 tableName = triggers.tableName.lower() 2751 get_time = def_get_time(tableName) 2752 2753 # set start and end times 2754 if not start and not end: 2755 times = [get_time(t) for t in triggers] 2756 if not start: 2757 start = min(times) 2758 if not end: 2759 end = max(times) 2760 2761 if not zero: 2762 zero = start 2763 2764 # set plot time unit whether it's used or not 2765 unit, timestr = time_unit(end-start) 2766 2767 # set ETG 2768 if not etg: 2769 if re.search('burst', tableName): 2770 etg = 'Burst' 2771 elif re.search('inspiral', tableName): 2772 etg = 'Inspiral' 2773 elif re.search('ringdown', tableName): 2774 etg = 'Ringdown' 2775 else: 2776 etg = 'Unknown' 2777 2778 # get limits 2779 calendar_time = kwargs.pop('calendar_time', False) 2780 if calendar_time: 2781 xlim = kwargs.pop('xlim', [gps2datenum(float(start)),\ 2782 gps2datenum(float(end))]) 2783 else: 2784 xlim = kwargs.pop('xlim', [float(start-zero)/unit, float(end-zero)/unit]) 2785 ylim = kwargs.pop('ylim', None) 2786 2787 # get axis scales 2788 logx = kwargs.pop('logx', False) 2789 logy = kwargs.pop('logy', False) 2790 2791 # get averaging time 2792 average = kwargs.pop('average',average) 2793 2794 # get colorbar options 2795 hidden_colorbar = kwargs.pop('hidden_colorbar', False) 2796 2797 # get savefig option 2798 bbox_inches = kwargs.pop('bbox_inches', 'tight') 2799 2800 # format ybins 2801 if not bins: 2802 bins = [[0,float('inf')]] 2803 ybins = [map(float, bin) for bin in bins] 2804 2805 # bin data 2806 tbins = {} 2807 rate = {} 2808 rms = {} 2809 for bin in ybins: 2810 if calendar_time: 2811 tbins[bin[0]] = list(gps2datenum(numpy.arange(float(start), float(end),\ 2812 average))) 2813 else: 2814 tbins[bin[0]] = list(numpy.arange(0,float(end-start), average)/unit) 2815 rate[bin[0]] = list(numpy.zeros(len(tbins[bin[0]]))) 2816 rms[bin[0]] = list(numpy.zeros(len(tbins[bin[0]]))) 2817 2818 for trig in triggers: 2819 x = int(float(getTrigAttribute(trig, 'time')-start)//average) 2820 y = getTrigAttribute(trig, bincolumn) 2821 z = getTrigAttribute(trig, rmscolumn) 2822 for bin in ybins: 2823 if bin[0] <= y < bin[1]: 2824 rms[bin[0]][x] += z*z 2825 rate[bin[0]][x] += 1 2826 break 2827 2828 # Normalize the RMS to get the mean not the sum 2829 for bin in ybins: 2830 for x in range(len(tbins[bin[0]])): 2831 if rate[bin[0]][x] : 2832 rms[bin[0]][x] = math.sqrt(rms[bin[0]][x]/rate[bin[0]][x]) 2833 2834 # if logscale includes zeros, pylab.scatter will break, so remove zeros 2835 if logy: 2836 for bin in ybins: 2837 removes = 0 2838 numtbins = len(tbins[bin[0]]) 2839 for rbin in xrange(0,numtbins): 2840 if rms[bin[0]][rbin-removes]==0: 2841 rms[bin[0]].pop(rbin-removes) 2842 tbins[bin[0]].pop(rbin-removes) 2843 removes+=1 2844 2845 # set labels 2846 etg = etg.replace('_', '\_') 2847 zero = LIGOTimeGPS('%.3f' % zero) 2848 if zero.nanoseconds==0: 2849 tlabel = datetime(*date.XLALGPSToUTC(LIGOTimeGPS(zero))[:6])\ 2850 .strftime("%B %d %Y, %H:%M:%S %ZUTC") 2851 else: 2852 tlabel = datetime(*date.XLALGPSToUTC(LIGOTimeGPS(zero.seconds))[:6])\ 2853 .strftime("%B %d %Y, %H:%M:%S %ZUTC") 2854 tlabel = tlabel.replace(' UTC', '.%.3s UTC' % zero.nanoseconds) 2855 xlabel = kwargs.pop('xlabel',\ 2856 'Time (%s) since %s (%s)' % (timestr, tlabel, zero)) 2857 ylabel = kwargs.pop('ylabel', 'RMS') 2858 title = kwargs.pop('title', '%s triggers binned by %s'\ 2859 % (etg, display_name(bincolumn))) 2860 if start and end: 2861 subtitle = '%s-%s' % (start, end) 2862 else: 2863 subtitle = " " 2864 subtitle = kwargs.pop('subtitle', subtitle) 2865 2866 # customise plot appearance 2867 set_rcParams() 2868 2869 # generms plot object 2870 plot = ScatterPlot(xlabel, ylabel, title, subtitle) 2871 2872 # plot rmss 2873 for bin in ybins: 2874 if logy: 2875 if len(rms[bin[0]])>0: 2876 plot.add_content(tbins[bin[0]], rms[bin[0]],\ 2877 label='-'.join(map(str, bin)), **kwargs) 2878 else: 2879 plot.add_content([1],[0.1], label='-'.join(map(str, bin)),\ 2880 visible=False) 2881 else: 2882 plot.add_content(tbins[bin[0]], rms[bin[0]], label='-'.join(map(str, bin)),\ 2883 **kwargs) 2884 2885 # finalise plot 2886 plot.finalize(logx=logx, logy=logy, hidden_colorbar=hidden_colorbar) 2887 2888 # set limits 2889 plot.ax.autoscale_view(tight=True, scalex=True, scaley=True) 2890 plot.ax.set_xlim(xlim) 2891 if ylim: 2892 plot.ax.set_ylim(ylim) 2893 2894 # normalize ticks 2895 set_ticks(plot.ax, calendar_time=calendar_time) 2896 2897 # save 2898 plot.savefig(outfile, bbox_inches=bbox_inches,\ 2899 bbox_extra_artists=plot.ax.texts) 2900 pylab.close(plot.fig)

2901

2902 # ============================================================================= 2903 # Plot segments 2904 2905 -def plot_segments(segdict, outfile, start=None, end=None, zero=None, 2906 keys=None, highlight_segments=None, **kwargs):

2907 2908 """ 2909 Plot the segments contained within the glue.segments.segmentlistdict 2910 segdict to the given path string outfile. The list keys can be given to 2911 guarantee the order of the segments on the y-axis. x-axis limits can be 2912 controlled using start, end and zero. The glue.segments.segmentlist object 2913 highlight_segments can be given to highlight a number of segments. 2914 2915 Arguments: 2916 2917 segdict : glue.segments.segmentlistdict 2918 2919 """ 2920 2921 if not start: 2922 minstart = lambda seglist: min(segdict[key][0][0] for key in segdict.keys()\ 2923 if segdict[key]) 2924 start = min(minstart(segdict[key]) for key in segdict.keys()) 2925 if not end: 2926 maxend = lambda seglist: max(segdict[key][0][1] for key in segdict.keys()\ 2927 if segdict[key]) 2928 end = max(maxend(segdict[key]) for key in segdict.keys()) 2929 if not zero: 2930 zero = start 2931 2932 # set plot time unit whether it's used or not 2933 unit, timestr = time_unit(end-start) 2934 2935 # set labels 2936 zero = LIGOTimeGPS('%.3f' % zero) 2937 if zero.nanoseconds==0: 2938 tlabel = datetime(*date.XLALGPSToUTC(LIGOTimeGPS(zero))[:6])\ 2939 .strftime("%B %d %Y, %H:%M:%S %ZUTC") 2940 else: 2941 tlabel = datetime(*date.XLALGPSToUTC(LIGOTimeGPS(zero.seconds))[:6])\ 2942 .strftime("%B %d %Y, %H:%M:%S %ZUTC") 2943 tlabel = tlabel.replace(' UTC', '.%.3s UTC' % zero.nanoseconds) 2944 xlabel = kwargs.pop('xlabel',\ 2945 'Time (%s) since %s (%s)' % (timestr, tlabel, zero)) 2946 ylabel = kwargs.pop('ylabel', "") 2947 title = kwargs.pop('title', '') 2948 subtitle = kwargs.pop('subtitle', '') 2949 2950 # get axis limits 2951 calendar_time = kwargs.pop('calendar_time', False) 2952 if calendar_time: 2953 xlim = kwargs.pop('xlim', [gps2datenum(float(start)),\ 2954 gps2datenum(float(end))]) 2955 else: 2956 xlim = kwargs.pop('xlim', [float(start-zero)/unit, float(end-zero)/unit]) 2957 2958 # get label param 2959 labels_inset = kwargs.pop('labels_inset', False) 2960 2961 # get savefig option 2962 bbox_inches = kwargs.pop('bbox_inches', 'tight') 2963 2964 # get colorbar option 2965 hidden_colorbar = kwargs.pop('hidden_colorbar', False) 2966 2967 if keys: 2968 # escape underscore, but don't do it twice 2969 keys = [key.replace('_','\_').replace('\\_','\_') for key in keys] 2970 segkeys = segdict.keys() 2971 newdict = segments.segmentlistdict() 2972 for key in segkeys: 2973 newkey = key.replace('_','\_').replace('\\\_','\_') 2974 newdict[newkey] = segdict[key] 2975 segdict = newdict 2976 2977 # set params 2978 set_rcParams() 2979 2980 plot = PlotSegmentsPlot(xlabel, ylabel, title, subtitle, t0=zero, unit=unit,\ 2981 calendar_time=calendar_time) 2982 plot.add_content(segdict, keys, **kwargs) 2983 plot.finalize(labels_inset=labels_inset, hidden_colorbar=hidden_colorbar) 2984 2985 # indicate last frame 2986 if highlight_segments: 2987 for seg in highlight_segments: 2988 plot.highlight_segment(seg) 2989 2990 # set x axis 2991 plot.ax.set_xlim(xlim) 2992 2993 plot.ax.grid(True,which='major') 2994 plot.ax.grid(True,which='majorminor') 2995 2996 set_ticks(plot.ax, calendar_time) 2997 2998 plot.savefig(outfile, bbox_inches=bbox_inches,\ 2999 bbox_extra_artists=plot.ax.texts) 3000 pylab.close(plot.fig)

3001

3002 # ============================================================================= 3003 # Helper functions 3004 # ============================================================================= 3005 3006 -def parse_plot_config(cp, section):

3007 3008 """ 3009 Parse ConfigParser.ConfigParser section for plot parameters. Sections should 3010 be name '[plot xcolumn-ycolumn-zcolumn]' e.g. 3011 '[plot time-peak_frequency-snr]'. Returns a pair of dicts with the 3012 following keys: 3013 3014 columns: 3015 3016 xcolumn : [ string | None ] 3017 column string to plot on x-axis 3018 ycolumn : [ string | None ] 3019 column string to plot on y-axis 3020 zcolumn : [ string | None ] 3021 column string to plot on z-axis 3022 3023 params: 3024 3025 xlim : list 3026 [xmin, xmax] pair for x-axis limits 3027 ylim : list 3028 [ymin, ymax] pair for y-axis limits 3029 zlim : list 3030 [zmin, zmax] pair for z-axis limits 3031 clim : list 3032 [cmin, cmax] pair for colorbar limits 3033 logx : bool 3034 True / False to plot log scale on x-axis 3035 logy : bool 3036 True / False to plot log scale on y-axis 3037 logz : bool 3038 True / False to plot log scale on z-axis 3039 """ 3040 3041 columns = {'xcolumn':None, 'ycolumn':None, 'zcolumn':None, 'rankcolumn':None} 3042 params = {} 3043 3044 plot = re.split('[\s-]', section)[1:] 3045 if len(plot)>=1: 3046 columns['xcolumn'] = plot[0] 3047 if len(plot)>=2: 3048 columns['ycolumn'] = plot[1] 3049 if len(plot)>=3: 3050 columns['zcolumn'] = plot[2] 3051 if len(plot)>=4: 3052 columns['rankcolumn'] = plot[3] 3053 3054 limits = ['xlim', 'ylim', 'zlim', 'clim', 'exponents', 'constants',\ 3055 'color_bins'] 3056 filters = ['poles', 'zeros'] 3057 bins = ['bins'] 3058 booleans = ['logx', 'logy', 'logz', 'cumulative', 'rate', 'detchar',\ 3059 'greyscale', 'zeroindicator', 'normalized', 'include_downtime',\ 3060 'calendar_time', 'fill', 'hidden_colorbar'] 3061 values = ['dcthresh','amplitude','num_bins','linewidth','average','s'] 3062 3063 # extract plot params as a dict 3064 params = {} 3065 for key,val in cp.items(section): 3066 val = val.rstrip('"').strip('"') 3067 if key in limits: 3068 params[key] = map(float, val.split(',')) 3069 elif key in filters: 3070 params[key] = map(complex, val.split(',')) 3071 elif key in bins: 3072 params[key] = map(lambda p: map(float, p.split(',')), val.split(';')) 3073 elif key in booleans: 3074 params[key] = cp.getboolean(section, key) 3075 elif key in values: 3076 params[key] = float(val) 3077 else: 3078 params[key] = val 3079 3080 return columns, params

3081

3082 # ============================================================================== 3083 # Plot color map 3084 3085 -def plot_color_map(data, outfile, data_limits=None, x_format='time',\ 3086 y_format='frequency', z_format='amplitude', zero=None,\ 3087 x_range=None, y_range=None, **kwargs):

3088 3089 """ 3090 Plots data in a 2d color map. 3091 """ 3092 3093 if not (isinstance(data, list) or isinstance(data, tuple)): 3094 data_sets = [data] 3095 data_limit_sets = [data_limits] 3096 else: 3097 data_sets = data 3098 data_limit_sets = data_limits 3099 if not len(data_sets) == len(data_limit_sets): 3100 raise AttributeError("You have given %d data sets and %d limit sets! Eh?"\ 3101 % (len(data_sets), len(data_limit_sets))) 3102 3103 # set data limits 3104 for i,data_limits in enumerate(data_limit_sets): 3105 if not data_limits: 3106 numrows, numcols = numpy.shape(data_sets[i]) 3107 if kwargs.has_key('xlim'): 3108 xlim = kwargs['xlim'] 3109 else: 3110 xlim = [0, numrows] 3111 if kwargs.has_key('ylim'): 3112 ylim = kwargs['ylim'] 3113 elif pylab.rcParams['image.origin'] == 'upper': 3114 ylim = [numrows, 0] 3115 else: 3116 ylim = [0, numrows] 3117 3118 data_limit_sets[i] = [xlim[0], xlim[1], ylim[0], ylim[1]] 3119 3120 # get limits 3121 xlim = kwargs.pop('xlim', None) 3122 ylim = kwargs.pop('ylim', None) 3123 zlim = kwargs.pop('zlim', None) 3124 clim = kwargs.pop('clim', None) 3125 calendar_time = kwargs.pop('calendar_time', False) 3126 3127 # get axis scales 3128 logx = kwargs.pop('logx', False) 3129 logy = kwargs.pop('logy', False) 3130 logz = kwargs.pop('logz', False) 3131 3132 # get savefig option 3133 bbox_inches = kwargs.pop('bbox_inches', 'tight') 3134 3135 # restrict data to meet limits if using log 3136 for i,(data, data_limits) in enumerate(zip(data_sets, data_limit_sets)): 3137 if logx and xlim: 3138 shape = numpy.shape(data) 3139 xmin, xmax = data_limits[0], data_limits[1] 3140 if x_range==None: 3141 x_range = numpy.logspace(numpy.log10(xmin), numpy.log10(xmax),\ 3142 num.shape[-1]) 3143 condition = (x_range>xlim[0]) & (x_range<=xlim[1]) 3144 newx = numpy.where(condition)[0] 3145 data2 = numpy.resize(data, (len(newx), shape[-2])) 3146 for j in xrange(newx): 3147 data2[:,j] = data[newx[j],:] 3148 data = data2.transpose() 3149 3150 if logy and ylim: 3151 shape = numpy.shape(data) 3152 ymin, ymax = data_limits[2], data_limits[3] 3153 if y_range==None: 3154 y_range = numpy.logspace(numpy.log10(ymin), numpy.log10(ymax),\ 3155 num=shape[-2]) 3156 condition = (y_range>ylim[0]) & (y_range<=ylim[1]) 3157 newy = numpy.where((condition))[0] 3158 data2 = numpy.resize(data, (len(newy), shape[-1])) 3159 for j in xrange(shape[-1]): 3160 data2[:,j] = data[:,j][condition] 3161 data = data2 3162 data_limits[2:] = ylim 3163 data_sets[i] = data 3164 3165 # get columnar params 3166 columns = list(map(str.lower, [x_format, y_format, z_format])) 3167 limits = [xlim, ylim, zlim] 3168 for i,lim in enumerate(limits): 3169 if lim: 3170 limits[i] = list(map(float, lim)) 3171 labels = ["", "", "", "", ""] 3172 3173 # get zero time for normalisation 3174 if 'time' in columns and not zero: 3175 i = columns.index('time') 3176 if limits[i]: 3177 zero = limits[i][0] 3178 else: 3179 zero = min(data_limits[0] for data_limits in data_limit_sets) 3180 3181 # format labels 3182 for i,(col,c) in enumerate(zip(columns, ['x', 'y', 'z'])): 3183 if col.lower() == 'time' and limits[i]: 3184 unit, timestr = time_unit(float(limits[i][1]-limits[i][0])) 3185 zero = LIGOTimeGPS('%.3f' % zero) 3186 if zero.nanoseconds==0: 3187 tlabel = datetime(*date.XLALGPSToUTC(LIGOTimeGPS(zero))[:6])\ 3188 .strftime("%B %d %Y, %H:%M:%S %ZUTC") 3189 else: 3190 tlabel = datetime(*date.XLALGPSToUTC(LIGOTimeGPS(zero.seconds))[:6])\ 3191 .strftime("%B %d %Y, %H:%M:%S %ZUTC") 3192 tlabel = tlabel.replace(' UTC', '.%.3s UTC' % zero.nanoseconds) 3193 labels[i] = kwargs.pop('%slabel' % c, 'Time (%s) since %s (%s)'\ 3194 % (timestr, tlabel, zero)) 3195 if calendar_time: 3196 limits[i] = gps2datenum(numpy.asarray(limits[i])) 3197 else: 3198 limits[i] = (numpy.asarray(limits[i])-float(zero))/unit 3199 for i,data_limits in enumerate(data_limit_sets): 3200 if calendar_time: 3201 data_limit_sets[i][0] = gps2datenum(float(data_limits[0])) 3202 data_limit_sets[i][1] = gps2datenum(float(data_limits[1])) 3203 else: 3204 data_limit_sets[i][0] = float(data_limits[0]-zero)/unit 3205 data_limit_sets[i][1] = float(data_limits[1]-zero)/unit 3206 else: 3207 labels[i] = kwargs.pop('%slabel' % c, display_name(col)) 3208 3209 labels[3] = kwargs.pop('title', "") 3210 labels[4] = kwargs.pop('subtitle', "") 3211 3212 # customise plot appearance 3213 set_rcParams() 3214 3215 # generate plot object 3216 plot = ColorMap(*labels) 3217 3218 # add data 3219 for i, (data, data_limits) in enumerate(zip(data_sets, data_limit_sets)): 3220 plot.add_content(data, data_limits, aspect='auto') 3221 3222 # finalize 3223 plot.finalize(logx=logx, logy=logy, logz=logz, clim=clim, origin='lower') 3224 3225 if limits[0] is not None and len(limits[0])==2: 3226 plot.ax.set_xlim(limits[0]) 3227 if limits[1] is not None and len(limits[1])==2: 3228 plot.ax.set_ylim(limits[1]) 3229 3230 # set global ticks 3231 set_ticks(plot.ax, calendar_time=calendar_time) 3232 3233 # save figure 3234 plot.savefig(outfile, bbox_inches=bbox_inches,\ 3235 bbox_extra_artists=plot.ax.texts) 3236 pylab.close(plot.fig)

3237

3238 # ============================================================================= 3239 # Significance drop plot (HVeto style) 3240 3241 -def plot_significance_drop(startsig, endsig, outfile, **kwargs):

3242 3243 """ 3244 Plot significance drop for each channel relative to the application of 3245 HVeto round veto segments. 3246 """ 3247 3248 # get channels 3249 channels = startsig.keys() 3250 for c in channels: 3251 if c not in endsig.keys(): 3252 raise AttributeError("Significance lists do not match.") 3253 channels.sort() 3254 3255 # find winner 3256 wch,wsig = max(startsig.items(), key=lambda x: x[1]) 3257 3258 # extract parameters 3259 kwargs.pop('xlim', None) 3260 ylim = kwargs.pop('ylim', (0, wsig+1)) 3261 xlabel = kwargs.pop('xlabel', "") 3262 ylabel = kwargs.pop('ylabel', "Significance") 3263 title = kwargs.pop('title', "Coincidence significance drop plot") 3264 subtitle = kwargs.pop('subtitle',\ 3265 "Winner: %s, significance: %s" % (wch, wsig)) 3266 3267 kwargs.setdefault('linestyle', '-') 3268 kwargs.setdefault('marker', 'o') 3269 color = kwargs.pop('color', None) 3270 3271 # customise plot appearance 3272 set_rcParams() 3273 pylab.rcParams.update({"figure.figsize":[24,6], "xtick.labelsize": 8}) 3274 3275 # get colorbar options 3276 hidden_colorbar = kwargs.pop('hidden_colorbar', False) 3277 3278 # get savefig option 3279 bbox_inches = kwargs.pop('bbox_inches', 'tight') 3280 3281 # generate plot object 3282 plot = DataPlot(xlabel, ylabel, title, subtitle) 3283 3284 kwargs.setdefault('markersize', 10) 3285 3286 # plot each channel's drop 3287 for i,c in enumerate(channels): 3288 s = startsig[c] 3289 e = endsig[c] 3290 col = color and color or s>e and 'b' or 'r' 3291 plot.add_content([i,i], [s,e], color=col, **kwargs) 3292 3293 # finalise plot object 3294 plot.finalize(logx=False, logy=False, hidden_colorbar=hidden_colorbar) 3295 3296 # set xticks to channel names and rotate 3297 plot.ax.set_xlim(-1, len(channels)) 3298 plot.ax.set_xticks(numpy.arange(0,len(channels))) 3299 plot.ax.set_xticklabels([c.replace('_','\_') for c in channels]) 3300 for i,t in enumerate(plot.ax.get_xticklabels()): 3301 t.set_rotation(315) 3302 #t.set_position(((i+1)/(len(channels)+1), 0.1)) 3303 t.set_verticalalignment('top') 3304 t.set_horizontalalignment('left') 3305 #t.set_bbox(dict(facecolor='white', alpha=0.5, edgecolor='none')) 3306 3307 # set ylim 3308 plot.ax.set_ylim(ylim) 3309 3310 # turn off x grid 3311 plot.ax.xaxis.grid(False) 3312 3313 plot.fig.subplots_adjust(bottom=0.3) 3314 3315 # save figure 3316 plot.savefig(outfile, bbox_inches=bbox_inches,\ 3317 bbox_extra_artists=plot.ax.texts) 3318 pylab.close(plot.fig)

3319

3320 # ============================================================================= 3321 # Plot sky positions 3322 3323 -def plot_sky_positions(skyTable, outfile, format='radians', zcolumn=None,\ 3324 projection=None, centre=None, detectors=[], lines={},\ 3325 range=[(0,0), (1,1)], **kwargs):

3326 3327 """ 3328 Plot latitude against longitude for the given ligolw table skyTable into the 3329 given outfile. Uses the mpl_toolkits basemap module to plot the sky sphere 3330 in a variety of projections, or simply a scatter plot if projection=None. 3331 Can plot lines and detector positions on top if given. 3332 3333 Arguments: 3334 3335 skyTable : glue.ligolw.table.Table 3336 ligolw table containing triggers or SkyPosition objects 3337 outfile : string 3338 string path for output plot 3339 3340 Keyword arguments: 3341 3342 zcolumn : string 3343 valid column of ligolw table to use for colorbar (optional). 3344 format : [ 'radians' | 'degrees' ] 3345 str identifying format of longtiude/ra, latitude/dec colums in table. 3346 Plot is always drawn in degrees. 3347 projection : str 3348 type of spherical projection to use, if any. See matplotlib Basemap 3349 documentation for details, recommended: 'ortho', 'hammer'. 3350 centre : tuple 3351 (longitude, latitude) pair on which to centre plot. Latitude centring 3352 only works for certain projections. 3353 detectors : list 3354 list of detector prefixes to plot, e.g. ['H1', 'L1']. 3355 lines : dict 3356 dict of name:table pairs from which to plot lines (+ marker) on top of 3357 points 3358 range : tuple 3359 ((xmin, ymin), (xmax, ymax)) tuples for lower left and upper right 3360 corners, in range 0-1, e.g. (0,0) is full lower left, (1,1) full upper 3361 right corners for full spherical projection. Only works with certain 3362 projections. 3363 3364 Unnamed keyword arguments: 3365 3366 logx : [ True | False ] 3367 boolean option to display x-axis in log scale. Not applicable when using 3368 projection. 3369 logy : [ True | False ] 3370 boolean option to display y-axis in log scale. Not applicable when using 3371 projection. 3372 logz : [ True | False ] 3373 boolean option to display z-axis in log scale. 3374 xlim : tuple 3375 (xmin, xmax) limits for x-axis (longitude). Triggers outside range are 3376 removed. 3377 ylim : tuple 3378 (ymin, ymax) limits for y-axis (latitude). Triggers outside range are 3379 removed. 3380 zlim : tuple 3381 (zmin, zmax) limits for z-axis. Triggers outside range are removed. 3382 clim : tuple 3383 (cmin, cmax) limits for color scale. Triggers outside range are moved 3384 onto boundary. 3385 xlabel : string 3386 label for x-axis 3387 ylabel : string 3388 label for y-axis 3389 zlabel : string 3390 label for z-axis 3391 title : string 3392 title for plot 3393 subtitle : string 3394 subtitle for plot 3395 3396 All other given arguments will be passed to matplotlib.axes.Axes.scatter. 3397 """ 3398 3399 format = format.lower() 3400 3401 # get labels 3402 if projection: 3403 xlabel = kwargs.pop("xlabel", "") 3404 ylabel = kwargs.pop("ylabel", "") 3405 else: 3406 xlabel = kwargs.pop("xlabel", "Longitude (degrees of arc)") 3407 ylabel = kwargs.pop("ylabel", "Latitude (degrees of arc)") 3408 if zcolumn: tmp = zcolumn 3409 else: tmp = "" 3410 zlabel = kwargs.pop("zlabel", display_name(tmp)) 3411 3412 # get ETG 3413 if re.search('burst', skyTable.tableName.lower()): 3414 etg = 'Burst' 3415 elif re.search('inspiral', skyTable.tableName.lower()): 3416 etg = 'Inspiral' 3417 elif re.search('ringdown', skyTable.tableName.lower()): 3418 etg = 'Ringdown' 3419 else: 3420 etg = 'Unknown' 3421 etg = etg.replace('_', '\_') 3422 3423 # get title 3424 title = kwargs.pop("title", "") 3425 subtitle = kwargs.pop("subtitle", "") 3426 3427 # get limits 3428 xlim = kwargs.pop('xlim', None) 3429 ylim = kwargs.pop('ylim', None) 3430 zlim = kwargs.pop('zlim', None) 3431 clim = kwargs.pop('clim', None) 3432 3433 # get logs 3434 logx = kwargs.pop('logx', False) 3435 logy = kwargs.pop('logy', False) 3436 logz = kwargs.pop('logz', False) 3437 3438 # get savefig option 3439 bbox_inches = kwargs.pop('bbox_inches', 'tight') 3440 3441 # get IFO data 3442 detData = [] 3443 detCol = [] 3444 for ifo in detectors: 3445 d = XLALTools.cached_detector[inject.prefix_to_name[ifo]] 3446 if format=='radians': 3447 detData.append((numpy.degrees(d.vertexLongitudeRadians),\ 3448 numpy.degrees(d.vertexLatitudeRadians))) 3449 else: 3450 detData.append((d.vertexLongitudeRadians, d.vertexLatitudeRadians)) 3451 if ifo in plotsegments.PlotSegmentsPlot.color_code.keys(): 3452 detCol.append(plotsegments.PlotSegmentsPlot.color_code[ifo]) 3453 else: 3454 detCol.append('cyan') 3455 3456 # get line data 3457 3458 lineData = [] 3459 lineCol = [] 3460 for line in lines: 3461 if isinstance(lines, dict): 3462 label = line 3463 line = lines[label] 3464 else: 3465 label = '_' 3466 # get column 3467 if re.search('multi_inspiral', line.tableName): 3468 loncol = 'ra' 3469 latcol = 'dec' 3470 else: 3471 loncol = 'longitude' 3472 latcol = 'latitude' 3473 # get data 3474 lon = [] 3475 lat = [] 3476 for p in line: 3477 if format=='radians': 3478 lon.append(numpy.degrees(getattr(p, loncol))) 3479 lat.append(numpy.degrees(getattr(p, latcol))) 3480 else: 3481 lon.append(getattr(p, loncol)) 3482 lat.append(getattr(p, latcol)) 3483 lineData.append((label, [lon, lat])) 3484 if len(lineData): 3485 for c in plotutils.default_colors(): 3486 if c=='b': continue 3487 lineCol.append(c) 3488 if len(lineCol)==len(lineData): break 3489 3490 # get point data 3491 columns = ['longitude', 'latitude'] 3492 if re.search('multi_inspiral', skyTable.tableName): 3493 columns[0] = 'ra' 3494 columns[1] = 'dec' 3495 if zcolumn: 3496 columns.append(zcolumn) 3497 3498 data = {} 3499 for col in columns: 3500 data[col] = [] 3501 for i,p in enumerate(skyTable): 3502 # test limits 3503 if xlim and not xlim[0] <= getTrigAttribute(p, columns[0]) <= xlim[1]: 3504 continue 3505 if ylim and not ylim[0] <= getTrigAttribute(p, columns[1]) <= ylim[1]: 3506 continue 3507 if zlim and not zlim[0] <= getTrigAttribute(p, columns[2]) <= zlim[1]: 3508 continue 3509 # get sky data 3510 for col in columns[0:2]: 3511 if format=='radians': 3512 data[col].append(numpy.degrees(getTrigAttribute(p, col))) 3513 else: 3514 data[col].append(getTrigAttribute(p, col)) 3515 # get z-axis data 3516 if zcolumn: 3517 data[columns[2]].append(getTrigAttribute(p, columns[2])) 3518 3519 # customise plot appearance 3520 set_rcParams() 3521 if projection not in [None, 'hammer', 'moll', 'robin', 'sinu', 'cyl', 'merc',\ 3522 'mill', 'gall']: 3523 pylab.rcParams.update({"figure.figsize":[8,6]}) 3524 3525 # generate plot 3526 if projection: 3527 # colorbar plot 3528 if zcolumn: 3529 plot = ColorbarSkyPositionsPlot(xlabel, ylabel, zlabel,\ 3530 title, subtitle) 3531 plot.add_content(data[columns[0]], data[columns[1]], data[zcolumn],\ 3532 **kwargs) 3533 for p,ifo,c in zip(detData, detectors, detCol): 3534 lon, lat = p 3535 plot.add_content([lon], [lat], [1], label=ifo, facecolor=c,\ 3536 edgecolor=c, **kwargs) 3537 for line,c in zip(lineData, lineCol): 3538 label, data = line 3539 plot.add_content(data[0], data[1], label=label, marker='+',\ 3540 edgecolor=c, s=4) 3541 plot.finalize(projection=projection, centre=centre, range=range,\ 3542 logz=logz, clim=clim) 3543 # normal plot 3544 else: 3545 plot = SkyPositionsPlot(xlabel, ylabel, title, subtitle) 3546 plot.add_content(data[columns[0]], data[columns[1]], label='_', **kwargs) 3547 for p,ifo,c in zip(detData, detectors, detCol): 3548 lon, lat = p 3549 plot.add_content([lon], [lat], label=ifo, facecolor=c, edgecolor=c,\ 3550 **kwargs) 3551 for line,c in zip(lineData, lineCol): 3552 label, data = line 3553 plot.add_content(data[0], data[1], label=label, marker='+',\ 3554 edgecolor=c, s=4) 3555 plot.finalize(projection=projection, centre=centre, range=range) 3556 3557 # generate square plot 3558 else: 3559 if zcolumn: 3560 plot = ColorbarScatterPlot(xlabel, ylabel, zlabel, title,\ 3561 subtitle) 3562 plot.add_content(data[columns[0]], data[columns[1]], data[columns[2]],\ 3563 **kwargs) 3564 plot.finalize(logx=logx, logy=logy, logz=logz, clim=clim) 3565 else: 3566 plot = ScatterPlot(xlabel, ylabel, title, subtitle) 3567 plot.add_content(data[columns[0]], data[columns[1]], **kwargs) 3568 plot.finalize(logx=logx, logy=logy) 3569 3570 # set axes 3571 plot.ax.autoscale_view(tight=True, scalex=True, scaley=True) 3572 plot.ax.set_xlim(xlim) 3573 plot.ax.set_ylim(ylim) 3574 3575 # save 3576 plot.savefig(outfile, bbox_inches=bbox_inches,\ 3577 bbox_extra_artists=plot.ax.texts) 3578 pylab.close(plot.fig)

3579

Source Code for Module pylal.dq.dqPlotUtils