Examples/python/scipy.py

0001 # This file is part of the ACTS project.
0002 #
0003 # Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 #
0005 # This Source Code Form is subject to the terms of the Mozilla Public
0006 # License, v. 2.0. If a copy of the MPL was not distributed with this
0007 # file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 """ROOT-free Gaussian fit backend for `ActsExamples::HistogramFitFunction`,
0010 using `scipy.optimize.curve_fit`. Suitable for use in the PyPI distribution,
0011 where the ROOT plugin (and `ActsPlugins::RootHistogramFit`) is not available.
0012 """
0013
0014 import numpy as np
0015 from scipy.optimize import curve_fit
0016
0017
0018 def _gaussian(x, amplitude, mean, sigma):
0019     return amplitude * np.exp(-0.5 * ((x - mean) / sigma) ** 2)
0020
0021
0022 def _gaussian_jac(x, amplitude, mean, sigma):
0023     """Analytic Jacobian of `_gaussian` w.r.t. (amplitude, mean, sigma).
0024
0025     Routes curve_fit's underlying MINPACK call through `_lmder` (exact
0026     derivatives) instead of `_lmdif` (forward-difference approximation).
0027     """
0028     z = (x - mean) / sigma
0029     g = np.exp(-0.5 * z * z)
0030     return np.stack(
0031         [g, amplitude * g * z / sigma, amplitude * g * z * z / sigma], axis=-1
0032     )
0033
0034
0035 def makeScipyHistogramFitFunction(maxfev=200, ftol=1e-6, xtol=1e-6):
0036     """Build a Gaussian-fit `HistogramFitFunction` backed by
0037     `scipy.optimize.curve_fit`.
0038
0039     The returned callable matches `ActsExamples::HistogramFitFunction`'s
0040     signature: `(hist, range) -> Optional[(mean, sigma, meanError,
0041     sigmaError)]`. Drops empty bins rather than weighting them at sigma=1,
0042     mirroring ROOT's "SQ0". The analytic Jacobian plus a bounded
0043     `maxfev`/`ftol`/`xtol` keeps `fitProfiles()` fast versus MINPACK's
0044     defaults.
0045
0046     @param maxfev Maximum number of function evaluations
0047     @param ftol Relative error desired in the sum of squares
0048     @param xtol Relative error desired in the approximate solution
0049     """
0050
0051     def fit(hist, rng):
0052         values = hist.histogram.values()
0053         edges = np.asarray(hist.histogram.axis(0).edges)
0054         centres = 0.5 * (edges[:-1] + edges[1:])
0055
0056         if rng is not None:
0057             xMin, xMax = rng
0058             mask = (centres >= xMin) & (centres <= xMax)
0059             centres = centres[mask]
0060             values = values[mask]
0061
0062         if np.count_nonzero(values) < 3 or values.sum() <= 0:
0063             return None
0064
0065         mean0 = np.average(centres, weights=np.clip(values, 0, None))
0066         sigma0 = max(
0067             np.sqrt(
0068                 np.average((centres - mean0) ** 2, weights=np.clip(values, 0, None))
0069             ),
0070             1e-6,
0071         )
0072         amplitude0 = values.max()
0073
0074         keep = values > 0
0075         fitCentres = centres[keep]
0076         fitValues = values[keep]
0077         errors = np.sqrt(fitValues)
0078
0079         try:
0080             with np.errstate(all="ignore"):
0081                 popt, pcov = curve_fit(
0082                     _gaussian,
0083                     fitCentres,
0084                     fitValues,
0085                     p0=[amplitude0, mean0, sigma0],
0086                     sigma=errors,
0087                     absolute_sigma=True,
0088                     jac=_gaussian_jac,
0089                     maxfev=maxfev,
0090                     ftol=ftol,
0091                     xtol=xtol,
0092                 )
0093         except RuntimeError:
0094             return None
0095
0096         if not np.all(np.isfinite(pcov)):
0097             return None
0098
0099         meanError = float(np.sqrt(pcov[1, 1]))
0100         sigmaError = float(np.sqrt(pcov[2, 2]))
0101         return (float(popt[1]), abs(float(popt[2])), meanError, sigmaError)
0102
0103     return fit