LALPulsar  6.1.0.1-89842e6
lalpulsar_knope_collate_results.py
Go to the documentation of this file.
1 # -*- coding: utf-8 -*-
2 #
3 # lalpulsar_knope_collate_results.py
4 #
5 # Copyright 2013, 2016
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 from __future__ import print_function, division
24 
25 # standard library imports
26 import argparse
27 import ast
28 import json
29 import datetime
30 import os
31 import re
32 import sys
33 from six.moves.configparser import ConfigParser
34 
35 # related third party imports
36 import numpy as np
37 
38 import matplotlib
39 
40 matplotlib.use("Agg")
41 
42 # local application/library specific imports
43 from lalpulsar import git_version
44 from lalpulsar.pulsarpputils import *
45 from lalpulsar.pulsarhtmlutils import *
46 
47 __author__ = "Matthew Pitkin <matthew.pitkin@ligo.org>"
48 __version__ = "git id {}".format(git_version.id)
49 __date__ = git_version.date
50 
51 
52 # create format for the output page
53 htmlpage = """
54 <!DOCTYPE html>
55 <html lang="en">
56 
57 <head>
58  <meta http-equiv="Content-Type" content="text/html; charset=utf-8"/>
59  <meta name="description" content="Known pulsar search results"/>
60  <meta charset="UTF-8">
61 
62  <title>Known pulsar search results</title>
63 
64  <link rel="stylesheet" type="text/css" href="{cssfile}"/>
65 </head>
66 
67 <body>
68 
69 <h1>Results table {nsources}</h1>
70 
71 <!-- The table of pulsars and results -->
72 <div class="tablediv">
73 {resultstable}
74 </div>
75 <br />
76 
77 <!-- Footnotes from the table -->
78 {footnotes}
79 <br />
80 
81 <!-- a footer -->
82 <div id="footer">
83 {footer}
84 </div>
85 
86 </body>
87 </html>
88 """
89 
90 
91 # main function
92 if __name__ == "__main__":
93  description = """This script will collate results pages from multiple pulsars into a signal table."""
94  epilog = """An example configuration file could contain the following:
95 
96 # a section containing output information
97 [output]
98 path = path_to_output_directory # the path to the output directory to contain the page [required]
99 
100 # a section containing input directory information
101 [input]
102 path = path_to_input_directory # the path to the directory in which the individual results page directories live [required]
103 
104 # a section containing general information for the output table
105 [general]
106 detectors = ['H1', 'L1', 'Joint'] # a list of detctors (including 'Joint') whose results will be output
107 sort_value = name # the pulsar parameter on which to sort the results (default: 'name') [Allowed: 'name', 'freq', 'ra', 'dec', 'sdlim', 'hul', 'ell', 'sdrat', 'dist']
108 sort_direction = ascending # sort in ascending or decending order (default: 'ascending')
109 results = ['h0ul', 'ell'] # a list of result parameters to output (default: ['h0ul'] - the h0 upper limit) [Allowed: 'h0ul', 'ell', 'sdrat', 'q22', 'bsn', 'bci', 'bcin', 'C21ul', 'C22ul']
110 parameters = ['f0rot', 'ra', 'dec'] # a list of pulsar parameters to output (default: ['f0rot'] - the pulsar's rotation frequency) [Allowed: 'f0rot', 'f0gw', 'f1rot', 'f1gw', 'sdlim', 'ra', 'dec', 'dist']
111 """
112 
113  parser = argparse.ArgumentParser(
114  description=description,
115  epilog=epilog,
116  formatter_class=argparse.RawDescriptionHelpFormatter,
117  )
118  parser.add_argument("inifile", help="The configuration (.ini) file")
119 
120  # parse input options
121  opts = parser.parse_args()
122 
123  inifile = opts.inifile
124 
125  cp = ConfigParser()
126  try:
127  cp.read(inifile)
128  except:
129  print(
130  "Error... cannot parse configuration file '%s'" % inifile, file=sys.stderr
131  )
132  sys.exit(-1)
133 
134  # get the output path
135  try:
136  outpath = cp.get("output", "path")
137  except:
138  print(
139  "Error... no output directory 'path' specified in [output] section.",
140  file=sys.stderr,
141  )
142  sys.exit(1)
143 
144  # try making directory
145  if not os.access(outpath, os.W_OK) and not os.path.isdir(
146  outpath
147  ): # check if directory already exists
148  try:
149  os.makedirs(outpath)
150  except:
151  print(
152  "Error... cannot make output directory '%s'." % outpath, file=sys.stderr
153  )
154  sys.exit(1)
155 
156  # get the directory containing all the individual pulsar result directories
157  try:
158  inpath = cp.get("input", "path")
159  except:
160  print(
161  "Error... no input directory 'path' specified in [input] section.",
162  file=sys.stderr,
163  )
164  sys.exit(1)
165 
166  if not os.path.isdir(inpath):
167  print("Error... input path '%s' does not exist." % inpath, file=sys.stderr)
168  sys.exit(1)
169 
170  # check how to sort the results (e.g. by name)
171  try:
172  sorttype = cp.get("general", "sort_value")
173  except:
174  sorttype = "name" # default to sort by pulsar name
175 
176  # sorting can be by name (equivalently right ascension), frequency, declination, spin-down limit, h0 upper limit, ellipticity upper limit, spin-down ratio or distance
177  if sorttype not in [
178  "name",
179  "freq",
180  "ra",
181  "dec",
182  "sdlim",
183  "h0ul",
184  "ell",
185  "sdrat",
186  "dist",
187  ]:
188  print(
189  "Error... sorting can only be on either 'name', 'freq', 'ra', 'dec', 'sdlim', 'hul', 'ell', 'sdrat' or 'dist'.",
190  file=sys.stderr,
191  )
192  sys.exit(1)
193 
194  # get the direction of the sorting ('ascending' or 'descending')
195  try:
196  sortdirection = cp.get("general", "sort_direction")
197  if sortdirection == "ascending":
198  reverse = False
199  elif sortdirection == "descending":
200  reverse = True
201  else: # default to ascending
202  reverse = False
203  except:
204  reverse = False # default to ascending
205 
206  # get list of detectors (and/or 'Joint') to use in the output table
207  try:
208  ifos = ast.literal_eval(cp.get("general", "detectors"))
209  except:
210  print("Error... could not parse list of detectors to include.", file=sys.stderr)
211  sys.exit(1)
212 
213  # get list of value to output (e.g. ['h0ul', 'ell'] for the h0 upper limit and ellipticity limit
214  try:
215  outputlims = [
216  ol.upper() for ol in ast.literal_eval(cp.get("general", "results"))
217  ] # covert to upper case
218  except: # as default output the h0 upper limit
219  outputlims = ["H0UL"]
220 
221  # get a list of the pulsar parameters to output
222  try:
223  outputvals = [
224  ov.upper() for ov in ast.literal_eval(cp.get("general", "parameters"))
225  ] # convert to upper case
226  except: # as default output the rotational frequency
227  outputvals = ["F0ROT"]
228 
229  # dictionary for output html and LaTeX pages
230  htmlinput = {}
231  latexinput = {}
232 
233  # get directories containing results pages
234  resdirs = os.listdir(inpath)
235  if len(resdirs) == 0:
236  print(
237  "Error... the input directory '%s' contains no other directories." % inpath,
238  file=sys.stderr,
239  )
240  sys.exit(1)
241 
242  resultsdata = {} # dictionary to hold results data
243  sourcedirs = [
244  os.path.join(inpath, rd)
245  for rd in resdirs
246  if os.path.isfile(os.path.join(inpath, rd, "{}.ini".format(rd)))
247  ]
248  totalsources = len(sourcedirs)
249  cursources = 0 # currently number of completed sources
250  for d in sourcedirs:
251  ld = os.listdir(d)
252  jsonfile = None
253  for fd in ld:
254  if ".json" == os.path.splitext(fd)[-1]:
255  jsonfile = os.path.join(d, fd)
256  break
257  if jsonfile == None: # no file found, so move on to next directory
258  continue
259 
260  # read in json file
261  try:
262  fp = open(jsonfile, "r")
263  pdata = json.load(fp)
264  fp.close()
265  except:
266  print(
267  "Warning... could not read JSON file '%s'. Skipping this directory."
268  % jsonfile
269  )
270  continue
271 
272  # check dictionary contains a 'Pulsar data' key
273  if "Pulsar data" not in pdata:
274  print(
275  "Warning... no 'Pulsar data' field in JSON file '%s'. Skipping this directory."
276  % jsonfile
277  )
278  continue
279 
280  # check dictionary contains 'PSR' key
281  if "PSR" not in pdata:
282  print(
283  "Warning... no 'PSR' pulsar name field in JSON file '%s'. Skipping this directory."
284  % jsonfile
285  )
286  continue
287 
288  # check that the request detectors are within the dictionary
289  ifopresent = True
290  for ifo in ifos:
291  if ifo not in pdata:
292  print(
293  "Warning... no key for detector '%s' in JSON file '%s'. Skipping this directory."
294  % (ifo, jsonfile)
295  )
296  ifopresent = False
297  break
298  if not ifopresent:
299  continue
300 
301  # add data into dictionary
302  resultsdata[pdata["PSR"]] = pdata
303  resultsdata[pdata["PSR"]]["path"] = os.path.relpath(
304  d, outpath
305  ) # relative path to result directory
306  cursources += 1
307 
308  if len(resultsdata) == 0:
309  print("Error... no reults were found!", file=sys.stderr)
310  sys.exit(1)
311 
312  # perform sorting
313  if "freq" == sorttype: # sort by frequency
314  sortlist = [
315  (resultsdata[pname]["Pulsar data"]["F0"], pname) for pname in resultsdata
316  ]
317  elif "dec" == sorttype: # sort by declination
318  sortlist = [
319  (resultsdata[pname]["Pulsar data"]["DEC"], pname) for pname in resultsdata
320  ]
321  elif "dist" == sorttype: # sort by distance
322  sortlist = [
323  (resultsdata[pname]["Pulsar data"]["DIST"], pname) for pname in resultsdata
324  ]
325  elif "sdlim" == sorttype: # sort by spin-down limit
326  sortlist = [
327  (resultsdata[pname]["Pulsar data"]["spin-down limit"], pname)
328  for pname in resultsdata
329  ]
330  elif (
331  "h0ul" == sorttype
332  ): # sort on h0 upper limit (get minimum h0 upper limit from all requested IFOs)
333  sortlist = [
334  (
335  min([resultsdata[pname][ifo]["Upper limits"]["H0"] for ifo in ifos]),
336  pname,
337  )
338  for pname in resultsdata
339  ]
340  elif (
341  "ell" == sorttype
342  ): # sort on ellipticity upper limit (get minimum upper limit from all requested IFOs)
343  sortlist = [
344  (
345  min([resultsdata[pname][ifo]["Upper limits"]["ELL"] for ifo in ifos]),
346  pname,
347  )
348  for pname in resultsdata
349  ]
350  elif (
351  "sdrat" == sorttype
352  ): # sort on ratio of result to spin-down limit (get minimum upper limit from all requested IFOs)
353  sortlist = [
354  (
355  min(
356  [
357  resultsdata[pname][ifo]["Upper limits"]["spin-down ratio"]
358  for ifo in ifos
359  ]
360  ),
361  pname,
362  )
363  for pname in resultsdata
364  ]
365  else: # sort by name (equivalent to sorting by RA) by default
366  sortlist = [(pname, pname) for pname in resultsdata]
367 
368  sortedlist = [
369  p[1] for p in sorted(sortlist, reverse=reverse)
370  ] # return a list of sorted names
371 
372  nifos = len(ifos) # number of IFOs (including "Joint")
373  numprepars = 1 + len(outputvals) # number of parameters in start of table
374  numlims = len(outputlims) # number of limit parameter to output
375 
376  # check if outputing coherent vs incoherent Bayes factors (Bci or Bcin) - these will only be used for if 'Joint' is present in the detectors
377  showbci = showbcin = 0
378  if "BCI" in outputlims:
379  if "Joint" in ifos:
380  showbci = 1
381  outputlims.pop(outputlims.index("BCI"))
382  numlims -= 1
383  if "BCIN" in outputlims:
384  if "Joint" in ifos:
385  showbcin = 1
386  outputlims.pop(outputlims.index("BCIN"))
387  numlims -= 1
388 
389  # start creating the results table
390  restable = htmltable()
391  ltable = latextable(ncolumns=(numprepars + numlims * nifos + showbci + showbcin))
392  ltable.set_caption("Limits on the gravitational wave amplitude for known pulsars")
393 
394  # first row
395  restable.addrow()
396  ltable.addrow(underline=True)
397  restable.adddata("", dataclass="bottomborder", header=True, colspan=numprepars)
398  ltable.adddata("~", multicolumn=numprepars)
399  for ifo in ifos:
400  ncolspan = numlims
401  if ifo == "Joint":
402  ncolspan += showbci + showbcin
403  restable.adddata(
404  ifo,
405  dataclass=ifo,
406  header=True,
407  datastyle="text-align:center; border-left:1px solid #000; border-bottom:1px solid #000",
408  colspan=ncolspan,
409  )
410  ltable.adddata(ifo, multicolumn=ncolspan)
411 
412  # second row
413  restable.addrow()
414  ltable.addrow(underline=True)
415  restable.adddata("Pulsar", dataclass="bottomborder", header=True)
416  ltable.adddata("Pulsar")
417  for prepar in outputvals:
418  restable.adddata(paramhtmldict[prepar], dataclass="bottomborder", header=True)
419  ltable.adddata(paramlatexdict[prepar])
420 
421  for ifo in ifos:
422  dataclass = "leftborder" # class for first value
423  datastyle = "border-bottom:1px solid #000" # style for first value
424  for ol in outputlims:
425  if ol[-2:] == "UL": # if value is an upper limit add credible region
426  cr = list(resultsdata.values())[0][ifo]["Upper limits"][
427  "credible region"
428  ]
429  restable.adddata(
430  paramhtmldict[ol].format(cr),
431  dataclass=dataclass,
432  datastyle=datastyle,
433  header=True,
434  )
435  ltable.adddata(paramlatexdict[ol].format(cr))
436  else:
437  restable.adddata(
438  paramhtmldict[ol],
439  dataclass=dataclass,
440  datastyle=datastyle,
441  header=True,
442  )
443  ltable.adddata(paramlatexdict[ol])
444  dataclass = "bottomborder"
445  datastyle = ""
446  if (
447  ifo == "Joint"
448  ): # check whether to show additional coherent vs incoherent Bayes factor values
449  if showbci:
450  restable.adddata(
451  paramhtmldict["BCI"],
452  dataclass=dataclass,
453  datastyle=datastyle,
454  header=True,
455  )
456  ltable.adddata(paramlatexdict["BCI"])
457  if showbcin:
458  restable.adddata(
459  paramhtmldict["BCIN"],
460  dataclass=dataclass,
461  datastyle=datastyle,
462  header=True,
463  )
464  ltable.adddata(paramlatexdict["BCIN"])
465 
466  # a dictionary to convert between parameter names
467  convdict = {
468  "RA": "RA",
469  "DEC": "DEC",
470  "H0UL": "H0",
471  "C21UL": "C21",
472  "C22UL": "C22",
473  "I21UL": "I21",
474  "I31UL": "I31",
475  "SDLIM": "spin-down limit",
476  "SDRAT": "spin-down ratio",
477  "BSN": "Signal vs Noise",
478  "BCI": "Coherent vs Incoherent",
479  "BCIN": "Coherent vs Incoherent or Noise",
480  }
481 
482  # checks for whegther footnotes are required
483  dagger = False
484  ddagger = False
485 
486  # loop through par files to produce plots
487  for pname in sortedlist:
488  restable.addrow()
489  ltable.addrow()
490 
491  pulsar = resultsdata[pname] # the dictionary of information for this pulsar
492 
493  # output pulsar names and pulsar parameters (include link to results page)
494  restable.adddata(
495  atag(os.path.join(pulsar["path"], pname + ".html"), linktext=pname).text
496  )
497  ltable.adddata(
498  re.sub(r"\-", "\\\\textminus", pname)
499  ) # substitute - ('hyphen') for \textminus in names
500 
501  for prepar in outputvals:
502  htmlsdtag = ""
503  latexsdtag = ""
504 
505  # set footnotes for pulsar's that have had a "corrected" spin-down
506  if "SDLIM" in prepar:
507  if (
508  pulsar["Pulsar data"]["P1_I"] != None
509  ): # spin-down has been corrected for intrinsic motion effects
510  if (
511  pulsar["Pulsar data"]["P1_I"] > 0.0
512  ): # double check that intrinsic period derivative is positive
513  htmlsdtag = htmltag("sup", tagtext="&dagger;").text
514  latexsdtag = r"$\dagger$"
515  dagger = True
516  elif pulsar["Pulsar data"]["ASSOC"] != None:
517  if "GC" in pulsar["Pulsar data"]["ASSOC"]:
518  htmlsdtag = htmltag("sup", tagtext="&Dagger;").text
519  latexsdtag = r"$\ddagger$"
520  ddagger = True
521 
522  if prepar in convdict:
523  pn = convdict[prepar]
524  else:
525  pn = prepar
526  if pn in pulsar["Pulsar data"]:
527  preval = pulsar["Pulsar data"][pn]
528  if preval == None:
529  prevalhtml = "*"
530  prevallatex = "*"
531  else:
532  disphtmlfunc = paramhtmldispfunc.__dict__[prepar]
533  displatexfunc = paramlatexdispfunc.__dict__[prepar]
534  prevalhtml = disphtmlfunc(preval)
535  prevallatex = displatexfunc(preval)
536  else:
537  prevalhtml = "*"
538  prevallatex = "*"
539 
540  restable.adddata(prevalhtml + htmlsdtag)
541  ltable.adddata(prevallatex + latexsdtag)
542 
543  for ifo in ifos:
544  dataclass = "leftborder"
545 
546  for limpar in outputlims:
547  if limpar in convdict:
548  ln = convdict[limpar]
549  else:
550  ln = limpar
551 
552  section = "Upper limits" # generally values will be in the 'Upper limits' section
553  if (
554  limpar == "BSN"
555  ): # if getting the Bayes factor look in 'Bayes factors' section
556  section = "Bayes factors"
557 
558  if ln in pulsar[ifo][section]:
559  limval = pulsar[ifo][section][ln]
560 
561  if limval == None:
562  limvalhtml = "*"
563  limvallatex = "*"
564  else:
565  if limpar == "BSN": # convert to log base 10
566  limval = limval / np.log(10.0)
567  disphtmlfunc = paramhtmldispfunc.__dict__[limpar]
568  displatexfunc = paramlatexdispfunc.__dict__[limpar]
569  limvalhtml = disphtmlfunc(limval)
570  limvallatex = displatexfunc(limval)
571  else:
572  limvalhtml = "*"
573  limvallatex = "*"
574 
575  restable.adddata(limvalhtml, dataclass=dataclass)
576  dataclass = ""
577  ltable.adddata(limvallatex)
578 
579  # check if needing to output coherent vs incoherent Bayes factors
580  if ifo == "Joint":
581  for bu, show in [("BCI", showbci), ("BCIN", showbcin)]:
582  if show:
583  bn = convdict[bu]
584  if bn in pulsar[ifo]["Bayes factors"]:
585  bval = pulsar[ifo]["Bayes factors"][bn]
586  if bval == None:
587  bvalhtml = "*"
588  bvallatex = "*"
589  else:
590  bval = bval / np.log(10.0)
591  disphtmlfunc = paramhtmldispfunc.__dict__[bu]
592  displatexfunc = paramlatexdispfunc.__dict__[bu]
593  bvalhtml = disphtmlfunc(bval)
594  bvallatex = displatexfunc(bval)
595  else:
596  bvalhtml = "*"
597  bvallatex = "*"
598 
599  restable.adddata(bvalhtml, dataclass=dataclass)
600  dataclass = ""
601  ltable.adddata(bvallatex)
602 
603  htmlinput["nsources"] = "(%d/%d sources)" % (cursources, totalsources)
604  htmlinput["resultstable"] = restable.tabletext
605 
606  # add footnotes
607  if dagger or ddagger:
608  htmlfootnotes = htmltag("div", tagclass="footnotes")
609  if dagger:
610  htmlfootnotes += htmltag("sup", tagtext="&dagger;").text
611  htmlfootnotes += " The spin-down limit was calculated using a spin-down "
612  htmlfootnotes += atag(
613  "http://www.atnf.csiro.au/research/pulsar/psrcat/psrcat_help.html?type=normal&highlight=p1_i#p1_i",
614  linktext="corrected",
615  ).text
616  htmlfootnotes += " for proper motion effects.\n"
617  htmlfootnotes += "<br />\n"
618  ltable.set_postamble(
619  "$\\dagger$ The pulsar's spin-down is corrected for proper motion effects. \\\\\n"
620  )
621  if ddagger:
622  htmlfootnotes += htmltag("sup", tagtext="&ddagger;").text
623  htmlfootnotes += " The spin-down limit was calculated using a characteristic spin-down age of 10<sup>9</sup> years.\n"
624  ltable.set_postamble(
625  ltable.postamble
626  + "$\\ddagger$ The pulsar's spin-down is calculated using a characteristic spin-down age of $10^9$ years.\n"
627  )
628  htmlinput["footnotes"] = htmlfootnotes.text
629  else:
630  htmlinput["footnotes"] = ""
631 
632  # create CSS
633  cssfile = os.path.join(outpath, "table.css")
634  fp = open(cssfile, "w")
635  fp.write(results_table_css)
636  fp.close()
637 
638  htmlinput["cssfile"] = os.path.basename(cssfile)
639 
640  # get time/date for file creation
641  now = datetime.datetime.now()
642 
643  # add footer containing author, date and command lines used for page
644  htmlinput["footer"] = "{} - {}<br><br>Command lines used:<br>{}<br>{}<br>".format(
645  __author__, now.strftime("%a %d %b %Y"), " ".join(sys.argv), __version__
646  )
647 
648  # create page
649  try:
650  htmlfile = os.path.join(outpath, "index.html")
651  fp = open(htmlfile, "w")
652  fp.write(htmlpage.format(**htmlinput))
653  fp.close()
654  except:
655  print("Error... there was a problem outputting the html page.", file=sys.stderr)
656  sys.exit(1)
657 
658  # output the LaTeX table
659  try:
660  latexfile = os.path.join(outpath, "resultstable.tex")
661  fp = open(latexfile, "w")
662  fp.write(ltable.tabletext)
663  fp.close()
664  except:
665  print(
666  "Error... there was a problem outputting the LaTeX table.", file=sys.stderr
667  )
668  sys.exit(1)
min
#define min(a, b)
Definition: CrossCorrToplist.c:71
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.pulsarhtmlutils.latextable
Class to make a return a LaTeX table.
Definition: pulsarhtmlutils.py:837
lalpulsar_knope_collate_results.displatexfunc
displatexfunc
Definition: lalpulsar_knope_collate_results.py:533
lalpulsar_knope_collate_results.disphtmlfunc
disphtmlfunc
Definition: lalpulsar_knope_collate_results.py:532