Introduce TMath::KNNDensity and deprecate TH1K

The `TH1K` class is deprecated and will be removed in 6.40. It did not
implement the `TH1` interface consistently, and limited the usability of
the k-neighbors method it implemented by closely coupling the algorithm
with the histogram class. Instead, one should use the new
`TMath::KNNDensity` function that implements the same mathematical
logic.

Author: guitargeek

README/ReleaseNotes/v638/index.md

@@ -10,6 +10,7 @@
 * The `rpath` build option is deprecated. It is now without effect.
   Relative RPATHs to the main ROOT libraries are unconditionally appended to all ROOT executables and libraries if the operating system supports it.
   If you want a ROOT build without RPATHs, use the canonical CMake variable `CMAKE_SKIP_INSTALL_RPATH=TRUE`.
+* The `TH1K` class is deprecated and will be removed in 6.40. It did not implement the `TH1` interface consistently, and limited the usability of the k-neighbors method it implemented by closely coupling the algorithm with the histogram class. Please use the new `TMath::KNNDensity` function that implements the same mathematical logic.
 
 ## Core Libraries
 * Behavior change: when selecting a template instantiation for a dictionary, all the template arguments have to be fully defined - the forward declarations are not enough any more. The error prompted by the dictionary generator will be `Warning: Unused class rule: MyTemplate<MyFwdDeclaredClass>`.

hist/hist/inc/TH1K.h

@@ -68,6 +68,6 @@ class TH1K : public TH1, public TArrayF {
    void    SetKOrd(Int_t k){fKOrd=k;}
 
    ClassDefOverride(TH1K,2)  //1-Dim Nearest Kth neighbour method
-};
+} R__DEPRECATED(6, 40, "Use TMath::KNNDensity");
 
 #endif

math/mathcore/inc/TMath.h

@@ -15,6 +15,8 @@
 #include "TMathBase.h"
 
 #include "TError.h"
+#include "ROOT/RSpan.hxx"
+
 #include <algorithm>
 #include <limits>
 #include <cmath>
@@ -574,6 +576,9 @@ struct Limits {
    Double_t Gamma(Double_t a,Double_t x);
    Double_t GammaDist(Double_t x, Double_t gamma, Double_t mu=0, Double_t beta=1);
    Double_t LnGamma(Double_t z);
+
+   void KNNDensity(std::span<const double> observations, std::span<const double> queries, std::span<double> result,
+                   int k, double dmin = 0.0);
 }
 
 ////////////////////////////////////////////////////////////////////////////////

math/mathcore/src/TMath.cxx

@@ -3196,6 +3196,92 @@ Double_t TMath::VavilovDenEval(Double_t rlam, Double_t *AC, Double_t *HC, Int_t
 }
 
 
