[TMath] Add Gradient and Laplacian methods for arrays

Author: pitkajuh

math/mathcore/inc/TMath.h

@@ -442,6 +442,10 @@ struct Limits {
    template <typename T> Long64_t  LocMax(Long64_t n, const T *a);
    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);
+
    // Hashing
    ULong_t Hash(const void *txt, Int_t ntxt);
    ULong_t Hash(const char *str);
@@ -1002,6 +1006,82 @@ Iterator TMath::LocMin(Iterator first, Iterator last) {
    return std::min_element(first, last);
 }
 
+////////////////////////////////////////////////////////////////////////////////
+/// \brief Calculate the one-dimensional gradient 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.
+///
+/// \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.
+/// \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)
+{
+   if (!f) {
+      ::Error("TMath::Gradient", "Input parameter f is empty.");
+      return nullptr;
+   } else if (n < 2) {
+      ::Error("TMath::Gradient", "Input parameter n=%lld is smaller than 2.", n);
+      return nullptr;
+   }
+   Long64_t i = 1;
+   T *result = new T[n];
+
+   // Forward difference
+   result[0] = (f[1] - f[0]) / h;
+
+   // Central difference
+   while (i < n - 1) {
+      result[i] = (f[i + 1] - f[i - 1]) / (2 * h);
+      i++;
+   }
+   // Backward difference
+   result[i] = (f[i] - f[i - 1]) / h;
+   return result;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+/// \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.
+///
+/// \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.
+/// \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)
+{
+   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);
+      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);
+
+   // Central difference
+   while (i < n - 1) {
+      result[i] = (f[i + 1] + f[i - 1] - 2 * f[i]) / (4 * h * h);
+      i++;
+   }
+   // Backward difference
+   result[i] = (2 * f[i] - 5 * f[i - 1] + 4 * f[i - 2] - f[i - 3]) / (4 * h * h);
+   return result;
+}
+
 ////////////////////////////////////////////////////////////////////////////////
 /// Returns index of array with the maximum element.
 /// If more than one element is maximum returns first found.

math/mathcore/test/testTMath.cxx

@@ -36,6 +36,32 @@ void testNormCross()
    << endl;
 }
 
+template <typename T>
+void 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};
+
+   // test results
+
+   for (Long64_t i = 0; i < n; 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);
+   }
+
+   delete [] gradient;
+   delete [] laplacian;
+
+}
 
 template <typename T, typename U>
 void testArrayFunctions()
@@ -174,6 +200,15 @@ void testTMath()
    testArrayFunctions<Long_t,Long64_t>();
    testArrayFunctions<Long64_t,Long64_t>();
 
+   cout << "\nArray derivative tests: " << endl;
+
+   testArrayDerivatives<Short_t>();
+   testArrayDerivatives<Int_t>();
+   testArrayDerivatives<Float_t>();
+   testArrayDerivatives<Double_t>();
+   testArrayDerivatives<Long_t>();
+   testArrayDerivatives<Long64_t>();
+
    cout << "\nIterator functions tests: " << endl;
 
    testIteratorFunctions<Short_t>();