1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26 """
27 Collection of functions to compute the efficiency and effective 4-volume
28 """
29
30 import sqlite3
31 import math
32 import pdb
33 from operator import itemgetter
34 import numpy
35 from scipy import special
36
37 from glue import segments
38 from glue.ligolw import table
39 from glue.ligolw import lsctables
40 from glue.ligolw import dbtables
41
42 from pylal import antenna
43 from pylal import ligolw_sqlutils as sqlutils
44 from pylal import ligolw_cbc_compute_durations as compute_dur
45
46
47
48
49
50
51
52
53
54
55
57 mchirp_DNS = (1.4+1.4) * (1./4)**(3.0/5.0)
58
59 return distance * (mchirp_DNS/mchirp)**(5.0/6.0)
60
61 -def decisive_dist(
62 h_dist, l_dist, v_dist,
63 mchirp, weight_dist, ifos):
64
65 dist_list = []
66 if 'H1' in ifos or 'H2' in ifos:
67 dist_list.append(h_dist)
68 if 'L1' in ifos:
69 dist_list.append(l_dist)
70 if 'V1' in ifos:
71 dist_list.append(v_dist)
72
73 if weight_dist:
74 return chirp_dist(sorted(dist_list)[1], mchirp)
75 else:
76 return sorted(dist_list)[1]
77
79 time = end_time + 1e-9*end_time_ns
80 return time
81
83 sqlquery = """
84 SELECT duration
85 FROM experiment_summary
86 JOIN experiment ON (
87 experiment_summary.experiment_id == experiment.experiment_id)
88 WHERE
89 datatype = :0
90 AND veto_def_name = :1
91 AND instruments = :2 """
92
93
94 total_dur = numpy.sum(connection.execute(sqlquery, (datatype, veto_cat, on_ifos)).fetchall() )
95
96 return total_dur
97
98
99
100
101
102
103
104
105
107
108 if dist_scale == "linear":
109 dist_bin_edges = numpy.arange(dist_bounds[0]-step, dist_bounds[1]+step, step)
110 elif dist_scale == "log":
111 log_limits = numpy.log10([dist_bounds[0], dist_bounds[1]])/numpy.log10(step)
112 dist_bin_edges = numpy.power(
113 step,
114 numpy.arange(log_limits[0]-1, log_limits[1]+1)
115 )
116
117 return dist_bin_edges
118
119
120 -def successful_injections(
121 connection,
122 tag,
123 on_ifos,
124 veto_cat,
125 dist_type = "distance",
126 weight_dist = False,
127 verbose = False):
128
129 """
130 My attempt to get a list of the simulations that actually made
131 it into some level of coincident time
132 """
133
134 xmldoc = dbtables.get_xml(connection)
135 connection.create_function('end_time_with_ns', 2, end_time_with_ns)
136
137
138 veto_segments = compute_dur.get_veto_segments(xmldoc, verbose)
139
140
141 sql_params_dict = {}
142 sqlquery = """
143 SELECT DISTINCT
144 simulation_id,
145 end_time_with_ns(geocent_end_time, geocent_end_time_ns),"""
146
147 if dist_type == "distance":
148 connection.create_function('distance_func', 2, chirp_dist)
149 sqlquery += """
150 distance_func(distance, sim_inspiral.mchirp)
151 FROM sim_inspiral """
152 elif dist_type == "decisive_distance":
153 connection.create_function('decisive_dist_func', 6, decisive_dist)
154 sql_params_dict['ifos'] = on_ifos
155 sql_params_dict['weight_dist'] = weight_dist
156 sqlquery += """
157 decisive_dist_func(
158 eff_dist_h, eff_dist_l, eff_dist_v,
159 sim_inspiral.mchirp, :weight_dist, :ifos)
160 FROM sim_inspiral """
161
162 if tag != 'ALL_INJ':
163
164 sqlquery += """
165 JOIN process_params ON (
166 process_params.process_id == sim_inspiral.process_id)
167 WHERE process_params.value = :usertag) """
168 sql_params_dict["usertag"] = tag
169 else:
170
171 tag = 'FULL_DATA'
172
173
174 ifo_segments = compute_dur.get_single_ifo_segments(
175 connection,
176 program_name = "inspiral",
177 usertag = tag)
178
179 zero_lag_dict = dict([(ifo, 0.0) for ifo in ifo_segments])
180
181 successful_inj = []
182
183 coinc_segs = compute_dur.get_coinc_segments(
184 ifo_segments - veto_segments[veto_cat],
185 zero_lag_dict)
186
187
188 for injection in connection.execute(sqlquery, sql_params_dict):
189 inj_segment = segments.segment(injection[1], injection[1])
190 if coinc_segs[on_ifos].intersects_segment( inj_segment ):
191 successful_inj.append( injection )
192
193 return successful_inj
194
195
196 -def found_injections(
197 connection,
198 tag,
199 on_ifos,
200 dist_type = "distance",
201 weight_dist = False,
202 verbose = False):
203
204 connection.create_function('end_time_with_ns', 2, end_time_with_ns)
205
206 sql_params_dict = {'ifos': on_ifos}
207 sqlquery = """
208 SELECT DISTINCT
209 sim_inspiral.simulation_id,
210 end_time_with_ns(geocent_end_time, geocent_end_time_ns), """
211
212 if dist_type == "distance":
213 connection.create_function('distance_func', 2, chirp_dist)
214 sqlquery += """
215 distance_func(distance, sim_inspiral.mchirp), """
216 elif dist_type == "decisive_distance":
217 connection.create_function('decisive_dist_func', 6, decisive_dist)
218 sql_params_dict['weight_dist'] = weight_dist
219 sqlquery += """
220 decisive_dist_func(
221 eff_dist_h, eff_dist_l, eff_dist_v,
222 sim_inspiral.mchirp, :weight_dist, :ifos), """
223
224 sqlquery += """
225 false_alarm_rate,
226 coinc_inspiral.snr
227 FROM
228 coinc_event_map AS coincs
229 JOIN coinc_event_map AS sims, coinc_inspiral, coinc_event, sim_inspiral ON (
230 coincs.coinc_event_id == sims.coinc_event_id
231 AND coinc_event.coinc_event_id == coincs.event_id
232 AND coinc_inspiral.coinc_event_id == coincs.event_id
233 AND sim_inspiral.simulation_id == sims.event_id)
234 JOIN process_params ON (
235 process_params.process_id == sim_inspiral.process_id)
236 WHERE
237 coincs.table_name = "coinc_event"
238 AND sims.table_name = "sim_inspiral"
239 AND coinc_event.instruments = :ifos """
240
241 if tag != 'ALL_INJ':
242 sqlquery += """
243 AND process_params.value = :usertag
244 """
245 sql_params_dict["tag"] = tag
246
247 injections = set(connection.execute(sqlquery, sql_params_dict).fetchall())
248
249
250 sqlquery = """
251 SELECT
252 end_time_with_ns(end_time, end_time_ns) AS trig_time,
253 snr AS trig_snr
254 FROM coinc_inspiral AS ci
255 JOIN experiment_map AS em, experiment_summary AS es ON (
256 ci.coinc_event_id == em.coinc_event_id
257 AND em.experiment_summ_id == es.experiment_summ_id )
258 WHERE es.datatype == 'all_data';
259 """
260 foreground = connection.executescript(sqlquery).fetchall()
261
262
263 inj_snr = numpy.array([inj[3] for inj in injections])
264 inj_time = numpy.array([inj[1] for inj in injections])
265 idx2remove = []
266 for time, snr in foreground:
267 indices = numpy.where(numpy.abs(inj_time - time) < 1.0)
268 if len(indices[0]):
269 idx = numpy.where(inj_snr[indices]/snr < 1.25)
270 if len(idx[0]):
271 idx2remove += list(indices[idx[0]])
272 for i in sorted(idx2remove, reverse=True):
273 del injections[i]
274
275
276 injections = sorted(injections, key=itemgetter(3), reverse=True)
277 found_inj = [inj[0:3] for inj in injections]
278 inj_fars = [inj[3] for inj in injections]
279 inj_snrs = [inj[4] for inj in injections]
280
281 return found_inj, inj_fars, inj_snrs
282
283
285 """
286 Calculate the optimal Bayesian credible interval for p(eff|k,n)
287 Posterior generated with binomial p(k|eff,n) and a uniform p(eff)
288 is the beta function: Beta(eff|k+1,n-k+1) where n is the number
289 of injected signals and k is the number of found signals.
290 """
291 eff_low = numpy.zeros(len(N))
292 eff_high = numpy.zeros(len(N))
293 for i, n in enumerate(N):
294 if n!= 0:
295
296 eff_cdf = special.betainc(K[i]+1, n-K[i]+1, eff_bin_edges)
297 pmf = ( eff_cdf[1:] - eff_cdf[:-1] )
298
299
300 a = numpy.argsort(pmf)[::-1]
301 sorted_cdf = numpy.cumsum(numpy.sort(pmf)[::-1])
302 j = numpy.argmin( numpy.abs(sorted_cdf - confidence) )
303
304 eff_low[i] = eff_bin_edges[:-1][ numpy.min(a[:(j+1)]) ]
305 eff_high[i] = eff_bin_edges[:-1][ numpy.max(a[:(j+1)]) ]
306
307 return eff_low, eff_high
308
309 -def detection_efficiency(
310 successful_inj,
311 found_inj,
312 found_fars,
313 far_list,
314 r,
315 confidence):
316 """
317 This function determines the peak efficiency for a given bin and associated
318 'highest density' confidence interval. The calculation is done for results
319 from each false-alarm-rate threshold
320 """
321
322 successful_inj = set(successful_inj) | set(found_inj)
323
324 successful_dist = [inj[2] for inj in successful_inj]
325 N, _ = numpy.histogram(successful_dist, bins = r)
326
327 significant_dist = [inj[2] for inj in found_inj]
328
329 eff_bin_edges = numpy.linspace(0, 1, 1e3+1)
330 eff = {
331 'mode': {},
332 'low': {},
333 'high': {}
334 }
335 for far in far_list:
336 for idx, coinc_far in enumerate(found_fars):
337 if coinc_far <= far:
338 new_start = idx
339 break
340
341 K, _ = numpy.histogram(significant_dist[new_start:], bins = r)
342 eff['mode'][far] = numpy.nan_to_num(numpy.float_(K)/N)
343
344
345 eff['low'][far], eff['high'][far] = binomial_confidence(K, N, eff_bin_edges, confidence)
346
347 return eff
348
349
350 -def rescale_dist(
351 on_ifos, dist_type, weight_dist,
352 phys_dist=None, param_dist=None
353 ):
354
355 N_signals = int(1e6)
356 trigTime = 0.0
357
358
359 if dist_type == 'decisive_distance':
360
361 ra = 360 * numpy.random.rand(N_signals)
362 dec = 180 * numpy.random.rand(N_signals) - 90
363
364 inclination = 180 * numpy.random.rand(N_signals)
365 polarization = 360 * numpy.random.rand(N_signals)
366
367 f_q = {}
368 for ifo in on_ifos:
369 f_q[ifo] = numpy.zeros(N_signals)
370 for index in range(N_signals):
371 _, _, _, f_q[ifo][index] = antenna.response(
372 trigTime,
373 ra[index], dec[index],
374 inclination[index], polarization[index],
375 'degree', ifo )
376
377 prob_d_d = {}
378 for j in range(len(phys_dist)-1):
379
380 volume = 4*numpy.pi/3 * numpy.random.uniform(
381 low = phys_dist[j]**3.0,
382 high = phys_dist[j+1]**3.0,
383 size = N_signals)
384 dist = numpy.power(volume*(3/(4*numpy.pi)), 1./3)
385
386
387 if dist_type == 'decisive_distance':
388 dist_eff = {}
389 for ifo in on_ifos:
390 dist_eff[ifo] = dist / f_q[ifo]
391 dist_dec = numpy.sort(dist_eff.values(), 0)[1]
392
393
394 if weight_dist:
395
396 mass1, mass2 = 0.13 * numpy.random.randn(2, N_signals) + 1.40
397 mchirp = numpy.power(mass1+mass2, -1./5) * numpy.power(mass1*mass2, 3./5)
398 if dist_type == 'decisive_distance':
399 dist_chirp = chirp_dist(dist_dec, mchirp)
400 if dist_type == 'distance':
401 dist_chirp = chirp_dist(dist, mchirp)
402 N_d, _ = numpy.histogram(dist_chirp, bins=param_dist)
403 else:
404 N_d, _ = numpy.histogram(dist_dec, bins=param_dist)
405
406 prob_d_d[phys_dist[j+1]] = numpy.float_(N_d)/numpy.sum(N_d)
407
408 return prob_d_d
409
411 """
412 This function creates a weighted average efficiency as a function of distance
413 by computing eff_wavg(D) = \sum_dc eff_mode(dc)p(dc|d). p(dc|d) is the probability
414 a signal has a parameterized distance dc if its physical distance is d.
415
416 The confidence interval for eff_wavg(d) is constructed from the quadrature sum
417 of the difference between the modes and the bounds, with each term again
418 weighted by p(dc|d).
419
420 This operation is done for each false-alarm-rate threshold.
421 """
422 eff_dist = {
423 'wavg': {},
424 'low': {},
425 'high': {}
426 }
427 for far, modes in measured_eff['mode'].items():
428 eff_dist['wavg'][far] = numpy.sum(modes * prob_dc_d.values(), 1)
429 eff_dist['low'][far] = numpy.sqrt(numpy.sum(
430 (measured_eff['low'][far] - modes)**2 * prob_dc_d.values(), 1)
431 )
432 eff_dist['high'][far] = numpy.sqrt(numpy.sum(
433 (measured_eff['high'][far] - modes)**2 * prob_dc_d.values(), 1)
434 )
435 return eff_dist
436
437
439 """
440 This function creates a weighted average efficiency within a given volume
441 by computing eff_wavg(D) = \sum_dc eff_mode(dc)p(dc|D). p(dc|D) is the
442 probability a signal has a parameterized distance dc if it falls within
443 physical distance D.
444
445 The confidence interval for eff_wavg(D) is constructed from the quadrature sum
446 of the difference between the modes and the bounds, with each term again
447 weighted by p(dc|D).
448
449 This operation is done for each false-alarm-rate threshold.
450 """
451
452 cumVol = numpy.cumsum(V_shell)
453 p_dc_d = numpy.array(prob_dc_d.values())
454 V_dc_D = numpy.cumsum((V_shell * p_dc_d.transpose()).transpose(), 0)
455 p_dc_D = (V_dc_D.transpose()/cumVol).transpose()
456
457 vol_eff = {
458 'wavg': {},
459 'low': {},
460 'high': {}
461 }
462 for far, modes in measured_eff['mode'].items():
463 vol_eff['wavg'][far] = numpy.sum(modes * p_dc_D, 1)
464 vol_eff['low'][far] = numpy.sqrt(numpy.sum(
465 (measured_eff['low'][far] - modes)**2 * p_dc_D, 1)
466 )
467 vol_eff['high'][far] = numpy.sqrt(numpy.sum(
468 (measured_eff['high'][far] - modes)**2 * p_dc_D, 1)
469 )
470
471 return vol_eff
472