1
2 """Contains functions to create various plots from lalapps_mvsc_player"""
3
4 __author__ = 'Tristan Miller <tmiller@caltech.edu>'
5 __prog__ = 'mvsc_plots'
6
7
8
9 import math,re
10 from matplotlib import pyplot
11 from matplotlib import cbook
12
13
14
15 -def patread(filename,station=None,readpat=False,readstr=False,\
16 headeronly=False,colsonly=False):
17 """Reads in a file from SprOutputWriterApp.
18
19 If readpat is false, reads a dat file. Otherwise, reads a pat file.
20 If headeronly is true, only outputs header in a list.
21 If cols only is true, only outputs header in dict.
22 If station is given, outputs two dicts, the first with more standardized
23 names (ie 'H1L1' => 't1t2')
24 If readstr is false, reads floats. Otherwise, reads strings (which keeps
25 full precision)."""
26
27 try:
28 f = open(filename)
29 except IOError:
30 print '***Error!*** Trouble opening file', filename
31 return
32
33 p = re.compile(r'\S+')
34
35 if readpat:
36 h = f.readline()
37
38 h = f.readline()
39 header = p.findall(h)
40
41 if headeronly:
42 f.close()
43 return header
44
45 cols = {}
46 for i in range(len(header)):
47 cols[header[i]] = i
48
49 if station:
50 iter = cols.iterkeys()
51 cols2 = {}
52
53 while True:
54 try:
55 key = iter.next()
56 origcopy = key
57
58 for j in range(0,6,2):
59 for i in range(0,len(station),2):
60 if key[j:j+2] == station[i:i+2]:
61 key = key[:j]+'t'+str(i/2+1)+key[j+2:]
62 break
63
64 cols2[key] = cols[origcopy]
65
66 except StopIteration:
67 break
68
69 if colsonly:
70 f.close()
71 if station:
72 return cols2,cols
73 else:
74 return cols
75
76
77
78 data = []
79 if readstr:
80 n = f.readline()
81 m = p.findall(n)
82 if m:
83 for i in range(len(m)):
84 data.append([m[i]])
85
86 while m:
87 n = f.readline()
88 m = p.findall(n)
89 if m:
90 for i in range(len(m)):
91 data[i].append(m[i])
92 else:
93 n = f.readline()
94 m = p.findall(n)
95 if m:
96 for i in range(len(m)):
97 data.append([float(m[i])])
98
99 while m:
100 n = f.readline()
101 m = p.findall(n)
102 if m:
103 for i in range(len(m)):
104 data[i].append(float(m[i]))
105
106 f.close()
107
108 print 'Finished reading file', filename
109
110 if station:
111 return data, cols2, cols
112 else:
113 return data,cols
114
115
116
118 """Reads the printout of SprBaggerDecisionTreeApp and returns info.
119
120 Info returned: List of FOM values
121 Table of number of tree splits
122 Number of timeslides in training set
123 Number of injections in training set
124 Number of timeslides in validation set
125 Number of injections in validation set"""
126
127 try:
128 f = open(filename)
129 except IOError:
130 print '***Error!*** Trouble opening file', filename
131 return
132
133 output = []
134 p1 = re.compile(r'Points in class [0(1):]{5} (\S+)')
135 counter = 0
136
137 while True:
138 n = f.readline()
139 m = p1.match(n)
140 if m:
141 output.append(int(m.group(1)))
142 counter += 1
143 elif counter >= 4:
144 break
145 elif not n:
146 print '***Error!*** Missing info in ', filename
147 return
148
149 p2 = re.compile(r'Validation FOM=(\S+)')
150 fom = []
151 flag = False
152 while True:
153 n = f.readline()
154 m = p2.match(n)
155 if m:
156 fom.append(float(m.group(1)))
157 flag = True
158 elif flag:
159 break
160 elif not n:
161 print '***Error!*** No FOM in ', filename
162 return
163
164 treesplits=[]
165 p3 = re.compile(r'Variable\s+(\S+)\s+Splits\s+(\S+)\s+Delta FOM\s+(\S+)')
166 while True:
167 n = f.readline()
168 m = p3.match(n)
169 if m:
170 treesplits.append([])
171 treesplits[-1].append(m.group(1))
172 treesplits[-1].append(int(m.group(2)))
173 treesplits[-1].append(float(m.group(3)))
174 elif not n:
175 break
176
177 f.close()
178 return fom,treesplits,output[0],output[1],output[2],output[3]
179
180
181
183 """Rewrites the output of SprOutputWriterApp because it couldn't be
184 bothered to save the numbers to full precision.
185
186 Returns the data (so that you don't have to reread the file later)"""
187
188 patdata,patcols = patread(patpath,readpat=True,readstr=True)
189
190 resultdata,resultcols = patread(resultpath)
191 header = patread(resultpath,headeronly=True)
192
193 try:
194 f = open(resultpath,'w')
195 except IOError:
196 print '***Error!*** Trouble opening file', filename
197 return
198
199 for i in range(len(header)):
200 f.write(header[i] + ' ')
201
202 baggercol = resultcols['Bagger']
203 for i in range(len(patdata[0])):
204 f.write('\n')
205
206 for j in range(3):
207 f.write(str(resultdata[j][i]) + ' ')
208
209 for j in range(len(patdata)-1):
210 f.write(patdata[j][i] + ' ')
211
212 f.write(str(resultdata[baggercol][i]))
213
214 f.close()
215
216 return resultdata
217
218
219
221 """Sorts data into injections and timeslides"""
222
223 injections = []
224 timeslides = []
225 for i in range(len(data)):
226 injections.append([])
227 timeslides.append([])
228
229 for i in range(len(data[0])):
230 if data[1][i] == 0:
231 for j in range(len(data)):
232 timeslides[j].append(data[j][i])
233 else:
234 for j in range(len(data)):
235 injections[j].append(data[j][i])
236
237 return injections,timeslides
238
239
240
241 -def ROC(timeslides,injections):
242 """Computes true and false positive rates.
243
244 May change the order of timeslides and injections."""
245
246 timeslides.sort()
247 injections.sort()
248 n = len(injections)
249 m = len(timeslides)
250
251 truepos = []
252 falsepos = []
253
254 for i in range(n):
255 if injections[i] != injections[i-1]:
256 for y in range(n):
257 if injections[y]>injections[i]:
258 truepos.append(1-(y-1)/float(n))
259 break
260 else:
261 truepos.append(0.0)
262
263
264
265
266 if timeslides[-1]<=injections[i]:
267 falsepos.append(1/float(m+1))
268 else:
269 for x in range(2,m+1):
270 if timeslides[-x]<=injections[i]:
271 falsepos.append((x-1)/float(m))
272 break
273 else:
274 falsepos.append(1.0)
275
276 return truepos, falsepos
277
278
279
281 """Calculates the Wilson interval, the confidence interval for a binomial
282 distribution.
283
284 Returns the appropriate upper and lower error bars.
285 The confidence level used is always 68%."""
286
287 n = float(n)
288 diff = math.sqrt(max(p*(1-p)/n + 0.25/n**2,0)) / (1+1/n)
289 mid = (-p/n + 0.5/n) / (1+1/n)
290 upper = mid+diff
291 lower = diff-mid
292
293 return upper,lower
294
295
296
298 """Transforms lists (x,y) into (x1,y1) such that plotting them will
299 create a stairs graph."""
300
301 x1 = []
302 y1 = [y[0]]
303 for i in range(len(x)-1):
304 x1.append(x[i])
305 x1.append(x[i])
306 y1.append(y[i+1])
307 y1.append(y[i+1])
308 x1.append(x[-1])
309
310 return x1,y1
311
312
313
314 -def top_events(data,cols,n,stations, mvsc_to_fan=None):
315 """Finds the n events with the highest MVSC values.
316
317 If multiple events tie for last, all will be included.
318 Must give data in string form to keep precision.
319 If given mvsc_to_fan (which is [list of mvsc cutoff values, list of
320 corresponding false alarm numbers]), then will add a column of FAN."""
321
322 baggercol = cols['Bagger']
323
324
325
326
327 data2 = []
328 for i in range(len(data[0])):
329 data2.append([])
330 for j in range(len(data)):
331 if j == baggercol:
332 data2[-1].append(float(data[j][i]))
333 else:
334 data2[-1].append( data[j][i] )
335
336 sorter = cbook.Sorter()
337 sorter(data2,baggercol)
338 mvsc_cutoff = data2[-min(len(data2),n)][baggercol]
339
340 events = []
341 index = []
342
343 for i in range(1,len(data2)+1):
344 if data2[-i][baggercol] >= mvsc_cutoff:
345 index.append(i)
346 events.append([])
347 for j in range(len(data2[0])):
348 events[-1].append(data2[-i][j])
349 else:
350 break
351
352 if mvsc_to_fan:
353
354 cols['FAN'] = len(cols)
355
356 for i in range(len(events)):
357 mvsc_cutoff = events[i][baggercol]
358
359 for j in range(1,len(mvsc_to_fan[0])+1):
360 if mvsc_to_fan[0][-j] <= mvsc_cutoff:
361 break
362
363 events[i].append(mvsc_to_fan[1][-j])
364
365 return events
366
367
368
369 -def ROCplot(data_inj,data_ts,cols,op_point = 1,ts_trig_ratio = 25):
370 """Creates an ROC plot from one file.
371
372 Returns the confidence interval of the resulting efficiency."""
373
374 timeslides = data_ts[cols['Bagger']][:]
375 injections = data_inj[cols['Bagger']][:]
376 timeslides_snr = []
377 injections_snr = []
378
379 effsnr1 = cols['t1get_effective_snr()']
380 effsnr2 = cols['t2get_effective_snr()']
381
382 for i in range(len(data_ts[0])):
383 timeslides_snr.append(data_ts[effsnr1][i]**2+ \
384 data_ts[effsnr2][i]**2)
385 for i in range(len(data_inj[0])):
386 injections_snr.append(data_inj[effsnr1][i]**2+ \
387 data_inj[effsnr2][i]**2)
388
389 truepos,falsepos = ROC(timeslides,injections)
390 truepos_snr,falsepos_snr = ROC(timeslides_snr,injections_snr)
391
392 binomerru = []
393 binomerrl = []
394 for i in range(len(truepos)):
395 upper,lower = wilson( truepos[i], len(injections) )
396 binomerru.append(upper)
397 binomerrl.append(lower)
398 binomerru.append(0)
399 binomerrl.append(0)
400 t = binomerru.pop()
401 t = binomerrl.pop()
402
403 binomerru_snr = []
404 binomerrl_snr = []
405 for i in range(len(truepos_snr)):
406 upper,lower = wilson( truepos_snr[i], len(injections) )
407 binomerru_snr.append(upper)
408 binomerrl_snr.append(lower)
409 binomerru_snr.append(0)
410 binomerrl_snr.append(0)
411 t = binomerru_snr.pop()
412 t = binomerrl_snr.pop()
413
414 falsepos2 = []
415 falsepos2_snr = []
416 for i in range(len(falsepos)):
417 falsepos2.append( falsepos[i]*float(len(timeslides))/ts_trig_ratio )
418 for i in range(len(falsepos_snr)):
419 falsepos2_snr.append( falsepos_snr[i]*float(len(timeslides))/ \
420 ts_trig_ratio )
421
422 falsepos1,truepos1 = stairs(falsepos2,truepos)
423 falsepos1_snr,truepos1_snr = stairs(falsepos2_snr,truepos_snr)
424
425 pyplot.figure()
426 pyplot.errorbar(falsepos1,truepos1,yerr=[binomerrl,binomerru] \
427 ,label='MVSC',color='green',linewidth=2)
428 pyplot.errorbar(falsepos1_snr,truepos1_snr,yerr=[binomerrl_snr, \
429 binomerru_snr], label='Effective SNR', \
430 color='blue',linewidth=2)
431 pyplot.semilogx()
432 pyplot.xlabel('False alarm number')
433 pyplot.ylabel('True positive rate')
434 pyplot.title('ROC curve')
435
436 op_point = max(op_point,1/(float(len(timeslides))))
437 pyplot.plot([op_point,op_point],[.85,1.01],'r',label='Operating point')
438
439 for i in range(len(falsepos)):
440 if falsepos2[i] <= op_point:
441 break
442
443 xmin,xmax = pyplot.xlim()
444 mid = truepos[i]
445 upper = mid + binomerru[2*i]
446 lower = mid - binomerrl[2*i]
447 pyplot.plot([xmin,xmax],[mid,mid],'r')
448 pyplot.plot([xmin,xmax],[lower,lower],'r:')
449 pyplot.plot([xmin,xmax],[upper,upper],'r:')
450
451 pyplot.xlim(xmin,xmax)
452
453 pyplot.ylim(.85,1.01)
454 pyplot.legend(loc='lower right')
455
456 return lower,upper
457
458
459
461 """Plots FOM vs bagger cycle."""
462
463 pyplot.figure()
464 pyplot.plot(fom)
465 pyplot.xlabel('Bagger cycle')
466 pyplot.ylabel('Figure of Merit')
467 pyplot.title('Figure of Merit vs number of Bagger cycles')
468 pyplot.xlim(-len(fom)*0.05,len(fom))
469
470
471
472 -def mvsc_vs_effsnr(data_inj,data_ts,cols, mvsc_cutoff=None,effsnr_cutoff=None,\
473 zerodata=None):
474 """Plots mvsc values vs the sum of the squares of effective snr."""
475
476 timeslides = data_ts[cols['Bagger']][:]
477 injections = data_inj[cols['Bagger']][:]
478 timeslides_snr = []
479 injections_snr = []
480
481 effsnr1 = cols['t1get_effective_snr()']
482 effsnr2 = cols['t2get_effective_snr()']
483
484 for i in range(len(data_ts[0])):
485 timeslides_snr.append(data_ts[effsnr1][i]**2+ \
486 data_ts[effsnr2][i]**2)
487 for i in range(len(data_inj[0])):
488 injections_snr.append(data_inj[effsnr1][i]**2+ \
489 data_inj[effsnr2][i]**2)
490
491 pyplot.figure()
492 pyplot.semilogy(timeslides,timeslides_snr,'kx',label='Timeslides')
493 pyplot.plot(injections,injections_snr,'r+',mec="red",label='Injections')
494 pyplot.xlabel('MVSC')
495 pyplot.ylabel('Combined effective SNR squared')
496 pyplot.title('MVSC vs Effective SNR')
497
498 if zerodata:
499 zerolag_snr = []
500 for i in range(len(zerodata[0])):
501 zerolag_snr.append(zerodata[cols['t1get_effective_snr()']][i]**2+ \
502 zerodata[cols['t2get_effective_snr()']][i]**2 )
503
504 pyplot.plot(zerodata[cols['Bagger']],zerolag_snr,'*g',mec='g',\
505 label='Zero Lag' )
506
507 if mvsc_cutoff:
508 ymin,ymax = pyplot.ylim()
509 pyplot.plot([mvsc_cutoff,mvsc_cutoff],[ymin,ymax],'b',\
510 label='Operating point cutoff')
511 if effsnr_cutoff:
512 pyplot.plot([-.1,1.1],[effsnr_cutoff,effsnr_cutoff],'b')
513
514 pyplot.xlim(-.1,1.1)
515 pyplot.legend(loc='upper left')
516
517
518
519 -def fraction_detected(data_inj,data_ts, cols, afar = None,mvsc_cutoff=None, \
520 distance=False ):
521 """Graphs the fraction detected vs snr or distance.
522
523 Returns the false alarm rate (afar), the mvsc cutoff and the effsnr
524 cutoff. Graphs vs snr by default, but if distance=True, will graph
525 against distance instead. Requires either the afar or mvsc_cutoff
526 in order to choose operating point."""
527
528
529
530 if distance:
531 minbin = 1.
532 stepfactor = 3
533 maxbin = 1000000000
534 else:
535
536 minbin = 60.
537 stepfactor = 1.2
538 maxbin = 1000
539
540 snrbins = [ minbin ]
541 i = 1
542 while snrbins[-1] <= maxbin:
543 snrbins.append( stepfactor**(i) * snrbins[0] )
544 i += 1
545
546
547
548 inj = []
549 ts = []
550 baggercol = cols['Bagger']
551 effsnr1 = cols['t1get_effective_snr()']
552 effsnr2 = cols['t2get_effective_snr()']
553
554 if distance:
555 speccol1 = cols['t1eff_distance']
556 speccol2 = cols['t2eff_distance']
557 for i in range(len(data_inj[0])):
558 inj.append( ( data_inj[cols['Bagger']][i], \
559 data_inj[effsnr1][i]**2 + data_inj[effsnr2][i]**2, \
560 (data_inj[speccol1][i]+ \
561 data_inj[speccol2][i])**3/8 ) )
562 else:
563 speccol1 = cols['t1snr']
564 speccol2 = cols['t2snr']
565 for i in range(len(data_inj[0])):
566 inj.append( ( data_inj[cols['Bagger']][i], \
567 data_inj[effsnr1][i]**2 + data_inj[effsnr2][i]**2, \
568 data_inj[speccol1][i]**2+ \
569 data_inj[speccol2][i]**2 ) )
570
571 for i in range(len(data_ts[0])):
572 ts.append( ( data_ts[baggercol][i], \
573 data_ts[effsnr1][i]**2 + data_ts[effsnr2][i]**2 ) )
574
575
576 sorter = cbook.Sorter()
577 ts_mvsc = sorter(ts,0)
578 ts_effsnr = sorter(ts,1)
579
580 if afar:
581
582
583 cutoff = int(afar*len(ts)) + 1
584 mvsc_cutoff = ts_mvsc[-cutoff][0]
585 elif mvsc_cutoff:
586 for i in range(1,len(ts_mvsc)+1):
587 if ts_mvsc[-i][0] < mvsc_cutoff:
588 break
589
590 cutoff = i
591 afar = float(cutoff - 1)/len(ts)
592 else:
593 print '***Error*** Must give fraction_detected either an afar value '+\
594 'or a mvsc_cutoff value'
595
596 effsnr_cutoff = ts_effsnr[-cutoff][1]
597
598
599 det_mvsc = []
600 notdet_mvsc = []
601 det_effsnr = []
602 notdet_effsnr = []
603 for i in range(len(inj)):
604 if inj[i][0] >= mvsc_cutoff:
605 det_mvsc.append(inj[i][2])
606 else:
607 notdet_mvsc.append(inj[i][2])
608
609 if inj[i][1] >= effsnr_cutoff:
610 det_effsnr.append(inj[i][2])
611 else:
612 notdet_effsnr.append(inj[i][2])
613
614
615 mvsc_frac = []
616 effsnr_frac = []
617 mvsc_erru = []
618 mvsc_errl = []
619 effsnr_erru = []
620 effsnr_errl = []
621 for i in range(len(snrbins)-1):
622 det_mvsc_total = 0
623 mvsc_total = 0
624
625 for j in range(len(det_mvsc)):
626 if (det_mvsc[j] > snrbins[i]) & (det_mvsc[j] <= snrbins[i+1]):
627 det_mvsc_total += 1
628 for j in range(len(notdet_mvsc)):
629 if (notdet_mvsc[j] > snrbins[i]) & \
630 (notdet_mvsc[j] <= snrbins[i+1]):
631 mvsc_total += 1
632
633 mvsc_total += det_mvsc_total
634
635 det_effsnr_total = 0
636 effsnr_total = 0
637
638 for j in range(len(det_effsnr)):
639 if (det_effsnr[j] > snrbins[i]) & (det_effsnr[j] <= snrbins[i+1]):
640 det_effsnr_total += 1
641 for j in range(len(notdet_effsnr)):
642 if (notdet_effsnr[j] > snrbins[i]) & \
643 (notdet_effsnr[j] <= snrbins[i+1]):
644 effsnr_total += 1
645
646 effsnr_total += det_effsnr_total
647
648
649 try:
650 mvsc_frac.append(float(det_mvsc_total)/mvsc_total)
651 upper,lower = wilson(float(det_mvsc_total)/mvsc_total, \
652 mvsc_total )
653 mvsc_erru.append(upper)
654 mvsc_errl.append(lower)
655 except ZeroDivisionError:
656 mvsc_frac.append(0.5)
657 mvsc_erru.append(0.5)
658 mvsc_errl.append(0.5)
659 try:
660 effsnr_frac.append(float(det_effsnr_total)/effsnr_total)
661 upper,lower = wilson( float(det_effsnr_total)/effsnr_total, \
662 effsnr_total )
663 effsnr_erru.append(upper)
664 effsnr_errl.append(lower)
665 except ZeroDivisionError:
666 effsnr_frac.append(0.5)
667 effsnr_erru.append(0.5)
668 effsnr_errl.append(0.5)
669
670 pyplot.figure()
671 pyplot.semilogx()
672 pyplot.errorbar(snrbins[:-1],effsnr_frac,marker='s',linestyle='-', \
673 yerr=[effsnr_errl,effsnr_erru],label='Effective SNR',mec='blue')
674 pyplot.errorbar(snrbins[:-1],mvsc_frac,marker='*',linestyle='-', \
675 yerr=[mvsc_errl,mvsc_erru],label='MVSC',mec='green')
676 pyplot.xlim(snrbins[0]/stepfactor**2,snrbins[-1]*stepfactor**2)
677 pyplot.ylim(-.1,1.1)
678
679 if distance:
680 pyplot.xlabel('Effective Volume')
681 pyplot.legend(loc='lower left')
682 else:
683 pyplot.xlabel('Combined SNR squared')
684 pyplot.legend(loc='lower right')
685
686 pyplot.ylabel('Fraction detected')
687 pyplot.title('Fraction of injections detected using different measures')
688
689 return afar,mvsc_cutoff,effsnr_cutoff
690
691
692
693 -def snr_vs_chisqr(data_inj,data_ts,cols,afar = 1.0/2000,zerodata=None ):
694 """Plots SNR vs Chi Squared, indicating which triggers were correctly
695 classified, and which were incorrectly classified.
696
697 afar is the allowed false alarm rate in the classification"""
698
699
700 inj = []
701 ts = []
702 baggercol = cols['Bagger']
703 chisq1 = cols['t1chisq']
704 chisq2 = cols['t2chisq']
705 snr1 = cols['t1snr']
706 snr2 = cols['t2snr']
707
708 for i in range(len(data_ts[0])):
709 ts.append( ( data_ts[baggercol][i], \
710 data_ts[chisq1][i]**2+ \
711 data_ts[chisq2][i]**2, \
712 data_ts[snr1][i]**2+ \
713 data_ts[snr2][i]**2 ) )
714 for i in range(len(data_inj[0])):
715 inj.append( ( data_inj[baggercol][i], \
716 data_inj[chisq1][i]**2+ \
717 data_inj[chisq2][i]**2, \
718 data_inj[snr1][i]**2+ \
719 data_inj[snr2][i]**2 ) )
720
721 cutoff = int(afar*len(ts)) + 1
722
723
724 sorter = cbook.Sorter()
725 ts_mvsc = sorter(ts,0)
726 mvsc_cutoff = ts_mvsc[-cutoff][0]
727
728
729 falsepos_chi = []
730 falsepos_snr = []
731 trueneg_chi = []
732 trueneg_snr = []
733 truepos_chi = []
734 truepos_snr = []
735 falseneg_chi = []
736 falseneg_snr = []
737 for i in range(len(ts)):
738 if ts[i][0] >= mvsc_cutoff:
739 falsepos_chi.append(ts[i][1])
740 falsepos_snr.append(ts[i][2])
741 else:
742 trueneg_chi.append(ts[i][1])
743 trueneg_snr.append(ts[i][2])
744 for i in range(len(inj)):
745 if inj[i][0] >= mvsc_cutoff:
746 truepos_chi.append(inj[i][1])
747 truepos_snr.append(inj[i][2])
748 else:
749 falseneg_chi.append(inj[i][1])
750 falseneg_snr.append(inj[i][2])
751
752 pyplot.figure()
753 pyplot.loglog(trueneg_snr,trueneg_chi,'xk',label='Filtered timeslides')
754 pyplot.plot(truepos_snr,truepos_chi,'+r',mec='r', \
755 label='Detected injections')
756 pyplot.plot(falseneg_snr,falseneg_chi,'oc',mec='c', \
757 label='Missed injections')
758 pyplot.plot(falsepos_snr,falsepos_chi,'sm',mec='m',label='False alarms')
759
760 if zerodata:
761 zero = []
762 for i in range(len(zerodata[0])):
763 zero.append( ( zerodata[cols['Bagger']][i],
764 zerodata[cols['t1chisq']][i]**2+ \
765 zerodata[cols['t2chisq']][i]**2, \
766 zerodata[cols['t1snr']][i]**2+ \
767 zerodata[cols['t2snr']][i]**2 ) )
768
769 zeropos_chi = []
770 zeropos_snr = []
771 zeroneg_chi = []
772 zeroneg_snr = []
773 for i in range(len(zero)):
774 if zero[i][0] >= mvsc_cutoff:
775 zeropos_chi.append(zero[i][1])
776 zeropos_snr.append(zero[i][2])
777 else:
778 zeroneg_chi.append(zero[i][1])
779 zeroneg_snr.append(zero[i][2])
780
781 pyplot.plot(zeroneg_snr,zeroneg_chi,'xb', \
782 label='Zero lag (noise)')
783 pyplot.plot(zeropos_snr,zeropos_chi,'*g',mec='g', \
784 label='Zero lag (signal)')
785
786 pyplot.legend(loc='upper left')
787 pyplot.ylabel('Combined Chisq squared')
788 pyplot.xlabel('Combined SNR squared')
789 pyplot.title('SNR vs Chi Squared')
790
791
792
793 -def FARplot(data_ts,cols,zerodata=None,ts_trig_ratio = 25):
794 """Graphs the cumulative number of detections vs MVSC threshold.
795
796 Returns the coordinates of plotted points so that the FAN can be extracted
797 from the MVSC value later."""
798
799 ts_mvsc = []
800
801 ts_mvsc = data_ts[cols['Bagger']][:]
802 ts_mvsc.sort()
803 n = len(ts_mvsc)
804
805 far = []
806 mvsc = []
807 mvsc_cutoff = 1
808 flag = True
809
810 for i in range(n):
811 if ts_mvsc[i] != ts_mvsc[i-1]:
812 far.append(float(n-i)/ts_trig_ratio)
813 mvsc.append(ts_mvsc[i])
814
815 if flag & (far[-1] < 1):
816 mvsc_cutoff = mvsc[-1]
817 flag = False
818
819 far2,mvsc2 = stairs(far,mvsc)
820
821 binomerru = []
822 binomerrl = []
823 for i in range(len(far)):
824 upper,lower = wilson( (far[i]/n)*ts_trig_ratio, n )
825 binomerru.append((upper*n)/ts_trig_ratio)
826 binomerrl.append((lower*n)/ts_trig_ratio)
827 binomerru.append(0)
828 binomerrl.append(0)
829 t = binomerru.pop()
830 t = binomerrl.pop()
831
832 pyplot.figure()
833 pyplot.errorbar(mvsc2,far2,yerr=[binomerrl,binomerru], \
834 label='Expected FAN (based on timeslides)')
835 pyplot.loglog()
836
837 if zerodata:
838 zero_mvsc = zerodata[cols['Bagger']][:]
839 zero_mvsc.sort()
840 m = len(zero_mvsc)
841
842 zerorate = []
843 zero_mvsc1 = []
844
845 for i in range(m):
846 if zero_mvsc[i] != zero_mvsc[i-1]:
847 zerorate.append(float(m-i))
848 zero_mvsc1.append(zero_mvsc[i])
849
850 zerorate2,zero_mvsc2 = stairs(zerorate,zero_mvsc1)
851
852 zbinomerru = []
853 zbinomerrl = []
854 for i in range(len(zerorate)):
855 upper,lower = wilson( zerorate[i]/m, m )
856 zbinomerru.append(upper*m)
857 zbinomerrl.append(lower*m)
858 zbinomerru.append(0)
859 zbinomerrl.append(0)
860 t = zbinomerru.pop()
861 t = zbinomerrl.pop()
862
863 pyplot.errorbar(zero_mvsc2,zerorate2,yerr=[zbinomerrl,zbinomerru],\
864 label='cumulative number of zero lags')
865
866 pyplot.legend(loc='lower left')
867
868 xmin,xmax = pyplot.xlim()
869 ymin,ymax = pyplot.ylim()
870 pyplot.plot([xmin,xmax],[1,1],'r')
871 pyplot.plot([mvsc_cutoff,mvsc_cutoff],[ymin,ymax],'r')
872 pyplot.xlim(xmin,xmax)
873 pyplot.ylim(ymin,ymax)
874
875 pyplot.xlabel('MVSC value cutoff')
876 pyplot.title('FAR plot')
877
878 return mvsc_cutoff,[mvsc,far]
879
880
881
882 -def IFANplot(data_ts,cols,zerodata,ts_trig_ratio=25):
883 """Plots the inverse false alarm number vs cumulative number of detections.
884
885 Returns false if fails, true if succeeds"""
886
887 ts_mvsc = data_ts[cols['Bagger']][:]
888 ts_mvsc.sort()
889 n = len(ts_mvsc)
890
891 zero_mvsc = zerodata[cols['Bagger']][:]
892 zero_mvsc.sort()
893 m = len(zero_mvsc)
894
895 cumnumber = []
896 cumnumber_erru = []
897 cumnumber_errl = []
898 ifan = []
899 ifan_erru = []
900 ifan_errl = []
901 for i in range(m):
902 if zero_mvsc[i] != zero_mvsc[i-1]:
903 cumnumber.append(float(m-i))
904 upper,lower = wilson(cumnumber[-1]/m,m)
905 cumnumber_erru.append(upper*m)
906 cumnumber_errl.append(lower*m)
907
908 for j in range(n):
909 if ts_mvsc[j] >= zero_mvsc[i]:
910 break
911
912 ifan.append(ts_trig_ratio/float(n-j))
913 upper1,lower1 = wilson(float(n-j)/n,n)
914 upper2 = lower1*n/(float(n-j)*(float(n-j)-lower1*n))
915 lower2 = upper1*n/(float(n-j)*(float(n-j)+upper1*n))
916 ifan_erru.append(upper2*ts_trig_ratio)
917 ifan_errl.append(lower2*ts_trig_ratio)
918
919 if len(cumnumber_errl) == 0:
920 return False
921
922 pyplot.figure()
923 pyplot.errorbar(ifan,cumnumber,xerr=[ifan_errl,ifan_erru], \
924 yerr=[cumnumber_errl,cumnumber_erru])
925 pyplot.loglog()
926
927 ymin,ymax = pyplot.ylim()
928 xmin,xmax = pyplot.xlim()
929 pyplot.plot([xmin,xmax],[1./xmin,1./xmax],'r',alpha=0.5)
930 pyplot.ylim(ymin,ymax)
931 pyplot.xlim(xmin,xmax)
932
933 pyplot.xlabel('Inverse False Alarm Number')
934 pyplot.ylabel('Cumulative Number')
935 pyplot.title('IFAN plot')
936
937 return True
938