diff --git a/math/mathcore/inc/TMath.h b/math/mathcore/inc/TMath.h index 134694e3cf88b..6b66deea74564 100644 --- a/math/mathcore/inc/TMath.h +++ b/math/mathcore/inc/TMath.h @@ -446,8 +446,10 @@ struct Limits { template Iterator LocMax(Iterator first, Iterator last); // Derivatives of an array - template T *Gradient(Long64_t n, T *f, double h = 1); - template T *Laplacian(Long64_t n, T *f, double h = 1); + template + T *Gradient(Long64_t n, T const *f, double h = 1); + template + 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::signaling_NaN(); @@ -934,7 +936,7 @@ inline Double_t TMath::Infinity() { template inline T TMath::Limits::Min() { - return (std::numeric_limits::min)(); //N.B. use this signature to avoid class with macro min() on Windows + return (std::numeric_limits::min)(); //N.B. use this signature to avoid clashes with macro min() on Windows } //////////////////////////////////////////////////////////////////////////////// @@ -942,7 +944,7 @@ inline T TMath::Limits::Min() { template inline T TMath::Limits::Max() { - return (std::numeric_limits::max)(); //N.B. use this signature to avoid class with macro max() on Windows + return (std::numeric_limits::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 -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 -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 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) { diff --git a/math/mathcore/test/testTMath.cxx b/math/mathcore/test/testTMath.cxx index acdc93ff7a134..5ddd07132c03a 100644 --- a/math/mathcore/test/testTMath.cxx +++ b/math/mathcore/test/testTMath.cxx @@ -1,14 +1,14 @@ +#include +#include + +#include #include #include #include #include #include -#include -#include - using std::cout, std::endl, std::vector, std::sort; -using namespace TMath; bool showVector = true; @@ -38,35 +38,45 @@ void testNormCross() } template -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 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 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 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(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 void testArrayFunctions() { + using namespace TMath; const U n = 10; const U k = 3; U index[n]; @@ -106,6 +116,7 @@ void testArrayFunctions() template void testIteratorFunctions() { + using namespace TMath; const Long64_t n = 10; vector 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 @@ -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(); - testArrayDerivatives(); - testArrayDerivatives(); - testArrayDerivatives(); - testArrayDerivatives(); - testArrayDerivatives(); + failure |= testArrayDerivatives(); + failure |= testArrayDerivatives(); + failure |= testArrayDerivatives(); + failure |= testArrayDerivatives(); + failure |= testArrayDerivatives(); 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; } diff --git a/roottest/root/math/CMakeLists.txt b/roottest/root/math/CMakeLists.txt index 9a26f3f797c44..1314db5ca11aa 100644 --- a/roottest/root/math/CMakeLists.txt +++ b/roottest/root/math/CMakeLists.txt @@ -1 +1,3 @@ ROOTTEST_ADD_TESTDIRS() + +ROOT_ADD_GTEST(MathCoreTests MathCoreTests.cxx LIBRARIES ROOT::MathCore ROOT::Hist) diff --git a/roottest/root/math/MathCoreTests.cxx b/roottest/root/math/MathCoreTests.cxx new file mode 100644 index 0000000000000..19e85e0eab039 --- /dev/null +++ b/roottest/root/math/MathCoreTests.cxx @@ -0,0 +1,77 @@ +#include +#include + +#include + +#include + +TEST(TMath, Gradient_Laplace) +{ + std::array parameters{2, 100., -1., -2., 0.1}; + std::array parameters1st{}; + std::array 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 vx; + std::array 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); + } +} \ No newline at end of file