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
28
29
30
31
32 """
33 This module contains classes and functions for post-processing the output
34 of the Bayesian parameter estimation codes.
35 """
36
37
38 import os
39 import sys
40 from math import cos,ceil,floor,sqrt,pi as pi_constant
41 import xml
42 from xml.dom import minidom
43 from operator import itemgetter
44
45
46 from lalinference.io import read_samples
47 import healpy as hp
48 import astropy.table
49 import lalinference.plot.cmap
50 import numpy as np
51 from numpy import fmod
52 import matplotlib
53 from matplotlib import pyplot as plt,cm as mpl_cm,lines as mpl_lines
54 from scipy import stats
55 from scipy import special
56 from scipy import signal
57 from scipy.optimize import newton
58 from scipy import interpolate
59 from numpy import linspace
60 import random
61 import socket
62 from itertools import combinations
63 from lalinference import LALInferenceHDF5PosteriorSamplesDatasetName as posterior_grp_name
64 import re
65
66 try:
67 import lalsimulation as lalsim
68 except ImportError:
69 print('Cannot import lalsimulation SWIG bindings')
70 raise
71
72 try:
73 from lalinference.imrtgr.nrutils import bbh_final_mass_non_spinning_Panetal, bbh_final_spin_non_spinning_Panetal, bbh_final_spin_non_precessing_Healyetal, bbh_final_mass_non_precessing_Healyetal, bbh_final_spin_projected_spin_Healyetal, bbh_final_mass_projected_spin_Healyetal, bbh_aligned_Lpeak_6mode_SHXJDK
74 except ImportError:
75 print('Cannot import lalinference.imrtgr.nrutils. Will suppress final parameter calculations.')
76
77 from matplotlib.ticker import FormatStrFormatter,ScalarFormatter,AutoMinorLocator
78
79 try:
80 hostname_short=socket.gethostbyaddr(socket.gethostname())[0].split('.',1)[1]
81 except:
82 hostname_short='Unknown'
83 if hostname_short=='ligo.caltech.edu' or hostname_short=='cluster.ldas.cit':
84 matplotlib.rcParams.update(
85 {'mathtext.fontset' : "custom",
86 'mathtext.fallback_to_cm' : True
87 })
88
89 from xml.etree.cElementTree import Element, SubElement, ElementTree, Comment, tostring, XMLParser
90
91
92 import pylal
93 import lal
94 from glue.ligolw import lsctables
95 from glue.ligolw import utils
96 from pylal import git_version
97
98 from _bayespputils import _skyhist_cart,_burnin
99
100 __author__="Ben Aylott <benjamin.aylott@ligo.org>, Ben Farr <bfarr@u.northwestern.edu>, Will M. Farr <will.farr@ligo.org>, John Veitch <john.veitch@ligo.org>, Vivien Raymond <vivien.raymond@ligo.org>"
101 __version__= "git id %s"%git_version.id
102 __date__= git_version.date
105 """Workaround for missing astropy.table.Table.replace_column method,
106 which was added in Astropy 1.1.
107
108 FIXME: remove this function when LALSuite depends on Astropy >= 1.1."""
109 index = table.colnames.index(old)
110 table.remove_column(old)
111 table.add_column(astropy.table.Column(new, name=old), index=index)
112
114 """Workaround for missing astropy.table.Table.as_array method,
115 which was added in Astropy 1.0.
116
117 FIXME: remove this function when LALSuite depends on Astropy >= 1.0."""
118 try:
119 return table.as_array()
120 except:
121 return table._data
122
123
124
125
126
127 logParams=['logl','loglh1','loglh2','logll1','loglv1','deltalogl','deltaloglh1','deltalogll1','deltaloglv1','logw','logprior','logpost','nulllogl','chain_log_evidence','chain_delta_log_evidence','chain_log_noise_evidence','chain_log_bayes_factor']
128
129 relativePhaseParams=[ a+b+'_relative_phase' for a,b in combinations(['h1','l1','v1'],2)]
130 snrParams=['snr','optimal_snr','matched_filter_snr'] + ['%s_optimal_snr'%(i) for i in ['h1','l1','v1']] + ['%s_cplx_snr_amp'%(i) for i in ['h1','l1','v1']] + ['%s_cplx_snr_arg'%(i) for i in ['h1', 'l1', 'v1']] + relativePhaseParams
131 calAmpParams=['calamp_%s'%(ifo) for ifo in ['h1','l1','v1']]
132 calPhaseParams=['calpha_%s'%(ifo) for ifo in ['h1','l1','v1']]
133 calParams = calAmpParams + calPhaseParams
134
135 massParams=['m1','m2','chirpmass','mchirp','mc','eta','q','massratio','asym_massratio','mtotal','mf']
136
137 spinParamsPrec=['a1','a2','phi1','theta1','phi2','theta2','costilt1','costilt2','costheta_jn','cosbeta','tilt1','tilt2','phi_jl','theta_jn','phi12','af']
138 spinParamsAli=['spin1','spin2','a1z','a2z']
139 spinParamsEff=['chi','effectivespin','chi_eff','chi_tot','chi_p']
140 spinParams=spinParamsPrec+spinParamsEff+spinParamsAli
141
142 cosmoParam=['m1_source','m2_source','mtotal_source','mc_source','redshift','mf_source']
143
144 ppEParams=['ppEalpha','ppElowera','ppEupperA','ppEbeta','ppElowerb','ppEupperB','alphaPPE','aPPE','betaPPE','bPPE']
145 tigerParams=['dchi%i'%(i) for i in range(8)] + ['dchi%il'%(i) for i in [5,6] ] + ['dxi%d'%(i+1) for i in range(6)] + ['dalpha%i'%(i+1) for i in range(5)] + ['dbeta%i'%(i+1) for i in range(3)] + ['dsigma%i'%(i+1) for i in range(4)]
146 bransDickeParams=['omegaBD','ScalarCharge1','ScalarCharge2']
147 massiveGravitonParams=['lambdaG']
148 tidalParams=['lambda1','lambda2','lam_tilde','dlam_tilde','lambdat','dlambdat']
149 energyParams=['e_rad', 'l_peak']
150 strongFieldParams=ppEParams+tigerParams+bransDickeParams+massiveGravitonParams+tidalParams+energyParams
151
152
153 distParams=['distance','distMPC','dist']
154 incParams=['iota','inclination','cosiota']
155 polParams=['psi','polarisation','polarization']
156 skyParams=['ra','rightascension','declination','dec']
157 phaseParams=['phase', 'phi0','phase_maxl']
158
159 timeParams=['time','time_mean']
160 endTimeParams=['l1_end_time','h1_end_time','v1_end_time']
161
162 statsParams=['logprior','logl','deltalogl','deltaloglh1','deltalogll1','deltaloglv1','deltaloglh2','deltaloglg1']
163 calibParams=['calpha_l1','calpha_h1','calpha_v1','calamp_l1','calamp_h1','calamp_v1']
164
165
166 confidenceLevels=[0.67,0.9,0.95,0.99]
167
168 greedyBinSizes={'mc':0.025,'m1':0.1,'m2':0.1,'mass1':0.1,'mass2':0.1,'mtotal':0.1,'mc_source':0.025,'m1_source':0.1,'m2_source':0.1,'mtotal_source':0.1,'eta':0.001,'q':0.01,'asym_massratio':0.01,'iota':0.01,'cosiota':0.02,'time':1e-4,'time_mean':1e-4,'distance':1.0,'dist':1.0,'redshift':0.01,'mchirp':0.025,'chirpmass':0.025,'spin1':0.04,'spin2':0.04,'a1z':0.04,'a2z':0.04,'a1':0.02,'a2':0.02,'phi1':0.05,'phi2':0.05,'theta1':0.05,'theta2':0.05,'ra':0.05,'dec':0.05,'chi':0.05,'chi_eff':0.05,'chi_tot':0.05,'chi_p':0.05,'costilt1':0.02,'costilt2':0.02,'thatas':0.05,'costheta_jn':0.02,'beta':0.05,'omega':0.05,'cosbeta':0.02,'ppealpha':1.0,'ppebeta':1.0,'ppelowera':0.01,'ppelowerb':0.01,'ppeuppera':0.01,'ppeupperb':0.01,'polarisation':0.04,'rightascension':0.05,'declination':0.05,'massratio':0.001,'inclination':0.01,'phase':0.05,'tilt1':0.05,'tilt2':0.05,'phi_jl':0.05,'theta_jn':0.05,'phi12':0.05,'flow':1.0,'phase_maxl':0.05,'calamp_l1':0.01,'calamp_h1':0.01,'calamp_v1':0.01,'calpha_h1':0.01,'calpha_l1':0.01,'calpha_v1':0.01,'logdistance':0.1,'psi':0.1,'costheta_jn':0.1,'mf':0.1,'mf_source':0.1,'af':0.02,'e_rad':0.1,'l_peak':0.1}
169 for s in snrParams:
170 greedyBinSizes[s]=0.05
171 for derived_time in ['h1_end_time','l1_end_time','v1_end_time','h1l1_delay','l1v1_delay','h1v1_delay']:
172 greedyBinSizes[derived_time]=greedyBinSizes['time']
173 for derived_phase in relativePhaseParams:
174 greedyBinSizes[derived_phase]=0.05
175 for param in tigerParams + bransDickeParams + massiveGravitonParams:
176 greedyBinSizes[param]=0.01
177 for param in tidalParams:
178 greedyBinSizes[param]=2.5
179
180 for loglname in statsParams:
181 greedyBinSizes[loglname]=0.1
182
183
184 __default_line_styles=['solid', 'dashed', 'dashdot', 'dotted']
185
186 __default_color_lst=['r','b','y','g','c','m']
187
188 __default_css_string="""
189 p,h1,h2,h3,h4,h5
190 {
191 font-family:"Trebuchet MS", Arial, Helvetica, sans-serif;
192 }
193
194 p
195 {
196 font-size:14px;
197 }
198
199 h1
200 {
201 font-size:20px;
202 }
203
204 h2
205 {
206 font-size:18px;
207 }
208
209 h3
210 {
211 font-size:16px;
212 }
213
214
215
216 table
217 {
218 font-family:"Trebuchet MS", Arial, Helvetica, sans-serif;
219 width:100%;
220 border-collapse:collapse;
221 }
222 td,th
223 {
224 font-size:12px;
225 border:1px solid #B5C1CF;
226 padding:3px 7px 2px 7px;
227 }
228 th
229 {
230 font-size:14px;
231 text-align:left;
232 padding-top:5px;
233 padding-bottom:4px;
234 background-color:#B3CEEF;
235 color:#ffffff;
236 }
237 #postable tr:hover
238 {
239 background: #DFF4FF;
240 }
241 #covtable tr:hover
242 {
243 background: #DFF4FF;
244 }
245 #statstable tr:hover
246 {
247 background: #DFF4FF;
248 }
249
250 img {
251 max-width: 510px;
252 max-height: 510px;
253 width:100%;
254 eight:100%;
255 }
256
257 .ppsection
258 {
259 border-bottom-style:double;
260 }
261
262 """
263 __default_javascript_string='''
264 //<![CDATA[
265 function toggle_visibility(tbid,lnkid)
266 {
267
268 if(document.all){document.getElementById(tbid).style.display = document.getElementById(tbid).style.display == 'block' ? 'none' : 'block';}
269
270 else{document.getElementById(tbid).style.display = document.getElementById(tbid).style.display == 'table' ? 'none' : 'table';}
271
272 document.getElementById(lnkid).value = document.getElementById(lnkid).value == '[-] Collapse' ? '[+] Expand' : '[-] Collapse';
273
274 }
275 //]]>
276
277 '''
278
279
280 from pylal.SimInspiralUtils import ExtractSimInspiralTableLIGOLWContentHandler
281 from glue.ligolw import lsctables
282 lsctables.use_in(ExtractSimInspiralTableLIGOLWContentHandler)
288 distributions={
289 'm1':'uniform',
290 'm2':'uniform',
291 'mc':None,
292 'eta':None,
293 'q':None,
294 'mtotal':'uniform',
295 'm1_source':None,
296 'm2_source':None,
297 'mtotal_source':None,
298 'mc_source':None,
299 'redshift':None,
300 'mf':None,
301 'mf_source':None,
302 'af':None,
303 'e_rad':None,
304 'l_peak':None,
305 'spin1':'uniform',
306 'spin2':'uniform',
307 'a1':'uniform',
308 'a2':'uniform',
309 'a1z':'uniform',
310 'a2z':'uniform',
311 'theta1':'uniform',
312 'theta2':'uniform',
313 'phi1':'uniform',
314 'phi2':'uniform',
315 'chi_eff':None,
316 'chi_tot':None,
317 'chi_p':None,
318 'tilt1':None,
319 'tilt2':None,
320 'costilt1':'uniform',
321 'costilt2':'uniform',
322 'iota':'np.cos',
323 'cosiota':'uniform',
324 'time':'uniform',
325 'time_mean':'uniform',
326 'dist':'x**2',
327 'ra':'uniform',
328 'dec':'np.cos',
329 'phase':'uniform',
330 'psi':'uniform',
331 'theta_jn':'np.sin',
332 'costheta_jn':'uniform',
333 'beta':None,
334 'cosbeta':None,
335 'phi_jl':None,
336 'phi12':None,
337 'logl':None,
338 'h1_end_time':None,
339 'l1_end_time':None,
340 'v1_end_time':None,
341 'h1l1_delay':None,
342 'h1v1_delay':None,
343 'l1v1_delay':None,
344 'lambdat' :None,
345 'dlambdat':None,
346 'lambda1' : 'uniform',
347 'lambda2': 'uniform',
348 'lam_tilde' : None,
349 'dlam_tilde': None,
350 'calamp_h1' : 'uniform',
351 'calamp_l1' : 'uniform',
352 'calpha_h1' : 'uniform',
353 'calpha_l1' : 'uniform',
354 'polar_eccentricity':'uniform',
355 'polar_angle':'uniform',
356 'alpha':'uniform'
357 }
358 try:
359 return distributions(name)
360 except:
361 return None
362
367 """
368 A lookup table for plot labels.
369 """
370 m1_names = ['mass1', 'm1']
371 m2_names = ['mass2', 'm2']
372 mc_names = ['mc','mchirp','chirpmass']
373 eta_names = ['eta','massratio','sym_massratio']
374 q_names = ['q','asym_massratio']
375 iota_names = ['iota','incl','inclination']
376 dist_names = ['dist','distance']
377 ra_names = ['rightascension','ra']
378 dec_names = ['declination','dec']
379 phase_names = ['phi_orb', 'phi', 'phase', 'phi0']
380 gr_test_names = ['dchi%d'%i for i in range(8)]+['dchil%d'%i for i in [5,6]]+['dxi%d'%(i+1) for i in range(6)]+['dalpha%d'%(i+1) for i in range(5)]+['dbeta%d'%(i+1) for i in range(3)]+['dsigma%d'%(i+1) for i in range(4)]
381
382 labels={
383 'm1':r'$m_1\,(\mathrm{M}_\odot)$',
384 'm2':r'$m_2\,(\mathrm{M}_\odot)$',
385 'mc':r'$\mathcal{M}\,(\mathrm{M}_\odot)$',
386 'eta':r'$\eta$',
387 'q':r'$q$',
388 'mtotal':r'$M_\mathrm{total}\,(\mathrm{M}_\odot)$',
389 'm1_source':r'$m_{1}^\mathrm{source}\,(\mathrm{M}_\odot)$',
390 'm2_source':r'$m_{2}^\mathrm{source}\,(\mathrm{M}_\odot)$',
391 'mtotal_source':r'$M_\mathrm{total}^\mathrm{source}\,(\mathrm{M}_\odot)$',
392 'mc_source':r'$\mathcal{M}^\mathrm{source}\,(\mathrm{M}_\odot)$',
393 'redshift':r'$z$',
394 'mf':r'$M_\mathrm{final}\,(\mathrm{M}_\odot)$',
395 'mf_source':r'$M_\mathrm{final}^\mathrm{source}\,(\mathrm{M}_\odot)$',
396 'af':r'$a_\mathrm{final}$',
397 'e_rad':r'$E_\mathrm{rad}\,(\mathrm{M}_\odot)$',
398 'l_peak':r'$L_\mathrm{peak}\,(10^{56}\,\mathrm{ergs}\,\mathrm{s}^{-1})$',
399 'spin1':r'$S_1$',
400 'spin2':r'$S_2$',
401 'a1':r'$a_1$',
402 'a2':r'$a_2$',
403 'a1z':r'$a_{1z}$',
404 'a2z':r'$a_{2z}$',
405 'theta1':r'$\theta_1\,(\mathrm{rad})$',
406 'theta2':r'$\theta_2\,(\mathrm{rad})$',
407 'phi1':r'$\phi_1\,(\mathrm{rad})$',
408 'phi2':r'$\phi_2\,(\mathrm{rad})$',
409 'chi_eff':r'$\chi_\mathrm{eff}$',
410 'chi_tot':r'$\chi_\mathrm{total}$',
411 'chi_p':r'$\chi_\mathrm{P}$',
412 'tilt1':r'$t_1\,(\mathrm{rad})$',
413 'tilt2':r'$t_2\,(\mathrm{rad})$',
414 'costilt1':r'$\mathrm{cos}(t_1)$',
415 'costilt2':r'$\mathrm{cos}(t_2)$',
416 'iota':r'$\iota\,(\mathrm{rad})$',
417 'cosiota':r'$\mathrm{cos}(\iota)$',
418 'time':r'$t_\mathrm{c}\,(\mathrm{s})$',
419 'time_mean':r'$<t>\,(\mathrm{s})$',
420 'dist':r'$d_\mathrm{L}\,(\mathrm{Mpc})$',
421 'ra':r'$\alpha$',
422 'dec':r'$\delta$',
423 'phase':r'$\phi\,(\mathrm{rad})$',
424 'psi':r'$\psi\,(\mathrm{rad})$',
425 'theta_jn':r'$\theta_\mathrm{JN}\,(\mathrm{rad})$',
426 'costheta_jn':r'$\mathrm{cos}(\theta_\mathrm{JN})$',
427 'beta':r'$\beta\,(\mathrm{rad})$',
428 'cosbeta':r'$\mathrm{cos}(\beta)$',
429 'phi_jl':r'$\phi_\mathrm{JL}\,(\mathrm{rad})$',
430 'phi12':r'$\phi_\mathrm{12}\,(\mathrm{rad})$',
431 'logl':r'$\mathrm{log}(\mathcal{L})$',
432 'h1_end_time':r'$t_\mathrm{H}$',
433 'l1_end_time':r'$t_\mathrm{L}$',
434 'v1_end_time':r'$t_\mathrm{V}$',
435 'h1l1_delay':r'$\Delta t_\mathrm{HL}$',
436 'h1v1_delay':r'$\Delta t_\mathrm{HV}$',
437 'l1v1_delay':r'$\Delta t_\mathrm{LV}$',
438 'lambdat' : r'$\tilde{\Lambda}$',
439 'dlambdat': r'$\delta \tilde{\Lambda}$',
440 'lambda1' : r'$\lambda_1$',
441 'lambda2': r'$\lambda_2$',
442 'lam_tilde' : r'$\tilde{\Lambda}$',
443 'dlam_tilde': r'$\delta \tilde{\Lambda}$',
444 'calamp_h1' : r'$\delta A_{H1}$',
445 'calamp_l1' : r'$\delta A_{L1}$',
446 'calpha_h1' : r'$\delta \phi_{H1}$',
447 'calpha_l1' : r'$\delta \phi_{L1}$',
448 'polar_eccentricity':r'$\epsilon_{polar}$',
449 'polar_angle':r'$\alpha_{polar}$',
450 'alpha':r'$\alpha_{polar}$',
451 'dchi0':r'$d\chi_0$',
452 'dchi1':r'$d\chi_1$',
453 'dchi2':r'$d\chi_2$',
454 'dchi3':r'$d\chi_3$',
455 'dchi4':r'$d\chi_4$',
456 'dchi5':r'$d\chi_5$',
457 'dchi5l':r'$d\chi_{5}^{(l)}$',
458 'dchi6':r'$d\chi_6$',
459 'dchi6l':r'$d\chi_{6}^{(l)}$',
460 'dchi7':r'$d\chi_7$',
461 'dxi1':r'$d\xi_1$',
462 'dxi2':r'$d\xi_2$',
463 'dxi3':r'$d\xi_3$',
464 'dxi4':r'$d\xi_4$',
465 'dxi5':r'$d\xi_5$',
466 'dxi6':r'$d\xi_6$',
467 'dalpha1':r'$d\alpha_1$',
468 'dalpha2':r'$d\alpha_2$',
469 'dalpha3':r'$d\alpha_3$',
470 'dalpha4':r'$d\alpha_4$',
471 'dalpha5':r'$d\alpha_5$',
472 'dbeta1':r'$d\beta_1$',
473 'dbeta2':r'$d\beta_2$',
474 'dbeta3':r'$d\beta_3$',
475 'dsigma1':r'$d\sigma_1$',
476 'dsigma2':r'$d\sigma_2$',
477 'dsigma3':r'$d\sigma_3$',
478 'dsigma4':r'$d\sigma_4$',
479 }
480
481
482 if param in m1_names:
483 param = 'm1'
484 elif param in m2_names:
485 param = 'm2'
486 elif param in mc_names:
487 param = 'mc'
488 elif param in eta_names:
489 param = 'eta'
490 elif param in q_names:
491 param = 'q'
492 elif param in iota_names:
493 param = 'iota'
494 elif param in dist_names:
495 param = 'dist'
496 elif param in ra_names:
497 param = 'ra'
498 elif param in dec_names:
499 param = 'dec'
500 elif param in phase_names:
501 param = 'phase'
502
503 try:
504 label = labels[param]
505 except KeyError:
506
507 label = param
508
509 return label
510
511
512
513
514
515 -class PosteriorOneDPDF(object):
516 """
517 A data structure representing one parameter in a chain of posterior samples.
518 The Posterior class generates instances of this class for pivoting onto a given
519 parameter (the Posterior class is per-Sampler oriented whereas this class represents
520 the same one parameter in successive samples in the chain).
521 """
522 - def __init__(self,name,posterior_samples,injected_value=None,injFref=None,trigger_values=None,prior=None):
523 """
524 Create an instance of PosteriorOneDPDF based on a table of posterior_samples.
525
526 @param name: A literal string name for the parameter.
527 @param posterior_samples: A 1D array of the samples.
528 @keyword injected_value: The injected or real value of the parameter.
529 @type injected_value: glue.ligolw.lsctables.SimInspiral
530 @keyword trigger_values: The trigger values of the parameter (dictionary w/ IFOs as keys).
531 @keyword prior: The prior value corresponding to each sample.
532 """
533 self.__name=name
534 self.__posterior_samples=np.array(posterior_samples)
535
536 self.__injFref=injFref
537 self.__injval=injected_value
538 self.__trigvals=trigger_values
539 self.__prior=prior
540
541 return
542
544 """
545 Container method. Defined as number of samples.
546 """
547 return len(self.__posterior_samples)
548
549 - def __getitem__(self,idx):
550 """
551 Container method . Returns posterior containing sample idx (allows slicing).
552 """
553 return PosteriorOneDPDF(self.__name, self.__posterior_samples[idx], injected_value=self.__injval, f_ref=self.__f_ref, trigger_values=self.__trigvals)
554
555 @property
557 """
558 Return the string literal name of the parameter.
559
560 @rtype: string
561 """
562 return self.__name
563
564 @property
566 """
567 Return the arithmetic mean for the marginal PDF on the parameter.
568
569 @rtype: number
570 """
571 return np.mean(self.__posterior_samples)
572
573 @property
575 """
576 Return the median value for the marginal PDF on the parameter.
577
578 @rtype: number
579 """
580 return np.median(self.__posterior_samples)
581
582 @property
584 """
585 Return the standard deviation of the marginal PDF on the parameter.
586
587 @rtype: number
588 """
589 try:
590 stdev = sqrt(np.var(self.__posterior_samples))
591 if not np.isfinite(stdev):
592 raise OverflowError
593 except OverflowError:
594 mean = np.mean(self.__posterior_samples)
595 stdev = mean * sqrt(np.var(self.__posterior_samples/mean))
596 return stdev
597
598 @property
600 """
601 Return the 'standard accuracy statistic' (stacc) of the marginal
602 posterior of the parameter.
603
604 stacc is a standard deviant incorporating information about the
605 accuracy of the waveform recovery. Defined as the mean of the sum
606 of the squared differences between the points in the PDF
607 (x_i - sampled according to the posterior) and the true value
608 (x_{\rm true}). So for a marginalized one-dimensional PDF:
609 stacc = \sqrt{\frac{1}{N}\sum_{i=1}^N (x_i-x_{\rm true})2}
610
611 @rtype: number
612 """
613 if self.__injval is None:
614 return None
615 else:
616 return np.sqrt(np.mean((self.__posterior_samples - self.__injval)**2.0))
617
618 @property
620 """
621 Return the injected value set at construction . If no value was set
622 will return None .
623
624 @rtype: undefined
625 """
626 return self.__injval
627
628 @property
629 - def trigvals(self):
630 """
631 Return the trigger values set at construction. If no value was set
632 will return None .
633
634 @rtype: undefined
635 """
636 return self.__trigvals
637
638
639 - def set_injval(self,new_injval):
640 """
641 Set the injected/real value of the parameter.
642
643 @param new_injval: The injected/real value to set.
644 @type new_injval: glue.ligolw.lsctables.SimInspiral
645 """
646
647 self.__injval=new_injval
648
649 - def set_trigvals(self,new_trigvals):
650 """
651 Set the trigger values of the parameter.
652
653 @param new_trigvals: Dictionary containing trigger values with IFO keys.
654 @type new_trigvals: glue.ligolw.lsctables.SnglInspiral
655 """
656
657 self.__trigvals=new_trigvals
658
659 @property
661 """
662 Return a 1D numpy.array of the samples.
663
664 @rtype: numpy.array
665 """
666 return self.__posterior_samples
667
668 - def delete_samples_by_idx(self,samples):
669 """
670 Remove samples from posterior, analagous to numpy.delete but opperates in place.
671
672 @param samples: A list containing the indexes of the samples to be removed.
673 @type samples: list
674 """
675 self.__posterior_samples=np.delete(self.__posterior_samples,samples).reshape(-1,1)
676
677 @property
678 - def gaussian_kde(self):
679 """
680 Return a SciPy gaussian_kde (representing a Gaussian KDE) of the samples.
681
682 @rtype: scipy.stats.kde.gaussian_kde
683 """
684 from numpy import seterr as np_seterr
685 from scipy import seterr as sp_seterr
686
687 np_seterr(under='ignore')
688 sp_seterr(under='ignore')
689 try:
690 return_value=stats.kde.gaussian_kde(np.transpose(self.__posterior_samples))
691 except:
692 exfile=open('exception.out','w')
693 np.savetxt(exfile,self.__posterior_samples)
694 exfile.close()
695 raise
696
697 return return_value
698
699 @property
701 """Returns the KL divergence between the prior and the posterior.
702 It measures the relative information content in nats. The prior is evaluated
703 at run time. It defaults to None. If None is passed, it just returns the information content
704 of the posterior."
705 """
706
707 def uniform(x):
708 return np.array([1./(np.max(x)-np.min(x)) for _ in x])
709
710 posterior, dx = np.histogram(self.samples,bins=36,normed=True)
711 from scipy.stats import entropy
712
713 prior = get_prior(self.name)
714 if prior is None:
715 raise ValueError
716 elif prior=='uniform':
717 prior+='(self.samples)'
718 elif 'x' in prior:
719 prior.replace('x','self.samples')
720 elif not(prior.startswith('np.')):
721 prior = 'np.'+prior
722 prior+='(self.samples)'
723 else:
724 raise ValueError
725
726 try:
727 prior_dist = eval(prior)
728 except:
729 raise ValueError
730
731 return entropy(posterior, qk=prior_dist)
732
733 - def prob_interval(self,intervals):
734 """
735 Evaluate probability intervals.
736
737 @param intervals: A list of the probability intervals [0-1]
738 @type intervals: list
739 """
740 list_of_ci=[]
741 samples_temp=np.sort(np.squeeze(self.samples))
742
743 for interval in intervals:
744 if interval<1.0:
745 samples_temp
746 N=np.size(samples_temp)
747
748 lower_idx=int(floor((N/2.0)*(1-interval)))
749 if lower_idx<0:
750 lower_idx=0
751
752 upper_idx=N-int(floor((N/2.0)*(1-interval)))
753 if upper_idx>N:
754 upper_idx=N-1
755
756 list_of_ci.append((float(samples_temp[lower_idx]),float(samples_temp[upper_idx])))
757 else:
758 list_of_ci.append((None,None))
759
760 return list_of_ci
761
762 -class Posterior(object):
763 """
764 Data structure for a table of posterior samples .
765 """
766 - def __init__(self,commonResultsFormatData,SimInspiralTableEntry=None,inj_spin_frame='OrbitalL', injFref=100,SnglInpiralList=None,name=None,description=None,votfile=None):
767 """
768 Constructor.
769
770 @param commonResultsFormatData: A 2D array containing the posterior
771 samples and related data. The samples chains form the columns.
772 @type commonResultsFormatData: custom type
773 @param SimInspiralTableEntry: A SimInspiralTable row containing the injected values.
774 @type SimInspiralTableEntry: glue.ligolw.lsctables.SimInspiral
775 @param SnglInspiralList: A list of SnglInspiral objects containing the triggers.
776 @type SnglInspiralList: list
777
778 """
779 common_output_table_header,common_output_table_raw =commonResultsFormatData
780 self._posterior={}
781 self._injFref=injFref
782 self._injection=SimInspiralTableEntry
783
784 self._triggers=SnglInpiralList
785 self._loglaliases=['deltalogl', 'posterior', 'logl','logL','likelihood']
786 self._logpaliases=['logp', 'logP','prior','logprior','Prior','logPrior']
787 self._votfile=votfile
788
789 common_output_table_header=[i.lower() for i in common_output_table_header]
790
791
792 self._injXMLFuncMap={
793 'mchirp':lambda inj:inj.mchirp,
794 'chirpmass':lambda inj:inj.mchirp,
795 'mc':lambda inj:inj.mchirp,
796 'mass1':lambda inj:inj.mass1,
797 'm1':lambda inj:inj.mass1,
798 'mass2':lambda inj:inj.mass2,
799 'm2':lambda inj:inj.mass2,
800 'mtotal':lambda inj:float(inj.mass1)+float(inj.mass2),
801 'eta':lambda inj:inj.eta,
802 'q':self._inj_q,
803 'asym_massratio':self._inj_q,
804 'massratio':lambda inj:inj.eta,
805 'sym_massratio':lambda inj:inj.eta,
806 'time': lambda inj:float(inj.get_end()),
807 'time_mean': lambda inj:float(inj.get_end()),
808 'end_time': lambda inj:float(inj.get_end()),
809 'phi0':lambda inj:inj.phi0,
810 'phi_orb': lambda inj: inj.coa_phase,
811 'dist':lambda inj:inj.distance,
812 'distance':lambda inj:inj.distance,
813 'ra':self._inj_longitude,
814 'rightascension':self._inj_longitude,
815 'long':self._inj_longitude,
816 'longitude':self._inj_longitude,
817 'dec':lambda inj:inj.latitude,
818 'declination':lambda inj:inj.latitude,
819 'lat':lambda inj:inj.latitude,
820 'latitude':lambda inj:inj.latitude,
821 'psi': lambda inj: np.mod(inj.polarization, np.pi),
822 'f_ref': lambda inj: self._injFref,
823 'polarisation':lambda inj:inj.polarization,
824 'polarization':lambda inj:inj.polarization,
825 'h1_end_time':lambda inj:float(inj.get_end('H')),
826 'l1_end_time':lambda inj:float(inj.get_end('L')),
827 'v1_end_time':lambda inj:float(inj.get_end('V')),
828 'lal_amporder':lambda inj:inj.amp_order}
829
830
831 for key, val in self._inj_spins(self._injection, frame=inj_spin_frame).items():
832 self._injXMLFuncMap[key] = val
833
834 for one_d_posterior_samples,param_name in zip(np.hsplit(common_output_table_raw,common_output_table_raw.shape[1]),common_output_table_header):
835
836 self._posterior[param_name]=PosteriorOneDPDF(param_name.lower(),one_d_posterior_samples,injected_value=self._getinjpar(param_name),injFref=self._injFref,trigger_values=self._gettrigpar(param_name))
837
838 if 'mchirp' in common_output_table_header and 'eta' in common_output_table_header \
839 and (not 'm1' in common_output_table_header) and (not 'm2' in common_output_table_header):
840 try:
841 print 'Inferring m1 and m2 from mchirp and eta'
842 (m1,m2)=mc2ms(self._posterior['mchirp'].samples, self._posterior['eta'].samples)
843 self._posterior['m1']=PosteriorOneDPDF('m1',m1,injected_value=self._getinjpar('m1'),trigger_values=self._gettrigpar('m1'))
844 self._posterior['m2']=PosteriorOneDPDF('m2',m2,injected_value=self._getinjpar('m2'),trigger_values=self._gettrigpar('m2'))
845 except KeyError:
846 print 'Unable to deduce m1 and m2 from input columns'
847
848
849 logLFound=False
850
851 for loglalias in self._loglaliases:
852
853 if loglalias in common_output_table_header:
854 try:
855 self._logL=self._posterior[loglalias].samples
856 except KeyError:
857 print "No '%s' column in input table!"%loglalias
858 continue
859 logLFound=True
860
861 if not logLFound:
862 raise RuntimeError("No likelihood/posterior values found!")
863 self._logP=None
864
865 for logpalias in self._logpaliases:
866 if logpalias in common_output_table_header:
867 try:
868 self._logP=self._posterior[logpalias].samples
869 except KeyError:
870 print "No '%s' column in input table!"%logpalias
871 continue
872 if not 'log' in logpalias:
873 self._logP=[np.log(i) for i in self._logP]
874
875 if name is not None:
876 self.__name=name
877
878 if description is not None:
879 self.__description=description
880
881 return
882
884 """
885 Add some usefule derived parameters (such as tilt angles, time delays, etc) in the Posterior object
886 """
887 injection=self._injection
888 pos=self
889
890 if 'mc' in pos.names:
891 mchirp_name = 'mc'
892 elif 'chirpmass' in pos.names:
893 mchirp_name = 'chirpmass'
894 else:
895 mchirp_name = 'mchirp'
896
897 if 'asym_massratio' in pos.names:
898 q_name = 'asym_massratio'
899 else:
900 q_name = 'q'
901
902 if 'sym_massratio' in pos.names:
903 eta_name= 'sym_massratio'
904 elif 'massratio' in pos.names:
905 eta_name= 'massratio'
906 else:
907 eta_name='eta'
908
909 if 'mass1' in pos.names and 'mass2' in pos.names :
910 pos.append_mapping(('m1','m2'),lambda x,y:(x,y) ,('mass1','mass2'))
911
912 if (mchirp_name in pos.names and eta_name in pos.names) and \
913 ('mass1' not in pos.names or 'm1' not in pos.names) and \
914 ('mass2' not in pos.names or 'm2' not in pos.names):
915
916 pos.append_mapping(('m1','m2'),mc2ms,(mchirp_name,eta_name))
917
918 if (mchirp_name in pos.names and q_name in pos.names) and \
919 ('mass1' not in pos.names or 'm1' not in pos.names) and \
920 ('mass2' not in pos.names or 'm2' not in pos.names):
921
922 pos.append_mapping(('m1','m2'),q2ms,(mchirp_name,q_name))
923 pos.append_mapping('eta',q2eta,(mchirp_name,q_name))
924
925 if ('m1' in pos.names and 'm2' in pos.names and not 'mtotal' in pos.names ):
926 pos.append_mapping('mtotal', lambda m1,m2: m1+m2, ('m1','m2') )
927
928 if('a_spin1' in pos.names): pos.append_mapping('a1',lambda a:a,'a_spin1')
929 if('a_spin2' in pos.names): pos.append_mapping('a2',lambda a:a,'a_spin2')
930 if('phi_spin1' in pos.names): pos.append_mapping('phi1',lambda a:a,'phi_spin1')
931 if('phi_spin2' in pos.names): pos.append_mapping('phi2',lambda a:a,'phi_spin2')
932 if('theta_spin1' in pos.names): pos.append_mapping('theta1',lambda a:a,'theta_spin1')
933 if('theta_spin2' in pos.names): pos.append_mapping('theta2',lambda a:a,'theta_spin2')
934
935 my_ifos=['h1','l1','v1']
936 for ifo1,ifo2 in combinations(my_ifos,2):
937 p1=ifo1+'_cplx_snr_arg'
938 p2=ifo2+'_cplx_snr_arg'
939 if p1 in pos.names and p2 in pos.names:
940 delta=np.mod(pos[p1].samples - pos[p2].samples + np.pi ,2.0*np.pi)-np.pi
941 pos.append(PosteriorOneDPDF(ifo1+ifo2+'_relative_phase',delta))
942
943
944
945
946
947
948
949
950
951
952 try:
953 if ('ra' in pos.names or 'rightascension' in pos.names) \
954 and ('declination' in pos.names or 'dec' in pos.names) \
955 and 'time' in pos.names:
956 from pylal import xlal,inject
957 from pylal import date
958 from pylal.date import XLALTimeDelayFromEarthCenter
959 from pylal.xlal.datatypes.ligotimegps import LIGOTimeGPS
960 import itertools
961 from numpy import array
962 detMap = {'H1': 'LHO_4k', 'H2': 'LHO_2k', 'L1': 'LLO_4k',
963 'G1': 'GEO_600', 'V1': 'VIRGO', 'T1': 'TAMA_300'}
964 if 'ra' in pos.names:
965 ra_name='ra'
966 else: ra_name='rightascension'
967 if 'dec' in pos.names:
968 dec_name='dec'
969 else: dec_name='declination'
970 ifo_times={}
971 my_ifos=['H1','L1','V1']
972 for ifo in my_ifos:
973 inj_time=None
974 if injection:
975 inj_time=float(injection.get_end(ifo[0]))
976 location=inject.cached_detector[detMap[ifo]].location
977 ifo_times[ifo]=array(map(lambda ra,dec,time: array([time[0]+XLALTimeDelayFromEarthCenter(location,ra[0],dec[0],LIGOTimeGPS(float(time[0])))]), pos[ra_name].samples,pos[dec_name].samples,pos['time'].samples))
978 loc_end_time=PosteriorOneDPDF(ifo.lower()+'_end_time',ifo_times[ifo],injected_value=inj_time)
979 pos.append(loc_end_time)
980 for ifo1 in my_ifos:
981 for ifo2 in my_ifos:
982 if ifo1==ifo2: continue
983 delay_time=ifo_times[ifo2]-ifo_times[ifo1]
984 if injection:
985 inj_delay=float(injection.get_end(ifo2[0])-injection.get_end(ifo1[0]))
986 else:
987 inj_delay=None
988 time_delay=PosteriorOneDPDF(ifo1.lower()+ifo2.lower()+'_delay',delay_time,inj_delay)
989 pos.append(time_delay)
990 except ImportError:
991 print 'Warning: Could not import lal python bindings, check you ./configured with --enable-swig-python'
992 print 'This means I cannot calculate time delays'
993
994
995 new_spin_params = ['tilt1','tilt2','theta_jn','beta']
996 if not set(new_spin_params).issubset(set(pos.names)):
997 old_params = ['f_ref',mchirp_name,'eta','iota','a1','theta1','phi1']
998 if 'a2' in pos.names: old_params += ['a2','theta2','phi2']
999 try:
1000 pos.append_mapping(new_spin_params, spin_angles, old_params)
1001 except KeyError:
1002 print "Warning: Cannot find spin parameters. Skipping spin angle calculations."
1003
1004
1005 if ('a1' in pos.names and 'tilt1' in pos.names and 'm1' in pos.names ) and ('a2' in pos.names and 'tilt2' in pos.names and 'm2' in pos.names):
1006 pos.append_mapping('chi_p', chi_precessing, ['a1', 'tilt1', 'm1', 'a2', 'tilt2', 'm2'])
1007
1008
1009
1010 if('distance' in pos.names):
1011 pos.append_mapping('redshift', calculate_redshift, 'distance')
1012
1013
1014 if ('m1' in pos.names) and ('redshift' in pos.names):
1015 pos.append_mapping('m1_source', source_mass, ['m1', 'redshift'])
1016
1017 if ('m2' in pos.names) and ('redshift' in pos.names):
1018 pos.append_mapping('m2_source', source_mass, ['m2', 'redshift'])
1019
1020 if ('mtotal' in pos.names) and ('redshift' in pos.names):
1021 pos.append_mapping('mtotal_source', source_mass, ['mtotal', 'redshift'])
1022
1023 if ('mc' in pos.names) and ('redshift' in pos.names):
1024 pos.append_mapping('mc_source', source_mass, ['mc', 'redshift'])
1025
1026
1027 if 'a1' in pos.names:
1028 if 'tilt1' in pos.names:
1029 pos.append_mapping('a1z', lambda a, tilt: a*np.cos(tilt), ('a1','tilt1'))
1030 else:
1031 pos.append_mapping('a1z', lambda x: x, 'a1')
1032 inj_az = None
1033 if injection is not None:
1034 inj_az = injection.spin1z
1035 pos['a1z'].set_injval(inj_az)
1036 pos.pop('a1')
1037 pos.append_mapping('a1', lambda x: np.abs(x), 'a1z')
1038
1039 if 'a2' in pos.names:
1040 if 'tilt2' in pos.names:
1041 pos.append_mapping('a2z', lambda a, tilt: a*np.cos(tilt), ('a2','tilt2'))
1042 else:
1043 pos.append_mapping('a2z', lambda x: x, 'a2')
1044 inj_az = None
1045 if injection is not None:
1046 inj_az = injection.spin2z
1047 pos['a2z'].set_injval(inj_az)
1048 pos.pop('a2')
1049 pos.append_mapping('a2', lambda x: np.abs(x), 'a2z')
1050
1051
1052 if ('m1' in pos.names and 'a1z' in pos.names) and ('m2' in pos.names and 'a2z' in pos.names):
1053 pos.append_mapping('chi_eff', lambda m1,a1z,m2,a2z: (m1*a1z + m2*a2z) / (m1 + m2), ('m1','a1z','m2','a2z'))
1054
1055
1056 if ('m1' in pos.names and 'a1' in pos.names and 'tilt1' in pos.names) and ('m2' in pos.names and 'a2' in pos.names and 'tilt2' in pos.names):
1057 pos.append_mapping('chi_tot', lambda m1,a1,m2,a2: (m1*a1 + m2*a2) / (m1 + m2), ('m1','a1','m2','a2'))
1058
1059
1060 if ('m1' in pos.names and 'a1' in pos.names and 'tilt1' in pos.names) and ('m2' in pos.names and 'a2' in pos.names and 'tilt2' in pos.names):
1061 pos.append_mapping('chi_p', chi_precessing, ['m1', 'a1', 'tilt1', 'm2', 'a2', 'tilt2'])
1062
1063
1064 if('distance' in pos.names):
1065 pos.append_mapping('redshift', calculate_redshift, 'distance')
1066
1067
1068 if ('m1' in pos.names) and ('redshift' in pos.names):
1069 pos.append_mapping('m1_source', source_mass, ['m1', 'redshift'])
1070
1071 if ('m2' in pos.names) and ('redshift' in pos.names):
1072 pos.append_mapping('m2_source', source_mass, ['m2', 'redshift'])
1073
1074 if ('mtotal' in pos.names) and ('redshift' in pos.names):
1075 pos.append_mapping('mtotal_source', source_mass, ['mtotal', 'redshift'])
1076
1077 if ('mc' in pos.names) and ('redshift' in pos.names):
1078 pos.append_mapping('mc_source', source_mass, ['mc', 'redshift'])
1079
1080
1081 new_tidal_params = ['lam_tilde','dlam_tilde']
1082 old_tidal_params = ['lambda1','lambda2','eta']
1083 if 'lambda1' in pos.names or 'lambda2' in pos.names:
1084 try:
1085 pos.append_mapping(new_tidal_params, symm_tidal_params, old_tidal_params)
1086 except KeyError:
1087 print "Warning: Cannot find tidal parameters. Skipping tidal calculations."
1088
1089
1090 old_spin_params = ['iota', 'theta1', 'phi1', 'theta2', 'phi2', 'beta']
1091 new_spin_params = ['theta_jn', 'phi_jl', 'tilt1', 'tilt2', 'phi12', 'a1', 'a2', 'm1', 'm2', 'f_ref']
1092 try:
1093 if pos['f_ref'].samples[0][0]==0.0:
1094 for name in ['flow','f_lower']:
1095 if name in pos.names:
1096 new_spin_params = ['theta_jn', 'phi_jl', 'tilt1', 'tilt2', 'phi12', 'a1', 'a2', 'm1', 'm2', name]
1097 except:
1098 print "No f_ref for SimInspiralTransformPrecessingNewInitialConditions()."
1099 if set(new_spin_params).issubset(set(pos.names)) and not set(old_spin_params).issubset(set(pos.names)):
1100 pos.append_mapping(old_spin_params, physical2radiationFrame, new_spin_params)
1101
1102
1103 if 'spin1' in pos.names:
1104 inj_a1 = inj_a2 = None
1105 if injection:
1106 inj_a1 = sqrt(injection.spin1x*injection.spin1x + injection.spin1y*injection.spin1y + injection.spin1z*injection.spin1z)
1107 inj_a2 = sqrt(injection.spin2x*injection.spin2x + injection.spin2y*injection.spin2y + injection.spin2z*injection.spin2z)
1108
1109 try:
1110 a1_samps = abs(pos['spin1'].samples)
1111 a1_pos = PosteriorOneDPDF('a1',a1_samps,injected_value=inj_a1)
1112 pos.append(a1_pos)
1113 except KeyError:
1114 print "Warning: problem accessing spin1 values."
1115
1116 try:
1117 a2_samps = abs(pos['spin2'].samples)
1118 a2_pos = PosteriorOneDPDF('a2',a2_samps,injected_value=inj_a2)
1119 pos.append(a2_pos)
1120 except KeyError:
1121 print "Warning: no spin2 values found."
1122
1123
1124 if ('m1' in pos.names) and ('m2' in pos.names):
1125 if ('tilt1' in pos.names) or ('tilt2' in pos.names):
1126 print "A precessing fit formula is not available yet. Using a non-precessing fit formula on the aligned-spin components."
1127 if ('a1z' in pos.names) and ('a2z' in pos.names):
1128 print "Using non-precessing fit formula [Healy at al (2014)] for final mass and spin (on masses and projected spin components)."
1129 try:
1130 pos.append_mapping('af', bbh_final_spin_non_precessing_Healyetal, ['m1', 'm2', 'a1z', 'a2z'])
1131 pos.append_mapping('mf', lambda m1, m2, chi1z, chi2z, chif: bbh_final_mass_non_precessing_Healyetal(m1, m2, chi1z, chi2z, chif=chif), ['m1', 'm2', 'a1z', 'a2z', 'af'])
1132 except Exception,e:
1133 print "Could not calculate final parameters. The error was: %s"%(str(e))
1134 elif ('a1' in pos.names) and ('a2' in pos.names):
1135 if ('tilt1' in pos.names) and ('tilt2' in pos.names):
1136 print "Projecting spin and using non-precessing fit formula [Healy at al (2014)] for final mass and spin."
1137 try:
1138 pos.append_mapping('af', bbh_final_spin_projected_spin_Healyetal, ['m1', 'm2', 'a1', 'a2', 'tilt1', 'tilt2'])
1139 pos.append_mapping('mf', bbh_final_mass_projected_spin_Healyetal, ['m1', 'm2', 'a1', 'a2', 'tilt1', 'tilt2', 'af'])
1140 except Exception,e:
1141 print "Could not calculate final parameters. The error was: %s"%(str(e))
1142 else:
1143 print "Using non-precessing fit formula [Healy at al (2014)] for final mass and spin (on masses and spin magnitudes)."
1144 try:
1145 pos.append_mapping('af', bbh_final_spin_non_precessing_Healyetal, ['m1', 'm2', 'a1', 'a2'])
1146 pos.append_mapping('mf', lambda m1, m2, chi1, chi2, chif: bbh_final_mass_non_precessing_Healyetal(m1, m2, chi1, chi2, chif=chif), ['m1', 'm2', 'a1', 'a2', 'af'])
1147 except Exception,e:
1148 print "Could not calculate final parameters. The error was: %s"%(str(e))
1149 else:
1150 print "Using non-spinning fit formula [Pan at al (2010)] for final mass and spin."
1151 try:
1152 pos.append_mapping('af', bbh_final_spin_non_spinning_Panetal, ['m1', 'm2'])
1153 pos.append_mapping('mf', bbh_final_mass_non_spinning_Panetal, ['m1', 'm2'])
1154 except Exception,e:
1155 print "Could not calculate final parameters. The error was: %s"%(str(e))
1156 if ('mf' in pos.names) and ('redshift' in pos.names):
1157 try:
1158 pos.append_mapping('mf_source', source_mass, ['mf', 'redshift'])
1159 except Exception,e:
1160 print "Could not calculate final source frame mass. The error was: %s"%(str(e))
1161
1162
1163 if ('mtotal_source' in pos.names) and ('mf_source' in pos.names):
1164 try:
1165 pos.append_mapping('e_rad', lambda mtot_s, mf_s: mtot_s-mf_s, ['mtotal_source', 'mf_source'])
1166 except Exception,e:
1167 print "Could not calculate radiated energy. The error was: %s"%(str(e))
1168
1169 if ('q' in pos.names) and ('a1z' in pos.names) and ('a2z' in pos.names):
1170 try:
1171 pos.append_mapping('l_peak', bbh_aligned_Lpeak_6mode_SHXJDK, ['q', 'a1z', 'a2z'])
1172 except Exception,e:
1173 print "Could not calculate peak luminosity. The error was: %s"%(str(e))
1174
1175 - def bootstrap(self):
1176 """
1177 Returns a new Posterior object that contains a bootstrap
1178 sample of self.
1179
1180 @rtype: Posterior
1181 """
1182 names=[]
1183 samples=[]
1184 for name,oneDpos in self._posterior.items():
1185 names.append(name)
1186 samples.append(oneDpos.samples)
1187
1188 samplesBlock=np.hstack(samples)
1189
1190 bootstrapSamples=samplesBlock[:,:]
1191 Nsamp=bootstrapSamples.shape[0]
1192
1193 rows=np.vsplit(samplesBlock,Nsamp)
1194
1195 for i in range(Nsamp):
1196 bootstrapSamples[i,:]=random.choice(rows)
1197
1198 return Posterior((names,bootstrapSamples),self._injection,self._triggers)
1199
1200 - def delete_samples_by_idx(self,samples):
1201 """
1202 Remove samples from all OneDPosteriors.
1203
1204 @param samples: The indexes of the samples to be removed.
1205 @type samples: list
1206 """
1207 for name,pos in self:
1208 pos.delete_samples_by_idx(samples)
1209 return
1210
1211 - def delete_NaN_entries(self,param_list):
1212 """
1213 Remove samples containing NaN in request params.
1214
1215 @param param_list: The parameters to be checked for NaNs.
1216 @type param_list: list
1217 """
1218 nan_idxs = np.array(())
1219 nan_dict = {}
1220 for param in param_list:
1221 nan_bool_array = np.isnan(self[param].samples).any(1)
1222 idxs = np.where(nan_bool_array == True)[0]
1223 if len(idxs) > 0:
1224 nan_dict[param]=len(idxs)
1225 nan_idxs = np.append(nan_idxs, idxs)
1226 total_samps = len(self)
1227 nan_samps = len(nan_idxs)
1228 if nan_samps is not 0:
1229 print "WARNING: removing %i of %i total samples due to NaNs:"% (nan_samps,total_samps)
1230 for param in nan_dict.keys():
1231 print "\t%i NaNs in %s."%(nan_dict[param],param)
1232 self.delete_samples_by_idx(nan_idxs)
1233 return
1234
1235 @property
1237 """Returns the Deviance Information Criterion estimated from the
1238 posterior samples. The DIC is defined as -2*(<log(L)> -
1239 Var(log(L))); smaller values are "better."
1240
1241 """
1242
1243 return -2.0*(np.mean(self._logL) - np.var(self._logL))
1244
1245 @property
1246 - def injection(self):
1247 """
1248 Return the injected values.
1249
1250 @rtype: glue.ligolw.lsctables.SimInspiral
1251 """
1252
1253 return self._injection
1254
1255 @property
1256 - def triggers(self):
1257 """
1258 Return the trigger values .
1259
1260 @rtype: list
1261 """
1262
1263 return self._triggers
1264
1265 - def _total_incl_restarts(self, samples):
1266 total=0
1267 last=samples[0]
1268 for x in samples[1:]:
1269 if x < last:
1270 total += last
1271 last = x
1272 total += samples[-1]
1273 return total
1274
1276 """
1277 Returns the number of cycles in the longest chain
1278
1279 @rtype: number
1280 """
1281 samps,header=self.samples()
1282 header=header.split()
1283 if not ('cycle' in header):
1284 raise RuntimeError("Cannot compute number of cycles in longest chain")
1285
1286 cycle_col=header.index('cycle')
1287 if 'chain' in header:
1288 chain_col=header.index('chain')
1289 chain_indexes=np.unique(samps[:,chain_col])
1290 max_cycle=0
1291 for ind in chain_indexes:
1292 chain_cycle_samps=samps[ samps[:,chain_col] == ind, cycle_col ]
1293 max_cycle=max(max_cycle, self._total_incl_restarts(chain_cycle_samps))
1294 return int(max_cycle)
1295 else:
1296 return int(self._total_incl_restarts(samps[:,cycle_col]))
1297
1298
1299 - def set_injection(self,injection):
1300 """
1301 Set the injected values of the parameters.
1302
1303 @param injection: A SimInspiralTable row object containing the injected parameters.
1304 @type injection: glue.ligolw.lsctables.SimInspiral
1305 """
1306 if injection is not None:
1307 self._injection=injection
1308 for name,onepos in self:
1309 new_injval=self._getinjpar(name)
1310 if new_injval is not None:
1311 self[name].set_injval(new_injval)
1312
1313 - def set_triggers(self,triggers):
1314 """
1315 Set the trigger values of the parameters.
1316
1317 @param triggers: A list of SnglInspiral objects.
1318 @type triggers: list
1319 """
1320 if triggers is not None:
1321 self._triggers=triggers
1322 for name,onepos in self:
1323 new_trigvals=self._gettrigpar(name)
1324 if new_trigvals is not None:
1325 self[name].set_trigvals(new_trigvals)
1326
1327
1328 - def _getinjpar(self,paramname):
1329 """
1330 Map parameter names to parameters in a SimInspiralTable .
1331 """
1332 if self._injection is not None:
1333 for key,value in self._injXMLFuncMap.items():
1334 if paramname.lower().strip() == key.lower().strip():
1335 try:
1336 return self._injXMLFuncMap[key](self._injection)
1337 except TypeError:
1338 return self._injXMLFuncMap[key]
1339 return None
1340
1341 - def _gettrigpar(self,paramname):
1342 """
1343 Map parameter names to parameters in a SnglInspiral.
1344 """
1345 vals = None
1346 if self._triggers is not None:
1347 for key,value in self._injXMLFuncMap.items():
1348 if paramname.lower().strip() == key.lower().strip():
1349 try:
1350 vals = dict([(trig.ifo,self._injXMLFuncMap[key](trig)) for trig in self._triggers])
1351 except TypeError:
1352 return self._injXMLFuncMap[key]
1353 except AttributeError:
1354 return None
1355 return vals
1356
1357 - def __getitem__(self,key):
1358 """
1359 Container method . Returns posterior chain,one_d_pos, with name one_d_pos.name.
1360 """
1361 return self._posterior[key.lower()]
1362
1363 - def __len__(self):
1364 """
1365 Container method. Defined as number of samples.
1366 """
1367 return len(self._logL)
1368
1369 - def __iter__(self):
1370 """
1371 Container method. Returns iterator from self.forward for use in
1372 for (...) in (...) etc.
1373 """
1374 return self.forward()
1375
1376 - def forward(self):
1377 """
1378 Generate a forward iterator (in sense of list of names) over Posterior
1379 with name,one_d_pos.
1380 """
1381 current_item = 0
1382 while current_item < self.dim:
1383 name=self._posterior.keys()[current_item]
1384 pos=self._posterior[name]
1385 current_item += 1
1386 yield name,pos
1387
1388 - def bySample(self):
1389 """
1390 Generate a forward iterator over the list of samples corresponding to
1391 the data stored within the Posterior instance. These are returned as
1392 ParameterSamples instances.
1393 """
1394 current_item=0
1395 pos_array,header=self.samples
1396 while current_item < len(self):
1397 sample_array=(np.squeeze(pos_array[current_item,:]))
1398 yield PosteriorSample(sample_array, header, header)
1399 current_item += 1
1400
1401
1402 @property
1404 """
1405 Return number of parameters.
1406 """
1407 return len(self._posterior.keys())
1408
1409 @property
1411 """
1412 Return list of parameter names.
1413 """
1414 nameslist=[]
1415 for key,value in self:
1416 nameslist.append(key)
1417 return nameslist
1418
1419 @property
1421 """
1422 Return dict {paramName:paramMean} .
1423 """
1424 meansdict={}
1425 for name,pos in self:
1426 meansdict[name]=pos.mean
1427 return meansdict
1428
1429 @property
1431 """
1432 Return dict {paramName:paramMedian} .
1433 """
1434 mediansdict={}
1435 for name,pos in self:
1436 mediansdict[name]=pos.median
1437 return mediansdict
1438
1439 @property
1441 """
1442 Return dict {paramName:paramStandardDeviation} .
1443 """
1444 stdsdict={}
1445 for name,pos in self:
1446 stdsdict[name]=pos.stdev
1447 return stdsdict
1448
1449 @property
1451 """
1452 Return qualified string containing the 'name' of the Posterior instance.
1453 """
1454 return self.__name
1455
1456 @property
1457 - def description(self):
1458 """
1459 Return qualified string containing a 'description' of the Posterior instance.
1460 """
1461 return self.__description
1462
1463 - def append(self,one_d_posterior):
1464 """
1465 Container method. Add a new OneDParameter to the Posterior instance.
1466 """
1467 self._posterior[one_d_posterior.name]=one_d_posterior
1468 return
1469
1470 - def pop(self,param_name):
1471 """
1472 Container method. Remove PosteriorOneDPDF from the Posterior instance.
1473 """
1474 return self._posterior.pop(param_name)
1475
1476 - def append_mapping(self, new_param_names, func, post_names):
1477 """
1478 Append posteriors pos1,pos2,...=func(post_names)
1479 """
1480
1481 import copy
1482
1483 if isinstance(post_names, str):
1484 old_post = copy.deepcopy(self[post_names])
1485 old_inj = old_post.injval
1486 old_trigs = old_post.trigvals
1487 if old_inj:
1488 new_inj = func(old_inj)
1489 else:
1490 new_inj = None
1491 if old_trigs:
1492 new_trigs = {}
1493 for IFO in old_trigs.keys():
1494 new_trigs[IFO] = func(old_trigs[IFO])
1495 else:
1496 new_trigs = None
1497
1498 samps = func(old_post.samples)
1499 new_post = PosteriorOneDPDF(new_param_names, samps, injected_value=new_inj, trigger_values=new_trigs)
1500 if new_post.samples.ndim is 0:
1501 print "WARNING: No posterior calculated for %s ..." % post.name
1502 else:
1503 self.append(new_post)
1504
1505 else:
1506 old_posts = [copy.deepcopy(self[post_name]) for post_name in post_names]
1507 old_injs = [post.injval for post in old_posts]
1508 old_trigs = [post.trigvals for post in old_posts]
1509 samps = func(*[post.samples for post in old_posts])
1510
1511 if isinstance(new_param_names, str):
1512 if None not in old_injs:
1513 inj = func(*old_injs)
1514 else:
1515 inj = None
1516 if None not in old_trigs:
1517 new_trigs = {}
1518 for IFO in old_trigs[0].keys():
1519 oldvals = [param[IFO] for param in old_trigs]
1520 new_trigs[IFO] = func(*oldvals)
1521 else:
1522 new_trigs = None
1523 new_post = PosteriorOneDPDF(new_param_names, samps, injected_value=inj, trigger_values=new_trigs)
1524 self.append(new_post)
1525
1526 else:
1527 if None not in old_injs:
1528 injs = func(*old_injs)
1529 else:
1530 injs = [None for name in new_param_names]
1531 if None not in old_trigs:
1532 new_trigs = [{} for param in range(len(new_param_names))]
1533 for IFO in old_trigs[0].keys():
1534 oldvals = [param[IFO] for param in old_trigs]
1535 newvals = func(*oldvals)
1536 for param,newval in enumerate(newvals):
1537 new_trigs[param][IFO] = newval
1538 else:
1539 new_trigs = [None for param in range(len(new_param_names))]
1540 if not samps: return()
1541 new_posts = [PosteriorOneDPDF(new_param_name,samp,injected_value=inj,trigger_values=new_trigs) for (new_param_name,samp,inj,new_trigs) in zip(new_param_names,samps,injs,new_trigs)]
1542 for post in new_posts:
1543 if post.samples.ndim is 0:
1544 print "WARNING: No posterior calculated for %s ..." % post.name
1545 else:
1546 self.append(post)
1547 return
1548
1549 - def _average_posterior(self, samples, post_name):
1550 """
1551 Returns the average value of the 'post_name' column of the
1552 given samples.
1553 """
1554 ap = 0.0
1555 for samp in samples:
1556 ap = ap + samp[post_name]
1557 return ap / len(samples)
1558
1559 - def _average_posterior_like_prior(self, samples, logl_name, prior_name, log_bias = 0):
1560 """
1561 Returns the average value of the posterior assuming that the
1562 'logl_name' column contains log(L) and the 'prior_name' column
1563 contains the prior (un-logged).
1564 """
1565 ap = 0.0
1566 for samp in samples:
1567 ap += np.exp(samp[logl_name]-log_bias)*samp[prior_name]
1568 return ap / len(samples)
1569
1570 - def _bias_factor(self):
1571 """
1572 Returns a sensible bias factor for the evidence so that
1573 integrals are representable as doubles.
1574 """
1575 return np.mean(self._logL)
1576
1577 - def di_evidence(self, boxing=64):
1578 """
1579 Returns the log of the direct-integration evidence for the
1580 posterior samples.
1581 """
1582 allowed_coord_names=["spin1", "spin2", "a1", "phi1", "theta1", "a2", "phi2", "theta2",
1583 "iota", "psi", "ra", "dec",
1584 "phi_orb", "phi0", "dist", "time", "mc", "mchirp", "chirpmass", "q"]
1585 samples,header=self.samples()
1586 header=header.split()
1587 coord_names=[name for name in allowed_coord_names if name in header]
1588 coordinatized_samples=[PosteriorSample(row, header, coord_names) for row in samples]
1589 tree=KDTree(coordinatized_samples)
1590
1591 if "prior" in header and "logl" in header:
1592 bf = self._bias_factor()
1593 return bf + np.log(tree.integrate(lambda samps: self._average_posterior_like_prior(samps, "logl", "prior", bf), boxing))
1594 elif "prior" in header and "likelihood" in header:
1595 bf = self._bias_factor()
1596 return bf + np.log(tree.integrate(lambda samps: self._average_posterior_like_prior(samps, "likelihood", "prior", bf), boxing))
1597 elif "post" in header:
1598 return np.log(tree.integrate(lambda samps: self._average_posterior(samps, "post"), boxing))
1599 elif "posterior" in header:
1600 return np.log(tree.integrate(lambda samps: self._average_posterior(samps, "posterior"), boxing))
1601 else:
1602 raise RuntimeError("could not find 'post', 'posterior', 'logl' and 'prior', or 'likelihood' and 'prior' columns in output to compute direct integration evidence")
1603
1605 """Returns an approximation to the log(evidence) obtained by
1606 fitting an ellipse around the highest-posterior samples and
1607 performing the harmonic mean approximation within the ellipse.
1608 Because the ellipse should be well-sampled, this provides a
1609 better approximation to the evidence than the full-domain HM."""
1610 allowed_coord_names=["spin1", "spin2", "a1", "phi1", "theta1", "a2", "phi2", "theta2",
1611 "iota", "psi", "ra", "dec",
1612 "phi_orb", "phi0", "dist", "time", "mc", "mchirp", "chirpmass", "q"]
1613 samples,header=self.samples()
1614 header=header.split()
1615
1616 n=int(0.05*samples.shape[0])
1617 if not n > 1:
1618 raise IndexError
1619
1620 coord_names=[name for name in allowed_coord_names if name in header]
1621 indexes=np.argsort(self._logL[:,0])
1622
1623 my_samples=samples[indexes[-n:], :]
1624 my_samples=np.array([PosteriorSample(sample,header,coord_names).coord() for sample in my_samples])
1625
1626 mu=np.mean(my_samples, axis=0)
1627 cov=np.cov(my_samples, rowvar=0)
1628
1629 d0=None
1630 for mysample in my_samples:
1631 d=np.dot(mysample-mu, np.linalg.solve(cov, mysample-mu))
1632 if d0 is None:
1633 d0 = d
1634 else:
1635 d0=max(d0,d)
1636
1637 ellipse_logl=[]
1638 ellipse_samples=[]
1639 for sample,logl in zip(samples, self._logL):
1640 coord=PosteriorSample(sample, header, coord_names).coord()
1641 d=np.dot(coord-mu, np.linalg.solve(cov, coord-mu))
1642
1643 if d <= d0:
1644 ellipse_logl.append(logl)
1645 ellipse_samples.append(sample)
1646
1647 if len(ellipse_samples) > 5*n:
1648 print 'WARNING: ellpise evidence region encloses significantly more samples than %d'%n
1649
1650 ellipse_samples=np.array(ellipse_samples)
1651 ellipse_logl=np.array(ellipse_logl)
1652
1653 ndim = len(coord_names)
1654 ellipse_volume=np.pi**(ndim/2.0)*d0**(ndim/2.0)/special.gamma(ndim/2.0+1)*np.sqrt(np.linalg.det(cov))
1655
1656 try:
1657 prior_index=header.index('prior')
1658 pmu=np.mean(ellipse_samples[:,prior_index])
1659 pstd=np.std(ellipse_samples[:,prior_index])
1660 if pstd/pmu > 1.0:
1661 print 'WARNING: prior variation greater than 100\% over elliptical volume.'
1662 approx_prior_integral=ellipse_volume*pmu
1663 except KeyError:
1664
1665 approx_prior_integral=ellipse_volume
1666
1667 ll_bias=np.mean(ellipse_logl)
1668 ellipse_logl = ellipse_logl - ll_bias
1669
1670 return np.log(approx_prior_integral) - np.log(np.mean(1.0/np.exp(ellipse_logl))) + ll_bias
1671
1673 """
1674 Returns the log of the harmonic mean evidence for the set of
1675 posterior samples.
1676 """
1677 bf = self._bias_factor()
1678 return bf + np.log(1/np.mean(1/np.exp(self._logL-bf)))
1679
1680 - def _posMaxL(self):
1681 """
1682 Find the sample with maximum likelihood probability. Returns value
1683 of likelihood and index of sample .
1684 """
1685 logl_vals=self._logL
1686 max_i=0
1687 max_logl=logl_vals[0]
1688 for i in range(len(logl_vals)):
1689 if logl_vals[i] > max_logl:
1690 max_logl=logl_vals[i]
1691 max_i=i
1692 return max_logl,max_i
1693
1694 - def _posMap(self):
1695 """
1696 Find the sample with maximum a posteriori probability. Returns value
1697 of posterior and index of sample .
1698 """
1699 logl_vals=self._logL
1700 if self._logP is not None:
1701 logp_vals=self._logP
1702 else:
1703 return None
1704
1705 max_i=0
1706 max_pos=logl_vals[0]+logp_vals[0]
1707 for i in range(len(logl_vals)):
1708 if logl_vals[i]+logp_vals[i] > max_pos:
1709 max_pos=logl_vals[i]+logp_vals[i]
1710 max_i=i
1711 return max_pos,max_i
1712
1713 - def _print_table_row(self,name,entries):
1714 """
1715 Print a html table row representation of
1716
1717 name:item1,item2,item3,...
1718 """
1719
1720 row_str='<tr><td>%s</td>'%name
1721 for entry in entries:
1722 row_str+='<td>%s</td>'%entry
1723 row_str+='</tr>'
1724 return row_str
1725
1726 @property
1728 """
1729 Return the maximum likelihood probability and the corresponding
1730 set of parameters.
1731 """
1732 maxLvals={}
1733 max_logl,max_i=self._posMaxL()
1734 for param_name in self.names:
1735 maxLvals[param_name]=self._posterior[param_name].samples[max_i][0]
1736
1737 return (max_logl,maxLvals)
1738
1739 @property
1741 """
1742 Return the maximum a posteriori probability and the corresponding
1743 set of parameters.
1744 """
1745 maxPvals={}
1746 max_pos,max_i=self._posMap()
1747 for param_name in self.names:
1748 maxPvals[param_name]=self._posterior[param_name].samples[max_i][0]
1749
1750 return (max_pos,maxPvals)
1751
1752
1753 - def samples(self):
1754 """
1755 Return an (M,N) numpy.array of posterior samples; M = len(self);
1756 N = dim(self) .
1757 """
1758 header_string=''
1759 posterior_table=[]
1760 for param_name,one_pos in self:
1761 column=np.array(one_pos.samples)
1762 header_string+=param_name+'\t'
1763 posterior_table.append(column)
1764 posterior_table=tuple(posterior_table)
1765 return np.column_stack(posterior_table),header_string
1766
1767 - def write_to_file(self,fname):
1768 """
1769 Dump the posterior table to a file in the 'common format'.
1770 """
1771 column_list=()
1772
1773 posterior_table,header_string=self.samples()
1774
1775 fobj=open(fname,'w')
1776
1777 fobj.write(header_string+'\n')
1778 np.savetxt(fobj,posterior_table,delimiter='\t')
1779 fobj.close()
1780
1781 return
1782
1783 - def gelman_rubin(self, pname):
1784 """
1785 Returns an approximation to the Gelman-Rubin statistic (see
1786 Gelman, A. and Rubin, D. B., Statistical Science, Vol 7,
1787 No. 4, pp. 457--511 (1992)) for the parameter given, accurate
1788 as the number of samples in each chain goes to infinity. The
1789 posterior samples must have a column named 'chain' so that the
1790 different chains can be separated.
1791 """
1792 from numpy import seterr as np_seterr
1793 np_seterr(all='raise')
1794
1795 if "chain" in self.names:
1796 chains=np.unique(self["chain"].samples)
1797 chain_index=self.names.index("chain")
1798 param_index=self.names.index(pname)
1799 data,header=self.samples()
1800 chainData=[data[ data[:,chain_index] == chain, param_index] for chain in chains]
1801 allData=np.concatenate(chainData)
1802 chainMeans=[np.mean(data) for data in chainData]
1803 chainVars=[np.var(data) for data in chainData]
1804 BoverN=np.var(chainMeans)
1805 W=np.mean(chainVars)
1806 sigmaHat2=W + BoverN
1807 m=len(chainData)
1808 VHat=sigmaHat2 + BoverN/m
1809 try:
1810 R = VHat/W
1811 except:
1812 print "Error when computer Gelman-Rubin R statistic for %s. This may be a fixed parameter"%pname
1813 R = np.nan
1814 return R
1815 else:
1816 raise RuntimeError('could not find necessary column header "chain" in posterior samples')
1817
1818 - def healpix_map(self, resol, nest=True):
1819 """Returns a healpix map in the pixel ordering that represents the
1820 posterior density (per square degree) on the sky. The pixels
1821 will be chosen to have at least the given resolution (in
1822 degrees).
1823
1824 """
1825
1826
1827 nside = 2
1828 while hp.nside2resol(nside, arcmin=True) > resol*60.0/2.0:
1829 nside *= 2
1830
1831 ras = self['ra'].samples.squeeze()
1832 decs = self['dec'].samples.squeeze()
1833
1834 phis = ras
1835 thetas = np.pi/2.0 - decs
1836
1837
1838 inds = hp.ang2pix(nside, thetas, phis, nest=False)
1839 counts = np.bincount(inds)
1840 if counts.shape[0] < hp.nside2npix(nside):
1841 counts = np.concatenate((counts, np.zeros(hp.nside2npix(nside) - counts.shape[0])))
1842
1843
1844 hpmap = hp.sphtfunc.smoothing(counts, sigma=resol*np.pi/180.0)
1845
1846 hpmap = hpmap / (np.sum(hpmap)*hp.nside2pixarea(nside, degrees=True))
1847
1848 if nest:
1849 hpmap = hp.reorder(hpmap, r2n=True)
1850
1851 return hpmap
1852
1853 - def __str__(self):
1854 """
1855 Define a string representation of the Posterior class ; returns
1856 a html formatted table of various properties of posteriors.
1857 """
1858 return_val='<table border="1" id="statstable"><tr><th/>'
1859 column_names=['maP','maxL','stdev','mean','median','stacc','injection value']
1860 IFOs = []
1861 if self._triggers is not None:
1862 IFOs = [trig.ifo for trig in self._triggers]
1863 for IFO in IFOs:
1864 column_names.append(IFO+' trigger values')
1865
1866 for column_name in column_names:
1867 return_val+='<th>%s</th>'%column_name
1868
1869 return_val+='</tr>'
1870
1871 for name,oned_pos in self:
1872
1873 max_logl,max_i=self._posMaxL()
1874 maxL=oned_pos.samples[max_i][0]
1875 max_post,max_j=self._posMap()
1876 maP=oned_pos.samples[max_j][0]
1877 mean=str(oned_pos.mean)
1878 stdev=str(oned_pos.stdev)
1879 median=str(np.squeeze(oned_pos.median))
1880 stacc=str(oned_pos.stacc)
1881 injval=str(oned_pos.injval)
1882 trigvals=oned_pos.trigvals
1883
1884 row = [maP,maxL,stdev,mean,median,stacc,injval]
1885 if self._triggers is not None:
1886 for IFO in IFOs:
1887 try:
1888 row.append(str(trigvals[IFO]))
1889 except TypeError:
1890 row.append(None)
1891 return_val+=self._print_table_row(name,row)
1892
1893 return_val+='</table>'
1894
1895 parser=XMLParser()
1896 parser.feed(return_val)
1897 Estr=parser.close()
1898
1899 elem=Estr
1900 rough_string = tostring(elem, 'utf-8')
1901 reparsed = minidom.parseString(rough_string)
1902 return_val=reparsed.toprettyxml(indent=" ")
1903 return return_val[len('<?xml version="1.0" ?>')+1:]
1904
1905 - def write_vot_info(self):
1906 """
1907 Writes the information stored in the VOTtree if there is one
1908 """
1909 target=VOT2HTML()
1910 parser=XMLParser(target=target)
1911 parser.feed(self._votfile)
1912 return parser.close()
1913
1914
1915
1916
1917 - def _inj_m1(self,inj):
1918 """
1919 Return the mapping of (mchirp,eta)->m1; m1>m2 i.e. return the greater of the mass
1920 components (m1) calculated from the chirp mass and the symmetric mass ratio.
1921
1922 @type inj: glue.ligolw.lsctables.SimInspiral
1923 @param inj: a custom type with the attributes 'mchirp' and 'eta'.
1924 @rtype: number
1925 """
1926 (mass1,mass2)=mc2ms(inj.mchirp,inj.eta)
1927 return mass1
1928
1929 - def _inj_m2(self,inj):
1930 """
1931 Return the mapping of (mchirp,eta)->m2; m1>m2 i.e. return the lesser of the mass
1932 components (m2) calculated from the chirp mass and the symmetric mass ratio.
1933
1934 @type inj: glue.ligolw.lsctables.SimInspiral
1935 @param inj: a custom type with the attributes 'mchirp' and 'eta'.
1936 @rtype: number
1937 """
1938 (mass1,mass2)=mc2ms(inj.mchirp,inj.eta)
1939 return mass2
1940
1941 - def _inj_q(self,inj):
1942 """
1943 Return the mapping of (mchirp,eta)->q; m1>m2 i.e. return the mass ratio q=m2/m1.
1944
1945 @type inj: glue.ligolw.lsctables.SimInspiral
1946 @param inj: a custom type with the attributes 'mchirp' and 'eta'.
1947 @rtype: number
1948 """
1949 (mass1,mass2)=mc2ms(inj.mchirp,inj.eta)
1950 return mass2/mass1
1951
1952 - def _inj_longitude(self,inj):
1953 """
1954 Return the mapping of longitude found in inj to the interval [0,2*pi).
1955
1956 @type inj: glue.ligolw.lsctables.SimInspiral
1957 @param inj: a custom type with the attribute 'longitude'.
1958 @rtype: number
1959 """
1960 if inj.longitude>2*pi_constant or inj.longitude<0.0:
1961 maplong=2*pi_constant*(((float(inj.longitude))/(2*pi_constant)) - floor(((float(inj.longitude))/(2*pi_constant))))
1962 print "Warning: Injected longitude/ra (%s) is not within [0,2\pi)! Angles are assumed to be in radians so this will be mapped to [0,2\pi). Mapped value is: %s."%(str(inj.longitude),str(maplong))
1963 return maplong
1964 else:
1965 return inj.longitude
1966
1967 - def _inj_spins(self, inj, frame='OrbitalL'):
1968 spins = {}
1969 f_ref = self._injFref
1970
1971 if not inj:
1972 spins = {}
1973
1974 else:
1975 axis = lalsim.SimInspiralGetFrameAxisFromString(frame)
1976
1977 m1, m2 = inj.mass1, inj.mass2
1978 mc, eta = inj.mchirp, inj.eta
1979
1980
1981 iota, s1x, s1y, s1z, s2x, s2y, s2z = \
1982 lalsim.SimInspiralInitialConditionsPrecessingApproxs(inj.inclination,
1983 inj.spin1x, inj.spin1y, inj.spin1z,
1984 inj.spin2x, inj.spin2y, inj.spin2z,
1985 m1*lal.MSUN_SI, m2*lal.MSUN_SI, f_ref, inj.coa_phase, axis)
1986
1987 a1, theta1, phi1 = cart2sph(s1x, s1y, s1z)
1988 a2, theta2, phi2 = cart2sph(s2x, s2y, s2z)
1989
1990 spins = {'a1':a1, 'theta1':theta1, 'phi1':phi1,
1991 'a2':a2, 'theta2':theta2, 'phi2':phi2,
1992 'iota':iota}
1993
1994
1995 if inj.spin1x == inj.spin1y == inj.spin2x == inj.spin2y == 0.:
1996 spins['a1z'] = inj.spin1z
1997 spins['a2z'] = inj.spin2z
1998
1999 L = orbital_momentum(f_ref, m1,m2, iota)
2000 S1 = np.hstack((s1x, s1y, s1z))
2001 S2 = np.hstack((s2x, s2y, s2z))
2002
2003 zhat = np.array([0., 0., 1.])
2004 aligned_comp_spin1 = array_dot(S1, zhat)
2005 aligned_comp_spin2 = array_dot(S2, zhat)
2006 chi = aligned_comp_spin1 + aligned_comp_spin2 + \
2007 np.sqrt(1. - 4.*eta) * (aligned_comp_spin1 - aligned_comp_spin2)
2008
2009 spins['spinchi'] = chi
2010
2011 S1 *= m1**2
2012 S2 *= m2**2
2013 J = L + S1 + S2
2014
2015 tilt1 = array_ang_sep(L, S1)
2016 tilt2 = array_ang_sep(L, S2)
2017 beta = array_ang_sep(J, L)
2018
2019 spins['tilt1'] = tilt1
2020 spins['tilt2'] = tilt2
2021 spins['beta'] = beta
2022
2023
2024
2025 phi0 = np.arctan2(J[1], J[0])
2026 phi0 = np.pi - phi0
2027
2028 J = ROTATEZ(phi0, J[0], J[1], J[2])
2029 L = ROTATEZ(phi0, L[0], L[1], L[2])
2030 S1 = ROTATEZ(phi0, S1[0], S1[1], S1[2])
2031 S2 = ROTATEZ(phi0, S2[0], S2[1], S2[2])
2032
2033
2034 theta_jn = array_polar_ang(J)
2035 spins['theta_jn'] = theta_jn
2036
2037 J = ROTATEY(theta_jn, J[0], J[1], J[2])
2038 L = ROTATEY(theta_jn, L[0], L[1], L[2])
2039 S1 = ROTATEY(theta_jn, S1[0], S1[1], S1[2])
2040 S2 = ROTATEY(theta_jn, S2[0], S2[1], S2[2])
2041
2042
2043 phi_jl = np.arctan2(L[1], L[0])
2044 phi_jl = np.pi - phi_jl
2045 spins['phi_jl'] = phi_jl
2046
2047
2048 J = ROTATEZ(phi_jl, J[0], J[1], J[2])
2049 L = ROTATEZ(phi_jl, L[0], L[1], L[2])
2050 S1 = ROTATEZ(phi_jl, S1[0], S1[1], S1[2])
2051 S2 = ROTATEZ(phi_jl, S2[0], S2[1], S2[2])
2052
2053 theta0 = array_polar_ang(L)
2054 J = ROTATEY(theta0, J[0], J[1], J[2])
2055 L = ROTATEY(theta0, L[0], L[1], L[2])
2056 S1 = ROTATEY(theta0, S1[0], S1[1], S1[2])
2057 S2 = ROTATEY(theta0, S2[0], S2[1], S2[2])
2058
2059
2060 phi1 = np.arctan2(S1[1], S1[0])
2061 phi2 = np.arctan2(S2[1], S2[0])
2062 if phi2 < phi1:
2063 phi12 = phi2 - phi1 + 2.*np.pi
2064 else:
2065 phi12 = phi2 - phi1
2066 spins['phi12'] = phi12
2067
2068 return spins
2069
2070 -class BurstPosterior(Posterior):
2071 """
2072 Data structure for a table of posterior samples .
2073 """
2074 - def __init__(self,commonResultsFormatData,SimBurstTableEntry=None,injFref=None,SnglBurstList=None,name=None,description=None,votfile=None):
2075 """
2076 Constructor.
2077
2078 @param commonResultsFormatData: A 2D array containing the posterior
2079 samples and related data. The samples chains form the columns.
2080 @type commonResultsFormatData: custom type
2081 @param SimInspiralTableEntry: A SimInspiralTable row containing the injected values.
2082 @type SimInspiralTableEntry: glue.ligolw.lsctables.SimInspiral
2083 @param SnglInspiralList: A list of SnglInspiral objects containing the triggers.
2084 @type SnglInspiralList: list
2085 """
2086 common_output_table_header,common_output_table_raw =commonResultsFormatData
2087 self._posterior={}
2088 self._injFref=injFref
2089 self._injection=SimBurstTableEntry
2090 self._triggers=SnglBurstList
2091 self._loglaliases=['posterior', 'logl','logL','likelihood', 'deltalogl']
2092 self._logpaliases=['logp', 'logP','prior','logprior','Prior','logPrior']
2093
2094 self._votfile=votfile
2095
2096 common_output_table_header=[i.lower() for i in common_output_table_header]
2097
2098
2099 self._injXMLFuncMap={
2100 'f0':lambda inj:inj.frequency,
2101 'frequency':lambda inj:inj.frequency,
2102 'centre_frequency':lambda inj:inj.frequency,
2103 'quality':lambda inj:inj.q,
2104 'hrss':lambda inj:inj.hrss,
2105 'loghrss':lambda inj:log(inj.hrss),
2106 'polar_angle':lambda inj:inj.pol_ellipse_angle,
2107 'pol_ellipse_angle':lambda inj:inj.pol_ellipse_angle,
2108 'pol_ellipse_e':lambda inj:inj.pol_ellipse_e,
2109 'alpha':lambda inj:inj.pol_ellipse_angle,
2110 'polar_eccentricity':lambda inj:inj.pol_ellipse_e,
2111 'eccentricity':lambda inj:inj.pol_ellipse_e,
2112 'time': lambda inj:float(inj.get_end()),
2113 'end_time': lambda inj:float(inj.get_end()),
2114 'ra':self._inj_longitude,
2115 'rightascension':self._inj_longitude,
2116 'long':self._inj_longitude,
2117 'longitude':self._inj_longitude,
2118 'dec':lambda inj:inj.dec,
2119 'declination':lambda inj:inj.dec,
2120 'lat':lambda inj:inj.dec,
2121 'latitude':lambda inj:inj.dec,
2122 'psi': lambda inj: np.mod(inj.psi, np.pi),
2123 'f_ref': lambda inj: self._injFref,
2124 'polarisation':lambda inj:inj.psi,
2125 'polarization':lambda inj:inj.psi,
2126 'duration':lambda inj:inj.duration,
2127 'h1_end_time':lambda inj:float(inj.get_end('H')),
2128 'l1_end_time':lambda inj:float(inj.get_end('L')),
2129 'v1_end_time':lambda inj:float(inj.get_end('V')),
2130 }
2131
2132 for one_d_posterior_samples,param_name in zip(np.hsplit(common_output_table_raw,common_output_table_raw.shape[1]),common_output_table_header):
2133
2134 self._posterior[param_name]=PosteriorOneDPDF(param_name.lower(),one_d_posterior_samples,injected_value=self._getinjpar(param_name),injFref=self._injFref,trigger_values=self._gettrigpar(param_name))
2135
2136 logLFound=False
2137
2138 for loglalias in self._loglaliases:
2139 if loglalias in common_output_table_header:
2140 try:
2141 self._logL=self._posterior[loglalias].samples
2142 except KeyError:
2143 print "No '%s' column in input table!"%loglalias
2144 continue
2145 logLFound=True
2146
2147 if not logLFound:
2148 raise RuntimeError("No likelihood/posterior values found!")
2149
2150 self._logP=None
2151 for logpalias in self._logpaliases:
2152 if logpalias in common_output_table_header:
2153 try:
2154 self._logP=self._posterior[logpalias].samples
2155 except KeyError:
2156 print "No '%s' column in input table!"%logpalias
2157 continue
2158 if not 'log' in logpalias:
2159 self._logP=[np.log(i) for i in self._logP]
2160 if name is not None:
2161 self.__name=name
2162
2163 if description is not None:
2164 self.__description=description
2165
2166 return
2167
2168
2169
2170
2171 - def _inj_longitude(self,inj):
2172 """
2173 Return the mapping of longitude found in inj to the interval [0,2*pi).
2174
2175 @type inj: glue.ligolw.lsctables.SimInspiral
2176 @param inj: a custom type with the attribute 'longitude'.
2177 @rtype: number
2178 """
2179 if inj.ra>2*pi_constant or inj.ra<0.0:
2180 maplong=2*pi_constant*(((float(inj.ra)/(2*pi_constant)) - floor(((float(inj.ra))/(2*pi_constant)))))
2181 print "Warning: Injected longitude/ra (%s) is not within [0,2\pi)! Angles are assumed to be in radians so this will be mapped to [0,2\pi). Mapped value is: %s."%(str(inj.ra),str(maplong))
2182 return maplong
2183 else:
2184 return inj.ra
2185
2187 """
2188 A kD-tree.
2189 """
2191 """
2192 Construct a kD-tree from a sequence of objects. Each object
2193 should return its coordinates using obj.coord().
2194 """
2195 if len(objects) == 0:
2196 raise RuntimeError("cannot construct kD-tree out of zero objects---you may have a repeated sample in your list")
2197 elif len(objects) == 1:
2198 self._objects = objects[:]
2199 coord=self._objects[0].coord()
2200 self._bounds = coord,coord
2201 elif self._same_coords(objects):
2202
2203 self._objects = [ objects[0] ]
2204 coord=self._objects[0].coord()
2205 self._bounds = coord,coord
2206 else:
2207 self._objects = objects[:]
2208 self._bounds = self._bounds_of_objects()
2209 low,high=self._bounds
2210 self._split_dim=self._longest_dimension()
2211 longest_dim = self._split_dim
2212 sorted_objects=sorted(self._objects, key=lambda obj: (obj.coord())[longest_dim])
2213 N = len(sorted_objects)
2214 bound=0.5*(sorted_objects[N/2].coord()[longest_dim] + sorted_objects[N/2-1].coord()[longest_dim])
2215 low = [obj for obj in self._objects if obj.coord()[longest_dim] < bound]
2216 high = [obj for obj in self._objects if obj.coord()[longest_dim] >= bound]
2217 if len(low)==0:
2218
2219
2220 low = [obj for obj in self._objects if obj.coord()[longest_dim]==bound]
2221 high = [obj for obj in self._objects if obj.coord()[longest_dim] > bound]
2222 self._left = KDTree(low)
2223 self._right = KDTree(high)
2224
2226 """
2227 True if and only if all the given objects have the same
2228 coordinates.
2229 """
2230 if len(objects) <= 1:
2231 return True
2232 coords = [obj.coord() for obj in objects]
2233 c0 = coords[0]
2234 for ci in coords[1:]:
2235 if not np.all(ci == c0):
2236 return False
2237 return True
2238
2240 """
2241 Bounds of the objects contained in the tree.
2242 """
2243 low=self._objects[0].coord()
2244 high=self._objects[0].coord()
2245 for obj in self._objects[1:]:
2246 low=np.minimum(low,obj.coord())
2247 high=np.maximum(high,obj.coord())
2248 return low,high
2249
2251 """
2252 Longest dimension of the tree bounds.
2253 """
2254 low,high = self._bounds
2255 widths = high-low
2256 return np.argmax(widths)
2257
2259 """
2260 Returns the objects in the tree.
2261 """
2262 return self._objects[:]
2263
2265 """
2266 Iterator over all the objects contained in the tree.
2267 """
2268 return self._objects.__iter__()
2269
2271 """
2272 Returns the left tree.
2273 """
2274 return self._left
2275
2277 """
2278 Returns the right tree.
2279 """
2280 return self._right
2281
2283 """
2284 Returns the dimension along which this level of the kD-tree
2285 splits.
2286 """
2287 return self._split_dim
2288
2290 """
2291 Returns the coordinates of the lower-left and upper-right
2292 corners of the bounding box for this tree: low_left, up_right
2293 """
2294 return self._bounds
2295
2297 """
2298 Returns the volume of the bounding box of the tree.
2299 """
2300 v = 1.0
2301 low,high=self._bounds
2302 for l,h in zip(low,high):
2303 v = v*(h - l)
2304 return v
2305
2307 """
2308 Returns the integral of f(objects) over the tree. The
2309 optional boxing parameter determines how deep to descend into
2310 the tree before computing f.
2311 """
2312
2313
2314
2315
2316
2317 def x(tree):
2318 return tree.volume()*f(tree._objects)
2319
2320 def y(a,b):
2321 return a+b
2322
2323 return self.operate(x,y,boxing=boxing)
2324
2326 """
2327 Operates on tree nodes exceeding boxing parameter depth.
2328 """
2329 if len(self._objects) <= boxing:
2330 return f(self)
2331 else:
2332
2333 return g(self._left.operate(f,g,boxing),self._right.operate(f,g,boxing))
2334
2337 """
2338 A kD-tree suitable for splitting parameter spaces and counting hypervolumes.
2339 Is modified from the KDTree class so that bounding boxes are stored. This means that
2340 there are no longer gaps in the hypervolume once the samples have been split into groups.
2341 """
2342 - def __init__(self, objects,boundingbox,dims=0):
2343 """
2344 Construct a kD-tree from a sequence of objects. Each object
2345 should return its coordinates using obj.coord().
2346 the obj should also store the bounds of the hypervolume its found in.
2347 for non-leaf objects we need the name of the dimension split and value at split.
2348 """
2349 self._dimension = dims
2350 self._bounds = boundingbox
2351 self._weight = 1
2352 if len(objects) == 0:
2353 raise RuntimeError("cannot construct kD-tree out of zero objects---you may have a repeated sample in your list")
2354 elif len(objects) == 1:
2355 self._objects = objects[:]
2356 elif self._same_coords(objects):
2357 self._weight = len(objects)
2358 self._objects = [ objects[0] ]
2359 coord=self._objects[0].coord()
2360 else:
2361 self._objects = objects[:]
2362 split_dim = self._dimension
2363 sorted_objects=sorted(self._objects, key=lambda obj: (obj.coord())[split_dim])
2364 N = len(sorted_objects)
2365 self._split_value = 0.5*(sorted_objects[N/2].coord()[split_dim] + sorted_objects[N/2-1].coord()[split_dim])
2366 bound = self._split_value
2367 low = [obj for obj in self._objects if obj.coord()[split_dim] < bound]
2368 high = [obj for obj in self._objects if obj.coord()[split_dim] >= bound]
2369 if len(low)==0:
2370
2371
2372 low = [obj for obj in self._objects if obj.coord()[split_dim] == bound]
2373 high = [obj for obj in self._objects if obj.coord()[split_dim] > bound]
2374 leftBoundingbox = []
2375 rightBoundingbox = []
2376 for i in self._bounds:
2377 leftBoundingbox.append(list(i))
2378 rightBoundingbox.append(list(i))
2379 leftBoundingbox[1][split_dim] = bound
2380 rightBoundingbox[0][split_dim] = bound
2381
2382
2383
2384 if (split_dim < (len(self._objects[0].coord()) - 1)):
2385 child_dim = split_dim + 1
2386 else:
2387 child_dim = 0
2388 self._left = KDTreeVolume(low,leftBoundingbox,dims = child_dim)
2389
2390 if (len(high) != 0):
2391 self._right = KDTreeVolume(high,rightBoundingbox,dims = child_dim)
2392 else:
2393 self._right = None
2394
2396 """
2397 True if and only if all the given objects have the same
2398 coordinates.
2399 """
2400 if len(objects) <= 1:
2401 return True
2402 coords = [obj.coord() for obj in objects]
2403 c0 = coords[0]
2404 for ci in coords[1:]:
2405 if not np.all(ci == c0):
2406 return False
2407 return True
2408
2410 """
2411 Returns the objects in the tree.
2412 """
2413 return self._objects[:]
2414
2416 """
2417 Iterator over all the objects contained in the tree.
2418 """
2419 return self._objects.__iter__()
2420
2422 """
2423 Returns the left tree.
2424 """
2425 return self._left
2426
2428 """
2429 Returns the right tree.
2430 """
2431 return self._right
2432
2434 """
2435 Returns the dimension along which this level of the kD-tree
2436 splits.
2437 """
2438 return self._split_dim
2439
2441 """
2442 Returns the coordinates of the lower-left and upper-right
2443 corners of the bounding box for this tree: low_left, up_right
2444 """
2445 return self._bounds
2446
2448 """
2449 Returns the volume of the bounding box of the tree.
2450 """
2451 v = 1.0
2452 low,high=self._bounds
2453 for l,h in zip(low,high):
2454 v = v*(h - l)
2455 return v
2456
2458 """
2459 Returns the integral of f(objects) over the tree. The
2460 optional boxing parameter determines how deep to descend into
2461 the tree before computing f.
2462 """
2463 def x(tree):
2464 return tree.volume()*f(tree._objects)
2465
2466 def y(a,b):
2467 return a+b
2468
2469 return self.operate(x,y,boxing=boxing)
2470
2472 """
2473 Operates on tree nodes exceeding boxing parameter depth.
2474 """
2475 if len(self._objects) <= boxing:
2476 return f(self)
2477 else:
2478 return g(self._left.operate(f,g,boxing),self._right.operate(f,g,boxing))
2479
2480 - def search(self,coordinates,boxing = 64):
2481 """
2482 takes a set of coordinates and searches down through the tree untill it gets
2483 to a box with less than 'boxing' objects in it and returns the box bounds,
2484 number of objects in the box, and the weighting.
2485 """
2486 if len(self._objects) <= boxing:
2487 return self._bounds,len(self._objects),self._weight
2488 elif coordinates[self._dimension] < self._split_value:
2489 return self._left.search(coordinates,boxing)
2490 else:
2491 return self._right.search(coordinates,boxing)
2492
2494 """
2495 copies tree structure, but with KDSkeleton as the new nodes.
2496 """
2497 boxN = boxing
2498 if len(self._objects) <= boxN:
2499 newNode = KDSkeleton(self.bounds(), left_child = None , right_child = None)
2500 if isArea:
2501 newNode.setImportance(len(self._objects),skyArea(self.bounds()))
2502 else:
2503 newNode.setImportance(len(self._objects),self.volume())
2504 return newNode
2505 else:
2506 if isArea:
2507 newNode = KDSkeleton(self.bounds, left_child = self._left.fillNewTree(boxN,isArea=True), right_child = self._right.fillNewTree(boxN,isArea=True))
2508 newNode.setImportance(len(self._objects),skyArea(self.bounds()))
2509 else:
2510 newNode = KDSkeleton(self.bounds, left_child = self._left.fillNewTree(boxN), right_child = self._right.fillNewTree(boxN))
2511 newNode.setImportance(len(self._objects),self.volume())
2512 newNode.setSplit(self._dimension,self._split_value)
2513 return newNode
2514
2516 """
2517 object to store the structure of a kd tree
2518 """
2519
2520 - def __init__(self, bounding_box, left_child = None, right_child = None):
2521 self._bounds = bounding_box
2522
2523 self._left = left_child
2524 self._right = right_child
2525 self._samples = 0
2526 self._splitValue = None
2527 self._splitDim = None
2528 self._importance = None
2529 self._volume = None
2530
2533
2536
2537 - def search(self,coordinates):
2538 """
2539 takes a set of coordinates and searches down through the tree untill it gets
2540 to a box with less than 'boxing' objects in it and returns the box bounds,
2541 number of objects in the box, and the weighting.
2542 """
2543 if self._left is None:
2544 return self._bounds, self._samples, self._importance
2545 elif coordinates[self._splitDim] < self._splitValue:
2546 return self._left.search(coordinates)
2547 else:
2548 return self._right.search(coordinates)
2549
2551 self._importance = sampleNumber/volume
2552 self._volume = volume
2553
2555 self._splitDim = dimension
2556 self._splitValue = value
2557
2558
2559 -class PosteriorSample(object):
2560 """
2561 A single parameter sample object, suitable for inclusion in a
2562 kD-tree.
2563 """
2564
2565 - def __init__(self, sample_array, headers, coord_names):
2566 """
2567 Given the sample array, headers for the values, and the names
2568 of the desired coordinates, construct a parameter sample
2569 object.
2570 """
2571 self._samples=sample_array[:]
2572 self._headers=headers
2573 if not (len(sample_array) == len(self._headers)):
2574 print "Header length = ", len(self._headers)
2575 print "Sample length = ", len(sample_array)
2576 raise RuntimeError("parameter and sample lengths do not agree")
2577 self._coord_names=coord_names
2578 self._coord_indexes=[self._headers.index(name) for name in coord_names]
2579
2580 - def __getitem__(self, key):
2581 """
2582 Return the element with the corresponding name.
2583 """
2584 key=key.lower()
2585 if key in self._headers:
2586 idx=self._headers.index(key)
2587 return self._samples[idx]
2588 else:
2589 raise KeyError("key not found in posterior sample: %s"%key)
2590
2592 """
2593 Return the coordinates for the parameter sample.
2594 """
2595 return self._samples[self._coord_indexes]
2596
2602 """
2603 Return analytic likelihood values.
2604 """
2605
2606 - def __init__(self, covariance_matrix_files, mean_vector_files):
2607 """
2608 Prepare analytic likelihood for the given parameters.
2609 """
2610
2611 if isinstance(covariance_matrix_files, str):
2612 covariance_matrix_files = [covariance_matrix_files]
2613 if isinstance(mean_vector_files, str):
2614 mean_vector_files = [mean_vector_files]
2615
2616 covarianceMatrices = [np.loadtxt(csvFile, delimiter=',') for csvFile in covariance_matrix_files]
2617 num_matrices = len(covarianceMatrices)
2618
2619 if num_matrices != len(mean_vector_files):
2620 raise RuntimeError('Must give a mean vector list for every covariance matrix')
2621
2622 param_line = open(mean_vector_files[0]).readline()
2623 self._params = [param.strip() for param in param_line.split(',')]
2624
2625 converter=lambda x: eval(x.replace('pi','%.32f'%pi_constant))
2626 self._modes = []
2627 for i in range(num_matrices):
2628 CM = covarianceMatrices[i]
2629 vecFile = mean_vector_files[i]
2630
2631 param_line = open(vecFile).readline()
2632 params = [param.strip() for param in param_line.split(',')]
2633 if set(params)!=set(self._params):
2634 raise RuntimeError('Parameters do not agree between mean vector files.')
2635
2636 sigmas = dict(zip(params,np.sqrt(CM.diagonal())))
2637 colNums = range(len(params))
2638 converters = dict(zip(colNums,[converter for i in colNums]))
2639 meanVectors = np.loadtxt(vecFile, delimiter=',', skiprows=1, converters=converters)
2640 try:
2641 for vec in meanVectors:
2642 means = dict(zip(params,vec))
2643 mode = [(param, stats.norm(loc=means[param],scale=sigmas[param])) for param in params]
2644 self._modes.append(dict(mode))
2645 except TypeError:
2646 means = dict(zip(params,meanVectors))
2647 mode = [(param, stats.norm(loc=means[param],scale=sigmas[param])) for param in params]
2648 self._modes.append(dict(mode))
2649 self._num_modes = len(self._modes)
2650
2651 - def pdf(self, param):
2652 """
2653 Return PDF function for parameter.
2654 """
2655 pdf = None
2656 if param in self._params:
2657 pdf = lambda x: (1.0/self._num_modes) * sum([mode[param].pdf(x) for mode in self._modes])
2658 return pdf
2659
2660 - def cdf(self, param):
2661 """
2662 Return PDF function for parameter.
2663 """
2664 cdf = None
2665 if param in self._params:
2666 cdf = lambda x: (1.0/self._num_modes) * sum([mode[param].cdf(x) for mode in self._modes])
2667 return cdf
2668
2669 @property
2671 """
2672 Return list of parameter names described by analytic likelihood function.
2673 """
2674 return self._params
2675
2676
2677
2678
2679
2680
2681
2682 -class htmlChunk(object):
2683 """
2684 A base class for representing web content using ElementTree .
2685 """
2686 - def __init__(self,tag,attrib=None,parent=None):
2687
2688 self._html=Element(tag)
2689 if attrib:
2690 for attribname,attribvalue in attrib.items():
2691 self._html.attrib[attribname]=attribvalue
2692 if parent:
2693 parent.append(self._html)
2694
2696 """
2697 Return a pretty-printed XML string of the htmlPage.
2698 """
2699 elem=self._html
2700 rough_string = tostring(elem)
2701 reparsed = minidom.parseString(rough_string)
2702 return reparsed.toprettyxml(indent=" ")
2703
2706
2707 - def write(self,string):
2708 parser=XMLParser()
2709 parser.feed(string)
2710 Estr=parser.close()
2711 self._html.append(Estr)
2712
2713 - def p(self,pstring):
2714 Ep=Element('p')
2715 Ep.text=pstring
2716 self._html.append(Ep)
2717 return Ep
2718
2719 - def h1(self,h1string):
2720 Ep=Element('h1')
2721 Ep.text=h1string
2722 self._html.append(Ep)
2723 return Ep
2724
2725 - def h5(self,h1string):
2726 Ep=Element('h5')
2727 Ep.text=h1string
2728 self._html.append(Ep)
2729 return Ep
2730
2731 - def h2(self,h2string):
2732 Ep=Element('h2')
2733 Ep.text=h2string
2734 self._html.append(Ep)
2735 return Ep
2736
2737 - def h3(self,h1string):
2738 Ep=Element('h3')
2739 Ep.text=h1string
2740 self._html.append(Ep)
2741 return Ep
2742
2744 Ebr=Element('br')
2745 self._html.append(Ebr)
2746 return Ebr
2747
2749 Ehr=Element('hr')
2750 self._html.append(Ehr)
2751 return Ehr
2752
2753 - def a(self,url,linktext):
2754 Ea=Element('a')
2755 Ea.attrib['href']=url
2756 Ea.text=linktext
2757 self._html.append(Ea)
2758 return Ea
2759
2760 - def tab(self,idtable=None):
2761 args={}
2762 if idtable is not None:
2763 args={'id':idtable}
2764
2765 Etab=Element('table',args)
2766 self._html.append(Etab)
2767 return Etab
2768
2770
2771 """
2772 Insert row in table tab.
2773 If given, label used as id for the table tag
2774 """
2775
2776 Etr=Element('tr')
2777 if label is not None:
2778 Etr.attrib['id']=label
2779 tab.append(Etr)
2780 return Etr
2781
2782 - def insert_td(self,row,td,label=None,legend=None):
2783 """
2784 Insert cell td into row row.
2785 Sets id to label, if given
2786 """
2787
2788 Etd=Element('td')
2789
2790 if type(td) is str:
2791 Etd.text=td
2792 else:
2793 td=tostring(td)
2794 td=minidom.parseString(td)
2795 td=td.toprettyxml(indent=" ")
2796 Etd.text=td
2797 if label is not None:
2798 Etd.attrib['id']=label
2799 if legend is not None:
2800 legend.a('#%s'%label,'%s'%label)
2801 legend.br()
2802 row.append(Etd)
2803 return Etd
2804
2806 self._html.append(element)
2807
2808
2809
2810 -class htmlPage(htmlChunk):
2811 """
2812 A concrete class for generating an XHTML(1) document. Inherits from htmlChunk.
2813 """
2814 - def __init__(self,title=None,css=None,javascript=None,toc=False):
2815 htmlChunk.__init__(self,'html',attrib={'xmlns':"http://www.w3.org/1999/xhtml"})
2816 self.doctype_str='<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">'
2817
2818 self._head=SubElement(self._html,'head')
2819 Etitle=SubElement(self._head,'title')
2820 self._body=SubElement(self._html,'body')
2821 self._css=None
2822 self._jscript=None
2823 if title is not None:
2824 Etitle.text=str(title)
2825 self._title=SubElement(self._body,'h1')
2826 self._title.text=title
2827 if javascript is not None:
2828 self._jscript=SubElement(self._head,'script')
2829 self._jscript.attrib['type']="text/javascript"
2830 self._jscript.text=str(javascript)
2831 if css is not None:
2832 self._css=SubElement(self._head,'style')
2833 self._css.attrib['type']="text/css"
2834 self._css.text=str(css)
2835
2836 - def __str__(self):
2837 return self.doctype_str+'\n'+self.toprettyxml()
2838
2839 - def add_section(self,section_name,legend=None):
2840 newSection=htmlSection(section_name)
2841 self._body.append(newSection._html)
2842 if legend is not None:
2843 legend.a('#%s'%section_name,'%s'%section_name)
2844 legend.br()
2845 return newSection
2846
2847 - def add_collapse_section(self,section_name,legend=None,innertable_id=None,start_closed=True):
2848 """
2849 Create a section embedded into a table that can be collapsed with a button
2850 """
2851 newSection=htmlCollapseSection(section_name,table_id=innertable_id,start_closed=start_closed)
2852 self._body.append(newSection._html)
2853 if legend is not None:
2854 legend.a('#%s'%section_name,'%s'%section_name)
2855 legend.br()
2856 return newSection
2857
2858 - def add_section_to_element(self,section_name,parent):
2859 """
2860 Create a section which is not appended to the body of html, but to the parent Element
2861 """
2862 newSection=htmlSection(section_name,htmlElement=parent,blank=True)
2863 parent.append(newSection._html)
2864 return newSection
2865
2866
2867 @property
2870
2871 @property
2874
2877 """
2878 Represents a block of html fitting within a htmlPage. Inherits from htmlChunk.
2879 """
2880 - def __init__(self,section_name,htmlElement=None,blank=False):
2881 if not blank:
2882 htmlChunk.__init__(self,'div',attrib={'class':'ppsection','id':section_name},parent=htmlElement)
2883 else:
2884 htmlChunk.__init__(self,'div',attrib={'style':'"color:#000000"','id':section_name},parent=htmlElement)
2885 self.h2(section_name)
2886
2889 """
2890 Represents a block of html fitting within a htmlPage. Inherits from htmlChunk.
2891 """
2892
2893 - def __init__(self,section_name,htmlElement=None,table_id=None,start_closed=True):
2894 htmlChunk.__init__(self,'div',attrib={'class':'ppsection','id':section_name},parent=htmlElement)
2895
2896 if table_id is None:
2897 table_id=random.randint(1,10000000)
2898 self.table_id=table_id
2899 self._start_closed=start_closed
2900
2901 - def write(self,string):
2902 k=random.randint(1,10000000)
2903 if self._start_closed:
2904 st='<table border="0" align="center" cellpadding="5" cellspacing="0"><tr bgcolor="#4682B4" height="50"><td width="5%%"><font size="4" face="tahoma" color="white"><a href="#"> Top</a></font></td><td width="45%%"><font size="4" face="tahoma" color="white"><strong>%s</strong></font></td><td bgcolor="#4682B4" align="center" width="50%%"><input id="lnk%s" type="button" value="[+] Expand" onclick="toggle_visibility(\'%s\',\'lnk%s\');"></input></td></tr><tr><td colspan="7">'%(self._html.attrib['id'],k,self.table_id,k)
2905 string=string.replace('table ', 'table style="display: none" ')
2906 else:
2907 st='<table border="0" align="center" cellpadding="5" cellspacing="0"><tr bgcolor="#4682B4" height="50"><td width="5%%"><font size="4" face="tahoma" color="white"><a href="#"> Top</a></font></td><td width="45%%"><font size="4" face="tahoma" color="white"><strong>%s</strong></font></td><td bgcolor="#4682B4" align="center" width="50%%"><input id="lnk%s" type="button" value="[-] Collapse" onclick="toggle_visibility(\'%s\',\'lnk%s\');"></input></td></tr><tr><td colspan="7">'%(self._html.attrib['id'],k,self.table_id,k)
2908 string=string.replace('table ', 'table style="display: table" ')
2909 st+=string
2910 st+='</td></tr></table>'
2911 htmlChunk.write(self,st)
2912
2918 """
2919 @deprecated: This is a pure python version of the C extension function
2920 pylal._bayespputils._skyhist_cart .
2921 """
2922
2923 N=len(skycarts)
2924 print 'operating on %d sky points'%(N)
2925 bins=np.zeros(N)
2926 for RAsample,decsample in sky_samples:
2927 sampcart=pol2cart(RAsample,decsample)
2928 maxdx=-1
2929 maxvalue=-1
2930 for i in xrange(0,N):
2931 dx=np.dot(sampcart,skycarts[i])
2932 if dx>maxvalue:
2933 maxdx=i
2934 maxvalue=dx
2935
2936 bins[maxdx]+=1
2937 return bins
2938
2940 """
2941 @deprecated: This is an old pure python version of the C extension function
2942 pylal._bayespputils._skyhist_cart .
2943 """
2944 N=len(skypoints)
2945 print 'operating on %d sky points' % (N)
2946 bins=zeros(N)
2947 j=0
2948 for sample in samples:
2949 seps=map(lambda s: ang_dist(sample[RAdim],sample[decdim],s[1],s[0]),skypoints)
2950 minsep=math.pi
2951 for i in range(0,N):
2952 if seps[i]<minsep:
2953 minsep=seps[i]
2954 mindx=i
2955 bins[mindx]=bins[mindx]+1
2956 j=j+1
2957 print 'Done %d/%d iterations, minsep=%f degrees'\
2958 %(j,len(samples),minsep*(180.0/3.1415926))
2959 return (skypoints,bins)
2960
2963 """
2964 Returns (injectionconf, toppoints), where injectionconf is the
2965 confidence level of the injection, contained in the injBin and
2966 toppoints is a list of (pointx, pointy, ptindex, frac), with
2967 pointx and pointy the (x,y) coordinates of the corresponding
2968 element of the points array, ptindex the index of the point in the
2969 array, and frac the cumulative fraction of points with larger
2970 posterior probability.
2971
2972 The hist argument should be a one-dimensional array that contains
2973 counts of sample points in each bin.
2974
2975 The points argument should be a 2-D array storing the sky location
2976 associated with each bin; the first index runs from 0 to NBins -
2977 1, while the second index runs from 0 to 1.
2978
2979 The injBin argument gives the bin index in which the injection is
2980 found.
2981
2982 The NSamples argument is used to normalize the histogram counts
2983 into fractional probability.
2984 """
2985
2986 histIndices=np.argsort(hist)[::-1]
2987
2988 toppoints=[]
2989 frac=0.0
2990 injConf=None
2991 for i in histIndices:
2992 frac+=float(hist[i])/float(NSamples)
2993 toppoints.append((points[i,0], points[i,1], i, frac))
2994 if i == injBin:
2995 injConf=frac
2996 print 'Injection found at confidence level %g'%injConf
2997
2998 return (injConf, toppoints)
2999
3000 -def _greedy_bin(greedyHist,greedyPoints,injection_bin_index,bin_size,Nsamples,confidence_levels):
3001 """
3002 An interal function representing the common, dimensionally-independent part of the
3003 greedy binning algorithms.
3004 """
3005
3006
3007 (injectionconfidence,toppoints)=_calculate_confidence_levels(greedyHist, greedyPoints, injection_bin_index, Nsamples)
3008
3009
3010 nBins=0
3011 confidence_levels.sort()
3012 reses={}
3013 toppoints=np.array(toppoints)
3014 for printcl in confidence_levels:
3015 nBins=np.searchsorted(toppoints[:,3], printcl) + 1
3016
3017 if nBins >= len(toppoints):
3018 nBins=len(toppoints)-1
3019
3020 accl=toppoints[nBins-1,3]
3021
3022 reses[printcl]=nBins*bin_size
3023
3024
3025 injection_area=None
3026 if injection_bin_index and injectionconfidence:
3027 i=list(np.nonzero(np.asarray(toppoints)[:,2]==injection_bin_index))[0]
3028 injection_area=bin_size*(i+1)
3029
3030 return toppoints,injectionconfidence,reses,injection_area
3031
3036
3038 size = int(len(items)*fraction)
3039 random.shuffle(items)
3040 return items[:size], items[size:]
3041
3043 if tree._left is None:
3044 tree.addSample()
3045 elif coordinates[tree._splitDim] < tree._splitValue:
3046 addSample(tree._left,coordinates)
3047 else:
3048 addSample(tree._right,coordinates)
3049
3056
3057 confidence_levels.sort()
3058
3059 class Harvester(list):
3060
3061 def __init__(self):
3062 list.__init__(self)
3063 self.unrho=0.
3064
3065 def __call__(self,tree):
3066 number_density=float(len(tree.objects()))/float(tree.volume())
3067 self.append([number_density,tree.volume(),tree.bounds()])
3068 self.unrho+=number_density
3069
3070 def close_ranks(self):
3071
3072 for i in range(len(self)):
3073 self[i][0]/=self.unrho
3074
3075 return sorted(self,key=itemgetter(0))
3076
3077 def h(a,b):
3078 pass
3079
3080 samples,header=posterior.samples()
3081 header=header.split()
3082 coord_names=["ra","dec","dist"]
3083 coordinatized_samples=[PosteriorSample(row, header, coord_names) for row in samples]
3084 tree=KDTree(coordinatized_samples)
3085
3086 a=Harvester()
3087 samples_per_bin=10
3088 tree.operate(a,h,boxing=samples_per_bin)
3089
3090 b=a.close_ranks()
3091 b.reverse()
3092
3093 acc_rho=0.
3094 acc_vol=0.
3095 cl_idx=0
3096 confidence_intervals={}
3097 for rho,vol,bounds in b:
3098 acc_rho+=rho
3099 acc_vol+=vol
3100
3101 if acc_rho>confidence_levels[cl_idx]:
3102 confidence_intervals[acc_rho]=acc_vol
3103 cl_idx+=1
3104 if cl_idx==len(confidence_levels):
3105 break
3106
3107 return confidence_intervals
3108
3110 """
3111 takes samples and applies a KDTree to them to return confidence levels
3112 returns confidence_intervals - dictionary of user_provided_CL:calculated_area
3113 b - ordered list of KD leaves
3114 injInfo - if injection values provided then returns
3115 [Bounds_of_inj_kd_leaf ,number_samples_in_box, weight_of_box,injection_CL ,injection_CL_area]
3116 Not quite sure that the repeated samples case is fixed, posibility of infinite loop.
3117 """
3118 confidence_levels.sort()
3119 from math import cos, pi
3120 class Harvester(list):
3121 """
3122 when called by kdtree.operate will be used to calculate the density of each bin (sky area)
3123 """
3124 def __init__(self):
3125 list.__init__(self)
3126 self.unrho=0.
3127
3128 def __call__(self,tree):
3129
3130 area = - (cos(pi/2. - tree.bounds()[0][1])-cos(pi/2. - tree.bounds()[1][1]))*(tree.bounds()[1][0] - tree.bounds()[0][0])
3131 number_density=float(len(tree.objects()))/float(area)
3132 self.append([number_density,len(tree.objects()),area,tree.bounds()])
3133 self.unrho+=number_density
3134
3135 def close_ranks(self):
3136
3137 for i in range(len(self)):
3138 self[i][0]/=self.unrho
3139
3140 return sorted(self,key=itemgetter(0))
3141
3142 def h(a,b):
3143 pass
3144
3145 peparser=PEOutputParser('common')
3146
3147 samples,header=posterior.samples()
3148 header=header.split()
3149 coord_names=["ra","dec"]
3150 initial_dimensions = [[0.,-pi/2.],[2.*pi, pi/2.]]
3151 coordinatized_samples=[PosteriorSample(row, header, coord_names) for row in samples]
3152 tree=KDTreeVolume(coordinatized_samples,initial_dimensions)
3153
3154 a=Harvester()
3155 tree.operate(a,h,boxing=samples_per_bin)
3156 totalSamples = len(tree.objects())
3157 b=a.close_ranks()
3158 b.reverse()
3159 samplecounter=0.0
3160 for entry in b:
3161 samplecounter += entry[1]
3162 entry[1] = float(samplecounter)/float(totalSamples)
3163
3164 acc_rho=0.
3165 acc_vol=0.
3166 cl_idx=0
3167
3168
3169 if posterior['ra'].injval is not None and posterior['dec'].injval is not None:
3170 injBound,injNum,injWeight = tree.search([posterior['ra'].injval,posterior['dec'].injval],boxing = samples_per_bin)
3171 injInfo = [injBound,injNum,injWeight]
3172 inj_area = - (cos(pi/2. - injBound[0][1])-cos(pi/2. - injBound[1][1]))*(injBound[1][0] - injBound[0][0])
3173 inj_number_density=float(injNum)/float(inj_area)
3174 inj_rho = inj_number_density / a.unrho
3175 else:
3176 injInfo = None
3177 inj_area = None
3178 inj_number_density=None
3179 inj_rho = None
3180
3181
3182 confidence_intervals={}
3183 for rho,confidence_level,vol,bounds in b:
3184 acc_vol+=vol
3185
3186 if confidence_level>confidence_levels[cl_idx]:
3187 print str(confidence_level)
3188 print acc_vol
3189 confidence_intervals[confidence_levels[cl_idx]]=acc_vol
3190 cl_idx+=1
3191 if cl_idx==len(confidence_levels):
3192 break
3193
3194 acc_vol = 0.
3195 for rho,sample_number,vol,bounds in b:
3196 acc_vol+=vol
3197 print 'total area: ' + str(acc_vol)
3198
3199
3200 inj_confidence = None
3201 inj_confidence_area = None
3202 if inj_rho is not None:
3203 acc_vol=0.
3204 for rho,confidence_level,vol,bounds in b:
3205 acc_vol+=vol
3206 if rho <= inj_rho:
3207 inj_confidence = confidence_level
3208 inj_confidence_area = acc_vol
3209 injInfo.append(inj_confidence)
3210 injInfo.append(inj_confidence_area)
3211 print 'inj ' +str(vol)
3212 break
3213 return confidence_intervals, b, injInfo
3214
3215 -def kdtree_bin(posterior,coord_names,confidence_levels,initial_boundingbox = None,samples_per_bin = 10):
3216 """
3217 takes samples and applies a KDTree to them to return confidence levels
3218 returns confidence_intervals - dictionary of user_provided_CL:calculated_volume
3219 b - ordered list of KD leaves
3220 initial_boundingbox - list of lists [upperleft_coords,lowerright_coords]
3221 injInfo - if injection values provided then returns
3222 [Bounds_of_inj_kd_leaf ,number_samples_in_box, weight_of_box,injection_CL ,injection_CL_volume]
3223 Not quite sure that the repeated samples case is fixed, posibility of infinite loop.
3224 """
3225 confidence_levels.sort()
3226 print confidence_levels
3227 class Harvester(list):
3228 """
3229 when called by kdtree.operate will be used to calculate the density of each bin
3230 """
3231 def __init__(self):
3232 list.__init__(self)
3233 self.unrho=0.
3234
3235 def __call__(self,tree):
3236 number_density=float(len(tree.objects()))/float(tree.volume())
3237 self.append([number_density,len(tree.objects()),tree.volume(),tree.bounds()])
3238 self.unrho+=number_density
3239
3240 def close_ranks(self):
3241
3242 for i in range(len(self)):
3243 self[i][0]/=self.unrho
3244
3245 return sorted(self,key=itemgetter(0))
3246
3247 def h(a,b):
3248 pass
3249
3250 peparser=PEOutputParser('common')
3251
3252 samples,header=posterior.samples()
3253 header=header.split()
3254 coordinatized_samples=[PosteriorSample(row, header, coord_names) for row in samples]
3255
3256
3257 if initial_boundingbox is None:
3258 low=coordinatized_samples[0].coord()
3259 high=coordinatized_samples[0].coord()
3260 for obj in coordinatized_samples[1:]:
3261 low=np.minimum(low,obj.coord())
3262 high=np.maximum(high,obj.coord())
3263 initial_boundingbox = [low,high]
3264
3265 tree=KDTreeVolume(coordinatized_samples,initial_boundingbox)
3266
3267 a=Harvester()
3268 tree.operate(a,h,boxing=samples_per_bin)
3269
3270 b=a.close_ranks()
3271 b.reverse()
3272 totalSamples = len(tree.objects())
3273 samplecounter=0.0
3274 for entry in b:
3275 samplecounter += entry[1]
3276 entry[1] = float(samplecounter)/float(totalSamples)
3277
3278 acc_rho=0.
3279 acc_vol=0.
3280 cl_idx=0
3281
3282
3283 def checkNone(listoParams):
3284 for param in listoParams:
3285 if posterior[param].injval is None:
3286 return False
3287 return True
3288
3289
3290 if checkNone(coord_names):
3291 injBound,injNum,injWeight = tree.search([posterior[x].injval for x in coord_names],boxing = samples_per_bin)
3292 injInfo = [injBound,injNum,injWeight]
3293
3294 inj_volume = 1.
3295 low = injBound[0]
3296 high = injBound[1]
3297 for aCoord,bCoord in zip(low,high):
3298 inj_volume = inj_volume*(bCoord - aCoord)
3299 inj_number_density=float(injNum)/float(inj_volume)
3300 inj_rho = inj_number_density / a.unrho
3301 print injNum,inj_volume,inj_number_density,a.unrho,injBound
3302 else:
3303 injInfo = None
3304 inj_area = None
3305 inj_number_density=None
3306 inj_rho = None
3307
3308
3309 confidence_intervals={}
3310 for rho,sample_number,vol,bounds in b:
3311 acc_vol+=vol
3312
3313 if sample_number>confidence_levels[cl_idx]:
3314 confidence_intervals[confidence_levels[cl_idx]]=(acc_vol,sample_number)
3315 cl_idx+=1
3316 if cl_idx==len(confidence_levels):
3317 break
3318
3319 acc_vol = 0.
3320 for rho,sample_number,vol,bounds in b:
3321 acc_vol+=vol
3322
3323
3324 inj_confidence = None
3325 inj_confidence_area = None
3326 if inj_rho is not None:
3327 print 'calculating cl'
3328 acc_vol=0.
3329 for rho,confidence_level,vol,bounds in b:
3330 acc_vol+=vol
3331 if rho <= inj_rho:
3332 inj_confidence = confidence_level
3333 inj_confidence_area = acc_vol
3334 injInfo.append(inj_confidence)
3335 injInfo.append(inj_confidence_area)
3336 break
3337
3338 return confidence_intervals, b, initial_boundingbox,injInfo
3339
3340 -def kdtree_bin2Step(posterior,coord_names,confidence_levels,initial_boundingbox = None,samples_per_bin = 10,injCoords = None,alternate = False, fraction = 0.5, skyCoords=False):
3341 """
3342 input: posterior class instance, list of confidence levels, optional choice of inital parameter space, samples per box in kdtree
3343 note initial_boundingbox is [[lowerbound of each param][upper bound of each param]], if not specified will just take limits of samples
3344 fraction is proportion of samples used for making the tree structure.
3345 returns: confidence_intervals, sorted list of kd objects, initial_boundingbox, injInfo
3346 where injInfo is [bounding box injection is found within, samples in said box, weighting of box (in case of repeated samples),inj_confidence, inj_confidence_area]
3347 """
3348 confidence_levels.sort()
3349
3350 samples,header=posterior.samples()
3351 numberSamples = len(samples)
3352 if alternate == False:
3353 samplesStructure, samplesFill = random_split(samples,fraction)
3354 else:
3355 samplesStructure = samples[:int(numberSamples*fraction)]
3356 samplesFill = samples[int(numberSamples*fraction):]
3357 samplesFillLen = len(samplesFill)
3358
3359 header=header.split()
3360 coordinatized_samples=[PosteriorSample(row, header, coord_names) for row in samplesStructure]
3361
3362 if skyCoords == True:
3363 initial_boundingbox = [[0,-pi_constant/2.],[2*pi_constant,pi_constant/2.]]
3364 if initial_boundingbox is None:
3365 low=coordinatized_samples[0].coord()
3366 high=coordinatized_samples[0].coord()
3367 for obj in coordinatized_samples[1:]:
3368 low=np.minimum(low,obj.coord())
3369 high=np.maximum(high,obj.coord())
3370 initial_boundingbox = [low,high]
3371 tree=KDTreeVolume(coordinatized_samples,initial_boundingbox)
3372 tree2fill = tree.fillNewTree(boxing=samples_per_bin, isArea = skyCoords)
3373
3374 columns = []
3375 for name in coord_names:
3376 columns.append(header.index(name))
3377
3378 for sample in samplesFill:
3379 tempSample=[]
3380 for column in columns:
3381 tempSample.append(sample[column])
3382 addSample(tree2fill,tempSample)
3383
3384 def getValues(tree,listing):
3385 if tree._left is None:
3386 listing.append([tree.bounds(),tree._importance,tree._samples,tree._volume])
3387 else:
3388 getValues(tree._left,listing)
3389 getValues(tree._right,listing)
3390
3391 listLeaves = []
3392 getValues(tree2fill,listLeaves)
3393
3394 clSamples = []
3395 for cl in confidence_levels:
3396 clSamples.append(samplesFillLen*cl)
3397
3398 sortedLeavesList = sorted(listLeaves, key=lambda importance: importance[1])
3399 sortedLeavesList.reverse()
3400 runningTotalSamples = 0
3401 for i in range(len(sortedLeavesList)):
3402 runningTotalSamples += sortedLeavesList[i][2]
3403 sortedLeavesList[i].append(float(runningTotalSamples)/samplesFillLen,)
3404
3405
3406 level = 0
3407 countSamples = 0
3408 volume = 0
3409 lencl = len(clSamples)
3410
3411 confidence_intervals={}
3412 interpConfAreas = {}
3413 countLeaves = 0
3414 for leaf in sortedLeavesList:
3415 countSamples += leaf[2]
3416 countLeaves += 1
3417 volume += leaf[3]
3418 if level < lencl and countSamples >= clSamples[level]:
3419 confidence_intervals[confidence_levels[level]]=(volume,float(countSamples)/samplesFillLen)
3420 interpConfAreas[confidence_levels[level]] = volume-leaf[3]*(countSamples-clSamples[level])/leaf[2]
3421 level +=1
3422
3423 if injCoords is not None:
3424 injBound,injNum,injImportance = tree2fill.search(injCoords)
3425 injInfo = [injBound,injNum,injImportance]
3426 else:
3427 injInfo = None
3428
3429
3430
3431 inj_confidence = None
3432 inj_confidence_area = None
3433 if injInfo is not None:
3434 acc_vol=0.
3435 acc_cl=0.
3436 for leaf in sortedLeavesList:
3437 acc_vol+=leaf[3]
3438 acc_cl+=leaf[2]
3439 if leaf[1] <= injImportance:
3440 inj_confidence = float(acc_cl)/samplesFillLen
3441 inj_confidence_area = acc_vol
3442 injInfo.append(inj_confidence)
3443 injInfo.append(inj_confidence_area)
3444 break
3445
3446 return sortedLeavesList, interpConfAreas, injInfo
3447
3448 def checkNone(listoParams):
3449 for param in listoParams:
3450 if posterior[param].injval is None:
3451 return False
3452 return True
3453
3455 """
3456 Determine the 2-parameter Bayesian Confidence Intervals using a greedy
3457 binning algorithm.
3458
3459 @param posterior: an instance of the Posterior class.
3460
3461 @param greedy2Params: a dict - {param1Name:param1binSize,param2Name:param2binSize} .
3462
3463 @param confidence_levels: A list of floats of the required confidence intervals [(0-1)].
3464 """
3465
3466
3467 par1_name,par2_name=greedy2Params.keys()
3468
3469
3470 par1pos=posterior[par1_name.lower()].samples
3471 par2pos=posterior[par2_name.lower()].samples
3472
3473
3474 par1_bin=greedy2Params[par1_name]
3475 par2_bin=greedy2Params[par2_name]
3476
3477
3478 par1_injvalue=posterior[par1_name.lower()].injval
3479 par2_injvalue=posterior[par2_name.lower()].injval
3480
3481
3482 par1pos_min=min(par1pos)[0]
3483 par2pos_min=min(par2pos)[0]
3484
3485 par1pos_max=max(par1pos)[0]
3486 par2pos_max=max(par2pos)[0]
3487
3488 par1pos_Nbins= int(ceil((par1pos_max - par1pos_min)/par1_bin))+1
3489
3490 par2pos_Nbins= int(ceil((par2pos_max - par2pos_min)/par2_bin))+1
3491
3492 greedyHist = np.zeros(par1pos_Nbins*par2pos_Nbins,dtype='i8')
3493 greedyPoints = np.zeros((par1pos_Nbins*par2pos_Nbins,2))
3494
3495
3496 par1_point=par1pos_min
3497 par2_point=par2pos_min
3498 for i in range(par2pos_Nbins):
3499
3500 par1_point=par1pos_min
3501 for j in range(par1pos_Nbins):
3502
3503 greedyPoints[j+par1pos_Nbins*i,0]=par1_point
3504 greedyPoints[j+par1pos_Nbins*i,1]=par2_point
3505 par1_point+=par1_bin
3506 par2_point+=par2_bin
3507
3508
3509
3510 injbin=None
3511 if par1_injvalue is not None and par2_injvalue is not None:
3512
3513 par1_binNumber=int(floor((par1_injvalue-par1pos_min)/par1_bin))
3514 par2_binNumber=int(floor((par2_injvalue-par2pos_min)/par2_bin))
3515
3516 injbin=int(par1_binNumber+par2_binNumber*par1pos_Nbins)
3517 elif par1_injvalue is None and par2_injvalue is not None:
3518 print "Injection value not found for %s!"%par1_name
3519
3520 elif par1_injvalue is not None and par2_injvalue is None:
3521 print "Injection value not found for %s!"%par2_name
3522
3523
3524 for par1_samp,par2_samp in zip(par1pos,par2pos):
3525 par1_samp=par1_samp[0]
3526 par2_samp=par2_samp[0]
3527 par1_binNumber=int(floor((par1_samp-par1pos_min)/par1_bin))
3528 par2_binNumber=int(floor((par2_samp-par2pos_min)/par2_bin))
3529 try:
3530 greedyHist[par1_binNumber+par2_binNumber*par1pos_Nbins]+=1
3531 except:
3532 raise RuntimeError("Problem binning samples: %i,%i,%i,%i,%i,%f,%f,%f,%f,%f,%f .")%(par1_binNumber,par2_binNumber,par1pos_Nbins,par2pos_Nbins,par1_binNumber+par2_binNumber*par1pos_Nbins,par1_samp,par1pos_min,par1_bin,par1_samp,par2pos_min,par2_bin)
3533
3534 toppoints,injection_cl,reses,injection_area=\
3535 _greedy_bin(
3536 greedyHist,
3537 greedyPoints,
3538 injbin,
3539 float(par1_bin*par2_bin),
3540 int(len(par1pos)),
3541 confidence_levels
3542 )
3543
3544 return toppoints,injection_cl,reses,injection_area
3545
3547 """
3548 Utility function to convert longitude,latitude on a unit sphere to
3549 cartesian co-ordinates.
3550 """
3551
3552 x=np.cos(lat)*np.cos(long)
3553 y=np.cos(lat)*np.sin(long)
3554 z=np.sin(lat)
3555 return np.array([x,y,z])
3556
3559 """
3560 Utiltiy function to convert r,theta,phi to cartesian co-ordinates.
3561 """
3562 x = r*np.sin(theta)*np.cos(phi)
3563 y = r*np.sin(theta)*np.sin(phi)
3564 z = r*np.cos(theta)
3565 return x,y,z
3566
3569 """
3570 Utility function to convert cartesian coords to r,theta,phi.
3571 """
3572 r = np.sqrt(x*x + y*y + z*z)
3573 theta = np.arccos(z/r)
3574 phi = np.fmod(2*pi_constant + np.arctan2(y,x), 2*pi_constant)
3575
3576 return r,theta,phi
3577
3580 """
3581 Greedy bins the sky posterior samples into a grid on the sky constructed so that
3582 sky boxes have roughly equal size (determined by skyres).
3583
3584 @param posterior: Posterior class instance containing ra and dec samples.
3585
3586 @param skyres: Desired approximate size of sky pixel on one side.
3587
3588 @param confidence_levels: List of desired confidence levels [(0-1)].
3589 """
3590
3591 from pylal import skylocutils
3592
3593 np.seterr(under='ignore')
3594
3595 if 'ra' in posterior.names and 'dec' in posterior.names:
3596 skypos=np.column_stack([posterior['ra'].samples,posterior['dec'].samples])
3597 raname='ra'
3598 decname='dec'
3599 elif 'rightascension' in posterior.names and 'declination' in posterior.names:
3600 skypos=np.column_stack([posterior['rightascension'].samples,posterior['declination'].samples])
3601 raname='rightascension'
3602 decname='declination'
3603 else:
3604 raise RuntimeError('could not find ra and dec or rightascention and declination in column names for sky position')
3605
3606 injvalues=None
3607
3608 sky_injpoint=(posterior[raname].injval,posterior[decname].injval)
3609
3610 skypoints=np.array(skylocutils.gridsky(float(skyres)))
3611 skycarts=map(lambda s: pol2cart(s[1],s[0]),skypoints)
3612 skyinjectionconfidence=None
3613
3614 shist=_skyhist_cart(np.array(skycarts),skypos)
3615
3616
3617 bins=skycarts
3618
3619
3620 injbin=None
3621 if None not in sky_injpoint:
3622 injhist=_skyhist_cart_slow(skycarts,np.array([sky_injpoint]))
3623 injbin=injhist.tolist().index(1)
3624 print 'Found injection in bin %d with co-ordinates %f,%f .'%(
3625 injbin,
3626 skypoints[injbin,0],
3627 skypoints[injbin,1]
3628 )
3629
3630 return _greedy_bin(shist,skypoints,injbin,float(skyres)*float(skyres),len(skypos),confidence_levels)
3631
3633 """Plots a sky map from a healpix map, optionally including an
3634 injected position.
3635
3636 :param hpmap: An array representing a healpix map (in nested
3637 ordering if ``nest = True``).
3638
3639 :param outdir: The output directory.
3640
3641 :param inj: If not ``None``, then ``[ra, dec]`` of the injection
3642 associated with the posterior map.
3643
3644 :param nest: Flag indicating the pixel ordering in healpix.
3645
3646 """
3647 import lalinference.plot as lp
3648
3649 fig = plt.figure(frameon=False, figsize=(8,6))
3650 ax = plt.subplot(111, projection='astro hours mollweide')
3651 ax.cla()
3652 ax.grid()
3653 lp.healpix_heatmap(hpmap, nest=nest, vmin=0.0, vmax=np.max(hpmap), cmap=plt.get_cmap('cylon'))
3654
3655 if inj is not None:
3656 ax.plot(inj[0], inj[1], '*', markerfacecolor='white', markeredgecolor='black', markersize=10)
3657
3658 plt.savefig(os.path.join(outdir, 'skymap.png'))
3659
3660 return fig
3661
3663 """Returns the area (in square degrees) for each confidence level with
3664 a greedy binning algorithm for the given healpix map.
3665
3666 """
3667
3668 hpmap = hpmap / np.sum(hpmap)
3669
3670 hpmap = np.sort(hpmap)[::-1]
3671 cum_hpmap = np.cumsum(hpmap)
3672
3673 pixarea = hp.nside2pixarea(hp.npix2nside(hpmap.shape[0]))
3674 pixarea = pixarea*(180.0/np.pi)**2
3675
3676 areas = []
3677 for cl in cls:
3678 npix = np.sum(cum_hpmap < cl)
3679 areas.append(npix*pixarea)
3680
3681 return np.array(areas)
3682
3684 """Returns the greedy p-value estimate for the given injection.
3685
3686 """
3687
3688 nside = hp.npix2nside(hpmap.shape[0])
3689 hpmap = hpmap / np.sum(hpmap)
3690
3691 injpix = hp.ang2pix(nside, np.pi/2.0-inj[1], inj[0], nest=nest)
3692 injvalue = hpmap[injpix]
3693
3694 return np.sum(hpmap[hpmap >= injvalue])
3695
3696
3697
3698 -def mc2ms(mc,eta):
3699 """
3700 Utility function for converting mchirp,eta to component masses. The
3701 masses are defined so that m1>m2. The rvalue is a tuple (m1,m2).
3702 """
3703 root = np.sqrt(0.25-eta)
3704 fraction = (0.5+root) / (0.5-root)
3705 invfraction = 1/fraction
3706
3707 m2= mc * np.power((1+fraction),0.2) / np.power(fraction,0.6)
3708
3709 m1= mc* np.power(1+invfraction,0.2) / np.power(invfraction,0.6)
3710 return (m1,m2)
3711
3712
3713
3714 -def q2ms(mc,q):
3715 """
3716 Utility function for converting mchirp,q to component masses. The
3717 masses are defined so that m1>m2. The rvalue is a tuple (m1,m2).
3718 """
3719 factor = mc * np.power(1+q, 1.0/5.0);
3720 m1 = factor * np.power(q, -3.0/5.0);
3721 m2 = factor * np.power(q, 2.0/5.0);
3722 return (m1,m2)
3723
3724
3725
3726 -def q2eta(mc,q):
3727 """
3728 Utility function for converting mchirp,q to eta. The
3729 rvalue is eta.
3730 """
3731 m1,m2 = q2ms(mc,q)
3732 eta = m1*m2/( (m1+m2)*(m1+m2) )
3733 return eta
3734
3735
3736
3737 -def mc2q(mc,eta):
3738 """
3739 Utility function for converting mchirp,eta to new mass ratio q (m2/m1).
3740 """
3741 m1,m2 = mc2ms(mc,eta)
3742 q = m2/m1
3743 return q
3744
3745
3746
3747 -def ang_dist(long1,lat1,long2,lat2):
3748 """
3749 Find the angular separation of (long1,lat1) and (long2,lat2), which are
3750 specified in radians.
3751 """
3752
3753 x1=np.cos(lat1)*np.cos(long1)
3754 y1=np.cos(lat1)*np.sin(long1)
3755 z1=np.sin(lat1)
3756 x2=np.cos(lat2)*np.cos(long2)
3757 y2=np.cos(lat2)*np.sin(long2)
3758 z2=np.sin(lat2)
3759 sep=math.acos(x1*x2+y1*y2+z1*z2)
3760 return(sep)
3761
3765 """
3766 Calculate dot products between vectors in rows of numpy arrays.
3767 """
3768 if vec1.ndim==1:
3769 product = (vec1*vec2).sum()
3770 else:
3771 product = (vec1*vec2).sum(axis=1).reshape(-1,1)
3772 return product
3773
3777 """
3778 Find angles between vectors in rows of numpy arrays.
3779 """
3780 vec1_mag = np.sqrt(array_dot(vec1, vec1))
3781 vec2_mag = np.sqrt(array_dot(vec2, vec2))
3782 return np.arccos(array_dot(vec1, vec2)/(vec1_mag*vec2_mag))
3783
3787 """
3788 Find polar angles of vectors in rows of a numpy array.
3789 """
3790 if vec.ndim==1:
3791 z = vec[2]
3792 else:
3793 z = vec[:,2].reshape(-1,1)
3794 norm = np.sqrt(array_dot(vec,vec))
3795 return np.arccos(z/norm)
3796
3800 """
3801 Compute general rotation matrices for a given angles and direction vectors.
3802 """
3803 cosa = np.cos(angle)
3804 sina = np.sin(angle)
3805 direction /= np.sqrt(array_dot(direction,direction))
3806
3807 try:
3808 nSamps = len(angle)
3809 R = np.array( [np.diag([i,i,i]) for i in cosa.flat] )
3810 R += np.array( [np.outer(direction[i],direction[i])*(1.0-cosa[i]) for i in range(nSamps)] )
3811 R += np.array( [np.array( [[ 0.0, -direction[i,2], direction[i,1]],
3812 [ direction[i,2], 0.0, -direction[i,0]],
3813 [-direction[i,1], direction[i,0], 0.0 ]] ) * sina[i] for i in range(nSamps)] )
3814
3815 except TypeError:
3816 R = np.diag([cosa,cosa,cosa])
3817 R += np.outer(direction,direction) * (1.0 - cosa)
3818 R += np.array( [[ 0.0, -direction[2], direction[1]],
3819 [ direction[2], 0.0, -direction[0]],
3820 [-direction[1], direction[0], 0.0 ]] ) * sina
3821 return R
3822
3826 """
3827 Rotate vectors using the given rotation matrices.
3828 """
3829 if vec.ndim == 1:
3830 newVec = np.dot(R[i],vec[i])
3831 else:
3832 newVec = np.array( [np.dot(R[i],vec[i]) for i in range(len(vec))] )
3833 return newVec
3834
3835
3836 -def ROTATEZ(angle, vx, vy, vz):
3837
3838 tmp1 = vx*np.cos(angle) - vy*np.sin(angle);
3839 tmp2 = vx*np.sin(angle) + vy*np.cos(angle);
3840 return np.asarray([tmp1,tmp2,vz])
3841
3843
3844 tmp1 = vx*np.cos(angle) + vz*np.sin(angle);
3845 tmp2 = - vx*np.sin(angle) + vz*np.cos(angle);
3846 return np.asarray([tmp1,vy,tmp2])
3847
3849 """
3850 Calculate orbital angular momentum vector.
3851 Note: The units of Lmag are different than what used in lalsimulation.
3852 Mc must be called in units of Msun here.
3853
3854 Note that if one wants to build J=L+S1+S2 with L returned by this function, S1 and S2
3855 must not get the Msun^2 factor.
3856 """
3857 eta = m1*m2/( (m1+m2)*(m1+m2) )
3858 Lmag = orbital_momentum_mag(fref, m1,m2,eta)
3859 Lx, Ly, Lz = sph2cart(Lmag, inclination, 0.0)
3860 return np.hstack((Lx,Ly,Lz))
3861
3864 v0 = np.power((m1+m2) *pi_constant * lal.MTSUN_SI * fref, 1.0/3.0)
3865
3866 PNFirst = (((m1+m2)**2)*eta)/v0
3867 PNSecond = 1+ (v0**2) * (3.0/2.0 +eta/6.0)
3868 Lmag= PNFirst*PNSecond
3869 return Lmag
3870
3872 """
3873 Calculate BH angular momentum vector.
3874 """
3875 Sx, Sy, Sz = sph2cart(m**2 * a, theta, phi)
3876 return np.hstack((Sx,Sy,Sz))
3877
3881 """
3882 Calculate best tidal parameters
3883 """
3884 lam_tilde = (8./13.)*((1.+7.*eta-31.*eta*eta)*(lambda1+lambda2) + np.sqrt(1.-4.*eta)*(1.+9.*eta-11.*eta*eta)*(lambda1-lambda2))
3885 dlam_tilde = (1./2.)*(np.sqrt(1.-4.*eta)*(1.-13272.*eta/1319.+8944.*eta*eta/1319.)*(lambda1+lambda2) + (1.-15910.*eta/1319.+32850.*eta*eta/1319.+3380.*eta*eta*eta/1319.)*(lambda1-lambda2))
3886 return lam_tilde, dlam_tilde
3887
3888 -def spin_angles(fref,mc,eta,incl,a1,theta1,phi1,a2=None,theta2=None,phi2=None):
3889 """
3890 Calculate physical spin angles.
3891 """
3892 singleSpin = None in (a2,theta2,phi2)
3893 m1, m2 = mc2ms(mc,eta)
3894 L = orbital_momentum(fref, m1,m2, incl)
3895 S1 = component_momentum(m1, a1, theta1, phi1)
3896 if not singleSpin:
3897 S2 = component_momentum(m2, a2, theta2, phi2)
3898 else:
3899 S2 = 0.0
3900 J = L + S1 + S2
3901 tilt1 = array_ang_sep(L,S1)
3902 if not singleSpin:
3903 tilt2 = array_ang_sep(L,S2)
3904 else:
3905 tilt2 = None
3906 theta_jn = array_polar_ang(J)
3907 beta = array_ang_sep(J,L)
3908 return tilt1, tilt2, theta_jn, beta
3909
3911 """
3912 Calculate the magnitude of the effective precessing spin
3913 following convention from Phys. Rev. D 91, 024043 -- arXiv:1408.1810
3914 note: the paper uses naming convention where m1 < m2
3915 (and similar for associated spin parameters) and q > 1
3916 """
3917 q_inv = m1/m2
3918 A1 = 2. + (3.*q_inv/2.)
3919 A2 = 2. + 3./(2.*q_inv)
3920 S1_perp = a1*np.sin(tilt1)*m1*m1
3921 S2_perp = a2*np.sin(tilt2)*m2*m2
3922 Sp = np.maximum(A1*S2_perp, A2*S1_perp)
3923 chi_p = Sp/(A2*m1*m1)
3924 return chi_p
3925
3927 """
3928 Calculate the redshift from the luminosity distance measurement using the
3929 Cosmology Calculator provided in LAL.
3930 By default assuming cosmological parameters from arXiv:1502.01589 - 'Planck 2015 results. XIII. Cosmological parameters'
3931 Using parameters from table 4, column 'TT+lowP+lensing+ext'
3932 This corresponds to Omega_M = 0.3065, Omega_Lambda = 0.6935, H_0 = 67.90 km s^-1 Mpc^-1
3933 Returns an array of redshifts
3934 """
3935 def find_z_root(z,dl,omega):
3936 return dl - lal.LuminosityDistance(omega,z)
3937
3938 omega = lal.CreateCosmologicalParameters(h,om,ol,w0,0.0,0.0)
3939 if isinstance(distance,float):
3940 z = np.array([newton(find_z_root,np.random.uniform(0.0,2.0),args = (distance,omega))])
3941 else:
3942 z = np.array([newton(find_z_root,np.random.uniform(0.0,2.0),args = (d,omega)) for d in distance[:,0]])
3943 return z.reshape(z.shape[0],1)
3944
3946 """
3947 Calculate source mass parameter for mass m as:
3948 m_source = m / (1.0 + z)
3949 For a parameter m.
3950 """
3951 return mass / (1.0 + redshift)
3952
3954 """
3955 Wrapper function for SimInspiralTransformPrecessingNewInitialConditions().
3956 Vectorizes function for use in append_mapping() methods of the posterior class.
3957 """
3958 try:
3959 import lalsimulation as lalsim
3960 except ImportError:
3961 print 'bayespputils.py: Cannot import lalsimulation SWIG bindings to calculate physical to radiation'
3962 print 'frame angles, did you remember to use --enable-swig-python when ./configuring lalsimulation?'
3963 return None
3964 from numpy import shape
3965 transformFunc = lalsim.SimInspiralTransformPrecessingNewInitialConditions
3966
3967
3968 m1_SI = m1*lal.MSUN_SI
3969 m2_SI = m2*lal.MSUN_SI
3970
3971
3972 ins = [theta_jn, phi_jl, tilt1, tilt2, phi12, a1, a2, m1_SI, m2_SI, fref]
3973 if len(shape(ins))>1:
3974
3975 try:
3976 for p,param in enumerate(ins):
3977 ins[p] = param.flatten()
3978 except:
3979 pass
3980
3981 try:
3982 results = np.array([transformFunc(t_jn, p_jl, t1, t2, p12, a1, a2, m1_SI, m2_SI, f) for (t_jn, p_jl, t1, t2, p12, a1, a2, m1_SI, m2_SI, f) in zip(*ins)])
3983 iota = results[:,0].reshape(-1,1)
3984 spin1x = results[:,1].reshape(-1,1)
3985 spin1y = results[:,2].reshape(-1,1)
3986 spin1z = results[:,3].reshape(-1,1)
3987 spin2x = results[:,4].reshape(-1,1)
3988 spin2y = results[:,5].reshape(-1,1)
3989 spin2z = results[:,6].reshape(-1,1)
3990 a1,theta1,phi1 = cart2sph(spin1x,spin1y,spin1z)
3991 a2,theta2,phi2 = cart2sph(spin2x,spin2y,spin2z)
3992
3993 mc = np.power(m1*m2,3./5.)*np.power(m1+m2,-1./5.)
3994 L = orbital_momentum(fref, m1,m2, iota)
3995 S1 = np.hstack([m1*m1*spin1x,m1*m1*spin1y,m1*m1*spin1z])
3996 S2 = np.hstack([m2*m2*spin2x,m2*m2*spin2y,m2*m2*spin2z])
3997 J = L + S1 + S2
3998 beta = array_ang_sep(J,L)
3999
4000 return iota, theta1, phi1, theta2, phi2, beta
4001 except:
4002
4003 return None
4004
4005 elif len(shape(ins))<=1:
4006
4007 try:
4008 for p,param in enumerate(ins):
4009 ins[p] = param
4010 except:
4011 pass
4012
4013 try:
4014 results = np.array(transformFunc(theta_jn, phi_jl, tilt1, tilt2, phi12, a1, a2, m1_SI, m2_SI, fref))
4015 iota = results[0]
4016 spin1x = results[1]
4017 spin1y = results[2]
4018 spin1z = results[3]
4019 spin2x = results[4]
4020 spin2y = results[5]
4021 spin2z = results[6]
4022 a1,theta1,phi1 = cart2sph(spin1x,spin1y,spin1z)
4023 a2,theta2,phi2 = cart2sph(spin2x,spin2y,spin2z)
4024
4025 mc = np.power(m1*m2,3./5.)*np.power(m1+m2,-1./5.)
4026 L = orbital_momentum(fref, m1,m2, iota)
4027 S1 = m1*m1*np.hstack([spin1x,spin1y,spin1z])
4028 S2 = m2*m2*np.hstack([spin2x,spin2y,spin2z])
4029 J = L + S1 + S2
4030 beta = array_ang_sep(J,L)
4031
4032 return iota, theta1, phi1, theta2, phi2, beta
4033
4034 except:
4035
4036 return None
4037
4041
4042 from scipy import seterr as sp_seterr
4043
4044 np.seterr(under='ignore')
4045 sp_seterr(under='ignore')
4046 pos_samps=onedpos.samples
4047 try:
4048 gkde=onedpos.gaussian_kde
4049 except np.linalg.linalg.LinAlgError:
4050 print 'Singular matrix in KDE. Skipping'
4051 else:
4052 ind=np.linspace(np.min(pos_samps),np.max(pos_samps),101)
4053 kdepdf=gkde.evaluate(ind)
4054 plt.plot(ind,kdepdf,color='green')
4055 return
4056
4058 plt.hist(pos_samps,kdepdf)
4059
4060
4061 -def plot_one_param_pdf(posterior,plot1DParams,analyticPDF=None,analyticCDF=None,plotkde=False):
4062 """
4063 Plots a 1D histogram and (gaussian) kernel density estimate of the
4064 distribution of posterior samples for a given parameter.
4065
4066 @param posterior: an instance of the Posterior class.
4067
4068 @param plot1DParams: a dict; {paramName:Nbins}
4069
4070 @param analyticPDF: an analytic probability distribution function describing the distribution.
4071
4072 @param analyticCDF: an analytic cumulative distribution function describing the distribution.
4073
4074 """
4075
4076
4077
4078 param=plot1DParams.keys()[0].lower()
4079 histbins=plot1DParams.values()[0]
4080
4081 pos_samps=posterior[param].samples
4082 injpar=posterior[param].injval
4083 trigvals=posterior[param].trigvals
4084
4085
4086 myfig=plt.figure(figsize=(6,4.5),dpi=150)
4087
4088 axes=plt.Axes(myfig,[0.15,0.15,0.6,0.76])
4089 myfig.add_axes(axes)
4090 majorFormatterX=ScalarFormatter(useMathText=True)
4091 majorFormatterX.format_data=lambda data:'%.6g'%(data)
4092 majorFormatterY=ScalarFormatter(useMathText=True)
4093 majorFormatterY.format_data=lambda data:'%.6g'%(data)
4094 majorFormatterX.set_scientific(True)
4095 majorFormatterY.set_scientific(True)
4096 offset=0.0
4097 if param.find('time')!=-1:
4098 offset=floor(min(pos_samps))
4099 pos_samps=pos_samps-offset
4100 if injpar is not None:
4101 injpar=injpar-offset
4102 ax1_name=param+' + %i'%(int(offset))
4103 else: ax1_name=param
4104
4105 (n, bins, patches)=plt.hist(pos_samps,histbins,normed='true',facecolor='grey')
4106 Nchars=max(map(lambda d:len(majorFormatterX.format_data(d)),axes.get_xticks()))
4107 if Nchars>8:
4108 Nticks=3
4109 elif Nchars>5:
4110 Nticks=4
4111 elif Nchars>4:
4112 Nticks=6
4113 else:
4114 Nticks=6
4115 locatorX=matplotlib.ticker.MaxNLocator(nbins=Nticks)
4116 xmin,xmax=plt.xlim()
4117 if param=='rightascension' or param=='ra':
4118 locatorX=RALocator(min=xmin,max=xmax)
4119 majorFormatterX=RAFormatter()
4120 if param=='declination' or param=='dec':
4121 locatorX=DecLocator(min=xmin,max=xmax)
4122 majorFormatterX=DecFormatter()
4123 axes.xaxis.set_major_formatter(majorFormatterX)
4124 axes.yaxis.set_major_formatter(majorFormatterY)
4125
4126 locatorX.view_limits(bins[0],bins[-1])
4127 axes.xaxis.set_major_locator(locatorX)
4128 if plotkde: plot_one_param_pdf_kde(myfig,posterior[param])
4129 histbinSize=bins[1]-bins[0]
4130 if analyticPDF:
4131 (xmin,xmax)=plt.xlim()
4132 x = np.linspace(xmin,xmax,2*len(bins))
4133 plt.plot(x, analyticPDF(x+offset), color='r', linewidth=2, linestyle='dashed')
4134 if analyticCDF:
4135
4136 max_samps=1000
4137 from numpy.random import permutation
4138 samps = permutation(pos_samps)[:max_samps,:].flatten()
4139 D,p = stats.kstest(samps+offset, analyticCDF, mode='asymp')
4140 plt.title("%s: ks p-value %.3f"%(param,p))
4141
4142 rbins=None
4143
4144 if injpar is not None:
4145
4146 delta_samps=max(pos_samps)-min(pos_samps)
4147 minrange=min(pos_samps)-0.05*delta_samps
4148 maxrange=max(pos_samps)+0.05*delta_samps
4149 if minrange<injpar and maxrange>injpar:
4150
4151 plt.axvline(injpar, color='r', linestyle='-.', linewidth=4)
4152
4153
4154
4155
4156
4157 if min(pos_samps)<injpar and max(pos_samps)>injpar:
4158 bins_to_inj=(injpar-bins[0])/histbinSize
4159 injbinh=int(floor(bins_to_inj))
4160 injbin_frac=bins_to_inj-float(injbinh)
4161
4162
4163 rbins=(sum(n[0:injbinh-1])+injbin_frac*n[injbinh])*histbinSize
4164
4165 if trigvals is not None:
4166 for IFO in [IFO for IFO in trigvals.keys()]:
4167 trigval = trigvals[IFO]
4168 if min(pos_samps)<trigval and max(pos_samps)>trigval:
4169 if IFO=='H1': color = 'r'
4170 elif IFO=='L1': color = 'g'
4171 elif IFO=='V1': color = 'm'
4172 else: color = 'c'
4173 plt.axvline(trigval, color=color, linestyle='-.')
4174
4175 plt.grid()
4176 plt.xlabel(plot_label(ax1_name))
4177 plt.ylabel('Probability Density')
4178
4179
4180 if(param.lower()=='ra' or param.lower()=='rightascension'):
4181 xmin,xmax=plt.xlim()
4182 plt.xlim(xmax,xmin)
4183
4184
4185
4186
4187
4188
4189
4190
4191
4192 return rbins,myfig
4193
4194
4195 -class RALocator(matplotlib.ticker.MultipleLocator):
4196 """
4197 RA tick locations with some intelligence
4198 """
4199 - def __init__(self,min=0.0,max=2.0*pi_constant):
4200 hour=pi_constant/12.0
4201 if(max-min)>12.0*hour:
4202 base=3.0*hour
4203 elif(max-min)>6.0*hour:
4204 base=2.0*hour
4205
4206 elif (max-min)>3.0*pi_constant/12.0:
4207 base=hour
4208 elif (max-min)>hour:
4209 base=hour/2.0
4210 else:
4211 base=hour/4.0
4212
4213 matplotlib.ticker.MultipleLocator.__init__(self,base=base)
4214
4215 -class DecLocator(matplotlib.ticker.MultipleLocator):
4216 """
4217 Dec tick locations with some intelligence
4218 """
4219 - def __init__(self, min=-pi_constant/2.0,max=pi_constant/2.0):
4220 deg=pi_constant/180.0
4221 if (max-min)>60*deg:
4222 base=30.0*deg
4223 elif (max-min)>20*deg:
4224 base=10*deg
4225 elif (max-min)>10*deg:
4226 base=5*deg
4227 else:
4228 base=deg
4229 matplotlib.ticker.MultipleLocator.__init__(self,base=base)
4230
4234
4238
4266
4290
4292 """
4293 round given angle in radians to integer hours, degrees, mins or secs
4294 accuracy can be 'hour'. 'deg', 'min', 'sec', 'all', all does nothing
4295 'arcmin', 'arcsec'
4296 """
4297 if accuracy=='all': return locs
4298 if accuracy=='hour': mult=24
4299 if accuracy=='deg': mult=360
4300 if accuracy=='min': mult=24*60
4301 if accuracy=='sec': mult=24*60*60
4302 if accuracy=='arcmin': mult=360*60
4303 if accuracy=='arcsec': mult=360*60*60
4304 mult=mult/(2.0*pi_constant)
4305 return round(rads*mult)/mult
4306
4308 secs=radians*12.0*3600/pi_constant
4309 hours, rem = divmod(secs, 3600 )
4310 mins,rem = divmod(rem, 60 )
4311 secs = rem
4312 if secs>=59.5:
4313 secs=secs-60
4314 mins=mins+1
4315 if mins>=59.5:
4316 mins=mins-60
4317 hours=hours+1
4318 if accuracy=='hour': return ur'%ih'%(hours)
4319 if accuracy=='min': return ur'%ih%im'%(hours,mins)
4320 if accuracy=='sec': return ur'%ih%im%2.0fs'%(hours,mins,secs)
4321 else:
4322 if abs(fmod(secs,60.0))>=0.5: return(getRAString(radians,accuracy='sec'))
4323 if abs(fmod(mins,60.0))>=0.5: return(getRAString(radians,accuracy='min'))
4324 else: return(getRAString(radians,accuracy='hour'))
4325
4327
4328 if matplotlib.rcParams['text.usetex']:
4329 degsymb='$^\circ$'
4330 minsymb="'"
4331 secsymb="''"
4332 else:
4333 degsymb=u'\u00B0'
4334 minsymb=u'\u0027'
4335 secsymb=u'\u2033'
4336 if(radians<0):
4337 radians=-radians
4338 sign=-1
4339 else: sign=+1
4340 deg,rem=divmod(radians,pi_constant/180.0)
4341 mins, rem = divmod(rem, pi_constant/(180.0*60.0))
4342 secs = rem * (180.0*3600.0)/pi_constant
4343
4344
4345
4346
4347
4348
4349 if (accuracy=='arcmin' or accuracy=='deg') and secs>30: mins=mins+1
4350 if accuracy=='deg' and mins>30: deg=deg+1
4351 if accuracy=='deg': return ur'%i'%(sign*deg)+degsymb
4352 if accuracy=='arcmin': return ur'%i%s%i%s'%(sign*deg,degsymb,mins,minsymb)
4353 if accuracy=='arcsec': return ur'%i%s%i%s%2.0f%s'%(sign*deg,degsymb,mins,minsymb,secs,secsymb)
4354 else:
4355
4356
4357
4358 return(getDecString(sign*radians,accuracy='deg'))
4359
4361 """
4362 Make a corner plot using the triangle module
4363 (See http://github.com/dfm/corner.py)
4364 @param posterior: The Posterior object
4365 @param levels: a list of confidence levels
4366 @param parnames: list of parameters to include
4367 """
4368 try:
4369 import corner
4370 except ImportError:
4371 try:
4372 import triangle as corner
4373 except ImportError:
4374 print 'Cannot load corner module. Try running\n\t$ pip install corner'
4375 return None
4376 parnames=filter(lambda x: x in posterior.names, parnames)
4377 labels = [plot_label(parname) for parname in parnames]
4378 data = np.hstack([posterior[p].samples for p in parnames])
4379 if posterior.injection:
4380 injvals=[posterior[p].injval for p in parnames]
4381 myfig=corner.corner(data,labels=labels,truths=injvals,quantiles=levels,plot_datapoints=False,bins=20)
4382 else:
4383 myfig=corner.corner(data,labels=labels,quantiles=levels,plot_datapoints=False,bins=20)
4384 return(myfig)
4385
4386
4387 -def plot_two_param_kde_greedy_levels(posteriors_by_name,plot2DkdeParams,levels,colors_by_name,line_styles=__default_line_styles,figsize=(4,3),dpi=250,figposition=[0.2,0.2,0.48,0.75],legend='right',hatches_by_name=None,Npixels=50):
4388 """
4389 Plots a 2D kernel density estimate of the 2-parameter marginal posterior.
4390
4391 @param posterior: an instance of the Posterior class.
4392
4393 @param plot2DkdeParams: a dict {param1Name:Nparam1Bins,param2Name:Nparam2Bins}
4394 """
4395
4396 from scipy import seterr as sp_seterr
4397 confidence_levels=levels
4398
4399
4400
4401 par2_name,par1_name=plot2DkdeParams.keys()
4402 xbin=plot2DkdeParams[par1_name]
4403 ybin=plot2DkdeParams[par2_name]
4404 levels= levels
4405 np.seterr(under='ignore')
4406 sp_seterr(under='ignore')
4407
4408 fig=plt.figure(1,figsize=figsize,dpi=dpi)
4409 plt.clf()
4410 axes=fig.add_axes(figposition)
4411 name_list=[]
4412
4413
4414 if len(line_styles)<len(levels):
4415 raise RuntimeError("Error: Need as many or more line styles to choose from as confidence levels to plot!")
4416
4417 CSlst=[]
4418 for name,posterior in posteriors_by_name.items():
4419 print 'Plotting '+name
4420 name_list.append(name)
4421 par1_injvalue=posterior[par1_name].injval
4422 par2_injvalue=posterior[par2_name].injval
4423
4424 par_trigvalues1=posterior[par1_name].trigvals
4425 par_trigvalues2=posterior[par2_name].trigvals
4426 xdat=posterior[par1_name].samples
4427 ydat=posterior[par2_name].samples
4428 a=np.squeeze(posterior[par1_name].samples)
4429 b=np.squeeze(posterior[par2_name].samples)
4430 offset=0.0
4431 if par1_name.find('time')!=-1:
4432 offset=floor(min(a))
4433 a=a-offset
4434 if par1_injvalue:
4435 par1_injvalue=par1_injvalue-offset
4436 ax1_name=par1_name+' + %i'%(int(offset))
4437 else: ax1_name=par1_name
4438
4439 if par2_name.find('time')!=-1:
4440 offset=floor(min(b))
4441 b=b-offset
4442 if par2_injvalue:
4443 par2_injvalue=par2_injvalue-offset
4444 ax2_name=par2_name+' + %i'%(int(offset))
4445 else: ax2_name=par2_name
4446
4447 samp=np.transpose(np.column_stack((xdat,ydat)))
4448
4449 try:
4450 kde=stats.kde.gaussian_kde(samp)
4451 den=kde(samp)
4452 except:
4453 return None
4454
4455
4456 Nx=Npixels
4457 Ny=Npixels
4458
4459
4460
4461 xsorted=np.sort(xdat,axis=None)
4462 ysorted=np.sort(ydat,axis=None)
4463 imin=int(0.01*len(xsorted))
4464 imax=int(0.99*len(xsorted))
4465 xmax=xsorted[imax]
4466 xmin=xsorted[imin]
4467 ymax=ysorted[imax]
4468 ymin=ysorted[imin]
4469 xmid=0.5*(xmin+xmax)
4470 ymid=0.5*(ymin+ymax)
4471 xmin=xmid+(xmin-xmid)/0.95
4472 xmax=xmid+(xmax-xmid)/0.95
4473 ymin=ymid+(ymin-ymid)/0.95
4474 ymax=ymid+(ymax-ymid)/0.95
4475 xax=np.linspace(xmin,xmax,Nx)
4476 yax=np.linspace(ymin,ymax,Ny)
4477
4478
4479
4480
4481
4482
4483
4484 x,y = np.meshgrid(xax,yax)
4485 grid_coords = np.row_stack( (x.flatten(),y.flatten()) )
4486 z = kde(grid_coords)
4487 z = z.reshape(Nx,Ny)
4488 densort=np.sort(den)[::-1]
4489 values=[]
4490 Npts=xdat.shape[0]
4491 zvalues=[]
4492 for level in levels:
4493 ilevel = int(Npts*level + 0.5)
4494 if ilevel >= Npts:
4495 ilevel = Npts-1
4496 zvalues.append(densort[ilevel])
4497 CS=plt.contour(x, y, z, zvalues,colors=[colors_by_name[name]],linestyles=line_styles )
4498 CSlst.append(CS)
4499
4500 if par1_injvalue is not None and par2_injvalue is not None:
4501 plt.plot([par1_injvalue],[par2_injvalue],'b*',scalex=False,scaley=False,markersize=12)
4502
4503 if par_trigvalues1 is not None and par_trigvalues2 is not None:
4504 par1IFOs = set([IFO for IFO in par_trigvalues1.keys()])
4505 par2IFOs = set([IFO for IFO in par_trigvalues2.keys()])
4506 IFOs = par1IFOs.intersection(par2IFOs)
4507 for IFO in IFOs:
4508 if IFO=='H1': color = 'r'
4509 elif IFO=='L1': color = 'g'
4510 elif IFO=='V1': color = 'm'
4511 else: color = 'c'
4512 plt.plot([par_trigvalues1[IFO]],[par_trigvalues2[IFO]],color=color,marker='o',scalex=False,scaley=False)
4513
4514 plt.xlabel(plot_label(par1_name))
4515 plt.ylabel(plot_label(par2_name))
4516 plt.grid()
4517
4518 if len(name_list)!=len(CSlst):
4519 raise RuntimeError("Error number of contour objects does not equal number of names! Use only *one* contour from each set to associate a name.")
4520
4521 full_name_list=[]
4522 dummy_lines=[]
4523 for plot_name in name_list:
4524 full_name_list.append(plot_name)
4525 if len(confidence_levels)>1:
4526 for ls_,cl in zip(line_styles[0:len(confidence_levels)],confidence_levels):
4527 dummy_lines.append(mpl_lines.Line2D(np.array([0.,1.]),np.array([0.,1.]),ls=ls_,color='k'))
4528 full_name_list.append('%s%%'%str(int(cl*100)))
4529
4530 fig_actor_lst = [cs.collections[0] for cs in CSlst]
4531 fig_actor_lst.extend(dummy_lines)
4532 if legend is not None:
4533 try:
4534 twodcontour_legend=plt.figlegend(tuple(fig_actor_lst), tuple(full_name_list), loc='right',framealpha=0.1)
4535 except:
4536 twodcontour_legend=plt.figlegend(tuple(fig_actor_lst), tuple(full_name_list), loc='right')
4537 for text in twodcontour_legend.get_texts():
4538 text.set_fontsize('small')
4539
4540 majorFormatterX=ScalarFormatter(useMathText=True)
4541 majorFormatterX.format_data=lambda data:'%.4g'%(data)
4542 majorFormatterY=ScalarFormatter(useMathText=True)
4543 majorFormatterY.format_data=lambda data:'%.4g'%(data)
4544 majorFormatterX.set_scientific(True)
4545 majorFormatterY.set_scientific(True)
4546 axes.xaxis.set_major_formatter(majorFormatterX)
4547 axes.yaxis.set_major_formatter(majorFormatterY)
4548 Nchars=max(map(lambda d:len(majorFormatterX.format_data(d)),axes.get_xticks()))
4549 if Nchars>8:
4550 Nticks=3
4551 elif Nchars>5:
4552 Nticks=4
4553 elif Nchars>4:
4554 Nticks=5
4555 else:
4556 Nticks=6
4557 locatorX=matplotlib.ticker.MaxNLocator(nbins=Nticks-1)
4558 if par1_name=='rightascension' or par1_name=='ra':
4559 (ramin,ramax)=plt.xlim()
4560 locatorX=RALocator(min=ramin,max=ramax)
4561 majorFormatterX=RAFormatter()
4562 if par1_name=='declination' or par1_name=='dec':
4563 (decmin,decmax)=plt.xlim()
4564 locatorX=DecLocator(min=decmin,max=decmax)
4565 majorFormatterX=DecFormatter()
4566 axes.xaxis.set_major_formatter(majorFormatterX)
4567 if par2_name=='rightascension' or par2_name=='ra':
4568 (ramin,ramax)=plt.ylim()
4569 locatorY=RALocator(ramin,ramax)
4570 axes.yaxis.set_major_locator(locatorY)
4571 majorFormatterY=RAFormatter()
4572 if par2_name=='declination' or par2_name=='dec':
4573 (decmin,decmax)=plt.ylim()
4574 locatorY=DecLocator(min=decmin,max=decmax)
4575 majorFormatterY=DecFormatter()
4576 axes.yaxis.set_major_locator(locatorY)
4577
4578 axes.yaxis.set_major_formatter(majorFormatterY)
4579
4580 axes.xaxis.set_major_locator(locatorX)
4581
4582 fix_axis_names(plt,par1_name,par2_name)
4583
4584 if(par1_name.lower()=='ra' or par1_name.lower()=='rightascension'):
4585 xmin,xmax=plt.xlim()
4586 if(xmin<0.0): xmin=0.0
4587 if(xmax>2.0*pi_constant): xmax=2.0*pi_constant
4588 plt.xlim(xmax,xmin)
4589
4590 return fig
4591
4594 """
4595 Plots a 2D kernel density estimate of the 2-parameter marginal posterior.
4596
4597 @param posterior: an instance of the Posterior class.
4598
4599 @param plot2DkdeParams: a dict {param1Name:Nparam1Bins,param2Name:Nparam2Bins}
4600 """
4601
4602 from scipy import seterr as sp_seterr
4603
4604 par1_name,par2_name=plot2DkdeParams.keys()
4605 Nx=plot2DkdeParams[par1_name]
4606 Ny=plot2DkdeParams[par2_name]
4607
4608 xdat=posterior[par1_name].samples
4609 ydat=posterior[par2_name].samples
4610
4611 par_injvalue1=posterior[par1_name].injval
4612 par_injvalue2=posterior[par2_name].injval
4613
4614 par_trigvalues1=posterior[par1_name].trigvals
4615 par_trigvalues2=posterior[par2_name].trigvals
4616
4617 np.seterr(under='ignore')
4618 sp_seterr(under='ignore')
4619
4620 myfig=plt.figure(1,figsize=(6,4),dpi=200)
4621 myfig.add_axes(plt.Axes(myfig,[0.20,0.20,0.75,0.7]))
4622 plt.clf()
4623
4624 xax=np.linspace(min(xdat),max(xdat),Nx)
4625 yax=np.linspace(min(ydat),max(ydat),Ny)
4626 x,y=np.meshgrid(xax,yax)
4627
4628 samp=np.transpose(np.column_stack((xdat,ydat)))
4629
4630 kde=stats.kde.gaussian_kde(samp)
4631
4632 grid_coords = np.append(x.reshape(-1,1),y.reshape(-1,1),axis=1)
4633
4634 z = kde(grid_coords.T)
4635 z = z.reshape(Nx,Ny)
4636
4637 values=[]
4638 for level in levels:
4639 ilevel = int(Npts*level + 0.5)
4640 if ilevel >= Npts:
4641 ilevel = Npts-1
4642 zvalues.append(densort[ilevel])
4643
4644 pp.contour(XS, YS, ZS, zvalues)
4645
4646 asp=xax.ptp()/yax.ptp()
4647
4648 plt.imshow(z,extent=(xax[0],xax[-1],yax[0],yax[-1]),aspect=asp,origin='lower')
4649 plt.colorbar()
4650
4651 if par_injvalue1 is not None and par_injvalue2 is not None:
4652 plt.plot([par_injvalue1],[par_injvalue2],'bo',scalex=False,scaley=False)
4653
4654 if par_trigvalues1 is not None and par_trigvalues2 is not None:
4655 par1IFOs = set([IFO for IFO in par_trigvalues1.keys()])
4656 par2IFOs = set([IFO for IFO in par_trigvalues2.keys()])
4657 IFOs = par1IFOs.intersection(par2IFOs)
4658 for IFO in IFOs:
4659 if IFO=='H1': color = 'r'
4660 elif IFO=='L1': color = 'g'
4661 elif IFO=='V1': color = 'm'
4662 else: color = 'c'
4663 plt.plot([par_trigvalues1[IFO]],[par_trigvalues2[IFO]],color=color,marker='o',scalex=False,scaley=False)
4664
4665 plt.xlabel(plot_label(par1_name))
4666 plt.ylabel(plot_label(par2_name))
4667 plt.grid()
4668
4669
4670 if(par1_name.lower()=='ra' or par1_name.lower()=='rightascension'):
4671 xmin,xmax=plt.xlim()
4672 plt.xlim(xmax,xmin)
4673
4674
4675
4676
4677
4678
4679
4680
4681
4682
4683
4684
4685
4686
4687
4688
4689
4690
4691
4692
4693
4694 return myfig
4695
4704
4705 -def histogram2D(posterior,greedy2Params,confidence_levels):
4706 """
4707 Returns a 2D histogram and edges for the two parameters passed in greedy2Params, plus the actual discrete confidence levels
4708 imposed by the finite number of samples.
4709 H,xedges,yedges,Hlasts = histogram2D(posterior,greedy2Params,confidence_levels)
4710 @param posterior: Posterior instance
4711 @param greedy2Params: a dict ;{param1Name:param1binSize,param2Name:param2binSize}
4712 @param confidence_levels: a list of the required confidence levels to plot on the contour map.
4713 """
4714
4715 par1_name,par2_name=greedy2Params.keys()
4716 par1_bin=greedy2Params[par1_name]
4717 par2_bin=greedy2Params[par2_name]
4718 par1_injvalue=posterior[par1_name.lower()].injval
4719 par2_injvalue=posterior[par2_name.lower()].injval
4720 a=np.squeeze(posterior[par1_name].samples)
4721 b=np.squeeze(posterior[par2_name].samples)
4722 par1pos_min=a.min()
4723 par2pos_min=b.min()
4724 par1pos_max=a.max()
4725 par2pos_max=b.max()
4726 par1pos_Nbins= int(ceil((par1pos_max - par1pos_min)/par1_bin))+1
4727 par2pos_Nbins= int(ceil((par2pos_max - par2pos_min)/par2_bin))+1
4728 H, xedges, yedges = np.histogram2d(a,b, bins=(par1pos_Nbins, par2pos_Nbins),normed=True)
4729 temp=np.copy(H)
4730 temp=temp.ravel()
4731 confidence_levels.sort()
4732 Hsum=0
4733 Hlasts=[]
4734 idxes=np.argsort(temp)
4735 j=len(idxes)-1
4736 for cl in confidence_levels:
4737 while float(Hsum/np.sum(H))<cl:
4738
4739 max_i=idxes[j]
4740 j-=1
4741 val = temp[max_i]
4742 Hlast=val
4743 Hsum+=val
4744 temp[max_i]=0
4745 Hlasts.append(Hlast)
4746 return (H,xedges,yedges,Hlasts)
4747
4748 -def plot_two_param_greedy_bins_contourf(posteriors_by_name,greedy2Params,confidence_levels,colors_by_name,figsize=(7,6),dpi=120,figposition=[0.3,0.3,0.5,0.5],legend='right',hatches_by_name=None):
4749 """
4750 @param posteriors_by_name A dictionary of posterior objects indexed by name
4751 @param greedy2Params: a dict ;{param1Name:param1binSize,param2Name:param2binSize}
4752 @param confidence_levels: a list of the required confidence levels to plot on the contour map.
4753
4754 """
4755 fig=plt.figure(1,figsize=figsize,dpi=dpi)
4756 plt.clf()
4757 fig.add_axes(figposition)
4758 CSlst=[]
4759 name_list=[]
4760 par1_name,par2_name=greedy2Params.keys()
4761 for name,posterior in posteriors_by_name.items():
4762 name_list.append(name)
4763 H,xedges,yedges,Hlasts=histogram2D(posterior,greedy2Params,confidence_levels+[0.99999999999999])
4764 extent= [xedges[0], yedges[-1], xedges[-1], xedges[0]]
4765 CS2=plt.contourf(yedges[:-1],xedges[:-1],H,Hlasts,extend='max',colors=[colors_by_name[name]] ,alpha=0.3 )
4766 CS=plt.contour(yedges[:-1],xedges[:-1],H,Hlasts,extend='max',colors=[colors_by_name[name]] )
4767 CSlst.append(CS)
4768
4769 plt.title("%s-%s confidence contours (greedy binning)"%(par1_name,par2_name))
4770 plt.xlabel(plot_label(par2_name))
4771 plt.ylabel(plot_label(par1_name))
4772 if len(name_list)!=len(CSlst):
4773 raise RuntimeError("Error number of contour objects does not equal number of names! Use only *one* contour from each set to associate a name.")
4774 full_name_list=[]
4775 dummy_lines=[]
4776 for plot_name in name_list:
4777 full_name_list.append(plot_name)
4778 if len(confidence_levels)>1:
4779 for cl in confidence_levels+[1]:
4780 dummy_lines.append(mpl_lines.Line2D(np.array([0.,1.]),np.array([0.,1.]),color='k'))
4781 full_name_list.append('%s%%'%str(int(cl*100)))
4782 fig_actor_lst = [cs.collections[0] for cs in CSlst]
4783 fig_actor_lst.extend(dummy_lines)
4784 if legend is not None:
4785 try:
4786 twodcontour_legend=plt.figlegend(tuple(fig_actor_lst), tuple(full_name_list), loc='right',framealpha=0.1)
4787 except:
4788 twodcontour_legend=plt.figlegend(tuple(fig_actor_lst), tuple(full_name_list), loc='right')
4789 for text in twodcontour_legend.get_texts():
4790 text.set_fontsize('small')
4791 fix_axis_names(plt,par1_name,par2_name)
4792 return fig
4793
4795 """
4796 Fixes names of axes
4797 """
4798 return
4799
4800
4801 if(par1_name.lower()=='ra' or par1_name.lower()=='rightascension'):
4802 ymin,ymax=plt.ylim()
4803 if(ymin<0.0): ylim=0.0
4804 if(ymax>2.0*pi_constant): ymax=2.0*pi_constant
4805 plt.ylim(ymax,ymin)
4806 if(par1_name.lower()=='ra' or par1_name.lower()=='rightascension'):
4807 locs, ticks = plt.yticks()
4808 newlocs, newticks = formatRATicks(locs)
4809 plt.yticks(newlocs,newticks)
4810 if(par1_name.lower()=='dec' or par1_name.lower()=='declination'):
4811 locs, ticks = plt.yticks()
4812 newlocs,newticks=formatDecTicks(locs)
4813 plt.yticks(newlocs,newticks)
4814
4815 if(par2_name.lower()=='ra' or par2_name.lower()=='rightascension'):
4816 xmin,xmax=plt.xlim()
4817 if(xmin<0.0): xmin=0.0
4818 if(xmax>2.0*pi_constant): xmax=2.0*pi_constant
4819 plt.xlim(xmax,xmin)
4820 if(par2_name.lower()=='ra' or par2_name.lower()=='rightascension'):
4821 locs, ticks = plt.xticks()
4822 newlocs, newticks = formatRATicks(locs)
4823 plt.xticks(newlocs,newticks,rotation=45)
4824 if(par2_name.lower()=='dec' or par2_name.lower()=='declination'):
4825 locs, ticks = plt.xticks()
4826 newlocs, newticks = formatDecTicks(locs)
4827 plt.xticks(newlocs,newticks,rotation=45)
4828 return plt
4829
4830 -def plot_two_param_greedy_bins_contour(posteriors_by_name,greedy2Params,confidence_levels,colors_by_name,line_styles=__default_line_styles,figsize=(4,3),dpi=250,figposition=[0.2,0.2,0.48,0.75],legend='right'):
4831 """
4832 Plots the confidence level contours as determined by the 2-parameter
4833 greedy binning algorithm.
4834
4835 @param posteriors_by_name: A dict containing Posterior instances referenced by some id.
4836
4837 @param greedy2Params: a dict ;{param1Name:param1binSize,param2Name:param2binSize}
4838
4839 @param confidence_levels: a list of the required confidence levels to plot on the contour map.
4840
4841 @param colors_by_name: A dict of colors cross-referenced to the above Posterior ids.
4842
4843 @param legend: Argument for legend placement or None for no legend ('right', 'upper left', 'center' etc)
4844
4845 """
4846
4847 fig=plt.figure(1,figsize=figsize,dpi=dpi)
4848 plt.clf()
4849
4850 axes=fig.add_axes(figposition)
4851
4852
4853 if len(line_styles)<len(confidence_levels):
4854 raise RuntimeError("Error: Need as many or more line styles to choose from as confidence levels to plot!")
4855
4856 CSlst=[]
4857 name_list=[]
4858 for name,posterior in posteriors_by_name.items():
4859
4860 name_list.append(name)
4861
4862 par1_name,par2_name=greedy2Params.keys()
4863
4864 par1_bin=greedy2Params[par1_name]
4865 par2_bin=greedy2Params[par2_name]
4866
4867
4868 par1_injvalue=posterior[par1_name.lower()].injval
4869 par2_injvalue=posterior[par2_name.lower()].injval
4870
4871
4872 par1_trigvalues=posterior[par1_name.lower()].trigvals
4873 par2_trigvalues=posterior[par2_name.lower()].trigvals
4874
4875 a=np.squeeze(posterior[par1_name].samples)
4876 b=np.squeeze(posterior[par2_name].samples)
4877
4878
4879 par1pos_min=a.min()
4880 par2pos_min=b.min()
4881
4882 par1pos_max=a.max()
4883 par2pos_max=b.max()
4884
4885 par1pos_Nbins= int(ceil((par1pos_max - par1pos_min)/par1_bin))+1
4886 par2pos_Nbins= int(ceil((par2pos_max - par2pos_min)/par2_bin))+1
4887
4888 if par1_name.find('time')!=-1:
4889 offset=floor(min(a))
4890 a=a-offset
4891 if par1_injvalue:
4892 par1_injvalue=par1_injvalue-offset
4893 ax1_name=par1_name+' + %i'%(int(offset))
4894 else: ax1_name=par1_name
4895
4896 if par2_name.find('time')!=-1:
4897 offset=floor(min(b))
4898 b=b-offset
4899 if par2_injvalue:
4900 par2_injvalue=par2_injvalue-offset
4901 ax2_name=par2_name+' + %i'%(int(offset))
4902 else: ax2_name=par2_name
4903
4904
4905 majorFormatterX=ScalarFormatter(useMathText=True)
4906 majorFormatterX.format_data=lambda data:'%.4g'%(data)
4907 majorFormatterY=ScalarFormatter(useMathText=True)
4908 majorFormatterY.format_data=lambda data:'%.4g'%(data)
4909 majorFormatterX.set_scientific(True)
4910 majorFormatterY.set_scientific(True)
4911 axes.xaxis.set_major_formatter(majorFormatterX)
4912 axes.yaxis.set_major_formatter(majorFormatterY)
4913
4914 H, xedges, yedges = np.histogram2d(a,b, bins=(par1pos_Nbins, par2pos_Nbins),normed=True)
4915
4916 extent = [xedges[0], yedges[-1], xedges[-1], xedges[0]]
4917
4918 temp=np.copy(H)
4919 temp=temp.ravel()
4920 confidence_levels.sort()
4921 Hsum=0
4922 Hlasts=[]
4923 idxes=np.argsort(temp)
4924 j=len(idxes)-1
4925 for cl in confidence_levels:
4926 while float(Hsum/np.sum(H))<cl:
4927
4928 max_i=idxes[j]
4929 j-=1
4930 val = temp[max_i]
4931 Hlast=val
4932 Hsum+=val
4933 temp[max_i]=0
4934 Hlasts.append(Hlast)
4935 CS=plt.contour(yedges[:-1],xedges[:-1],H,Hlasts,colors=[colors_by_name[name]],linestyles=line_styles)
4936 plt.grid()
4937 if(par1_injvalue is not None and par2_injvalue is not None):
4938 plt.plot([par2_injvalue],[par1_injvalue],'b*',scalex=False,scaley=False,markersize=12)
4939 if(par1_trigvalues is not None and par2_trigvalues is not None):
4940 par1IFOs = set([IFO for IFO in par1_trigvalues.keys()])
4941 par2IFOs = set([IFO for IFO in par2_trigvalues.keys()])
4942 IFOs = par1IFOs.intersection(par2IFOs)
4943 if IFO=='H1': color = 'r'
4944 elif IFO=='L1': color = 'g'
4945 elif IFO=='V1': color = 'm'
4946 else: color = 'c'
4947 plt.plot([par2_trigvalues[IFO]],[par1_trigvalues[IFO]],color=color,marker='*',scalex=False,scaley=False)
4948 CSlst.append(CS)
4949
4950 Nchars=max(map(lambda d:len(majorFormatterX.format_data(d)),axes.get_xticks()))
4951 if Nchars>8:
4952 Nticks=3
4953 elif Nchars>5:
4954 Nticks=4
4955 elif Nchars>4:
4956 Nticks=5
4957 else:
4958 Nticks=6
4959 locatorX=matplotlib.ticker.MaxNLocator(nbins=Nticks-1)
4960 if par2_name=='rightascension' or par2_name=='ra':
4961 (ramin,ramax)=plt.xlim()
4962 locatorX=RALocator(min=ramin,max=ramax)
4963 majorFormatterX=RAFormatter()
4964 if par2_name=='declination' or par2_name=='dec':
4965 (decmin,decmax)=plt.xlim()
4966 locatorX=DecLocator(min=decmin,max=decmax)
4967 majorFormatterX=DecFormatter()
4968 axes.xaxis.set_major_formatter(majorFormatterX)
4969 if par1_name=='rightascension' or par1_name=='ra':
4970 (ramin,ramax)=plt.ylim()
4971 locatorY=RALocator(ramin,ramax)
4972 axes.yaxis.set_major_locator(locatorY)
4973 majorFormatterY=RAFormatter()
4974 if par1_name=='declination' or par1_name=='dec':
4975 (decmin,decmax)=plt.ylim()
4976 locatorY=DecLocator(min=decmin,max=decmax)
4977 majorFormatterY=DecFormatter()
4978 axes.yaxis.set_major_locator(locatorY)
4979
4980 axes.yaxis.set_major_formatter(majorFormatterY)
4981
4982 axes.xaxis.set_major_locator(locatorX)
4983
4984
4985 plt.xlabel(plot_label(ax2_name))
4986 plt.ylabel(plot_label(ax1_name))
4987
4988 if len(name_list)!=len(CSlst):
4989 raise RuntimeError("Error number of contour objects does not equal number of names! Use only *one* contour from each set to associate a name.")
4990 full_name_list=[]
4991 dummy_lines=[]
4992
4993 for plot_name in name_list:
4994 full_name_list.append(plot_name)
4995 if len(confidence_levels)>1:
4996 for ls_,cl in zip(line_styles[0:len(confidence_levels)],confidence_levels):
4997 dummy_lines.append(mpl_lines.Line2D(np.array([0.,1.]),np.array([0.,1.]),ls=ls_,color='k'))
4998 full_name_list.append('%s%%'%str(int(cl*100)))
4999
5000 fig_actor_lst = [cs.collections[0] for cs in CSlst]
5001
5002 fig_actor_lst.extend(dummy_lines)
5003
5004 if legend is not None:
5005 try:
5006 twodcontour_legend=plt.figlegend(tuple(fig_actor_lst), tuple(full_name_list), loc='right',framealpha=0.1)
5007 except:
5008 twodcontour_legend=plt.figlegend(tuple(fig_actor_lst), tuple(full_name_list), loc='right')
5009 for text in twodcontour_legend.get_texts():
5010 text.set_fontsize('small')
5011
5012
5013
5014
5015
5016
5017
5018
5019
5020
5021
5022
5023
5024
5025
5026
5027
5028 if(par2_name.lower()=='ra' or par2_name.lower()=='rightascension'):
5029 xmin,xmax=plt.xlim()
5030 if(xmin<0.0): xmin=0.0
5031 if(xmax>2.0*pi_constant): xmax=2.0*pi_constant
5032 plt.xlim(xmax,xmin)
5033
5034
5035
5036
5037
5038
5039
5040
5041
5042 return fig
5043
5046 """
5047 Histograms of the ranked pixels produced by the 2-parameter greedy
5048 binning algorithm colured by their confidence level.
5049
5050 @param toppoints: Nx2 array of 2-parameter posterior samples.
5051
5052 @param posterior: an instance of the Posterior class.
5053
5054 @param greedy2Params: a dict ;{param1Name:param1binSize,param2Name:param2binSize}
5055 """
5056
5057 from scipy import seterr as sp_seterr
5058
5059 np.seterr(under='ignore')
5060 sp_seterr(under='ignore')
5061
5062
5063 par1_name,par2_name=greedy2Params.keys()
5064
5065 par1_bin=greedy2Params[par1_name]
5066 par2_bin=greedy2Params[par2_name]
5067
5068 a=np.squeeze(posterior[par1_name].samples)
5069 b=np.squeeze(posterior[par2_name].samples)
5070
5071
5072 par1_injvalue=posterior[par1_name.lower()].injval
5073 par2_injvalue=posterior[par2_name.lower()].injval
5074
5075
5076 par1pos_min=a.min()
5077 par2pos_min=b.min()
5078
5079 par1pos_max=a.max()
5080 par2pos_max=b.max()
5081
5082 par1pos_Nbins= int(ceil((par1pos_max - par1pos_min)/par1_bin))+1
5083 par2pos_Nbins= int(ceil((par2pos_max - par2pos_min)/par2_bin))+1
5084
5085
5086 if par1_name.find('time')!=-1:
5087 offset=floor(min(a))
5088 a=a-offset
5089 if par1_injvalue:
5090 par1_injvalue=par1_injvalue-offset
5091 ax1_name=par1_name+' + %i'%(int(offset))
5092 else: ax1_name=par1_name
5093
5094 if par2_name.find('time')!=-1:
5095 offset=floor(min(b))
5096 b=b-offset
5097 if par2_injvalue:
5098 par2_injvalue=par2_injvalue-offset
5099 ax2_name=par2_name+' + %i'%(int(offset))
5100 else: ax2_name=par2_name
5101
5102
5103
5104 par1_trigvalues=posterior[par1_name.lower()].trigvals
5105 par2_trigvalues=posterior[par2_name.lower()].trigvals
5106
5107 myfig=plt.figure()
5108 axes=plt.Axes(myfig,[0.3,0.3,0.95-0.3,0.90-0.3])
5109 myfig.add_axes(axes)
5110
5111
5112 plt.xlabel(plot_label(ax2_name))
5113 plt.ylabel(plot_label(ax1_name))
5114
5115
5116 bins=(50,50)
5117
5118 majorFormatterX=ScalarFormatter(useMathText=True)
5119 majorFormatterX.format_data=lambda data:'%.4g'%(data)
5120 majorFormatterY=ScalarFormatter(useMathText=True)
5121 majorFormatterY.format_data=lambda data:'%.4g'%(data)
5122 majorFormatterX.set_scientific(True)
5123 majorFormatterY.set_scientific(True)
5124 axes.xaxis.set_major_formatter(majorFormatterX)
5125 axes.yaxis.set_major_formatter(majorFormatterY)
5126 H, xedges, yedges = np.histogram2d(a,b, bins,normed=False)
5127
5128
5129
5130 temp=np.copy(H)
5131 temp=temp.flatten()
5132
5133 Hsum=0
5134 Hsum_actual=np.sum(H)
5135
5136 idxes=np.argsort(temp)
5137 j=len(idxes)-1
5138 while Hsum<Hsum_actual:
5139
5140 max_i=idxes[j]
5141 j-=1
5142 val = temp[max_i]
5143 Hsum+=int(val)
5144 temp[max_i]=0
5145
5146
5147 H.flat[max_i]=1-float(Hsum)/float(Hsum_actual)
5148
5149 extent = [yedges[0], yedges[-1], xedges[0], xedges[-1]]
5150 plt.imshow(np.flipud(H), axes=axes, aspect='auto', extent=extent, interpolation='nearest',cmap='gray_r')
5151 plt.gca().autoscale_view()
5152 plt.colorbar()
5153
5154
5155
5156 Nchars=max(map(lambda d:len(majorFormatterX.format_data(d)),axes.get_xticks()))
5157 if Nchars>8:
5158 Nticks=3
5159 elif Nchars>5:
5160 Nticks=4
5161 elif Nchars>4:
5162 Nticks=5
5163 else:
5164 Nticks=6
5165 locatorX=matplotlib.ticker.MaxNLocator(nbins=Nticks-1)
5166 (xmin,xmax)=plt.xlim()
5167 (ymin,ymax)=plt.ylim()
5168 if par2_name=='rightascension' or par2_name=='ra':
5169 locatorX=RALocator(min=xmin,max=xmax)
5170 majorFormatterX=RAFormatter()
5171 if par2_name=='declination' or par2_name=='dec':
5172 locatorX=DecLocator(min=xmin,max=xmax)
5173 majorFormatterX=DecFormatter()
5174 if par1_name=='rightascension' or par1_name=='ra':
5175 locatorY=RALocator(min=ymin,max=ymax)
5176 axes.yaxis.set_major_locator(locatorY)
5177 majorFormatterY=RAFormatter()
5178 if par1_name=='declination' or par1_name=='dec':
5179 locatorY=DecLocator(min=ymin,max=ymax)
5180 axes.yaxis.set_major_locator(locatorY)
5181 majorFormatterY=DecFormatter()
5182
5183 axes.xaxis.set_major_formatter(majorFormatterX)
5184 axes.yaxis.set_major_formatter(majorFormatterY)
5185
5186 axes.xaxis.set_major_locator(locatorX)
5187
5188 if par1_injvalue is not None and par2_injvalue is not None:
5189 plt.plot([par2_injvalue],[par1_injvalue],'bo',scalex=False,scaley=False)
5190
5191 if par1_trigvalues is not None and par2_trigvalues is not None:
5192 par1IFOs = set([IFO for IFO in par1_trigvalues.keys()])
5193 par2IFOs = set([IFO for IFO in par2_trigvalues.keys()])
5194 IFOs = par1IFOs.intersection(par2IFOs)
5195 if IFO=='H1': color = 'r'
5196 elif IFO=='L1': color = 'g'
5197 elif IFO=='V1': color = 'm'
5198 else: color = 'c'
5199 plt.plot([par2_trigvalues[IFO]],[par1_trigvalues[IFO]],color=color,marker='o',scalex=False,scaley=False)
5200
5201
5202
5203
5204
5205
5206
5207
5208
5209
5210
5211
5212
5213
5214 if(par2_name.lower()=='ra' or par2_name.lower()=='rightascension'):
5215 xmin,xmax=plt.xlim()
5216 if(xmin)<0.0: xmin=0.0
5217 if(xmax>2.0*pi_constant): xmax=2.0*pi_constant
5218 plt.xlim(xmax,xmin)
5219
5220
5221
5222
5223
5224
5225
5226
5227
5228 return myfig
5229
5231 """
5232 Determine the 1-parameter Bayesian Confidence Interval using a greedy
5233 binning algorithm.
5234
5235 @param posterior: an instance of the posterior class.
5236
5237 @param greedy1Param: a dict; {paramName:paramBinSize}.
5238
5239 @param confidence_levels: A list of floats of the required confidence intervals [(0-1)].
5240 """
5241
5242 paramName=greedy1Param.keys()[0]
5243 par_bin=greedy1Param.values()[0]
5244 par_samps=posterior[paramName.lower()].samples
5245
5246 parpos_min=min(par_samps)[0]
5247 parpos_max=max(par_samps)[0]
5248
5249 par_point=parpos_min
5250
5251 parpos_Nbins= int(ceil((parpos_max - parpos_min)/par_bin))+1
5252
5253 greedyPoints=np.zeros((parpos_Nbins,2))
5254
5255 greedyHist=np.zeros(parpos_Nbins,dtype='i8')
5256
5257
5258 for i in range(parpos_Nbins):
5259 greedyPoints[i,0]=par_point
5260 greedyPoints[i,1]=par_point
5261 par_point+=par_bin
5262
5263 for par_samp in par_samps:
5264 par_samp=par_samp[0]
5265 par_binNumber=int(floor((par_samp-parpos_min)/par_bin))
5266 try:
5267 greedyHist[par_binNumber]+=1
5268 except IndexError:
5269 print "IndexError: bin number: %i total bins: %i parsamp: %f "\
5270 %(par_binNumber,parpos_Nbins,par_samp)
5271
5272
5273 injbin=None
5274 par_injvalue=posterior[paramName].injval
5275 if par_injvalue:
5276 par_binNumber=floor((par_injvalue-parpos_min)/par_bin)
5277 injbin=par_binNumber
5278
5279 toppoints,injectionconfidence,reses,injection_area=_greedy_bin(greedyHist,greedyPoints,injbin,float(par_bin),int(len(par_samps)),confidence_levels)
5280 cl_intervals=[]
5281 confidence_levels.sort()
5282 for cl in confidence_levels:
5283 ind=np.nonzero(toppoints[:,-1]<cl)
5284
5285 if len(ind[0]) > 1:
5286 cl_intervals.append((np.min(toppoints[ind,0]),np.max(toppoints[ind,0])))
5287
5288 else:
5289
5290 cl_intervals.append((toppoints[ind[0],0],toppoints[ind[0],0]))
5291
5292 return toppoints,injectionconfidence,reses,injection_area,cl_intervals
5293
5296 """
5297 Calculates the smallest contigious 1-parameter confidence interval for a
5298 set of given confidence levels.
5299
5300 @param posterior: an instance of the Posterior class.
5301
5302 @param contInt1Params: a dict {paramName:paramBinSize}.
5303
5304 @param confidence_levels: Required confidence intervals.
5305
5306 """
5307 oneDContCL={}
5308 oneDContInj={}
5309
5310 paramName=contInt1Params.keys()[0]
5311 par_bin=contInt1Params.values()[0]
5312
5313 par_injvalue=posterior[paramName].injval
5314
5315 par_samps=posterior[paramName].samples
5316
5317 parpos_min=min(par_samps)
5318 parpos_max=max(par_samps)
5319
5320 par_point=parpos_min
5321 parpos_Nbins= int(ceil((parpos_max - parpos_min)/par_bin))+1
5322
5323 greedyHist=np.zeros(parpos_Nbins,dtype='i8')
5324
5325 for par_samp in par_samps:
5326 par_binNumber=int(floor((par_samp-parpos_min)/par_bin))
5327 try:
5328 greedyHist[par_binNumber]+=1
5329 except IndexError:
5330 print "IndexError: bin number: %i total bins: %i parsamp: %f bin: %f - %f"\
5331 %(
5332 par_binNumber,
5333 parpos_Nbins,
5334 par_samp,
5335 greedyPoints[par_binNumber-1,0],
5336 greedyPoints[par_binNumber-1,0]+par_bin
5337 )
5338
5339 injbin=None
5340
5341 if par_injvalue:
5342 par_binNumber=floor((par_injvalue-parpos_min)/par_bin)
5343 injbin=par_binNumber
5344
5345 j=0
5346
5347 len_par_samps=len(par_samps)
5348
5349 injinterval=None
5350
5351
5352 while j < len(confidence_levels):
5353 confidence_level=confidence_levels[j]
5354
5355 max_left=0
5356 max_right=0
5357
5358 for i in range(len(greedyHist)):
5359
5360 max_frac=None
5361 left=0
5362 right=i
5363
5364
5365 while right<len(greedyHist):
5366 Npoints=sum(greedyHist[left:right])
5367 frac=float(Npoints)/float(len_par_samps)
5368
5369
5370 if max_frac is None:
5371 max_frac=frac
5372 max_left=left
5373 max_right=right
5374 else:
5375 if frac>max_frac:
5376 max_frac=frac
5377 max_left=left
5378 max_right=right
5379
5380 left+=1
5381 right+=1
5382
5383 if injbin is not None and injinterval is None:
5384 if injbin in range(max_left,max_right):
5385 injinterval=(max_right-max_left)*par_bin
5386 oneDContInj['interval']=injinterval
5387 oneDContInj['confidence']=1-frac
5388 if max_frac > confidence_level:
5389 break
5390
5391 max_frac=None
5392
5393 if max_frac is None:
5394 print "Cant determine intervals at %f confidence!"%confidence_level
5395 else:
5396
5397 oneDContCL['left']=max_left*par_bin
5398 oneDContCL['right']=max_right*par_bin
5399 oneDContCL['width']=(max_right-max_left)*par_bin
5400 k=j
5401 while k+1<len(confidence_levels) :
5402 if confidence_levels[k+1]<max_frac:
5403 j+=1
5404 k+=1
5405 j+=1
5406
5407 return oneDContCL,oneDContInj
5408
5409 -def burnin(data,spin_flag,deltaLogP,outputfile):
5410
5411 pos,bayesfactor=_burnin(data,spin_flag,deltaLogP,outputfile)
5412
5413 return pos,bayesfactor
5414
5419
5422 """Returns an estimate of the autocorrelation function of a given
5423 series. Returns only the positive-lag portion of the ACF,
5424 normalized so that the zero-th element is 1."""
5425 x=series-np.mean(series)
5426 y=np.conj(x[::-1])
5427
5428 acf=np.fft.ifftshift(signal.fftconvolve(y,x,mode='full'))
5429
5430 N=series.shape[0]
5431
5432 acf = acf[0:N]
5433
5434 return acf/acf[0]
5435
5438 """Attempts to find a self-consistent estimate of the
5439 autocorrelation length of a given series.
5440
5441 If C(tau) is the autocorrelation function (normalized so C(0) = 1,
5442 for example from the autocorrelation procedure in this module),
5443 then the autocorrelation length is the smallest s such that
5444
5445 1 + 2*C(1) + 2*C(2) + ... + 2*C(M*s) < s
5446
5447 In words: the autocorrelation length is the shortest length so
5448 that the sum of the autocorrelation function is smaller than that
5449 length over a window of M times that length.
5450
5451 The maximum window length is restricted to be len(series)/K as a
5452 safety precaution against relying on data near the extreme of the
5453 lags in the ACF, where there is a lot of noise. Note that this
5454 implies that the series must be at least M*K*s samples long in
5455 order to get a reliable estimate of the ACL.
5456
5457 If no such s can be found, raises ACLError; in this case it is
5458 likely that the series is too short relative to its true
5459 autocorrelation length to obtain a consistent ACL estimate."""
5460
5461 if acf is None:
5462 acf=autocorrelation(series)
5463 acf[1:] *= 2.0
5464
5465 imax=int(acf.shape[0]/K)
5466
5467
5468 cacf=np.cumsum(acf)
5469 s=np.arange(1, cacf.shape[0]+1)/float(M)
5470
5471
5472
5473 estimates=np.flatnonzero(cacf[:imax] < s[:imax])
5474
5475 if estimates.shape[0] > 0:
5476
5477
5478 return s[estimates[0]]
5479 else:
5480
5481 raise ACLError('autocorrelation length too short for consistent estimate')
5482
5485 """
5486 Compute the effective sample size, calculating the ACL using only
5487 the second half of the samples to avoid ACL overestimation due to
5488 chains equilibrating after adaptation.
5489 """
5490 N = len(samples)
5491 acf = autocorrelation(samples[N/2:])
5492 try:
5493 acl = autocorrelation_length_estimate(samples[N/2:], acf=acf)
5494 except ACLError:
5495 acl = N
5496 Neffective = floor(N/acl)
5497 acl *= Nskip
5498 return (Neffective, acl, acf)
5499
5502 triggers=None
5503
5504 from glue.ligolw import ligolw
5505 coincXML = utils.load_filename(xml_file, contenthandler = lsctables.use_in(ligolw.LIGOLWContentHandler))
5506 coinc = lsctables.CoincTable.get_table(coincXML)
5507 coincMap = lsctables.CoincMapTable.get_table(coincXML)
5508 snglInsps = lsctables.SnglInspiralTable.get_table(coincXML)
5509
5510 if (trignum>len(coinc)):
5511 raise RuntimeError("Error: You asked for trigger %d, but %s contains only %d triggers" %(trignum,coincfile,len(tiggers)))
5512 else:
5513 coincEventID = coinc.getColumnByName('coinc_event_id')[trignum]
5514 eventIDs = [row.event_id for row in coincMap if row.coinc_event_id == coincEventID]
5515 triggers = [row for row in snglInsps if row.event_id in eventIDs]
5516 return triggers
5517
5523 """
5524 Given a list of files, threshold value, and a desired
5525 number of outputs posterior samples, return the skip number to
5526 achieve the desired number of posterior samples.
5527 """
5528 if nDownsample is None:
5529 print "Max ACL(s):"
5530 splineParams=["spcal_npts", "spcal_active","constantcal_active"]
5531 for i in np.arange(25):
5532 for k in lal.cached_detector_by_prefix:
5533 splineParams.append(k+'_spcal_freq_'+str(i))
5534 splineParams.append(k+'_spcal_logfreq_'+str(i))
5535
5536 nonParams = ["logpost", "post", "cycle", "timestamp", "snrh1", "snrl1", "snrv1",
5537 "margtime","margtimephi","margtime","time_max","time_min",
5538 "time_mean", "time_maxl","sky_frame","psdscaleflag","logdeltaf","flow","f_ref",
5539 "lal_amporder","lal_pnorder","lal_approximant","tideo","spino","signalmodelflag",
5540 "temperature","nifo","nlocaltemps","ntemps","randomseed","samplerate","segmentlength","segmentstart",
5541 "t0", "phase_maxl", "azimuth", "cosalpha", "lal_amporder"] + logParams + snrParams + splineParams
5542 fixedParams = [p for p in samples.colnames if all(x==samples[p][0] for x in samples[p])]
5543 print "Fixed parameters: "+str(fixedParams)
5544 nonParams.extend(fixedParams)
5545 params = [p for p in samples.colnames if p.lower() not in nonParams]
5546 stride=np.diff(samples['cycle'])[0]
5547 results = np.array([np.array(effectiveSampleSize(samples[param])[:2]) for param in params])
5548 nEffs = results[:,0]
5549 nEffective = min(nEffs)
5550 ACLs = results[:,1]
5551 maxACLind = np.argmax(ACLs)
5552 maxACL = ACLs[maxACLind]
5553
5554 print "%i (%s)." %(stride*maxACL,params[maxACLind])
5555
5556 nskip = 1
5557 if nDownsample is not None:
5558 if len(samples) > nDownsample:
5559 nskip *= floor(len(samples)/nDownsample)
5560 nskip = int(nskip)
5561 else:
5562 nEff = nEffective
5563 if nEff > 1:
5564 if len(samples) > nEff:
5565 nskip = int(ceil(len(samples)/nEff))
5566 else:
5567 nskip = np.nan
5568 return nskip
5569
5571 """
5572 A parser for the output of Bayesian parameter estimation codes.
5573
5574 TODO: Will be abstract class when LDG moves over to Python >2.6,
5575 inherited by each method .
5576 """
5578 if inputtype is 'mcmc_burnin':
5579 self._parser=self._mcmc_burnin_to_pos
5580 elif inputtype is 'ns':
5581 self._parser=self._ns_to_pos
5582 elif inputtype is 'common':
5583 self._parser=self._common_to_pos
5584 elif inputtype is 'fm':
5585 self._parser=self._followupmcmc_to_pos
5586 elif inputtype is "inf_mcmc":
5587 self._parser=self._infmcmc_to_pos
5588 elif inputtype is "xml":
5589 self._parser=self._xml_to_pos
5590 elif inputtype == 'hdf5':
5591 self._parser = self._hdf5_to_pos
5592 elif inputtype == 'hdf5s':
5593 self._parser = self._hdf5s_to_pos
5594 else:
5595 raise ValueError('Invalid value for "inputtype": %r' % inputtype)
5596
5597 - def parse(self,files,**kwargs):
5598 """
5599 Parse files.
5600 """
5601 return self._parser(files,**kwargs)
5602
5603 - def _infmcmc_to_pos(self,files,outdir=None,deltaLogP=None,fixedBurnins=None,nDownsample=None,oldMassConvention=False,**kwargs):
5604 """
5605 Parser for lalinference_mcmcmpi output.
5606 """
5607 if not (fixedBurnins is None):
5608 if not (deltaLogP is None):
5609 print "Warning: using deltaLogP criteria in addition to fixed burnin"
5610 if len(fixedBurnins) == 1 and len(files) > 1:
5611 print "Only one fixedBurnin criteria given for more than one output. Applying this to all outputs."
5612 fixedBurnins = np.ones(len(files),'int')*fixedBurnins[0]
5613 elif len(fixedBurnins) != len(files):
5614 raise RuntimeError("Inconsistent number of fixed burnin criteria and output files specified.")
5615 print "Fixed burning criteria: ",fixedBurnins
5616 else:
5617 fixedBurnins = np.zeros(len(files))
5618 logPThreshold=-np.inf
5619 if not (deltaLogP is None):
5620 logPThreshold= - deltaLogP
5621 print "Eliminating any samples before log(Post) = ", logPThreshold
5622 nskips=self._find_ndownsample(files, logPThreshold, fixedBurnins, nDownsample)
5623 if nDownsample is None:
5624 print "Downsampling to take only uncorrelated posterior samples from each file."
5625 if len(nskips) == 1 and np.isnan(nskips[0]):
5626 print "WARNING: All samples in chain are correlated. Downsampling to 10000 samples for inspection!!!"
5627 nskips=self._find_ndownsample(files, logPThreshold, fixedBurnins, 10000)
5628 else:
5629 for i in range(len(nskips)):
5630 if np.isnan(nskips[i]):
5631 print "%s eliminated since all samples are correlated."
5632 else:
5633 print "Downsampling by a factor of ", nskips[0], " to achieve approximately ", nDownsample, " posterior samples"
5634 if outdir is None:
5635 outdir=''
5636 runfileName=os.path.join(outdir,"lalinfmcmc_headers.dat")
5637 postName="posterior_samples.dat"
5638 runfile=open(runfileName, 'w')
5639 outfile=open(postName, 'w')
5640 try:
5641 self._infmcmc_output_posterior_samples(files, runfile, outfile, logPThreshold, fixedBurnins, nskips, oldMassConvention)
5642 finally:
5643 runfile.close()
5644 outfile.close()
5645 return self._common_to_pos(open(postName,'r'))
5646
5647
5648 - def _infmcmc_output_posterior_samples(self, files, runfile, outfile, logPThreshold, fixedBurnins, nskips=None, oldMassConvention=False):
5649 """
5650 Concatenate all the samples from the given files into outfile.
5651 For each file, only those samples past the point where the
5652 log(post) > logPThreshold are concatenated after eliminating
5653 fixedBurnin.
5654 """
5655 nRead=0
5656 outputHeader=False
5657 acceptedChains=0
5658 if nskips is None:
5659 nskips = np.ones(len(files),'int')
5660 for infilename,i,nskip,fixedBurnin in zip(files,range(1,len(files)+1),nskips,fixedBurnins):
5661 infile=open(infilename,'r')
5662 try:
5663 print "Writing header of %s to %s"%(infilename,runfile.name)
5664 runInfo,header=self._clear_infmcmc_header(infile)
5665 runfile.write('Chain '+str(i)+':\n')
5666 runfile.writelines(runInfo)
5667 print "Processing file %s to %s"%(infilename,outfile.name)
5668 write_fref = False
5669 if 'f_ref' not in header:
5670 write_fref = True
5671 f_ref=self._find_infmcmc_f_ref(runInfo)
5672 if oldMassConvention:
5673
5674
5675
5676 header=[self._swaplabel12(label) for label in header]
5677 if not outputHeader:
5678 for label in header:
5679 outfile.write(label)
5680 outfile.write(" ")
5681 if write_fref:
5682 outfile.write("f_ref")
5683 outfile.write(" ")
5684 outfile.write("chain")
5685 outfile.write("\n")
5686 outputHeader=header
5687 iterindex=header.index("cycle")
5688 logpindex=header.index("logpost")
5689 output=False
5690 for line in infile:
5691 line=line.lstrip()
5692 lineParams=line.split()
5693 iter=int(lineParams[iterindex])
5694 logP=float(lineParams[logpindex])
5695 if (iter > fixedBurnin) and (logP >= logPThreshold):
5696 output=True
5697 if output:
5698 if nRead % nskip == 0:
5699 for label in outputHeader:
5700
5701
5702
5703
5704
5705 outfile.write(lineParams[header.index(label)])
5706 outfile.write("\t")
5707 if write_fref:
5708 outfile.write(f_ref)
5709 outfile.write("\t")
5710 outfile.write(str(i))
5711 outfile.write("\n")
5712 nRead=nRead+1
5713 if output: acceptedChains += 1
5714 finally:
5715 infile.close()
5716 print "%i of %i chains accepted."%(acceptedChains,len(files))
5717
5719 if label[-1] == '1':
5720 return label[0:-1] + '2'
5721 elif label[-1] == '2':
5722 return label[0:-1] + '1'
5723 else:
5724 return label[:]
5725
5727 """
5728 Given a list of files, reads them, finding the maximum log(post)
5729 """
5730 maxLogP = -np.inf
5731 for inpname in files:
5732 infile=open(inpname, 'r')
5733 try:
5734 runInfo,header=self._clear_infmcmc_header(infile)
5735 logpindex=header.index("logpost")
5736 for line in infile:
5737 line=line.lstrip().split()
5738 logP=float(line[logpindex])
5739 if logP > maxLogP:
5740 maxLogP=logP
5741 finally:
5742 infile.close()
5743 print "Found max log(post) = ", maxLogP
5744 return maxLogP
5745
5747 """
5748 Given a list of files, threshold value, and a desired
5749 number of outputs posterior samples, return the skip number to
5750 achieve the desired number of posterior samples.
5751 """
5752 nfiles = len(files)
5753 ntots=[]
5754 nEffectives = []
5755 if nDownsample is None: print "Max ACL(s):"
5756 for inpname,fixedBurnin in zip(files,fixedBurnins):
5757 infile = open(inpname, 'r')
5758 try:
5759 runInfo,header = self._clear_infmcmc_header(infile)
5760 header = [name.lower() for name in header]
5761 logpindex = header.index("logpost")
5762 iterindex = header.index("cycle")
5763 deltaLburnedIn = False
5764 fixedBurnedIn = False
5765 adapting = True
5766 lines=[]
5767 ntot=0
5768 for line in infile:
5769 line = line.lstrip().split()
5770 iter = int(line[iterindex])
5771 logP = float(line[logpindex])
5772 if iter > fixedBurnin:
5773 fixedBurnedIn = True
5774
5775 elif fixedBurnedIn:
5776 fixedBurnedIn = False
5777 ntot = 0
5778 lines = []
5779 if logP > logPthreshold:
5780 deltaLburnedIn = True
5781 if iter > 0:
5782 adapting = False
5783 if fixedBurnedIn and deltaLburnedIn and not adapting:
5784 ntot += 1
5785 lines.append(line)
5786 ntots.append(ntot)
5787 if nDownsample is None:
5788 try:
5789 splineParams=["spcal_npts", "spcal_active","constantcal_active"]
5790 for i in np.arange(5):
5791 for k in ['h1','l1']:
5792 splineParams.append(k+'_spcal_freq'+str(i))
5793 splineParams.append(k+'_spcal_logfreq'+str(i))
5794
5795 nonParams = ["logpost", "cycle", "timestamp", "snrh1", "snrl1", "snrv1",
5796 "margtime","margtimephi","margtime","time_max","time_min",
5797 "time_mean", "time_maxl","sky_frame","psdscaleflag","logdeltaf","flow","f_ref",
5798 "lal_amporder","lal_pnorder","lal_approximant","tideo","spino","signalmodelflag",
5799 "temperature","nifo","nlocaltemps","ntemps","randomseed","samplerate","segmentlength","segmentstart",
5800 "t0", "phase_maxl", "azimuth", "cosalpha"] + logParams + snrParams + splineParams
5801 nonParamsIdxs = [header.index(name) for name in nonParams if name in header]
5802 samps = np.array(lines).astype(float)
5803 fixedIdxs = np.where(np.amin(samps,axis=0)-np.amax(samps,axis=0) == 0.0)[0]
5804 nonParamsIdxs.extend(fixedIdxs)
5805 paramIdxs = [i for i in range(len(header)) if i not in nonParamsIdxs]
5806 stride=samps[1,iterindex] - samps[0,iterindex]
5807 results = np.array([np.array(effectiveSampleSize(samps[:,i])[:2]) for i in paramIdxs])
5808 nEffs = results[:,0]
5809 nEffectives.append(min(nEffs))
5810 ACLs = results[:,1]
5811 maxACLind = np.argmax(ACLs)
5812 maxACL = ACLs[maxACLind]
5813
5814 maxACLind = paramIdxs[maxACLind]
5815 print "%i (%s) for chain %s." %(stride*maxACL,header[maxACLind],inpname)
5816 except:
5817 nEffectives.append(None)
5818 print "Error computing effective sample size of %s!"%inpname
5819
5820 finally:
5821 infile.close()
5822 nskips = np.ones(nfiles)
5823 ntot = sum(ntots)
5824 if nDownsample is not None:
5825 if ntot > nDownsample:
5826 nskips *= int(floor(ntot/nDownsample))
5827 else:
5828 for i in range(nfiles):
5829 nEff = nEffectives[i]
5830 ntot = ntots[i]
5831 if nEff > 1:
5832 if ntot > nEff:
5833 nskips[i] = int(ceil(ntot/nEff))
5834 else:
5835 nskips[i] = None
5836 return nskips
5837
5839 """
5840 Searches through header to determine reference frequency of waveforms.
5841 If no fRef given, calls _find_infmcmc_f_lower to get the lower frequency
5842 bound, which is the default reference frequency for LALInference.
5843 """
5844 fRef = None
5845 runInfoIter = iter(runInfo)
5846 for line in runInfoIter:
5847 headers=line.lstrip().lower().split()
5848 try:
5849 fRefColNum = headers.index('f_ref')
5850 info = runInfoIter.next().lstrip().lower().split()
5851 fRef = info[-1]
5852 break
5853 except ValueError:
5854 continue
5855
5856
5857
5858 if not fRef:
5859 runInfoIter = iter(runInfo)
5860 for line in runInfoIter:
5861 headers=line.lstrip().lower().split()
5862 try:
5863 if headers[0]=="command":
5864 try:
5865 fRefInd = headers.index('--fref')+1
5866 fRef = headers[fRefInd]
5867 except ValueError:
5868 pass
5869 break
5870 except IndexError:
5871 continue
5872
5873
5874 if not fRef:
5875 fRef = self._find_infmcmc_f_lower(runInfo)
5876
5877 return fRef
5878
5880 """
5881 Searches through header to determine starting frequency of waveforms.
5882 Assumes same for all IFOs.
5883 """
5884 runInfo = iter(runInfo)
5885 for line in runInfo:
5886 headers=line.lstrip().lower().split()
5887 try:
5888 flowColNum = headers.index('f_low')
5889 IFOinfo = runInfo.next().lstrip().lower().split()
5890 f_lower = IFOinfo[flowColNum]
5891 break
5892 except ValueError:
5893 continue
5894 return f_lower
5895
5897 """
5898 Reads lalinference_mcmcmpi file given, returning the run info and
5899 common output header information.
5900 """
5901 runInfo = []
5902 for line in infile:
5903 runInfo.append(line)
5904 headers=line.lstrip().lower().split()
5905 try:
5906 headers.index('cycle')
5907 break
5908 except ValueError:
5909 continue
5910 else:
5911 raise RuntimeError("couldn't find line with 'cycle' in LALInferenceMCMC input")
5912 return runInfo[:-1],headers
5913
5914
5916 """
5917 Parser for SPINspiral output .
5918 """
5919 raise NotImplementedError
5920 if deltaLogP is not None:
5921 pos,bayesfactor=burnin(data,spin,deltaLogP,"posterior_samples.dat")
5922 return self._common_to_pos(open("posterior_samples.dat",'r'))
5923
5924 - def _ns_to_pos(self,files,Nlive=None,Npost=None,posfilename='posterior_samples.dat'):
5925 """
5926 Parser for nested sampling output.
5927 files : list of input NS files
5928 Nlive : Number of live points
5929 Npost : Desired number of posterior samples
5930 posfilename : Posterior output file name (default: 'posterior_samples.dat')
5931 """
5932 try:
5933 from lalapps.nest2pos import draw_N_posterior_many,draw_posterior_many
5934 except ImportError:
5935 print "Need lalapps.nest2pos to convert nested sampling output!"
5936 raise
5937
5938 if Nlive is None:
5939 raise RuntimeError("Need to specify number of live points in positional arguments of parse!")
5940
5941
5942
5943
5944 it = iter(files)
5945
5946
5947 parsfilename = (it.next()).strip('.gz')+'_params.txt'
5948
5949 if os.path.isfile(parsfilename):
5950 print 'Looking for '+parsfilename
5951
5952 if os.access(parsfilename,os.R_OK):
5953
5954 with open(parsfilename,'r') as parsfile:
5955 outpars=parsfile.readline()+'\n'
5956 else:
5957 raise RuntimeError('Cannot open parameters file %s!'%(parsfilename))
5958
5959 else:
5960 outpars='mchirp \t eta \t time \t phi0 \t dist \t RA \t \
5961 dec \t psi \t iota \t logl \n'
5962
5963
5964 parsvec=outpars.split()
5965 logLcol=-1
5966 for i in range(len(parsvec)):
5967 if parsvec[i].lower()=='logl':
5968 logLcol=i
5969 if logLcol==-1:
5970 print 'Error! Could not find logL column in parameter list: %s'%(outpars)
5971 raise RuntimeError
5972
5973 inarrays=map(np.loadtxt,files)
5974 if Npost is None:
5975 pos=draw_posterior_many(inarrays,[Nlive for f in files],logLcols=[logLcol for f in files])
5976 else:
5977 pos=draw_N_posterior_many(inarrays,[Nlive for f in files],Npost,logLcols=[logLcol for f in files])
5978
5979 with open(posfilename,'w') as posfile:
5980
5981 posfile.write(outpars)
5982
5983 for row in pos:
5984 for i in row:
5985 posfile.write('%10.12e\t' %(i))
5986 posfile.write('\n')
5987
5988 with open(posfilename,'r') as posfile:
5989 return_val=self._common_to_pos(posfile)
5990
5991 return return_val
5992
5994 """
5995 Parser for followupMCMC output.
5996 """
5997 return self._common_to_pos(open(files[0],'r'),delimiter=',')
5998
5999
6001 """
6002 Parser for MultiNest output.
6003 """
6004 return self._common_to_pos(open(files[0],'r'))
6005
6007 """
6008 Parser for VOTable XML Using
6009 """
6010 from xml.etree import ElementTree as ET
6011 xmlns='http://www.ivoa.net/xml/VOTable/v1.1'
6012 try:
6013 register_namespace=ET.register_namespace
6014 except AttributeError:
6015 def register_namespace(prefix,uri):
6016 ET._namespace_map[uri]=prefix
6017 register_namespace('vot',xmlns)
6018 tree = ET.ElementTree()
6019
6020 tree.parse(infile)
6021
6022 tables = tree.findall('.//{%s}TABLE'%(xmlns))
6023 for table in tables:
6024 if table.get('utype')=='lalinference:results:posteriorsamples':
6025 return(self._VOTTABLE2pos(table))
6026 for table in tables:
6027 if table.get('utype')=='lalinference:results:nestedsamples':
6028 nsresource=[node for node in tree.findall('{%s}RESOURCE'%(xmlns)) if node.get('utype')=='lalinference:results'][0]
6029 return(self._VOTTABLE2pos(vo_nest2pos(nsresource)))
6030 raise RuntimeError('Cannot find "Posterior Samples" TABLE element in XML input file %s'%(infile))
6031
6033 """
6034 Parser for a VOT TABLE element with FIELDs and TABLEDATA elements
6035 """
6036 from xml.etree import ElementTree as ET
6037 xmlns='http://www.ivoa.net/xml/VOTable/v1.1'
6038 try:
6039 register_namespace=ET.register_namespace
6040 except AttributeError:
6041 def register_namespace(prefix,uri):
6042 ET._namespace_map[uri]=prefix
6043 register_namespace('vot',xmlns)
6044
6045 header=[]
6046 for field in table.findall('./{%s}FIELD'%(xmlns)):
6047 header.append(field.attrib['name'])
6048 if(len(header)==0):
6049 raise RuntimeError('Unable to find FIELD nodes for table headers in XML table')
6050 data=table.findall('./{%s}DATA'%(xmlns))
6051 tabledata=data[0].find('./{%s}TABLEDATA'%(xmlns))
6052 llines=[]
6053 for row in tabledata:
6054 llines.append(np.array(map(lambda a:float(a.text),row)))
6055 flines=np.array(llines)
6056 for i in range(0,len(header)):
6057 if header[i].lower().find('log')!=-1 and header[i].lower() not in logParams and re.sub('log', '', header[i].lower()) not in [h.lower() for h in header]:
6058 print 'exponentiating %s'%(header[i])
6059
6060 flines[:,i]=np.exp(flines[:,i])
6061
6062 header[i]=re.sub('log', '', header[i], flags=re.IGNORECASE)
6063 if header[i].lower().find('sin')!=-1 and re.sub('sin', '', header[i].lower()) not in [h.lower() for h in header]:
6064 print 'asining %s'%(header[i])
6065 flines[:,i]=np.arcsin(flines[:,i])
6066 header[i]=re.sub('sin', '', header[i], flags=re.IGNORECASE)
6067 if header[i].lower().find('cos')!=-1 and re.sub('cos', '', header[i].lower()) not in [h.lower() for h in header]:
6068 print 'acosing %s'%(header[i])
6069 flines[:,i]=np.arccos(flines[:,i])
6070 header[i]=re.sub('cos', '', header[i], flags=re.IGNORECASE)
6071 header[i]=header[i].replace('(','')
6072 header[i]=header[i].replace(')','')
6073 print 'Read columns %s'%(str(header))
6074 return header,flines
6075
6076 - def _hdf5s_to_pos(self, infiles, fixedBurnins=None, deltaLogP=None, nDownsample=None, tablename=None, **kwargs):
6077 from astropy.table import vstack
6078
6079 if fixedBurnins is None:
6080 fixedBurnins = np.zeros(len(infiles))
6081
6082 if len(infiles) > 1:
6083 multiple_chains = True
6084
6085 chains = []
6086 for i, [infile, fixedBurnin] in enumerate(zip(infiles, fixedBurnins)):
6087 chain = self._hdf5_to_table(infile, fixedBurnin=fixedBurnin, deltaLogP=deltaLogP, nDownsample=nDownsample, multiple_chains=multiple_chains, tablename=tablename, **kwargs)
6088 chain.add_column(astropy.table.Column(i*np.ones(len(chain)), name='chain'))
6089 chains.append(chain)
6090
6091
6092 if deltaLogP is not None:
6093 logPThreshold = -np.inf
6094 for chain in chains:
6095 if len(chain) > 0:
6096 logPThreshold = max([logPThreshold, max(chain['logpost'])- deltaLogP])
6097 print("Eliminating any samples before log(L) = {}".format(logPThreshold))
6098
6099 for i, chain in enumerate(chains):
6100 if deltaLogP is not None:
6101 above_threshold = np.arange(len(chain))[chain['logpost'] > logPThreshold]
6102 burnin_idx = above_threshold[0] if len(above_threshold) > 0 else len(chain)
6103 else:
6104 burnin_idx = 0
6105 chains[i] = chain[burnin_idx:]
6106
6107 samples = vstack(chains)
6108
6109
6110 if nDownsample is not None:
6111 nskip = find_ndownsample(samples, nDownsample)
6112 samples = samples[::nskip]
6113
6114 return samples.colnames, as_array(samples).view(float).reshape(-1, len(samples.columns))
6115
6116 - def _hdf5_to_table(self, infile, deltaLogP=None, fixedBurnin=None, nDownsample=None, multiple_chains=False, tablename=None, **kwargs):
6117 """
6118 Parse a HDF5 file and return an array of posterior samples ad list of
6119 parameter names. Equivalent to '_common_to_pos' and work in progress.
6120 """
6121 if not tablename:
6122 samples = read_samples(infile, tablename=lalinference.LALInferenceHDF5PosteriorSamplesDatasetName)
6123 else:
6124 samples = read_samples(infile, tablename=tablename)
6125 params = samples.colnames
6126
6127 for param in params:
6128 param_low = param.lower()
6129 if param_low.find('log') != -1 and param_low not in logParams and re.sub('log', '', param_low) not in [p.lower() for p in params]:
6130 print('exponentiating %s' % param)
6131 new_param = re.sub('log', '', param, flags=re.IGNORECASE)
6132 samples[new_param] = np.exp(samples[param])
6133 del samples[param]
6134 param = new_param
6135 if param_low.find('sin') != -1 and re.sub('sin', '', param_low) not in [p.lower() for p in params]:
6136 print('asining %s' % param)
6137 new_param = re.sub('sin', '', param, flags=re.IGNORECASE)
6138 samples[new_param] = np.arcsin(samples[param])
6139 del samples[param]
6140 param = new_param
6141 if param_low.find('cos') != -1 and re.sub('cos', '', param_low) not in [p.lower() for p in params]:
6142 print('acosing %s' % param)
6143 new_param = re.sub('cos', '', param, flags=re.IGNORECASE)
6144 samples[new_param] = np.arccos(samples[param])
6145 del samples[param]
6146 param = new_param
6147
6148 if param != param.replace('(', ''):
6149 samples.rename_column(param, param.replace('(', ''))
6150 if param != param.replace(')', ''):
6151 samples.rename_column(param, param.replace(')', ''))
6152
6153
6154 replace_column(samples, param, samples[param].astype(float))
6155
6156 params = samples.colnames
6157 print('Read columns %s' % str(params))
6158
6159
6160 if 'cycle' in params:
6161 if not (fixedBurnin is None):
6162 if not (deltaLogP is None):
6163 print "Warning: using deltaLogP criteria in addition to fixed burnin"
6164 print "Fixed burning criteria: ",fixedBurnin
6165 else:
6166 fixedBurnin = 0
6167
6168 burned_in_cycles = np.arange(len(samples))[samples['cycle'] > fixedBurnin]
6169 burnin_idx = burned_in_cycles[0] if len(burned_in_cycles) > 0 else len(samples)
6170 samples = samples[burnin_idx:]
6171
6172 logPThreshold=-np.inf
6173 if len(samples) > 0 and not (deltaLogP is None):
6174 logPThreshold = max(samples['logpost'])- deltaLogP
6175 print "Eliminating any samples before log(post) = ", logPThreshold
6176 burnin_idx = np.arange(len(samples))[samples['logpost'] > logPThreshold][0]
6177 samples = samples[burnin_idx:]
6178
6179 if len(samples) > 0:
6180 nskip = find_ndownsample(samples, nDownsample)
6181 if nDownsample is None:
6182 print "Downsampling to take only uncorrelated posterior samples from each file."
6183 if np.isnan(nskip) and not multiple_chains:
6184 print "WARNING: All samples in chain are correlated. Downsampling to 10000 samples for inspection!!!"
6185 nskip = find_ndownsample(samples, 10000)
6186 samples = samples[::nskip]
6187 else:
6188 if np.isnan(nskip):
6189 print "WARNING: All samples in {} are correlated.".format(infile)
6190 samples = samples[-1:]
6191 else:
6192 print "Downsampling by a factor of ", nskip, " to achieve approximately ", nDownsample, " posterior samples"
6193 samples = samples[::nskip]
6194
6195 return samples
6196
6197 - def _hdf5_to_pos(self, infile, fixedBurnins=None, deltaLogP=None, nDownsample=None, tablename=None, **kwargs):
6201
6203 """
6204 Parse a file in the 'common format' and return an array of posterior
6205 samples and list of parameter names. Will apply inverse functions to
6206 columns with names containing sin,cos,log.
6207 """
6208
6209 [headerfile,delimiter]=info
6210
6211 if headerfile==None:
6212 formatstr=infile.readline().lstrip()
6213 else:
6214 hf=open(headerfile,'r')
6215 formatstr=hf.readline().lstrip()
6216 hf.close()
6217
6218 formatstr=formatstr.replace('#','')
6219 formatstr=formatstr.replace('"','')
6220
6221 header=formatstr.split(delimiter)
6222 header[-1]=header[-1].rstrip('\n')
6223 nparams=len(header)
6224 llines=[]
6225 dec=re.compile(r'^[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?$|^inf$')
6226
6227 for line_number,line in enumerate(infile):
6228 sline=line.split(delimiter)
6229 if sline[-1] == '\n':
6230 del(sline[-1])
6231 proceed=True
6232 if len(sline)<1:
6233 print 'Ignoring empty line in input file: %s'%(sline)
6234 proceed=False
6235 elif len(sline)!=nparams:
6236 sys.stderr.write('WARNING: Malformed row %i, read %i elements but there is meant to be %i\n'%(line_number,len(sline),nparams))
6237 proceed=False
6238
6239 for elemn,st in enumerate(sline):
6240 s=st.replace('\n','')
6241 if dec.search(s) is None:
6242 print 'Warning! Ignoring non-numeric data after the header: %s. Row = %i,Element=%i'%(s,line_number,elemn)
6243 proceed=False
6244 elif s is '\n':
6245 proceed=False
6246
6247 if proceed:
6248 llines.append(map(float,sline))
6249
6250 flines=np.array(llines)
6251
6252 if not flines.any():
6253 raise RuntimeError("ERROR: no lines read in!")
6254
6255
6256 for i in range(0,len(header)):
6257 if header[i].lower().find('log')!=-1 and header[i].lower() not in logParams and re.sub('log', '', header[i].lower()) not in [h.lower() for h in header]:
6258 print 'exponentiating %s'%(header[i])
6259
6260 flines[:,i]=np.exp(flines[:,i])
6261
6262 header[i]=re.sub('log', '', header[i], flags=re.IGNORECASE)
6263 if header[i].lower().find('sin')!=-1 and re.sub('sin', '', header[i].lower()) not in [h.lower() for h in header]:
6264 print 'asining %s'%(header[i])
6265 flines[:,i]=np.arcsin(flines[:,i])
6266 header[i]=re.sub('sin', '', header[i], flags=re.IGNORECASE)
6267 if header[i].lower().find('cos')!=-1 and re.sub('cos', '', header[i].lower()) not in [h.lower() for h in header]:
6268 print 'acosing %s'%(header[i])
6269 flines[:,i]=np.arccos(flines[:,i])
6270 header[i]=re.sub('cos', '', header[i], flags=re.IGNORECASE)
6271 header[i]=header[i].replace('(','')
6272 header[i]=header[i].replace(')','')
6273 print 'Read columns %s'%(str(header))
6274 return header,flines
6275
6278 result={}
6279 lines=fo.split('\n')
6280 chain_line=False
6281 for line in lines:
6282
6283 if '[1]' in line:
6284 key=line.replace('[1]','').strip(' ').strip('"')
6285 result[key]={}
6286 out=result[key]
6287 continue
6288 if result is not {}:
6289 if 'chain' in line:
6290 chain_line=True
6291 continue
6292 if chain_line:
6293 chain_line=False
6294 key=line.strip('"').split()[1]
6295 out[key]=[]
6296 out2=out[key]
6297 else:
6298 try:
6299 newline=line.strip('"').split()
6300 if newline is not []:
6301 out2.append(line.strip('"').split())
6302 except:
6303 pass
6304
6305 return result
6306
6309 """
6310 Parse a VO Table RESOURCE containing nested sampling output and
6311 return a VOTable TABLE element with posterior samples in it.
6312 This can be added to an existing tree by the user.
6313 Nlive will be read from the nsresource, unless specified
6314 """
6315 from xml.etree import ElementTree as ET
6316 import copy
6317 from math import log, exp
6318 xmlns='http://www.ivoa.net/xml/VOTable/v1.1'
6319 try:
6320 register_namespace=ET.register_namespace
6321 except AttributeError:
6322 def register_namespace(prefix,uri):
6323 ET._namespace_map[uri]=prefix
6324 register_namespace('vot',xmlns)
6325
6326 postable=ET.Element("{%s}TABLE"%(xmlns),attrib={'name':'Posterior Samples','utype':'lalinference:results:posteriorsamples'})
6327 i=0
6328 nstables=[resource for resource in nsresource.findall("./{%s}TABLE"%(xmlns)) if resource.get("utype")=="lalinference:results:nestedsamples"]
6329
6330 nstable=nstables[0]
6331 if Nlive is None:
6332 runstateResource = [resource for resource in nsresource.findall("./{%s}RESOURCE"%(xmlns)) if resource.get("utype")=="lalinference:state"][0]
6333 algTable = [table for table in runstateResource.findall("./{%s}TABLE"%(xmlns)) if table.get("utype")=="lalinference:state:algorithmparams"][0]
6334 Nlive = int ([param for param in algTable.findall("./{%s}PARAM"%(xmlns)) if param.get("name")=='Nlive'][0].get('value'))
6335 print 'Found Nlive %i'%(Nlive)
6336 if Nlive is None:
6337 raise RuntimeError("Cannot find number of live points in XML table, please specify")
6338 logLcol = None
6339 for fieldnode in nstable.findall('./{%s}FIELD'%xmlns):
6340 if fieldnode.get('name') == 'logL':
6341 logLcol=i
6342 i=i+1
6343 postable.append(copy.deepcopy(fieldnode))
6344 for paramnode in nstable.findall('./{%s}PARAM'%(xmlns)):
6345 postable.append(copy.deepcopy(paramnode))
6346 if logLcol is None:
6347 RuntimeError("Unable to find logL column")
6348 posdataNode=ET.Element("{%s}DATA"%(xmlns))
6349 postabledataNode=ET.Element("{%s}TABLEDATA"%(xmlns))
6350 postable.append(posdataNode)
6351 posdataNode.append(postabledataNode)
6352 nstabledata=nstable.find('./{%s}DATA/{%s}TABLEDATA'%(xmlns,xmlns))
6353 logw=log(1.0 - exp(-1.0/float(Nlive)))
6354 weights=[]
6355 for row in nstabledata:
6356 logL=float(row[logLcol].text)
6357 weights.append(logL-logw)
6358 logw=logw-1.0/float(Nlive)
6359 mw=max(weights)
6360 weights = [w - mw for w in weights]
6361 for (row,weight) in zip(nstabledata,weights):
6362 if weight > log(random.random()):
6363 postabledataNode.append(copy.deepcopy(row))
6364 return postable
6365
6366 xmlns='http://www.ivoa.net/xml/VOTable/v1.1'
6370 self.html=htmlSection("VOTable information")
6371 self.skiptable=0
6372
6373 - def start(self,tag,attrib):
6374 if tag=='{%s}TABLE'%(xmlns):
6375 if attrib['utype']=='lalinference:results:nestedsamples'\
6376 or attrib['utype']=='lalinference:results:posteriorsamples':
6377 self.skiptable=1
6378 else:
6379 self.skiptable=0
6380 self.tableouter=htmlChunk('div')
6381 self.tableouter.h2(attrib['name'])
6382 try:
6383 self.tableouter.p(attrib['utype'])
6384 except KeyError:
6385 pass
6386 self.fixedparams=htmlChunk('table',attrib={'class':'statstable'},parent=self.tableouter)
6387 self.table=htmlChunk('table',attrib={'class':'statstable'},parent=self.tableouter)
6388 self.tabheader=htmlChunk('tr',parent=self.table)
6389 if tag=='{%s}FIELD'%(xmlns):
6390 self.field=htmlChunk('th',{'name':attrib['name']},parent=self.tabheader)
6391 if tag=='{%s}TR'%(xmlns):
6392 self.tabrow=htmlChunk('tr',parent=self.table)
6393 if tag=='{%s}TD'%(xmlns):
6394 self.td=htmlChunk('td',parent=self.tabrow)
6395 if tag=='{%s}PARAM'%(xmlns):
6396 pnode=htmlChunk('tr',parent=self.fixedparams)
6397 namenode=htmlChunk('td',parent=pnode)
6398 namenode.p(attrib['name'])
6399 valnode=htmlChunk('td',parent=pnode)
6400 valnode.p(attrib['value'])
6401
6402 - def end(self,tag):
6403 if tag=='{%s}TABLE'%(xmlns):
6404 if not self.skiptable:
6405 self.html.append(self.tableouter._html)
6406 if tag=='{%s}FIELD'%(xmlns):
6407 self.field.p(self.data)
6408 if tag=='{%s}TD'%(xmlns):
6409 self.td.p(self.data)
6410
6411 - def data(self,data):
6413
6418 """Returns (high - low), the width of the given confidence
6419 bounds."""
6420
6421 return cl_bound[1] - cl_bound[0]
6422
6424 """Returns the number of samples within the given confidence
6425 bounds."""
6426
6427 return np.sum((samples >= cl_bound[0]) & (samples <= cl_bound[1]))
6428
6430 """Returns a tuple (relative_change, fractional_uncertainty,
6431 percentile_uncertainty) giving the uncertainty in confidence
6432 intervals from multiple posteriors.
6433
6434 The uncertainty in the confidence intervals is the difference in
6435 length between the widest interval, formed from the smallest to
6436 largest values among all the cl_bounds, and the narrowest
6437 interval, formed from the largest-small and smallest-large values
6438 among all the cl_bounds. Note that neither the smallest nor the
6439 largest confidence intervals necessarily correspond to one of the
6440 cl_bounds.
6441
6442 The relative change relates the confidence interval uncertainty to
6443 the expected value of the parameter, the fractional uncertainty
6444 relates this length to the length of the confidence level from the
6445 combined posteriors, and the percentile uncertainty gives the
6446 change in percentile over the combined posterior between the
6447 smallest and largest confidence intervals.
6448
6449 @param cl The confidence level (between 0 and 1).
6450
6451 @param cl_bounds A list of (low, high) pairs giving the confidence
6452 interval associated with each posterior.
6453
6454 @param posteriors A list of PosteriorOneDPDF objects giving the
6455 posteriors."""
6456
6457 Ns=[p.samples.shape[0] for p in posteriors]
6458 Nsamplers=len(Ns)
6459
6460
6461
6462 all_samples = np.squeeze(np.concatenate([p.samples for p in posteriors], axis=0))
6463 weights = np.squeeze(np.concatenate([p.samples*0.0+1.0/(Nsamplers*N) for (N,p) in zip(Ns,posteriors)], axis=0))
6464
6465 isort=np.argsort(all_samples)
6466
6467 all_samples = all_samples[isort]
6468 weights = weights[isort]
6469
6470 param_mean = np.average(all_samples, weights=weights)
6471
6472 N=all_samples.shape[0]
6473
6474 alpha = (1.0 - cl)/2.0
6475
6476 wttotal = np.cumsum(weights)
6477 ilow = np.nonzero(wttotal >= alpha)[0][0]
6478 ihigh = np.nonzero(wttotal >= 1.0-alpha)[0][0]
6479
6480 all_cl_bound = (all_samples[ilow], all_samples[ihigh])
6481
6482 low_bounds = np.array([l for (l,h) in cl_bounds])
6483 high_bounds = np.array([h for (l,h) in cl_bounds])
6484
6485 largest_cl_bound = (np.min(low_bounds), np.max(high_bounds))
6486 smallest_cl_bound = (np.max(low_bounds), np.min(high_bounds))
6487
6488 if smallest_cl_bound[1] < smallest_cl_bound[0]:
6489
6490 smallest_cl_bound = (0.0, 0.0)
6491
6492 ci_uncertainty = _cl_width(largest_cl_bound) - _cl_width(smallest_cl_bound)
6493
6494 relative_change = ci_uncertainty/param_mean
6495
6496 frac_uncertainty = ci_uncertainty/_cl_width(all_cl_bound)
6497
6498 quant_uncertainty = float(_cl_count(largest_cl_bound, all_samples) - _cl_count(smallest_cl_bound, all_samples))/float(N)
6499
6500 return (relative_change, frac_uncertainty, quant_uncertainty)
6501
6898
6899
6900 -def plot_psd(psd_files,outpath=None,f_min=30.):
6901 myfig2=plt.figure(figsize=(15,15),dpi=500)
6902 ax=plt.subplot(1,1,1)
6903 colors={'H1':'r','L1':'g','V1':'m','I1':'k','J1':'y'}
6904
6905 if outpath is None:
6906 outpath=os.getcwd()
6907 tmp=[]
6908 for f in psd_files:
6909 if not os.path.isfile(f):
6910 print "PSD file %s has not been found and won't be plotted\n"%f
6911 else:
6912 tmp.append(f)
6913 if tmp==[]:
6914 return None
6915 else:
6916 psd_files=tmp
6917
6918 freqs = {}
6919 for f in psd_files:
6920 data=np.loadtxt(f)
6921 freq=data[:,0]
6922 data=data[:,1]
6923 idx=f.find('-PSD.dat')
6924 ifo=f[idx-2:idx]
6925 freqs[ifo.lower()] = freq
6926 fr=[]
6927 da=[]
6928 for (f,d) in zip(freq,data):
6929 if f>f_min and d!=0.0 and np.isfinite(d):
6930 fr.append(f)
6931 da.append(d)
6932 plt.loglog(fr,da,colors[ifo],label=ifo,alpha=0.5,linewidth=3)
6933 plt.xlim([min(fr),max(fr)])
6934 plt.xlabel("Frequency [Hz]",fontsize=26)
6935 plt.ylabel("PSD",fontsize=26)
6936 plt.legend(loc='best')
6937 plt.grid(which='both')
6938 try:
6939 plt.tight_layout()
6940 myfig2.savefig(os.path.join(outpath,'PSD.png'),bbox_inches='tight')
6941 except:
6942 myfig2.savefig(os.path.join(outpath,'PSD.png'))
6943 myfig2.clf()
6944
6945 return freqs
6946
6947 cred_level = lambda cl, x: np.sort(x, axis=0)[int(cl*len(x))]
6950 """Return location of lower or upper confidence levels
6951 Args:
6952 x: List of samples.
6953 cl: Confidence level to return the bound of.
6954 lower: If ``True``, return the lower bound, otherwise return the upper bound.
6955 """
6956 if lower:
6957 return cred_level((1.-cl)/2, x)
6958 else:
6959 return cred_level((1.+cl)/2, x)
6960
6969
6970 -def plot_spline_pos(logf, ys, nf=100, level=0.9, color='k', label=None, xform=None):
6971 """Plot calibration posterior estimates for a spline model in log space.
6972 Args:
6973 logf: The (log) location of spline control points.
6974 ys: List of posterior draws of function at control points ``logf``
6975 nx: Number of points to evaluate spline at for plotting.
6976 level: Credible level to fill in.
6977 color: Color to plot with.
6978 label: Label for plot.
6979 xform: Function to transform the spline into plotted values.
6980 """
6981 f = np.exp(logf)
6982 fs = np.linspace(f.min(), f.max(), nf)
6983
6984 data = np.zeros((ys.shape[0], nf))
6985
6986 if xform is None:
6987 zs = ys
6988 else:
6989 zs = xform(ys)
6990
6991 mu = np.mean(zs, axis=0)
6992 lower_cl = mu - cred_interval(zs, level, lower=True)
6993 upper_cl = cred_interval(zs, level, lower=False) - mu
6994 plt.errorbar(np.exp(logf), mu, yerr=[lower_cl, upper_cl], fmt='.', color=color, lw=4, alpha=0.5, capsize=0)
6995
6996 for i, samp in enumerate(ys):
6997 temp = interpolate.spline(logf, samp, np.log(fs))
6998 if xform is None:
6999 data[i] = temp
7000 else:
7001 data[i] = xform(temp)
7002
7003 line, = plt.plot(fs, np.mean(data, axis=0), color=color, label=label)
7004 color = line.get_color()
7005 plt.fill_between(fs, cred_interval(data, level), cred_interval(data, level, lower=False), color=color, alpha=.1, linewidth=0.1)
7006 plt.xlim(f.min()-.5, f.max()+50)
7007
7009 fig, [ax1, ax2] = plt.subplots(2, 1, figsize=(15, 15), dpi=500)
7010
7011 font_size = 32
7012 if outpath is None:
7013 outpath=os.getcwd()
7014
7015 params = pos.names
7016 ifos = np.unique([param.split('_')[0] for param in params if 'spcal_freq' in param])
7017 for ifo in ifos:
7018 if ifo=='h1': color = 'r'
7019 elif ifo=='l1': color = 'g'
7020 elif ifo=='v1': color = 'm'
7021 else: color = 'c'
7022
7023
7024 freq_params = np.sort([param for param in params if
7025 '{0}_spcal_freq'.format(ifo) in param])
7026
7027 logfreqs = np.log([pos[param].median for param in freq_params])
7028
7029
7030 plt.sca(ax1)
7031 amp_params = np.sort([param for param in params if
7032 '{0}_spcal_amp'.format(ifo) in param])
7033 if len(amp_params) > 0:
7034 amp = 100*np.column_stack([pos[param].samples for param in amp_params])
7035 plot_spline_pos(logfreqs, amp, color=color, level=level, label="{0} (mean, {1}%)".format(ifo.upper(), int(level*100)))
7036
7037
7038 plt.sca(ax2)
7039 phase_params = np.sort([param for param in params if
7040 '{0}_spcal_phase'.format(ifo) in param])
7041 if len(phase_params) > 0:
7042 phase = np.column_stack([pos[param].samples for param in phase_params])
7043
7044 plot_spline_pos(logfreqs, phase, color=color, level=level, label="{0} (mean, {1}%)".format(ifo.upper(), int(level*100)), xform=spline_angle_xform)
7045
7046 ax1.tick_params(labelsize=.75*font_size)
7047 ax2.tick_params(labelsize=.75*font_size)
7048 try:
7049 plt.legend(loc='upper right', prop={'size':.75*font_size}, framealpha=0.1)
7050 except:
7051 plt.legend(loc='upper right', prop={'size':.75*font_size})
7052 ax1.set_xscale('log')
7053 ax2.set_xscale('log')
7054
7055 ax2.set_xlabel('Frequency (Hz)', fontsize=font_size)
7056 ax1.set_ylabel('Amplitude (%)', fontsize=font_size)
7057 ax2.set_ylabel('Phase (deg)', fontsize=font_size)
7058
7059 outp = os.path.join(outpath, 'calibration.png')
7060 try:
7061 fig.tight_layout()
7062 fig.savefig(outp, bbox_inches='tight')
7063 except:
7064 fig.savefig(outp)
7065 plt.close(fig)
7066
7106 def startElement(self,name,attrs):
7107 if attrs.has_key('Name') and table.Table.TableName(attrs['Name'])==self.tabname:
7108 self.tableElementName=name
7109
7110 ligolw.LIGOLWContentHandler.startElement(self,name,attrs)
7111 self.intable=True
7112 elif self.intable:
7113 ligolw.LIGOLWContentHandler.startElement(self,name,attrs)
7114 def endElement(self,name):
7115 if self.intable: ligolw.LIGOLWContentHandler.endElement(self,name)
7116 if self.intable and name==self.tableElementName: self.intable=False
7117
7118 lsctables.use_in(LIGOLWContentHandlerExtractSimBurstTable)
7119
7120
7121 srate=4096.0
7122 seglen=10.
7123 length=srate*seglen
7124 deltaT=1/srate
7125 deltaF = 1.0 / (length* deltaT);
7126
7127
7128 pad=0.4
7129 timeToFreqFFTPlan = CreateForwardREAL8FFTPlan(int(length), 1 );
7130 freqToTimeFFTPlan = CreateReverseREAL8FFTPlan(int(length), 1 );
7131 window=CreateTukeyREAL8Window(int(length),2.0*pad*srate/length);
7132
7133 segStart=939936910.000
7134 strainFinj= CreateCOMPLEX16FrequencySeries("strainF",segStart,0.0,deltaF,DimensionlessUnit,int(length/2. +1));
7135 strainTinj=CreateREAL8TimeSeries("strainT",segStart,0.0,1.0/srate,DimensionlessUnit,int(length));
7136 strainFrec= CreateCOMPLEX16FrequencySeries("strainF",segStart,0.0,deltaF,DimensionlessUnit,int(length/2. +1));
7137 strainTrec=CreateREAL8TimeSeries("strainT",segStart,0.0,1.0/srate,DimensionlessUnit,int(length));
7138 GlobREAL8time=None
7139 f_min=25
7140 f_ref=100
7141 f_max=srate/2.0
7142
7143 plot_fmax=2048
7144 plot_fmin=0.01
7145 plot_tmin=1e11
7146 plot_tmax=-1e11
7147
7148 inj_strains=dict((i,{"T":{'x':None,'strain':None},"F":{'x':None,'strain':None}}) for i in ifos)
7149 rec_strains=dict((i,{"T":{'x':None,'strain':None},"F":{'x':None,'strain':None}}) for i in ifos)
7150
7151 inj_domain=None
7152 rec_domain=None
7153 font_size=26
7154 if simburst is not None:
7155 skip=0
7156 try:
7157 xmldoc = utils.load_filename(simburst,contenthandler=LIGOLWContentHandlerExtractSimBurstTable)
7158 tbl = lsctables.SimBurstTable.get_table(xmldoc)
7159 if event>0:
7160 tbl=tbl[event]
7161 else:
7162 tbl=tbl[0]
7163 except:
7164 print "Cannot read event %s from table %s. Won't plot injected waveform \n"%(event,simburst)
7165 skip=1
7166 if not skip:
7167 REAL8time=tbl.time_geocent_gps+1e-9*tbl.time_geocent_gps_ns
7168 GPStime=LIGOTimeGPS(REAL8time)
7169 GlobREAL8time=(REAL8time)
7170 strainTinj.epoch=LIGOTimeGPS(round(GlobREAL8time,0)-seglen/2.)
7171 strainFinj.epoch=LIGOTimeGPS(round(GlobREAL8time,0)-seglen/2.)
7172 f0=tbl.frequency
7173 q=tbl.q
7174 dur=tbl.duration
7175 hrss=tbl.hrss
7176 polar_e_angle=tbl.pol_ellipse_angle
7177 polar_e_ecc=tbl.pol_ellipse_e
7178
7179 BurstExtraParams=None
7180 wf=str(tbl.waveform)
7181
7182 injapproximant=lalinf.GetBurstApproximantFromString(wf)
7183 ra=tbl.ra
7184 dec=tbl.dec
7185 psi=tbl.psi
7186
7187 if SimBurstImplementedFDApproximants(injapproximant):
7188 inj_domain='F'
7189 [plus,cross]=SimBurstChooseFDWaveform(deltaF, deltaT, f0, q,dur, f_min, f_max,hrss,polar_e_angle ,polar_e_ecc,BurstExtraParams, injapproximant)
7190 elif SimBurstImplementedTDApproximants(injapproximant):
7191 inj_domain='T'
7192 [plus,cross]=SimBurstChooseTDWaveform(deltaT, f0, q,dur, f_min, f_max,hrss,polar_e_angle ,polar_e_ecc,BurstExtraParams, injapproximant)
7193 else:
7194 print "\nThe approximant %s doesn't seem to be recognized by lalinference!\n Skipping WF plots\n"%injapproximant
7195 return None
7196
7197 for ifo in ifos:
7198 (fp,fc,fa,qv)=ant.response(REAL8time,ra,dec,0.0,psi,'radians',ifo)
7199 if inj_domain=='T':
7200
7201 tCinstrain=np.floor(REAL8time-float(strainTinj.epoch))/deltaT
7202
7203
7204 tSinstrain=int( (REAL8time-fabs(float(plus.epoch)) - fabs(float(strainTinj.epoch)))/deltaT)
7205 rem=(REAL8time-fabs(float(plus.epoch)) - fabs(float(strainTinj.epoch)))/deltaT-tSinstrain
7206
7207
7208 for k in np.arange(tSinstrain):
7209 strainTinj.data.data[k]=0.0
7210
7211 for k in np.arange(plus.data.length):
7212 strainTinj.data.data[k+tSinstrain]=((fp*plus.data.data[k]+fc*cross.data.data[k]))
7213
7214 for k in np.arange(strainTinj.data.length- (tSinstrain +plus.data.length)):
7215 strainTinj.data.data[k+tSinstrain+plus.data.length]=0.0
7216 for k in np.arange(strainTinj.data.length):
7217 strainTinj.data.data[k]*=window.data.data[k]
7218 np.savetxt('file.out',zip(np.array([strainTinj.epoch + j*deltaT for j in arange(strainTinj.data.length)]),np.array([strainTinj.data.data[j] for j in arange(strainTinj.data.length)])))
7219
7220 inj_strains[ifo]["T"]['strain']=np.array([strainTinj.data.data[j] for j in arange(strainTinj.data.length)])
7221 inj_strains[ifo]["T"]['x']=np.array([strainTinj.epoch + j*deltaT for j in arange(strainTinj.data.length)])
7222
7223 for j in arange(strainFinj.data.length):
7224 strainFinj.data.data[j]=0.0
7225 REAL8TimeFreqFFT(strainFinj,strainTinj,timeToFreqFFTPlan);
7226 twopit=2.*np.pi*(rem*deltaT)
7227 for k in arange(strainFinj.data.length):
7228 re = cos(twopit*deltaF*k)
7229 im = -sin(twopit*deltaF*k)
7230 strainFinj.data.data[k]*= (re + 1j*im);
7231
7232 inj_strains[ifo]["F"]['strain']=np.array([strainFinj.data.data[k] for k in arange(int(strainFinj.data.length))])
7233 inj_strains[ifo]["F"]['x']=np.array([strainFinj.f0+ k*strainFinj.deltaF for k in arange(int(strainFinj.data.length))])
7234 elif inj_domain=='F':
7235 for k in np.arange(strainFinj.data.length):
7236 if k<plus.data.length:
7237 strainFinj.data.data[k]=((fp*plus.data.data[k]+fc*cross.data.data[k]))
7238 else:
7239 strainFinj.data.data[k]=0.0
7240 twopit=2.*np.pi*(REAL8time-float(strainFinj.epoch))
7241 for k in arange(strainFinj.data.length):
7242 re = cos(twopit*deltaF*k)
7243 im = -sin(twopit*deltaF*k)
7244 strainFinj.data.data[k]*= (re + 1j*im);
7245
7246 inj_strains[ifo]["F"]['strain']=np.array([strainFinj.data.data[k] for k in arange(int(strainFinj.data.length))])
7247 inj_strains[ifo]["F"]['x']=np.array([strainFinj.f0+ k*strainFinj.deltaF for k in arange(int(strainFinj.data.length))])
7248
7249
7250 if f0 is not None and f0 is not np.nan:
7251 if q is not None and q is not np.nan:
7252 sigmaF=f0/q
7253 if f0-6.*sigmaF>plot_fmin:
7254 plot_fmin=f0-6.*sigmaF
7255 if f0+6.*sigmaF<plot_fmax:
7256 plot_fmax=f0+6.*sigmaF
7257 sigmaT=q/(2.*pi*f0)
7258 if REAL8time-6.*sigmaT<plot_tmin:
7259 plot_tmin=REAL8time-6.*sigmaT
7260 if REAL8time+6.*sigmaT>plot_tmax:
7261 plot_tmax=REAL8time+6.*sigmaT
7262
7263 if dur is not None and dur is not np.nan:
7264 sigmaF=1./sqrt(2.)/pi/dur
7265 if 0+6.*sigmaF<plot_fmax:
7266 plot_fmax=0+6.*sigmaF
7267 plot_fmin=0.0
7268 sigmaT=dur/sqrt(2.)
7269 if REAL8time-6.*sigmaT<plot_tmin:
7270 plot_tmin=REAL8time-6.*sigmaT
7271 if REAL8time+6.*sigmaT>plot_tmax:
7272 plot_tmax=REAL8time+6.*sigmaT
7273
7274
7275 if pos is not None:
7276
7277
7278 _,which=pos._posMap()
7279
7280 if 'time' in pos.names:
7281 REAL8time=pos['time'].samples[which][0]
7282 elif 'time_maxl' in pos.names:
7283 REAL8time=pos['time_maxl'].samples[which][0]
7284 elif 'time_mean' in pos.names:
7285 REAL8time=pos['time_mean'].samples[which][0]
7286 elif 'time_min' in pos.names and 'time_max' in pos.names:
7287 REAL8time=pos['time_min'].samples[which][0]+0.5*(pos['time_max'].samples[which][0]-pos['time_min'].samples[which][0])
7288 else:
7289 print "ERROR: could not find any time parameter in the posterior file. Not plotting the WF...\n"
7290 return None
7291
7292
7293 skip=0
7294
7295 try:
7296 approximant=int(pos['LAL_APPROXIMANT'].samples[which][0])
7297 except:
7298 skip=1
7299 if skip==0:
7300 GPStime=LIGOTimeGPS(REAL8time)
7301 if GlobREAL8time is None:
7302 GlobREAL8time=REAL8time
7303 strainTrec.epoch=LIGOTimeGPS(round(GlobREAL8time,0)-seglen/2.)
7304 strainFrec.epoch=LIGOTimeGPS(round(GlobREAL8time,0)-seglen/2.)
7305 if "duration" in pos.names:
7306 dur=pos["duration"].samples[which][0]
7307 else:
7308 dur=np.nan
7309 if "quality" in pos.names:
7310 q=pos['quality'].samples[which][0]
7311 else:
7312 q=np.nan
7313 if 'frequency' in pos.names:
7314 f0=pos['frequency'].samples[which][0]
7315 else:
7316 f0=np.nan
7317 try:
7318 hrss=pos['hrss'].samples[which][0]
7319 except:
7320 hrss=exp(pos['loghrss'].samples[which][0])
7321 if np.isnan(q) and not np.isnan(dur):
7322 q=sqrt(2)*pi*dur
7323 alpha=None
7324 if 'alpha' in pos.names:
7325 alpha=pos['alpha'].samples[which][0]
7326 polar_e_angle=alpha
7327 polar_e_ecc=pos['polar_eccentricity'].samples[which][0]
7328 elif 'polar_ellipse_angle' in pos.names:
7329 polar_e_angle=pos['polar_ellipse_angle'].samples[which][0]
7330 polar_e_ecc=pos['polar_eccentricity'].samples[which][0]
7331
7332 BurstExtraParams=None
7333
7334
7335
7336 if SimBurstImplementedFDApproximants(approximant):
7337 rec_domain='F'
7338 [plus,cross]=SimBurstChooseFDWaveform(deltaF, deltaT, f0, q,dur, f_min, f_max,hrss,polar_e_angle ,polar_e_ecc,BurstExtraParams, approximant)
7339 elif SimBurstImplementedTDApproximants(approximant):
7340 rec_domain='T'
7341 [plus,cross]=SimBurstChooseTDWaveform(deltaT, f0, q,dur, f_min, f_max,hrss,polar_e_angle ,polar_e_ecc,BurstExtraParams, approximant)
7342 else:
7343 print "The approximant %s doesn't seem to be recognized by lalinference!\n Skipping WF plots\n"%approximant
7344 return None
7345 ra=pos['ra'].samples[which][0]
7346 dec=pos['dec'].samples[which][0]
7347 psi=pos['psi'].samples[which][0]
7348 fs={}
7349 for ifo in ifos:
7350 (fp,fc,fa,qv)=ant.response(REAL8time,ra,dec,0.0,psi,'radians',ifo)
7351 if rec_domain=='T':
7352
7353 tCinstrain=np.floor(REAL8time-float(strainTrec.epoch))/deltaT
7354
7355 tSinstrain=int( (REAL8time-fabs(float(plus.epoch)) - fabs(float(strainTrec.epoch)))/deltaT)
7356
7357
7358 rem=(REAL8time-fabs(float(plus.epoch)) - fabs(float(strainTrec.epoch)))/deltaT-tSinstrain
7359
7360
7361
7362 for k in np.arange(tSinstrain):
7363 strainTrec.data.data[k]=0.0
7364
7365 for k in np.arange(plus.data.length):
7366 strainTrec.data.data[k+tSinstrain]=((fp*plus.data.data[k]+fc*cross.data.data[k]))
7367
7368 for k in np.arange(strainTrec.data.length- (tSinstrain +plus.data.length)):
7369 strainTrec.data.data[k+tSinstrain+plus.data.length]=0.0
7370 for k in np.arange(strainTrec.data.length):
7371 strainTrec.data.data[k]*=window.data.data[k]
7372
7373 rec_strains[ifo]["T"]['strain']=np.array([strainTrec.data.data[j] for j in arange(strainTrec.data.length)])
7374 rec_strains[ifo]["T"]['x']=np.array([strainTrec.epoch + j*deltaT for j in arange(strainTrec.data.length)])
7375
7376 for j in arange(strainFrec.data.length):
7377 strainFrec.data.data[j]=0.0
7378 REAL8TimeFreqFFT(strainFrec,strainTrec,timeToFreqFFTPlan);
7379 twopit=2.*np.pi*(rem*deltaT)
7380 for k in arange(strainFrec.data.length):
7381 re = cos(twopit*deltaF*k)
7382 im = -sin(twopit*deltaF*k)
7383 strainFrec.data.data[k]*= (re + 1j*im);
7384
7385 rec_strains[ifo]["F"]['strain']=np.array([strainFrec.data.data[k] for k in arange(int(strainFrec.data.length))])
7386 rec_strains[ifo]["F"]['x']=np.array([strainFrec.f0+ k*strainFrec.deltaF for k in arange(int(strainFrec.data.length))])
7387 elif rec_domain=='F':
7388 for k in np.arange(strainFrec.data.length):
7389 if k<plus.data.length:
7390 strainFrec.data.data[k]=((fp*plus.data.data[k]+fc*cross.data.data[k]))
7391 else:
7392 strainFrec.data.data[k]=0.0
7393 twopit=2.*np.pi*(REAL8time-float(strainFrec.epoch))
7394 for k in arange(strainFrec.data.length):
7395 re = cos(twopit*deltaF*k)
7396 im = -sin(twopit*deltaF*k)
7397 strainFrec.data.data[k]*= (re + 1j*im);
7398
7399 rec_strains[ifo]["F"]['strain']=np.array([strainFrec.data.data[k] for k in arange(int(strainFrec.data.length))])
7400 rec_strains[ifo]["F"]['x']=np.array([strainFrec.f0+ k*strainFrec.deltaF for k in arange(int(strainFrec.data.length))])
7401
7402
7403 if f0 is not None and f0 is not np.nan:
7404 if q is not None and q is not np.nan:
7405 sigmaF=f0/q
7406 if f0-6.*sigmaF>plot_fmin:
7407 plot_fmin=f0-6.*sigmaF
7408 if f0+6.*sigmaF<plot_fmax:
7409 plot_fmax=f0+6.*sigmaF
7410 sigmaT=q/(2.*pi*f0)
7411 if REAL8time-6.*sigmaT<plot_tmin:
7412 plot_tmin=REAL8time-6.*sigmaT
7413 if REAL8time+6.*sigmaT>plot_tmax:
7414 plot_tmax=REAL8time+6.*sigmaT
7415
7416 if dur is not None and dur is not np.nan:
7417 sigmaF=1./sqrt(2.)/pi/dur
7418 if 0+6.*sigmaF<plot_fmax:
7419 plot_fmax=0+6.*sigmaF
7420 plot_fmin=0.0
7421 sigmaT=dur/sqrt(2.)
7422 if REAL8time-6.*sigmaT<plot_tmin:
7423 plot_tmin=REAL8time-6.*sigmaT
7424 if REAL8time+6.*sigmaT>plot_tmax:
7425 plot_tmax=REAL8time+6.*sigmaT
7426
7427 myfig=plt.figure(1,figsize=(10,7))
7428
7429 rows=len(ifos)
7430 cols=2
7431
7432
7433
7434 global_domain="F"
7435 if rec_domain is not None and inj_domain is not None:
7436 if rec_domain=="T" and inj_domain=="T":
7437 global_domain="T"
7438 elif rec_domain is not None:
7439 if rec_domain=="T":
7440 global_domain="T"
7441 elif inj_domain is not None:
7442 if inj_domain=="T":
7443 global_domain="T"
7444
7445 A,axes=plt.subplots(nrows=rows,ncols=cols,sharex=False,sharey=False)
7446 plt.setp(A,figwidth=10,figheight=7)
7447 for (r,i) in zip(np.arange(rows),ifos):
7448 for c in np.arange(cols):
7449 ax=axes[r]
7450 if type(ax)==np.ndarray:
7451 ax=ax[c]
7452 else:
7453 ax=axes[c]
7454 if rec_strains[i]["T"]['strain'] is not None or rec_strains[i]["F"]['strain'] is not None:
7455 if c==0:
7456 if global_domain=="T":
7457 ax.plot(rec_strains[i]["T"]['x'],rec_strains[i]["T"]['strain'],colors_rec[i],label='%s maP'%i,linewidth=5)
7458 ax.set_xlim([plot_tmin,plot_tmax])
7459
7460 else:
7461 data=rec_strains[i]["F"]['strain']
7462 f=rec_strains[i]["F"]['x']
7463 mask=np.logical_and(f>=plot_fmin,f<=plot_fmax)
7464 ys=data
7465 ax.plot(f[mask],real(ys[mask]),'-',color=colors_rec[i],label='%s maP'%i,linewidth=5)
7466 ax.set_xlim([plot_fmin,plot_fmax])
7467 else:
7468 data=rec_strains[i]["F"]['strain']
7469 f=rec_strains[i]["F"]['x']
7470 mask=np.logical_and(f>=plot_fmin,f<=plot_fmax)
7471 ys=data
7472 ax.loglog(f[mask],absolute(ys[mask]),'--',color=colors_rec[i],linewidth=5)
7473 ax.grid(True,which='both')
7474 ax.set_xlim([plot_fmin,plot_fmax])
7475 if inj_strains[i]["T"]['strain'] is not None or inj_strains[i]["F"]['strain'] is not None:
7476 if c==0:
7477 if global_domain=="T":
7478 ax.plot(inj_strains[i]["T"]['x'],inj_strains[i]["T"]['strain'],colors_inj[i],label='%s inj'%i,linewidth=2)
7479 ax.set_xlim([plot_tmin,plot_tmax])
7480 else:
7481 data=inj_strains[i]["F"]['strain']
7482 f=inj_strains[i]["F"]['x']
7483 mask=np.logical_and(f>=plot_fmin,f<=plot_fmax)
7484 ys=data
7485 ax.plot(f[mask],real(ys[mask]),'-',color=colors_inj[i],label='%s inj'%i,linewidth=2)
7486 ax.set_xlim([plot_fmin,plot_fmax])
7487 else:
7488 data=inj_strains[i]["F"]['strain']
7489 f=inj_strains[i]["F"]['x']
7490 mask=np.logical_and(f>=plot_fmin,f<=plot_fmax)
7491 ys=data
7492 ax.loglog(f[mask],absolute(ys[mask]),'--',color=colors_inj[i],linewidth=2)
7493 ax.grid(True,which='both')
7494 ax.set_xlim([plot_fmin,plot_fmax])
7495 if r==0:
7496 if c==0:
7497 if global_domain=="T":
7498 ax.set_title(r"$h(t)$",fontsize=font_size)
7499 else:
7500 ax.set_title(r"$\Re[h(f)]$",fontsize=font_size)
7501 else:
7502 ax.set_title(r"$|h(f)|$",fontsize=font_size)
7503 elif r==rows-1:
7504 if c==0:
7505 if global_domain=="T":
7506 ax.set_xlabel("time [s]",fontsize=font_size)
7507 else:
7508 ax.set_xlabel("frequency [Hz]",fontsize=font_size)
7509 else:
7510 ax.set_xlabel("frequency [Hz]",fontsize=font_size)
7511
7512 ax.legend(loc='best')
7513 ax.grid(True)
7514
7515
7516 A.savefig(os.path.join(path,'WF_DetFrame.png'),bbox_inches='tight')
7517 return inj_strains, rec_strains
7518
7519 -def make_1d_table(html,legend,label,pos,pars,noacf,GreedyRes,onepdfdir,sampsdir,savepdfs,greedy,analyticLikelihood,nDownsample):
7520
7521 from numpy import unique, sort
7522 global confidenceLevels
7523 confidence_levels=confidenceLevels
7524
7525 out={}
7526 if pars==[]:
7527 return out
7528 if set(pos.names)-set(pars)==set(pos.names):
7529 return out
7530
7531
7532 tabid='onedmargtable_'+label.lower()
7533 html_ompdf=html.add_collapse_section('1D marginal posterior PDFs (%s)'%label,legend=legend,innertable_id=tabid)
7534
7535 if not noacf:
7536 html_ompdf_write= '<table id="%s"><tr><th>Histogram and Kernel Density Estimate</th><th>Samples used</th><th>Autocorrelation</th></tr>'%tabid
7537 else:
7538 html_ompdf_write= '<table id="%s"><tr><th>Histogram and Kernel Density Estimate</th><th>Samples used</th></tr>'%tabid
7539
7540 Nskip=0
7541 if 'chain' in pos.names:
7542 data,header=pos.samples()
7543 par_index=pos.names.index('cycle')
7544 chain_index=pos.names.index("chain")
7545 chains=unique(pos["chain"].samples)
7546 chainCycles = [sort(data[ data[:,chain_index] == chain, par_index ]) for chain in chains]
7547 chainNcycles = []
7548 chainNskips = []
7549 for cycles in chainCycles:
7550 if len(cycles) > 1:
7551 chainNcycles.append(cycles[-1] - cycles[0])
7552 chainNskips.append(cycles[1] - cycles[0])
7553 else:
7554 chainNcycles.append(1)
7555 chainNskips.append(1)
7556 elif 'cycle' in pos.names:
7557 cycles = sort(pos['cycle'].samples)
7558 if len(cycles) > 1:
7559 Ncycles = cycles[-1]-cycles[0]
7560 Nskip = cycles[1]-cycles[0]
7561 else:
7562 Ncycles = 1
7563 Nskip = 1
7564
7565 printed=0
7566 for par_name in pars:
7567 par_name=par_name.lower()
7568 try:
7569 pos[par_name.lower()]
7570 except KeyError:
7571
7572 continue
7573 try:
7574 par_bin=GreedyRes[par_name]
7575 except KeyError:
7576 print "Bin size is not set for %s, skipping binning."%par_name
7577 continue
7578
7579
7580 binParams={par_name:par_bin}
7581 injection_area=None
7582 injection_area=None
7583 if greedy:
7584 if printed==0:
7585 print "Using greedy 1-d binning credible regions\n"
7586 printed=1
7587 toppoints,injectionconfidence,reses,injection_area,cl_intervals=greedy_bin_one_param(pos,binParams,confidence_levels)
7588 else:
7589 if printed==0:
7590 print "Using 2-step KDE 1-d credible regions\n"
7591 printed=1
7592 if pos[par_name].injval is None:
7593 injCoords=None
7594 else:
7595 injCoords=[pos[par_name].injval]
7596 _,reses,injstats=kdtree_bin2Step(pos,[par_name],confidence_levels,injCoords=injCoords)
7597 if injstats is not None:
7598 injectionconfidence=injstats[3]
7599 injection_area=injstats[4]
7600
7601 print "Generating 1D plot for %s."%par_name
7602 out[par_name]=reses
7603
7604 pdf=cdf=None
7605 if analyticLikelihood:
7606 pdf = analyticLikelihood.pdf(par_name)
7607 cdf = analyticLikelihood.cdf(par_name)
7608
7609 oneDPDFParams={par_name:50}
7610 try:
7611 rbins,plotFig=plot_one_param_pdf(pos,oneDPDFParams,pdf,cdf,plotkde=False)
7612 except:
7613 print "Failed to produce plot for %s."%par_name
7614 continue
7615
7616 figname=par_name+'.png'
7617 oneDplotPath=os.path.join(onepdfdir,figname)
7618 plotFig.savefig(oneDplotPath)
7619 if(savepdfs): plotFig.savefig(os.path.join(onepdfdir,par_name+'.pdf'))
7620 plt.close(plotFig)
7621
7622 if rbins:
7623 print "r of injected value of %s (bins) = %f"%(par_name, rbins)
7624
7625
7626 myfig=plt.figure(figsize=(4,3.5),dpi=200)
7627 pos_samps=pos[par_name].samples
7628 if not ("chain" in pos.names):
7629
7630
7631 plt.plot(pos_samps,'k.', markersize=5, alpha=0.5, linewidth=0.0, figure=myfig)
7632 maxLen=len(pos_samps)
7633 else:
7634
7635
7636
7637 data,header=pos.samples()
7638 par_index=pos.names.index(par_name)
7639 chain_index=pos.names.index("chain")
7640 chains=unique(pos["chain"].samples)
7641 chainData=[data[ data[:,chain_index] == chain, par_index ] for chain in chains]
7642 chainDataRanges=[range(len(cd)) for cd in chainData]
7643 maxLen=max([len(cd) for cd in chainData])
7644 for rng, data in zip(chainDataRanges, chainData):
7645 plt.plot(rng, data, marker='.', markersize=1, alpha=0.5, linewidth=0.0,figure=myfig)
7646 plt.title("Gelman-Rubin R = %g"%(pos.gelman_rubin(par_name)))
7647
7648
7649
7650
7651
7652
7653
7654 injpar=pos[par_name].injval
7655
7656 if injpar is not None:
7657
7658 minrange=min(pos_samps)-0.05*(max(pos_samps)-min(pos_samps))
7659 maxrange=max(pos_samps)+0.05*(max(pos_samps)-min(pos_samps))
7660 if minrange<injpar and maxrange>injpar:
7661 plt.axhline(injpar, color='r', linestyle='-.',linewidth=4)
7662 myfig.savefig(os.path.join(sampsdir,figname.replace('.png','_samps.png')))
7663 if(savepdfs): myfig.savefig(os.path.join(sampsdir,figname.replace('.png','_samps.pdf')))
7664 plt.close(myfig)
7665 acfail=0
7666 if not (noacf):
7667 acffig=plt.figure(figsize=(4,3.5),dpi=200)
7668 if not ("chain" in pos.names):
7669 data=pos_samps[:,0]
7670 try:
7671 (Neff, acl, acf) = effectiveSampleSize(data, Nskip)
7672 lines=plt.plot(acf, 'k.', marker='.', markersize=1, alpha=0.5, linewidth=0.0, figure=acffig)
7673
7674 if nDownsample is None:
7675 plt.title('Autocorrelation Function')
7676 elif 'cycle' in pos.names:
7677 last_color = lines[-1].get_color()
7678 plt.axvline(acl/Nskip, linestyle='-.', color=last_color)
7679 plt.title('ACL = %i N = %i'%(acl,Neff))
7680 except FloatingPointError:
7681
7682 acfail=1
7683 pass
7684 else:
7685 try:
7686 acls = []
7687 Nsamps = 0.0;
7688 for rng, data, Nskip, Ncycles in zip(chainDataRanges, chainData, chainNskips, chainNcycles):
7689 (Neff, acl, acf) = effectiveSampleSize(data, Nskip)
7690 acls.append(acl)
7691 Nsamps += Neff
7692 lines=plt.plot(acf,'k.', marker='.', markersize=1, alpha=0.5, linewidth=0.0, figure=acffig)
7693
7694 if nDownsample is not None:
7695 last_color = lines[-1].get_color()
7696 plt.axvline(acl/Nskip, linestyle='-.', color=last_color)
7697 if nDownsample is None:
7698 plt.title('Autocorrelation Function')
7699 else:
7700 plt.title('ACL = %i N = %i'%(max(acls),Nsamps))
7701 except FloatingPointError:
7702
7703 acfail=1
7704 pass
7705
7706 acffig.savefig(os.path.join(sampsdir,figname.replace('.png','_acf.png')))
7707 if(savepdfs): acffig.savefig(os.path.join(sampsdir,figname.replace('.png','_acf.pdf')))
7708 plt.close(acffig)
7709
7710 if not noacf:
7711 if not acfail:
7712 acfhtml='<td width="30%"><img width="100%" src="1Dsamps/'+figname.replace('.png', '_acf.png')+'"/></td>'
7713 else:
7714 acfhtml='<td>ACF generation failed!</td>'
7715 html_ompdf_write+='<tr><td width="30%"><img width="100%" src="1Dpdf/'+figname+'"/></td><td width="30%"><img width="100%" src="1Dsamps/'+figname.replace('.png','_samps.png')+'"/></td>'+acfhtml+'</tr>'
7716 else:
7717 html_ompdf_write+='<tr><td width="30%"><img width="100%" src="1Dpdf/'+figname+'"/></td><td width="30%"><img width="100%" src="1Dsamps/'+figname.replace('.png','_samps.png')+'"/></td></tr>'
7718
7719 html_ompdf_write+='</table>'
7720 html_ompdf.write(html_ompdf_write)
7721
7722 return out
7723