37 from six.moves
import cPickle
as pickle
39 from time
import strftime
43 from numpy
import (exp, cos, sin, size, cov, unique, hsplit, log, squeeze)
47 from matplotlib
import pyplot
as plt
51 inches_per_pt = 1.0/72.27
52 golden_mean = (2.236-1.0)/2.0
53 fig_width = fig_width_pt*inches_per_pt
54 fig_height = fig_width*golden_mean
55 fig_size = [fig_width,fig_height]
56 matplotlib.rcParams.update(
57 {
'axes.labelsize': 16,
59 'legend.fontsize': 16,
60 'xtick.labelsize': 16,
61 'ytick.labelsize': 16,
63 'figure.figsize': fig_size,
64 'font.family':
"serif",
65 'font.serif': [
'Times',
'Palatino',
'New Century Schoolbook',
'Bookman',
'Computer Modern Roman',
'Times New Roman',
'Liberation Serif'],
66 'font.weight':
'normal',
73 from lalinference
import bayespputils
as bppu
74 from lalinference
import git_version
76 from ligo.lw
import ligolw
77 from ligo.lw
import lsctables
78 from ligo.lw
import utils
80 __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>"
81 __version__=
"git id %s"%git_version.id
82 __date__= git_version.date
88 USER = os.environ[
'USER']
89 HOST = socket.getfqdn()
90 address=address.split(
',')
92 SUBJECT=
"LALInference result is ready at "+HOST+
"!"
94 fslocation=os.path.abspath(path)
95 webpath=
'posplots.html'
96 url =
guess_url(os.path.join(fslocation,webpath))
97 TEXT=
"Hi "+USER+
",\nYou have a new parameter estimation result on "+HOST+
".\nYou can view the result at "+url+
"\n"
98 cmd=
'echo "%s" | mail -s "%s" "%s"'%(TEXT,SUBJECT,
', '.join(address))
99 proc = subprocess.Popen(cmd, stdout=subprocess.PIPE,stderr=subprocess.STDOUT, shell=
True)
100 (out, err) = proc.communicate()
105 Pickle/serialize 'obj' into 'fname'.
107 filed=open(fname,
'w')
108 pickle.dump(obj,filed)
112 filed=open(fname,
'w')
113 for key,value
in dict.items():
114 filed.write(
"%s %s\n"%(
str(key),
str(value)) )
126 for arg
in parser.rargs:
128 if arg[:2] ==
"--" and len(arg) > 2:
131 if arg[:1] ==
"-" and len(arg) > 1
and not floatable(arg):
135 del parser.rargs[:len(args)]
137 if getattr(parser.values, opt.dest):
138 oldargs = getattr(parser.values, opt.dest)
141 setattr(parser.values, opt.dest, args)
145 out=bppu.htmlChunk(
'div',parent=parent)
147 row=out.insert_row(tab)
149 out.insert_td(row,
str(key))
150 row2=out.insert_row(tab)
151 for val
in d.values():
152 out.insert_td(row2,
str(val))
158 sec=bppu.htmlSection(h5grp.name,htmlElement=parent)
161 if(isinstance(h5grp[group],h5py.Group)):
167 outdir,data,oneDMenus,twoDGreedyMenu,GreedyRes,
168 confidence_levels,twoDplots,
170 injfile=None,eventnum=None,
171 trigfile=None,trignum=None,
174 dievidence=False,boxing=64,difactor=1.0,
178 bayesfactornoise=None,bayesfactorcoherent=None,
182 ns_flag=False,ns_Nlive=None,
184 ss_flag=False,ss_spin_flag=False,
186 li_flag=False,deltaLogP=None,fixedBurnins=None,nDownsample=None,oldMassConvention=False,
190 inj_spin_frame='OrbitalL',
196 RconvergenceTests=False,
200 covarianceMatrices=None,
209 This is a demonstration script for using the functionality/data structures
210 contained in lalinference.bayespputils . It will produce a webpage from a file containing
211 posterior samples generated by the parameter estimation codes with 1D/2D plots
212 and stats from the marginal posteriors for each parameter/set of parameters.
214 if eventnum
is not None and injfile
is None:
215 print(
"You specified an event number but no injection file. Ignoring!")
217 if trignum
is not None and trigfile
is None:
218 print(
"You specified a trigger number but no trigger file. Ignoring!")
220 if trignum
is None and trigfile
is not None:
221 print(
"You specified a trigger file but no trigger number. Taking first entry (the case for GraceDB events).")
225 raise RuntimeError(
'You must specify an input data file')
228 raise RuntimeError(
"You must specify an output directory.")
230 if not os.path.isdir(outdir):
234 peparser=bppu.PEOutputParser(
'fm')
235 commonResultsObj=peparser.parse(data)
237 elif ns_flag
and not ss_flag:
238 peparser=bppu.PEOutputParser(
'ns')
239 commonResultsObj=peparser.parse(data,Nlive=ns_Nlive)
241 elif ss_flag
and not ns_flag:
242 peparser=bppu.PEOutputParser(
'mcmc_burnin')
243 commonResultsObj=peparser.parse(data,spin=ss_spin_flag,deltaLogP=deltaLogP)
246 peparser=bppu.PEOutputParser(
'inf_mcmc')
247 commonResultsObj=peparser.parse(data,outdir=outdir,deltaLogP=deltaLogP,fixedBurnins=fixedBurnins,nDownsample=nDownsample,oldMassConvention=oldMassConvention)
249 elif ss_flag
and ns_flag:
250 raise RuntimeError(
"Undefined input format. Choose only one of:")
252 elif '.hdf' in data[0]
or '.h5' in data[0]:
254 peparser = bppu.PEOutputParser(
'hdf5s')
255 commonResultsObj=peparser.parse(data,deltaLogP=deltaLogP,fixedBurnins=fixedBurnins,nDownsample=nDownsample)
257 fixedBurnins = fixedBurnins
if fixedBurnins
is not None else None
258 peparser = bppu.PEOutputParser(
'hdf5')
259 commonResultsObj=peparser.parse(data[0],deltaLogP=deltaLogP,fixedBurnins=fixedBurnins,nDownsample=nDownsample)
261 peparser=bppu.PEOutputParser(
'common')
262 commonResultsObj=peparser.parse(open(data[0],
'r'),info=[header,
None])
264 if os.path.isfile(data[0]+
"_header.txt"):
266 shutil.copy2(data[0]+
"_header.txt", os.path.join(outdir,
'nest_headers.txt'))
271 ps,samps = commonResultsObj
273 f_refIdx = ps.index(
'f_ref')
274 injFref = unique(samps[:,f_refIdx])
276 print(
"ERROR: Expected f_ref to be constant for all samples. Can't tell which value was injected!")
286 if injfile
and eventnum
is not None:
287 print(
'Looking for event %i in %s\n'%(eventnum,injfile))
288 xmldoc = utils.load_filename(injfile,contenthandler=lsctables.use_in(ligolw.LIGOLWContentHandler))
289 siminspiraltable=lsctables.SimInspiralTable.get_table(xmldoc)
290 injection=siminspiraltable[eventnum]
294 if trigfile
is not None and trignum
is not None:
295 triggers = bppu.readCoincXML(trigfile, trignum)
299 if bayesfactornoise
is not None:
300 bfile=open(bayesfactornoise,
'r')
303 if(len(BSN.split())!=1):
306 if bayesfactorcoherent
is not None:
307 bfile=open(bayesfactorcoherent,
'r')
312 if snrfactor
is not None:
313 if not os.path.isfile(snrfactor):
314 print(
"No snr file provided or wrong path to snr file\n")
318 snrfile=open(snrfactor,
'r')
319 snrs=snrfile.readlines()
324 snrstring=snrstring +
" "+
str(snr[0:-1])+
" ,"
325 snrstring=snrstring[0:-1]
329 pos = bppu.Posterior(commonResultsObj,SimInspiralTableEntry=injection,inj_spin_frame=inj_spin_frame, injFref=injFref,SnglInspiralList=triggers)
332 analyticLikelihood =
None
333 if covarianceMatrices
and meanVectors:
334 analyticLikelihood = bppu.AnalyticLikelihood(covarianceMatrices, meanVectors)
337 oneDMenu = analyticLikelihood.names
340 for i
in range(len(oneDMenu)):
341 for j
in range(i+1,len(oneDMenu)):
342 twoDGreedyMenu.append([oneDMenu[i],oneDMenu[j]])
343 twoDplots = twoDGreedyMenu
345 if eventnum
is None and injfile
is not None:
347 injections = lsctables.SimInspiralTable.get_table(utils.load_filename(injfile, contenthandler=lsctables.use_in(ligolw.LIGOLWContentHandler)))
349 if(len(injections)<1):
351 print(
'Warning: Cannot find injection with end time %f' %(pos[
'time'].mean))
353 print(
"Warning: No 'time' column!")
357 injection = itertools.ifilter(
lambda a: abs(float(a.get_end()) - pos[
'time'].mean) < 0.1, injections).next()
358 pos.set_injection(injection)
360 print(
"Warning: No 'time' column!")
362 pos.extend_posterior()
365 for i
in oneDMenus.keys():
366 oneDMenu+=oneDMenus[i]
368 functions = {
'cos':cos,
'sin':sin,
'exp':exp,
'log':log}
369 for pos_name
in oneDMenu:
370 if pos_name
not in pos.names:
371 for func
in functions.keys():
372 old_pos_name = pos_name.replace(func,
'')
373 if pos_name.find(func)==0
and old_pos_name
in pos.names:
374 print(
"Taking %s of %s ..."% (func,old_pos_name))
375 pos.append_mapping(pos_name,functions[func],old_pos_name)
378 requested_params = set(pos.names).intersection(set(oneDMenu))
379 pos.delete_NaN_entries(requested_params)
382 if analyticLikelihood:
383 dievidence_names = [
'post',
'posterior',
'logl',
'prior',
'likelihood',
'cycle',
'chain']
384 [pos.pop(param)
for param
in pos.names
if param
not in analyticLikelihood.names
and param
not in dievidence_names]
388 print(
"Number of posterior samples: %i"%len(pos))
391 print(
str(pos.means))
394 print(
str(pos.medians))
397 max_pos,max_pos_co=pos.maxL
401 posfilename=os.path.join(outdir,
'posterior_samples.dat')
402 pos.write_to_file(posfilename)
408 html=bppu.htmlPage(
'Posterior PDFs',css=bppu.__default_css_string,javascript=bppu.__default_javascript_string)
411 html_meta=html.add_section(
'Summary')
412 table=html_meta.tab()
413 row=html_meta.insert_row(table,label=
'thisrow')
414 td=html_meta.insert_td(row,
'',label=
'Samples')
415 SampsStats=html.add_section_to_element(
'Samples',td)
416 SampsStats.p(
'Produced from '+
str(len(pos))+
' posterior samples.')
417 if 'chain' in pos.names:
418 acceptedChains = unique(pos[
'chain'].samples)
419 acceptedChainText =
'%i of %i chains accepted: %i'%(len(acceptedChains),len(data),acceptedChains[0])
420 if len(acceptedChains) > 1:
421 for chain
in acceptedChains[1:]:
422 acceptedChainText +=
', %i'%(chain)
423 SampsStats.p(acceptedChainText)
424 if 'cycle' in pos.names:
425 SampsStats.p(
'Longest chain has '+
str(pos.longest_chain_cycles())+
' cycles.')
426 filenames=
'Samples read from %s'%(data[0])
428 for fname
in data[1:]:
429 filenames+=
', '+
str(fname)
430 SampsStats.p(filenames)
431 td=html_meta.insert_td(row,
'',label=
'SummaryLinks')
432 legend=html.add_section_to_element(
'Sections',td)
435 if '.h5' in data[0]
or '.hdf' in data[0]:
436 html_hdf=html.add_section(
'Metadata',legend=legend)
438 with h5py.File(data[0],
'r')
as h5grp:
442 if bayesfactorcoherent
is not None or bayesfactornoise
is not None:
443 html_model=html.add_section(
'Model selection',legend=legend)
444 if bayesfactornoise
is not None:
445 html_model.p(
'log Bayes factor ( coherent vs gaussian noise) = %s, Bayes factor=%f'%(BSN,exp(float(BSN))))
446 if bayesfactorcoherent
is not None:
447 html_model.p(
'log Bayes factor ( coherent vs incoherent OR noise ) = %s, Bayes factor=%f'%(BCI,exp(float(BCI))))
450 html_model=html.add_section(
'Direct Integration Evidence',legend=legend)
451 log_ev = log(difactor) + pos.di_evidence(boxing=boxing)
453 evfilename=os.path.join(outdir,
"evidence.dat")
454 evout=open(evfilename,
"w")
457 evout.write(
str(log_ev))
459 print(
"Computing direct integration evidence = %g (log(Evidence) = %g)"%(ev, log_ev))
460 html_model.p(
'Direct integration evidence is %g, or log(Evidence) = %g. (Boxing parameter = %d.)'%(ev,log_ev,boxing))
461 if 'logl' in pos.names:
462 log_ev=pos.harmonic_mean_evidence()
463 html_model.p(
'Compare to harmonic mean evidence of %g (log(Evidence) = %g).'%(exp(log_ev),log_ev))
467 html_model=html.add_section(
'Elliptical Evidence',legend=legend)
468 log_ev = pos.elliptical_subregion_evidence()
470 evfilename=os.path.join(outdir,
'ellevidence.dat')
471 evout=open(evfilename,
'w')
472 evout.write(
str(ev) +
' ' +
str(log_ev))
474 print(
'Computing elliptical region evidence = %g (log(ev) = %g)'%(ev, log_ev))
475 html_model.p(
'Elliptical region evidence is %g, or log(Evidence) = %g.'%(ev, log_ev))
477 if 'logl' in pos.names:
478 log_ev=pos.harmonic_mean_evidence()
479 html_model.p(
'Compare to harmonic mean evidence of %g (log(Evidence = %g))'%(exp(log_ev), log_ev))
481 print(
'Warning: Sample size too small to compute elliptical evidence!')
484 if snrfactor
is not None:
485 html_snr=html.add_section(
'Signal to noise ratio(s)',legend=legend)
486 html_snr.p(
'%s'%snrstring)
489 html_dic = html.add_section(
'Deviance Information Criterion', legend=legend)
490 html_dic.p(
'DIC = %.1f'%pos.DIC)
491 dicout = open(os.path.join(outdir,
'dic.dat'),
'w')
493 dicout.write(
'%.1f\n'%pos.DIC)
500 html_stats=html.add_collapse_section(
'Summary statistics',legend=legend,innertable_id=tabid,start_closed=
False)
501 html_stats.write(
str(pos))
502 statfilename=os.path.join(outdir,
"summary_statistics.dat")
503 statout=open(statfilename,
"w")
504 statout.write(
"\tmaP\tmaxL\tKL\tstdev\tmean\tmedian\tstacc\tinjection\tvalue\n")
506 warned_about_kl =
False
507 for statname,statoned_pos
in pos:
509 statmax_pos,max_i=pos._posMaxL()
510 statmaxL=statoned_pos.samples[max_i][0]
512 statKL = statoned_pos.KL()
514 if not warned_about_kl:
515 print(
"Unable to compute KL divergence")
516 warned_about_kl =
True
519 statmax_pos,max_j=pos._posMap()
520 statmaxP=statoned_pos.samples[max_j][0]
521 statmean=
str(statoned_pos.mean)
522 statstdev=
str(statoned_pos.stdev)
523 statmedian=
str(squeeze(statoned_pos.median))
524 statstacc=
str(statoned_pos.stacc)
525 statinjval=
str(statoned_pos.injval)
527 statarray=[
str(i)
for i
in [statname,statmaxP,statmaxL,statKL,statstdev,statmean,statmedian,statstacc,statinjval]]
528 statout.write(
"\t".join(statarray))
538 sky_injection_cl=
None
541 html_wf=html.add_collapse_section(
'Sky Localization and Waveform',innertable_id=tabid)
543 table=html_wf.tab(idtable=tabid)
544 row=html_wf.insert_row(table,label=
'SkyandWF')
545 skytd=html_wf.insert_td(row,
'',label=
'SkyMap',legend=legend)
546 html_sky=html.add_section_to_element(
'SkyMap',skytd)
548 if skyres
is not None and \
549 (
'ra' in pos.names
and 'dec' in pos.names):
551 if pos[
'dec'].injval
is not None and pos[
'ra'].injval
is not None:
552 inj_position=[pos[
'ra'].injval,pos[
'dec'].injval]
556 hpmap = pos.healpix_map(float(skyres), nest=
True)
557 bppu.plot_sky_map(hpmap, outdir, inj=inj_position, nest=
True)
559 if inj_position
is not None:
560 html_sky.p(
'Injection found at p = %g'%bppu.skymap_inj_pvalue(hpmap, inj_position, nest=
True))
562 html_sky.write(
'<a href="skymap.png" target="_blank"><img src="skymap.png"/></a>')
564 html_sky_write=
'<table border="1" id="statstable"><tr><th>Confidence region</th><th>size (sq. deg)</th></tr>'
566 areas = bppu.skymap_confidence_areas(hpmap, confidence_levels)
567 for cl, area
in zip(confidence_levels, areas):
568 html_sky_write+=
'<tr><td>%g</td><td>%g</td></tr>'%(cl, area)
569 html_sky_write+=(
'</table>')
571 html_sky.write(html_sky_write)
573 html_sky.write(
'<b>No skymap generated!</b>')
574 wfdir=os.path.join(outdir,
'Waveform')
575 if not os.path.isdir(wfdir):
578 wfpointer= bppu.plot_waveform(pos=pos,siminspiral=injfile,event=eventnum,path=wfdir)
579 except Exception
as e:
581 print(
"Could not create WF plot. The error was: %s\n"%
str(e))
582 wftd=html_wf.insert_td(row,
'',label=
'Waveform',legend=legend)
583 wfsection=html.add_section_to_element(
'Waveforms',wftd)
584 if wfpointer
is not None:
585 wfsection.write(
'<a href="Waveform/WF_DetFrame.png" target="_blank"><img src="Waveform/WF_DetFrame.png"/></a>')
587 print(
"Could not create WF plot.\n")
588 wfsection.write(
"<b>No Waveform generated!</b>")
590 wftd=html_wf.insert_td(row,
'',label=
'PSDs',legend=legend)
591 wfsection=html.add_section_to_element(
'PSDs',wftd)
592 if psd_files
is not None:
593 psd_files=list(psd_files.split(
','))
594 psddir=os.path.join(outdir,
'PSDs')
595 if not os.path.isdir(psddir):
598 if 'flow' in pos.names:
599 f_low = pos[
'flow'].samples.min()
602 bppu.plot_psd(psd_files,outpath=psddir, f_min=f_low)
603 wfsection.write(
'<a href="PSDs/PSD.png" target="_blank"><img src="PSDs/PSD.png"/></a>')
604 except Exception
as e:
605 print(
"Could not create PSD plot. The error was: %s\n"%
str(e))
606 wfsection.write(
"<b>PSD plotting failed</b>")
608 wfsection.write(
"<b>No PSD files provided</b>")
611 if np.any([
'spcal_amp' in param
for param
in pos.names])
or np.any([
'spcal_phase' in param
for param
in pos.names]):
612 wftd=html_wf.insert_td(row,
'',label=
'Calibration',legend=legend)
613 wfsection=html.add_section_to_element(
'Calibration',wftd)
614 bppu.plot_calibration_pos(pos, outpath=outdir)
615 wfsection.write(
'<a href="calibration.png" target="_blank"><img src="calibration.png"/></a>')
618 if set([
'a1',
'a2',
'tilt1',
'tilt2']).issubset(pos.names):
619 wftd=html_wf.insert_td(row,
'',label=
'DiskPlot',legend=legend)
620 wfsection=html.add_section_to_element(
'DiskPlot',wftd)
621 lalinference.plot.make_disk_plot(pos, outpath=outdir)
622 wfsection.write(
'<a href="comp_spin_pos.png" target="_blank"><img src="comp_spin_pos.png"/></a>')
627 onepdfdir=os.path.join(outdir,
'1Dpdf')
628 if not os.path.isdir(onepdfdir):
629 os.makedirs(onepdfdir)
631 sampsdir=os.path.join(outdir,
'1Dsamps')
632 if not os.path.isdir(sampsdir):
633 os.makedirs(sampsdir)
636 for i
in oneDMenus.keys():
637 rss=bppu.make_1d_table(html,legend,i,pos,oneDMenus[i],noacf,GreedyRes,onepdfdir,sampsdir,savepdfs,greedy,analyticLikelihood,nDownsample)
643 tabid=
'onedconftable'
644 html_ogci=html.add_collapse_section(
'1D confidence intervals (greedy binning)',legend=legend,innertable_id=tabid)
645 html_ogci_write=
'<table id="%s" border="1"><tr><th/>'%tabid
646 clasciiout=
"#parameter \t"
647 confidence_levels.sort()
648 for cl
in confidence_levels:
649 html_ogci_write+=
'<th>%f</th>'%cl
650 clasciiout+=
"%s\t"%(
'%.02f'%cl)
652 html_ogci_write+=
'<th>Injection Confidence Level</th>'
653 html_ogci_write+=
'<th>Injection Confidence Interval</th>'
654 clasciiout+=
"Injection_Confidence_Level\t"
655 clasciiout+=
"Injection_Confidence_Interval"
657 html_ogci_write+=
'</tr>'
660 for par_name
in oneDMenu:
661 par_name=par_name.lower()
663 pos[par_name.lower()]
668 par_bin=GreedyRes[par_name]
670 print(
"Bin size is not set for %s, skipping binning."%par_name)
672 binParams={par_name:par_bin}
676 print(
"Using greedy 1-d binning credible regions\n")
678 toppoints,injectionconfidence,reses,injection_area,cl_intervals=bppu.greedy_bin_one_param(pos,binParams,confidence_levels)
681 print(
"Using 2-step KDE 1-d credible regions\n")
683 if pos[par_name].injval
is None:
686 injCoords=[pos[par_name].injval]
687 _,reses,injstats=bppu.kdtree_bin2Step(pos,[par_name],confidence_levels,injCoords=injCoords)
688 if injstats
is not None:
689 injectionconfidence=injstats[3]
690 injection_area=injstats[4]
692 BCItableline=
'<tr><td>%s</td>'%(par_name)
693 clasciiout+=
"%s\t"%par_name
694 cls=list(reses.keys())
698 BCItableline+=
'<td>%f</td>'%reses[cl]
699 clasciiout+=
"%f\t"%reses[cl]
700 if injection
is not None:
701 if injectionconfidence
is not None and injection_area
is not None:
703 BCItableline+=
'<td>%f</td>'%injectionconfidence
704 BCItableline+=
'<td>%f</td>'%injection_area
705 clasciiout+=
"%f\t"%injectionconfidence
706 clasciiout+=
"%f"%injection_area
709 BCItableline+=
'<td/>'
710 BCItableline+=
'<td/>'
713 BCItableline+=
'</tr>'
716 html_ogci_write+=BCItableline
718 html_ogci_write+=
'</table>'
719 html_ogci.write(html_ogci_write)
724 cornerdir=os.path.join(outdir,
'corner')
725 if not os.path.isdir(cornerdir):
726 os.makedirs(cornerdir)
727 massParams=[
'mtotal',
'm1',
'm2',
'mc']
728 distParams=[
'distance',
'distMPC',
'dist',
'distance_maxl']
729 incParams=[
'iota',
'inclination',
'theta_jn']
730 polParams=[
'psi',
'polarisation',
'polarization']
731 skyParams=[
'ra',
'rightascension',
'declination',
'dec']
733 spininducedquadParams = [
'dquadmon1',
'dquadmon2',
'dquadmons',
'dquadmona']
734 spinParams=[
'spin1',
'spin2',
'a1',
'a2',
'a1z',
'a2z',
'phi1',
'theta1',
'phi2',
'theta2',
'chi',
'effectivespin',
'chi_eff',
'chi_tot',
'chi_p',
'beta',
'tilt1',
'tilt2',
'phi_jl',
'theta_jn',
'phi12']
735 sourceParams=[
'm1_source',
'm2_source',
'mtotal_source',
'mc_source',
'redshift']
736 intrinsicParams=massParams+spinParams
737 extrinsicParams=incParams+distParams+polParams+skyParams
738 sourceFrameParams=sourceParams+distParams
740 myfig=bppu.plot_corner(pos,[0.05,0.5,0.95],parnames=intrinsicParams)
747 html_corner+=
'<tr><td width="100%"><a href="corner/intrinsic.png" target="_blank"><img width="70%" src="corner/intrinsic.png"/></a></td></tr>'
748 myfig.savefig(os.path.join(cornerdir,
'intrinsic.png'))
749 myfig.savefig(os.path.join(cornerdir,
'intrinsic.pdf'))
752 myfig=bppu.plot_corner(pos,[0.05,0.5,0.95],parnames=extrinsicParams)
757 myfig.savefig(os.path.join(cornerdir,
'extrinsic.png'))
758 myfig.savefig(os.path.join(cornerdir,
'extrinsic.pdf'))
759 html_corner+=
'<tr><td width="100%"><a href="corner/extrinsic.png" target="_blank"><img width="70%" src="corner/extrinsic.png"/></a></td></tr>'
762 myfig=bppu.plot_corner(pos,[0.05,0.5,0.95],parnames=sourceFrameParams)
767 myfig.savefig(os.path.join(cornerdir,
'sourceFrame.png'))
768 myfig.savefig(os.path.join(cornerdir,
'sourceFrame.pdf'))
769 html_corner+=
'<tr><td width="100%"><a href="corner/sourceFrame.png" target="_blank"><img width="70%" src="corner/sourceFrame.png"/></a></td></tr>'
773 html_corner=
'<table id="%s" border="1">'%tabid+html_corner
774 html_corner+=
'</table>'
776 html_co=html.add_collapse_section(
'Corner plots',legend=legend,innertable_id=tabid)
777 html_co.write(html_corner)
779 fout=open(os.path.join(outdir,
'confidence_levels.txt'),
'w')
780 fout.write(clasciiout)
791 margdir=os.path.join(outdir,
'2Dkde')
792 if not os.path.isdir(margdir):
795 twobinsdir=os.path.join(outdir,
'2Dbins')
796 if not os.path.isdir(twobinsdir):
797 os.makedirs(twobinsdir)
799 greedytwobinsdir=os.path.join(outdir,
'greedy2Dbins')
800 if not os.path.isdir(greedytwobinsdir):
801 os.makedirs(greedytwobinsdir)
806 html_tcig=html.add_collapse_section(
'2D confidence intervals (greedy binning)',legend=legend,innertable_id=tabid)
808 html_tcig_write=
'<table id="%s" border="1"><tr><th/>'%tabid
809 confidence_levels.sort()
810 for cl
in confidence_levels:
811 html_tcig_write+=
'<th>%f</th>'%cl
813 html_tcig_write+=
'<th>Injection Confidence Level</th>'
814 html_tcig_write+=
'<th>Injection Confidence Interval</th>'
815 html_tcig_write+=
'</tr>'
819 twodkdeplots_flag=twodkdeplots
820 if twodkdeplots_flag:
822 html_tcmp=html.add_collapse_section(
'2D Marginal PDFs',legend=legend,innertable_id=tabid)
824 html_tcmp_write=
'<table border="1" id="%s">'%tabid
826 tabid=
'2dgreedytable'
827 html_tgbh=html.add_collapse_section(
'2D Greedy Bin Histograms',legend=legend,innertable_id=tabid)
828 html_tgbh_write=
'<table border="1" id="%s">'%tabid
833 for par1_name,par2_name
in twoDGreedyMenu:
834 par1_name=par1_name.lower()
835 par2_name=par2_name.lower()
837 if par1_name == par2_name:
continue
839 pos[par1_name.lower()]
844 pos[par2_name.lower()]
850 par1_bin=GreedyRes[par1_name]
852 print(
"Bin size is not set for %s, skipping %s/%s binning."%(par1_name,par1_name,par2_name))
855 par2_bin=GreedyRes[par2_name]
857 print(
"Bin size is not set for %s, skipping %s/%s binning."%(par2_name,par1_name,par2_name))
861 par1_pos=pos[par1_name].samples
862 par2_pos=pos[par2_name].samples
863 if (
size(unique(par1_pos))<2
or size(unique(par2_pos))<2):
868 greedy2Params={par1_name:par1_bin,par2_name:par2_bin}
871 toppoints,injection_cl,reses,injection_area=\
872 bppu.greedy_bin_two_param(pos,greedy2Params,confidence_levels)
874 print(
"BCI %s-%s:"%(par1_name,par2_name))
878 BCItableline=
'<tr><td>%s-%s</td>'%(par1_name,par2_name)
879 cls=list(reses.keys())
883 BCItableline+=
'<td>%f</td>'%reses[cl]
885 if injection
is not None:
886 if injection_cl
is not None:
887 BCItableline+=
'<td>%f</td>'%injection_cl
888 BCItableline+=
'<td>'+
str(injection_area)+
'</td>'
891 BCItableline+=
'<td/>'
892 BCItableline+=
'<td/>'
894 BCItableline+=
'</tr>'
897 html_tcig_write+=BCItableline
903 greedy2ContourPlot=bppu.plot_two_param_kde_greedy_levels({
'Result':pos},greedy2Params,[0.67,0.9,0.95],{
'Result':
'k'})
904 greedy2contourpath=os.path.join(greedytwobinsdir,
'%s-%s_greedy2contour.png'%(par1_name,par2_name))
905 if greedy2ContourPlot
is not None:
906 greedy2ContourPlot.savefig(greedy2contourpath)
907 if(savepdfs): greedy2ContourPlot.savefig(greedy2contourpath.replace(
'.png',
'.pdf'))
908 plt.close(greedy2ContourPlot)
910 greedy2HistFig=bppu.plot_two_param_greedy_bins_hist(pos,greedy2Params,confidence_levels)
911 greedy2histpath=os.path.join(greedytwobinsdir,
'%s-%s_greedy2.png'%(par1_name,par2_name))
912 greedy2HistFig.savefig(greedy2histpath)
913 if(savepdfs): greedy2HistFig.savefig(greedy2histpath.replace(
'.png',
'.pdf'))
914 plt.close(greedy2HistFig)
916 greedyFile = open(os.path.join(twobinsdir,
'%s_%s_greedy_stats.txt'%(par1_name,par2_name)),
'w')
920 greedyFile.write(
"%lf %lf\n"%(cl,reses[cl]))
923 if [par1_name,par2_name]
in twoDplots
or [par2_name,par1_name]
in twoDplots :
924 print(
'Generating %s-%s greedy hist plot'%(par1_name,par2_name))
926 par1_pos=pos[par1_name].samples
927 par2_pos=pos[par2_name].samples
929 if (
size(unique(par1_pos))<2
or size(unique(par2_pos))<2):
931 head,figname=os.path.split(greedy2histpath)
932 head,figname_c=os.path.split(greedy2contourpath)
934 html_tgbh_write+=
'<tr>'
935 html_tgbh_write+=
'<td width="30%"><img width="100%" src="greedy2Dbins/'+figname+
'"/>[<a href="greedy2Dbins/'+figname_c+
'">contour</a>]</td>'
938 html_tgbh_write+=
'</tr>'
943 if twodkdeplots_flag
is True:
944 if [par1_name,par2_name]
in twoDplots
or [par2_name,par1_name]
in twoDplots :
945 print(
'Generating %s-%s plot'%(par1_name,par2_name))
947 par1_pos=pos[par1_name].samples
948 par2_pos=pos[par2_name].samples
950 if (
size(unique(par1_pos))<2
or size(unique(par2_pos))<2):
953 plot2DkdeParams={par1_name:50,par2_name:50}
954 myfig=bppu.plot_two_param_kde(pos,plot2DkdeParams)
956 figname=par1_name+
'-'+par2_name+
'_2Dkernel.png'
957 twoDKdePath=os.path.join(margdir,figname)
960 html_tcmp_write+=
'<tr>'
961 html_tcmp_write+=
'<td width="30%"><img width="100%" src="2Dkde/'+figname+
'"/></td>'
964 html_tcmp_write+=
'</tr>'
968 myfig.savefig(twoDKdePath)
969 if(savepdfs): myfig.savefig(twoDKdePath.replace(
'.png',
'.pdf'))
972 print(
'Unable to generate 2D kde levels for %s-%s'%(par1_name,par2_name))
976 html_tcig_write+=
'</table>'
977 html_tcig.write(html_tcig_write)
979 if twodkdeplots_flag
is True:
982 html_tcmp_write+=
'<td/>'
986 html_tcmp_write+=
'</tr>'
987 html_tcmp_write+=
'</table>'
988 html_tcmp.write(html_tcmp_write)
990 html_tcmp.a(
"2Dkde/",
'All 2D marginal PDFs (kde)')
993 while row_count_gb!=0:
994 html_tgbh_write+=
'<td/>'
998 html_tgbh_write+=
'</tr>'
999 html_tgbh_write+=
'</table>'
1000 html_tgbh.write(html_tgbh_write)
1002 html_tgbh.a(
"greedy2Dbins/",
'All 2D Greedy Bin Histograms')
1004 if RconvergenceTests
is True:
1005 convergenceResults=bppu.convergenceTests(pos,gelman=
False)
1007 if convergenceResults
is not None:
1009 html_conv_test=html.add_collapse_section(
'Convergence tests',legend=legend,innertable_id=tabid)
1011 for test,test_data
in convergenceResults.items():
1015 html_conv_test.h3(test)
1017 html_conv_table_rows={}
1018 html_conv_table_header=
''
1019 for chain,chain_data
in test_data.items():
1020 html_conv_table_header+=
'<th>%s</th>'%chain
1023 for data
in chain_data:
1026 html_conv_table_rows[data[0]]+=
'<td>'+data[1]+
'</td>'
1028 html_conv_table_rows[data[0]]=
'<td>'+data[1]+
'</td>'
1030 html_conv_table=
'<table id="%s"><tr><th>Chain</th>'%tabid+html_conv_table_header+
'</tr>'
1031 for row_name,row
in html_conv_table_rows.items():
1032 html_conv_table+=
'<tr><td>%s</td>%s</tr>'%(row_name,row)
1033 html_conv_table+=
'</table>'
1034 html_conv_test.write(html_conv_table)
1035 if data_found
is False:
1036 html_conv_test.p(
'No convergence diagnostics generated!')
1040 html_stats_cov=html.add_collapse_section(
'Covariance matrix',legend=legend,innertable_id=tabid)
1041 pos_samples,table_header_string=pos.samples()
1045 cov_matrix=cov(pos_samples,rowvar=0,bias=1)
1048 table_header_list=table_header_string.split()
1049 cov_table_string=
'<table border="1" id="%s"><tr><th/>'%tabid
1050 for header
in table_header_list:
1051 cov_table_string+=
'<th>%s</th>'%header
1052 cov_table_string+=
'</tr>'
1054 cov_column_list=hsplit(cov_matrix,cov_matrix.shape[1])
1056 for cov_column,cov_column_name
in zip(cov_column_list,table_header_list):
1057 cov_table_string+=
'<tr><th>%s</th>'%cov_column_name
1058 for cov_column_element
in cov_column:
1059 cov_table_string+=
'<td>%.3e</td>'%(cov_column_element[0])
1060 cov_table_string+=
'</tr>'
1061 cov_table_string+=
'</table>'
1062 html_stats_cov.write(cov_table_string)
1065 print(
'Unable to compute the covariance matrix.')
1068 html_footer=html.add_section(
'')
1069 html_footer.p(
'Produced using cbcBayesPostProc.py at '+strftime(
"%Y-%m-%d %H:%M:%S")+
' .')
1072 for arg
in sys.argv:
1075 html_footer.p(
'Command line: %s'%cc_args)
1076 html_footer.p(git_version.verbose_msg)
1079 resultspage=open(os.path.join(outdir,
'posplots.html'),
'w')
1080 resultspage.write(
str(html))
1086 USAGE=
'''%prog [options] datafile.dat [datafile2.dat ...]
1087 Generate a web page displaying results of parameter estimation based on the contents
1088 of one or more datafiles containing samples from one of the bayesian algorithms (MCMC, nested sampling).
1089 Options specify which extra statistics to compute and allow specification of additional information.
1092 if __name__==
'__main__':
1094 from optparse
import OptionParser
1095 parser=OptionParser(USAGE)
1096 parser.add_option(
"-o",
"--outpath", dest=
"outpath",help=
"make page and plots in DIR", metavar=
"DIR")
1097 parser.add_option(
"-d",
"--data",dest=
"data",action=
"callback",callback=multipleFileCB,help=
"datafile")
1099 parser.add_option(
"-i",
"--inj",dest=
"injfile",help=
"SimInsipral injection file",metavar=
"INJ.XML",default=
None)
1100 parser.add_option(
"-t",
"--trig",dest=
"trigfile",help=
"Coinc XML file",metavar=
"COINC.XML",default=
None)
1101 parser.add_option(
"--skyres",dest=
"skyres",help=
"Sky resolution to use to calculate sky box size",default=
None)
1102 parser.add_option(
"--eventnum",dest=
"eventnum",action=
"store",default=
None,help=
"event number in SimInspiral file of this signal",type=
"int",metavar=
"NUM")
1103 parser.add_option(
"--trignum",dest=
"trignum",action=
"store",default=
None,help=
"trigger number in CoincTable",type=
"int",metavar=
"NUM")
1104 parser.add_option(
"--bsn",action=
"store",default=
None,help=
"Optional file containing the bayes factor signal against noise",type=
"string")
1105 parser.add_option(
"--bci",action=
"store",default=
None,help=
"Optional file containing the bayes factor coherent signal model against incoherent signal model.",type=
"string")
1106 parser.add_option(
"--snr",action=
"store",default=
None,help=
"Optional file containing the SNRs of the signal in each IFO",type=
"string")
1107 parser.add_option(
"--dievidence",action=
"store_true",default=
False,help=
"Calculate the direct integration evidence for the posterior samples")
1108 parser.add_option(
"--boxing",action=
"store",default=64,help=
"Boxing parameter for the direct integration evidence calculation",type=
"int",dest=
"boxing")
1109 parser.add_option(
"--evidenceFactor",action=
"store",default=1.0,help=
"Overall factor (normalization) to apply to evidence",type=
"float",dest=
"difactor",metavar=
"FACTOR")
1111 parser.add_option(
'--ellipticEvidence', action=
'store_true', default=
False,help=
'Estimate the evidence by fitting ellipse to highest-posterior points.', dest=
'ellevidence')
1113 parser.add_option(
"--plot-2d", action=
"store_true", default=
False,help=
"Make individual 2-D plots.")
1114 parser.add_option(
"--header",action=
"store",default=
None,help=
"Optional file containing the header line for posterior samples",type=
"string")
1116 parser.add_option(
"--ns",action=
"store_true",default=
False,help=
"(inspnest) Parse input as if it was output from parallel nested sampling runs.")
1117 parser.add_option(
"--Nlive",action=
"store",default=
None,help=
"(inspnest) Number of live points used in each parallel nested sampling run.",type=
"int")
1118 parser.add_option(
"--xflag",action=
"store_true",default=
False,help=
"(inspnest) Convert x to iota.")
1120 parser.add_option(
"--ss",action=
"store_true",default=
False,help=
"(SPINspiral) Parse input as if it was output from SPINspiral.")
1121 parser.add_option(
"--spin",action=
"store_true",default=
False,help=
"(SPINspiral) Specify spin run (15 parameters). ")
1123 parser.add_option(
"--lalinfmcmc",action=
"store_true",default=
False,help=
"(LALInferenceMCMC) Parse input from LALInferenceMCMC.")
1124 parser.add_option(
"--inj-spin-frame",default=
'OrbitalL', help=
"The reference frame used for the injection (default: OrbitalL)")
1125 parser.add_option(
"--downsample",action=
"store",default=
None,help=
"(LALInferenceMCMC) approximate number of samples to record in the posterior",type=
"int")
1126 parser.add_option(
"--deltaLogL",action=
"store",default=
None,help=
"(LALInferenceMCMC) Difference in logL to use for convergence test. (DEPRECATED)",type=
"float")
1127 parser.add_option(
"--deltaLogP",action=
"store",default=
None,help=
"(LALInferenceMCMC) Difference in logpost to use for burnin criteria.",type=
"float")
1128 parser.add_option(
"--fixedBurnin",dest=
"fixedBurnin",action=
"callback",callback=multipleFileCB,help=
"(LALInferenceMCMC) Fixed number of iteration for burnin.")
1129 parser.add_option(
"--oldMassConvention",action=
"store_true",default=
False,help=
"(LALInferenceMCMC) if activated, m2 > m1; otherwise m1 > m2 in PTMCMC.output.*.00")
1131 parser.add_option(
"--fm",action=
"store_true",default=
False,help=
"(followupMCMC) Parse input as if it was output from followupMCMC.")
1133 parser.add_option(
"--no-acf", action=
"store_true", default=
False, dest=
"noacf")
1135 parser.add_option(
"--twodkdeplots", action=
"store_true", default=
False, dest=
"twodkdeplots")
1137 parser.add_option(
"--RconvergenceTests", action=
"store_true", default=
False, dest=
"RconvergenceTests")
1138 parser.add_option(
"--nopdfs",action=
"store_false",default=
True,dest=
"nopdfs")
1139 parser.add_option(
"-c",
"--covarianceMatrix",dest=
"covarianceMatrices",action=
"append",default=
None,help=
"CSV file containing covariance (must give accompanying mean vector CSV. Can add more than one matrix.")
1140 parser.add_option(
"-m",
"--meanVectors",dest=
"meanVectors",action=
"append",default=
None,help=
"Comma separated list of locations of the multivariate gaussian described by the correlation matrix. First line must be list of params in the order used for the covariance matrix. Provide one list per covariance matrix.")
1141 parser.add_option(
"--email",action=
"store",default=
None,type=
"string",metavar=
"user@ligo.org",help=
"Send an e-mail to the given address with a link to the finished page.")
1142 parser.add_option(
"--archive",action=
"store",default=
None,type=
"string",metavar=
"results.tar.gz",help=
"Create the given tarball with all results")
1143 parser.add_option(
"--psdfiles",action=
"store",default=
None,type=
"string",metavar=
"H1,L1,V1",help=
"comma separater list of ASCII files with PSDs, one per IFO")
1144 parser.add_option(
"--kdecredibleregions",action=
"store_true",default=
False,help=
"If given, will use 2-step KDE trees to estimate 1-d credible regions [default false: use greedy binning]")
1145 parser.add_option(
"--noplot-source-frame", action=
"store_true", default=
False,help=
"Don't make 1D plots of source-frame masses")
1146 (opts,args)=parser.parse_args()
1150 datafiles=datafiles+args
1152 datafiles=datafiles + opts.data
1154 if opts.fixedBurnin:
1156 if len(opts.fixedBurnin) == 1:
1157 fixedBurnins = [
int(opts.fixedBurnin[0])
for df
in datafiles]
1159 fixedBurnins = [
int(fixedBurnin)
for fixedBurnin
in opts.fixedBurnin]
1162 from lalinference.bayespputils import massParams,spinParams,cosmoParam,strongFieldParams,distParams,incParams,polParams,skyParams,phaseParams,timeParams,endTimeParams,statsParams,calibParams,snrParams,tidalParams,fourPiecePolyParams,spectralParams
1165 oneDMenus={
'Masses':
None,
'SourceFrame':
None,
'Timing':
None,
'Extrinsic':
None,
'Spins':
None,
'StrongField':
None,
'Others':
None}
1167 oneDMenus[
'Masses']= massParams
1168 oneDMenus[
'Extrinsic']=incParams+distParams+polParams+skyParams+phaseParams
1169 oneDMenus[
'Spins']= spinParams
1170 oneDMenus[
'Timing']=timeParams+endTimeParams
1171 oneDMenus[
'StrongField']= strongFieldParams
1172 oneDMenus[
'Others']=snrParams+statsParams+calibParams
1173 oneDMenus[
'SourceFrame']= cosmoParam
1175 if opts.noplot_source_frame:
1176 oneDMenus[
'SourceFrame']= []
1178 ifos_menu=[
'h1',
'l1',
'v1']
1179 from itertools
import combinations
1180 for ifo1,ifo2
in combinations(ifos_menu,2):
1181 oneDMenus[
'Timing'].append(ifo1+ifo2+
'_delay')
1185 for mp1,mp2
in combinations(massParams,2):
1186 twoDGreedyMenu.append([mp1, mp2])
1187 for mp
in massParams:
1188 for d
in distParams:
1189 twoDGreedyMenu.append([mp,d])
1190 for mp
in massParams:
1191 for sp
in spinParams:
1192 twoDGreedyMenu.append([mp,sp])
1193 for mp
in massParams:
1194 for dchi
in bppu.tigerParams:
1195 twoDGreedyMenu.append([mp,dchi])
1196 for mp
in massParams:
1197 for lvp
in bppu.lorentzInvarianceViolationParams:
1198 twoDGreedyMenu.append([mp,lvp])
1199 for sp
in spinParams:
1200 for lvp
in bppu.lorentzInvarianceViolationParams:
1201 twoDGreedyMenu.append([sp,lvp])
1202 for dchi
in bppu.tigerParams:
1203 for lvp
in bppu.lorentzInvarianceViolationParams:
1204 twoDGreedyMenu.append([dchi, lvp])
1205 for eg
in bppu.energyParams:
1206 for lvp
in bppu.lorentzInvarianceViolationParams:
1207 twoDGreedyMenu.append([eg, lvp])
1208 for dp
in distParams:
1209 for lvp
in bppu.lorentzInvarianceViolationParams:
1210 twoDGreedyMenu.append([dp, lvp])
1211 for dp
in distParams:
1212 for sp
in snrParams:
1213 twoDGreedyMenu.append([dp,sp])
1214 for dp
in distParams:
1215 for ip
in incParams:
1216 twoDGreedyMenu.append([dp,ip])
1217 for dp
in distParams:
1218 for sp
in skyParams:
1219 twoDGreedyMenu.append([dp,sp])
1220 for dp
in distParams:
1221 for sp
in spinParams:
1222 twoDGreedyMenu.append([dp,sp])
1223 for ip
in incParams:
1224 for sp
in skyParams:
1225 twoDGreedyMenu.append([ip,sp])
1226 for ip
in incParams:
1227 for sp
in spinParams:
1228 twoDGreedyMenu.append([ip,sp])
1229 for phip
in phaseParams:
1230 twoDGreedyMenu.append([ip,phip])
1231 for psip
in polParams:
1232 twoDGreedyMenu.append([ip,psip])
1233 for sp1
in skyParams:
1234 for sp2
in skyParams:
1235 if not (sp1 == sp2):
1236 twoDGreedyMenu.append([sp1, sp2])
1237 for sp1,sp2
in combinations(spinParams,2):
1238 twoDGreedyMenu.append([sp1, sp2])
1239 for dc1,dc2
in combinations(bppu.tigerParams,2):
1240 twoDGreedyMenu.append([dc1,dc2])
1241 for mp
in massParams:
1242 for tp
in tidalParams:
1244 twoDGreedyMenu.append([mp, tp])
1245 for tp
in fourPiecePolyParams:
1247 twoDGreedyMenu.append([mp, tp])
1248 for tp
in spectralParams:
1250 twoDGreedyMenu.append([mp, tp])
1251 for sp1,sp2
in combinations(snrParams,2):
1252 twoDGreedyMenu.append([sp1,sp2])
1253 twoDGreedyMenu.append([
'lambda1',
'lambda2'])
1254 twoDGreedyMenu.append([
'lam_tilde',
'dlam_tilde'])
1255 twoDGreedyMenu.append([
'lambdat',
'dlambdat'])
1256 twoDGreedyMenu.append([
'logp1',
'gamma1',
'gamma2',
'gamma3'])
1257 twoDGreedyMenu.append([
'SDgamma0',
'SDgamma1',
'SDgamma2',
'SDgamma3'])
1258 for psip
in polParams:
1259 for phip
in phaseParams:
1260 twoDGreedyMenu.append([psip,phip])
1261 for sp
in skyParams:
1262 twoDGreedyMenu.append([psip,sp])
1263 for sp
in spinParams:
1264 twoDGreedyMenu.append([psip,sp])
1266 for i
in calibParams[3:]:
1267 twoDGreedyMenu.append([i,
'distance'])
1271 greedyBinSizes=bppu.greedyBinSizes
1275 for dt1,dt2
in combinations([
'h1_end_time',
'l1_end_time',
'v1_end_time'],2):
1276 twoDGreedyMenu.append([dt1,dt2])
1277 for dt1,dt2
in combinations( [
'h1l1_delay',
'l1v1_delay',
'h1v1_delay'],2):
1278 twoDGreedyMenu.append([dt1,dt2])
1280 if opts.deltaLogL
and not opts.deltaLogP:
1281 print(
"DEPRECATION WARNING: --deltaLogL has been replaced by --deltaLogP. Using the posterior to define burnin criteria")
1282 deltaLogP = opts.deltaLogL
1284 deltaLogP = opts.deltaLogP
1286 confidenceLevels=bppu.confidenceLevels
1289 twoDplots=twoDGreedyMenu
1291 opts.outpath,datafiles,oneDMenus,twoDGreedyMenu,
1292 greedyBinSizes,confidenceLevels,twoDplots,
1294 injfile=opts.injfile,eventnum=opts.eventnum,
1295 trigfile=opts.trigfile,trignum=opts.trignum,
1298 dievidence=opts.dievidence,boxing=opts.boxing,difactor=opts.difactor,
1300 ellevidence=opts.ellevidence,
1302 bayesfactornoise=opts.bsn,bayesfactorcoherent=opts.bci,
1306 ns_flag=opts.ns,ns_Nlive=opts.Nlive,
1308 ss_flag=opts.ss,ss_spin_flag=opts.spin,
1310 li_flag=opts.lalinfmcmc,deltaLogP=deltaLogP,fixedBurnins=fixedBurnins,nDownsample=opts.downsample,oldMassConvention=opts.oldMassConvention,
1314 inj_spin_frame=opts.inj_spin_frame,
1318 twodkdeplots=opts.twodkdeplots,
1320 RconvergenceTests=opts.RconvergenceTests,
1322 savepdfs=opts.nopdfs,
1324 covarianceMatrices=opts.covarianceMatrices,
1326 meanVectors=opts.meanVectors,
1330 psd_files=opts.psdfiles,
1331 greedy=not(opts.kdecredibleregions)
1334 if opts.archive
is not None:
1336 subprocess.call([
"tar",
"cvzf",opts.archive,opts.outpath])
1344 except Exception
as e:
1345 print(
'Unable to send notification email')
1346 print(
"The error was %s\n"%
str(e))
def oneD_dict_to_file(dict, fname)
def cbcBayesPostProc(outdir, data, oneDMenus, twoDGreedyMenu, GreedyRes, confidence_levels, twoDplots, injfile=None, eventnum=None, trigfile=None, trignum=None, skyres=None, dievidence=False, boxing=64, difactor=1.0, ellevidence=False, bayesfactornoise=None, bayesfactorcoherent=None, snrfactor=None, ns_flag=False, ns_Nlive=None, ss_flag=False, ss_spin_flag=False, li_flag=False, deltaLogP=None, fixedBurnins=None, nDownsample=None, oldMassConvention=False, fm_flag=False, inj_spin_frame='OrbitalL', noacf=False, twodkdeplots=False, RconvergenceTests=False, savepdfs=True, covarianceMatrices=None, meanVectors=None, header=None, psd_files=None, greedy=True ## If true will use greedy bin for 1-d credible regions. Otherwise use 2-steps KDE)
This is a demonstration script for using the functionality/data structures contained in lalinference....
def pickle_to_file(obj, fname)
Pickle/serialize 'obj' into 'fname'.
def email_notify(address, path)
def dict2html(d, parent=None)
def extract_hdf5_metadata(h5grp, parent=None)
def guess_url(fslocation)
Try to work out the web address of a given path.