Sky map overlap integral can now use MOC (multi-ordered) GW sky maps; fixes #19

This MR adds the ability to handle multi-ordered (MOC) sky maps, compared to be only able to use flat sky maps before. In addition, can handle single pixel secondary sky maps (as would be applicable for Swift, since the localization is on the order or less than a single pixel).

This MR enables the ability to handle the following combinations of sky maps:

1. MOC-MOC
2. MOC-flat
3. flat-flat (previous method)
4. MOC-ra/dec
5. flat-ra/dec

Bin scripts have been updated with these new options.

Appropriate unittests were added to verify these methods give identical results.

These changes were also tested and verified with Gracedb Events (e.g. https://gracedb-playground.ligo.org/superevents/S220317ch/view/).

bin/raven_skymap_overlap

@@ -32,24 +32,62 @@ __author__ = "Brandon Piotrzkowski <brandon.piotrzkowski@ligo.org>"
 # Global imports.
 import numpy as np
 from ligo.raven import search
-from ligo.skymap.io import fits
+from ligo.skymap.io import read_sky_map
 
-from optparse import Option, OptionParser
+import argparse
 
 # Command line options.
-parser = OptionParser(
-    description = __doc__,
-    usage = "%prog [options]")
-opts, args = parser.parse_args()
-print(args)
+# arguments
+parser = argparse.ArgumentParser(description='New Nemo User Account Creator')
+parser.add_argument('-i', '--inputfiles', nargs='+', dest="inputfiles", type=str, required=True, metavar=('se_skymap_filename', 'ext_skymap_filename'),
+                    help="Full path to input file of usernames, one username per line.")
+parser.add_argument('-m', '--se_moc', dest="se_moc", action="store_true", help="Assume the superevent sky map is multi-ordered (MOC).")
+parser.add_argument('-M', '--ext_moc', dest="ext_moc", action="store_true", help="Assume the external sky map is multi-ordered (MOC).")
+parser.add_argument('-r', '--se_ring', dest="se_ring", action="store_true", help="Assume the superevent map uses RING ordering rather than nested.")
+parser.add_argument('-R', '--ext_ring', dest="ext_ring", action="store_true", help="Assume the external sky map uses RING ordering rather than nested.")
+parser.add_argument('-c', '--ra_dec', nargs=2, dest="ra_dec", type=float, metavar=('ra', 'dec'),
+                    help="RA and dec of external sky map if a single pixel sky map.")
+args = parser.parse_args()
 
-# Check only two skymaps are given
-if len(args) != 2:
-    raise AssertionError('only provide 2 sky map filenames')
-   
-se_skymap, header = fits.read_sky_map(args[0], moc=False)
-ext_skymap, header = fits.read_sky_map(args[1], moc=False)
 
+# Check either one or two skymaps are given
+if not (0 < len(args.inputfiles) <= 2):
+    raise AssertionError('Please provide one or two sky map filenames')   
 
-result = search.skymap_overlap_integral(se_skymap, ext_skymap)
+if args.ra_dec:
+    ra, dec = args.ra_dec[0], args.ra_dec[1]
+else:
+    ra, dec = None, None
+
+try:
+    if args.ext_moc:
+        ext_table = read_sky_map(args.inputfiles[1], moc=args.ext_moc)
+        ext_skymap = ext_table['PROBDENSITY']
+        ext_uniq = ext_table['UNIQ']
+    else:
+        ext_skymap, header = read_sky_map(args.inputfiles[1], moc=args.ext_moc, nest=not args.ext_ring)
+        ext_uniq = []
+    if args.se_moc:
+        se_table = read_sky_map(args.inputfiles[0], moc=args.se_moc)
+        se_skymap = se_table['PROBDENSITY']
+        se_uniq = se_table['UNIQ']
+    else:
+        se_skymap, header = read_sky_map(args.inputfiles[0], moc=args.se_moc, nest=not args.se_ring)
+        se_uniq = []
+except IndexError:
+    if args.se_moc:
+        se_table = read_sky_map(args.inputfiles[0], moc=args.se_moc)
+        se_skymap = se_table['PROBDENSITY']
+        se_uniq = se_table['UNIQ']
+    else:
+        se_skymap, header = read_sky_map(args.inputfiles[0], moc=args.se_moc, nest=not args.se_ring)
+        se_uniq = []
+    ext_skymap = []
+    ext_uniq = []
+
+
+result = search.skymap_overlap_integral(se_skymap, ext_skymap,
+                                        se_skymap_uniq=se_uniq, ext_skymap_uniq=ext_uniq,
+                                        ra=ra, dec=dec,
+                                        se_nested=not args.se_ring, ext_nested=not args.ext_ring)
 print("Skymap overlap integral: {}".format(result))