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