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, 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
  79 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
  93
  94 /*! \brief Test fixture for table comparision with analytical/numerical functions */
  95 template <typename T>
  96 class SplineTableTest : public SplineTableTestBase
  97 {
  98     public:
  99         SplineTableTest() : tolerance_(T::defaultTolerance) {}
 100
 101         /*! \brief Set a new tolerance to be used in table function comparison
 102          *
 103          *  \param tol New tolerance to use
 104          */
 105         void
 106         setTolerance(real tol) { tolerance_ = tol; }
 107
 108         //! \cond internal
 109         /*! \internal \brief
 110          * Assertion predicate formatter for comparing table with function/derivative
 111          */
 112         template<int numFuncInTable = 1, int funcIndex = 0>
 113         void
 114         testSplineTableAgainstFunctions(const std::string                    &desc,
 115                                         const std::function<double(double)>  &refFunc,
 116                                         const std::function<double(double)>  &refDer,
 117                                         const T                              &table,
 118                                         const std::pair<real, real>          &testRange);
 119         //! \endcond
 120
 121     private:
 122         real        tolerance_;    //!< Tolerance to use
 123 };
 124
 125 template <class T>
 126 template<int numFuncInTable, int funcIndex>
 127 void
 128 SplineTableTest<T>::testSplineTableAgainstFunctions(const std::string                    &desc,
 129                                                     const std::function<double(double)>  &refFunc,
 130                                                     const std::function<double(double)>  &refDer,
 131                                                     const T                              &table,
 132                                                     const std::pair<real, real>          &testRange)
 133 {
 134     real                   dx = (testRange.second - testRange.first) / s_testPoints_;
 135
 136     FloatingPointTolerance funcTolerance(relativeToleranceAsFloatingPoint(0.0, tolerance_));
 137
 138     for (real x = testRange.first; x < testRange.second; x += dx)
 139     {
 140         real h                = std::sqrt(GMX_REAL_EPS);
 141         real secondDerivative = (refDer(x+h)-refDer(x))/h;
 142
 143         real testFuncValue;
 144         real testDerValue;
 145
 146         table.template evaluateFunctionAndDerivative<numFuncInTable, funcIndex>(x, &testFuncValue, &testDerValue);
 147
 148         // Check that we get the same values from function/derivative-only methods
 149         real tmpFunc, tmpDer;
 150
 151         table.template evaluateFunction<numFuncInTable, funcIndex>(x, &tmpFunc);
 152
 153         table.template evaluateDerivative<numFuncInTable, funcIndex>(x, &tmpDer);
 154
 155         // Before we even start to think about errors related to the table interpolation
 156         // accuracy, we want to test that the interpolations are consistent whether we
 157         // call the routine that evaluates both the function and derivative or only one
 158         // of them.
 159         // Note that for these tests the relevant tolerance is NOT the default one
 160         // provided based on the requested accuracy of the table, but a tolerance related
 161         // to the floating-point precision used. For now we only allow deviations up
 162         // to 4 ulp (one for the FMA order, and then some margin).
 163         FloatingPointTolerance  consistencyTolerance(ulpTolerance(4));
 164
 165         FloatingPointDifference evaluateFuncDiff(tmpFunc, testFuncValue);
 166         if (!consistencyTolerance.isWithin(evaluateFuncDiff))
 167         {
 168             ADD_FAILURE()
 169             << "Interpolation inconsistency for table " << desc << std::endl
 170             << numFuncInTable << " function(s) in table, testing index " << funcIndex << std::endl
 171             << "First failure at x = " << x << std::endl
 172             << "Function value when evaluating function & derivative: " << testFuncValue << std::endl
 173             << "Function value when evaluating only function:         " << tmpFunc << std::endl;
 174             return;
 175         }
 176
 177         FloatingPointDifference evaluateDerDiff(tmpDer, testDerValue);
 178         if (!consistencyTolerance.isWithin(evaluateDerDiff))
 179         {
 180             ADD_FAILURE()
 181             << "Interpolation inconsistency for table " << desc << std::endl
 182             << numFuncInTable << " function(s) in table, testing index " << funcIndex << std::endl
 183             << "First failure at x = " << x << std::endl
 184             << "Derivative value when evaluating function & derivative: " << testDerValue << std::endl
 185             << "Derivative value when evaluating only derivative:       " << tmpDer << 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()
 237             << "Failing comparison with function for table " << desc << std::endl
 238             << numFuncInTable << " function(s) in table, testing index " << funcIndex << std::endl
 239             << "Test range is ( " << testRange.first << " , " << testRange.second << " ) " << 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
 261 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
 272 coulombDerivative(double r)
 273 {
 274     return -1.0/(r*r);
 275 }
 276
 277 /*! \brief Function similar to power-6 Lennard-Jones dispersion
 278  *
 279  *  \param r argument
 280  *  \return r^-6
 281  */
 282 double
 283 lj6Function(double r)
 284 {
 285     return std::pow(r, -6.0);
 286 }
 287
 288 /*! \brief Derivative (not force) of the power-6 Lennard-Jones dispersion
 289  *
 290  *  \param r argument
 291  *  \return -6.0*r^-7
 292  */
 293 double
 294 lj6Derivative(double r)
 295 {
 296     return -6.0*std::pow(r, -7.0);
 297 }
 298
 299 /*! \brief Function similar to power-12 Lennard-Jones repulsion
 300  *
 301  *  \param r argument
 302  *  \return r^-12
 303  */
 304 double
 305 lj12Function(double r)
 306 {
 307     return std::pow(r, -12.0);
 308 }
 309
 310 /*! \brief Derivative (not force) of the power-12 Lennard-Jones repulsion
 311  *
 312  *  \param r argument
 313  *  \return -12.0*r^-13
 314  */
 315 double
 316 lj12Derivative(double r)
 317 {
 318     return -12.0*std::pow(r, -13.0);
 319 }
 320
 321 /*! \brief The sinc function, sin(r)/r
 322  *
 323  *  \param r argument
 324  *  \return sin(r)/r
 325  */
 326 double
 327 sincFunction(double r)
 328 {
 329     return std::sin(r)/r;
 330 }
 331
 332 /*! \brief Derivative of the sinc function
 333  *
 334  *  \param r argument
 335  *  \return derivative of sinc, (r*cos(r)-sin(r))/r^2
 336  */
 337 double
 338 sincDerivative(double r)
 339 {
 340     return (r*std::cos(r)-std::sin(r))/(r*r);
 341 }
 342
 343 /*! \brief Function for the direct-space PME correction to 1/r
 344  *
 345  *  \param r argument
 346  *  \return PME correction function, erf(r)/r
 347  */
 348 double
 349 pmeCorrFunction(double r)
 350 {
 351     if (r == 0)
 352     {
 353         return 2.0/std::sqrt(M_PI);
 354     }
 355     else
 356     {
 357         return std::erf(r)/r;
 358     }
 359 }
 360
 361 /*! \brief Derivative of the direct-space PME correction to 1/r
 362  *
 363  *  \param r argument
 364  *  \return Derivative of the PME correction function.
 365  */
 366 double
 367 pmeCorrDerivative(double r)
 368 {
 369     if (r == 0)
 370     {
 371         return 0;
 372     }
 373     else
 374     {
 375         return (2.0*std::exp(-r*r)/std::sqrt(3.14159265358979323846)*r-erf(r))/(r*r);
 376     }
 377 }
 378
 379 /*! \brief Typed-test list. We test QuadraticSplineTable and CubicSplineTable
 380  */
 381 typedef ::testing::Types<QuadraticSplineTable, CubicSplineTable> SplineTableTypes;
 382 TYPED_TEST_CASE(SplineTableTest, SplineTableTypes);
 383
 384
 385 TYPED_TEST(SplineTableTest, HandlesIncorrectInput)
 386 {
 387     // negative range
 388     EXPECT_THROW_GMX(TypeParam( {{"LJ12", lj12Function, lj12Derivative}}, {-1.0, 0.0}), gmx::InvalidInputError);
 389
 390     // Too small range
 391     EXPECT_THROW_GMX(TypeParam( {{"LJ12", lj12Function, lj12Derivative}}, {1.0, 1.00001}), gmx::InvalidInputError);
 392
 393     // bad tolerance
 394     EXPECT_THROW_GMX(TypeParam( {{"LJ12", lj12Function, lj12Derivative}}, {1.0, 2.0}, 1e-20), gmx::ToleranceError);
 395
 396     // Range is so close to 0.0 that table would require >1e6 points
 397     EXPECT_THROW_GMX(TypeParam( {{"LJ12", lj12Function, lj12Derivative}}, {1e-4, 2.0}), gmx::ToleranceError);
 398
 399     // mismatching function/derivative
 400     EXPECT_THROW_GMX(TypeParam( { {"BadLJ12", lj12Derivative, lj12Function}}, {1.0, 2.0}), gmx::InconsistentInputError);
 401 }
 402
 403
 404 #ifndef NDEBUG
 405 TYPED_TEST(SplineTableTest, CatchesOutOfRangeValues)
 406 {
 407     TypeParam      table( {{"LJ12", lj12Function, lj12Derivative}}, {0.2, 1.0});
 408     real           x, func, der;
 409
 410     x = -GMX_REAL_EPS;
 411     EXPECT_THROW_GMX(table.evaluateFunctionAndDerivative(x, &func, &der), gmx::RangeError);
 412
 413     x = 1.0;
 414     EXPECT_THROW_GMX(table.evaluateFunctionAndDerivative(x, &func, &der), gmx::RangeError);
 415 }
 416 #endif
 417
 418
 419 TYPED_TEST(SplineTableTest, Sinc)
 420 {
 421     // Sinc hits some sensitive parts of the table construction code which means
 422     // we will not have full relative accuracy close to the zeros in the
 423     // derivative. Since this is intentially a pathological function we reduce
 424     // the interval slightly for now.
 425     std::pair<real, real>  range(0.1, 3.1);
 426
 427     TypeParam              sincTable( {{"Sinc", sincFunction, sincDerivative}}, range);
 428
 429     TestFixture::testSplineTableAgainstFunctions("Sinc", sincFunction, sincDerivative, sincTable, range);
 430 }
 431
 432
 433 TYPED_TEST(SplineTableTest, LJ12)
 434 {
 435     std::pair<real, real>  range(0.2, 2.0);
 436
 437     TypeParam              lj12Table( {{"LJ12", lj12Function, lj12Derivative}}, range);
 438
 439     TestFixture::testSplineTableAgainstFunctions("LJ12", lj12Function, lj12Derivative, lj12Table, range);
 440 }
 441
 442
 443 TYPED_TEST(SplineTableTest, PmeCorrection)
 444 {
 445     std::pair<real, real>  range(0.0, 4.0);
 446     real                   tolerance = 1e-5;
 447
 448     TypeParam              pmeCorrTable( {{"PMECorr", pmeCorrFunction, pmeCorrDerivative}}, range, tolerance);
 449
 450     TestFixture::setTolerance(tolerance);
 451     TestFixture::testSplineTableAgainstFunctions("PMECorr", pmeCorrFunction, pmeCorrDerivative, pmeCorrTable, range);
 452 }
 453
 454
 455
 456 TYPED_TEST(SplineTableTest, HandlesIncorrectNumericalInput)
 457 {
 458     // Lengths do not match
 459     std::vector<double>   functionValues(10);
 460     std::vector<double>   derivativeValues(20);
 461     EXPECT_THROW_GMX(TypeParam( {{"EmptyVectors", functionValues, derivativeValues, 0.001}},
 462                                 {1.0, 2.0}), gmx::InconsistentInputError);
 463
 464     // Upper range is 2.0, spacing 0.1. This requires at least 21 points. Make sure we get an error for 20.
 465     functionValues.resize(20);
 466     derivativeValues.resize(20);
 467     EXPECT_THROW_GMX(TypeParam( {{"EmptyVectors", functionValues, derivativeValues, 0.1}},
 468                                 {1.0, 2.0}), gmx::InconsistentInputError);
 469
 470     // Create some test data
 471     functionValues.clear();
 472     derivativeValues.clear();
 473
 474     std::vector<double>   badDerivativeValues;
 475     double                spacing = 1e-3;
 476
 477     for (std::size_t i = 0; i < 1001; i++)
 478     {
 479         double x    = i * spacing;
 480         double func = (x >= 0.1) ? lj12Function(x)   : 0.0;
 481         double der  = (x >= 0.1) ? lj12Derivative(x) : 0.0;
 482
 483         functionValues.push_back(func);
 484         derivativeValues.push_back(der);
 485         badDerivativeValues.push_back(-der);
 486     }
 487
 488     // Derivatives not consistent with function
 489     EXPECT_THROW_GMX(TypeParam( {{"NumericalBadLJ12", functionValues, badDerivativeValues, spacing}},
 490                                 {0.2, 1.0}), gmx::InconsistentInputError);
 491
 492     // Spacing 1e-3 is not sufficient for r^-12 in range [0.1,1.0]
 493     // Make sure we get a tolerance error
 494     EXPECT_THROW_GMX(TypeParam( {{"NumericalLJ12", functionValues, derivativeValues, spacing}},
 495                                 {0.2, 1.0}), gmx::ToleranceError);
 496 }
 497
 498
 499 TYPED_TEST(SplineTableTest, NumericalInputPmeCorr)
 500 {
 501     std::pair<real, real>  range(0.0, 4.0);
 502     std::vector<double>    functionValues;
 503     std::vector<double>    derivativeValues;
 504
 505     double                 inputSpacing = 1e-3;
 506     real                   tolerance    = 1e-5;
 507
 508     // We only need data up to the argument 4.0, but add 1% margin
 509     for (std::size_t i = 0; i < range.second*1.01/inputSpacing; i++)
 510     {
 511         double x    = i * inputSpacing;
 512
 513         functionValues.push_back(pmeCorrFunction(x));
 514         derivativeValues.push_back(pmeCorrDerivative(x));
 515     }
 516
 517     TypeParam  pmeCorrTable( {{"NumericalPMECorr", functionValues, derivativeValues, inputSpacing}},
 518                              range, tolerance);
 519
 520     TestFixture::setTolerance(tolerance);
 521     TestFixture::testSplineTableAgainstFunctions("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}}, 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, table, range);
 532     TestFixture::template testSplineTableAgainstFunctions<2, 1>("LJ12", lj12Function, lj12Derivative, table, range);
 533
 534     // Test the methods that evaluated both functions for one value
 535     real     x        = 0.5 * (range.first + range.second);
 536     real     refFunc0 = lj6Function(x);
 537     real     refDer0  = lj6Derivative(x);
 538     real     refFunc1 = lj12Function(x);
 539     real     refDer1  = lj12Derivative(x);
 540
 541     real     tstFunc0, tstDer0, tstFunc1, tstDer1;
 542     real     tmpFunc0, tmpFunc1, tmpDer0, tmpDer1;
 543
 544     // test that we reproduce the reference functions
 545     table.evaluateFunctionAndDerivative(x, &tstFunc0, &tstDer0, &tstFunc1, &tstDer1);
 546
 547     real funcErr0 = std::abs(tstFunc0-refFunc0) / std::abs(refFunc0);
 548     real funcErr1 = std::abs(tstFunc1-refFunc1) / std::abs(refFunc1);
 549     real derErr0  = std::abs(tstDer0-refDer0) / std::abs(refDer0);
 550     real derErr1  = std::abs(tstDer1-refDer1) / std::abs(refDer1);
 551
 552     // Use asserts, since the following ones compare to these values.
 553     ASSERT_LT(funcErr0, TypeParam::defaultTolerance);
 554     ASSERT_LT(derErr0, TypeParam::defaultTolerance);
 555     ASSERT_LT(funcErr1, TypeParam::defaultTolerance);
 556     ASSERT_LT(derErr1, TypeParam::defaultTolerance);
 557
 558     // Test that function/derivative-only interpolation methods work
 559     table.evaluateFunction(x, &tmpFunc0, &tmpFunc1);
 560     table.evaluateDerivative(x, &tmpDer0, &tmpDer1);
 561     EXPECT_EQ(tstFunc0, tmpFunc0);
 562     EXPECT_EQ(tstFunc1, tmpFunc1);
 563     EXPECT_EQ(tstDer0, tmpDer0);
 564
 565     // Test that scrambled order interpolation methods work
 566     table.template evaluateFunctionAndDerivative<2, 1, 0>(x, &tstFunc1, &tstDer1, &tstFunc0, &tstDer0);
 567     EXPECT_EQ(tstFunc0, tmpFunc0);
 568     EXPECT_EQ(tstFunc1, tmpFunc1);
 569     EXPECT_EQ(tstDer0, tmpDer0);
 570     EXPECT_EQ(tstDer1, tmpDer1);
 571
 572     // Test scrambled order for function/derivative-only methods
 573     table.template evaluateFunction<2, 1, 0>(x, &tmpFunc1, &tmpFunc0);
 574     table.template evaluateDerivative<2, 1, 0>(x, &tmpDer1, &tmpDer0);
 575     EXPECT_EQ(tstFunc0, tmpFunc0);
 576     EXPECT_EQ(tstFunc1, tmpFunc1);
 577     EXPECT_EQ(tstDer0, tmpDer0);
 578     EXPECT_EQ(tstDer1, tmpDer1);
 579 }
 580
 581 TYPED_TEST(SplineTableTest, ThreeFunctions)
 582 {
 583     std::pair<real, real>  range(0.2, 2.0);
 584
 585     TypeParam              table( {{"Coulomb", coulombFunction, coulombDerivative}, {"LJ6", lj6Function, lj6Derivative}, {"LJ12", lj12Function, lj12Derivative}}, range);
 586
 587     // Test entire range for each function
 588     TestFixture::template testSplineTableAgainstFunctions<3, 0>("Coulomb", coulombFunction, coulombDerivative, table, range);
 589     TestFixture::template testSplineTableAgainstFunctions<3, 1>("LJ6", lj6Function, lj6Derivative, table, range);
 590     TestFixture::template testSplineTableAgainstFunctions<3, 2>("LJ12", lj12Function, lj12Derivative, table, range);
 591
 592     // Test the methods that evaluated both functions for one value
 593     real     x        = 0.5 * (range.first + range.second);
 594     real     refFunc0 = coulombFunction(x);
 595     real     refDer0  = coulombDerivative(x);
 596     real     refFunc1 = lj6Function(x);
 597     real     refDer1  = lj6Derivative(x);
 598     real     refFunc2 = lj12Function(x);
 599     real     refDer2  = lj12Derivative(x);
 600
 601     real     tstFunc0, tstDer0, tstFunc1, tstDer1, tstFunc2, tstDer2;
 602     real     tmpFunc0, tmpFunc1, tmpFunc2, tmpDer0, tmpDer1, tmpDer2;
 603
 604     // test that we reproduce the reference functions
 605     table.evaluateFunctionAndDerivative(x, &tstFunc0, &tstDer0, &tstFunc1, &tstDer1, &tstFunc2, &tstDer2);
 606
 607     real funcErr0 = std::abs(tstFunc0-refFunc0) / std::abs(refFunc0);
 608     real derErr0  = std::abs(tstDer0-refDer0) / std::abs(refDer0);
 609     real funcErr1 = std::abs(tstFunc1-refFunc1) / std::abs(refFunc1);
 610     real derErr1  = std::abs(tstDer1-refDer1) / std::abs(refDer1);
 611     real funcErr2 = std::abs(tstFunc2-refFunc2) / std::abs(refFunc2);
 612     real derErr2  = std::abs(tstDer2-refDer2) / std::abs(refDer2);
 613
 614     // Use asserts, since the following ones compare to these values.
 615     ASSERT_LT(funcErr0, TypeParam::defaultTolerance);
 616     ASSERT_LT(derErr0, TypeParam::defaultTolerance);
 617     ASSERT_LT(funcErr1, TypeParam::defaultTolerance);
 618     ASSERT_LT(derErr1, TypeParam::defaultTolerance);
 619     ASSERT_LT(funcErr2, TypeParam::defaultTolerance);
 620     ASSERT_LT(derErr2, TypeParam::defaultTolerance);
 621
 622     // Test that function/derivative-only interpolation methods work
 623     table.evaluateFunction(x, &tmpFunc0, &tmpFunc1, &tmpFunc2);
 624     table.evaluateDerivative(x, &tmpDer0, &tmpDer1, &tmpDer2);
 625     EXPECT_EQ(tstFunc0, tmpFunc0);
 626     EXPECT_EQ(tstFunc1, tmpFunc1);
 627     EXPECT_EQ(tstFunc2, tmpFunc2);
 628     EXPECT_EQ(tstDer0, tmpDer0);
 629     EXPECT_EQ(tstDer1, tmpDer1);
 630     EXPECT_EQ(tstDer2, tmpDer2);
 631
 632     // Test two functions out of three
 633     table.template evaluateFunctionAndDerivative<3, 0, 1>(x, &tmpFunc0, &tmpDer0, &tmpFunc1, &tmpDer1);
 634     EXPECT_EQ(tstFunc0, tmpFunc0);
 635     EXPECT_EQ(tstFunc1, tmpFunc1);
 636     EXPECT_EQ(tstDer0, tmpDer0);
 637     EXPECT_EQ(tstDer1, tmpDer1);
 638
 639     // two out of three, function/derivative-only
 640     table.template evaluateFunction<3, 0, 1>(x, &tmpFunc0, &tmpFunc1);
 641     table.template evaluateDerivative<3, 0, 1>(x, &tmpDer0, &tmpDer1);
 642     EXPECT_EQ(tstFunc0, tmpFunc0);
 643     EXPECT_EQ(tstFunc1, tmpFunc1);
 644     EXPECT_EQ(tstDer0, tmpDer0);
 645     EXPECT_EQ(tstDer1, tmpDer1);
 646
 647     // Test that scrambled order interpolation methods work
 648     table.template evaluateFunctionAndDerivative<3, 2, 1, 0>(x, &tstFunc2, &tstDer2, &tstFunc1, &tstDer1, &tstFunc0, &tstDer0);
 649     EXPECT_EQ(tstFunc0, tmpFunc0);
 650     EXPECT_EQ(tstFunc1, tmpFunc1);
 651     EXPECT_EQ(tstFunc2, tmpFunc2);
 652     EXPECT_EQ(tstDer0, tmpDer0);
 653     EXPECT_EQ(tstDer1, tmpDer1);
 654     EXPECT_EQ(tstDer2, tmpDer2);
 655
 656     // Test scrambled order for function/derivative-only methods
 657     table.template evaluateFunction<3, 2, 1, 0>(x, &tmpFunc2, &tmpFunc1, &tmpFunc0);
 658     table.template evaluateDerivative<3, 2, 1, 0>(x, &tmpDer2, &tmpDer1, &tmpDer0);
 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
 667 #if GMX_SIMD_HAVE_REAL
 668 TYPED_TEST(SplineTableTest, Simd)
 669 {
 670     std::pair<real, real>  range(0.2, 1.0);
 671     TypeParam              table( {{"LJ12", lj12Function, lj12Derivative}}, range);
 672
 673     // We already test that the SIMD operations handle the different elements
 674     // correctly in the SIMD module, so here we only test that interpolation
 675     // works for a single value in the middle of the interval
 676
 677     real     x       = 0.5 * (range.first + range.second);
 678     real     refFunc = lj12Function(x);
 679     real     refDer  = lj12Derivative(x);
 680     SimdReal tstFunc, tstDer;
 681     real     funcErr, derErr;
 682     alignas(GMX_SIMD_ALIGNMENT) real alignedMem[GMX_SIMD_REAL_WIDTH];
 683
 684     table.evaluateFunctionAndDerivative(SimdReal(x), &tstFunc, &tstDer);
 685
 686     store(alignedMem, tstFunc);
 687     funcErr = std::abs(alignedMem[0]-refFunc) / std::abs(refFunc);
 688
 689     store(alignedMem, tstDer);
 690     derErr  = std::abs(alignedMem[0]-refDer ) / std::abs(refDer );
 691
 692     EXPECT_LT(funcErr, TypeParam::defaultTolerance);
 693     EXPECT_LT(derErr, TypeParam::defaultTolerance);
 694 }
 695
 696 TYPED_TEST(SplineTableTest, SimdTwoFunctions)
 697 {
 698     std::pair<real, real>  range(0.2, 2.0);
 699
 700     TypeParam              table( {{"LJ6", lj6Function, lj6Derivative}, {"LJ12", lj12Function, lj12Derivative}}, range);
 701
 702     // We already test that the SIMD operations handle the different elements
 703     // correctly in the SIMD module, so here we only test that interpolation
 704     // works for a single value in the middle of the interval
 705
 706     real     x        = 0.5 * (range.first + range.second);
 707     real     refFunc0 = lj6Function(x);
 708     real     refDer0  = lj6Derivative(x);
 709     real     refFunc1 = lj12Function(x);
 710     real     refDer1  = lj12Derivative(x);
 711     SimdReal tstFunc0, tstDer0;
 712     SimdReal tstFunc1, tstDer1;
 713     real     funcErr0, derErr0;
 714     real     funcErr1, derErr1;
 715     alignas(GMX_SIMD_ALIGNMENT) real alignedMem[GMX_SIMD_REAL_WIDTH];
 716
 717     table.evaluateFunctionAndDerivative(SimdReal(x), &tstFunc0, &tstDer0, &tstFunc1, &tstDer1);
 718
 719     store(alignedMem, tstFunc0);
 720     funcErr0 = std::abs(alignedMem[0]-refFunc0) / std::abs(refFunc0);
 721
 722     store(alignedMem, tstDer0);
 723     derErr0  = std::abs(alignedMem[0]-refDer0 ) / std::abs(refDer0 );
 724
 725     store(alignedMem, tstFunc1);
 726     funcErr1 = std::abs(alignedMem[0]-refFunc1) / std::abs(refFunc1);
 727
 728     store(alignedMem, tstDer1);
 729     derErr1  = std::abs(alignedMem[0]-refDer1 ) / std::abs(refDer1 );
 730
 731     EXPECT_LT(funcErr0, TypeParam::defaultTolerance);
 732     EXPECT_LT(derErr0, TypeParam::defaultTolerance);
 733     EXPECT_LT(funcErr1, TypeParam::defaultTolerance);
 734     EXPECT_LT(derErr1, TypeParam::defaultTolerance);
 735 }
 736 #endif
 737
 738 #if GMX_SIMD_HAVE_REAL && !defined NDEBUG
 739 TYPED_TEST(SplineTableTest, CatchesOutOfRangeValuesSimd)
 740 {
 741     std::pair<real, real>                   range(0.2, 1.0);
 742     TypeParam                               table( {{"LJ12", lj12Function, lj12Derivative}}, range);
 743     SimdReal                                x, func, der;
 744
 745     AlignedArray<real, GMX_SIMD_REAL_WIDTH> alignedMem;
 746
 747     alignedMem.fill(range.first);
 748     // Make position 1 incorrect if width>=2, otherwise position 0
 749     // range.first-GMX_REAL_EPS is not invalid. See comment in table.
 750     alignedMem[ (GMX_SIMD_REAL_WIDTH >= 2) ? 1 : 0] = -GMX_REAL_EPS;
 751     x = load<SimdReal>(alignedMem);
 752
 753     EXPECT_THROW_GMX(table.evaluateFunctionAndDerivative(x, &func, &der), gmx::RangeError);
 754
 755     // Make position 1 incorrect if width>=2, otherwise position 0
 756     alignedMem[ (GMX_SIMD_REAL_WIDTH >= 2) ? 1 : 0] = range.second;
 757     x = load<SimdReal>(alignedMem);
 758
 759     EXPECT_THROW_GMX(table.evaluateFunctionAndDerivative(x, &func, &der), gmx::RangeError);
 760 }
 761
 762 TYPED_TEST(SplineTableTest, AcceptsInRangeValuesSimd)
 763 {
 764     std::pair<real, real>  range(0.2, 1.0);
 765     TypeParam              table( {{"LJ12", lj12Function, lj12Derivative}}, range);
 766     SimdReal               x, func, der;
 767
 768     alignas(GMX_SIMD_ALIGNMENT) real alignedMem[GMX_SIMD_REAL_WIDTH];
 769
 770     // Test all values between 0 and range.second
 771     for (std::size_t i = 0; i < GMX_SIMD_REAL_WIDTH; i++)
 772     {
 773         alignedMem[i] = range.second*(1.0-GMX_REAL_EPS)*i/(GMX_SIMD_REAL_WIDTH-1);
 774     }
 775     x = load<SimdReal>(alignedMem);
 776
 777     EXPECT_NO_THROW_GMX(table.evaluateFunctionAndDerivative(x, &func, &der));
 778 }
 779
 780 #endif
 781
 782 } // namespace
 783
 784 } // namespace test
 785
 786 } // namespace gmx