Adding 2D Gaussian filtering

README.md

@@ -62,6 +62,18 @@ features requests thought up.
 
 ### Update History
 
+##### 0.9.8
+* Adding 2D Gaussian smoothing for the contour surfaces.
+* Renaming ``contourf`` and ``contourf_alpha`` to ``shade`` and ``shade_alpha``.
+* Updating some of the example plots.
+
+##### 0.9.7
+* Updating package setup scripts.
+
+##### 0.9.6
+* Updating package setup scripts.
+
+
 ##### 0.9.5
 * Adding markdown paper.
 

chainconsumer/chain.py

@@ -101,7 +101,7 @@ def add_chain(self, chain, parameters=None, name=None, weights=None, posterior=N
         return self
 
     def configure_general(self, bins=None, flip=True, rainbow=None, colours=None,
-                          serif=True, plot_hists=True, max_ticks=5, kde=False, smooth=10):  # pragma: no cover
+                          serif=True, plot_hists=True, max_ticks=5, kde=False, smooth=3):  # pragma: no cover
         r""" Configure the general plotting parameters common across the bar
         and contour plots. If you do not call this explicitly, the :func:`plot`
         method will invoke this method automatically.
@@ -158,7 +158,7 @@ def configure_general(self, bins=None, flip=True, rainbow=None, colours=None,
         self.parameters_general["rainbow"] = rainbow
         self.parameters_general["plot_hists"] = plot_hists
         self.parameters_general["kde"] = kde
-        if not smooth:
+        if not smooth or kde:
             smooth = None
         self.parameters_general["smooth"] = smooth
         if colours is None:
@@ -169,8 +169,8 @@ def configure_general(self, bins=None, flip=True, rainbow=None, colours=None,
         self._configured_general = True
         return self
 
-    def configure_contour(self, sigmas=None, cloud=None, contourf=None,
-                          contourf_alpha=1.0):  # pragma: no cover
+    def configure_contour(self, sigmas=None, cloud=None, shade=None,
+                          shade_alpha=None):  # pragma: no cover
         """ Configure the default variables for the contour plots. If you do not call this
         explicitly, the :func:`plot` method will invoke this method automatically.
 
