Repository: cta-observatory/ctapipe
Pull request #309 (2017-01-20 12:05:35)

Cleanup  utils.Histogram
- Histogram.hist is now Histogram.data (less confusing)
- examples in documentation
- Histogram.from_fits() now works correctly again and is a static member function (can be used as a constructor)
- removed fromFits option in main constructor


Patch:
diff --git a/ctapipe/reco/shower_max.py b/ctapipe/reco/shower_max.py
index dfda9c2d..02972f66 100644
--- a/ctapipe/reco/shower_max.py
+++ b/ctapipe/reco/shower_max.py
@@ -35,7 +35,7 @@ class ShowerMaxEstimator:
             thickness.append(float(line.split()[col_thickness]))
 
         self.atmosphere = Histogram(axisNames=["altitude"])
-        self.atmosphere.hist = thickness * u.g * u.cm ** -2
+        self.atmosphere.data = thickness * u.g * u.cm ** -2
         self.atmosphere.bin_lower_edges = [np.array(altitude) * u.km]
 
     def interpolate(self, arg, outlierValue=0., order=3):
@@ -52,11 +52,11 @@ class ShowerMaxEstimator:
         coordinate = (argv - bin_u_edge) / (bin_u_edge - bin_l_edge) \
             * (bin_u - bin_l) + bin_u
 
