[6.40] [math] Add test and improve doxygen of finite difference methods

math/mathcore/inc/TMath.h

@@ -446,8 +446,10 @@ struct Limits {
    template <typename Iterator> Iterator LocMax(Iterator first, Iterator last);
 
    // Derivatives of an array
-   template <typename T> T *Gradient(Long64_t n, T *f, double h = 1);
-   template <typename T> T *Laplacian(Long64_t n, T *f, double h = 1);
+   template <typename T>
+   T *Gradient(Long64_t n, T const *f, double h = 1);
+   template <typename T>
+   T *Laplacian(Long64_t n, T const *f, double h = 1);
 
    // Hashing
    ULong_t Hash(const void *txt, Int_t ntxt);
@@ -612,7 +614,7 @@ inline Double_t TMath::Tan(Double_t x)
    { return tan(x); }
 
 ////////////////////////////////////////////////////////////////////////////////
-/// Returns the hyperbolic sine of `x.
+/// Returns the hyperbolic sine of `x`.
 
 inline Double_t TMath::SinH(Double_t x)
    { return sinh(x); }
@@ -916,7 +918,7 @@ inline Double_t TMath::QuietNaN() {
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-/// Returns a signaling NaN as defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Signaling_NaN).
+/// Returns a signaling NaN [as defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Signaling_NaN).
 
 inline Double_t TMath::SignalingNaN() {
    return std::numeric_limits<Double_t>::signaling_NaN();
@@ -934,15 +936,15 @@ inline Double_t TMath::Infinity() {
 
 template<typename T>
 inline T TMath::Limits<T>::Min() {
-   return (std::numeric_limits<T>::min)();    //N.B. use this signature to avoid class with macro min() on Windows
+   return (std::numeric_limits<T>::min)();    //N.B. use this signature to avoid clashes with macro min() on Windows
 }
 
 ////////////////////////////////////////////////////////////////////////////////
 /// Returns minimum double representation.
 
 template<typename T>
 inline T TMath::Limits<T>::Max() {
-   return (std::numeric_limits<T>::max)();  //N.B. use this signature to avoid class with macro max() on Windows
+   return (std::numeric_limits<T>::max)();  //N.B. use this signature to avoid clashes with macro max() on Windows
 }
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -1014,19 +1016,21 @@ Iterator TMath::LocMin(Iterator first, Iterator last) {
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-/// \brief Calculate the one-dimensional gradient of an array with length n.
+/// \brief Calculate the one-dimensional finite gradient of an array with length n.
+/// It is assumed that the values in the array are spaced uniformly in "x", which
+/// would amount to taking the derivative w.r.t. x with an infinitesimal step size `h`.
+///
 /// The first value in the returned array is a forward difference,
 /// the next n-2 values are central differences, and the last is a backward difference.
 ///
-/// \note Function leads to undefined behavior if n does not match the length of f
 /// \param n the number of points in the array
-/// \param f the array of points.
-/// \param h the step size. The default step size is 1.
+/// \param f the array of points.  It is assumed that these are spaced uniformly in "x".
+/// \param h the distance between the points in `f`.
 /// \return an array of size n with the gradient. Returns nullptr if n < 2 or f empty. Ownership is transferred to the
 /// caller.
 
 template <typename T>
-T *TMath::Gradient(const Long64_t n, T *f, const double h)
+T *TMath::Gradient(const Long64_t n, T const *f, const double h)
 {
    if (!f) {
       ::Error("TMath::Gradient", "Input parameter f is empty.");
@@ -1052,40 +1056,42 @@ T *TMath::Gradient(const Long64_t n, T *f, const double h)
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-/// \brief Calculate the Laplacian of an array with length n.
-/// The first value in the returned array is a forward difference,
-/// the next n-2 values are central differences, and the last is a backward difference.
+/// \brief Calculate second-order finite differences of an array with length n at second-order accuracy.
+/// It is assumed that the values in the array are spaced uniformly in "x", where
+/// the spacing is represented by the argument `h`.
+///
+/// The first value in the returned array is a forward difference at 2nd order accuracy,
+/// the next n-2 values are central differences, and the last is the backward difference.
 ///
-/// \note Function leads to undefined behavior if n does not match the length of f
 /// \param n the number of points in the array
-/// \param f the array of points.
-/// \param h the step size. The default step size is 1.
+/// \param f the array of points. It is assumed that these are spaced uniformly in "x".
+/// \param h the distance between the points in `f`.
 /// \return an array of size n with the laplacian. Returns nullptr if n < 4 or f empty. Ownership is transferred to the
 /// caller.
 
 template <typename T>
-T *TMath::Laplacian(const Long64_t n, T *f, const double h)
+T *TMath::Laplacian(const Long64_t n, T const *f, const double h)
 {
    if (!f) {
       ::Error("TMath::Laplacian", "Input parameter f is empty.");
       return nullptr;
    } else if (n < 4) {
-      ::Error("TMath::Laplacian", "Input parameter n=%lld is smaller than 4.", n);
+      ::Error("TMath::Laplacian", "Need at least four elements, got %lld", n);
       return nullptr;
    }
    Long64_t i = 1;
    T *result = new T[n];
 
    // Forward difference
-   result[0] = (4 * f[2] + 2 * f[0] - 5 * f[1] - f[3]) / (4 * h * h);
+   result[0] = (4 * f[2] + 2 * f[0] - 5 * f[1] - f[3]) / (h * h);
 
    // Central difference
    while (i < n - 1) {
-      result[i] = (f[i + 1] + f[i - 1] - 2 * f[i]) / (4 * h * h);
+      result[i] = (f[i + 1] + f[i - 1] - 2 * f[i]) / (h * h);
       i++;
    }
    // Backward difference
-   result[i] = (2 * f[i] - 5 * f[i - 1] + 4 * f[i - 2] - f[i - 3]) / (4 * h * h);
+   result[i] = (2 * f[i] - 5 * f[i - 1] + 4 * f[i - 2] - f[i - 3]) / (h * h);
    return result;
 }
 
@@ -1516,7 +1522,7 @@ template <typename T> Double_t TMath::ModeHalfSample(Long64_t n, const T *a, con
       const size_t N = std::ceil(n * 0.5);
       const size_t start = jMin;
       const size_t stop = start + n - N + 1; // +1 since we use < and not <=
-      // Find sequentally what v_range is smallest by sliding the half-window
+      // Find sequentially what v_range is smallest by sliding the half-window
       for (size_t i = start; i < stop; i++) {
          Double_t range = values[i + N - 1] - values[i];
          if (range < min_v_range) {

math/mathcore/test/testTMath.cxx

@@ -1,14 +1,14 @@
+#include <TMath.h>
+#include <TError.h>
+
+#include <array>
 #include <iostream>
 #include <vector>
 #include <string>
 #include <cstring>
 #include <typeinfo>
 
-#include <TMath.h>
-#include <TError.h>
-
 using std::cout, std::endl, std::vector, std::sort;
-using namespace TMath;
 
 bool showVector = true;
 
@@ -38,35 +38,45 @@ void testNormCross()
 }
 
 template <typename T>
-void testArrayDerivatives()
+bool testArrayDerivatives()
 {
-   const Long64_t n = 10;
-   const double h = 0.1;
-   T sa[n] = {18, 47, 183, 98, 56, 74, 28, 75, 10, 89};
-   T *gradient = TMath::Gradient(n, sa, h);
-   T *laplacian = TMath::Laplacian(n, sa, h);
-
-   const T gradienta[n] = {290, 825, 255, -635, -120, -140, 5, -90, 70, 790};
-   const T laplaciana[n] = {10875, 2675, -5525, 1075, 1500, -1600, 2325, -2800, 3600, 10000};
+   bool failure = false;
+   constexpr Long64_t n = 10;
+   constexpr double h = 0.1;
+   constexpr std::array<T, n> sa = {18, 47, 183, 98, 56, 74, 28, 75, 10, 89};
+   T *gradient = TMath::Gradient(sa.size(), sa.data(), h);
+   T *laplacian = TMath::Laplacian(sa.size(), sa.data(), h);
+
+   // see https://en.wikipedia.org/wiki/Finite_difference_coefficient
+   constexpr std::array<int, n> gradienta = {290, 825, 255, -635, -120, -140, 5, -90, 70, 790};
+   // central differences are e.g. (sa[x-1] - 2* sa[x] + sa[x+1]) / h^2
+   constexpr std::array<int, n> laplaciana = {43500, 10700, -22100, 4300, 6000, -6400, 9300, -11200, 14400, 40000};
 
    // test results
-
    for (Long64_t i = 0; i < n; i++) {
-      if (gradient[i] != gradienta[i])
+      if (gradient[i] != gradienta[i]) {
          Error("testArrayDerivatives", "For Gradient, different values found at i = %lld", i);
-
-      if (laplacian[i] != laplaciana[i])
-         Error("testArrayDerivatives", "For Laplacian, different values found at i = %lld", i);
+         failure = true;
+      }
+
+      // We have to allow a bit of rounding errors, because TMath internally computes in floating-point:
+      if (abs(static_cast<int>(laplacian[i]) - laplaciana[i]) > 1) {
+         Error("testArrayDerivatives", "For Laplacian, different values found at i = %lld.", i);
+         std::cerr << laplacian[i] << " " << laplaciana[i] << "\n";
+         failure = true;
+      }
    }
 
    delete [] gradient;
    delete [] laplacian;
 
+   return failure;
 }
 
 template <typename T, typename U>
 void testArrayFunctions()
 {
+   using namespace TMath;
    const U n = 10;
    const U k = 3;
    U index[n];
@@ -106,6 +116,7 @@ void testArrayFunctions()
 template <typename T>
 void testIteratorFunctions()
 {
+   using namespace TMath;
    const Long64_t n = 10;
    vector<Int_t> index(n);
    Long64_t is;
@@ -148,9 +159,7 @@ void testPoints(T x, T y)
 
    T dx[4] = {0, 0, 2, 2};
    T dy[4] = {0, 2, 2, 0};
-   cout << "Point(" << x << "," << y << ") IsInside?: "
-        << IsInside( x, y, n, dx, dy) << endl;
-
+   cout << "Point(" << x << "," << y << ") IsInside?: " << TMath::IsInside(x, y, n, dx, dy) << endl;
 }
 
 template <typename T>
@@ -160,7 +169,7 @@ void testPlane()
    T dp2[3] = {1,0,0};
    T dp3[3] = {0,1,0};
    T dn[3];
-   Normal2Plane(dp1, dp2, dp3, dn);
+   TMath::Normal2Plane(dp1, dp2, dp3, dn);
    cout << "Normal: ("
         << dn[0] << ", "
         << dn[1] << ", "
@@ -178,7 +187,8 @@ void testBreitWignerRelativistic()
 
   for (Int_t i=0;i<=nPoints;i++) {
     Double_t currentX = xMinimum+i*xStepSize;
-    cout << "BreitWignerRelativistic(" << currentX << "," << median << "," << gamma << ") = " << BreitWignerRelativistic(currentX,median,gamma) << endl;
+    cout << "BreitWignerRelativistic(" << currentX << "," << median << "," << gamma
+         << ") = " << TMath::BreitWignerRelativistic(currentX, median, gamma) << endl;
   }
 }
 
@@ -244,8 +254,12 @@ void testHalfSampleMode()
    R__ASSERT(TMath::ModeHalfSample(testdata7_n, testdata7, weightdata7) == -1);
 }
 
-void testTMath()
+int main()
 {
+   // TODO: Warning, many tests don't yet return any success or failure information,
+   // so don't trust that they will be red.
+   bool failure = false;
+
    cout << "Starting tests on TMath..." << endl;
 
    cout << "\nNormCross tests: " << endl;
@@ -265,12 +279,11 @@ void testTMath()
 
    cout << "\nArray derivative tests: " << endl;
 
-   testArrayDerivatives<Short_t>();
-   testArrayDerivatives<Int_t>();
-   testArrayDerivatives<Float_t>();
-   testArrayDerivatives<Double_t>();
-   testArrayDerivatives<Long_t>();
-   testArrayDerivatives<Long64_t>();
+   failure |= testArrayDerivatives<Int_t>();
+   failure |= testArrayDerivatives<Float_t>();
+   failure |= testArrayDerivatives<Double_t>();
+   failure |= testArrayDerivatives<Long_t>();
+   failure |= testArrayDerivatives<Long64_t>();
 
    cout << "\nIterator functions tests: " << endl;
 
@@ -298,11 +311,6 @@ void testTMath()
 
    cout << "\nHalfSampleMode tests: " << endl;
    testHalfSampleMode();
-}
-
-int main()
-{
-   testTMath();
 
-   return 0;
+   return failure ? 1 : 0;
 }

roottest/root/math/CMakeLists.txt

@@ -1 +1,3 @@
 ROOTTEST_ADD_TESTDIRS()
+
+ROOT_ADD_GTEST(MathCoreTests MathCoreTests.cxx LIBRARIES ROOT::MathCore ROOT::Hist)

roottest/root/math/MathCoreTests.cxx

@@ -0,0 +1,77 @@
+#include <TMath.h>
+#include <TF1.h>
+
+#include <gtest/gtest.h>
+
+#include <array>
+
+TEST(TMath, Gradient_Laplace)
+{
+   std::array<double, 5> parameters{2, 100., -1., -2., 0.1};
+   std::array<double, 4> parameters1st{};
+   std::array<double, 3> parameters2nd{};
+   for (unsigned int i = 0; i < 4; ++i) {
+      parameters1st[i] = parameters[i + 1] * (i + 1);
+   }
+   for (unsigned int i = 0; i < 3; ++i) {
+      parameters2nd[i] = parameters1st[i + 1] * (i + 1);
+   }
+
+   TF1 poly("poly", "pol4", -10, 10);
+   TF1 poly1st("derivative", "pol3", -10, 10);
+   TF1 poly2nd("secondDerivative", "pol2", -10, 10);
+   poly.SetParameters(parameters.data());
+   poly1st.SetParameters(parameters1st.data());
+   poly2nd.SetParameters(parameters2nd.data());
+
+   constexpr std::size_t nPoint = 10000;
+   std::array<double, nPoint> vx;
+   std::array<double, nPoint> vPoly;
+   for (unsigned int i = 0; i < nPoint; ++i) {
+      const auto x = -10. + 20. / nPoint * i;
+      vx[i] = x;
+      vPoly[i] = poly.Eval(x);
+   }
+
+   auto grad = TMath::Gradient(nPoint, vPoly.data(), 20. / nPoint);
+   auto lap = TMath::Laplacian(nPoint, vPoly.data(), 20. / nPoint);
+
+   auto relativeDiff = [](double val, double ref) {
+      if (ref == 0.) {
+         return std::fabs(val - ref);
+      }
+      return std::fabs(val - ref) / ref;
+   };
+
+   // Check forward/backward differences
+   EXPECT_LT(relativeDiff(grad[0], poly1st.Eval(-10)), 0.001);
+   EXPECT_LT(relativeDiff(grad[nPoint - 1], poly1st.Eval(vx[nPoint - 1])), 0.001);
+   EXPECT_LT(relativeDiff(lap[0], poly2nd.Eval(-10)), 0.001);
+   EXPECT_LT(relativeDiff(lap[nPoint - 1], poly2nd.Eval(vx[nPoint - 1])), 0.001);
+
+   {
+      double squaredDiff_grad = 0.;
+      double maxRelDiff_grad = 0.;
+      double squaredDiff_laplace = 0.;
+      double maxRelDiff_laplace = 0.;
+      // The points on the edges are forward/backward differences, so they will diverge more
+      // Therefore, run the comparison only on the centre
+      for (unsigned int i = 1; i < nPoint - 1; ++i) {
+         const auto x = vx[i];
+         const double diff = poly1st.Eval(x) - grad[i];
+         squaredDiff_grad += diff * diff;
+         maxRelDiff_grad = std::max(relativeDiff(grad[i], poly1st.Eval(x)), maxRelDiff_grad);
+
+         const double diff2 = poly2nd.Eval(x) - lap[i];
+         squaredDiff_laplace += diff2 * diff2;
+         maxRelDiff_laplace = std::max(relativeDiff(lap[i], poly2nd.Eval(x)), maxRelDiff_laplace);
+      }
+
+      // Central differences
+      EXPECT_LT(maxRelDiff_grad, 0.01);
+      EXPECT_LT(std::sqrt(squaredDiff_grad), 0.01);
+
+      EXPECT_LT(maxRelDiff_laplace, 0.01);
+      EXPECT_LT(std::sqrt(squaredDiff_laplace), 0.01);
+   }
+}