src/gromacs/tables/tests/splinetable.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2015,2016,2017,2018,2019,2020,2021, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Tests for simple math functions.eval
  38  *
  39  * \author Erik Lindahl <erik.lindahl@gmail.com>
  40  * \ingroup module_tables
  41  */
  42 #include "gmxpre.h"
  43
  44 #include <cmath>
  45
  46 #include <algorithm>
  47 #include <functional>
  48 #include <utility>
  49
  50 #include <gtest/gtest.h>
  51
  52 #include "gromacs/math/units.h"
  53 #include "gromacs/math/utilities.h"
  54 #include "gromacs/options/basicoptions.h"
  55 #include "gromacs/options/ioptionscontainer.h"
  56 #include "gromacs/simd/simd.h"
  57 #include "gromacs/tables/cubicsplinetable.h"
  58 #include "gromacs/tables/quadraticsplinetable.h"
  59
  60 #include "testutils/testasserts.h"
  61 #include "testutils/testoptions.h"
  62
  63
  64 namespace gmx
  65 {
  66
  67 namespace test
  68 {
  69
  70 namespace
  71 {
  72
  73 class SplineTableTestBase : public ::testing::Test
  74 {
  75 public:
  76     static int s_testPoints_; //!< Number of points to use. Public so we can set it as option
  77 };
  78
  79 int SplineTableTestBase::s_testPoints_ = 100;
  80
  81 /*! \cond */
  82 /*! \brief Command-line option to adjust the number of points used to test SIMD math functions. */
  83 GMX_TEST_OPTIONS(SplineTableTestOptions, options)
  84 {
  85     options->addOption(::gmx::IntegerOption("npoints")
  86                                .store(&SplineTableTestBase::s_testPoints_)
  87                                .description("Number of points to test for spline table functions"));
  88 }
  89 /*! \endcond */
  90
  91
  92 /*! \brief Test fixture for table comparision with analytical/numerical functions */
  93 template<typename T>
  94 class SplineTableTest : public SplineTableTestBase
  95 {
  96 public:
  97     SplineTableTest() : tolerance_(T::defaultTolerance) {}
  98
  99     /*! \brief Set a new tolerance to be used in table function comparison
 100      *
 101      *  \param tol New tolerance to use
 102      */
 103     void setTolerance(real tol) { tolerance_ = tol; }
 104
 105     //! \cond internal
 106     /*! \internal \brief
 107      * Assertion predicate formatter for comparing table with function/derivative
 108      */
 109     template<int numFuncInTable = 1, int funcIndex = 0>
 110     void testSplineTableAgainstFunctions(const std::string&                   desc,
 111                                          const std::function<double(double)>& refFunc,
 112                                          const std::function<double(double)>& refDer,
 113                                          const T&                             table,
 114                                          const std::pair<real, real>&         testRange);
 115     //! \endcond
 116
 117 private:
 118     real tolerance_; //!< Tolerance to use
 119 };
 120
 121 template<class T>
 122 template<int numFuncInTable, int funcIndex>
 123 void SplineTableTest<T>::testSplineTableAgainstFunctions(const std::string&                   desc,
 124                                                          const std::function<double(double)>& refFunc,
 125                                                          const std::function<double(double)>& refDer,
 126                                                          const T&                             table,
 127                                                          const std::pair<real, real>& testRange)
 128 {
 129     real dx = (testRange.second - testRange.first) / s_testPoints_;
 130
 131     FloatingPointTolerance funcTolerance(relativeToleranceAsFloatingPoint(0.0, tolerance_));
 132
 133     for (real x = testRange.first; x < testRange.second; x += dx) // NOLINT(clang-analyzer-security.FloatLoopCounter)
 134     {
 135         real h                = std::sqrt(GMX_REAL_EPS);
 136         real secondDerivative = (refDer(x + h) - refDer(x)) / h;
 137
 138         real testFuncValue;
 139         real testDerValue;
 140
 141         table.template evaluateFunctionAndDerivative<numFuncInTable, funcIndex>(
 142                 x, &testFuncValue, &testDerValue);
 143
 144         // Check that we get the same values from function/derivative-only methods
 145         real tmpFunc, tmpDer;
 146
 147         table.template evaluateFunction<numFuncInTable, funcIndex>(x, &tmpFunc);
 148
 149         table.template evaluateDerivative<numFuncInTable, funcIndex>(x, &tmpDer);
 150
 151         // Before we even start to think about errors related to the table interpolation
 152         // accuracy, we want to test that the interpolations are consistent whether we
 153         // call the routine that evaluates both the function and derivative or only one
 154         // of them.
 155         // Note that for these tests the relevant tolerance is NOT the default one
 156         // provided based on the requested accuracy of the table, but a tolerance related
 157         // to the floating-point precision used. For now we only allow deviations up
 158         // to 4 ulp (one for the FMA order, and then some margin).
 159         FloatingPointTolerance consistencyTolerance(ulpTolerance(4));
 160
 161         FloatingPointDifference evaluateFuncDiff(tmpFunc, testFuncValue);
 162         if (!consistencyTolerance.isWithin(evaluateFuncDiff))
 163         {
 164             ADD_FAILURE() << "Interpolation inconsistency for table " << desc << std::endl
 165                           << numFuncInTable << " function(s) in table, testing index " << funcIndex
 166                           << std::endl
 167                           << "First failure at x = " << x << std::endl
 168                           << "Function value when evaluating function & derivative: " << testFuncValue
 169                           << std::endl
 170                           << "Function value when evaluating only function:         " << tmpFunc
 171                           << std::endl;
 172             return;
 173         }
 174
 175         FloatingPointDifference evaluateDerDiff(tmpDer, testDerValue);
 176         if (!consistencyTolerance.isWithin(evaluateDerDiff))
 177         {
 178             ADD_FAILURE() << "Interpolation inconsistency for table " << desc << std::endl
 179                           << numFuncInTable << " function(s) in table, testing index " << funcIndex
 180                           << std::endl
 181                           << "First failure at x = " << x << std::endl
 182                           << "Derivative value when evaluating function & derivative: " << testDerValue
 183                           << std::endl
 184                           << "Derivative value when evaluating only derivative:       " << tmpDer
 185                           << std::endl;
 186             return;
 187         }
 188
 189         // Next, we should examine that the table is exact enough relative
 190         // to the requested accuracy in the interpolation.
 191         //
 192         // There are two sources of errors that we need to account for when checking the values,
 193         // and we only fail the test if both of these tolerances are violated:
 194         //
 195         // 1) First, we have the normal relative error of the test vs. reference value. For this
 196         //    we use the normal testutils relative tolerance checking.
 197         //    However, there is an additional source of error: When we calculate the forces we
 198         //    use average higher derivatives over the interval to improve accuracy, but this
 199         //    also means we won't reproduce values at table points exactly. This is usually not
 200         //    an issue since the tolerances we have are much larger, but when the reference derivative
 201         //    value is exactly zero the relative error will be infinite. To account for this, we
 202         //    extract the spacing from the table and evaluate the reference derivative at a point
 203         //    this much larger too, and use the largest of the two values as the reference
 204         //    magnitude for the derivative when setting the relative tolerance.
 205         //    Note that according to the table function definitions, we should be allowed to evaluate
 206         //    it one table point beyond the range (this is done already for construction).
 207         //
 208         // 2) Second, due to the loss-of-accuracy when calculating the index through rtable
 209         //    there is an internal absolute tolerance that we can calculate.
 210         //    The source of this error is the subtraction eps=rtab-[rtab], which leaves an
 211         //    error proportional to eps_machine*rtab=eps_machine*x*tableScale.
 212         //    To lowest order, the term in the function and derivative values (respectively) that
 213         //    are proportional to eps will be the next-higher derivative multiplied by the spacing.
 214         //    This means the truncation error in the value is derivative*x*eps_machine, and in the
 215         //    derivative the error is 2nd_derivative*x*eps_machine.
 216
 217         real refFuncValue    = refFunc(x);
 218         real refDerValue     = refDer(x);
 219         real nextRefDerValue = refDer(x + table.tableSpacing());
 220
 221         real derMagnitude = std::max(std::abs(refDerValue), std::abs(nextRefDerValue));
 222
 223         // Since the reference magnitude will change over each interval we need to re-evaluate
 224         // the derivative tolerance inside the loop.
 225         FloatingPointTolerance derTolerance(relativeToleranceAsFloatingPoint(derMagnitude, tolerance_));
 226
 227         FloatingPointDifference funcDiff(refFuncValue, testFuncValue);
 228         FloatingPointDifference derDiff(refDerValue, testDerValue);
 229
 230         real allowedAbsFuncErr = std::abs(refDerValue) * x * GMX_REAL_EPS;
 231         real allowedAbsDerErr  = std::abs(secondDerivative) * x * GMX_REAL_EPS;
 232
 233         if ((!funcTolerance.isWithin(funcDiff) && funcDiff.asAbsolute() > allowedAbsFuncErr)
 234             || (!derTolerance.isWithin(derDiff) && derDiff.asAbsolute() > allowedAbsDerErr))
 235         {
 236             ADD_FAILURE() << "Failing comparison with function for table " << desc << std::endl
 237                           << numFuncInTable << " function(s) in table, testing index " << funcIndex
 238                           << std::endl
 239                           << "Test range is ( " << testRange.first << " , " << testRange.second
 240                           << " ) " << std::endl
 241                           << "Tolerance             = " << tolerance_ << std::endl
 242                           << "First failure at    x = " << x << std::endl
 243                           << "Reference function    = " << refFuncValue << std::endl
 244                           << "Test table function   = " << testFuncValue << std::endl
 245                           << "Allowed abs func err. = " << allowedAbsFuncErr << std::endl
 246                           << "Reference derivative  = " << refDerValue << std::endl
 247                           << "Test table derivative = " << testDerValue << std::endl
 248                           << "Allowed abs der. err. = " << allowedAbsDerErr << std::endl
 249                           << "Actual abs der. err.  = " << derDiff.asAbsolute() << std::endl;
 250             return;
 251         }
 252     }
 253 }
 254
 255
 256 /*! \brief Function similar to coulomb electrostatics
 257  *
 258  *  \param r argument
 259  *  \return r^-1
 260  */
 261 double coulombFunction(double r)
 262 {
 263     return 1.0 / r;
 264 }
 265
 266 /*! \brief Derivative (not force) of coulomb electrostatics
 267  *
 268  *  \param r argument
 269  *  \return -r^-2
 270  */
 271 double coulombDerivative(double r)
 272 {
 273     return -1.0 / (r * r);
 274 }
 275
 276 /*! \brief Function similar to power-6 Lennard-Jones dispersion
 277  *
 278  *  \param r argument
 279  *  \return r^-6
 280  */
 281 double lj6Function(double r)
 282 {
 283     return std::pow(r, -6.0);
 284 }
 285
 286 /*! \brief Derivative (not force) of the power-6 Lennard-Jones dispersion
 287  *
 288  *  \param r argument
 289  *  \return -6.0*r^-7
 290  */
 291 double lj6Derivative(double r)
 292 {
 293     return -6.0 * std::pow(r, -7.0);
 294 }
 295
 296 /*! \brief Function similar to power-12 Lennard-Jones repulsion
 297  *
 298  *  \param r argument
 299  *  \return r^-12
 300  */
 301 double lj12Function(double r)
 302 {
 303     return std::pow(r, -12.0);
 304 }
 305
 306 /*! \brief Derivative (not force) of the power-12 Lennard-Jones repulsion
 307  *
 308  *  \param r argument
 309  *  \return -12.0*r^-13
 310  */
 311 double lj12Derivative(double r)
 312 {
 313     return -12.0 * std::pow(r, -13.0);
 314 }
 315
 316 /*! \brief The sinc function, sin(r)/r
 317  *
 318  *  \param r argument
 319  *  \return sin(r)/r
 320  */
 321 double sincFunction(double r)
 322 {
 323     return std::sin(r) / r;
 324 }
 325
 326 /*! \brief Derivative of the sinc function
 327  *
 328  *  \param r argument
 329  *  \return derivative of sinc, (r*cos(r)-sin(r))/r^2
 330  */
 331 double sincDerivative(double r)
 332 {
 333     return (r * std::cos(r) - std::sin(r)) / (r * r);
 334 }
 335
 336 /*! \brief Function for the direct-space PME correction to 1/r
 337  *
 338  *  \param r argument
 339  *  \return PME correction function, erf(r)/r
 340  */
 341 double pmeCorrFunction(double r)
 342 {
 343     if (r == 0)
 344     {
 345         return 2.0 / std::sqrt(M_PI);
 346     }
 347     else
 348     {
 349         return std::erf(r) / r;
 350     }
 351 }
 352
 353 /*! \brief Derivative of the direct-space PME correction to 1/r
 354  *
 355  *  \param r argument
 356  *  \return Derivative of the PME correction function.
 357  */
 358 double pmeCorrDerivative(double r)
 359 {
 360     if (r == 0)
 361     {
 362         return 0;
 363     }
 364     else
 365     {
 366         return (2.0 * std::exp(-r * r) / std::sqrt(3.14159265358979323846) * r - erf(r)) / (r * r);
 367     }
 368 }
 369
 370 /*! \brief Typed-test list. We test QuadraticSplineTable and CubicSplineTable
 371  */
 372 typedef ::testing::Types<QuadraticSplineTable, CubicSplineTable> SplineTableTypes;
 373 TYPED_TEST_CASE(SplineTableTest, SplineTableTypes);
 374
 375
 376 TYPED_TEST(SplineTableTest, HandlesIncorrectInput)
 377 {
 378     // negative range
 379     EXPECT_THROW_GMX(TypeParam({ { "LJ12", lj12Function, lj12Derivative } }, { -1.0, 0.0 }),
 380                      gmx::InvalidInputError);
 381
 382     // Too small range
 383     EXPECT_THROW_GMX(TypeParam({ { "LJ12", lj12Function, lj12Derivative } }, { 1.0, 1.00001 }),
 384                      gmx::InvalidInputError);
 385
 386     // bad tolerance
 387     EXPECT_THROW_GMX(TypeParam({ { "LJ12", lj12Function, lj12Derivative } }, { 1.0, 2.0 }, 1e-20),
 388                      gmx::ToleranceError);
 389
 390     // Range is so close to 0.0 that table would require >1e6 points
 391     EXPECT_THROW_GMX(TypeParam({ { "LJ12", lj12Function, lj12Derivative } }, { 1e-4, 2.0 }),
 392                      gmx::ToleranceError);
 393
 394     // mismatching function/derivative
 395     EXPECT_THROW_GMX(TypeParam({ { "BadLJ12", lj12Derivative, lj12Function } }, { 1.0, 2.0 }),
 396                      gmx::InconsistentInputError);
 397 }
 398
 399
 400 #ifndef NDEBUG
 401 TYPED_TEST(SplineTableTest, CatchesOutOfRangeValues)
 402 {
 403     TypeParam table({ { "LJ12", lj12Function, lj12Derivative } }, { 0.2, 1.0 });
 404     real      x, func, der;
 405
 406     x = -GMX_REAL_EPS;
 407     EXPECT_THROW_GMX(table.evaluateFunctionAndDerivative(x, &func, &der), gmx::RangeError);
 408
 409     x = 1.0;
 410     EXPECT_THROW_GMX(table.evaluateFunctionAndDerivative(x, &func, &der), gmx::RangeError);
 411 }
 412 #endif
 413
 414
 415 TYPED_TEST(SplineTableTest, Sinc)
 416 {
 417     // Sinc hits some sensitive parts of the table construction code which means
 418     // we will not have full relative accuracy close to the zeros in the
 419     // derivative. Since this is intentially a pathological function we reduce
 420     // the interval slightly for now.
 421     std::pair<real, real> range(0.1, 3.1);
 422
 423     TypeParam sincTable({ { "Sinc", sincFunction, sincDerivative } }, range);
 424
 425     TestFixture::testSplineTableAgainstFunctions("Sinc", sincFunction, sincDerivative, sincTable, range);
 426 }
 427
 428
 429 TYPED_TEST(SplineTableTest, LJ12)
 430 {
 431     std::pair<real, real> range(0.2, 2.0);
 432
 433     TypeParam lj12Table({ { "LJ12", lj12Function, lj12Derivative } }, range);
 434
 435     TestFixture::testSplineTableAgainstFunctions("LJ12", lj12Function, lj12Derivative, lj12Table, range);
 436 }
 437
 438
 439 TYPED_TEST(SplineTableTest, PmeCorrection)
 440 {
 441     std::pair<real, real> range(0.0, 4.0);
 442     real                  tolerance = 1e-5;
 443
 444     TypeParam pmeCorrTable({ { "PMECorr", pmeCorrFunction, pmeCorrDerivative } }, range, tolerance);
 445
 446     TestFixture::setTolerance(tolerance);
 447     TestFixture::testSplineTableAgainstFunctions(
 448             "PMECorr", pmeCorrFunction, pmeCorrDerivative, pmeCorrTable, range);
 449 }
 450
 451
 452 TYPED_TEST(SplineTableTest, HandlesIncorrectNumericalInput)
 453 {
 454     // Lengths do not match
 455     std::vector<double> functionValues(10);
 456     std::vector<double> derivativeValues(20);
 457     EXPECT_THROW_GMX(
 458             TypeParam({ { "EmptyVectors", functionValues, derivativeValues, 0.001 } }, { 1.0, 2.0 }),
 459             gmx::InconsistentInputError);
 460
 461     // Upper range is 2.0, spacing 0.1. This requires at least 21 points. Make sure we get an error for 20.
 462     functionValues.resize(20);
 463     derivativeValues.resize(20);
 464     EXPECT_THROW_GMX(TypeParam({ { "EmptyVectors", functionValues, derivativeValues, 0.1 } }, { 1.0, 2.0 }),
 465                      gmx::InconsistentInputError);
 466
 467     // Create some test data
 468     functionValues.clear();
 469     derivativeValues.clear();
 470
 471     std::vector<double> badDerivativeValues;
 472     double              spacing = 1e-3;
 473
 474     for (std::size_t i = 0; i < 1001; i++)
 475     {
 476         double x    = i * spacing;
 477         double func = (x >= 0.1) ? lj12Function(x) : 0.0;
 478         double der  = (x >= 0.1) ? lj12Derivative(x) : 0.0;
 479
 480         functionValues.push_back(func);
 481         derivativeValues.push_back(der);
 482         badDerivativeValues.push_back(-der);
 483     }
 484
 485     // Derivatives not consistent with function
 486     EXPECT_THROW_GMX(TypeParam({ { "NumericalBadLJ12", functionValues, badDerivativeValues, spacing } },
 487                                { 0.2, 1.0 }),
 488                      gmx::InconsistentInputError);
 489
 490     // Spacing 1e-3 is not sufficient for r^-12 in range [0.1,1.0]
 491     // Make sure we get a tolerance error
 492     EXPECT_THROW_GMX(
 493             TypeParam({ { "NumericalLJ12", functionValues, derivativeValues, spacing } }, { 0.2, 1.0 }),
 494             gmx::ToleranceError);
 495 }
 496
 497
 498 TYPED_TEST(SplineTableTest, NumericalInputPmeCorr)
 499 {
 500     std::pair<real, real> range(0.0, 4.0);
 501     std::vector<double>   functionValues;
 502     std::vector<double>   derivativeValues;
 503
 504     double inputSpacing = 1e-3;
 505     real   tolerance    = 1e-5;
 506
 507     // We only need data up to the argument 4.0, but add 1% margin
 508     for (std::size_t i = 0; i < range.second * 1.01 / inputSpacing; i++)
 509     {
 510         double x = i * inputSpacing;
 511
 512         functionValues.push_back(pmeCorrFunction(x));
 513         derivativeValues.push_back(pmeCorrDerivative(x));
 514     }
 515
 516     TypeParam pmeCorrTable(
 517             { { "NumericalPMECorr", functionValues, derivativeValues, inputSpacing } }, range, tolerance);
 518
 519     TestFixture::setTolerance(tolerance);
 520     TestFixture::testSplineTableAgainstFunctions(
 521             "NumericalPMECorr", pmeCorrFunction, pmeCorrDerivative, pmeCorrTable, range);
 522 }
 523
 524 TYPED_TEST(SplineTableTest, TwoFunctions)
 525 {
 526     std::pair<real, real> range(0.2, 2.0);
 527
 528     TypeParam table({ { "LJ6", lj6Function, lj6Derivative }, { "LJ12", lj12Function, lj12Derivative } },
 529                     range);
 530
 531     // Test entire range for each function. This will use the method that interpolates a single function
 532     TestFixture::template testSplineTableAgainstFunctions<2, 0>(
 533             "LJ6", lj6Function, lj6Derivative, table, range);
 534     TestFixture::template testSplineTableAgainstFunctions<2, 1>(
 535             "LJ12", lj12Function, lj12Derivative, table, range);
 536
 537     // Test the methods that evaluated both functions for one value
 538     real x        = 0.5 * (range.first + range.second);
 539     real refFunc0 = lj6Function(x);
 540     real refDer0  = lj6Derivative(x);
 541     real refFunc1 = lj12Function(x);
 542     real refDer1  = lj12Derivative(x);
 543
 544     real tstFunc0, tstDer0, tstFunc1, tstDer1;
 545     real tmpFunc0, tmpFunc1, tmpDer0, tmpDer1;
 546
 547     // test that we reproduce the reference functions
 548     table.evaluateFunctionAndDerivative(x, &tstFunc0, &tstDer0, &tstFunc1, &tstDer1);
 549
 550     real funcErr0 = std::abs(tstFunc0 - refFunc0) / std::abs(refFunc0);
 551     real funcErr1 = std::abs(tstFunc1 - refFunc1) / std::abs(refFunc1);
 552     real derErr0  = std::abs(tstDer0 - refDer0) / std::abs(refDer0);
 553     real derErr1  = std::abs(tstDer1 - refDer1) / std::abs(refDer1);
 554
 555     // Use asserts, since the following ones compare to these values.
 556     ASSERT_LT(funcErr0, TypeParam::defaultTolerance);
 557     ASSERT_LT(derErr0, TypeParam::defaultTolerance);
 558     ASSERT_LT(funcErr1, TypeParam::defaultTolerance);
 559     ASSERT_LT(derErr1, TypeParam::defaultTolerance);
 560
 561     // Test that function/derivative-only interpolation methods work
 562     table.evaluateFunction(x, &tmpFunc0, &tmpFunc1);
 563     table.evaluateDerivative(x, &tmpDer0, &tmpDer1);
 564     EXPECT_EQ(tstFunc0, tmpFunc0);
 565     EXPECT_EQ(tstFunc1, tmpFunc1);
 566     EXPECT_EQ(tstDer0, tmpDer0);
 567
 568     // Test that scrambled order interpolation methods work
 569     table.template evaluateFunctionAndDerivative<2, 1, 0>(x, &tstFunc1, &tstDer1, &tstFunc0, &tstDer0);
 570     EXPECT_EQ(tstFunc0, tmpFunc0);
 571     EXPECT_EQ(tstFunc1, tmpFunc1);
 572     EXPECT_EQ(tstDer0, tmpDer0);
 573     EXPECT_EQ(tstDer1, tmpDer1);
 574
 575     // Test scrambled order for function/derivative-only methods
 576     table.template evaluateFunction<2, 1, 0>(x, &tmpFunc1, &tmpFunc0);
 577     table.template evaluateDerivative<2, 1, 0>(x, &tmpDer1, &tmpDer0);
 578     EXPECT_EQ(tstFunc0, tmpFunc0);
 579     EXPECT_EQ(tstFunc1, tmpFunc1);
 580     EXPECT_EQ(tstDer0, tmpDer0);
 581     EXPECT_EQ(tstDer1, tmpDer1);
 582 }
 583
 584 TYPED_TEST(SplineTableTest, ThreeFunctions)
 585 {
 586     std::pair<real, real> range(0.2, 2.0);
 587
 588     TypeParam table({ { "Coulomb", coulombFunction, coulombDerivative },
 589                       { "LJ6", lj6Function, lj6Derivative },
 590                       { "LJ12", lj12Function, lj12Derivative } },
 591                     range);
 592
 593     // Test entire range for each function
 594     TestFixture::template testSplineTableAgainstFunctions<3, 0>(
 595             "Coulomb", coulombFunction, coulombDerivative, table, range);
 596     TestFixture::template testSplineTableAgainstFunctions<3, 1>(
 597             "LJ6", lj6Function, lj6Derivative, table, range);
 598     TestFixture::template testSplineTableAgainstFunctions<3, 2>(
 599             "LJ12", lj12Function, lj12Derivative, table, range);
 600
 601     // Test the methods that evaluated both functions for one value
 602     real x        = 0.5 * (range.first + range.second);
 603     real refFunc0 = coulombFunction(x);
 604     real refDer0  = coulombDerivative(x);
 605     real refFunc1 = lj6Function(x);
 606     real refDer1  = lj6Derivative(x);
 607     real refFunc2 = lj12Function(x);
 608     real refDer2  = lj12Derivative(x);
 609
 610     real tstFunc0, tstDer0, tstFunc1, tstDer1, tstFunc2, tstDer2;
 611     real tmpFunc0, tmpFunc1, tmpFunc2, tmpDer0, tmpDer1, tmpDer2;
 612
 613     // test that we reproduce the reference functions
 614     table.evaluateFunctionAndDerivative(x, &tstFunc0, &tstDer0, &tstFunc1, &tstDer1, &tstFunc2, &tstDer2);
 615
 616     real funcErr0 = std::abs(tstFunc0 - refFunc0) / std::abs(refFunc0);
 617     real derErr0  = std::abs(tstDer0 - refDer0) / std::abs(refDer0);
 618     real funcErr1 = std::abs(tstFunc1 - refFunc1) / std::abs(refFunc1);
 619     real derErr1  = std::abs(tstDer1 - refDer1) / std::abs(refDer1);
 620     real funcErr2 = std::abs(tstFunc2 - refFunc2) / std::abs(refFunc2);
 621     real derErr2  = std::abs(tstDer2 - refDer2) / std::abs(refDer2);
 622
 623     // Use asserts, since the following ones compare to these values.
 624     ASSERT_LT(funcErr0, TypeParam::defaultTolerance);
 625     ASSERT_LT(derErr0, TypeParam::defaultTolerance);
 626     ASSERT_LT(funcErr1, TypeParam::defaultTolerance);
 627     ASSERT_LT(derErr1, TypeParam::defaultTolerance);
 628     ASSERT_LT(funcErr2, TypeParam::defaultTolerance);
 629     ASSERT_LT(derErr2, TypeParam::defaultTolerance);
 630
 631     // Test that function/derivative-only interpolation methods work
 632     table.evaluateFunction(x, &tmpFunc0, &tmpFunc1, &tmpFunc2);
 633     table.evaluateDerivative(x, &tmpDer0, &tmpDer1, &tmpDer2);
 634     EXPECT_EQ(tstFunc0, tmpFunc0);
 635     EXPECT_EQ(tstFunc1, tmpFunc1);
 636     EXPECT_EQ(tstFunc2, tmpFunc2);
 637     EXPECT_EQ(tstDer0, tmpDer0);
 638     EXPECT_EQ(tstDer1, tmpDer1);
 639     EXPECT_EQ(tstDer2, tmpDer2);
 640
 641     // Test two functions out of three
 642     table.template evaluateFunctionAndDerivative<3, 0, 1>(x, &tmpFunc0, &tmpDer0, &tmpFunc1, &tmpDer1);
 643     EXPECT_EQ(tstFunc0, tmpFunc0);
 644     EXPECT_EQ(tstFunc1, tmpFunc1);
 645     EXPECT_EQ(tstDer0, tmpDer0);
 646     EXPECT_EQ(tstDer1, tmpDer1);
 647
 648     // two out of three, function/derivative-only
 649     table.template evaluateFunction<3, 0, 1>(x, &tmpFunc0, &tmpFunc1);
 650     table.template evaluateDerivative<3, 0, 1>(x, &tmpDer0, &tmpDer1);
 651     EXPECT_EQ(tstFunc0, tmpFunc0);
 652     EXPECT_EQ(tstFunc1, tmpFunc1);
 653     EXPECT_EQ(tstDer0, tmpDer0);
 654     EXPECT_EQ(tstDer1, tmpDer1);
 655
 656     // Test that scrambled order interpolation methods work
 657     table.template evaluateFunctionAndDerivative<3, 2, 1, 0>(
 658             x, &tstFunc2, &tstDer2, &tstFunc1, &tstDer1, &tstFunc0, &tstDer0);
 659     EXPECT_EQ(tstFunc0, tmpFunc0);
 660     EXPECT_EQ(tstFunc1, tmpFunc1);
 661     EXPECT_EQ(tstFunc2, tmpFunc2);
 662     EXPECT_EQ(tstDer0, tmpDer0);
 663     EXPECT_EQ(tstDer1, tmpDer1);
 664     EXPECT_EQ(tstDer2, tmpDer2);
 665
 666     // Test scrambled order for function/derivative-only methods
 667     table.template evaluateFunction<3, 2, 1, 0>(x, &tmpFunc2, &tmpFunc1, &tmpFunc0);
 668     table.template evaluateDerivative<3, 2, 1, 0>(x, &tmpDer2, &tmpDer1, &tmpDer0);
 669     EXPECT_EQ(tstFunc0, tmpFunc0);
 670     EXPECT_EQ(tstFunc1, tmpFunc1);
 671     EXPECT_EQ(tstFunc2, tmpFunc2);
 672     EXPECT_EQ(tstDer0, tmpDer0);
 673     EXPECT_EQ(tstDer1, tmpDer1);
 674     EXPECT_EQ(tstDer2, tmpDer2);
 675 }
 676
 677 #if GMX_SIMD_HAVE_REAL
 678 TYPED_TEST(SplineTableTest, Simd)
 679 {
 680     std::pair<real, real> range(0.2, 1.0);
 681     TypeParam             table({ { "LJ12", lj12Function, lj12Derivative } }, range);
 682
 683     // We already test that the SIMD operations handle the different elements
 684     // correctly in the SIMD module, so here we only test that interpolation
 685     // works for a single value in the middle of the interval
 686
 687     real                             x       = 0.5 * (range.first + range.second);
 688     real                             refFunc = lj12Function(x);
 689     real                             refDer  = lj12Derivative(x);
 690     SimdReal                         tstFunc, tstDer;
 691     real                             funcErr, derErr;
 692     alignas(GMX_SIMD_ALIGNMENT) real alignedMem[GMX_SIMD_REAL_WIDTH];
 693
 694     table.evaluateFunctionAndDerivative(SimdReal(x), &tstFunc, &tstDer);
 695
 696     store(alignedMem, tstFunc);
 697     funcErr = std::abs(alignedMem[0] - refFunc) / std::abs(refFunc);
 698
 699     store(alignedMem, tstDer);
 700     derErr = std::abs(alignedMem[0] - refDer) / std::abs(refDer);
 701
 702     EXPECT_LT(funcErr, TypeParam::defaultTolerance);
 703     EXPECT_LT(derErr, TypeParam::defaultTolerance);
 704 }
 705
 706 TYPED_TEST(SplineTableTest, SimdTwoFunctions)
 707 {
 708     std::pair<real, real> range(0.2, 2.0);
 709
 710     TypeParam table({ { "LJ6", lj6Function, lj6Derivative }, { "LJ12", lj12Function, lj12Derivative } },
 711                     range);
 712
 713     // We already test that the SIMD operations handle the different elements
 714     // correctly in the SIMD module, so here we only test that interpolation
 715     // works for a single value in the middle of the interval
 716
 717     real                             x        = 0.5 * (range.first + range.second);
 718     real                             refFunc0 = lj6Function(x);
 719     real                             refDer0  = lj6Derivative(x);
 720     real                             refFunc1 = lj12Function(x);
 721     real                             refDer1  = lj12Derivative(x);
 722     SimdReal                         tstFunc0, tstDer0;
 723     SimdReal                         tstFunc1, tstDer1;
 724     real                             funcErr0, derErr0;
 725     real                             funcErr1, derErr1;
 726     alignas(GMX_SIMD_ALIGNMENT) real alignedMem[GMX_SIMD_REAL_WIDTH];
 727
 728     table.evaluateFunctionAndDerivative(SimdReal(x), &tstFunc0, &tstDer0, &tstFunc1, &tstDer1);
 729
 730     store(alignedMem, tstFunc0);
 731     funcErr0 = std::abs(alignedMem[0] - refFunc0) / std::abs(refFunc0);
 732
 733     store(alignedMem, tstDer0);
 734     derErr0 = std::abs(alignedMem[0] - refDer0) / std::abs(refDer0);
 735
 736     store(alignedMem, tstFunc1);
 737     funcErr1 = std::abs(alignedMem[0] - refFunc1) / std::abs(refFunc1);
 738
 739     store(alignedMem, tstDer1);
 740     derErr1 = std::abs(alignedMem[0] - refDer1) / std::abs(refDer1);
 741
 742     EXPECT_LT(funcErr0, TypeParam::defaultTolerance);
 743     EXPECT_LT(derErr0, TypeParam::defaultTolerance);
 744     EXPECT_LT(funcErr1, TypeParam::defaultTolerance);
 745     EXPECT_LT(derErr1, TypeParam::defaultTolerance);
 746 }
 747 #endif
 748
 749 #if GMX_SIMD_HAVE_REAL && !defined NDEBUG
 750 TYPED_TEST(SplineTableTest, CatchesOutOfRangeValuesSimd)
 751 {
 752     std::pair<real, real> range(0.2, 1.0);
 753     TypeParam             table({ { "LJ12", lj12Function, lj12Derivative } }, range);
 754     SimdReal              x, func, der;
 755
 756     AlignedArray<real, GMX_SIMD_REAL_WIDTH> alignedMem;
 757
 758     alignedMem.fill(range.first);
 759     // Make position 1 incorrect if width>=2, otherwise position 0
 760     // range.first-GMX_REAL_EPS is not invalid. See comment in table.
 761     alignedMem[(GMX_SIMD_REAL_WIDTH >= 2) ? 1 : 0] = -GMX_REAL_EPS;
 762     x                                              = load<SimdReal>(alignedMem);
 763
 764     EXPECT_THROW_GMX(table.evaluateFunctionAndDerivative(x, &func, &der), gmx::RangeError);
 765
 766     // Make position 1 incorrect if width>=2, otherwise position 0
 767     alignedMem[(GMX_SIMD_REAL_WIDTH >= 2) ? 1 : 0] = range.second;
 768     x                                              = load<SimdReal>(alignedMem);
 769
 770     EXPECT_THROW_GMX(table.evaluateFunctionAndDerivative(x, &func, &der), gmx::RangeError);
 771 }
 772
 773 TYPED_TEST(SplineTableTest, AcceptsInRangeValuesSimd)
 774 {
 775     std::pair<real, real> range(0.2, 1.0);
 776     TypeParam             table({ { "LJ12", lj12Function, lj12Derivative } }, range);
 777     SimdReal              x, func, der;
 778
 779     alignas(GMX_SIMD_ALIGNMENT) real alignedMem[GMX_SIMD_REAL_WIDTH];
 780
 781     // Test all values between 0 and range.second
 782     for (std::size_t i = 0; i < GMX_SIMD_REAL_WIDTH; i++)
 783     {
 784         alignedMem[i] = range.second * (1.0 - GMX_REAL_EPS) * i / (GMX_SIMD_REAL_WIDTH - 1);
 785     }
 786     x = load<SimdReal>(alignedMem);
 787
 788     EXPECT_NO_THROW_GMX(table.evaluateFunctionAndDerivative(x, &func, &der));
 789 }
 790
 791 #endif
 792
 793 } // namespace
 794
 795 } // namespace test
 796
 797 } // namespace gmx