-        return ndimage.map_coordinates(self.atmosphere.hist,
+        return ndimage.map_coordinates(self.atmosphere.data,
                                        [[coordinate]],
                                        order=order,
                                        cval=outlierValue)[0]\
-            * self.atmosphere.hist.unit
+            * self.atmosphere.data.unit
 
     def find_shower_max_height(self, energy, h_first_int, gamma_alt):
         """ 
@@ -96,7 +96,7 @@ class ShowerMaxEstimator:
         t_shower_max = t_first_int + c
 
         # now find the height with the wanted thickness
-        for ii, thick1 in enumerate(self.atmosphere.hist):
+        for ii, thick1 in enumerate(self.atmosphere.data):
             if t_shower_max > thick1:
                 height1 = self.atmosphere.bin_lower_edges[0][ii]
                 height2 = self.atmosphere.bin_lower_edges[0][ii - 1]
diff --git a/ctapipe/utils/fitshistogram.py b/ctapipe/utils/fitshistogram.py
index dbc6c753..a1ccf030 100644
--- a/ctapipe/utils/fitshistogram.py
+++ b/ctapipe/utils/fitshistogram.py
@@ -22,41 +22,64 @@ class Histogram:
 
     Parameters
     ----------
-    bins: array_like(int)
+    nbins: array_like(int)
         list of number of bins for each dimension
     ranges: list(tuple)
         list of (min,max) values for each dimension
     name: str
         name of histogram (will be used as FITS extension name when written
         to a file
-    axisNames: list(str)
+    axis_names: list(str)
         name of each axis
-    initFromFITS: None or `astropy.io.fits.ImageHDU`
-        if not None, construct the Histogram object from the given
-        FITS HDU
+
+    Examples
+    --------
+
+    >>> hist = Histogram(nbins=(10,10), ranges=[[-1,1], [-1,1]])
+    >>> data = np.random.normal(shape=(2*100)) # make 100 random 2D events
+    >>> hist.fill(data)
+
+
+    Get a point in the histogram (can also get multiple values at once by
+    passing an array)
+
+    >>> val = hist.get_value([0.1,-0.5])
+    >>> vals = hist.get_value([[0.1,-0.5], [0.9,0.9]])
+
+    Get the full data array and do things with it:
+
+    >>> meanx = hist.data.mean(axis=0)
+
+    Write it to a FITS image file:
+
+    >>> hist.to_fits().writeto("output.fits")
+
+    Read it from FITS image file:
+
+    >>> hist2 = Histogram.from_fits("output.fits")
 
     """
 
     def __init__(self, nbins=None, ranges=None, name="Histogram",
-                 axisNames=None, initFromFITS=None):
-        """ Initialize an unfilled histogram (need to call either
-        fill() or specify initFromFITS to put data into it)
+                 axis_names=None):
+        """ Initialize an unfilled histogram (need to call fill()  put data into it)
+
+        see also
+        --------
+        The `Histogram.from_fits()` constructor
         """
 
-        self.hist = np.zeros(nbins)
-        self._binLowerEdges = None  # TODO: should be a property, get only
+        self.data = np.zeros(nbins)
+        self._binLowerEdges = None
         self._nbins = np.array([nbins]).flatten()
         self._ranges = np.array(ranges, ndmin=2)
-        self.valueScale = None
-        self.valueZero = None
+        self.value_scale = None
+        self.value_zero = None
         self.name = name
         self._ctypes = None
-        self.axisNames = axisNames
+        self.axis_names = axis_names
         self._numsamples = 0
 
-        if (initFromFITS):
-            self.from_fits(initFromFITS)
-
         # sanity check on inputs:
 
         if self.ndims < 1:
@@ -65,17 +88,17 @@ class Histogram:
             raise ValueError("Dimensions of ranges {0} don't match bins {1}"
                              .format(len(self._ranges), self.ndims))
 
-        if self.axisNames is not None:  # ensure the array is size ndims
-            self.axisNames = np.array(self.axisNames)
-            self.axisNames.resize(self.ndims)
+        if self.axis_names is not None:  # ensure the array is size ndims
+            self.axis_names = np.array(self.axis_names)
+            self.axis_names.resize(self.ndims)
         else:
-            self.axisNames = ["axis{}".format(x) for x in range(self.ndims)]
+            self.axis_names = ["axis{}".format(x) for x in range(self.ndims)]
         
     def __str__(self,):
         return ("Histogram(name='{name}', axes={axnames}, "
                 "nbins={nbins}, ranges={ranges})"
                 .format(name=self.name, ranges=self._ranges,
-                        nbins=self._nbins, axnames=self.axisNames))
+                        nbins=self._nbins, axnames=self.axis_names))
     
     @property
     def bin_lower_edges(self):
@@ -107,9 +130,9 @@ class Histogram:
         """
         returns the number of outlier points (the number of input
         datapoints - the sum of the histogram). This assumes the data
-        of the histogram is unmodified (and is "counts")
+        of the histogram is unmodified (and represents "counts").
         """
-        return self._numsamples - self.hist.sum()
+        return self._numsamples - self.data.sum()
 
     def fill(self, datapoints, **kwargs):
         """
@@ -126,11 +149,10 @@ class Histogram:
             creating the histogram
         """
 
-        hist, binLowerEdges = np.histogramdd(datapoints,
-                                             bins=self._nbins,
-                                             range=self._ranges, **kwargs)
+        hist, __ = np.histogramdd(datapoints, bins=self._nbins,
+                                  range=self._ranges, **kwargs)
 
-        self.hist += hist
+        self.data += hist
         self._numsamples += len(datapoints)
 
     def bin_centers(self, index):
@@ -149,7 +171,7 @@ class Histogram:
         myhist.to_fits().writeto("outputfile.fits")
 
         """
-        ohdu = fits.ImageHDU(data=self.hist.transpose())
+        ohdu = fits.ImageHDU(data=self.data.transpose())
         ohdu.name = self.name
         ndim = len(self._nbins)
 
@@ -173,30 +195,31 @@ class Histogram:
             bin0pix = 0.5  # lower-left corner of first bin
             bin0coord = self._ranges[dim][0]
 
-            name = self.axisNames[dim]
+            name = self.axis_names[dim]
 
             ctype = name[0:4] + "-   "
-            if (self._ctypes != None):
+            if self._ctypes is not None:
                 ctype = self._ctypes[dim]
 
-            ohdu.header.update("CTYPE%d" % (dim + 1), ctype, name)
-            ohdu.header.update("CDELT%d" % (dim + 1), delta)
-            ohdu.header.update("CRVAL%d" % (dim + 1), bin0coord)
-            ohdu.header.update("CRPIX%d" % (dim + 1), bin0pix)
-            ohdu.header.update("CNAME%d" % (dim + 1), self.axisNames[dim])
+            ohdu.header.set("CTYPE%d" % (dim + 1), ctype, name)
+            ohdu.header.set("CDELT%d" % (dim + 1), delta)
+            ohdu.header.set("CRVAL%d" % (dim + 1), bin0coord)
+            ohdu.header.set("CRPIX%d" % (dim + 1), bin0pix)
+            ohdu.header.set("CNAME%d" % (dim + 1), self.axis_names[dim])
 
-        if (self.valueScale):
-            ohdu.header.update("BSCALE", 1.0 / float(self.valueScale))
+        if self.value_scale:
+            ohdu.header.set("BSCALE", 1.0 / float(self.value_scale))
 
-        if (self.valueZero):
-            ohdu.header.update("BZERO", float(self.valueZero))
+        if self.value_zero:
+            ohdu.header.set("BZERO", float(self.value_zero))
 
-        ohdu.header.update("NSAMP", self._numsamples,
-                           "Number of samples originally filled")
+        ohdu.header.set("NSAMP", self._numsamples, "Number of samples "
+                                                   "originally filled")
 
         return ohdu
 
-    def from_fits(self, inputFITS):
+    @staticmethod
+    def from_fits(input_fits):
         """
         load a FITS histogram file
 
@@ -204,62 +227,67 @@ class Histogram:
         - `inputfits`: filename or fits HDU
         """
 
-        if (type(inputFITS) == str):
-            hdu = fits.open(inputFITS)[1]
+        hist = Histogram()
+
+        if type(input_fits) == str:
+            hdu = fits.open(input_fits)[1]
         else:
-            hdu = inputFITS
+            hdu = input_fits
         
-        self.hist = hdu.data.transpose()
-        self._nbins = self.hist.shape
+        hist.data = hdu.data.transpose()
+        hist._nbins = hist.data.shape
 
         wcs = WCS(hdu.header)
-        ndim = len(self._nbins)
+        ndim = len(hist._nbins)
 
         edges = []
-        self._ranges = []
-        axisnames = []
-        self._ctypes = []
+        hist._ranges = []
+        axis_names = []
+        hist._ctypes = []
 
         for dim in range(ndim):
             # note that histogramdd returns edges for 0-N+1 (including
             # the lower edge of the non-existant next bin), so we do
             # the same here to keep things the same
-            ax = np.zeros((self._nbins[dim] + 1, ndim))
-            ax[:, dim] = np.arange(self._nbins[dim] + 1) + 0.5
+            ax = np.zeros((hist._nbins[dim] + 1, ndim))
+            ax[:, dim] = np.arange(hist._nbins[dim] + 1) + 0.5
             edges.append(wcs.wcs_pix2world(ax, 1)[:, dim])
-            self._ranges.append((edges[dim][0], edges[dim][-1]))
-            self._ctypes.append(hdu.header["CTYPE%d" % (dim + 1)])
-            if (hdu.header.get("CNAME%d" % (dim + 1))):
-                axisnames.append(hdu.header["CNAME%d" % (dim + 1)])
+            hist._ranges.append((edges[dim][0], edges[dim][-1]))
+            hist._ctypes.append(hdu.header["CTYPE%d" % (dim + 1)])
+            if hdu.header.get("CNAME%d" % (dim + 1)):
+                axis_names.append(hdu.header["CNAME%d" % (dim + 1)])
             else:
-                axisnames.append(self._ctypes[dim][0:4])
+                axis_names.append(hist._ctypes[dim][0:4])
 
-        self.axisNames = np.array(axisnames)
+        hist.axis_names = np.array(axis_names)
 
-        self._binLowerEdges = edges
-        self._ranges = np.array(self._ranges)
+        hist._binLowerEdges = edges
+        hist._ranges = np.array(hist._ranges)
 
-        if(hdu.header.get("BSCALE")):
-            self.valueScale = hdu.header["BSCALE"]
+        if hdu.header.get("BSCALE"):
+            hist.value_scale = hdu.header["BSCALE"]
 
-        if(hdu.header.get("BSCALE")):
-            self.valueZero = hdu.header["BZERO"]
+        if hdu.header.get("BSCALE"):
+            hist.value_zero = hdu.header["BZERO"]
 
-        if(hdu.header.get("NSAMP")):
-            self._numsamples += int(hdu.header["NSAMP"])
+        if hdu.header.get("NSAMP"):
+            hist._numsamples += int(hdu.header["NSAMP"])
 
-    def get_value(self, coords, outlierValue=None):
+        return hist
+
+    def get_value(self, coords, outlier_value=None):
         """ Returns the values of the histogram at the given world
         coordinate(s) 
 
-        Arguments:
-
-        - `coords`: list/array of M coordinates of dimension N (where
-          the N is the histogram dimension)
+        Parameters
+        ----------
 
-        - `outlierValue`: value for outliers, if None, coordinates
-          outside the edges of the histogram will be given the edge
-          value
+        coords: array_like
+            array of M coordinates of dimension N (where
+            the N is the histogram dimension)
+        outlier_value: float or None
+          value for outliers, if None, coordinates outside the edges of the
+          histogram will be given the edge value
 
         """
         # if (self._binLowerEdges == None):
@@ -275,28 +303,26 @@ class Histogram:
                                      self.bin_lower_edges[ii][1:])
                          for ii in range(ndims)])
 
