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