Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALPulsar 7.1.1.1-5e288d3
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
lalpulsar_knope_result_page.py
Go to the documentation of this file.
1##python
2# -*- coding: utf-8 -*-
3#
4# lalpulsar_knope_result_page.py
5#
6# Copyright 2015
7# Matthew Pitkin <matthew.pitkin@ligo.org>,
8#
9# This program is free software; you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation; either version 2 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program; if not, write to the Free Software
21# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
22# MA 02110-1301, USA.
23
24## \file
25## \ingroup lalpulsar_bin_HeterodyneSearch
26"""
27The KNOwn Pulsar pipelinE - lalpulsar_knope_result_page
28
29Script from creating results pages for pulsars from the known pulsar search
30"""
31
32from __future__ import print_function, division
33
34import argparse
35from configparser import ConfigParser
36import sys
37import ast
38import numpy as np
39import re
40import copy
41import os
42import fnmatch
43import ast
44import datetime
45import json
46from scipy import stats
47import h5py
48import itertools
49
50import matplotlib
51
52matplotlib.use("Agg")
53
54from lalpulsar.pulsarpputils import *
55from lalpulsar.pulsarhtmlutils import *
56from lalinference.bayespputils import Posterior, PosteriorOneDPDF
57from lalinference import git_version
58
59__author__ = "Matthew Pitkin <matthew.pitkin@ligo.org>"
60__version__ = "git id %s" % git_version.id
61__date__ = git_version.date
62
63# try importing scotchcorner
64try:
65 from scotchcorner import scotchcorner
66except ImportError:
67 print(
68 "Could no import scotchcorner: make sure scotchcorner is installed (e.g. 'pip install scotchcorner') and in the PYTHONPATH",
69 file=sys.stderr,
70 )
71 sys.exit(1)
72
73# create format for the output page
74htmlpage = """
75<!DOCTYPE html>
76<html lang="en">
77
78<head>
79 <meta http-equiv="Content-Type" content="text/html; charset=utf-8"/>
80 <meta name="description" content="PSR {psrname}"/>
81 <meta charset="UTF-8">
82
83 <title>{title}</title>
84
85 <link rel="stylesheet" type="text/css" href="{cssfile}"/>
86
87</head>
88
89<body>
90
91<!- add in javascript to allow text toggling -->
92<script language="javascript">
93function toggle(id) {{
94 var ele = document.getElementById(id);
95 if(ele.style.display == "block"){{
96 ele.style.display = "none";
97 }}
98 else
99 ele.style.display = "block";
100}}
101</script>
102
103<!-- Page title -->
104<h1>{h1title}</h1>
105
106<!-- Links to the parts of the page -->
107<div class="pagelinks">
108{linkstable}
109</div>
110
111<!-- table of pulsar parameters from EM data -->
112<table>
113<tr>
114<td>
115<div class="pulsartable" id="pulsartable">
116{pulsartable}
117</div>
118</td>
119
120<td>
121<!-- table of derived gravitational wave upper limits, SNRs and evidence ratios -->
122<div class="limitstable" id="limitstable">
123{limitstable}
124</div>
125</td>
126</tr>
127</table>
128
129<br />
130
131<!-- plots of 1D and 2D marginalised posteriors for GW amplitude and phase parameters -->
132<div class="selectedposteriors" id="selectedposteriors">
133{selectedposteriors}
134</div>
135
136<br />
137
138<!-- corner plot of all parameters (hidden by default) -->
139<div class="jointposteriors" id="jointposteriors">
140{jointposteriors}
141</div>
142
143<br />
144
145<!-- (if background runs have been produced): if running a multi-detector analysis the plot
146 coherent vs. incoherent Bayes factors (with SNR giving colour of points) of background
147 and foreground, or if single detector plot signal vs. noise Bayes factor against SNR. -->
148<div class="evidenceplots" id="evidenceplots">
149{evidenceplots}
150</div>
151
152<br />
153
154<!-- plots of pre-processed data (heterodyned/spectrally interpolated time series) -->
155<div class="dataplots" id="dataplots">
156{dataplots}
157</div>
158
159<br />
160
161<!-- plots of the posterior samples for all parameters -->
162<div class="posteriorsamples" id="posteriorsamples">
163{posteriorsamples}
164</div>
165
166<br />
167
168<!-- table of posterior statistics: means, max. likelihoods, confidence intervals -->
169<div class="posteriorstats" id="posteriorstats">
170{posteriorstats}
171</div>
172
173<br />
174
175<!-- a footer -->
176<div id="footer">
177{footer}
178</div>
179
180</body>
181</html>
182"""
183
184
185def set_spin_down(p1_I, assoc, f0, f1, n=5.0):
186 """
187 Set the spin-down of the source based on the intrinsic period derivative (p1_I) corrected for any proper motion/
188 globular cluster acceleration if available, or if not give AND the pulsar is in a globular cluster base the
189 spin-down on assuming an age of 10^9 years (defaulting to the source being a gravitar, with n=5).
190 Otherwise just return the unadjusted spin-down.
191 """
192
193 if p1_I != None and p1_I > 0.0:
194 return -p1_I * f0**2 # convert period derivative into frequency derivative
195 elif assoc != None:
196 if "GC" in assoc: # check if a globular cluster pulsar
197 return -f0 / ((n - 1.0) * 1.0e9 * 365.25 * 86400.0)
198 else:
199 return f1
200 else:
201 return f1
202
203
204def create_psr_table(par):
205 """
206 Create a html table of some information from the pulsar parameter file.
207 """
208 table = htmltable() # create table class
209
210 # list of parameters to display (in this order)
211 paramdisplist = [
212 "RAJ",
213 "DECJ",
214 "F0",
215 "F1",
216 "F2",
217 "PEPOCH",
218 "A1",
219 "E",
220 "EPS1",
221 "EPS2",
222 "OM",
223 "T0",
224 "TASC",
225 "PB",
226 ]
227
228 for param in paramdisplist:
229 pa = par[param]
230 if pa != None:
231 table.addrow()
232 dispfunc = paramhtmldispfunc.__dict__[param]
233 table.adddata(paramhtmldict[param])
234 table.adddata(dispfunc(str(pa)))
235
236 return table.tabletext
237
238
239class create_data_table:
240 def __init__(
241 self, datafiles, outputdir, figformats=["png"], asdtime=86400, harmonics=[2.0]
242 ):
243 """
244 Initialise with a dictionary keyed in detector names containing paths to the equivalent pre-processed
245 data file(s) for that detector (if harmonics contains more than one frequency harmonic this there should
246 be a list of files, with one for each harmonic and given in the same order as the harmonics lists).
247 figformats gives a list of the output figure formats (defaulting to .png), and asd time is the time to use
248 for each FFT when creating the spectrograms/ASDs, which defaults to 86400 seconds (1 day).
249 """
250
251 self._datatable = htmltag(
252 "h2", "Analysis statistics", newline=True
253 ).text # the output html table
254
255 if not isinstance(harmonics, list):
256 self._harmonics = list(harmonics)
257 else:
258 self._harmonics = harmonics
259 self._nharmonics = len(self._harmonics)
260 self._ifos = list(datafiles.keys()) # get list of detectors
261 self._datafiles = (
262 {}
263 ) # dictionary (keyed to harmonics) of dictionaries (keyed to detector) of data files
264
265 for i, ifo in enumerate(datafiles):
266 if isinstance(datafiles[ifo], list):
267 if len(datafiles[ifo]) != self._nharmonics:
268 print(
269 "Error... number of pre-processed data files is not the same as the number of harmonics.",
270 file=sys.stderr,
271 )
272 sys.exit(1)
273 else:
274 if i == 0:
275 for h in self._harmonics:
276 self._datafiles[h] = {}
277
278 for j, df in enumerate(datafiles[ifo]):
279 if not os.path.isfile(df):
280 print(
281 "Error... pre-processed data file '%s' for '%s' does not exist."
282 % (df, ifo),
283 file=sys.stderr,
284 )
285 sys.exit(1)
286 else:
287 self._datafiles[self._harmonics[j]][ifo] = df
288 else:
289 if self._nharmonics == 1:
290 self._datafiles[self._harmonics[0]] = datafiles
291 for df in list(self._datafiles[self._harmonics[0]].values()):
292 if not os.path.isfile(df):
293 print(
294 "Error... pre-processed data file '%s' for '%s' does not exist."
295 % (df, ifo),
296 file=sys.stderr,
297 )
298 sys.exit(1)
299 break
300 else:
301 print(
302 "Error... number of pre-processed data files is not the same as the number of harmonics.",
303 file=sys.stderr,
304 )
305 sys.exit(1)
306
307 self._outputdir = outputdir
308 if not os.path.isdir(self._outputdir):
309 print(
310 "Error... output path '%s' for data plots does not exist"
311 % self._outputdir,
312 file=sys.stderr,
313 )
314 sys.exit(1)
315
316 self._figformats = figformats # figure output formats
317 if "png" not in figformats:
318 print(
319 "Error... must include 'png' in the output figure formats",
320 file=sys.stderr,
321 )
322 sys.exit(1)
323
324 # get plots of absolute values, amplitude spectral density and spectrograms of the pre-processed data
325 self._asdtime = asdtime
326 self._asds = {}
327 self.create_plots()
328
329 self._starts = {} # start times of data files
330 self._ends = {} # end times of data files
331 self._length = {} # length (seconds) of data files
332 self._duty_factor = {} # duty factors of data
333
334 # get analysis data statistics (start time, end time, duty factor)
335 self.get_stats()
336
337 # create table of results for all detectors
338 self.create_table()
339
340 def create_plots(self):
341 # create the plots of the absolute values, amplitude spectral density and spectrograms of the
342 # pre-processed data.
343 self._Bkplots = {}
344 self._asdplots = {}
345 self._fscanplots = {}
346
347 for h in self._harmonics:
348 self._Bkplots[h] = {}
349 self._asdplots[h] = {}
350 self._fscanplots[h] = {}
351 self._asds[h] = {}
352
353 Bkfigs, asdfigs, fscanfigs, asdlist, _ = plot_Bks_ASDs(
354 self._datafiles[h],
355 delt=self._asdtime,
356 plotpsds=True,
357 plotfscan=True,
358 removeoutlier=5.0,
359 )
360
361 if self._nharmonics == 1:
362 harmtext = "" # don't add any harmonics suffix to the file name
363 else:
364 harmtext = "_%d" % int(h)
365
366 # output the plots
367 for i, ifo in enumerate(self._ifos):
368 self._asds[h][ifo] = asdlist[
369 i
370 ] # set the estimate of the amplitude spectral density for that IFO
371
372 for ftype in self._figformats:
373 for fig, fignameprefix in zip(
374 [Bkfigs[i], fscanfigs[i], asdfigs[i]], ["Bk", "fscan", "ASD"]
375 ):
376 # try outputting figures
377 figname = os.path.join(
378 outdir, "%s_%s%s.%s" % (fignameprefix, ifo, harmtext, ftype)
379 )
380 try:
381 fig.savefig(figname)
382 except:
383 print(
384 "Error... problem creating figure '%s'." % figname,
385 file=sys.stderr,
386 )
387 sys.exit(1)
388
389 if ftype == "png": # add figure files
390 self._Bkplots[h][ifo] = os.path.join(
391 outdir, "Bk_%s%s.%s" % (ifo, harmtext, ftype)
392 )
393 self._asdplots[h][ifo] = os.path.join(
394 outdir, "ASD_%s%s.%s" % (ifo, harmtext, ftype)
395 )
396 self._fscanplots[h][ifo] = os.path.join(
397 outdir, "fscan_%s%s.%s" % (ifo, harmtext, ftype)
398 )
399
400 @property
401 def Bkplots(self):
402 # return dictionary of plots of absolute values of the pre-processed data for each detector
403 return self._Bkplots
404
405 @property
406 def asdplots(self):
407 # return dictionary of plot of the amplitude specatral densities for each detector
408 return self._asdplots
409
410 @property
411 def fscanplots(self):
412 # return dictionary of spectrogram plots of the pre-processed data for each detector
413 return self._fscanplots
414
415 @property
416 def asds(self):
417 # return dictionary of noise amplitude spectral density estimates for each detector
418 return self._asds
419
420 def get_stats(self):
421 # get analysis statistics (data start, end and duty cycle) for each detector
422 for i, ifo in enumerate(self._ifos):
423 # use np.loadtxt to open file and get start and end times
424 try:
425 bkd = np.loadtxt(
426 self._datafiles[self._harmonics[0]][ifo], comments=["#", "%"]
427 )
428 except:
429 print(
430 "Error... could not load file '%s'."
431 % self._datafiles[self._harmonics[0]][ifo],
432 file=sys.stderr,
433 )
434 sys.exit(1)
435
436 self._starts[ifo] = bkd[0, 0] # start time
437 self._ends[ifo] = bkd[-1, 0] # end time
438 dt = np.min(np.diff(bkd[:, 0])) # minimum time step in data
439 self._length[ifo] = float(len(bkd)) * dt
440 self._duty_factor[ifo] = (
441 100.0 * self._length[ifo] / (self._ends[ifo] - self._starts[ifo])
442 )
443
444 def analysis_stats_td(self, ifo):
445 # return the inside of the html <td> element giving the statistics for a particular detector
446 table = htmltable()
447 table.addrow()
448 table.adddata("&nbsp;")
449 table.adddata(ifo, dataclass=ifo)
450 tdvals = [
451 ("Start (GPS)", str(int(self._starts[ifo]))),
452 ("End (GPS)", str(int(self._ends[ifo]))),
453 ("Length (sec)", str(int(self._length[ifo]))),
454 ("Duty factor (%)", "%.lf" % self._duty_factor[ifo]),
455 ]
456 for item in tdvals:
457 table.addrow()
458 table.adddata(item[0])
459 table.adddata(item[1])
460
461 return table.tabletext
462
463 def create_table(self):
464 table = htmltable()
465
466 for ifo in self._ifos:
467 for h in self._harmonics:
468 table.addrow()
469 if self._nharmonics != 1:
470 table.adddata("Data for %df" % int(h), colspan=2, header=True)
471 table.addrow()
472
473 table.adddata(
474 self.analysis_stats_td(ifo), datastyle="text-align: center;"
475 )
476 plotlink = atag(
477 os.path.basename(self._Bkplots[h][ifo]),
478 '<img class="dataplot" src="{}"/>'.format(
479 os.path.basename(self._Bkplots[h][ifo])
480 ),
481 ).text
482 table.adddata(plotlink)
483 table.addrow()
484 plotlink = atag(
485 os.path.basename(self._asdplots[h][ifo]),
486 '<img class="asdplot" src="{}"/>'.format(
487 os.path.basename(self._asdplots[h][ifo])
488 ),
489 ).text
490 table.adddata(plotlink)
491 plotlink = atag(
492 os.path.basename(self._fscanplots[h][ifo]),
493 '<img class="dataplot" src="{}"/>'.format(
494 os.path.basename(self._fscanplots[h][ifo])
495 ),
496 ).text
497 table.adddata(plotlink)
498
499 self._datatable += table.tabletext
500
501 def __str__(self): # string method
502 return self._datatable
503
504
505class posteriors:
506 """
507 Get sample posteriors and created a set of functions for outputting tables, plots and posterior statistics
508 """
509
510 def __init__(
511 self,
512 postfiles,
513 outputdir,
514 ifos=None,
515 harmonics=[2],
516 modeltype="waveform",
517 biaxial=False,
518 usegwphase=False,
519 parfile=None,
520 priorfile=None,
521 subtracttruths=False,
522 showcontours=False,
523 ):
524 """
525 Initialise with a dictionary keyed in detector names containing paths to the equivalent posterior samples
526 file for that detector.
527 """
528 self._outputdir = outputdir
529 if not os.path.isdir(self._outputdir):
530 print(
531 "Error... output path '%s' for data plots does not exist"
532 % self._outputdir,
533 file=sys.stderr,
534 )
535 sys.exit(1)
536
537 if ifos is None: # get list of detectors from postfiles dictionary
538 self._ifos = list(postfiles.keys()) # get list of detectors
539 else:
540 if isinstance(ifos, list):
541 self._ifos = ifos
542 else:
543 self._ifos = [ifos]
544 # check ifos are in postfiles dictionary
545 for ifo in self._ifos:
546 if ifo not in postfiles:
547 print(
548 "Error... posterior files for detector '%s' not given" % ifo,
549 file=sys.stderr,
550 )
551 sys.exit(1)
552
553 self._postfiles = postfiles
554 self._posteriors = {} # dictionary of posterior objects
555 self._posterior_stats = {} # dictionary if posteriors statistics
556 self._signal_evidence = {} # dictionary of signal evidence values
557 self._noise_evidence = {} # dictionary of noise evidence values
558 self._maxL = {} # dictionary of maximum log likelihood values
559 self._Bsn = {} # dictionary of signal vs noise Bayes factors
560 self._signal_evidence = {} # dictionary of signal evidences
561 self._noise_evidence = {} # dictionary of noise evidences
562 self._Bci = None # coherent versus incoherent Bayes factor
563 self._Bcin = None # coherent versus incoherent or noise Bayes factor
564 self._optimal_snrs = {} # dictionary of optimal matched filter SNRs
565 self._parameters = [] # list of the source parameters in the posterior files
566 self._parfile = parfile # set the TEMPO(2) parameter file used
567 self._injection_parameters = None # injection/heterodyne parameters
568 self._priorfile = priorfile # the prior file used for the analysis
569 self._prior_parameters = {} # the prior parameters
570 self._injection_credible_regions = (
571 {}
572 ) # dictionary of minimal credible regions within which an injected parameter is found
573 self._harmonics = harmonics # list of frequency harmonics
574 self._modeltype = modeltype # the model type ('waveform' or 'source')
575 self._biaxial = (
576 biaxial # whether the source is a biaxial star (rather than triaxial)
577 )
578 self._usegwphase = (
579 usegwphase # whether to use GW phase rather than rotational phase
580 )
581 self._subtract_truths = subtracttruths # set whether to subtract true/heterodyned values of phase parameters from the distributions (so true/heterodyned value is at zero)
582 self._showcontours = showcontours # set whether to show probability contours on 2d posterior plots
583
584 # check if parameter file has been given
585 if self._parfile is not None:
586 # try and read it
587 try:
588 self._injection_parameters = psr_par(self._parfile)
589 except:
590 print(
591 "Error... cannot read injection parameter file '%s'."
592 % self._parfile,
593 file=sys.stderr,
594 )
595 sys.exit(1)
596
597 if self._usegwphase: # change initial phase if required
598 if hasattr(self._injection_parameters, "PHI0"):
599 phi0val = 2.0 * self._injection_parameters["PHI0"]
600 setattr(self._injection_parameters, "PHI0", phi0val)
601
602 # if RA_RAD and DEC_RAD are set then copy them into RA and DEC attributes
603 if hasattr(self._injection_parameters, "RA_RAD"):
604 setattr(
605 self._injection_parameters,
606 "RA",
607 self._injection_parameters["RA_RAD"],
608 )
609 if hasattr(self._injection_parameters, "DEC_RAD"):
610 setattr(
611 self._injection_parameters,
612 "DEC",
613 self._injection_parameters["DEC_RAD"],
614 )
615
616 # if DIST is set, then use the original value in KPC as the posterior samples will be in KPC
617 if hasattr(self._injection_parameters, "DIST"):
618 setattr(
619 self._injection_parameters,
620 "DIST",
621 self._injection_parameters["DIST_ORIGINAL"],
622 )
623
624 if "Joint" in self._ifos: # put 'Joint' at the end
625 j = self._ifos.pop(self._ifos.index("Joint"))
626 self._ifos.append(j)
627
628 # check if prior file has been given
629 if self._priorfile is not None:
630 # try and read it
631 try:
632 pf = open(self._priorfile, "r")
633 except:
634 print(
635 "Error... cannot read prior file '%s'." % self._priorfile,
636 file=sys.stderr,
637 )
638 sys.exit(1)
639
640 for line in pf.readlines(): # read in priors
641 priorlinevals = line.split()
642 if len(priorlinevals) < 4:
643 print(
644 "Error... there must be at least four values on each line of the prior file '%s'."
645 % self._priorfile,
646 file=sys.stderr,
647 )
648 sys.exit(1)
649
650 if priorlinevals[1] not in [
651 "uniform",
652 "fermidirac",
653 "gaussian",
654 "gmm",
655 "loguniform",
656 ]:
657 print(
658 "Error... the prior for '%s' must be either 'uniform', 'loguniform', 'gmm', 'fermidirac', or 'gaussian'."
659 % priorlinevals[0],
660 file=sys.stderr,
661 )
662 sys.exit(1)
663
664 if priorlinevals[1] in [
665 "uniform",
666 "fermidirac",
667 "gaussian",
668 "loguniform",
669 ]:
670 if len(priorlinevals) != 4:
671 print(
672 "Error... there must be four values on each line of the prior file '%s'."
673 % self._priorfile,
674 file=sys.stderr,
675 )
676 sys.exit(1)
677 ranges = np.array(
678 [float(priorlinevals[2]), float(priorlinevals[3])]
679 ) # set ranges
680
681 if (
682 self._usegwphase and priorlinevals[0].lower() == "phi0"
683 ): # adjust phi0 priors if using GW phase
684 ranges = 2.0 * ranges
685 elif priorlinevals[1] == "gmm":
686 try:
687 ranges = {}
688 ranges["nmodes"] = priorlinevals[2]
689 ranges["means"] = ast.literal_eval(priorlinevals[3]) # means
690 ranges["covs"] = ast.literal_eval(
691 priorlinevals[4]
692 ) # covariance matrices
693 ranges["weights"] = ast.literal_eval(
694 priorlinevals[5]
695 ) # weights
696 plims = []
697 splitpars = [
698 p.lower() for p in priorlinevals[0].split(":")
699 ] # split GMM parameter names
700 npars = len(splitpars)
701 for lims in priorlinevals[6:-1]: # limits for each parameter
702 plims.append(np.array(ast.literal_eval(lims)))
703 ranges["limits"] = plims
704 ranges["npars"] = npars
705 if self._usegwphase and "phi0" in splitpars:
706 phi0idx = splitpars.index("phi0")
707 # adjust phi0 priors (multiply limits, means and covariances by 2.) if using GW phase
708 for ci in range(ranges["nmodes"]):
709 ranges["means"][ci][phi0idx] = (
710 2.0 * ranges["means"][ci][phi0idx]
711 )
712 for ck in range(
713 npars
714 ): # this will correctly multiply the variances by 4 when ck == phi0idx
715 ranges["covs"][ci][phi0idx][ck] = (
716 2.0 * ranges["covs"][ci][phi0idx][ck]
717 )
718 ranges["covs"][ci][ck][phi0idx] = (
719 2.0 * ranges["covs"][ci][ck][phi0idx]
720 )
721 ranges["limits"][phi0idx] = 2.0 * ranges["limits"][phi0idx]
722 for sidx, spar in enumerate(splitpars):
723 ranges["pindex"] = sidx
724 self._prior_parameters[spar] = {"gmm": ranges}
725 continue
726 except:
727 # if can't read in GMM prior correctly then just continue
728 continue
729
730 self._prior_parameters[priorlinevals[0]] = {priorlinevals[1]: ranges}
731
732 # check posterior files exist
733 for ifo in postfiles:
734 if not os.path.isfile(postfiles[ifo]):
735 print(
736 "Error... posterior samples file '%s' for '%s' does not exist."
737 % (postfiles[ifo], ifo),
738 file=sys.stderr,
739 )
740 sys.exit(1)
741 else: # convert into posterior class
742 pos, sigev, noiseev = pulsar_nest_to_posterior(postfiles[ifo])
743 self._posteriors[ifo] = pos
744 self._signal_evidence[ifo] = sigev
745 self._noise_evidence[ifo] = noiseev
746 self._Bsn[ifo] = sigev - noiseev
747
748 theseparams = []
749 for param in pos.names:
750 thispar = True
751 # remove parameters with 'logl', 'logw' or 'logprior' in the name
752 for testpar in ["logl", "logw", "logprior"]:
753 if testpar in param.lower():
754 thispar = False
755 break
756 if thispar:
757 theseparams.append(param)
758
759 if len(self._parameters) == 0:
760 self._parameters = copy.copy(theseparams)
761 else: # make sure different detectors do not have different sets of parameters
762 for param in theseparams:
763 if param not in self._parameters:
764 print(
765 "Error... parameter '%s' is not defined in posteriors samples for '%s'."
766 % (param, ifo),
767 file=sys.stderr,
768 )
769 sys.exit(1)
770
771 # rotate phi0 and psi into the 0->pi and 0->pi/2 ranges respectively if required (see Eqn. 45 of http://arxiv.org/abs/1501.05832)
772 if modeltype == "source":
773 phi0samples = None
774 psisamples = None
775 if "phi0" in pos.names:
776 phi0samples = pos["phi0"].samples
777 if "psi" in pos.names:
778 psisamples = pos["psi"].samples
779
780 # rotate psi values by pi/2 increments into the 0->pi/2 range
781 if psisamples is not None:
782 for i in range(len(psisamples)):
783 nrots = np.abs(
784 np.floor(psisamples[i] / (np.pi / 2.0))
785 ) # number of rotations to return to range
786 psisamples[i] = np.mod(psisamples[i], np.pi / 2.0)
787 if phi0samples is not None: # rotate phi0 appropriately
788 phi0samples[i] += nrots * (np.pi / 2.0)
789 psisnew = PosteriorOneDPDF("psi", psisamples)
790 pos.pop("psi")
791 pos.append(psisnew)
792
793 # make sure phi0 values are between 0->pi
794 if phi0samples is not None:
795 phi0samples = np.mod(phi0samples, np.pi)
796 phi0new = PosteriorOneDPDF("phi0", phi0samples)
797 pos.pop("phi0")
798 pos.append(phi0new)
799
800 if (
801 self._usegwphase
802 ): # try switching phi0 to 2*phi0 if working with l=m=2 gravitational wave initial phase (e.g. for hardware injections)
803 if "phi0" in pos.names:
804 phi0new = PosteriorOneDPDF("phi0", 2.0 * pos["phi0"].samples)
805 pos.pop("phi0")
806 pos.append(phi0new)
807
808 # if just working with results at 2f, but C22 and phi22 parameters have been used, convert back to h0 and phi0
809 if len(self._harmonics) == 1:
810 if self._harmonics[0] == 2.0:
811 if "c22" in pos.names and self._modeltype == "waveform":
812 h0new = PosteriorOneDPDF("h0", 2.0 * pos["c22"].samples)
813 pos.pop("c22")
814 pos.append(h0new)
815 ih0 = self._parameters.index("h0")
816 self._parameters.pop(ih0)
817 self._parameters.insert(ih0, "C22")
818 if self._parfile is not None:
819 if not hasattr(self._injection_parameters, "H0"):
820 setattr(
821 self._injection_parameters,
822 "H0",
823 2.0 * self._injection_parameters["C22"],
824 )
825
826 if "phi22" in pos.names and modeltype == "waveform":
827 phi0new = PosteriorOneDPDF(
828 "phi0", 0.5 * pos["phi22"].samples
829 )
830 pos.pop("phi0")
831 pos.append(phi0new)
832 iphi0 = self._parameters.index("phi0")
833 self._parameters.pop(iphi0)
834 self._parameters.insert(iphi0, "PHI0")
835 if self._parfile is not None:
836 if not hasattr(self._injection_parameters, "PHI0"):
837 setattr(
838 self._injection_parameters,
839 "PHI0",
840 0.5 * self._injection_parameters["PHI22"],
841 )
842
843 # for a biaxial star (and if using the waveform model) set phi22 = 2*phi21
844 if len(self._harmonics) == 2:
845 if (
846 1.0 in self._harmonics
847 and 2.0 in self._harmonics
848 and self._biaxial
849 and self._modeltype == "waveform"
850 ):
851 if "phi21" in pos.names:
852 phi22 = PosteriorOneDPDF(
853 "phi22", 2.0 * pos["phi21"].samples
854 )
855 pos.append(phi22)
856 if self._parfile is not None:
857 if hasattr(self._injection_parameters, "PHI21"):
858 setattr(
859 self._injection_parameters,
860 "PHI22",
861 2.0 * self._injection_parameters["PHI21"],
862 )
863
864 self._posteriors[ifo] = pos
865
866 self._get_snrs() # get the optimal SNR of the signal
867
868 self._get_bayes_factors() # get the Bayes factors
869
870 @property
871 def parameters(self):
872 # return a list of the source parameter names
873 return self._parameters
874
875 @property
876 def bsn(self):
877 # return a dictionary of the signal vs noise log Bayes factors
878 return self._Bsn
879
880 @property
881 def snrs(self):
882 # return dictionary of snrs
883 return self._optimal_snrs
884
885 @property
886 def bci(self):
887 # return coherent vs incoherent Bayes factor
888 return self._Bci
889
890 @property
891 def bcin(self):
892 # return coherent vs incoherent or noise Bayes factor
893 return self._Bcin
894
895 def h0_ul(self, ifo):
896 # return the h0 upper limit for a given detector
897 if ifo in self._h0ul:
898 return self._h0ul[ifo]
899 else:
900 return None
901
902 def ellipticity_ul(self, ifo):
903 # return the ellipticity upper limit
904 if ifo in self._ellipticity:
905 return self._ellipticity[ifo]
906 else:
907 return None
908
909 def q22_ul(self, ifo):
910 # return the Q22 quadrupole moment upper limit
911 if ifo in self._q22:
912 return self._q22[ifo]
913 else:
914 return None
915
916 def C21_ul(self, ifo):
917 # return the C21 upper limit
918 if ifo in self._C21:
919 return self._C21[ifo]
920 else:
921 return None
922
923 def C22_ul(self, ifo):
924 # return the C22 upper limit
925 if ifo in self._C22:
926 return self._C22[ifo]
927 else:
928 return None
929
930 def I21_ul(self, ifo):
931 # return the I21 upper limit
932 if ifo in self._I21:
933 return self._I21[ifo]
934 else:
935 return None
936
937 def I31_ul(self, ifo):
938 # return the I31 upper limit
939 if ifo in self._I31:
940 return self._I31[ifo]
941 else:
942 return None
943
944 def sdlim_ratio(self, ifo):
945 # return the spin-down limit ratio
946 if ifo in self._sdratio:
947 return self._sdratio[ifo]
948 else:
949 return None
950
951 def _get_snrs(self):
952 # get the returned match filter SNRs
953 for ifo in self._ifos:
954 self._optimal_snrs[ifo] = self.get_snr(
955 os.path.dirname(self._postfiles[ifo])
956 )
957
958 def get_snr(self, pdir):
959 # get SNR from files in pdir
960 snr = 0.0
961 # get files with SNR in the name in the posterior file directory
962 snrfiles = [sf for sf in os.listdir(pdir) if "SNR" in sf]
963 if len(snrfiles) < 1:
964 print("Error... no SNR files are given for '%s'" % ifo, file=sys.stderr)
965 sys.exit(1)
966
967 for (
968 snrfile
969 ) in (
970 snrfiles
971 ): # average SNR values from mulitple files (e.g. if nested sampling has been run multiple times for the same signal)
972 fp = open(os.path.join(pdir, snrfile), "r")
973 lines = [line.strip() for line in fp.readlines()]
974
975 if "# Recovered SNR" not in lines:
976 print(
977 "Error... no recovered SNRs are given in the SNR file '%s'."
978 % snrfile,
979 file=sys.stderr,
980 )
981 sys.exit(1)
982
983 # just return the final (either coherent or single detector-single frequency value of the SNR)
984 linevals = lines[-1].split()
985 snr += float(linevals[-1])
986 return snr / len(snrfiles)
987
988 def _get_bayes_factors(self):
989 # get the Bayes factors (actually odds ratios with equal priors for all hypotheses) for the signal
990 nifos = 0
991 ifosn = [] # list of signal and noise evidences for each detector
992 for ifo in self._ifos:
993 (
994 self._Bsn[ifo],
995 self._signal_evidence[ifo],
996 self._noise_evidence[ifo],
997 self._maxL[ifo],
998 ) = self.get_bayes_factor(self._postfiles[ifo])
999
1000 if ifo != "Joint":
1001 nifos += 1
1002 ifosn.append(
1003 {"s": self._signal_evidence[ifo], "n": self._noise_evidence[ifo]}
1004 )
1005
1006 # get all combinations of (incoherent) noise and signal hypotheses
1007 combs = [
1008 list(i) for i in itertools.product(["s", "n"], repeat=nifos)
1009 ] # see e.g. http://stackoverflow.com/q/14931769/1862861
1010 incoherentcombs = -np.inf
1011 incoherentsig = 0.0 # incoherent signal in all detectors
1012 for comb in combs:
1013 # don't include the all noise hypotheses (as we have that already as self._noise_evidence['Joint'])
1014 if comb.count("n") != len(
1015 comb
1016 ): # see e.g. http://stackoverflow.com/a/3844948/1862861
1017 combsum = 0.0
1018 for i, cval in enumerate(comb):
1019 combsum += ifosn[i][cval]
1020 incoherentcombs = np.logaddexp(incoherentcombs, combsum)
1021 if comb.count("s") == len(comb):
1022 incoherentsig = combsum
1023
1024 # get the coherent vs incoherent noise odds ratio (assuming all hypotheses have equal priors)
1025 if len(self._ifos) > 2 and "Joint" in self._ifos:
1026 self._Bci = self._signal_evidence["Joint"] - incoherentsig
1027 self._Bcin = self._signal_evidence["Joint"] - np.logaddexp(
1028 incoherentcombs, self._noise_evidence["Joint"]
1029 )
1030
1031 def get_bayes_factor(self, postfile):
1032 # return the Bayes factor extracted from a posterior file
1033 try:
1034 fe = os.path.splitext(postfile)[-1].lower() # file extension
1035 if fe == ".h5" or fe == ".hdf": # HDF5 file
1036 # open hdf5 file
1037 f = h5py.File(postfile, "r")
1038 a = f["lalinference"]["lalinference_nest"]
1039 evdata = (
1040 a.attrs["log_bayes_factor"],
1041 a.attrs["log_evidence"],
1042 a.attrs["log_noise_evidence"],
1043 a.attrs["log_max_likelihood"],
1044 )
1045 f.close()
1046 else: # try old/legacy file format
1047 B = np.loadtxt(postfile.replace(".gz", "") + "_B.txt")
1048 evdata = tuple(B.tolist())
1049 except:
1050 print(
1051 "Error... could not extract evidences from '%s'." % postfile,
1052 file=sys.stderr,
1053 )
1054 sys.exit(1)
1055
1056 return (
1057 evdata # line contains Bsn, signal evidence, noise evidence, max likelihood
1058 )
1059
1060 def snr(self, ifo):
1061 # return the SNR for a given detector
1062 if ifo in self._optimal_snrs:
1063 return self._optimal_snrs[ifo]
1064 else:
1065 return None
1066
1067 def create_joint_posterior_plot(
1068 self,
1069 parameters,
1070 bins=20,
1071 ifo=None,
1072 truths=None,
1073 credintervals=[0.9],
1074 filename=None,
1075 figformats=["png"],
1076 ratio=3,
1077 figlimits=None,
1078 contourlimits=None,
1079 jointsamples=True,
1080 whichtruth=None,
1081 scatter_kwargs={},
1082 ):
1083 # create a plot with the 1D and 2D joint posteriors (if ifo is None then use all ifos in class)
1084 if ifo != None and ifo in self._ifos:
1085 plotifos = [ifo]
1086 else:
1087 plotifos = self._ifos
1088
1089 # if a joint detector posterior is present, make that first so other data is plotted on top
1090 if "Joint" in plotifos:
1091 j = plotifos.pop(plotifos.index("Joint"))
1092 plotifos.insert(0, j)
1093
1094 coldict = {
1095 "H1": "red",
1096 "H2": "cyan",
1097 "L1": "green",
1098 "V1": "blue",
1099 "G1": "magenta",
1100 "Joint": "black",
1101 }
1102
1103 if len(parameters) < 2 or len(parameters) > len(self._parameters):
1104 print(
1105 "Error... can only plot posterior distributions for more than one parameters",
1106 file=sys.stderr,
1107 )
1108 sys.exit(1)
1109
1110 thistruth = []
1111 for param in parameters:
1112 if param not in self._parameters:
1113 print(
1114 "Error... the requested parameter '%s' is not recognised" % param,
1115 file=sys.stderr,
1116 )
1117 sys.exit(1)
1118
1119 if (
1120 self._injection_parameters is not None and truths is None
1121 ): # we have injection values
1122 thistruth.append(self._injection_parameters[param.upper()])
1123
1124 if self._injection_parameters is not None and truths is None:
1125 truths = thistruth
1126
1127 if truths is not None:
1128 if isinstance(truths, list): # just a list of the true parameter values
1129 newtruths = {}
1130 if len(truths) != len(parameters):
1131 print(
1132 "Error... number of true values must be the same as number of parameters.",
1133 file=sys.stderr,
1134 )
1135 sys.exit(1)
1136 else:
1137 # convert to dictionary
1138 for ifo in plotifos:
1139 newtruths[ifo] = truths
1140 truths = newtruths
1141 if isinstance(
1142 truths, dict
1143 ): # dictionary keyed to IFO with lists of true parameters for each IFO (this is really designed for the SNR vs Bayes factor plots)
1144 for ifo in plotifos:
1145 if ifo not in truths:
1146 print(
1147 "Error... problem in truths values. Detector '%s' not recognised."
1148 % ifo,
1149 file=sys.stderr,
1150 )
1151 sys.exit(1)
1152 else:
1153 if len(truths[ifo]) != len(parameters):
1154 print(
1155 "Error... number of true values must be the same as number of parameters.",
1156 file=sys.stderr,
1157 )
1158 sys.exit(1)
1159 else:
1160 truths = {}
1161 for ifo in plotifos:
1162 truths[ifo] = None # set all to be None
1163
1164 # list of truth values to subtract
1165 subtract_truths = []
1166 amppars = [
1167 "h0",
1168 "c21",
1169 "c22",
1170 "i21",
1171 "i31",
1172 "cosiota",
1173 "iota",
1174 "phi0",
1175 "phi21",
1176 "phi22",
1177 "lambda",
1178 "costheta",
1179 "theta",
1180 "psi",
1181 ] # a list of "amplitude" parameters for which subtract truths won't be applied
1182 if truths is not None and self._subtract_truths is not None:
1183 if parameters[0] not in amppars:
1184 subtract_truths.append(0) # add the true/heterodyned parameter
1185
1186 # use first ifo and get the required posterior samples
1187 x = self._posteriors[plotifos[0]][parameters[0]].samples
1188 labels = []
1189 if parameters[0].upper() in paramlatexdict:
1190 labels.append(paramlatexdict[parameters[0].upper()])
1191 else:
1192 labels.append(parameters[0])
1193 for param in parameters[1:]:
1194 x = np.hstack((x, self._posteriors[plotifos[0]][param].samples))
1195 if param.upper() in paramlatexdict:
1196 labels.append(paramlatexdict[param.upper()])
1197 else:
1198 labels.append(param)
1199 if truths is not None and self._subtract_truths is not None:
1200 if param not in amppars:
1201 subtract_truths.append(parameters.index(param))
1202
1203 if len(subtract_truths) == 0:
1204 subtract_truths = None
1205
1206 # set styles to different detectors
1207 if plotifos[0] == "Joint":
1208 histops = {
1209 "histtype": "stepfilled",
1210 "color": "darkslategrey",
1211 "edgecolor": coldict["Joint"],
1212 "linewidth": 1.5,
1213 }
1214 contourops = {"colors": coldict[plotifos[0]]}
1215 showcontours = self._showcontours
1216 if whichtruth == "Joint" or whichtruth == "all":
1217 truthops = {"color": "black", "markeredgewidth": 2}
1218 else:
1219 truthops = {}
1220 showpoints = jointsamples
1221 else:
1222 showcontours = self._showcontours
1223 contourops = {"colors": "dark" + coldict[plotifos[0]]}
1224 showpoints = True
1225 if len(plotifos) == 1: # if just one detector use a filled histogram
1226 histops = {
1227 "histtype": "stepfilled",
1228 "color": coldict[plotifos[0]],
1229 "edgecolor": coldict[plotifos[0]],
1230 "linewidth": 1.5,
1231 }
1232 truthops = {"color": "black", "markeredgewidth": 2}
1233 else:
1234 histops = {
1235 "histtype": "step",
1236 "color": coldict[plotifos[0]],
1237 "edgecolor": coldict[plotifos[0]],
1238 "linewidth": 1,
1239 }
1240 if whichtruth == plotifos[0]:
1241 truthops = {"color": "black", "markeredgewidth": 2}
1242 elif whichtruth == "all":
1243 truthops = {
1244 "color": "dark" + coldict[plotifos[0]],
1245 "markeredgewidth": 2,
1246 }
1247 else:
1248 truthops = {}
1249
1250 sc = scotchcorner(
1251 x,
1252 bins=bins,
1253 ratio=ratio,
1254 labels=labels,
1255 truths=truths[plotifos[0]],
1256 datatitle=plotifos[0],
1257 showlims="both",
1258 hist_kwargs=histops,
1259 showcontours=showcontours,
1260 contour_levels=credintervals,
1261 contour_kwargs=contourops,
1262 truths_kwargs=truthops,
1263 contour_limits=contourlimits,
1264 limits=figlimits,
1265 show_level_labels=False,
1266 showpoints=showpoints,
1267 scatter_kwargs=scatter_kwargs,
1268 subtract_truths=subtract_truths,
1269 )
1270
1271 # now add the rest to the plots
1272 if len(plotifos) > 1:
1273 for k, ifo in enumerate(plotifos[1:]):
1274 histops = {
1275 "histtype": "step",
1276 "color": coldict[ifo],
1277 "edgecolor": coldict[ifo],
1278 "linewidth": 1,
1279 }
1280 x = self._posteriors[ifo][parameters[0]].samples
1281 for param in parameters[1:]:
1282 x = np.hstack((x, self._posteriors[ifo][param].samples))
1283 showcontours = self._showcontours
1284 contourops = {"colors": "dark" + coldict[plotifos[k + 1]]}
1285 if whichtruth == plotifos[k + 1]:
1286 truthops = {"color": "black", "markeredgewidth": 2}
1287 elif whichtruth == "all":
1288 truthops = {
1289 "color": "dark" + coldict[plotifos[k + 1]],
1290 "markeredgewidth": 2,
1291 }
1292 else:
1293 truthops = {}
1294 sc.add_data(
1295 x,
1296 hist_kwargs=histops,
1297 datatitle=ifo,
1298 truths=truths[ifo],
1299 showcontours=showcontours,
1300 contour_kwargs=contourops,
1301 contour_levels=credintervals,
1302 show_level_labels=False,
1303 truths_kwargs=truthops,
1304 scatter_kwargs=scatter_kwargs,
1305 contour_limits=contourlimits,
1306 limits=figlimits,
1307 )
1308
1309 # add priors plots if required
1310 if self._priorfile is not None and len(self._prior_parameters) > 0:
1311 for priorparam in self._prior_parameters:
1312 if priorparam.lower() in parameters:
1313 # get axes for that parameter
1314 thisax = sc.get_axis(labels[parameters.index(priorparam.lower())])
1315
1316 atruth = None
1317 if (
1318 truths is not None
1319 and self._subtract_truths is not None
1320 and priorparam.lower() not in amppars
1321 ):
1322 atruth = truths[plotifos[0]][
1323 parameters.index(priorparam.lower())
1324 ]
1325
1326 # check if this is the vertical histogram or not
1327 vertaxrange = sc.histvert[-1].get_ylim()
1328 yl = thisax.get_ylim()
1329 if (
1330 yl[0] == vertaxrange[0] and yl[1] == vertaxrange[1]
1331 ): # vertical histogram
1332 self.plot_prior(
1333 thisax,
1334 priorparam,
1335 self._prior_parameters,
1336 truth=atruth,
1337 orientation="vertical",
1338 )
1339 else:
1340 self.plot_prior(
1341 thisax,
1342 priorparam,
1343 self._prior_parameters,
1344 truth=atruth,
1345 orientation="horizontal",
1346 )
1347
1348 # output the plots
1349 if "png" not in figformats and "svg" not in figformats:
1350 print(
1351 "Error... must include 'png' and/or 'svg' in the output figure formats",
1352 file=sys.stderr,
1353 )
1354 sys.exit(1)
1355
1356 if filename == None: # use default file names
1357 if len(parameters) == len(self._parameters):
1358 outfilepre = os.path.join(self._outputdir, "all_posteriors")
1359 else:
1360 outfilepre = os.path.join(self._outputdir, "vs".join(parameters))
1361 else:
1362 outfilepre = os.path.join(self._outputdir, filename)
1363
1364 outfiles = []
1365 for ftype in figformats:
1366 outfile = outfilepre + "." + ftype
1367 try:
1368 sc.fig.subplots_adjust(
1369 left=0.18, bottom=0.18
1370 ) # adjust size to accommodate axes labels
1371 sc.savefig(outfile)
1372 except:
1373 print(
1374 "Error... could not output posterior plot file '%s'." % outfile,
1375 file=sys.stderr,
1376 )
1377 sys.exit(1)
1378 outfiles.append(outfile)
1379
1380 return outfiles # list of output figure file names
1381
1382 def plot_prior(
1383 self, ax, param, prior, orientation="horizontal", truth=0.0, npoints=100
1384 ):
1385 # plot the prior distribution (with truth subtracted if non-zero)
1386 priortype = list(prior[param].keys())[0]
1387 priorrange = list(prior[param].values())[0]
1388
1389 if truth is None:
1390 truth = 0.0
1391
1392 if orientation == "horizontal":
1393 valrange = np.linspace(ax.get_xlim()[0], ax.get_xlim()[1], npoints)
1394 elif orientation == "vertical":
1395 valrange = np.linspace(ax.get_ylim()[0], ax.get_ylim()[1], npoints)
1396 else:
1397 print(
1398 "Error... axis orientation not recognised. Must be 'horizontal' or 'vertical'.",
1399 file=sys.stderr,
1400 )
1401 sys.exit(1)
1402
1403 # get the prior distributions
1404 if priortype == "uniform":
1405 vals = stats.uniform.pdf(
1406 valrange, priorrange[0] - truth, priorrange[1] - priorrange[0]
1407 )
1408 elif priortype == "gaussian":
1409 # crude (not taking account of wrap-around) shift of phi0 and psi into 0->pi and 0->pi/2 ranges)
1410 if param.lower() == "psi":
1411 priorrange[0] = np.mod(priorrange[0], np.pi / 2.0)
1412 if param.lower() == "phi0":
1413 if self._usegwphase:
1414 priorrange[0] = np.mod(priorrange[0], 2.0 * pi)
1415 else:
1416 priorrange[0] = np.mod(priorrange[0], pi)
1417
1418 vals = stats.norm.pdf(valrange, priorrange[0] - truth, priorrange[1])
1419 elif (
1420 priortype == "fermidirac"
1421 ): # don't subtract truth from Fermi-Dirac as it should only be use for amplitude parameters anyway
1422 sigma = priorrange[0]
1423 r = priorrange[1]
1424 mu = sigma * r
1425 vals = 1.0 / (
1426 (sigma * np.log(1.0 + np.exp(r)))
1427 * (np.exp((valrange - mu) / sigma) + 1.0)
1428 )
1429 elif priortype == "loguniform":
1430 vals = np.zeros(len(valrange))
1431 indices = (valrange >= priorrange[0]) & (valrange <= priorrange[1])
1432 vals[indices] = 1.0 / (
1433 valrange[indices] * np.log(priorrange[1] / priorrange[0])
1434 )
1435 elif priortype == "gmm":
1436 vals = np.zeros(len(valrange))
1437 nmodes = priorrange["nmodes"] # number of modes in GMM
1438 pindex = priorrange["pindex"] # parameter index in GMM prior
1439 mmeans = np.array([mm[pindex] for mm in priorrange["means"]]) # means
1440 mcovs = priorrange["covs"] # covariance matrices
1441 mstddevs = []
1442 for cov in mcovs:
1443 mstddevs.append(np.sqrt(cov[pindex][pindex]))
1444 mweights = priorrange["weights"]
1445 # check bounds
1446 bmin = -np.inf
1447 bmax = np.inf
1448 plims = priorrange["limits"][pindex]
1449 if np.isfinite(plims[0]):
1450 bmin = plims[0]
1451 if np.isfinite(plims[1]):
1452 bmax = plims[1]
1453 indices = (valrange >= bmin) & (valrange <= bmax)
1454
1455 # crude (not taking account of wrap-around) shift of phi0 and psi into 0->pi and 0->pi/2 ranges)
1456 if param.lower() == "psi":
1457 mmeans = np.mod(mmeans, np.pi / 2.0)
1458 if param.lower() == "phi0":
1459 if self._usegwphase:
1460 mmeans = np.mod(mmeans, 2.0 * pi)
1461 else:
1462 mmeans = np.mod(mmeans, pi)
1463
1464 # create Gaussian mixture model
1465 for i in range(nmodes):
1466 vals[indices] += mweights[i] * stats.norm.pdf(
1467 valrange[indices], mmeans[i] - truth, mstddevs[i]
1468 )
1469 # normalise
1470 vals = vals / np.trapz(vals, valrange)
1471 else:
1472 print(
1473 "Error... prior type '%s' not recognised." % priortype, file=sys.stderr
1474 )
1475 sys.exit(1)
1476
1477 # make plots
1478 if orientation == "horizontal":
1479 ax.plot(
1480 valrange, vals, linestyle="--", color="lightslategray", linewidth=1.5
1481 )
1482 if orientation == "vertical":
1483 ax.plot(
1484 vals, valrange, linestyle="--", color="lightslategray", linewidth=1.5
1485 )
1486
1487 # create standard joint plots
1488 def create_joint_plots_table(self, allparams=False, title="Joint distributions"):
1489 # create a table containing a set of standard 2D or 3D joint posterior plots (depending on the analysis setup)
1490 header = htmltag("h2", title, newline=True)
1491
1492 table = htmltable()
1493 table.addrow()
1494 paramlist = []
1495
1496 # set limits for credible interval contours plots using maximum allowed ranges (see e.g. Tables 1 and 2 of http://arxiv.org/abs/1508.00416v2)
1497 limits = {}
1498 for param in self._parameters:
1499 if param == "h0":
1500 # only set lower limit on amplitudes if posteriors are close (less than 3 sigma) to zero
1501 limits[param] = ()
1502 for ifo in self._ifos:
1503 if (
1504 self._posteriors[ifo][param].mean
1505 - 3.0 * self._posteriors[ifo][param].stdev
1506 < 0.0
1507 ):
1508 limits[param] = (0.0, None) # must be greater than 0.
1509 break
1510 elif param == "cosiota":
1511 limits[param] = (-1.0, 1.0) # cos(iota)
1512 elif param == "phi0":
1513 if self._usegwphase:
1514 limits[param] = (0.0, 2.0 * np.pi)
1515 else:
1516 limits[param] = (0.0, np.pi)
1517 elif param == "psi":
1518 if self._biaxial and self._modeltype == "source":
1519 limits[param] = (0.0, np.pi)
1520 else:
1521 limits[param] = (0.0, np.pi / 2.0)
1522 elif param == "c21" or param == "c22":
1523 if self._biaxial and self._modeltype == "waveform":
1524 limits[param] = ()
1525 else:
1526 limits[param] = ()
1527 for ifo in self._ifos:
1528 if (
1529 self._posteriors[ifo][param].mean
1530 - 3.0 * self._posteriors[ifo][param].stdev
1531 < 0.0
1532 ):
1533 limits[param] = (0.0, None) # must be greater than 0.
1534 break
1535 elif param == "lambda":
1536 limits[param] = (0.0, np.pi)
1537 elif param == "costheta":
1538 limits[param] = (0.0, 1.0)
1539 else: # default to empty tuple
1540 limits[param] = () # empty tuple
1541
1542 if allparams: # plot all parameters
1543 paramlist = [self._parameters]
1544 else:
1545 if len(self._harmonics) == 1:
1546 if harmonics[0] == 2.0: # for emission at twice the rotation frequency
1547 # plot h0 versus cos(iota) and psi verus phi0
1548 paramlist = [["h0", "cosiota"], ["phi0", "psi"]]
1549 else:
1550 if self._modeltype == "source":
1551 print(
1552 "Error... do not know 'source' parameterisation for emission at purely the rotation frequency",
1553 file=sys.stderr,
1554 )
1555 sys.exit(1)
1556 else:
1557 # plot C21 versus cosiota and psi versus phi21
1558 paramlist = [["c21", "cosiota"], ["phi21", "psi"]]
1559 else:
1560 if self._modeltype == "source":
1561 if self._biaxial:
1562 # plot I31 vs cosiota vs phi0 and psi vs lambda vs theta (I21 is zero)
1563 paramlist = [
1564 ["i31", "cosiota", "phi0"],
1565 ["psi", "lambda", "costheta"],
1566 ]
1567 else:
1568 # plot I21 vs I31 vs cosiota and phi0 vs psi vs lambda vs theta
1569 paramlist = [
1570 ["i21", "i31", "cosiota"],
1571 ["phi0", "psi", "lambda", "costheta"],
1572 ]
1573 else:
1574 # plot C21 vs C22 vs cosiota and psi vs phi21 vs phi22
1575 paramlist = [["c21", "c22", "cosiota"], ["psi", "phi21", "phi22"]]
1576
1577 notpresent = False
1578 for parampairs in paramlist:
1579 contourlimits = []
1580 figlimits = []
1581 for p in parampairs:
1582 if p not in self._parameters:
1583 notpresent = True
1584 else:
1585 contourlimits.append(limits[p])
1586 figlimits.append(limits[p])
1587 if notpresent:
1588 break
1589
1590 pf = self.create_joint_posterior_plot(
1591 parampairs,
1592 figformats=["png"],
1593 ratio=2,
1594 figlimits=figlimits,
1595 contourlimits=contourlimits,
1596 jointsamples=False,
1597 )
1598 if allparams:
1599 tagclass = "jointplot"
1600 else:
1601 tagclass = "posplot"
1602 table.adddata(
1603 atag(
1604 os.path.basename(pf[0]),
1605 '<img class="{}" src="{}"/>'.format(
1606 tagclass, os.path.basename(pf[0])
1607 ),
1608 ).text
1609 )
1610
1611 return header.text + table.tabletext
1612
1613 def create_sample_plot_table(self, figformats=["png"]):
1614 # create the table of posterior sample plots
1615 if "png" not in figformats:
1616 print(
1617 "Error... must include 'png' in the output figure formats",
1618 file=sys.stderr,
1619 )
1620 sys.exit(1)
1621
1622 header = htmltag(
1623 "h2", "Posterior samples", newline=True
1624 ) # the output html table
1625
1626 table = htmltable()
1627
1628 # get the plots
1629 for param in self._parameters:
1630 chainfig = plot_posterior_chain(
1631 self._posteriors.values(), param, self._ifos, grr=None
1632 )
1633 chainfile = os.path.join(self._outputdir, param + "_postchain")
1634
1635 for ftype in figformats:
1636 thischainfile = chainfile + "." + ftype
1637 try:
1638 chainfig.savefig(thischainfile)
1639 except:
1640 print(
1641 "Error... could not output posterior samples plot for '%s'."
1642 % param,
1643 file=sys.stderr,
1644 )
1645 sys.exit(1)
1646
1647 if ftype == "png":
1648 table.addrow()
1649 table.adddata(
1650 atag(
1651 os.path.basename(thischainfile),
1652 '<img class="chainplot" src="{}"/>'.format(
1653 os.path.basename(thischainfile)
1654 ),
1655 ).text
1656 )
1657
1658 self._samples_table = header.text + table.tabletext
1659 return self._samples_table # output html for table
1660
1661 def create_stats_table(self, credints=[95]):
1662 # create table of posterior sample statistics, including credible intervals given by the list
1663 self._credint = credints
1664 for ci in self._credint:
1665 if ci >= 100.0 or ci <= 0.0:
1666 print(
1667 "Error... credible interval '%s' is outside the range between 0%% and 100%%."
1668 % str(ci),
1669 file=sys.stderr,
1670 )
1671 sys.exit(1)
1672
1673 header = htmltag("h2", "Posterior statistics", newline=True)
1674
1675 table = htmltable()
1676
1677 cols = ["Detector", "max. poterior", "mean", "median", "&sigma;"]
1678 for ci in self._credint:
1679 cols.append("%d%% credible interval" % int(ci))
1680 if self._injection_parameters != None:
1681 cols.insert(0, "Inj. value") # prepend the injection value
1682 cols.insert(0, "&nbsp;") # add empty cell as header line for the parameter name
1683
1684 # create header row
1685 table.addrow()
1686 for col in cols:
1687 table.adddata(col, header=True)
1688
1689 # loop through parameters header row to table
1690 for i, param in enumerate(self._parameters):
1691 pu = param.upper()
1692 dispkwargs = (
1693 {}
1694 ) # any additional arguments required by value display function
1695
1696 # parameter html to display
1697 if pu in paramhtmldict:
1698 pdisp = paramhtmldict[pu]
1699
1700 # some different display styles (compared to the defaults) for certain parameters
1701 if pu in ["RA", "DEC", "RAJ", "DECJ"]:
1702 dispkwargs = {
1703 "stype": "rads"
1704 } # display as rad rather than hh/dd:mm:ss string
1705 if pu == "F0":
1706 dispkwargs = {"dp": 2} # just display with 2 decimal places
1707 else:
1708 pdisp = param
1709
1710 # parameter value function to display
1711 if pu in paramhtmldispfunc.__dict__:
1712 dispfunc = paramhtmldispfunc.__dict__[pu]
1713 else:
1714 dispfunc = paramhtmldispfunc.__dict__["DEFAULTSTR"]
1715
1716 for j, ifo in enumerate(self._ifos):
1717 table.addrow()
1718 if j == 0: # add parameter name and injection value first)
1719 table.adddata(pdisp, rowspan=len(self._ifos))
1720 if self._injection_parameters != None:
1721 if hasattr(self._injection_parameters, param.upper()):
1722 table.adddata(
1723 dispfunc(
1724 str(self._injection_parameters[param.upper()]),
1725 **dispkwargs,
1726 ),
1727 rowspan=len(self._ifos),
1728 )
1729 else:
1730 table.adddata(
1731 "*", rowspan=len(self._ifos)
1732 ) # no value for this parameter
1733
1734 if i == 0:
1735 self._injection_credible_regions[
1736 ifo
1737 ] = (
1738 {}
1739 ) # dictionary for injection minimal credible regions to be input for each parameter
1740 maxL, maxLparams = self._posteriors[ifo].maxL
1741 table.adddata(ifo, dataclass=ifo)
1742 table.adddata(
1743 dispfunc(str(maxLparams[param]), **dispkwargs)
1744 ) # maximum likelihood
1745 table.adddata(
1746 dispfunc(str(self._posteriors[ifo].means[param]), **dispkwargs)
1747 ) # mean value
1748 table.adddata(
1749 dispfunc(str(self._posteriors[ifo].medians[param]), **dispkwargs)
1750 ) # median value
1751 tdispkwargs = dispkwargs
1752 if param.upper() in ["T0", "TASC"]:
1753 tdispkwargs = {"stype": "diff"} # don't display values in MJD
1754 table.adddata(
1755 dispfunc(str(self._posteriors[ifo].stdevs[param]), **tdispkwargs)
1756 ) # standard deviations
1757
1758 for k, ci in enumerate(credints):
1759 paramval = None
1760 if k == 0 and self._injection_parameters != None:
1761 if hasattr(self._injection_parameters, param.upper()):
1762 paramval = self._injection_parameters[param.upper()]
1763 low, high, cr = self.credible_interval(ifo, param, ci, paramval)
1764
1765 if k == 0 and self._injection_parameters != None:
1766 if hasattr(self._injection_parameters, param.upper()):
1767 self._injection_credible_regions[ifo][param] = cr
1768
1769 table.adddata(
1770 "({}, {})".format(
1771 dispfunc(str(low), **dispkwargs),
1772 dispfunc(str(high), **dispkwargs),
1773 )
1774 )
1775
1776 self._stats_section = header.text + table.tabletext
1777 return self._stats_section
1778
1779 def create_limits_table(self, freq, sdlim=None, dist=None, ul=95):
1780 # create a table with upper limits on amplitudes and Bayes factors
1781 self._limits_table = ""
1782
1783 table = htmltable()
1784
1785 self._h0ul = {}
1786 self._C21 = {}
1787 self._C22 = {}
1788 self._I21 = {}
1789 self._I31 = {}
1790 self._sdratio = {}
1791 self._ellipticity = {}
1792 self._q22 = {}
1793
1794 # add header row
1795 table.addrow()
1796
1797 table.adddata("&nbsp;", header=True) # nothing in the first column
1798
1799 if len(self._harmonics) == 1 and self._harmonics[0] == 2:
1800 table.adddata(paramhtmldict["H0UL"].format(ul), header=True)
1801 table.adddata(paramhtmldict["ELL"], header=True)
1802 table.adddata(paramhtmldict["Q22"], header=True)
1803 table.adddata(paramhtmldict["SDRAT"], header=True)
1804 elif (
1805 len(self._harmonics) == 1
1806 and self._harmonics[0] == 1
1807 and self._modeltype == "waveform"
1808 ):
1809 table.adddata(paramhtmldict["C21UL"].format(ul), header=True)
1810 elif len(self._harmonics) == 2 and self._modeltype == "waveform":
1811 table.adddata(paramhtmldict["C21UL"].format(ul), header=True)
1812 table.adddata(paramhtmldict["C22UL"].format(ul), header=True)
1813 table.adddata(
1814 paramhtmldict["H0UL"].format(ul) + " (from C<sub>22</sub>)", header=True
1815 )
1816 table.adddata(paramhtmldict["ELL"] + " (from C<sub>22</sub>)", header=True)
1817 table.adddata(paramhtmldict["Q22"] + " (from C<sub>22</sub>)", header=True)
1818 table.adddata("ratio", header=True)
1819 elif len(self._harmonics) == 2 and self._modeltype == "source":
1820 if not self._biaixal:
1821 table.adddata(
1822 paramhtmldict["I21UL"].format(ul), header=True
1823 ) # not needed for a biaxial source
1824 table.adddata(paramhtmldict["I31UL"].format(ul), header=True)
1825 else:
1826 print(
1827 "Error... do not know how to output results table for this situation.",
1828 file=sys.stderr,
1829 )
1830 sys.exit(1)
1831
1832 # add SNR
1833 table.adddata(paramhtmldict["SNR"], header=True)
1834
1835 # add max log likelihood
1836 table.adddata(paramhtmldict["MAXL"], header=True)
1837
1838 # add Bayes factor (signal vs noise)
1839 table.adddata(paramhtmldict["BSN"], header=True)
1840
1841 # check whether to add coherent vs incoherent Bayes factor (and coherent vs (incoherent or noise) Bayes factor)
1842 if len(self._ifos) > 2 and "Joint" in self._ifos:
1843 table.adddata(paramhtmldict["BCI"], header=True)
1844 table.adddata(paramhtmldict["BCIN"], header=True)
1845
1846 for ifo in self._ifos:
1847 table.addrow()
1848 table.adddata("{}".format(ifo), dataclass=ifo)
1849
1850 if self._modeltype == "waveform" and len(self._harmonics) == 2:
1851 for p in ["c21", "c22"]:
1852 Cul = upper_limit_greedy(
1853 self._posteriors[ifo][p].samples, upperlimit=(ul / 100.0)
1854 )
1855 table.adddata("{}".format(exp_str(Cul)), dataclass=ifo)
1856 if p == "c21":
1857 self._C21[ifo] = Cul
1858 if p == "c22":
1859 self._C22[ifo] = Cul
1860
1861 self._h0ul[ifo] = self._C22[ifo] * 2.0
1862 table.adddata(
1863 "{}".format(exp_str(self._h0ul[ifo])), dataclass=ifo
1864 ) # convert to h0
1865
1866 if (len(self._harmonics) == 1 and self._harmonics[0] == 2) or (
1867 len(self._harmonics) == 2 and self._modeltype == "waveform"
1868 ):
1869 # get h0 upper limit
1870 if len(self._harmonics) == 1 and self._harmonics[0] == 2:
1871 h0ul = upper_limit_greedy(
1872 self._posteriors[ifo]["h0"].samples, upperlimit=(ul / 100.0)
1873 )
1874 self._h0ul[ifo] = h0ul
1875 table.adddata("{}".format(exp_str(h0ul)), dataclass=ifo)
1876
1877 # ellipticity
1878 if dist == None:
1879 self._ellipticity[ifo] = None
1880 table.adddata("*", dataclass=ifo)
1881 else:
1882 self._ellipticity[ifo] = h0_to_ellipticity(
1883 self._h0ul[ifo], freq, dist
1884 )
1885 table.adddata(
1886 "{}".format(exp_str(self._ellipticity[ifo])), dataclass=ifo
1887 )
1888
1889 # quadrupole
1890 if dist == None:
1891 self._q22[ifo] = None
1892 table.adddata("*", dataclass=ifo)
1893 else:
1894 self._q22[ifo] = h0_to_quadrupole(self._h0ul[ifo], freq, dist)
1895 table.adddata("{}".format(exp_str(self._q22[ifo])), dataclass=ifo)
1896
1897 # get spin down limit ratio
1898 if sdlim == None:
1899 self._sdratio[ifo] = None
1900 table.adddata("*", dataclass=ifo)
1901 else:
1902 self._sdratio[ifo] = self._h0ul[ifo] / sdlim
1903 table.adddata("%.02f" % self._sdratio[ifo], dataclass=ifo)
1904
1905 if len(self._harmonics) == 2 and self._modeltype == "source":
1906 if not self._biaxial:
1907 self._I21[ifo] = upper_limit_greedy(
1908 self._posteriors[ifo]["I21"].samples, upperlimit=(ul / 100.0)
1909 )
1910 table.adddata("{}".format(exp_str(self._I21[ifo])), dataclass=ifo)
1911
1912 self._I31[ifo] = upper_limit_greedy(
1913 self._posteriors[ifo]["I31"].samples, upperlimit=(ul / 100.0)
1914 )
1915 table.adddata("{}".format(exp_str(self._I31[ifo])), dataclass=ifo)
1916
1917 # add SNR
1918 if ifo in self._optimal_snrs:
1919 table.adddata("%.1f" % self._optimal_snrs[ifo], dataclass=ifo)
1920 else:
1921 table.adddata("*", dataclass=ifo)
1922
1923 # add max log likelihood values
1924 table.adddata(
1925 "%.f" % (self._maxL[ifo] / np.log(10.0)), dataclass=ifo
1926 ) # convert to log base 10
1927 table.adddata(
1928 "%.1f" % (self._Bsn[ifo] / np.log(10.0)), dataclass=ifo
1929 ) # convert to log base 10
1930
1931 # add coherent vs incoherent and coherent vs (incoherent or noise) for Joint ifo)
1932 if len(self._ifos) > 2 and "Joint" in self._ifos:
1933 if ifo == "Joint":
1934 table.adddata(
1935 "%.1f" % (self._Bci / np.log(10.0)), dataclass=ifo
1936 ) # convert to log base 10
1937 table.adddata(
1938 "%.1f" % (self._Bcin / np.log(10.0)), dataclass=ifo
1939 ) # convert to log base 10
1940 else:
1941 table.adddata("*", dataclass=ifo)
1942 table.adddata("*", dataclass=ifo)
1943
1944 self._limits_table += table.tabletext
1945 return self._limits_table
1946
1947 def credible_interval(self, ifo, param, ci=95, paramval=None):
1948 # get the given percentage (minimum) credible interval from samples for detector given by ifo and param
1949 # (for injections, where we have parameter values (paramval) get the corresponding smallest credible
1950 # region that contains the parameter value)
1951 samples = self._posteriors[ifo][param].samples.squeeze()
1952 samples.sort()
1953 lowbound, highbound, _ = self._ci_loop(samples, ci)
1954
1955 cival = None
1956 if paramval != None:
1957 cifound = False
1958 # loop over different credible intervals until finding the one that contains the injection
1959 for cit in range(1, 101):
1960 l, h, cr = self._ci_loop(samples, cit)
1961 if paramval >= l and paramval <= h:
1962 cifound = True
1963 cival = cit
1964 break
1965 if not cifound:
1966 cival = 100
1967
1968 return lowbound, highbound, cival
1969
1970 def _ci_loop(self, sortedsamples, ci):
1971 lowbound = sortedsamples[0]
1972 highbound = sortedsamples[-1]
1973 cregion = highbound - lowbound
1974 lsam = len(sortedsamples)
1975 cllen = int(lsam * float(ci) / 100.0)
1976 for j in range(lsam - cllen):
1977 if sortedsamples[j + cllen] - sortedsamples[j] < cregion:
1978 lowbound = sortedsamples[j]
1979 highbound = sortedsamples[j + cllen]
1980 cregion = highbound - lowbound
1981
1982 return lowbound, highbound, cregion
1983
1984
1985class create_background(posteriors):
1986 """
1987 Get information (evidence ratios and SNRs) from any the background analyses
1988 """
1989
1990 def __init__(
1991 self,
1992 backgrounddirs,
1993 snrs,
1994 Bsn,
1995 outputdir,
1996 Bci=None,
1997 Bcin=None,
1998 showcontours=True,
1999 ):
2000 # initialise with a dictionary (keyed to detectors) of directories containing the background analyses,
2001 # a dictionary of signal vs noise Bayes factors, a coherent vs incoherent Bayes factor and a coherent
2002 # vs incoherent or noise Bayes factor (all created from the posterior class)
2003 self._ifos_ifos = list(backgrounddirs.keys()) # detectors
2004 self._backgrounddirs = backgrounddirs
2005 self._dir_lists = (
2006 {}
2007 ) # dictionary with a list of background results directories for each detector
2008 self._Bci_fore = Bci # the "foreground" Bayes factor for coherent vs incoherent
2009 self._Bcin_fore = (
2010 Bcin # the "foreground" Bayes factor for coherent vs incoherent or noise
2011 )
2012
2013 self._Bsn_Bsn = (
2014 {}
2015 ) # dictionary of "background" signal vs noise Bayes factor lists for each detector
2016 self._Bci_Bci = [] # list of "background" coherent vs incoherent Bayes factors
2017 self._Bcin_Bcin = (
2018 []
2019 ) # list of "background" coherent vs incoherent or noise Bayes factors
2020 self._optimal_snrs_optimal_snrs = (
2021 {}
2022 ) # dictionary of "background" SNR lists for each detector
2023 self._signal_evidence_signal_evidence = {}
2024 self._noise_evidence_noise_evidence = {}
2025 self._maxL_maxL = {}
2026
2027 # figure and contour limits
2028 self._figlimits = [(0.0, None), ()] # assumes SNR is first value
2029 self._contourlimits = [(0.0, None), ()]
2030
2031 # set some default arguments for create_joint_posterior_plot
2032 self._plot_kwargs = {
2033 "ratio": 2,
2034 "whichtruth": "all",
2035 "scatter_kwargs": {"alpha": 0.5},
2036 "figlimits": self._figlimits,
2037 "contourlimits": self._contourlimits,
2038 }
2039
2040 self._Bsn_prob = (
2041 {}
2042 ) # probability of the "foreground value" given a Gaussian KDE applied to the background
2043 self._Bci_prob = None
2044 self._Bcin_prob = None
2045
2046 self._injection_parameters_injection_parameters = None # set to None
2047
2048 self._outputdir_outputdir = outputdir # set output directory
2049
2050 dirlen = []
2051 for i, ifo in enumerate(self._ifos_ifos):
2052 if ifo not in Bsn or ifo not in snrs:
2053 print(
2054 "Error... Bsn/SNR dictionary is not consistent with the background directories dictionary.",
2055 file=sys.stderr,
2056 )
2057 sys.exit(1)
2058
2059 self._Bsn_Bsn[ifo] = [] # create empty list
2060 self._optimal_snrs_optimal_snrs[ifo] = [] # create empty list
2061 self._signal_evidence_signal_evidence[ifo] = []
2062 self._noise_evidence_noise_evidence[ifo] = []
2063 self._maxL_maxL[ifo] = []
2064
2065 # get directory lists
2066 self._dir_lists[ifo] = [
2067 os.path.join(self._backgrounddirs[ifo], d)
2068 for d in os.listdir(self._backgrounddirs[ifo])
2069 if os.path.isdir(os.path.join(self._backgrounddirs[ifo], d))
2070 ]
2071 dirlen.append(len(self._dir_lists[ifo]))
2072 if dirlen[i] == 0:
2073 print(
2074 "Error... no background results directories were present for '%s'."
2075 % ifo,
2076 file=sys.stderr,
2077 )
2078 sys.exit(1)
2079
2080 # check there are the same number of background runs for each ifo
2081 if i > 0:
2082 if dirlen[i] != dirlen[0]:
2083 print(
2084 "Error... number of background results directories not then same for different detectors."
2085 % ifo,
2086 file=sys.stderr,
2087 )
2088 sys.exit(1)
2089
2090 self._Bsn_fore = Bsn # the "foreground" Bayes factor for signal vs noise
2091 self._snrs_fore = snrs # the "foreground" optimal snrs
2092
2093 # get SNRs
2094 self._get_snrs_get_snrs()
2095
2096 # get Bayes factors
2097 self._get_bayes_factors_get_bayes_factors()
2098
2099 # get probabilities the background distribution being greater than the foreground value
2100 self._get_bsn_prob()
2101 if len(self._ifos_ifos) > 2 and "Joint" in self._ifos_ifos:
2102 self._get_bci_prob()
2103
2104 @property
2105 def bci_prob(self):
2106 return self._Bci_prob
2107
2108 @property
2109 def bcin_prob(self):
2110 return self._Bcin_prob
2111
2112 def bsn_prob(self, ifo):
2113 return self._Bsn_prob[ifo]
2114
2115 def _get_snrs(self):
2116 # get the returned match filter SNRs
2117 for ifo in self._ifos_ifos: # loop over detectors
2118 # loop over background directories
2119 for pdir in self._dir_lists[ifo]:
2120 self._optimal_snrs_optimal_snrs[ifo].append(self.get_snr(pdir))
2121
2122 def _get_bayes_factors(self):
2123 # get the Bayes factors for the signal
2124 incoherent = np.zeros(
2125 len(self._dir_lists[self._ifos_ifos[0]])
2126 ) # the incoherent evidence
2127
2128 for ifo in self._ifos_ifos:
2129 for i, pdir in enumerate(self._dir_lists[ifo]):
2130 pfiles = os.listdir(pdir)
2131 Bsn = None
2132 for pfile in pfiles: # get HDF5 file with extenstion '.hdf' or '.h5'
2133 if fnmatch.fnmatch(
2134 pfile, "posterior_samples*.hdf"
2135 ) or fnmatch.fnmatch(pfile, "posterior_samples*.h5"):
2136 Bsn, sigev, noiseev, maxL = self.get_bayes_factor(
2137 os.path.join(pdir, pfile)
2138 )
2139
2140 self._Bsn_Bsn[ifo].append(Bsn)
2141 self._signal_evidence_signal_evidence[ifo].append(sigev)
2142 self._noise_evidence_noise_evidence[ifo].append(noiseev)
2143 self._maxL_maxL[ifo].append(maxL)
2144 break
2145 if Bsn == None:
2146 print(
2147 "Error... no HDF5 file with 'posterior_samples' in the name was found in '%s'."
2148 % pdir,
2149 file=sys.stderr,
2150 )
2151 sys.exit(1)
2152
2153 if ifo != "Joint":
2154 incoherent[i] += self._signal_evidence_signal_evidence[ifo][i]
2155
2156 # get the coherent vs incoherent noise evidence
2157 if len(self._ifos_ifos) > 2 and "Joint" in self._ifos_ifos:
2158 for i in range(len(self._dir_lists[self._ifos_ifos[0]])):
2159 self._Bci_Bci.append(self._signal_evidence_signal_evidence["Joint"][i] - incoherent[i])
2160 self._Bcin_Bcin.append(
2161 self._signal_evidence_signal_evidence["Joint"][i]
2162 - np.logaddexp(incoherent[i], self._noise_evidence_noise_evidence["Joint"][i])
2163 )
2164
2165 def _get_bsn_prob(self):
2166 # get the probability of the background distribution of signal vs noise Bayes factors
2167 # being greater than the foreground value for each detector by using a Gaussian KDE on the
2168 # 1D samples from the background
2169 for ifo in self._ifos_ifos:
2170 kernel = stats.gaussian_kde(self._Bsn_Bsn[ifo])
2171 self._Bsn_prob[ifo] = kernel.integrate_box_1d(
2172 self._Bsn_fore[ifo], np.inf
2173 ) # integrate between foreground value and infinity
2174
2175 def _get_bci_prob(self):
2176 kernel = stats.gaussian_kde(self._Bci_Bci)
2177 self._Bci_prob = kernel.integrate_box_1d(self._Bci_fore, np.inf)
2178 kernel = stats.gaussian_kde(self._Bcin_Bcin)
2179 self._Bcin_prob = kernel.integrate_box_1d(self._Bcin_fore, np.inf)
2180
2181 def bsn_plot(self, credint=[0.5, 0.95]):
2182 # create 2D scatter plots of SNR vs Bsn (with KDE estimated probability contours - if present just for Joint)
2183 # and 1D histograms for each - use the create_joint_posterior_plot's function
2184
2185 self._posteriors_posteriors = {}
2186 truths = {}
2187
2188 # create fake posterior objects (requires logL to be there) for each detector
2189 for ifo in self._ifos_ifos:
2190 self._posteriors_posteriors[ifo] = Posterior((["tmp", "logL"], np.zeros((100, 2))))
2191
2192 # add the SNR and Bsn values
2193 snrs = PosteriorOneDPDF("snr", np.array([self._optimal_snrs_optimal_snrs[ifo]]).T)
2194 self._posteriors_posteriors[ifo].append(snrs)
2195
2196 bsn = PosteriorOneDPDF(
2197 "bsn", np.array([self._Bsn_Bsn[ifo]]).T / np.log(10.0)
2198 ) # convert into base 10 log
2199 self._posteriors_posteriors[ifo].append(bsn)
2200
2201 # remove temporary varaibles
2202 self._posteriors_posteriors[ifo].pop("tmp")
2203 self._posteriors_posteriors[ifo].pop("logl")
2204
2205 self._parameters_parameters = ["snr", "bsn"] # set parameters
2206 truths[ifo] = [self._snrs_fore[ifo], self._Bsn_fore[ifo] / np.log(10.0)]
2207
2208 return self.create_joint_posterior_plot(
2209 ["snr", "bsn"],
2210 bins=int(np.log2(len(snrs.samples))),
2211 truths=truths,
2212 credintervals=credint,
2213 filename="bsn",
2214 **self._plot_kwargs,
2215 )
2216
2217 def bci_plot(self, credint=[0.5, 0.95], which="bci"):
2218 self._posteriors_posteriors = {}
2219 if which == "bcin":
2220 truths = [self._snrs_fore["Joint"], self._Bcin_fore / np.log(10.0)]
2221 self._parameters_parameters = ["snr", "bcin"]
2222 else: # default to bci
2223 truths = [self._snrs_fore["Joint"], self._Bci_fore / np.log(10.0)]
2224 self._parameters_parameters = ["snr", "bci"]
2225
2226 if (
2227 which != "bci"
2228 ): # force which to be 'bci' (as there are only the two options)
2229 which = "bci"
2230
2231 self._posteriors_posteriors["Joint"] = Posterior((["tmp", "logL"], np.zeros((100, 2))))
2232
2233 # add the SNR and Bci values
2234 snrs = PosteriorOneDPDF("snr", np.array([self._optimal_snrs_optimal_snrs["Joint"]]).T)
2235 self._posteriors_posteriors["Joint"].append(snrs)
2236
2237 if which == "bcin":
2238 bci = PosteriorOneDPDF("bcin", np.array([self._Bcin_Bcin]).T / np.log(10.0))
2239 else:
2240 bci = PosteriorOneDPDF("bci", np.array([self._Bci_Bci]).T / np.log(10.0))
2241 self._posteriors_posteriors["Joint"].append(bci)
2242
2243 curifos = self._ifos_ifos # save current detectors
2244 self._ifos_ifos = ["Joint"] # set for just the joint detector
2245 bciplot = self.create_joint_posterior_plot(
2246 self._parameters_parameters,
2247 bins=int(np.log2(len(snrs.samples))),
2248 truths=truths,
2249 credintervals=credint,
2250 filename=which,
2251 **self._plot_kwargs,
2252 )
2253 self._ifos_ifos = curifos # reset detectors (in case they are needed later)
2254
2255 return bciplot
2256
2257 def create_background_table(self):
2258 # create table with background plots in them
2259 background_plots_table = htmltag(
2260 "h2", "Evidence distributions", newline=True
2261 ).text
2262
2263 table = htmltable()
2264 table.addrow()
2265 cols = 1 # number of columns
2266
2267 # add Bsn plot
2268 bsnplot = self.bsn_plot()
2269 table.adddata(
2270 atag(
2271 os.path.basename(bsnplot[0]),
2272 '<img class="backgroundplot" src="{}"/>'.format(
2273 os.path.basename(bsnplot[0])
2274 ),
2275 ).text
2276 )
2277
2278 if self._Bci_fore != None:
2279 bciplot = self.bci_plot(which="bci")
2280 table.adddata(
2281 atag(
2282 os.path.basename(bciplot[0]),
2283 '<img class="backgroundplot" src="{}"/>'.format(
2284 os.path.basename(bciplot[0])
2285 ),
2286 ).text
2287 )
2288 cols += 1
2289
2290 if self._Bcin_fore != None:
2291 bcinplot = self.bci_plot(which="bcin")
2292 table.adddata(
2293 atag(
2294 os.path.basename(bcinplot[0]),
2295 '<img class="backgroundplot" src="{}"/>'.format(
2296 os.path.basename(bcinplot[0])
2297 ),
2298 ).text
2299 )
2300 cols += 1
2301
2302 # add probabilities for the background distribution being greater than the foreground value
2303 innertable = htmltable(tableclass="evidencetable")
2304 innertable.addrow()
2305 innertable.adddata("&nbsp;")
2306 innertable.adddata(
2307 "Probability of background distribution being greater than foreground",
2308 colspan=(len(self._ifos_ifos) + cols - 1),
2309 header=True,
2310 )
2311 innertable.addrow()
2312 innertable.adddata("Detector", header=True)
2313 for ifo in self._ifos_ifos:
2314 innertable.adddata(ifo, dataclass=ifo)
2315 if self._Bci_fore != None:
2316 innertable.adddata("B<sub>CvI</sub>")
2317 if self._Bcin_fore != None:
2318 innertable.adddata("B<sub>CvIN</sub>")
2319
2320 innertable.addrow()
2321 innertable.adddata("&nbsp;")
2322 for ifo in self._ifos_ifos:
2323 innertable.adddata(exp_str(self.bsn_prob(ifo)), dataclass=ifo)
2324
2325 if self._Bci_fore != None:
2326 innertable.adddata(exp_str(float(self.bci_prob), 1))
2327
2328 if self._Bcin_fore != None:
2329 innertable.adddata(exp_str(float(self.bcin_prob), 1))
2330
2331 table.addrow()
2332 table.adddata(innertable.tabletext, colspan=cols)
2333
2334 background_plots_table += table.tabletext
2335
2336 return background_plots_table
2337
2338
2339if __name__ == "__main__":
2340 description = """This script will create a results page for a single pulsar from the known pulsar analysis pipeline. A configuration (.ini) file is required."""
2341 epilog = """An example configuration file could contain the following:
2342
2343# a section for general analysis information
2344[general]
2345parfile = 'path_to_par_file' # the path to the pulsar parameter file
2346detectors = ['H1', 'L1'] # a list of the detectors to use
2347with_joint = True # a boolean stating whether to also add the joint multi-detector analysis
2348joint_only = False # a boolean stating whether to only output the joint multi-detector analysis
2349with_background = False # a boolean stating whether to include background analyses
2350injection = False # a boolean stating whether this pulsar is a software/hardware injection
2351upper_limit = 95 # the percentage credible upper limit value
2352credible_interval = [95] # a list of the percentage credible intervals for output statistics
2353use_gw_phase = False # a boolean stating whether to assume the initial phase parameter is the rotational, or gravitational wave (l=m=2), phase (e.g. if looking a hardware injections)
2354harmonics = [2] # a list of the frequency harmonics used in this analysis
2355model_type = waveform # either 'waveform' (default) or 'source' specify the parameterisation
2356biaxial = False # set whether the signal searched for was from a biaxial source
2357priorfile = 'path_to_prior' # the path to the prior file used in the analysis (if given priors will be plotted)
2358
2359# a section for parameter estimation files
2360[parameter_estimation]
2361posteriors = {'H1': 'path_to_H1_posteriors', 'L1': 'path_to_L1_posteriors', 'Joint': 'path_to_Joint_posteriors'} # a dictionary (keyed on detectors) pointing to the locations of posterior sample files
2362background = {'H1': 'path_to_H1_background_dir', 'L1': 'path_to_L1_background_dir', 'Joint': 'path_to_Joint_backgroud_dir'} # a dictionary (keyed on detectors) pointing to directories containing the background posterior files
2363
2364# a section for pre-processed data information
2365[data]
2366files = {'H1': 'path_to_H1_data', 'L1': 'path_to_L1_data'} # a dictionary (keyed on detectors) pointing to the locations (or lists of locations for multiple harmonics) of pre-processed data files
2367
2368# a section for the output location for this pulsar
2369[output]
2370path = 'path_to_output_base_directory' # the path to the base directory in which the results page will be created
2371indexpage = 'path_to_index_page' # an optional path (relative to the base directory) to the index page containing results from multiple pulsars
2372
2373# a section for plotting options
2374[plotting]
2375all_posteriors = False # a boolean stating whether to show joint posterior plots of all parameters (default: False)
2376subtract_truths = False # a boolean stating whether to subtract the true/heterodyned value from any phase parameters to centre the plot at zero for that value
2377show_contours = False # a boolean stating whether to show probabilty contours on 2D posterior plots (default: False)
2378eps_output = False # a boolean stating whether to also output eps versions of figures (png figures will automatically be produced)
2379pdf_output = False # a boolean stating whether to also output pdf versions of figures (png figures will automatically be produced)
2380
2381"""
2382
2383 parser = argparse.ArgumentParser(
2384 description=description,
2385 epilog=epilog,
2386 formatter_class=argparse.RawDescriptionHelpFormatter,
2387 )
2388 parser.add_argument("inifile", help="The configuration (.ini) file")
2389
2390 # parse input options
2391 opts = parser.parse_args()
2392
2393 inifile = opts.inifile
2394
2395 # open and parse config file
2396 cp = ConfigParser()
2397 try:
2398 cp.read(inifile)
2399 except:
2400 print(
2401 "Error... cannot parse configuration file '%s'" % inifile, file=sys.stderr
2402 )
2403 sys.exit(1)
2404
2405 # get the output directory
2406 try:
2407 outdir = cp.get("output", "path")
2408 except:
2409 print(
2410 "Error... no output directory 'path' specified in [output] section.",
2411 file=sys.stderr,
2412 )
2413 sys.exit(1)
2414
2415 # create directory if required
2416 if not os.access(outdir, os.W_OK) and not os.path.isdir(
2417 outdir
2418 ): # check if directory already exists
2419 try:
2420 os.makedirs(outdir)
2421 except:
2422 print(
2423 "Error... cannot make output directory '%s'." % outdir, file=sys.stderr
2424 )
2425 sys.exit(1)
2426
2427 # get the index page location
2428 try:
2429 indexpage = cp.get("output", "indexpage")
2430 except:
2431 indexpage = None
2432
2433 # get the pulsar parameter file
2434 parfile = None
2435 try:
2436 parfile = cp.get("general", "parfile")
2437 except:
2438 print(
2439 "Error... Pulsar parameter 'parfile' must specified in the [general] section.",
2440 file=sys.stderr,
2441 )
2442 sys.exit(1)
2443
2444 # read in data from par file
2445 try:
2446 par = psr_par(parfile)
2447 except:
2448 print(
2449 "Error... Pulsar parameter (.par) file '%s' could not be opened!" % parfile,
2450 file=sys.stderr,
2451 )
2452 sys.exit(1)
2453
2454 # get pulsar PSRJ name
2455 pname = par["PSRJ"]
2456 if not pname:
2457 print(
2458 "Error... no PSRJ name set in pulsar parameter file '%s'." % parfile,
2459 file=sys.stderr,
2460 )
2461
2462 # check/get required parameters
2463 f0 = par["F0"] # pulsar rotation frequency (Hz)
2464 if not f0:
2465 print(
2466 "Error... no 'F0' value in the parameter file '%s'." % parfile,
2467 file=sys.stderr,
2468 )
2469 sys.exit(1)
2470
2471 f1 = par["F1"] # pulsar first time derivative of rotational frequency (Hz/s)
2472 if not f1:
2473 f1 = 0.0 # set to zero if not given
2474
2475 # get the upper limit credible interval
2476 try:
2477 upperlim = ast.literval_eval(cp.get("general", "upper_limit"))
2478 except: # default to 95%
2479 upperlim = 95
2480
2481 # get credible intervals for output statistics
2482 try:
2483 credints = ast.literval_eval(cp.get("general", "credible_interval"))
2484 except: # default to 95%
2485 credints = [95]
2486
2487 # check whether looking at an injection or not
2488 try:
2489 injection = cp.getboolean("general", "injection")
2490 except:
2491 injection = (
2492 False # if nothing is given then assume that this is not an injection
2493 )
2494
2495 # check whether to use the rotational, or gravitational wave phase (e.g. for hardware injections), in output plots
2496 try:
2497 usegwphase = cp.getboolean("general", "use_gw_phase")
2498 except:
2499 usegwphase = False
2500
2501 # check whther a prior file is given
2502 try:
2503 priorfile = cp.get("general", "priorfile")
2504 except:
2505 priorfile = None
2506
2507 # attempt to get pulsar distance, proper motion corrected age and any association (e.g. GC from the ATNF catalogue)
2508 dist = p1_I = assoc = sdlim = f1sd = None
2509 atnfurl = None
2510 if not injection:
2511 # set ATNF URL where pulsar can be found
2512 atnfurl = (
2513 "http://www.atnf.csiro.au/people/pulsar/psrcat/proc_form.php?version="
2514 + ATNF_VERSION
2515 )
2516 atnfurl += "&startUserDefined=true&pulsar_names=" + re.sub(r"\+", "%2B", pname)
2517 atnfurl += "&ephemeris=long&submit_ephemeris=Get+Ephemeris&state=query"
2518
2519 # try getting information already parsed from ATNF catalogue by lalpulsar_knope pipeline setup
2520 jsonfile = os.path.join(outdir, pname + ".json")
2521 tryatnf = True
2522 if os.path.isfile(jsonfile):
2523 try:
2524 fp = open(jsonfile, "r")
2525 info = json.load(fp)
2526 fp.close()
2527
2528 # extract distance, intrinsic period derivative, pulsar association and required URL
2529 dist = info["Pulsar data"]["DIST"]
2530 p1_I = info["Pulsar data"]["P1_I"]
2531 assoc = info["Pulsar data"]["ASSOC"]
2532 tryatnf = False
2533 except:
2534 print(
2535 "Warning... could not read in JSON file '%s'." % jsonfile,
2536 file=sys.stderr,
2537 )
2538
2539 if tryatnf: # try getting ATNF info now
2540 pinfo = get_atnf_info(pname)
2541 if pinfo is not None:
2542 dist, p1_I, assoc, atnfurlref = pinfo # unpack values
2543
2544 # if distance is in the par file use that instead
2545 if par["DIST"]:
2546 dist = par["DIST"] / KPC # convert back into kpc
2547
2548 # set the corrected spin-down value (based on intrinsic period derivative or using conservative age (10^9) value from GC pulsars)
2549 f1sd = set_spin_down(p1_I, assoc, f0, f1)
2550
2551 # get spin-down limit
2552 if f1sd is not None and dist is not None:
2553 sdlim = spin_down_limit(f0, f1sd, dist)
2554
2555 # check whether to only include results from a joint (multi-detector) analysis
2556 try:
2557 jointonly = cp.getboolean("general", "joint_only")
2558 except:
2559 jointonly = False
2560
2561 # check which frequency harmonics were used in this analysis
2562 try:
2563 harmonics = ast.literal_eval(cp.get("general", "harmonics"))
2564 except:
2565 harmonics = [2] # default to twice the rotation frequency
2566
2567 # check whether the waveform or source model is being used
2568 try:
2569 modeltype = cp.get("general", "model_type")
2570 except:
2571 modeltype = "waveform"
2572
2573 if modeltype not in ["waveform", "source"]:
2574 print(
2575 "Error... unknown 'model type' '%s' specified." % modeltype, file=sys.stderr
2576 )
2577 sys.exit(1)
2578
2579 try:
2580 biaxial = cp.getboolean("general", "biaxial")
2581 except:
2582 biaxial = False
2583
2584 # check whether a background analysis has also been performed
2585 try:
2586 with_background = cp.getboolean("general", "with_background")
2587 except:
2588 with_background = False
2589
2590 if with_background:
2591 try:
2592 backgrounddir = ast.literal_eval(
2593 cp.get("parameter_estimation", "background")
2594 )
2595 except:
2596 with_background = False
2597
2598 # check whether to plot all joint posteriors
2599 try:
2600 allposteriors = cp.getboolean("plotting", "all_posteriors")
2601 except:
2602 allposteriors = False # default to False
2603
2604 # check whether to subtract true/heterodyned values from distributions
2605 try:
2606 subtracttruths = cp.getboolean("plotting", "subtract_truths")
2607 except:
2608 subtracttruths = False
2609
2610 # check whether to show probability contours on 2D posterior plots
2611 try:
2612 showcontours = cp.getboolean("plotting", "show_contours")
2613 except:
2614 showcontours = False
2615
2616 figformat = ["png"] # default to outputting png versions of figures
2617 # check whether to (also) output figures as eps
2618 try:
2619 witheps = cp.getboolean("plotting", "eps_output")
2620 except:
2621 witheps = False
2622
2623 if witheps:
2624 figformat.append("eps")
2625
2626 # check whether to (also) output figures as pdf
2627 try:
2628 withpdf = cp.getboolean("plotting", "pdf_output")
2629 except:
2630 withpdf = False
2631
2632 if withpdf:
2633 figformat.append("pdf")
2634
2635 # get the detectors to use
2636 ifos = [] # list of detectors to use
2637 withjoint = False
2638 preprocdat = None
2639 datatable = None
2640 if not jointonly:
2641 try:
2642 ifos = ast.literal_eval(cp.get("general", "detectors"))
2643 except:
2644 print(
2645 "Error... could not parse list of 'detectors' in the [general] section.",
2646 file=sys.stderr,
2647 )
2648 sys.exit(1)
2649
2650 if not isinstance(ifos, list):
2651 print(
2652 "Error... the 'detectors' value in the [general] section must be a list.",
2653 file=sys.stderr,
2654 )
2655 sys.exit(1)
2656
2657 if len(ifos) < 1:
2658 print(
2659 "Error... the 'detectors' value in the [general] section must contain at least one detector name.",
2660 file=sys.stderr,
2661 )
2662 sys.exit(1)
2663
2664 # check whether to include the joint (multi-detector) analysis
2665 if len(ifos) > 1: # only required if there's been more than one detector
2666 try:
2667 withjoint = cp.getboolean("general", "with_joint")
2668 except:
2669 withjoint = True # default to including the joint (multi-detector) analysis results
2670
2671 # get paths to pre-processed data files
2672 try:
2673 preprocdat = ast.literal_eval(cp.get("data", "files"))
2674 except:
2675 print(
2676 "Error... could not parse dictionary of 'files' in the [data] section.",
2677 file=sys.stderr,
2678 )
2679 sys.exit(1)
2680
2681 if not isinstance(preprocdat, dict):
2682 print(
2683 "Error... the 'files' value in the [data] section must be a dictionary.",
2684 file=sys.stderr,
2685 )
2686 sys.exit(1)
2687
2688 # check there is a value for each detector
2689 for ifo in ifos:
2690 if ifo not in preprocdat: # check if detector is in dictionary
2691 print(
2692 "Error... no pre-processed data file is given for '%s'." % ifo,
2693 file=sys.stderr,
2694 )
2695 sys.exit(1)
2696
2697 # make table with all stats and plots
2698 datatable = create_data_table(
2699 preprocdat, outdir, figformats=figformat, harmonics=harmonics
2700 )
2701
2702 # add 'Joint' on to the list of detectors if required
2703 if withjoint or jointonly:
2704 ifos.append("Joint")
2705
2706 # get the posterior sample files
2707 try:
2708 postfiles = ast.literal_eval(cp.get("parameter_estimation", "posteriors"))
2709 except:
2710 print(
2711 "Error... could not parse dictionary of 'posteriors' in the [parameter_estimation] section.",
2712 file=sys.stderr,
2713 )
2714 sys.exit(1)
2715
2716 if not isinstance(postfiles, dict):
2717 print(
2718 "Error... the 'posteriors' value in the [parameter_estimation] section must be a dictionary.",
2719 file=sys.stderr,
2720 )
2721 sys.exit(1)
2722
2723 for ifo in ifos:
2724 if ifo not in postfiles: # check if detector is in dictionary
2725 print(
2726 "Error... no posterior file is given for '%s'." % ifo, file=sys.stderr
2727 )
2728 sys.exit(1)
2729 else: # check file exists
2730 if not os.path.isfile(postfiles[ifo]):
2731 print(
2732 "Error... posterior file '%s' for '%s' does not exist."
2733 % (postfiles[ifo], ifo),
2734 file=sys.stderr,
2735 )
2736 sys.exit(1)
2737
2738 # create list of links to the various parts of the page
2739 linktable = htmltag("div", tagstyle="text-align: left; float: left")
2740 linkstable = htmltable()
2741 linkstable.addrow()
2742
2743 # dictionary for output html page
2744 htmlinput = {}
2745
2746 htmlinput["psrname"] = pname
2747 htmlinput["title"] = pname
2748 if injection:
2749 titlename = "INJ " + pname
2750 else:
2751 titlename = "PSR " + pname
2752 if atnfurl != None:
2753 htmlinput["h1title"] = atag(atnfurl, linktext=titlename).text
2754 else:
2755 htmlinput["h1title"] = titlename
2756
2757 # create table of pulsar info from par file
2758 psrtable = create_psr_table(par)
2759 htmlinput["pulsartable"] = psrtable
2760
2761 # get posterior class (containing samples, sample plots and posterior plots)
2762 postinfo = posteriors(
2763 postfiles,
2764 outdir,
2765 ifos=ifos,
2766 harmonics=harmonics,
2767 modeltype=modeltype,
2768 biaxial=biaxial,
2769 parfile=parfile,
2770 usegwphase=usegwphase,
2771 subtracttruths=subtracttruths,
2772 priorfile=priorfile,
2773 showcontours=showcontours,
2774 )
2775
2776 # create table of upper limits, SNR and evidence ratios
2777 htmlinput["limitstable"] = postinfo.create_limits_table(
2778 f0, sdlim=sdlim, dist=dist, ul=upperlim
2779 )
2780
2781 htmlinput["selectedposteriors"] = postinfo.create_joint_plots_table(
2782 title="Selected parameters"
2783 )
2784
2785 if allposteriors:
2786 htmlinput["jointposteriors"] = postinfo.create_joint_plots_table(allparams=True)
2787 linkstable.adddata(
2788 "Posteriors ("
2789 + atag("#selectedposteriors", linktext="Selected").text
2790 + ", "
2791 + atag("#jointposteriors", linktext="All").text
2792 + ")",
2793 dataclass="rightborder",
2794 )
2795 else:
2796 htmlinput["jointposteriors"] = ""
2797 linkstable.adddata(
2798 atag("#selectedposteriors", linktext="Posteriors").text,
2799 dataclass="rightborder",
2800 )
2801
2802 if with_background:
2803 bginfo = create_background(
2804 backgrounddir,
2805 postinfo.snrs,
2806 postinfo.bsn,
2807 outdir,
2808 Bci=postinfo.bci,
2809 Bcin=postinfo.bcin,
2810 )
2811
2812 # add background evidence plots
2813 htmlinput["evidenceplots"] = bginfo.create_background_table()
2814 linkstable.adddata(
2815 atag("#evidenceplots", linktext="Backgrounds").text, dataclass="rightborder"
2816 )
2817 else:
2818 htmlinput["evidenceplots"] = ""
2819
2820 # add data plots
2821 if datatable is not None:
2822 htmlinput["dataplots"] = str(datatable)
2823 linkstable.adddata(
2824 atag("#dataplots", linktext="Data Plots").text, dataclass="rightborder"
2825 )
2826
2827 # add posterior samples plots
2828 htmlinput["posteriorsamples"] = postinfo.create_sample_plot_table(
2829 figformats=figformat
2830 )
2831 linkstable.adddata(
2832 atag("#posteriorsamples", linktext="Samples").text, dataclass="rightborder"
2833 )
2834
2835 # add statistics
2836 htmlinput["posteriorstats"] = postinfo.create_stats_table(credints=credints)
2837 linkstable.adddata(atag("#posteriorstats", linktext="Statistics").text)
2838
2839 linktable.set_tagtext(linkstable.tabletext)
2840
2841 # set index page link
2842 if indexpage != None:
2843 indexlink = htmltag(
2844 "div", tagstyle="text-align: left; float: right; padding-right: 8px"
2845 )
2846 indexlink.set_tagtext(atag(indexpage, linktext="Index Page").text)
2847 indexlinktxt = indexlink.text
2848 else:
2849 indexlinktxt = ""
2850
2851 htmlinput["linkstable"] = linktable.text + indexlinktxt
2852
2853 # output CSS file
2854 cssfile = os.path.join(outdir, "resultspage.css")
2855 fp = open(cssfile, "w")
2856 fp.write(result_page_css)
2857 fp.close()
2858
2859 htmlinput["cssfile"] = os.path.basename(cssfile)
2860
2861 # get time/date for file creation
2862 now = datetime.datetime.now()
2863
2864 # add footer containing author, date and command lines used for page
2865 htmlinput["footer"] = "{} - {}<br><br>Command lines used:<br>{}<br>{}<br>".format(
2866 __author__, now.strftime("%a %d %b %Y"), " ".join(sys.argv), __version__
2867 )
2868
2869 # create page
2870 try:
2871 htmlfile = os.path.join(outdir, pname + ".html")
2872 fp = open(htmlfile, "w")
2873 fp.write(htmlpage.format(**htmlinput))
2874 fp.close()
2875 except:
2876 print("Error... there was a problem outputting the html page.", file=sys.stderr)
2877 sys.exit(1)
2878
2879 # output JSON file with results information
2880 info = {} # a dictionary with the analysis information
2881
2882 info["PSR"] = pname # set the pulsar name
2883
2884 # data about the pulsar
2885 info["Pulsar data"] = {}
2886 info["Pulsar data"]["F0"] = f0
2887 info["Pulsar data"]["F0ROT"] = f0
2888 info["Pulsar data"]["F0GW"] = 2.0 * f0 # assuming l=m=2 mode emission
2889 info["Pulsar data"]["F1"] = f1
2890 info["Pulsar data"]["F1ROT"] = f1
2891 info["Pulsar data"]["F1GW"] = 2.0 * f1 # assuming l=m=2 mode emission
2892 info["Pulsar data"]["F1SD"] = f1sd # acceleration corrected f1
2893 info["Pulsar data"][
2894 "P1_I"
2895 ] = p1_I # intrinsic period derivative (corrected for proper motion/GC acceleration effects)
2896 info["Pulsar data"]["DIST"] = dist
2897 info["Pulsar data"]["RA"] = par["RA_RAD"] # right ascension in radians
2898 info["Pulsar data"]["DEC"] = par["DEC_RAD"] # declination in radians
2899 info["Pulsar data"]["START"] = par["START"]
2900 info["Pulsar data"]["FINISH"] = par["FINISH"]
2901 info["Pulsar data"]["BINARY"] = par["BINARY"]
2902 info["Pulsar data"]["PEPOCH"] = par["PEPOCH"]
2903 info["Pulsar data"]["POSEPOCH"] = par["POSEPOCH"]
2904 info["Pulsar data"]["EPHEM"] = par["EPHEM"]
2905 info["Pulsar data"]["UNITS"] = par["UNITS"]
2906 info["Pulsar data"]["ASSOC"] = assoc # any association for the pulsar
2907 info["Pulsar data"]["spin-down limit"] = sdlim
2908 info["Pulsar data"]["par file"] = parfile
2909 info["Pulsar data"]["ATNF URL"] = atnfurl
2910
2911 # data about the upper limits
2912 for ifo in ifos:
2913 info[ifo] = {}
2914 info[ifo]["Upper limits"] = {}
2915 info[ifo]["Upper limits"]["credible region"] = upperlim
2916 info[ifo]["Upper limits"]["H0"] = postinfo.h0_ul(ifo)
2917 info[ifo]["Upper limits"]["ELL"] = postinfo.ellipticity_ul(ifo)
2918 info[ifo]["Upper limits"]["Q22"] = postinfo.q22_ul(ifo)
2919 info[ifo]["Upper limits"]["C21"] = postinfo.C21_ul(ifo)
2920 info[ifo]["Upper limits"]["C22"] = postinfo.C22_ul(ifo)
2921 info[ifo]["Upper limits"]["I21"] = postinfo.I21_ul(ifo)
2922 info[ifo]["Upper limits"]["I31"] = postinfo.I21_ul(ifo)
2923 info[ifo]["Upper limits"]["spin-down ratio"] = postinfo.sdlim_ratio(ifo)
2924
2925 # evidence ratios
2926 info[ifo]["Bayes factors"] = {}
2927 info[ifo]["Bayes factors"]["Signal vs Noise"] = postinfo.bsn[ifo]
2928
2929 if "Joint" == ifo:
2930 info[ifo]["Bayes factors"]["Coherent vs Incoherent"] = postinfo.bci
2931 info[ifo]["Bayes factors"][
2932 "Coherent vs Incoherent or Noise"
2933 ] = postinfo.bcin
2934
2935 info[ifo]["SNR"] = postinfo.snr(ifo)
2936
2937 # amplitude noise spectal densities
2938 if ifo != "Joint":
2939 info[ifo]["Amplitude spectral density"] = {}
2940
2941 if len(harmonics) == 1:
2942 info[ifo]["Amplitude spectral density"]["Spectrum"] = datatable.asds[
2943 harmonics[0]
2944 ][ifo].tolist()
2945 info[ifo]["Amplitude spectral density"]["Mean"] = np.mean(
2946 datatable.asds[harmonics[0]][ifo]
2947 )
2948 info[ifo]["Amplitude spectral density"]["Median"] = np.median(
2949 datatable.asds[harmonics[0]][ifo]
2950 )
2951 info[ifo]["Amplitude spectral density"]["Maximum"] = np.max(
2952 datatable.asds[harmonics[0]][ifo]
2953 )
2954 else:
2955 for h in harmonics:
2956 info[ifo]["Amplitude spectral density"][
2957 "%df harmonic" % int(h)
2958 ] = {}
2959 info[ifo]["Amplitude spectral density"]["%df harmonic" % int(h)][
2960 "Spectrum"
2961 ] = datatable.asds[h][ifo].tolist()
2962 info[ifo]["Amplitude spectral density"]["%df harmonic" % int(h)][
2963 "Mean"
2964 ] = np.mean(datatable.asds[h][ifo])
2965 info[ifo]["Amplitude spectral density"]["%df harmonic" % int(h)][
2966 "Median"
2967 ] = np.median(datatable.asds[h][ifo])
2968 info[ifo]["Amplitude spectral density"]["%df harmonic" % int(h)][
2969 "Maximum"
2970 ] = np.max(datatable.asds[h][ifo])
2971
2972 jsonfile = os.path.join(outdir, pname + ".json")
2973 fp = open(jsonfile, "w")
2974 json.dump(info, fp, indent=2)
2975 fp.close()
2976
2977 sys.exit(0)
lalpulsar.pulsarhtmlutils.htmltable
Class to make and return a html table.
Definition: pulsarhtmlutils.py:736
lalpulsar.pulsarhtmlutils.htmltag
A class to create a html tag.
Definition: pulsarhtmlutils.py:607
lalpulsar.pulsarpputils.psr_par
Definition: pulsarpputils.py:392
lalpulsar_knope_result_page.create_background
Get information (evidence ratios and SNRs) from any the background analyses.
Definition: lalpulsar_knope_result_page.py:1988
lalpulsar_knope_result_page.create_background._outputdir
_outputdir
Definition: lalpulsar_knope_result_page.py:2048
lalpulsar_knope_result_page.create_background._backgrounddirs
_backgrounddirs
Definition: lalpulsar_knope_result_page.py:2004
lalpulsar_knope_result_page.create_background._contourlimits
_contourlimits
Definition: lalpulsar_knope_result_page.py:2029
lalpulsar_knope_result_page.create_background._noise_evidence
_noise_evidence
Definition: lalpulsar_knope_result_page.py:2024
lalpulsar_knope_result_page.create_background._signal_evidence
_signal_evidence
Definition: lalpulsar_knope_result_page.py:2023
lalpulsar_knope_result_page.create_background._maxL
_maxL
Definition: lalpulsar_knope_result_page.py:2025
lalpulsar_knope_result_page.create_background.bci_prob
def bci_prob(self)
Definition: lalpulsar_knope_result_page.py:2105
lalpulsar_knope_result_page.create_background._Bsn_fore
_Bsn_fore
Definition: lalpulsar_knope_result_page.py:2090
lalpulsar_knope_result_page.create_background._get_bci_prob
def _get_bci_prob(self)
Definition: lalpulsar_knope_result_page.py:2175
lalpulsar_knope_result_page.create_background._plot_kwargs
_plot_kwargs
Definition: lalpulsar_knope_result_page.py:2032
lalpulsar_knope_result_page.create_background._Bsn_prob
_Bsn_prob
Definition: lalpulsar_knope_result_page.py:2040
lalpulsar_knope_result_page.create_background._dir_lists
_dir_lists
Definition: lalpulsar_knope_result_page.py:2005
lalpulsar_knope_result_page.create_background._Bci_prob
_Bci_prob
Definition: lalpulsar_knope_result_page.py:2043
lalpulsar_knope_result_page.create_background._Bcin_prob
_Bcin_prob
Definition: lalpulsar_knope_result_page.py:2044
lalpulsar_knope_result_page.create_background._parameters
_parameters
Definition: lalpulsar_knope_result_page.py:2205
lalpulsar_knope_result_page.create_background._figlimits
_figlimits
Definition: lalpulsar_knope_result_page.py:2028
lalpulsar_knope_result_page.create_background._Bci
_Bci
Definition: lalpulsar_knope_result_page.py:2016
lalpulsar_knope_result_page.create_background._posteriors
_posteriors
Definition: lalpulsar_knope_result_page.py:2185
lalpulsar_knope_result_page.create_background.bsn_plot
def bsn_plot(self, credint=[0.5, 0.95])
Definition: lalpulsar_knope_result_page.py:2181
lalpulsar_knope_result_page.create_background.create_background_table
def create_background_table(self)
Definition: lalpulsar_knope_result_page.py:2257
lalpulsar_knope_result_page.create_background._ifos
_ifos
Definition: lalpulsar_knope_result_page.py:2003
lalpulsar_knope_result_page.create_background._injection_parameters
_injection_parameters
Definition: lalpulsar_knope_result_page.py:2046
lalpulsar_knope_result_page.create_background._Bcin
_Bcin
Definition: lalpulsar_knope_result_page.py:2017
lalpulsar_knope_result_page.create_background._Bsn
_Bsn
Definition: lalpulsar_knope_result_page.py:2013
lalpulsar_knope_result_page.create_background.bci_plot
def bci_plot(self, credint=[0.5, 0.95], which="bci")
Definition: lalpulsar_knope_result_page.py:2217
lalpulsar_knope_result_page.create_background._snrs_fore
_snrs_fore
Definition: lalpulsar_knope_result_page.py:2091
lalpulsar_knope_result_page.create_background.__init__
def __init__(self, backgrounddirs, snrs, Bsn, outputdir, Bci=None, Bcin=None, showcontours=True)
Initialise with a dictionary keyed in detector names containing paths to the equivalent posterior sam...
Definition: lalpulsar_knope_result_page.py:1999
lalpulsar_knope_result_page.create_background._Bcin_fore
_Bcin_fore
Definition: lalpulsar_knope_result_page.py:2009
lalpulsar_knope_result_page.create_background._get_bsn_prob
def _get_bsn_prob(self)
Definition: lalpulsar_knope_result_page.py:2165
lalpulsar_knope_result_page.create_background._optimal_snrs
_optimal_snrs
Definition: lalpulsar_knope_result_page.py:2020
lalpulsar_knope_result_page.create_background._get_snrs
def _get_snrs(self)
Definition: lalpulsar_knope_result_page.py:2115
lalpulsar_knope_result_page.create_background._get_bayes_factors
def _get_bayes_factors(self)
Definition: lalpulsar_knope_result_page.py:2122
lalpulsar_knope_result_page.create_background._Bci_fore
_Bci_fore
Definition: lalpulsar_knope_result_page.py:2008
lalpulsar_knope_result_page.create_background.bsn_prob
def bsn_prob(self, ifo)
Definition: lalpulsar_knope_result_page.py:2112
lalpulsar_knope_result_page.create_background.bcin_prob
def bcin_prob(self)
Definition: lalpulsar_knope_result_page.py:2109
lalpulsar_knope_result_page.create_data_table
Definition: lalpulsar_knope_result_page.py:239
lalpulsar_knope_result_page.create_data_table._Bkplots
_Bkplots
Definition: lalpulsar_knope_result_page.py:343
lalpulsar_knope_result_page.create_data_table._ends
_ends
Definition: lalpulsar_knope_result_page.py:330
lalpulsar_knope_result_page.create_data_table._ifos
_ifos
Definition: lalpulsar_knope_result_page.py:260
lalpulsar_knope_result_page.create_data_table.create_plots
def create_plots(self)
Definition: lalpulsar_knope_result_page.py:340
lalpulsar_knope_result_page.create_data_table._datatable
_datatable
Definition: lalpulsar_knope_result_page.py:251
lalpulsar_knope_result_page.create_data_table._asds
_asds
Definition: lalpulsar_knope_result_page.py:326
lalpulsar_knope_result_page.create_data_table.fscanplots
def fscanplots(self)
Definition: lalpulsar_knope_result_page.py:411
lalpulsar_knope_result_page.create_data_table.__str__
def __str__(self)
Definition: lalpulsar_knope_result_page.py:501
lalpulsar_knope_result_page.create_data_table._nharmonics
_nharmonics
Definition: lalpulsar_knope_result_page.py:259
lalpulsar_knope_result_page.create_data_table._length
_length
Definition: lalpulsar_knope_result_page.py:331
lalpulsar_knope_result_page.create_data_table.analysis_stats_td
def analysis_stats_td(self, ifo)
Definition: lalpulsar_knope_result_page.py:444
lalpulsar_knope_result_page.create_data_table.Bkplots
def Bkplots(self)
Definition: lalpulsar_knope_result_page.py:401
lalpulsar_knope_result_page.create_data_table.get_stats
def get_stats(self)
Definition: lalpulsar_knope_result_page.py:420
lalpulsar_knope_result_page.create_data_table.asdplots
def asdplots(self)
Definition: lalpulsar_knope_result_page.py:406
lalpulsar_knope_result_page.create_data_table._starts
_starts
Definition: lalpulsar_knope_result_page.py:329
lalpulsar_knope_result_page.create_data_table._asdtime
_asdtime
Definition: lalpulsar_knope_result_page.py:325
lalpulsar_knope_result_page.create_data_table._duty_factor
_duty_factor
Definition: lalpulsar_knope_result_page.py:332
lalpulsar_knope_result_page.create_data_table.__init__
def __init__(self, datafiles, outputdir, figformats=["png"], asdtime=86400, harmonics=[2.0])
Initialise with a dictionary keyed in detector names containing paths to the equivalent pre-processed...
Definition: lalpulsar_knope_result_page.py:249
lalpulsar_knope_result_page.create_data_table._asdplots
_asdplots
Definition: lalpulsar_knope_result_page.py:344
lalpulsar_knope_result_page.create_data_table._outputdir
_outputdir
Definition: lalpulsar_knope_result_page.py:307
lalpulsar_knope_result_page.create_data_table._harmonics
_harmonics
Definition: lalpulsar_knope_result_page.py:256
lalpulsar_knope_result_page.create_data_table._figformats
_figformats
Definition: lalpulsar_knope_result_page.py:316
lalpulsar_knope_result_page.create_data_table._fscanplots
_fscanplots
Definition: lalpulsar_knope_result_page.py:345
lalpulsar_knope_result_page.create_data_table.asds
def asds(self)
Definition: lalpulsar_knope_result_page.py:416
lalpulsar_knope_result_page.create_data_table.create_table
def create_table(self)
Definition: lalpulsar_knope_result_page.py:463
lalpulsar_knope_result_page.create_data_table._datafiles
_datafiles
Definition: lalpulsar_knope_result_page.py:261
lalpulsar_knope_result_page.posteriors
Get sample posteriors and created a set of functions for outputting tables, plots and posterior stati...
Definition: lalpulsar_knope_result_page.py:508
lalpulsar_knope_result_page.posteriors.snr
def snr(self, ifo)
Definition: lalpulsar_knope_result_page.py:1060
lalpulsar_knope_result_page.posteriors.bcin
def bcin(self)
Definition: lalpulsar_knope_result_page.py:891
lalpulsar_knope_result_page.posteriors._prior_parameters
_prior_parameters
Definition: lalpulsar_knope_result_page.py:569
lalpulsar_knope_result_page.posteriors._posteriors
_posteriors
Definition: lalpulsar_knope_result_page.py:554
lalpulsar_knope_result_page.posteriors.__init__
def __init__(self, postfiles, outputdir, ifos=None, harmonics=[2], modeltype="waveform", biaxial=False, usegwphase=False, parfile=None, priorfile=None, subtracttruths=False, showcontours=False)
Initialise with a dictionary keyed in detector names containing paths to the equivalent posterior sam...
Definition: lalpulsar_knope_result_page.py:527
lalpulsar_knope_result_page.posteriors._ellipticity
_ellipticity
Definition: lalpulsar_knope_result_page.py:1791
lalpulsar_knope_result_page.posteriors._h0ul
_h0ul
Definition: lalpulsar_knope_result_page.py:1785
lalpulsar_knope_result_page.posteriors.bsn
def bsn(self)
Definition: lalpulsar_knope_result_page.py:876
lalpulsar_knope_result_page.posteriors.create_joint_plots_table
def create_joint_plots_table(self, allparams=False, title="Joint distributions")
Definition: lalpulsar_knope_result_page.py:1488
lalpulsar_knope_result_page.posteriors._postfiles
_postfiles
Definition: lalpulsar_knope_result_page.py:553
lalpulsar_knope_result_page.posteriors._sdratio
_sdratio
Definition: lalpulsar_knope_result_page.py:1790
lalpulsar_knope_result_page.posteriors._q22
_q22
Definition: lalpulsar_knope_result_page.py:1792
lalpulsar_knope_result_page.posteriors._modeltype
_modeltype
Definition: lalpulsar_knope_result_page.py:574
lalpulsar_knope_result_page.posteriors.C22_ul
def C22_ul(self, ifo)
Definition: lalpulsar_knope_result_page.py:923
lalpulsar_knope_result_page.posteriors._stats_section
_stats_section
Definition: lalpulsar_knope_result_page.py:1776
lalpulsar_knope_result_page.posteriors._limits_table
_limits_table
Definition: lalpulsar_knope_result_page.py:1781
lalpulsar_knope_result_page.posteriors.create_sample_plot_table
def create_sample_plot_table(self, figformats=["png"])
Definition: lalpulsar_knope_result_page.py:1613
lalpulsar_knope_result_page.posteriors._I21
_I21
Definition: lalpulsar_knope_result_page.py:1788
lalpulsar_knope_result_page.posteriors._samples_table
_samples_table
Definition: lalpulsar_knope_result_page.py:1658
lalpulsar_knope_result_page.posteriors.create_stats_table
def create_stats_table(self, credints=[95])
Definition: lalpulsar_knope_result_page.py:1661
lalpulsar_knope_result_page.posteriors._get_bayes_factors
def _get_bayes_factors(self)
Definition: lalpulsar_knope_result_page.py:988
lalpulsar_knope_result_page.posteriors._parameters
_parameters
Definition: lalpulsar_knope_result_page.py:565
lalpulsar_knope_result_page.posteriors._ci_loop
def _ci_loop(self, sortedsamples, ci)
Definition: lalpulsar_knope_result_page.py:1970
lalpulsar_knope_result_page.posteriors._priorfile
_priorfile
Definition: lalpulsar_knope_result_page.py:568
lalpulsar_knope_result_page.posteriors.get_snr
def get_snr(self, pdir)
Definition: lalpulsar_knope_result_page.py:958
lalpulsar_knope_result_page.posteriors.C21_ul
def C21_ul(self, ifo)
Definition: lalpulsar_knope_result_page.py:916
lalpulsar_knope_result_page.posteriors._noise_evidence
_noise_evidence
Definition: lalpulsar_knope_result_page.py:557
lalpulsar_knope_result_page.posteriors._I31
_I31
Definition: lalpulsar_knope_result_page.py:1789
lalpulsar_knope_result_page.posteriors._outputdir
_outputdir
Definition: lalpulsar_knope_result_page.py:528
lalpulsar_knope_result_page.posteriors._injection_parameters
_injection_parameters
Definition: lalpulsar_knope_result_page.py:567
lalpulsar_knope_result_page.posteriors.create_limits_table
def create_limits_table(self, freq, sdlim=None, dist=None, ul=95)
Definition: lalpulsar_knope_result_page.py:1779
lalpulsar_knope_result_page.posteriors.bci
def bci(self)
Definition: lalpulsar_knope_result_page.py:886
lalpulsar_knope_result_page.posteriors._maxL
_maxL
Definition: lalpulsar_knope_result_page.py:558
lalpulsar_knope_result_page.posteriors.plot_prior
def plot_prior(self, ax, param, prior, orientation="horizontal", truth=0.0, npoints=100)
Definition: lalpulsar_knope_result_page.py:1384
lalpulsar_knope_result_page.posteriors.ellipticity_ul
def ellipticity_ul(self, ifo)
Definition: lalpulsar_knope_result_page.py:902
lalpulsar_knope_result_page.posteriors.parameters
def parameters(self)
Definition: lalpulsar_knope_result_page.py:871
lalpulsar_knope_result_page.posteriors._optimal_snrs
_optimal_snrs
Definition: lalpulsar_knope_result_page.py:564
lalpulsar_knope_result_page.posteriors._showcontours
_showcontours
Definition: lalpulsar_knope_result_page.py:582
lalpulsar_knope_result_page.posteriors._parfile
_parfile
Definition: lalpulsar_knope_result_page.py:566
lalpulsar_knope_result_page.posteriors.create_joint_posterior_plot
def create_joint_posterior_plot(self, parameters, bins=20, ifo=None, truths=None, credintervals=[0.9], filename=None, figformats=["png"], ratio=3, figlimits=None, contourlimits=None, jointsamples=True, whichtruth=None, scatter_kwargs={})
Definition: lalpulsar_knope_result_page.py:1082
lalpulsar_knope_result_page.posteriors._Bcin
_Bcin
Definition: lalpulsar_knope_result_page.py:563
lalpulsar_knope_result_page.posteriors._Bsn
_Bsn
Definition: lalpulsar_knope_result_page.py:559
lalpulsar_knope_result_page.posteriors._get_snrs
def _get_snrs(self)
Definition: lalpulsar_knope_result_page.py:951
lalpulsar_knope_result_page.posteriors._C22
_C22
Definition: lalpulsar_knope_result_page.py:1787
lalpulsar_knope_result_page.posteriors.h0_ul
def h0_ul(self, ifo)
Definition: lalpulsar_knope_result_page.py:895
lalpulsar_knope_result_page.posteriors.sdlim_ratio
def sdlim_ratio(self, ifo)
Definition: lalpulsar_knope_result_page.py:944
lalpulsar_knope_result_page.posteriors._C21
_C21
Definition: lalpulsar_knope_result_page.py:1786
lalpulsar_knope_result_page.posteriors._subtract_truths
_subtract_truths
Definition: lalpulsar_knope_result_page.py:581
lalpulsar_knope_result_page.posteriors.q22_ul
def q22_ul(self, ifo)
Definition: lalpulsar_knope_result_page.py:909
lalpulsar_knope_result_page.posteriors.I21_ul
def I21_ul(self, ifo)
Definition: lalpulsar_knope_result_page.py:930
lalpulsar_knope_result_page.posteriors.I31_ul
def I31_ul(self, ifo)
Definition: lalpulsar_knope_result_page.py:937
lalpulsar_knope_result_page.posteriors.snrs
def snrs(self)
Definition: lalpulsar_knope_result_page.py:881
lalpulsar_knope_result_page.posteriors.get_bayes_factor
def get_bayes_factor(self, postfile)
Definition: lalpulsar_knope_result_page.py:1031
lalpulsar_knope_result_page.posteriors._signal_evidence
_signal_evidence
Definition: lalpulsar_knope_result_page.py:556
lalpulsar_knope_result_page.posteriors._injection_credible_regions
_injection_credible_regions
Definition: lalpulsar_knope_result_page.py:570
lalpulsar_knope_result_page.posteriors._ifos
_ifos
Definition: lalpulsar_knope_result_page.py:538
lalpulsar_knope_result_page.posteriors._posterior_stats
_posterior_stats
Definition: lalpulsar_knope_result_page.py:555
lalpulsar_knope_result_page.posteriors._usegwphase
_usegwphase
Definition: lalpulsar_knope_result_page.py:578
lalpulsar_knope_result_page.posteriors._biaxial
_biaxial
Definition: lalpulsar_knope_result_page.py:575
lalpulsar_knope_result_page.posteriors.credible_interval
def credible_interval(self, ifo, param, ci=95, paramval=None)
Definition: lalpulsar_knope_result_page.py:1947
lalpulsar_knope_result_page.posteriors._Bci
_Bci
Definition: lalpulsar_knope_result_page.py:562
lalpulsar_knope_result_page.posteriors._harmonics
_harmonics
Definition: lalpulsar_knope_result_page.py:573
lalpulsar_knope_result_page.posteriors._credint
_credint
Definition: lalpulsar_knope_result_page.py:1663
lalpulsar.pulsarhtmlutils.exp_str
def exp_str(f, p=1, otype="html")
Definition: pulsarhtmlutils.py:1002
lalpulsar.pulsarpputils.pulsar_nest_to_posterior
def pulsar_nest_to_posterior(postfile, nestedsamples=False, removeuntrig=True)
This function will import a posterior sample file created by lalapps_nest2pos (or a nested sample fil...
Definition: pulsarpputils.py:2830
lalpulsar.pulsarpputils.upper_limit_greedy
def upper_limit_greedy(pos, upperlimit=0.95, nbins=100)
Definition: pulsarpputils.py:959
lalpulsar.pulsarpputils.h0_to_quadrupole
def h0_to_quadrupole(h0, freq, dist)
Definition: pulsarpputils.py:739
lalpulsar.pulsarpputils.spin_down_limit
def spin_down_limit(freq, fdot, dist)
Definition: pulsarpputils.py:722
lalpulsar.pulsarpputils.plot_posterior_chain
def plot_posterior_chain(poslist, param, ifos, grr=None, withhist=0, mplparams=False)
Definition: pulsarpputils.py:984
lalpulsar.pulsarpputils.get_atnf_info
def get_atnf_info(psr)
Get the pulsar (psr) distance (DIST in kpc), proper motion corrected period derivative (P1_I) and any...
Definition: pulsarpputils.py:3568
lalpulsar.pulsarpputils.plot_Bks_ASDs
def plot_Bks_ASDs(Bkdata, delt=86400, plotpsds=True, plotfscan=False, removeoutlier=None, mplparams=False)
Definition: pulsarpputils.py:1565
lalpulsar.pulsarpputils.h0_to_ellipticity
def h0_to_ellipticity(h0, freq, dist)
Definition: pulsarpputils.py:732
lalpulsar_MakeSFTDAG.int
int
Definition: lalpulsar_MakeSFTDAG.py:413
lalpulsar_MakeSFTDAG.str
str
Definition: lalpulsar_MakeSFTDAG.py:411
lalpulsar_PiecewiseSearch.float
float
Definition: lalpulsar_PiecewiseSearch.py:122
lalpulsar_PiecewiseSearch.format
format
Definition: lalpulsar_PiecewiseSearch.py:409
lalpulsar_knope_result_page.create_psr_table
def create_psr_table(par)
Create a html table of some information from the pulsar parameter file.
Definition: lalpulsar_knope_result_page.py:207
lalpulsar_knope_result_page.set_spin_down
def set_spin_down(p1_I, assoc, f0, f1, n=5.0)
Set the spin-down of the source based on the intrinsic period derivative (p1_I) corrected for any pro...
Definition: lalpulsar_knope_result_page.py:191
lalpulsar_knope_result_page.file
file
Definition: lalpulsar_knope_result_page.py:69
test_flat_prior.range
range
Definition: test_flat_prior.py:105