-        maxbin = np.array(self.hist.shape)
+        maxbin = np.array(self.data.shape)
 
         # deal with out-of-range values:
-        if (outlierValue == None):
+        if outlier_value is None:
             # extrapolate (simply for now, just takes edge value)
             bins[bins < 0] = 0
             for ii in range(ndims):
                 bins[ii][bins[ii] >= maxbin[ii]] = maxbin[ii] - 1
         else:
             if (bins >= maxbin).any() or (bins < 0).any():
-                return outlierValue
+                return outlier_value
 
-        return self.hist[tuple(bins)]
+        return self.data[tuple(bins)]
 
-    def draw_2d(self, dims=[0, 1], **kwargs):
+    def draw_2d(self, dims=(0, 1), **kwargs):
         """draw the histogram using pcolormesh() (only works for 2D
         histograms currently)
 
         Parameters
         ----------
-        self: type
-            description
         dims: (int,int)
             indices of which dimensions to draw
         kwargs: 
@@ -304,7 +330,7 @@ class Histogram:
         """
         from matplotlib import pyplot
 
-        if self.hist.ndim < 2:
+        if self.data.ndim < 2:
             raise ValueError("Too few dimensions")
 
         if len(dims) != 2:
@@ -312,29 +338,30 @@ class Histogram:
 
         pyplot.pcolormesh(self.bin_lower_edges[dims[0]],
                           self.bin_lower_edges[dims[1]],
-                          self.hist, **kwargs)
+                          self.data, **kwargs)
         pyplot.title(self.name)