@@ -185,10 +185,10 @@ def configure_contour(self, sigmas=None, cloud=None, contourf=None,
             and [1, 2] for multiple chains.
         cloud : bool, optional
             If set, overrides the default behaviour and plots the cloud or not
-        contourf : bool, optional
+        shade : bool, optional
             If set, overrides the default behaviour and plots filled contours or not
-        contourf_alpha : float, optional
-            Filled contour alpha value override.
+        shade_alpha : float, optional
+            Filled contour alpha value override. Default is 1.0
         """
         num_chains = len(self.chains)
 
@@ -205,10 +205,15 @@ def configure_contour(self, sigmas=None, cloud=None, contourf=None,
             cloud = False
         self.parameters_contour["cloud"] = cloud
 
-        if contourf is None:
-            contourf = num_chains == 1
-        self.parameters_contour["contourf"] = contourf
-        self.parameters_contour["contourf_alpha"] = contourf_alpha
+        if shade is None:
+            shade = num_chains <= 2
+        self.parameters_contour["shade"] = shade
+        if shade_alpha is None:
+            if num_chains == 1:
+                shade_alpha = 1.0
+            else:
+                shade_alpha = np.sqrt(1 / num_chains)
+        self.parameters_contour["shade_alpha"] = shade_alpha
 
         self._configured_contour = True
 
@@ -223,13 +228,15 @@ def configure_bar(self, summary=None, shade=None):  # pragma: no cover
             Will not work if you have multiple chains
         shade : bool, optional
             If set to true, shades in confidence regions in under histogram. By default
-            this happens if you have a single chain, but is disabled if you are comparing
-            multiple chains.
+            this happens if you less than 3 chains, but is disabled if you are comparing
+            more chains.
         """
         if summary is not None:
             summary = summary and len(self.chains) == 1
         self.parameters_bar["summary"] = summary
-        self.parameters_bar["shade"] = shade if shade is not None else len(self.chains) == 1
+        if shade is None:
+            shade = len(self.chains) <= 2
+        self.parameters_bar["shade"] = shade
         self._configured_bar = True
         return self
 
@@ -753,14 +760,13 @@ def _plot_bars(self, ax, parameter, chain_row, weights, colour, bins=25, flip=Fa
 
         kde = self.parameters_general["kde"]
         smooth = self.parameters_general["smooth"]
-        do_smooth = smooth is not None and smooth > 0
-        if not do_smooth:
-            smooth = 1
-        bins = np.linspace(extents[0], extents[1], smooth * bins + 1)
+        bins, smooth = self._get_smoothed_bins(smooth, bins)
+
+        bins = np.linspace(extents[0], extents[1], bins)
         hist, edges = np.histogram(chain_row, bins=bins, normed=True, weights=weights)
         edge_center = 0.5 * (edges[:-1] + edges[1:])
-        if do_smooth:
-            hist = gaussian_filter(hist, int(smooth / 2), mode='constant')
+        if smooth:
+            hist = gaussian_filter(hist, smooth, mode='constant')
         if kde:
             assert np.all(weights == 1.0), "You can only use KDE if your weights are all one. " \
                                            "If you would like weights, please vote for this issue: " \
@@ -774,7 +780,7 @@ def _plot_bars(self, ax, parameter, chain_row, weights, colour, bins=25, flip=Fa
                 ax.plot(xs, pdf.evaluate(xs), color=colour)
             interpolator = pdf.evaluate
         else:
-            if do_smooth:
+            if smooth:
                 if flip:
                     ax.plot(hist, edge_center, color=colour)
                 else:
@@ -786,7 +792,8 @@ def _plot_bars(self, ax, parameter, chain_row, weights, colour, bins=25, flip=Fa
                     orientation = "vertical"
                 ax.hist(edge_center, weights=hist, bins=edges, histtype="step",
                         color=colour, orientation=orientation)
-            interpolator = interp1d(edge_center, hist, kind="nearest")
+            interp_type = "linear" if smooth else "nearest"
+            interpolator = interp1d(edge_center, hist, kind=interp_type)
 
         if self.parameters_bar["shade"] and fit_values is not None:
             lower = fit_values[0]
@@ -816,27 +823,26 @@ def _plot_contour(self, ax, x, y, w, px, py, colour, bins=25, truth=None):  # pr
 
         levels = 1.0 - np.exp(-0.5 * self.parameters_contour["sigmas"] ** 2)
         smooth = self.parameters_general["smooth"]
-        do_smooth = smooth is not None and smooth > 0
-        if not do_smooth:
-            smooth = 1
+        bins, smooth = self._get_smoothed_bins(smooth, bins, marginalsied=False)
+
         colours = self._scale_colours(colour, len(levels))
         colours2 = [self._scale_colour(colours[0], 0.7)] + \
                    [self._scale_colour(c, 0.8) for c in colours[:-1]]
 
-        hist, x_bins, y_bins = np.histogram2d(x, y, bins=(bins * smooth + 1), weights=w)
+        hist, x_bins, y_bins = np.histogram2d(x, y, bins=bins, weights=w)
         x_centers = 0.5 * (x_bins[:-1] + x_bins[1:])
         y_centers = 0.5 * (y_bins[:-1] + y_bins[1:])
-        if do_smooth:
-            hist = gaussian_filter(hist, int(smooth / 2), mode='constant')
+        if smooth:
+            hist = gaussian_filter(hist, smooth, mode='constant')
         hist[hist == 0] = 1E-16
         vals = self._convert_to_stdev(hist.T)
         if self.parameters_contour["cloud"]:
-            skip = max(1, x.size / 80000)
+            skip = max(1, int(x.size / 50000))
             ax.scatter(x[::skip], y[::skip], s=10, alpha=0.4, c=colours[1],
                        marker=".", edgecolors="none")
-        if self.parameters_contour["contourf"]:
+        if self.parameters_contour["shade"]:
             ax.contourf(x_centers, y_centers, vals, levels=levels, colors=colours,
-                        alpha=self.parameters_contour["contourf_alpha"])
+                        alpha=self.parameters_contour["shade_alpha"])
         ax.contour(x_centers, y_centers, vals, levels=levels, colors=colours2)
 
         if truth is not None:
@@ -984,20 +990,25 @@ def _convert_to_stdev(self, sigma):  # pragma: no cover
 
         return sigma_cumsum[i_unsort].reshape(shape)
 
+    def _get_smoothed_bins(self, smooth, bins, marginalsied=True):
+        if smooth is None or not smooth or smooth == 0:
+            return bins, 0
+        else:
+            return ((3 if marginalsied else 2) * smooth * bins), smooth
+
     def _get_parameter_summary(self, data, weights, parameter, chain_index, desired_area=0.6827):
         if not self._configured_general:
             self.configure_general()
 
         smooth = self.parameters_general["smooth"]
-        do_smooth = smooth is not None and smooth > 0
-        if not do_smooth:
-            smooth = 1
         bins = self.parameters_general['bins'][chain_index]
-        hist, edges = np.histogram(data, bins=(smooth * bins + 1), normed=True, weights=weights)
+        bins, smooth = self._get_smoothed_bins(smooth, bins)
+
+        hist, edges = np.histogram(data, bins=bins, normed=True, weights=weights)
         edge_centers = 0.5 * (edges[1:] + edges[:-1])
         xs = np.linspace(edge_centers[0], edge_centers[-1], 10000)
-        if do_smooth:
-            hist = gaussian_filter(hist, int(smooth / 2), mode='constant')
+        if smooth:
+            hist = gaussian_filter(hist, smooth, mode='constant')
 
         if self.parameters_general["kde"]:
             kde_xs = np.linspace(edge_centers[0], edge_centers[-1], max(100, int(bins)))

examples/customisations/plot_cloud_sigma.py

@@ -7,6 +7,10 @@
 
 In this example we display more sigma levels, turn on the point cloud, and
 disable the parameter summaries on the top of the marginalised distributions.
+
+Note that because of the very highly correlated distribution we have, it is
+useful to increase the number of bins the plots are generated with, to capture the
+thinness of the correlation.
 """
 
 import numpy as np
@@ -19,7 +23,7 @@
     data = multivariate_normal(normal(size=3), 0.5 * (cov + cov.T), size=100000)
 
     c = ChainConsumer().add_chain(data, parameters=["$x$", "$y$", "$z$"])
-    c.configure_bar(summary=False)
+    c.configure_bar(summary=False).configure_general(bins=1.4)
     c.configure_contour(cloud=True, sigmas=np.linspace(0, 2, 10))
     fig = c.plot()
 

examples/customisations/plot_colours_shade.py

@@ -5,7 +5,7 @@
 
 Choose custom colours and plot multiple chains with shading.
 
-Normally when plotting more than one chain, shading is removed so
+Normally when plotting more than two chains, shading is removed so
 you can clearly see the outlines. However, you can turn shading back
 on and modify the shade opacity if you prefer colourful plots.
 
@@ -22,13 +22,15 @@
 if __name__ == "__main__":
     np.random.seed(2)
     cov = normal(size=(2, 2)) + np.identity(2)
-    data = multivariate_normal(normal(size=2), 0.5 * (cov + cov.T), size=100000)
+    d1 = multivariate_normal(normal(size=2), 0.5 * (cov + cov.T), size=100000)
     cov = normal(size=(2, 2)) + np.identity(2)
-    data2 = multivariate_normal(normal(size=2), 0.5 * (cov + cov.T), size=100000)
+    d2 = multivariate_normal(normal(size=2), 0.5 * (cov + cov.T), size=100000)
+    cov = normal(size=(2, 2)) + np.identity(2)
+    d3 = multivariate_normal(normal(size=2), 0.5 * (cov + cov.T), size=100000)
 
-    c = ChainConsumer().add_chain(data, parameters=["$x$", "$y$"]).add_chain(data2)
-    c.configure_general(colours=["#B32222", "#D1D10D"])
-    c.configure_contour(contourf=True, contourf_alpha=0.5)
+    c = ChainConsumer().add_chain(d1, parameters=["$x$", "$y$"]).add_chain(d2).add_chain(d3)
+    c.configure_general(colours=["#B32222", "#D1D10D", "#455A64"])
+    c.configure_contour(shade=True, shade_alpha=0.2)
     c.configure_bar(shade=True)
     fig = c.plot()
 

examples/customisations/plot_dont_flip.py

@@ -16,7 +16,8 @@
 from chainconsumer import ChainConsumer
 
 if __name__ == "__main__":
-    data = np.random.multivariate_normal([0.0, 4.0], [[1.0, 0.7], [0.7, 1.5]], size=100000)
+    np.random.seed(0)
+    data = np.random.multivariate_normal([0.0, 4.0], [[1.0, 0.7], [0.7, 1.5]], size=1000000)
 
     c = ChainConsumer().add_chain(data, parameters=["$x_1$", "$x_2$"])
     c.configure_general(flip=False, max_ticks=10)

examples/customisations/plot_fewer_parameters.py

@@ -7,6 +7,9 @@
 
 Here we only plot the first four parameters. You could also simply pass the number four,
 which means the *first* four parameters.
+
+For fun, we also plot everything in green. Note you don't need to give multiple colours,
+the shading is all computed from the colour hex code.
 """
 
 import numpy as np
@@ -19,7 +22,7 @@
     cov = random(size=(6, 6))
     data = multivariate_normal(normal(size=6), 0.5 * (cov + cov.T), size=200000)
     parameters = ["$x$", "$y$", "$z$", "$a$", "$b$", "$c$"]
-    c = ChainConsumer().add_chain(data, parameters=parameters)
+    c = ChainConsumer().add_chain(data, parameters=parameters).configure_general(colours=["#388E3C"])
     fig = c.plot(parameters=parameters[:4], figsize="page")
 
     fig.set_size_inches(2.5 + fig.get_size_inches())  # Resize fig for doco. You don't need this.

examples/customisations/plot_no_smooth.py

@@ -5,7 +5,7 @@
 
 We can turn off the default gaussian filter on marginalised distributions.
 
-This can be done by setting ``smooth`` to either ``0`` or ``None``.
+This can be done by setting ``smooth`` to either ``0``, ``None`` or ``False``.
 Note that the parameter summaries also have smoothing turned off, and
 thus summaries may change.
 

examples/customisations/plot_rainbow_serif_bins.py

@@ -9,6 +9,12 @@
 enable it before that point if you wish. You can also pick how many bins you want to
 display your data with.
 
+You can see that in this example, we pick too many bins and would not get good
+summaries. If you simply want more (or less) bins than the default estimate,
+if you input a float instead of an integer, the number of bins will simply scale
+by that amount. For example, if the estimated picks 20 bins, and you set ``bins=1.5``
+your plots and summaries would be calculated with 30 bins.
+
 """
 import numpy as np
 from numpy.random import normal, random, multivariate_normal

examples/plot_introduction.py

@@ -8,13 +8,20 @@
 We give truth values, parameter labels and set the figure size to
 fit one column of a two column document.
 
+Note that the parameter summaries are all calculated from the chosen bin size and take
+into account if the data is being smoothed or not. It is thus important to consider
+whether you want smoothing enabled or (depending on your surfaces) more or less
+bins than automatically estimated. See the extended customisation examples for
+more information.
+
 """
 
 import numpy as np
 from chainconsumer import ChainConsumer
 
 if __name__ == "__main__":
-    data = np.random.multivariate_normal([0.0, 4.0], [[1.0, 0.7], [0.7, 1.5]], size=100000)
+    np.random.seed(0)
+    data = np.random.multivariate_normal([0.0, 4.0], [[1.0, 0.7], [0.7, 1.5]], size=1000000)
 
     c = ChainConsumer()
     c.add_chain(data, parameters=["$x_1$", "$x_2$"])

test_chain.py

@@ -28,7 +28,7 @@ def test_summary_no_smooth(self):
         tolerance = 5e-2
         consumer = ChainConsumer()
         consumer.add_chain(self.data)
-        consumer.configure_general(smooth=0)
+        consumer.configure_general(smooth=0, bins=2.4)
         summary = consumer.get_summary()
         actual = np.array(list(summary[0].values())[0])
         expected = np.array([3.5, 5.0, 6.5])