Source code for ixpeobssim.core.hist

#!/usr/bin/env python
# Copyright (C) 2018--2020, the ixpeobssim team.
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 3 of the License, or
# (at your option) any later version.
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# GNU General Public License for more details.
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

"""Histogram facilities.

from __future__ import print_function, division

import numbers

from import fits
import matplotlib
import numpy
import scipy.stats
import scipy.signal

from ixpeobssim.core.fitting import fit_histogram, USE_ABSOLUTE_SIGMA
from ixpeobssim.instrument.gpd import gpd_map_binning, GPD_PHYSICAL_HALF_SIDE_X,\
from ixpeobssim.utils.matplotlib_ import plt, setup_gca, draggable_colorbar
from ixpeobssim.utils.logging_ import logger

# pylint: disable=invalid-name, too-many-arguments, attribute-defined-outside-init, too-many-instance-attributes

[docs] class xHistogramBase: """Base class for an n-dimensional histogram. This interface to histograms is profoundly different for the minimal numpy/matplotlib approach, where histogramming methods return bare vectors of bin edges and counts. The main underlying ideas are * we keep track of the bin contents, the bin entries and the sum of the weights squared (the latter for the purpose of calculating the errors); * we control the axis label and the plotting styles; * we provide two separate interfaces, fill() and set_content(), to fill the histogram from either unbinned or binned data; * we support the basic arithmetics (addition, subtraction and multiplication by a scalar); * we support full data persistence (I/O) in FITS format. Note that this base class is not meant to be instantiated directly, and the interfaces to concrete histograms of specific dimensionality are defined in the sub-classes. Parameters ---------- binning : n-tuple of array the bin edges on the different axes. labels : n-tuple of strings the text labels for the different axes. """ PLOT_OPTIONS = dict() def __init__(self, binning, labels): """Constructor. """ assert len(labels) == len(binning) + 1 self.binning = tuple(binning) self.labels = list(labels) self.shape = tuple(len(bins) - 1 for bins in self.binning) self.num_axes = len(self.binning) self.content = numpy.zeros(shape=self.shape, dtype=float) self.entries = numpy.zeros(shape=self.shape, dtype=float) self.sumw2 = numpy.zeros(shape=self.shape, dtype=float)
[docs] def fill(self, *data, weights=None): """Fill the histogram from unbinned data. Note this method is returning the histogram instance, so that the function call can be chained. """ if weights is None: weights = numpy.ones(data[0].shape, dtype=float) elif isinstance(weights, numbers.Number): weights = numpy.full(data[0].shape, weights, dtype=float) data = numpy.vstack(data).T content, _ = numpy.histogramdd(data, bins=self.binning, weights=weights) entries, _ = numpy.histogramdd(data, bins=self.binning) sumw2, _ = numpy.histogramdd(data, bins=self.binning, weights=weights**2.) self.content += content self.entries += entries self.sumw2 += sumw2 return self
[docs] def set_content(self, content, entries=None, errors=None): """Set the bin contents programmatically from binned data. Note this method is returning the histogram instance, so that the function call can be chained. """ assert content.shape == self.shape self.content = content if entries is not None: assert entries.shape == self.shape self.entries = entries if errors is not None: self.set_errors(errors) return self
[docs] def errors(self): """Return the bin errors. """ return numpy.sqrt(self.sumw2)
[docs] def set_errors(self, errors): """ """ assert errors.shape == self.shape self.sumw2 = errors**2.
[docs] def bin_centers(self, axis=0): """Return the bin centers for a specific axis. """ return 0.5 * (self.binning[axis][1:] + self.binning[axis][:-1])
[docs] def bin_widths(self, axis=0): """Return the bin widths for a specific axis. """ return numpy.diff(self.binning[axis])
def _axis_complement(self, axis=0): """Return a tuple of integers identifying all the histogram axes but the one passed as an arguments. If, e.g, num_axes is 3, calling _axis_complement(1) will return (0, 2). This is used to sum and/or average on all the histogram dimensions except for a given one, see, e.g., the mean() and rms() methods below. """ return tuple(ax for ax in range(self.num_axes) if ax != axis)
[docs] def mean(self, axis=0): """Calculate the (binned) mean eastimate along a given axis. """ content = numpy.sum(self.content, axis=self._axis_complement(axis)) return numpy.sum(self.bin_centers(axis) * content) / self.sum()
[docs] def rms(self, axis=0): """Calculate the (binned) rms eastimate along a given axis. .. warning:: Mind this is not using anything clever (e.g. the Welford algorithm) and is not particularly numerically stable. """ bin_centers = self.bin_centers(axis) content = numpy.sum(self.content, axis=self._axis_complement(axis)) _sum = numpy.sum(bin_centers * content) / self.sum() _sum2 = numpy.sum(bin_centers**2. * content) / self.sum() return numpy.sqrt(_sum2 - _sum**2.)
[docs] @staticmethod def bisect(binning, values, side='left'): """Return the indices corresponding to a given array of values for a given binning. """ return numpy.searchsorted(binning, values, side) - 1
[docs] def find_bin(self, *coords): """Find the bin corresponding to a given set of "physical" coordinates on the histogram axes. This returns a tuple of integer indices that can be used to address the histogram content. """ return tuple([self.bisect(binning, value) for binning, value in zip(self.binning, coords)])
[docs] def find_bin_value(self, *coords): """Find the histogram content corresponding to a given set of "physical" coordinates on the histogram axes. """ return self.content[self.find_bin(*coords)]
[docs] def num_entries(self): """Return the number of entries in the histogram. """ return self.entries.sum()
[docs] def sum(self, axis=None): """return the sum of weights in the histogram. """ return self.content.sum(axis)
[docs] def empty_copy(self): """Create an empty copy of a histogram. """ return self.__class__(*self.binning, *self.labels)
[docs] def copy(self): """Create a full copy of a histogram. """ hist = self.empty_copy() hist.set_content(self.content.copy(), self.entries.copy()) return hist
def __add__(self, other): """Histogram addition. """ hist = self.empty_copy() hist.set_content(self.content + other.content, self.entries + other.entries, numpy.sqrt(self.sumw2 + other.sumw2)) return hist def __sub__(self, other): """Histogram subtraction. """ hist = self.empty_copy() hist.set_content(self.content - other.content, self.entries + other.entries, numpy.sqrt(self.sumw2 + other.sumw2)) return hist def __mul__(self, value): """Histogram multiplication by a scalar. Args ---- value : array_like The scale factor for the multiplication---must be either a scalar or an array of the same shape of the histogram content. """ hist = self.empty_copy() hist.set_content(self.content * value, self.entries, self.errors() * value) return hist def __rmul__(self, value): """Histogram multiplication by a scalar. """ return self.__mul__(value)
[docs] def set_axis_label(self, axis, label): """Set the label for a given axis. """ self.labels[axis] = label
def _plot(self, **kwargs): """No-op plot() method, to be overloaded by derived classes. """ raise NotImplementedError('_plot() not implemented for %s' % self.__class__.__name__)
[docs] def plot(self, **kwargs): """Plot the histogram. """ for key, value in self.PLOT_OPTIONS.items(): kwargs.setdefault(key, value) self._plot(**kwargs) setup_gca(xmin=self.binning[0][0], xmax=self.binning[0][-1], xlabel=self.labels[0], ylabel=self.labels[1])
[docs] @staticmethod def label_keyword(axis): """Header keyword for the axis labels. """ return 'LABEL%d' % axis
[docs] @staticmethod def entries_hdu_name(): """Extension name for the bin entries. """ return 'ENTRIES'
[docs] @staticmethod def sumw2_hdu_name(): """Extension name for the bin entries. """ return 'SUMW2'
[docs] @staticmethod def binning_hdu_name(axis): """Extension name for the axis binnings. """ return 'BINNING%d' % axis
[docs] @staticmethod def binning_col_name(): """Column name for the axis binnings. """ return 'EDGES'
[docs] @classmethod def from_file(cls, file_path): """ Load the histogram from a FITS file. Note that we transpose the image data back at read time---see the comment in the save() method. """'Loading histogram from %s...', file_path) with as hdu_list: content = hdu_list['Primary'].data.T entries = hdu_list[cls.entries_hdu_name()].data.T errors = numpy.sqrt(hdu_list[cls.sumw2_hdu_name()].data.T) num_axes = len(content.shape) edges = [hdu_list[cls.binning_hdu_name(i)].data[cls.binning_col_name()] \ for i in range(num_axes)] labels = [hdu_list['Primary'].header[cls.label_keyword(i)] for i in range(num_axes + 1)] hist = cls(*edges, *labels) hist.set_content(content, entries, errors) return hist
[docs] def save(self, file_path, overwrite=True, **header_keywords): """ Save the histogram (with edges) to a FITS file. Note that all the image data are transposed so that the thing can be correctly visualized with standard FITS viewers---at least in two dimensions. """'Saving histogram to %s...', file_path) # Save the histogram content in a series of images. hdu_list = [ fits.PrimaryHDU(self.content.T), fits.ImageHDU(self.entries.T, name=self.entries_hdu_name()), fits.ImageHDU(self.sumw2.T, name=self.sumw2_hdu_name()) ] for key, value in header_keywords.items(): hdu_list[0].header.set(key, value) for i, label in enumerate(self.labels): hdu_list[0].header.set(self.label_keyword(i), label) for i, binning in enumerate(self.binning): col = fits.Column(name=self.binning_col_name(), array=binning, format='E') hdu = fits.BinTableHDU.from_columns([col]) = self.binning_hdu_name(i) hdu_list.append(hdu) hdu_list = fits.HDUList(hdu_list) hdu_list.writeto(file_path, overwrite=overwrite)
[docs] class xHistogram1d(xHistogramBase): """Container class for one-dimensional histograms. """ PLOT_OPTIONS = dict(lw=1.25, alpha=0.4, histtype='stepfilled') def __init__(self, xbins, xlabel='', ylabel='Entries/bin'): """Constructor. """ xHistogramBase.__init__(self, (xbins, ), [xlabel, ylabel]) def _plot(self, **kwargs): """Overloaded make_plot() method. """ plt.hist(self.bin_centers(0), self.binning[0], weights=self.content, **kwargs)
[docs] def errorbar_data(self): """Return the x, y, dy arrays that can be used to build a scatter plot (with errors) from the histogram. """ return self.bin_centers(0), self.content, self.errors()
[docs] def errorbar(self, **kwargs): """Plot the histogram as a scatter plot. """ kwargs.setdefault('fmt', 'o') x, y, dy = self.errorbar_data() plt.errorbar(x, y, dy, **kwargs) setup_gca(xlabel=self.labels[0], ylabel=self.labels[1])
[docs] def scatter_plot(self): """Turn the histogram into a scatter plot. .. warning:: This is to be removed. """ return xScatterPlot(self.bin_centers(0), self.content, xlabel=self.labels[0])
[docs] def gaussian_kde_smooth(self, bandwidth=2.): """Create a copy of the histogram where the weights are smoothed with a gaussian kernel density estimatore. Args ---- bandwidth : float The sigma of the gaussian kernel, in (fractional) number of bins. """ n = int(5. * bandwidth + 0.5) kernel_domain = numpy.linspace(-n, n, 2 * n + 1) kernel = scipy.stats.norm(0, bandwidth).pdf(kernel_domain) hist = self.copy() hist.content = numpy.convolve(self.content, kernel, 'same') return hist
[docs] def fit(self, model, p0=None, sigma=None, xmin=-numpy.inf, xmax=numpy.inf, absolute_sigma=USE_ABSOLUTE_SIGMA, check_finite=True, method=None, verbose=True, **kwargs): """Fit the histogram with a model. """ return fit_histogram(model, self, p0, sigma, xmin, xmax, absolute_sigma, check_finite, method, verbose, **kwargs)
[docs] class xScatterPlot: """Small class encapsulating a scatter plot. Technically speaking, this would not belong here, as a scatter plot is not strictly related to any of the histogram classes, but a 1-dimensional histogram with errors can techincally be turned into a scatter plot, so we introduce the concept here for completeness. .. warning:: Consider removing this class. """ def __init__(self, x, y, dy=None, dx=None, xlabel=None, ylabel=None): """Constructor. """ self.x = x self.y = y self.dy = dy self.dx = dx self.labels = (xlabel, ylabel)
[docs] def plot(self, **kwargs): """Plot the scatter plot. """ kwargs.setdefault('fmt', 'o') kwargs.setdefault('markeredgecolor', 'black') plt.errorbar(self.x, self.y, self.dy, self.dx, **kwargs) setup_gca(xlabel=self.labels[0], ylabel=self.labels[1])
[docs] class xHistogram2d(xHistogramBase): """Container class for two-dimensional histograms. """ PLOT_OPTIONS = dict(cmap=plt.get_cmap('hot')) def __init__(self, xbins, ybins, xlabel='', ylabel='', zlabel='Entries/bin'): """Constructor. """ xHistogramBase.__init__(self, (xbins, ybins), [xlabel, ylabel, zlabel]) def _plot(self, logz=False, **kwargs): """Overloaded make_plot() method. """ x, y = (v.flatten() for v in numpy.meshgrid(self.bin_centers(0), self.bin_centers(1))) bins = self.binning w = self.content.T.flatten() if logz: # Hack for a deprecated functionality in matplotlib 3.3.0 # Parameters norm and vmin/vmax should not be used simultaneously # If logz is requested, we intercent the bounds when created the norm # and refrain from passing vmin/vmax downstream. vmin = kwargs.pop('vmin', None) vmax = kwargs.pop('vmax', None) kwargs.setdefault('norm', matplotlib.colors.LogNorm(vmin, vmax)) plt.hist2d(x, y, bins, weights=w, **kwargs) colorbar = draggable_colorbar() if self.labels[2] is not None: colorbar.set_label(self.labels[2])
[docs] def gaussian_kde_smooth(self, bandwidth=(2., 2.)): """Create a copy of the histogram where the weights are smoothed with a gaussian kernel density estimatore. .. warning:: This is essentially untested. Args ---- bandwidth : 2-element tuple float The sigma of the gaussian kernel, in (fractional) number of bins. """ nx, ny = [int(5. * bw + 0.5) for bw in bandwidth] dx = numpy.linspace(-nx, nx, 2 * nx + 1) dy = numpy.linspace(-ny, ny, 2 * ny + 1) kernel_domain = numpy.dstack(numpy.meshgrid(dx, dy)) sigma = numpy.array([[bandwidth[0], 0.], [0., bandwidth[1]]]) kernel = scipy.stats.multivariate_normal((0., 0.), sigma).pdf(kernel_domain) hist = self.copy() hist.content = scipy.signal.convolve2d(self.content, kernel, 'same') return hist
[docs] def hslice(self, bin_): """Return the horizontal slice for a given bin. """ hist = xHistogram1d(self.binning[0], self.labels[0]) hist.set_content(self.content[:, bin_], self.entries[:, bin_]) return hist
[docs] def hslices(self): """Return a list of all the horizontal slices. """ return tuple(self.hslice(bin_) for bin_ in range(self.shape[1]))
[docs] def hbisect(self, y): """Return the horizontal slice corresponding to a given y value. """ return self.hslice(self.bisect(self.binning[1], y))
[docs] def vslice(self, bin_): """Return the vertical slice for a given bin. """ hist = xHistogram1d(self.binning[1], self.labels[1]) hist.set_content(self.content[bin_], self.entries[bin_]) return hist
[docs] def vslices(self): """Return a list of all the vertical slices. """ return tuple(self.vslice(bin) for bin in range(self.shape[0]))
[docs] def vbisect(self, x): """Return the vertical slice corresponding to a given y value. """ return self.vslice(self.bisect(self.binning[0], x))
[docs] class xModulationCube2d(xHistogram2d): """Specialized class for modulations cubes. """ def __init__(self, xbins, ybins): """Constructor. """ labels = ['Energy [keV]', '$\\phi$ [rad]', 'Entries/bin'] xHistogram2d.__init__(self, xbins, ybins, *labels)
[docs] class xGpdMap2d(xHistogram2d): """2-dimensional GPD map. """ def __init__(self, nside=10, zlabel='Entries/bin'): """Constructor. """ edges = gpd_map_binning(GPD_PHYSICAL_HALF_SIDE_X, GPD_PHYSICAL_HALF_SIDE_Y, nside) labels = ['x [mm]', 'y [mm]', zlabel] xHistogram2d.__init__(self, *edges, *labels)
[docs] class xHistogram3d(xHistogramBase): """Container class for three-dimensional histograms. """ def __init__(self, xbins, ybins, zbins, xlabel='', ylabel='', zlabel='', wlabel='Entries/bin'): """Constructor. """ xHistogramBase.__init__(self, (xbins, ybins, zbins), [xlabel, ylabel, zlabel, wlabel])
[docs] class xGpdMap3d(xHistogram3d): """Three-dimensional histogram where the first two axes represent the GPD active area in Physical coordinates. """ def __init__(self, nside, zbins, zlabel='', wlabel='Entries/bin'): """Constructor. """ edges = gpd_map_binning(GPD_PHYSICAL_HALF_SIDE_X, GPD_PHYSICAL_HALF_SIDE_Y, nside) xHistogram3d.__init__(self, *edges, zbins, 'x [mm]', 'y [mm]', zlabel, wlabel)