-        pyplot.xlabel(self.axisNames[dims[0]])
-        pyplot.ylabel(self.axisNames[dims[1]])
+        pyplot.xlabel(self.axis_names[dims[0]])
+        pyplot.ylabel(self.axis_names[dims[1]])
 
     def draw_1d(self, dim=0, **kwargs):
         from matplotlib import pyplot
 
         # todo fix this to work properly with dim argument!
-        pyplot.plot(self.bin_centers(dim), self.hist, drawstyle='steps-mid',
+        pyplot.plot(self.bin_centers(dim), self.data, drawstyle='steps-mid',
                     **kwargs)
 
     def interpolate(self, nbins):
         """
         Change the shape of the histogram using an n-dimensional
-        interpolation function.
+        interpolation function (via `ndimage.map_coordinates`).
 
-        Arguments:
-
-        - `bins`: a tuple of the new number of bins to interpolate
-          this histogram over (e.g. if the original histogram was
-          (100,100), setting bins to (200,200) would provide double
-          the resolution, with the interviening bins interpolated.
+        Parameters
+        ----------
+        nbins: tuple of int
+            a tuple of the new number of bins to interpolate
+            this histogram over (e.g. if the original histogram was
+            (100,100), setting bins to (200,200) would provide double
+            the resolution, with the interviening bins interpolated.
         """
 
         oldbins = self._nbins
@@ -345,5 +372,11 @@ class Histogram:
                            for X in range(len(nbins))])
 
         self._nbins = nbins
