Index: trunk/npstat/swig/npstat_fitting.py =================================================================== --- trunk/npstat/swig/npstat_fitting.py (revision 0) +++ trunk/npstat/swig/npstat_fitting.py (revision 851) @@ -0,0 +1,28 @@ +""" +This module provides a number of useful wrapper classes and functions +related to fitting probability densities with Minuit2. It requires +"minuit2py" package. +""" + +__author__ = "Igor Volobouev (i.volobouev@ttu.edu)" +__version__ = "1.0" +__date__ = "Oct 30 2022" + +import numpy as np +import npstat as ns + +# The function to minimize should be represented by +# a python callable. If the gradient is not calculated +# explicitly, the callable should just return the value +# of the function. Subsequently, the callable should +# be used as an argument of the "Minuit2PySimpleFCN" +# class constructor. +class NegativeLogLikelihood: + def __init__(self, sample, distroConstructor, minDensity=2.2251e-308): + self.sample = sample + self.distroConstructor = distroConstructor + self.minDensity = minDensity + def __call__(self, *args): + distro = self.distroConstructor(*args) + densities = ns.scanDensity1D(distro, self.sample) + return -np.sum(np.log(np.maximum(densities, self.minDensity))) Index: trunk/npstat/swig/npstat_plotting.py =================================================================== --- trunk/npstat/swig/npstat_plotting.py (revision 850) +++ trunk/npstat/swig/npstat_plotting.py (revision 851) @@ -1,134 +1,144 @@ +""" +This module provides a number of useful wrapper classes and functions +related to plotting histograms and densities implemented in the "NPStat" +nonparametric statistics package (Python 3 version). +""" + +__author__ = "Igor Volobouev (i.volobouev@ttu.edu)" +__version__ = "1.0" +__date__ = "Oct 1 2022" + import npstat import numpy as np from npstat_utils import BoxND, arrayValueTypeAsNumpy, arrayNDFromNumpy def _remove_keys(keys, dictionary): for key in keys: if key in dictionary: del dictionary[key] def colormeshHisto2D(ax, h, **options): # We will flip the y axis in order to use the "normal" plotting # convention instead of the matrix convention assert h.dim() == 2, "Not a 2-d histogram" xedges = h.axis(0).binEdges() yedges = h.axis(1).binEdges(True) bins = npstat.arrayNDToNumpy(h.binContents(), False) xv, yv = np.meshgrid(xedges, yedges, indexing='ij') # Get our options before passing **options to the pcolormesh xlabel = options.get('xlabel', h.axis(0).label()) ylabel = options.get('ylabel', h.axis(1).label()) title = options.get('title', h.title()) _remove_keys(('xlabel', 'ylabel', 'title'), options) mesh = ax.pcolormesh(xv, yv, np.flip(bins, 1), **options) ax.set_xlabel(xlabel) ax.set_ylabel(ylabel) ax.set_title(title) return mesh def colormeshDensity2D(ax, distro, bBox, shape, **options): # We will flip the y axis in order to use the "normal" plotting # convention instead of the matrix convention assert distro.dim() == 2, "Not a 2-d distribution" assert len(bBox) == 2, "Not a 2-d bounding box" assert len(shape) == 2, "Not a 2-d shape" xedges = np.linspace(bBox[0].min(), bBox[0].max(), shape[0]+1) yedges = np.linspace(bBox[1].max(), bBox[1].min(), shape[1]+1) # The following makes a numpy array scan = npstat.scanDensityND(distro, bBox, shape) xv, yv = np.meshgrid(xedges, yedges, indexing='ij') xlabel = options.get('xlabel', 'X') ylabel = options.get('ylabel', 'Y') title = options.get('title', None) _remove_keys(('xlabel', 'ylabel', 'title'), options) mesh = ax.pcolormesh(xv, yv, np.flip(scan, 1), **options) ax.set_xlabel(xlabel) ax.set_ylabel(ylabel) if title is not None: ax.set_title(title) return mesh def surfacePlotDensity2D(ax, distro, bBox, shape, **options): # We will flip the y axis in order to use the "normal" plotting # convention instead of the matrix convention assert distro.dim() == 2, "Not a 2-d distribution" assert len(bBox) == 2, "Not a 2-d bounding box" assert len(shape) == 2, "Not a 2-d shape" xedges = np.linspace(bBox[0].min(), bBox[0].max(), shape[0]) yedges = np.linspace(bBox[1].max(), bBox[1].min(), shape[1]) # The following makes a numpy array scan = npstat.scanDensityND(distro, bBox, shape, True) xv, yv = np.meshgrid(xedges, yedges, indexing='ij') xlabel = options.get('xlabel', 'X') ylabel = options.get('ylabel', 'Y') zlabel = options.get('zlabel', 'Density') title = options.get('title', None) _remove_keys(('xlabel', 'ylabel', 'zlabel', 'title'), options) surf = ax.plot_surface(xv, yv, np.flip(scan, 1), **options) ax.set_xlabel(xlabel) ax.set_ylabel(ylabel) ax.set_zlabel(zlabel) if title is not None: ax.set_title(title) return surf def colormeshNumpyArray(ax, arr, bBox, **options): assert arr.ndim == 2, "Not a 2-d array" assert len(bBox) == 2, "Not a 2-d bounding box" shape = arr.shape xedges = np.linspace(bBox[0].min(), bBox[0].max(), shape[0]+1) flip = options.get('flipY', False) or options.get('flipy', False) if flip: yedges = np.linspace(bBox[1].max(), bBox[1].min(), shape[1]+1) else: yedges = np.linspace(bBox[1].min(), bBox[1].max(), shape[1]+1) xv, yv = np.meshgrid(xedges, yedges, indexing='ij') xlabel = options.get('xlabel', 'X') ylabel = options.get('ylabel', 'Y') title = options.get('title', None) _remove_keys(('xlabel', 'ylabel', 'title', 'flipY', 'flipy'), options) if flip: mesh = ax.pcolormesh(xv, yv, np.flip(arr, 1), **options) else: mesh = ax.pcolormesh(xv, yv, arr, **options) ax.set_xlabel(xlabel) ax.set_ylabel(ylabel) if title is not None: ax.set_title(title) return mesh def plotHistoOutline1D(ax, h, **options): xlabel = options.get('xlabel', h.axis(0).label()) ylabel = options.get('ylabel', h.accumulatedDataLabel()) title = options.get('title', h.title()) color = options.get('color', 'k') width = h.axis(0).length() xmin = options.get('xmin', h.axis(0).min() - 0.05*width) xmax = options.get('xmax', h.axis(0).max() + 0.05*width) x, y = npstat.histoOutline1D(h) lowerbound = min(y) upperbound = max(y) if lowerbound > 0.0: lowerbound = 0.0 yrange = upperbound - lowerbound if yrange > 0.0: up = upperbound + yrange*0.1 down = lowerbound - yrange*0.05 else: up = 1.0 down = -0.25 ymin = options.get('ymin', down) ymax = options.get('ymax', up) _remove_keys(('xlabel', 'ylabel', 'title', 'color', 'xmin', 'xmax', 'ymin', 'ymax'), options) ax.plot(x, y, color, **options) ax.set_xlim(xmin, xmax) ax.set_ylim(ymin, ymax) ax.set_xlabel(xlabel) ax.set_ylabel(ylabel) ax.set_title(title) def uncertaintyBandPolygon(x, ymin, ymax): newx = np.concatenate((x, np.flip(x))) newy = np.concatenate((ymin, np.flip(ymax))) return newx, newy