Introduce TMath::KNNDensity and deprecate TH1K

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

@@ -1,7 +1,7 @@
 import ROOT
 import gc
 import os
 import unittest
 
 
 def _leak(obj):
@@ -16,7 +16,7 @@
     def test_tstyle_memory_management(self):
         """Regression test for https://github.com/root-project/root/issues/16918"""
 
         h1 = ROOT.TH1F("h1", "", 100, 0, 10)
 
         style = ROOT.TStyle("NewSTYLE", "")
         groot = ROOT.ROOT.GetROOT()
@@ -27,19 +27,19 @@
         """Regression test for https://github.com/root-project/root/issues/16942"""
         # The test is just that the memory regulation works correctly and the
         # application does not segfault
         f2 = ROOT.TF2("f2", "sin(x)*sin(y)/x/y")
 
     def test_tf3_memory_regulation(self):
         """Make sure TF3 is properly managed by the memory regulation logic"""
         # The test is just that the memory regulation works correctly and the
         # application does not segfault
         f3 = ROOT.TF3("f3","[0] * sin(x) + [1] * cos(y) + [2] * z",0,10,0,10,0,10)
 
     def test_tcolor_memory_regulation(self):
         """Make sure TColor is properly managed by the memory regulation logic"""
         # The test is just that the memory regulation works correctly and the
         # application does not segfault
         c = ROOT.TColor(42, 42, 42)
 
     def test_ttree_clone_in_file_context(self):
         """Test that CloneTree() doesn't give the ownership to Python when
@@ -49,8 +49,8 @@
 
         ttree = ROOT.TTree("tree", "tree")
 
         with ROOT.TFile(filename, "RECREATE") as infile:
             ttree_clone = ttree.CloneTree()
 
         os.remove(filename)
 
@@ -88,8 +88,7 @@
             "TH1I": ("h", "h", 10, 0, 10),
             "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 <[email protected]>
+
+
+\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

math/mathcore/test/CMakeLists.txt

@@ -95,6 +95,7 @@ ROOT_ADD_GTEST(testRootFinder testRootFinder.cxx  LIBRARIES ${Libraries})
 
 ROOT_ADD_GTEST(testKahan testKahan.cxx LIBRARIES Core MathCore)
 ROOT_ADD_GTEST(testDelaunay2D testDelaunay2D.cxx LIBRARIES Core MathCore)
+ROOT_ADD_GTEST(testKNNDensity testKNNDensity.cxx LIBRARIES Core MathCore Hist)
 
 if(clad)
   ROOT_ADD_GTEST(CladDerivatorTests CladDerivatorTests.cxx LIBRARIES Core MathCore)

math/mathcore/test/testKNNDensity.cxx

@@ -0,0 +1,54 @@
+#include "TMath.h"
+#include "TH1K.h"
+
+#include "gtest/gtest.h"
+
+TEST(DNNDensity, BasicTests)
+{
+   auto check = [](std::vector<double> values, int k, double dist, int k_actual) {
+      double query = 0.;
+      std::span<const double> queryspan{&query, 1};
+
+      double output = 0.;
+      std::span<double> outspan{&output, 1};
+
+      // To compute the reference value
+      std::size_t n = values.size();
+      const double factor = n * 0.5 / (n + 1.0);
+
+      TH1K hist{"hist", "hist", 19, -9.5, 9.5};
+      hist.SetKOrd(k);
+
+      double dmin = 1.;
+      TMath::KNNDensity(values, queryspan, outspan, k, k == 0 ? dmin : 0.);
+      for (auto &x : values) {
+         hist.Fill(x);
+      }
+      double ref = factor == 0. ? 0. : factor * k_actual * (1. / std::abs(dist));
+      double th1kval = hist.GetBinContent(10);
+
+      std::cout << output << std::endl;
+      std::cout << ref << std::endl;
+      std::cout << th1kval << std::endl;
+      std::cout << std::endl;
+   };
+
+   check({}, /*k*/ 0, /*dist*/ 0., /*k_actual*/ 0);
+   check({}, /*k*/ 1, /*dist*/ 0., /*k_actual*/ 0);
+   check({}, /*k*/ 2, /*dist*/ 0., /*k_actual*/ 0);
+
+   std::vector<double> values{-3.1, -2.4, -1.5, 0., 0.9, 1.9, 2.9, 4.4};
+   std::vector<double> valuesabs = values;
+   std::sort(valuesabs.begin(), valuesabs.end(), [](double a, double b) { return std::abs(a) < std::abs(b); });
+
+   for (std::size_t i = 0; i < values.size() - 1; ++i) {
+      check(values, /*k=*/i + 1, /*d*/ valuesabs[std::max(i, std::size_t(1))], /*k_actual*/ i + 1);
+   }
+
+   std::cout << "_" << std::endl;
+
+   // check({0., 0.1, }, 0, 1,  2);
+   // check({3.5, 2.5, 1.5, 0., 1., 2., 3., 4.}, 0, 1,  2);
+
+   check({0., 0.2, 0.4, 0.6, 0.8, 1., 1.2, 1.4}, 0, 1.2, 7);
+}

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!

tutorials/hist/index.md

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