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>`.

bindings/pyroot/pythonizations/python/ROOT/_pythonization/_th1.py

@@ -233,7 +233,6 @@ def _FillWithNumpyArray(self, *args):
     "TH1L",
     "TH1F",
     "TH1D",
-    "TH1K",
     "TProfile",
 ]
 

bindings/pyroot/pythonizations/python/ROOT/_pythonization/_uhi.py

@@ -566,7 +566,6 @@ def _get_sum_of_weights_squared(self) -> np.typing.NDArray[Any]:  # noqa: F821
     "TH1C": _values_by_copy,
     "TH2C": _values_by_copy,
     "TH3C": _values_by_copy,
-    "TH1K": _values_by_copy,
     "TH2K": _values_by_copy,
     "TH3K": _values_by_copy,
     "TProfile": _values_by_copy,

bindings/pyroot/pythonizations/test/memory.py

@@ -89,7 +89,6 @@ def test_objects_ownership_with_subdir(self):
             "TH1L": ("h", "h", 10, 0, 10),
             "TH1F": ("h", "h", 10, 0, 10),
             "TH1D": ("h", "h", 10, 0, 10),
-            "TH1K": ("h", "h", 10, 0, 10),
             "TProfile": ("h", "h", 10, 0, 10),
             "TH2C": ("h", "h", 10, 0, 10, 10, 0, 10),
             "TH2S": ("h", "h", 10, 0, 10, 10, 0, 10),

bindings/pyroot/pythonizations/test/uhi_indexing.py

@@ -12,7 +12,7 @@ def _special_setting(hist):
     # For these classes, SetBinContent works differently than for other classes,
     # as in it does not set the bin content to the specified value, but does some other calculations
     # for that, these classes will be tested differently
-    return isinstance(hist, (ROOT.TProfile, ROOT.TProfile2D, ROOT.TProfile2Poly, ROOT.TProfile3D, ROOT.TH1K))
+    return isinstance(hist, (ROOT.TProfile, ROOT.TProfile2D, ROOT.TProfile2Poly, ROOT.TProfile3D))
 
 
 def _get_index_for_dimension(hist, index):

gui/webgui6/src/TWebCanvas.cxx

@@ -36,7 +36,6 @@
 #include "TF2.h"
 #include "TH1.h"
 #include "TH2.h"
-#include "TH1K.h"
 #include "THStack.h"
 #include "TMultiGraph.h"
 #include "TEnv.h"
@@ -935,27 +934,6 @@ void TWebCanvas::CreatePadSnapshot(TPadWebSnapshot &paddata, TPad *pad, Long64_t
          CreatePadSnapshot(paddata.NewSubPad(), (TPad *)obj, version, nullptr);
       } else if (!process_primitives) {
          continue;
-      } else if (obj->InheritsFrom(TH1K::Class())) {
-         flush_master();
-         TH1K *hist = static_cast<TH1K *>(obj);
-
-         Int_t nbins = hist->GetXaxis()->GetNbins();
-
-         TH1D *h1 = new TH1D("__dummy_name__", hist->GetTitle(), nbins, hist->GetXaxis()->GetXmin(), hist->GetXaxis()->GetXmax());
-         h1->SetDirectory(nullptr);
-         h1->SetName(hist->GetName());
-         hist->TAttLine::Copy(*h1);
-         hist->TAttFill::Copy(*h1);
-         hist->TAttMarker::Copy(*h1);
-         for (Int_t n = 1; n <= nbins; ++n)
-             h1->SetBinContent(n, hist->GetBinContent(n));
-
-         TIter fiter(hist->GetListOfFunctions());
-         while (auto fobj = fiter())
-            h1->GetListOfFunctions()->Add(fobj->Clone());
-
-         paddata.NewPrimitive(obj, iter.GetOption()).SetSnapshot(TWebSnapshot::kObject, h1, kTRUE);
-
       } else if (obj->InheritsFrom(TH1::Class())) {
          flush_master();
 

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

hist/hist/src/TH1K.cxx

@@ -28,6 +28,9 @@ In this method :
              that DistanceToNearestKthNeighbour > BinWidth and K >=3
 
 This class has been implemented by Victor Perevoztchikov <perev@bnl.gov>
+
+
+\deprecated 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.
 */
 
 ClassImp(TH1K);

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

@@ -3195,6 +3195,87 @@ Double_t TMath::VavilovDenEval(Double_t rlam, Double_t *AC, Double_t *HC, Int_t
    return v;
 }
 
+/// \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 TMath::KNNDensity(std::span<const double> observations, std::span<const double> queries, std::span<double> result,
+                       int k, double dmin)
+{
+   int n = observations.size();
+   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 == 0.0 ? 0.0 : factor * k_actual / ff;
+   }
+}
 
 //explicitly instantiate template functions from VecCore
 #ifdef R__HAS_VECCORE