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