.gitignore

@@ -8,6 +8,7 @@ GW*/output/*.fits
 GW*/output/*.txt
 GW*/output/*.png
 .ipynb_checkpoints
+*.json
 *.err
 *.png
 *.svg
@@ -15,3 +16,4 @@ GW*/output/*.png
 *.dag.*
 site
 env
+__pycache__

docs/index.md

@@ -68,12 +68,12 @@ These events *do* have SNR time series, so the test exercises BAYESTAR's code pa
 
 ### Results
 
-| Event     | Baseline                          | Latest                        | Comparison                        |
-| --------- | --------------------------------- | ----------------------------- | --------------------------------- |
-| GW170814  | ![baseline][GW170814/baseline]    | ![latest][GW170814/latest]    | ![latest][GW170814/compare]       |
-|           | ![baseline][GW170814/3d/baseline] | ![latest][GW170814/3d/latest] | ![latest][GW170814/3d/compare]    |
-| GW170817  | ![baseline][GW170817/baseline]    | ![latest][GW170817/latest]    | ![latest][GW170817/compare]       |
-|           | ![baseline][GW170817/3d/baseline] | ![latest][GW170817/3d/latest] | ![latest][GW170817/3d/compare]    |
+| Event         | Baseline                                              | Latest                                            | [JS Divergence](https://en.wikipedia.org/wiki/Jensen%E2%80%93Shannon_divergence)  | Image Comparison                                                  |
+| ------------- | ----------------------------------------------------- | ------------------------------------------------- | --------------------------------------------------------------------------------- | ----------------------------------------------------------------- |
+{% for event in 'GW170814 GW170817'.split() -%}
+| {{ event }}   | ![baseline]({{ event }}/baseline/bayestar.png)        | ![latest]({{ event }}/output/bayestar.png)        | {{ format_float_scientific(json(event + '/output/divergence.json')['2d'], 2) }}   | ![difference]({{ event }}/output/bayestar-failed-diff.png)        |
+|               | ![baseline]({{ event }}/baseline/bayestar.volume.png) | ![latest]({{ event }}/output/bayestar.volume.png) | {{ format_float_scientific(json(event + '/output/divergence.json')['3d'], 2) }}   | ![difference]({{ event }}/output/bayestar.volume-failed-diff.png) |
+{% endfor %}
 
 [f2y/15/searched_prob]:             first2years/comparison/2015/searched_prob.svg
 [f2y/16/searched_prob]:             first2years/comparison/2016/searched_prob.svg
@@ -96,16 +96,3 @@ These events *do* have SNR time series, so the test exercises BAYESTAR's code pa
 [f2y/kde/16/searched_area_hist]:    first2years/kde/comparison/2016/searched_area_hist.svg
 [f2y/kde/15/offset_hist]:           first2years/kde/comparison/2015/offset_hist.svg
 [f2y/kde/16/offset_hist]:           first2years/kde/comparison/2016/offset_hist.svg
-
-[GW170814/3d/baseline]:             GW170814/baseline/bayestar.volume.png
-[GW170814/3d/compare]:              GW170814/output/bayestar.volume-failed-diff.png
-[GW170814/3d/latest]:               GW170814/output/bayestar.volume.png
-[GW170814/baseline]:                GW170814/baseline/bayestar.png
-[GW170814/compare]:                 GW170814/output/bayestar-failed-diff.png
-[GW170814/latest]:                  GW170814/output/bayestar.png
-[GW170817/3d/baseline]:             GW170817/baseline/bayestar.volume.png
-[GW170817/3d/compare]:              GW170817/output/bayestar.volume-failed-diff.png
-[GW170817/3d/latest]:               GW170817/output/bayestar.volume.png
-[GW170817/baseline]:                GW170817/baseline/bayestar.png
-[GW170817/compare]:                 GW170817/output/bayestar-failed-diff.png
-[GW170817/latest]:                  GW170817/output/bayestar.png

htcondor/divergence.py

+import argparse
+import json
+
+from astropy_healpix import HEALPix
+from astropy.table import join
+from ligo.skymap.io import read_sky_map
+from ligo.skymap.healpix_tree import HEALPIX_MACHINE_ORDER, HEALPIX_MACHINE_NSIDE
+from ligo.skymap.moc import uniq2nest
+import numpy as np
+from scipy.stats import norm
+
+HPX = HEALPix(nside=HEALPIX_MACHINE_NSIDE, order='nested')
+
+
+def plog2p(p):
+    """Calculate p * log2(p), but return 0 if p is 0."""
+    result = np.empty_like(p)
+    zero = p == 0
+    result[zero] = 0
+    p_nonzero = p[~zero]
+    result[~zero] = p_nonzero * np.log2(p_nonzero)
+    return result
+
+
+def js_divergence(p, q, dx):
+    """Calculate the Jenson-Shannon divergence.
+
+    This is the symmetrized Kullback-Leibler (KL) divergence computed using
+    base 2 logarithms and divided by 2. It is symmetric and has a value between
+    0 (for identical distributions) and 1."""
+    integrand = plog2p(p) + plog2p(q) - 2 * plog2p(0.5 * (p + q))
+    return 0.5 * (integrand * dx).sum()
+
+
+parser = argparse.ArgumentParser()
+parser.add_argument('skymaps', nargs=2)
+args = parser.parse_args()
+
+# Load sky map, convert from UNIQ indexing to NESTED indexing at machine order
+tables = []
+for arg in args.skymaps:
+    table = read_sky_map(arg, moc=True)
+    order, ipix = uniq2nest(table.columns.pop('UNIQ'))
+    ipix <<= 2 * (HEALPIX_MACHINE_ORDER - order)
+    table['NESTED'] = ipix
+    tables.append(table)
+
+# Resample both sky maps to a common grid
+table = join(
+    *tables, keys='NESTED', join_type='outer', metadata_conflicts='silent')
+table.sort('NESTED')
+for colname, column in table.columns.items():
+    if colname != 'NESTED':
+        for slice in np.ma.clump_masked(column):
+            if slice.start > 0:
+                column[slice] = column[slice.start - 1]
+
+# Calculate area per pixel
+dA = np.diff(table['NESTED'], append=HPX.npix) * HPX.pixel_area
+
+# Calculate JS divergence for 2D sky map
+divergence_2d = js_divergence(
+    table['PROBDENSITY_1'], table['PROBDENSITY_2'], dA)
+
+# Calculate JS divergence for 3D sky map
+max_r = 6 * max(table.meta['distmean'] for table in tables)
+n_r = 1000
+d_r = max_r / n_r
+r = d_r * np.arange(1, n_r)
+p = norm(loc=table['DISTMU_1'], scale=table['DISTSIGMA_1']).pdf(r[:, np.newaxis]) * table['PROBDENSITY_1'] * table['DISTNORM_1']
+q = norm(loc=table['DISTMU_2'], scale=table['DISTSIGMA_2']).pdf(r[:, np.newaxis]) * table['PROBDENSITY_2'] * table['DISTNORM_2']
+# Calculate volume of each voxel, defined as the region within the
+# HEALPix pixel and contained within the two centric spherical shells
+# with radii (r - d_r / 2) and (r + d_r / 2).
+dV = (np.square(r[:, np.newaxis]) + np.square(d_r) / 12) * d_r * dA
+divergence_3d = js_divergence(p, q, dV)
+
+# Done
+print(json.dumps({'2d': divergence_2d, '3d': divergence_3d}))

htcondor/historical.sh

@@ -7,4 +7,5 @@ ligo-skymap-plot-volume baseline/*.fits* -o baseline/bayestar.volume.png
 ligo-skymap-plot-volume output/*.fits* -o output/bayestar.volume.png --align-to baseline/*.fits*
 python -c 'from matplotlib.testing.compare import compare_images; print(compare_images("baseline/bayestar.png", "output/bayestar.png", 0))'
 python -c 'from matplotlib.testing.compare import compare_images; print(compare_images("baseline/bayestar.volume.png", "output/bayestar.volume.png", 0))'
+python ../htcondor/divergence.py baseline/*.fits* output/*.fits* > output/divergence.json
 fitsheader output/*.fits* > output/bayestar.txt

main.py

+import json as _json
+
+import numpy as np
+
+
+def define_env(env):
+    env.macro(np.format_float_scientific)
+
+    @env.macro
+    def json(filename):
+        with open(filename) as f:
+            return _json.load(f)

mkdocs.yml

@@ -9,3 +9,4 @@ plugins:
         - "first2years/kde/input/*"
         - "first2years/kde/output/*"
         - "*.err"
+  - macros

requirements.txt

 ligo.skymap
 mkdocs
 mkdocs-exclude
+mkdocs-macros-plugin
 mkdocs-material