-        self.hist = ndimage.map_coordinates(self.hist, coords)
+        self.data = ndimage.map_coordinates(self.data, coords)
         self._binLowerEdges = None  # need to be recalculated
+
+    @property
+    def hist(self):
+        """ for backward compatibility. Use `Histogram.data` for read/write
+        access"""
+        return self.data
\ No newline at end of file


Test patch:
diff --git a/ctapipe/utils/tests/test_fitshistogram.py b/ctapipe/utils/tests/test_fitshistogram.py
index c41e8aa1..a1197d56 100644
--- a/ctapipe/utils/tests/test_fitshistogram.py
+++ b/ctapipe/utils/tests/test_fitshistogram.py
@@ -1,6 +1,17 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 import numpy as np
-from ..fitshistogram import Histogram
+import pytest
+from ctapipe.utils.fitshistogram import Histogram
+
+def compare_histograms(hist1: Histogram, hist2: Histogram):
+    """ check that 2 histograms are identical in value """
+    assert hist1.ndims == hist2.ndims
+    assert (hist1.axis_names == hist2.axis_names).all()
+    assert (hist1.data == hist2.data).all
+
+    for ii in range(hist1.ndims):
+        assert np.isclose(hist1.bin_lower_edges[ii],
+                          hist2.bin_lower_edges[ii]).all()
 
 
 def test_histogram_str():
@@ -57,15 +68,50 @@ def test_outliers():
     """
     H = Histogram(nbins=[5, 10], ranges=[[-2.5, 2.5], [-1, 1]])
     H.fill(np.array([[1, 1], ]))
-    val1 = H.get_value((100, 100), outlierValue=-10000)
-    val2 = H.get_value((-100, 0), outlierValue=None)
+    val1 = H.get_value((100, 100), outlier_value=-10000)
+    val2 = H.get_value((-100, 0), outlier_value=None)
     assert val1 == -10000
     assert val2 == 0
 
 
-def test_histogram_write_fits():
+def get_temp_fits_filename():
+
+    return fn
+
+@pytest.fixture(scope='session')
+def histogram_file(tmpdir_factory):
+    """ a fixture that fetches a temporary output dir/file for a test
+    histogram"""
+    return str(tmpdir_factory.mktemp('data').join('histogram_test.fits'))
+
+
+def test_histogram_fits(histogram_file):
     """
     Write to fits,read back, and check
     """
-    # TODO: implement
-    pass
+
+    hist = Histogram(nbins=[5, 11], ranges=[[-2.5, 2.5], [-1, 1]])
+    hist.fill(np.array([[0, 0],
+                        [0, 0.5]]))
+
+    hist.to_fits().writeto(histogram_file, clobber=True)
+    newhist = Histogram.from_fits(histogram_file)
+
+    # check that the values are the same
+    compare_histograms(hist, newhist)
+
+
+def test_histogram_interpolate():
+    H = Histogram(nbins=[5, 11], ranges=[[-2.5, 2.5], [-1, 1]])
+    H.fill(np.array([[0, 0],
+                     [0,0.5]]))
+
+    for testpoint in [(0,0), (0,1), (1,0), (3,3)]:
+        val0 = H.get_value(testpoint)
+        H.interpolate((10,22))
+        val1 = H.get_value(testpoint)
+        H.interpolate((5,11))
+        val2 = H.get_value(testpoint)
+
+        # at least check the interpolation is undoable
+        assert np.isclose(val0[0], val2[0])
\ No newline at end of file