1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19 """
20 This modules provides some wrapper functions for the classes defined in pylal.plotutils.
21 """
22
23
24
25
26
27 from __future__ import division
28 import re
29 import copy
30 import numpy
31 import datetime
32 from pylal import git_version
33 from pylal import date
34 from pylal import plotutils
35 from pylal.datatypes import LIGOTimeGPS
36 from glue.ligolw import lsctables
37 from glue.ligolw import table
38
39 pylab = plotutils.pylab
40 matplotlib = pylab.matplotlib
41
42 __author__ = "Duncan M. Macleod <duncan.macleod@ligo.org>"
43 __version__ = git_version.id
44 __date__ = git_version.date
45
46
47 jet = copy.deepcopy(matplotlib.cm.jet)
48
49 jet._segmentdata['red'] = ([x == jet._segmentdata["red"][-1] and (1,0.8,0.8)\
50 or x for x in jet._segmentdata['red'] ])
51 jet._segmentdata['blue'] = ([x == jet._segmentdata["blue"][0] and (0,0.7,0.7)\
52 or x for x in jet._segmentdata['blue'] ])
53
54
55
56
57
59 """
60 Extract column from the given glue.ligolw.table lsctable as numpy array.
61 Tries to use a 'get_col() function if available, otherwise uses
62 getColumnByName(), treating 'time' as a special case for the known
63 Burst/Inspiral/Ringdown tables.
64 """
65
66 column = str(column).lower()
67 obj_type = str(type(lsctable))
68
69
70 if hasattr(lsctable, 'get_%s' % column):
71 return numpy.asarray(getattr(lsctable, 'get_%s' % column)())
72
73
74 if column == 'time':
75
76 if re.match("sim_inspiral", lsctable.tableName, re.I):
77 timecolumn = "geocent_end_time"
78 elif re.match("sim_burst", lsctable.tableName, re.I):
79 timecolumn = "time_geocent_gps"
80 elif re.match("sim_ringdown", lsctable.tableName, re.I):
81 timecolumn = "geocent_start_time"
82 elif re.search("inspiral", lsctable.tableName, re.I):
83 timecolumn = "end_time"
84 elif re.search("burst", lsctable.tableName, re.I):
85 timecolumn = "peak_time"
86 elif re.search("(ringdown|stochastic)", lsctable.tableName, re.I):
87 timecolumn = "start_time"
88
89
90 return numpy.asarray(lsctable.getColumnByName(timecolumn)) + \
91 numpy.asarray(lsctable.getColumnByName("%s_ns"%timecolumn))*1e-9
92
93 return numpy.asarray(lsctable.getColumnByName(column))
94
95
96
97
98
99 -def plottable(lsctable, outfile, xcolumn="time", ycolumn="snr",\
100 colorcolumn=None, rankcolumn=None, t0=None, **kwargs):
101 """
102 Plot any column of a valid LigoLW table against any other, coloured by any
103 other, or plot all the same if you're that way inclined. Multiple tables
104 can be given for any non-coloured plot.
105
106 "time" as a column means the output of get_peak for Burst tables, get_end
107 for Inspiral tables, and get_start for ringdown tables
108
109 Arguments:
110
111 tables : [ dict | glue.ligolw.table.Table ]
112 dict of ("name", table) pairs, or single LigoLW tables
113 outfile : string
114 string path for output plot
115
116 Keyword arguments:
117
118 xcolumn : string
119 valid column of triggers table to plot on x-axis
120 ycolumn : string
121 valid column of triggers table to plot on y-axis
122 zcolumn : string
123 valid column of triggers table to use for colorbar (optional).
124 """
125
126 if isinstance(lsctable, table.Table):
127 tables = {"_":lsctable}
128 else:
129 tables = lsctable
130 tablenames = tables.keys()
131 tables = [tables[n] for n in tablenames]
132
133
134 xlim = kwargs.pop('xlim', None)
135 ylim = kwargs.pop('ylim', None)
136 zlim = kwargs.pop('zlim', None)
137 colorlim = kwargs.pop('colorlim', None)
138 if zlim and not colorlim:
139 colorlim = zlim
140
141
142 columns = list(map(str.lower, [xcolumn, ycolumn]))
143 if colorcolumn:
144 columns.append(colorcolumn.lower())
145 if rankcolumn:
146 columns.append(rankcolumn.lower())
147 else:
148 columns.append(colorcolumn.lower())
149
150
151 limits = [xlim, ylim, zlim, None]
152
153
154 if "time" in columns and not t0:
155 if xlim:
156 t0 = float(xlim[0])
157 else:
158 t0 = numpy.infty
159 for t in tables:
160 timedata = get_column(t, "time").astype(float)
161 if len(timedata):
162 t0 = min(t0, timedata.min())
163 if numpy.isinf(t0):
164 t0 = 0
165 t0 = int(t0)
166
167
168
169
170
171
172 data = list()
173 for i,name in enumerate(tablenames):
174 data.append(list())
175 for c in columns:
176 c = c.lower()
177 if ((c not in tables[i].columnnames and c != "time") and not
178 hasattr(tables[i], "get_%s" % c.lower())):
179 raise RuntimeError("Column '%s' not found in table '%s'."\
180 % (c,name))
181 data[-1].append(get_column(tables[i], c.lower()).astype(float))
182 if not len(data[-1][-1].shape)\
183 or data[-1][-1].shape[0] != len(tables[i]):
184 raise AttributeError("No data loaded for column \"%s\" and "\
185 "table \"%s\"" % (c, name))
186
187
188 for i in range(len(columns)):
189 if not limits[i]:
190 mins = [data[j][i].min() for j in range(len(tablenames))\
191 if len(data[j][i].shape) and data[j][i].shape[0] != 0]
192 if len(mins):
193 lmin = min(mins)
194 lmax = max(data[j][i].max() for j in range(len(tablenames))\
195 if len(data[j][i]))
196 limits[i] = [lmin, lmax]
197
198
199 if "time" in columns:
200 idx = columns.index("time")
201 if limits[idx]:
202 unit, timestr = plotutils.time_axis_unit(limits[idx][1]\
203 - limits[idx][0])
204 else:
205 unit, timestr = plotutils.time_axis_unit(1)
206
207
208 for i,name in enumerate(tablenames):
209 for j,column in enumerate(columns):
210 if column == "time":
211 data[i][j] = (data[i][j] - t0) / unit
212 if i==0 and limits[j]:
213 limits[j] = [float(limits[j][0] - t0) / unit,\
214 float(limits[j][1] - t0) / unit]
215
216
217 plotted = True
218 for j,column in enumerate(columns):
219 if limits[j]:
220 plotted = plotted & (limits[j][0] <= data[i][j])\
221 & (limits[j][1] >= data[i][j])
222
223
224 if not isinstance(plotted, bool):
225 for j,d in enumerate(data[i]):
226 data[i][j] = d[plotted]
227
228
229
230
231
232 loudest = None
233 if len(columns) == 4 and len(tablenames) == 1 and len(data[0][0]) > 0:
234 idx = data[0][3].argmax()
235 loudest = [data[0][j][idx] for j in range(len(columns))]
236
237
238
239
240
241 label = {}
242 for i,(label,column) in enumerate(zip(["xlabel", "ylabel", "colorlabel"],\
243 columns)):
244
245 l = kwargs.pop(label, None)
246 if l is None:
247
248 if column == "time" and limits[i]:
249 zerostr = datetime.datetime(\
250 *date.XLALGPSToUTC(LIGOTimeGPS(t0))[:6])\
251 .strftime("%B %d %Y, %H:%M:%S %ZUTC")
252 l = "Time (%s) since %s (%s)" % (timestr, zerostr, t0)
253
254 else:
255 l = plotutils.display_name(column)
256 if label == "xlabel":
257 xlabel = l
258 elif label == "ylabel":
259 ylabel = l
260 else:
261 colorlabel = l
262
263 title = kwargs.pop("title", "")
264 subtitle = kwargs.pop("subtitle", None)
265 if subtitle is None and loudest:
266 subtitle = "Loudest event by %s:" % plotutils.display_name(columns[-1])
267 for j,c in enumerate(columns):
268 if j!= 0 and c == columns[j-1]: continue
269 lstr = loudest[j]
270 if c == "time":
271 lstr = lstr * unit + t0
272 subtitle += " %s=%.2f" % (plotutils.display_name(c), lstr)
273 else:
274 loudest = None
275
276
277
278
279
280
281 logx = kwargs.pop("logx", False)
282 logy = kwargs.pop("logy", False)
283 logcolor = kwargs.pop("logcolor", False)
284
285
286 loc = kwargs.pop("loc", 0)
287
288
289 hidden_colorbar = kwargs.pop("hidden_colorbar", False)
290
291
292 bbox_inches = kwargs.pop("bbox_inches", "tight")
293
294
295 dqstyle = (kwargs.pop("detchar", False) or
296 kwargs.pop('detchar_style', False))
297 dqthresh = (kwargs.pop('dcthreshold', False) or
298 kwargs.pop('detchar_style_threshold', 10))
299
300
301 greyscale = kwargs.pop("greyscale", False)
302 if greyscale and not kwargs.has_key("cmap"):
303 kwargs["cmap"] =\
304 pylab.matplotlib.colors.LinearSegmentedColormap("clrs",\
305 pylab.matplotlib.cm.hot_r._segmentdata)
306 elif not kwargs.has_key("cmap"):
307 kwargs["cmap"] = jet
308
309
310
311
312
313 tablenames = map(plotutils.display_name, tablenames)
314
315 if len(columns) == 2:
316 plot = plotutils.ScatterPlot(xlabel, ylabel,title, subtitle)
317 for i in range(len(tablenames)):
318 plot.add_content(data[i][0], data[i][1], label=tablenames[i],\
319 **kwargs)
320 plot.finalize(loc=loc)
321 if hidden_colorbar:
322 plotutils.add_colorbar(plot.ax, visible=False)
323 else:
324 if dqstyle:
325 plot = plotutils.DQScatterPlot(xlabel, ylabel, colorlabel,\
326 title, subtitle)
327 else:
328 plot = plotutils.ColorbarScatterPlot(xlabel, ylabel, colorlabel,\
329 title, subtitle)
330 plot.add_content(data[0][0], data[0][1], data[0][2],\
331 label=tablenames[0], **kwargs)
332 kwargs.pop("cmap", None)
333 if dqstyle:
334 plot.finalize(logcolor=logcolor, clim=colorlim, loc=loc,\
335 threshold=dqthresh)
336 else:
337 plot.finalize(logcolor=logcolor, clim=colorlim, loc=loc)
338
339
340 if loudest:
341 plot.ax.plot([loudest[0]], [loudest[1]], marker="*", color="gold",\
342 markersize=15)
343
344
345 if logx:
346 plot.ax.set_xscale("log")
347 if logy:
348 plot.ax.set_yscale("log")
349 plot.ax.autoscale_view(tight=True, scalex=True, scaley=True)
350
351 if limits[0]:
352 plot.ax.set_xlim(limits[0])
353 if limits[1]:
354 plot.ax.set_ylim(limits[1])
355
356 plot.ax.grid(True, which="both")
357 if xcolumn == "time":
358 plotutils.set_time_ticks(plot.ax)
359 plotutils.set_minor_ticks(plot.ax)
360
361
362 if greyscale:
363 plot.ax.patch.set_facecolor("#E8E8E8")
364 plot.savefig(outfile, bbox_inches=bbox_inches,\
365 bbox_extra_artists=plot.ax.texts)
366 plot.close()
367
368
369
370
371
372 -def plothistogram(lsctable, outfile, column="snr", numbins=100,\
373 colorcolumn=None, colorbins=10, normalize=None,\
374 cumulative=None, bar=False, **kwargs):
375 """
376 Plot any column of a valid LigoLW table against any other, coloured by any
377 other, or plot all the same if you're that way inclined. Multiple tables
378 can be given for any non-coloured plot.
379
380 "time" as a column means the output of get_peak for Burst tables, get_end
381 for Inspiral tables, and get_start for ringdown tables
382
383 Arguments:
384
385 table : [ dict | glue.ligolw.table.Table ]
386 dict of ("name", table) pairs, or single LigoLW table
387 outfile : string
388 string path for output plot
389
390 Keyword arguments:
391 xcolumn : string
392 valid column of triggers table to plot on x-axis
393 ycolumn : string
394 valid column of triggers table to plot on y-axis
395 zcolumn : string
396 valid column of triggers table to use for colorbar (optional).
397 """
398
399 if isinstance(lsctable, table.Table):
400 tables = {"_":lsctable}
401 else:
402 tables = lsctable
403 tablenames = sorted(tables.keys())
404 tables = [tables[n] for n in tablenames]
405
406
407 xlim = kwargs.pop("xlim", None)
408 ylim = kwargs.pop("ylim", None)
409 zlim = kwargs.pop("zlim", None)
410 clim = kwargs.pop("clim", zlim)
411
412
413 logx = kwargs.pop("logx", False)
414 logy = kwargs.pop("logy", False)
415 logcolor = kwargs.pop("logcolor", False)
416
417
418 hidden_colorbar = kwargs.pop("hidden_colorbar", False)
419
420
421 bbox_inches = kwargs.pop("bbox_inches", "tight")
422
423
424 fill = kwargs.pop("fill", False)
425
426
427 rate = kwargs.pop("rate", False)
428 if normalize:
429 rate = True
430 if not hasattr(normalize, "__iter__"):
431 normalize = [normalize]*len(tablenames)
432 else:
433 normalize = [1]*len(tablenames)
434 loc = kwargs.pop("loc", 0)
435
436
437 xlabel = kwargs.pop("xlabel", plotutils.display_name(column))
438 if rate and cumulative:
439 ylabel = kwargs.pop("ylabel", "Cumulative rate (Hz)")
440 elif rate:
441 ylabel = kwargs.pop("ylabel", "Rate (Hz)")
442 elif not rate and cumulative:
443 ylabel = kwargs.pop("ylabel", "Cumulative number")
444 else:
445 ylabel = kwargs.pop("ylabel", "Number")
446
447 title = kwargs.pop("title", "")
448 subtitle = kwargs.pop("subtitle","")
449 colorlabel = kwargs.pop("colorlabel", colorcolumn\
450 and plotutils.display_name(colorcolumn) or "")
451
452
453
454
455
456 if colorcolumn:
457 colordata = get_column(tables[0], colorcolumn).astype(float)
458 if not clim and colordata.size != 0:
459 clim = [colordata.min(), colordata.max()]
460 elif not clim:
461 clim = [1, 10]
462 if isinstance(colorbins, int):
463 numcolorbins = colorbins
464 if logcolor:
465 colorbins = numpy.logspace(*numpy.log10(clim), num=numcolorbins)
466 else:
467 colorbins = numpy.linspace(*clim, num=numcolorbins)
468 else:
469 numcolorbins = len(colorbins)
470
471 data = list()
472 for i,lsctable in enumerate(tables):
473 d = get_column(lsctable, column).astype(float)
474 if i==0 and colorcolumn:
475 data.append(list())
476 for j in range(numcolorbins):
477 data[i].append(d[d>=colorbins[j]])
478 else:
479 data.append(d)
480
481
482
483
484
485 if colorcolumn:
486 raise NotImplementedError("This has not been implemented in "+\
487 "plotutils yet, here's your chance!")
488 plot = plotutils.ColorbarLineHistogram(xlabel, ylabel, colorlabel,\
489 title, subtitle)
490 for dataset,colorbin in zip(data[0], colorbins):
491 plot.add_content(dataset, normalize=normalize[0],\
492 colorvalue=colorbin, label=tablenames[0], **kwargs)
493 for i,name in enumerate(tablenames[1:]):
494 plot.add_content(data[i+1], normalize=normalize[i+1],\
495 label=name, **kwargs)
496 plot.finalize(logcolor=logcolor, colorlabel=colorlabel, loc=loc)
497 else:
498 color = kwargs.pop("color", None)
499 if isinstance(color, str):
500 color = color.split(',')
501 if len(color) == 1:
502 color = [color[0]]*len(serieslist)
503 plot = plotutils.LineHistogram(xlabel, ylabel, title, subtitle)
504 for i,name in enumerate(tablenames):
505 if color:
506 kwargs["color"] = color[i]
507 plot.add_content(data[i], normalize=normalize[i], label=name,\
508 **kwargs)
509 plot.finalize(loc=loc, logx=logx, logy=logy, bar=bar, num_bins=numbins,\
510 cumulative=cumulative)
511
512 plot.ax.autoscale_view(tight=True, scalex=True, scaley=True)
513
514 if xlim:
515 plot.ax.set_xlim(xlim)
516 if ylim:
517 plot.ax.set_ylim(ylim)
518
519
520 plot.savefig(outfile, bbox_inches=bbox_inches,\
521 bbox_extra_artists=plot.ax.texts)
522 plot.close()
523
524
525
526
527
528 -def plotrate(lsctable, outfile, stride=60, column="peak_frequency", bins=[],\
529 t0=None, **kwargs):
530 """
531 Plot rate versus time for the given ligolw table triggers, binned by the
532 given bincolumn using the bins list.
533
534 Arguments:
535
536 lsctable : glue.ligolw.table.Table
537 LIGOLW table containing a list of triggers
538 outfile : string
539 string path for output plot
540 """
541
542 timedata = get_column(lsctable, "time")
543 ratedata = get_column(lsctable, column)
544
545 if timedata.size:
546 timedata, ratedata = map(numpy.asarray,\
547 zip(*sorted(zip(timedata,ratedata),\
548 key=lambda (t,r): t)))
549
550
551
552 xlim = kwargs.pop("xlim", None)
553 if xlim:
554 start,end = xlim
555 elif timedata.size:
556 start = timedata.min()
557 end = timedata.max()
558 xlim = start, end
559 else:
560 start = 0
561 end = 1
562
563
564 logx = kwargs.pop("logx", False)
565 logy = kwargs.pop("logy", False)
566 ylim = kwargs.pop("ylim", False)
567 hidden_colorbar = kwargs.pop("hidden_colorbar", False)
568 bbox_inches = kwargs.pop("bbox_inches", "tight")
569 loc = kwargs.pop("loc", 0)
570
571
572 xbins = numpy.arange(float(start), float(end), stride)
573 if not bins:
574 bins = [[0, float("inf")]]
575 ybins = [map(float, bin_) for bin_ in bins]
576
577
578 xextent = (xbins[0], xbins[-1]+stride)
579 yextent = (ybins[0][0], ybins[-1][-1])
580 inanybin = (timedata > xextent[0]) & (timedata <= xextent[1])\
581 & (ratedata > yextent[0]) & (ratedata <= yextent[1])
582 timedata = timedata[inanybin]
583 ratedata = ratedata[inanybin]
584
585 ydata = numpy.zeros((len(xbins), len(ybins)))
586 for i,xbin in enumerate(xbins):
587
588 try:
589 idx = (timedata >= xbin+stride).nonzero()[0][0]
590
591 ratebin = numpy.asarray(sorted(ratedata[:idx]))
592 break_ = False
593
594 ratebin = numpy.asarray(sorted(ratedata[:idx]))
595 except IndexError:
596 idx = None
597 ratebin = numpy.asarray(sorted(ratedata))
598 break_ = True
599
600 for j,ybin in enumerate(ybins):
601
602 try:
603 yidx = (ratebin >= ybin[1]).nonzero()[0][0]
604 except IndexError:
605 ydata[i][j] = len(ratebin)/stride
606 break
607 ydata[i][j] = len(ratebin[:yidx])/stride
608 ratebin = ratebin[yidx:]
609 if break_: break
610 timedata = timedata[idx:]
611 xdata = xbins + stride/2
612
613
614 if t0:
615 xdata -= t0
616 start -= t0
617 end -= t0
618 xlabel = kwargs.pop("xlabel", None)
619 if not xlabel:
620 if not t0:
621 t0 = start
622 start -= t0
623 end -= t0
624 xdata -= t0
625 unit, timestr = plotutils.time_axis_unit(end-start)
626 start /= unit
627 end /= unit
628 xdata = xdata/unit
629 t0 = LIGOTimeGPS(t0)
630 if t0.nanoseconds==0:
631 xlabel = datetime.datetime(*date.XLALGPSToUTC(LIGOTimeGPS(t0))[:6])\
632 .strftime("%B %d %Y, %H:%M:%S %ZUTC")
633 else:
634 xlabel = datetime.datetime(*date.XLALGPSToUTC(\
635 LIGOTimeGPS(t0.seconds))[:6])\
636 .strftime("%B %d %Y, %H:%M:%S %ZUTC")
637 xlabel = "Time (%s) since %s (%s)"\
638 % (timestr,
639 xlabel.replace(" UTC", ".%.3s UTC" % t0.nanoseconds),\
640 t0)
641 ylabel = kwargs.pop("ylabel", "Rate (Hz)")
642 title = kwargs.pop("title", "")
643 subtitle = kwargs.pop("subtitle", "")
644
645
646
647
648
649 plot = plotutils.ScatterPlot(xlabel, ylabel, title, subtitle)
650
651
652 if logy:
653 ydata = numpy.ma.masked_where(ydata==0, ydata, copy=False)
654
655 for i,ybin in enumerate(ybins):
656 if list(ybin) == [0, numpy.inf]:
657 label = "_"
658 else:
659 label = "%s-%s" % (ybin[0], ybin[1])
660 plot.add_content(xdata, ydata[:,i], label=label, **kwargs)
661
662 plot.finalize(loc=loc)
663 if hidden_colorbar:
664 plotutils.add_colorbar(plot.ax, visible=False)
665 if logx:
666 plot.ax.set_xscale("log")
667 if logy:
668 plot.ax.set_yscale("log")
669 plot.ax.autoscale_view(tight=True, scalex=True, scaley=True)
670 plot.ax.set_xlim(start, end)
671 if ylim:
672 plot.ax.set_ylim(ylim)
673
674
675 plot.ax.grid(True, which="both")
676 plotutils.set_time_ticks(plot.ax)
677 plotutils.set_minor_ticks(plot.ax)
678
679
680 plot.savefig(outfile, bbox_inches=bbox_inches,\
681 bbox_extra_artists=plot.ax.texts)
682 plot.close()
683