+/// \brief Computes a 1D k-nearest neighbor density estimate for a set of query points.
+///
+/// For each query point, this function estimates the local density based on
+/// the distance to the k-th nearest neighbor. A minimum distance threshold
+/// `dmin` is used to avoid division by zero and numerical instability when a
+/// query point coincides with an observation or neighbors are extremely close.
+///
+/// \param observations Data points used for density estimation.
+/// \param queries Points at which to estimate the density.
+/// \param result Output for the density estimates. Must have the same size as `queries`.
+/// \param k Number of nearest neighbors to consider. If k >= number of observations,
+///          it is automatically reduced to n-1.
+/// \param dmin Minimum allowed neighbor distance to prevent division by zero. Default is zero.
+///
+/// The density estimate for each query point x is:
+///
+/// \f[
+/// \hat{p}(x) = \frac{n}{2(n+1)} \frac{k_{\text{actual}}}{d_\(k_{\text{actual}}\)(x)},
+/// \f]
+///
+/// where \(d_\(k_{\text{actual}}\)(x)\) is the distance and \(k_{\text{actual}}\) may be different from \(k\)
+/// when it needs to be greater than \(k\) to keep the \(d > d_{\text{min}} constraint.
+/// If no neighbor is outside `dmin`, the minimum allowed distance will be used as the distance itself.
+///
+/// \note This function implements the math of the former `TH1K` class. There are some differences:
+///
+///   1. The `TH1K` class additionally multiplied all density estimates with the bin width.
+///   2. If the `k` parameter was not specified by the user, the `TH1K` class used the bin width as `dmin`.
+///   3. If the `k` parameter was explicitly specified, the `dmin` parameter was set to 1.e-6.
+void KNNDensity(std::span<const double> observations, std::span<const double> queries, std::span<double> result, int k,
+                double dmin)
+{
+   int n = observations.size();
+   if (n == 0) {
+      Error("KNNDensity", "No observations provided");
+      return;
+   }
+   if (k < 0) {
+      Error("KNNDensity", "k must be positive");
+      return;
+   }
+   if (k >= static_cast<int>(n))
+      k = static_cast<int>(n - 1);
+
+   if (result.size() != queries.size()) {
+      Error("KNNDensity", "Result span size must match query span size");
+      return;
+   }
+
+   // Make a sorted copy of the observations for binary search
+   std::vector<double> sorted_obs(observations.begin(), observations.end());
+   std::sort(sorted_obs.begin(), sorted_obs.end());
+
+   // constant factor in the output
+   const double factor = n * 0.5 / (n + 1.0);
+
+   for (std::size_t i = 0; i < queries.size(); ++i) {
+      double x = queries[i];
+
+      // Find index in sorted_obs
+      int jr = std::distance(sorted_obs.begin(), std::lower_bound(sorted_obs.begin(), sorted_obs.end(), x));
+      int jl = jr - 1;
+
+      double ff = 0.0;
+      int k_actual = 0;
+
+      for (; k_actual + 1 <= k || ff <= dmin; ++k_actual) {
+         double fl = (jl >= 0) ? std::abs(sorted_obs[jl] - x) : std::numeric_limits<double>::max();
+         double fr = (jr < n) ? std::abs(sorted_obs[jr] - x) : std::numeric_limits<double>::max();
+
+         if (jl < 0 && jr >= n)
+            break;
+
+         ff = std::min(fl, fr);
+         fl < fr ? --jl : ++jr;
+      }
+
+      if (ff == 0.0)
+         ff = dmin; // avoid division by zero
+
+      // Modified kNN density formula
+      result[i] = factor * k_actual / ff;
+   }
+}
+
+
 //explicitly instantiate template functions from VecCore
 #ifdef R__HAS_VECCORE
 #include <Math/Types.h>

roottest/graphics/CMakeLists.txt

@@ -112,7 +112,6 @@ set(TEST_DESCRIPTIONS
     "h2_cut a hist"
     "h2proj a hist"
     "histpalettecolor a hist"
-    "hksimple a hist"
     "hlHisto1 a hist"
     "hlHisto2 a hist"
     #"hlHisto4 a hist" TODO: Fails on all platforms!

roottest/graphics/macros/hist/hksimple.C

@@ -80,7 +80,6 @@ The examples below showcase the same functionalities through 3 different impleme
 | hist009_TH1_normalize.C            |                                  |                               | Normalizing a Histogram. |
 | hist010_TH1_two_scales.C           | hist010_TH1_two_scales.py        | hist010_TH1_two_scales_uhi.py|Draw two histograms on one canvas using different y-axis scales. |
 | hist011_TH1_legend_autoplaced.C    |                                  |                               | Automatic placing of the legend. |
-| hist012_TH1_hksimple.C             |                                  |                               | Dynamic filling of TH1K histograms. |
 | hist013_TH1_rebin.C                |                                  |                               | Create a variable bin-width histogram and change bin sizes. |
 | hist014_TH1_cumulative.C           |                                  |                               | Illustrate use of the TH1::GetCumulative method. |
 | hist015_TH1_read_and_draw.C        | hist015_TH1_read_and_draw.py     | hist015_TH1_read_and_draw_uhi.py|Read a 1D histogram from a ROOT File and draw it. |