Use stable logsumexp in grid plots

cosmosis/plotting/grid_plots.py

@@ -23,13 +23,21 @@
 from .plotter import Plotter
 import itertools
 import scipy.optimize
+from scipy.special import logsumexp
 try:
 	from cosmosis import output as output_module
 except ImportError:
 	print("Running without cosmosis: no pretty section names or running on ini files")
 
 
 class GridPlotter(Plotter):
+	@staticmethod
+	def _logsumexp(values):
+		"""Stable log(sum(exp(values))) for grid marginalization."""
+		if np.size(values) == 0:
+			return -np.inf
+		return logsumexp(values)
+
 	#def __init__(self, *args, **kwargs):
 		# super(GridPlotter, self).__init__(*args, **kwargs)
 		# convert the loaded chain data sets into grids
@@ -69,7 +77,7 @@ def _plot_1d(self, name1):
 		#marginalize
 		for k,v1 in enumerate(vals1):
 			w = np.where(cols1==v1)
-			like_sum[k] = np.log(np.exp(like[w]).sum())
+			like_sum[k] = self._logsumexp(like[w])
 		like = like_sum.flatten()
 
 		#linearly interpolate
@@ -149,7 +157,7 @@ def _plot_2d(self, name1, name2, log_like=True):
 		for k,(v1, v2) in enumerate(itertools.product(vals1, vals2)):
 			w = np.where((cols1==v1)&(cols2==v2))
 			i,j = np.unravel_index(k, like_sum.shape)
-			like_sum[i,j] = np.log(np.exp(like[w]).sum())
+			like_sum[i,j] = self._logsumexp(like[w])
 		like = like_sum.flatten()
 
 		#Normalize the log-likelihood to peak=0
@@ -210,4 +218,3 @@ def toggle_edge(edges, x1,y1,x2,y2):
 		sm = pylab.cm.ScalarMappable(cmap=colormap, norm=norm)
 		sm._A = [] #hack from StackOverflow to make this work
 		pylab.colorbar(sm, label='Likelihood')
-

cosmosis/test/test_grid_plots.py

@@ -0,0 +1,16 @@
+import numpy as np
+
+from cosmosis.plotting.grid_plots import GridPlotter
+
+
+def test_grid_logsumexp_handles_extreme_log_likelihoods():
+    values = np.array([1000.0, 999.0, -1000.0])
+    result = GridPlotter._logsumexp(values)
+    expected = 1000.0 + np.log1p(np.exp(-1.0))
+    np.testing.assert_allclose(result, expected)
+    assert np.isfinite(result)
+
+
+def test_grid_logsumexp_empty_group_is_negative_infinity():
+    result = GridPlotter._logsumexp(np.array([]))
+    assert result == -np.inf