src/gromacs/tables/cubicsplinetable.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2016,2017,2018,2019,2020, 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
  36 /*! \internal \file
  37  * \brief
  38  * Implements classes for cubic spline table functions
  39  *
  40  * \author Erik Lindahl <erik.lindahl@gmail.com>
  41  * \ingroup module_tables
  42  */
  43 #include "gmxpre.h"
  44
  45 #include "cubicsplinetable.h"
  46
  47 #include <cmath>
  48
  49 #include <algorithm>
  50 #include <functional>
  51 #include <initializer_list>
  52 #include <utility>
  53 #include <vector>
  54
  55 #include "gromacs/tables/tableinput.h"
  56 #include "gromacs/utility/alignedallocator.h"
  57 #include "gromacs/utility/arrayref.h"
  58 #include "gromacs/utility/exceptions.h"
  59 #include "gromacs/utility/real.h"
  60
  61 #include "splineutil.h"
  62
  63 namespace gmx
  64 {
  65
  66 namespace
  67 {
  68
  69 /*! \brief Calculate table elements from function/derivative data
  70  *
  71  * \param functionValue0   Function value for the present table index
  72  * \param functionValue1   Function value for the next table index
  73  * \param derivativeValue0 Derivative value for the present table index
  74  * \param derivativeValue1 Derivative value for the next table index
  75  * \param spacing          Distance between table points
  76  * \param Y                Function value for table index
  77  * \param F                Component to multiply with offset eps
  78  * \param G                Component to multiply with eps^2
  79  * \param H                Component to multiply with eps^3
  80  */
  81 void calculateCubicSplineCoefficients(double  functionValue0,
  82                                       double  functionValue1,
  83                                       double  derivativeValue0,
  84                                       double  derivativeValue1,
  85                                       double  spacing,
  86                                       double* Y,
  87                                       double* F,
  88                                       double* G,
  89                                       double* H)
  90 {
  91     *Y = functionValue0;
  92     *F = spacing * derivativeValue0;
  93     *G = 3.0 * (functionValue1 - functionValue0) - spacing * (derivativeValue1 + 2.0 * derivativeValue0);
  94     *H = -2.0 * (functionValue1 - functionValue0) + spacing * (derivativeValue1 + derivativeValue0);
  95 }
  96
  97 /*! \brief Perform cubic spline interpolation in interval from function/derivative
  98  *
  99  * \param      functionValue0               Function value for the present table index
 100  * \param      functionValue1               Function value for the next table index
 101  * \param      derivativeValue0             Derivative value for the present table index
 102  * \param      derivativeValue1             Derivative value for the next table index
 103  * \param      spacing                      Distance between table points
 104  * \param      eps                          Offset from lower table point for evaluation
 105  * \param[out] interpolatedFunctionValue    Output function value
 106  * \param[out] interpolatedDerivativeValue  Output derivative value
 107  */
 108 void cubicSplineInterpolationFromFunctionAndDerivative(double  functionValue0,
 109                                                        double  functionValue1,
 110                                                        double  derivativeValue0,
 111                                                        double  derivativeValue1,
 112                                                        double  spacing,
 113                                                        double  eps,
 114                                                        double* interpolatedFunctionValue,
 115                                                        double* interpolatedDerivativeValue)
 116 {
 117     double Y, F, G, H;
 118
 119     calculateCubicSplineCoefficients(
 120             functionValue0, functionValue1, derivativeValue0, derivativeValue1, spacing, &Y, &F, &G, &H);
 121
 122     double Fp = fma(fma(H, eps, G), eps, F);
 123
 124     *interpolatedFunctionValue   = fma(Fp, eps, Y);
 125     *interpolatedDerivativeValue = fma(eps, fma(2.0 * eps, H, G), Fp) / spacing;
 126 }
 127
 128
 129 /*! \brief Construct the data for a single cubic table from analytical functions
 130  *
 131  * \param[in]  function             Analytical functiojn
 132  * \param[in]  derivative           Analytical derivative
 133  * \param[in]  range                Upper/lower limit of region to tabulate
 134  * \param[in]  spacing              Distance between table points
 135  * \param[out] yfghTableData        Output cubic spline table with Y,F,G,H entries
 136  */
 137 void fillSingleCubicSplineTableData(const std::function<double(double)>& function,
 138                                     const std::function<double(double)>& derivative,
 139                                     const std::pair<real, real>&         range,
 140                                     double                               spacing,
 141                                     std::vector<real>*                   yfghTableData)
 142 {
 143     int endIndex = static_cast<int>(range.second / spacing + 2);
 144
 145     yfghTableData->resize(4 * endIndex);
 146
 147     double      maxMagnitude      = 0.0001 * GMX_REAL_MAX;
 148     bool        functionIsInRange = true;
 149     std::size_t lastIndexInRange  = endIndex - 1;
 150
 151     for (int i = endIndex - 1; i >= 0; i--)
 152     {
 153         double x = i * spacing;
 154         double tmpFunctionValue;
 155         double tmpDerivativeValue;
 156         double nextHigherFunction;
 157         double nextHigherDerivative;
 158         double Y, F, G, H;
 159
 160         if (range.first > 0 && i == 0)
 161         {
 162             // Avoid x==0 if it is not in the range, since it can lead to
 163             // singularities even if the value for i==1 was within or required magnitude
 164             functionIsInRange = false;
 165         }
 166
 167         if (functionIsInRange)
 168         {
 169             tmpFunctionValue     = function(x);
 170             tmpDerivativeValue   = derivative(x);
 171             nextHigherFunction   = ((i + 1) < endIndex) ? function(x + spacing) : 0.0;
 172             nextHigherDerivative = ((i + 1) < endIndex) ? derivative(x + spacing) : 0.0;
 173
 174             if (std::abs(tmpFunctionValue) > maxMagnitude || std::abs(tmpDerivativeValue) > maxMagnitude)
 175             {
 176                 functionIsInRange = false; // Once this happens, it never resets to true again
 177             }
 178         }
 179
 180         if (functionIsInRange)
 181         {
 182             calculateCubicSplineCoefficients(tmpFunctionValue,
 183                                              nextHigherFunction,
 184                                              tmpDerivativeValue,
 185                                              nextHigherDerivative,
 186                                              spacing,
 187                                              &Y,
 188                                              &F,
 189                                              &G,
 190                                              &H);
 191             lastIndexInRange--;
 192         }
 193         else
 194         {
 195             double lastIndexY = (*yfghTableData)[4 * lastIndexInRange];
 196             double lastIndexF = (*yfghTableData)[4 * lastIndexInRange + 1];
 197
 198             Y = lastIndexY + lastIndexF * (i - lastIndexInRange);
 199             F = lastIndexF;
 200             G = 0.0;
 201             H = 0.0;
 202         }
 203
 204         (*yfghTableData)[4 * i]     = Y;
 205         (*yfghTableData)[4 * i + 1] = F;
 206         (*yfghTableData)[4 * i + 2] = G;
 207         (*yfghTableData)[4 * i + 3] = H;
 208     }
 209 }
 210
 211
 212 /*! \brief Construct the data for a single cubic table from vector data
 213  *
 214  * \param[in]  function             Input vector with function data
 215  * \param[in]  derivative           Input vector with derivative data
 216  * \param[in]  inputSpacing         Distance between points in input vectors
 217  * \param[in]  range                Upper/lower limit of region to tabulate
 218  * \param[in]  spacing              Distance between table points
 219  * \param[out] yfghTableData        Output cubic spline table with Y,F,G,H entries
 220  */
 221 void fillSingleCubicSplineTableData(ArrayRef<const double>       function,
 222                                     ArrayRef<const double>       derivative,
 223                                     double                       inputSpacing,
 224                                     const std::pair<real, real>& range,
 225                                     double                       spacing,
 226                                     std::vector<real>*           yfghTableData)
 227 {
 228     int endIndex = static_cast<int>(range.second / spacing + 2);
 229
 230     std::vector<double> tmpFunction(endIndex);
 231     std::vector<double> tmpDerivative(endIndex);
 232
 233     double      maxMagnitude      = 0.0001 * GMX_REAL_MAX;
 234     bool        functionIsInRange = true;
 235     std::size_t lastIndexInRange  = endIndex - 1;
 236
 237     // Interpolate function and derivative values in positions needed for output
 238     for (int i = endIndex - 1; i >= 0; i--)
 239     {
 240         double x     = i * spacing;
 241         double xtab  = x / inputSpacing;
 242         int    index = static_cast<int>(xtab);
 243         double eps   = xtab - index;
 244
 245         if (range.first > 0 && i == 0)
 246         {
 247             // Avoid x==0 if it is not in the range, since it can lead to
 248             // singularities even if the value for i==1 was within or required magnitude
 249             functionIsInRange = false;
 250         }
 251
 252         if (functionIsInRange
 253             && (std::abs(function[index]) > maxMagnitude || std::abs(derivative[index]) > maxMagnitude))
 254         {
 255             functionIsInRange = false; // Once this happens, it never resets to true again
 256         }
 257
 258         if (functionIsInRange)
 259         {
 260             cubicSplineInterpolationFromFunctionAndDerivative(function[index],
 261                                                               function[index + 1],
 262                                                               derivative[index],
 263                                                               derivative[index + 1],
 264                                                               inputSpacing,
 265                                                               eps,
 266                                                               &(tmpFunction[i]),
 267                                                               &(tmpDerivative[i]));
 268             lastIndexInRange--;
 269         }
 270         else
 271         {
 272             double lastIndexFunction   = tmpFunction[lastIndexInRange];
 273             double lastIndexDerivative = tmpDerivative[lastIndexInRange];
 274             tmpFunction[i] = lastIndexFunction + lastIndexDerivative * (i - lastIndexInRange) * spacing;
 275             tmpDerivative[i] = lastIndexDerivative;
 276         }
 277     }
 278
 279     yfghTableData->resize(4 * endIndex);
 280
 281     for (int i = 0; i < endIndex; i++)
 282     {
 283         double Y, F, G, H;
 284
 285         double nextFunction   = ((i + 1) < endIndex) ? tmpFunction[i + 1] : 0.0;
 286         double nextDerivative = ((i + 1) < endIndex) ? tmpDerivative[i + 1] : 0.0;
 287
 288         calculateCubicSplineCoefficients(
 289                 tmpFunction[i], nextFunction, tmpDerivative[i], nextDerivative, spacing, &Y, &F, &G, &H);
 290         (*yfghTableData)[4 * i]     = Y;
 291         (*yfghTableData)[4 * i + 1] = F;
 292         (*yfghTableData)[4 * i + 2] = G;
 293         (*yfghTableData)[4 * i + 3] = H;
 294     }
 295 }
 296
 297 } // namespace
 298
 299
 300 #if GMX_DOUBLE
 301 const real CubicSplineTable::defaultTolerance = 1e-10;
 302 #else
 303 const real CubicSplineTable::defaultTolerance = 10.0 * GMX_FLOAT_EPS;
 304 #endif
 305
 306 CubicSplineTable::CubicSplineTable(std::initializer_list<AnalyticalSplineTableInput> analyticalInputList,
 307                                    const std::pair<real, real>&                      range,
 308                                    real                                              tolerance) :
 309     numFuncInTable_(analyticalInputList.size()),
 310     range_(range)
 311 {
 312     // Sanity check on input values
 313     if (range_.first < 0.0 || (range_.second - range_.first) < 0.001)
 314     {
 315         GMX_THROW(InvalidInputError(
 316                 "Range to tabulate cannot include negative values and must span at least 0.001"));
 317     }
 318
 319     if (tolerance < GMX_REAL_EPS)
 320     {
 321         GMX_THROW(ToleranceError("Table tolerance cannot be smaller than GMX_REAL_EPS"));
 322     }
 323
 324     double minQuotient = GMX_REAL_MAX;
 325
 326     // loop over all functions to find smallest spacing
 327     for (const auto& thisFuncInput : analyticalInputList)
 328     {
 329         try
 330         {
 331             internal::throwUnlessDerivativeIsConsistentWithFunction(
 332                     thisFuncInput.function, thisFuncInput.derivative, range_);
 333         }
 334         catch (gmx::GromacsException& ex)
 335         {
 336             ex.prependContext("Error generating cubic spline table for function '"
 337                               + thisFuncInput.desc + "'");
 338             throw;
 339         }
 340         // Calculate the required table spacing h. The error we make with a third order polynomial
 341         // (second order for derivative) will be described by the fourth-derivative correction term.
 342         //
 343         // This means we can compute the required spacing as h =
 344         // 0.5*cbrt(72*sqrt(3)*tolerance**min(f'/f'''')), where f'/f'''' is the first and fourth
 345         // derivative of the function, respectively. Since we already have an analytical form of the
 346         // derivative, we reduce the numerical errors by calculating the quotient of the function
 347         // and third derivative of the input-derivative-analytical function instead.
 348
 349         double thisMinQuotient = internal::findSmallestQuotientOfFunctionAndThirdDerivative(
 350                 thisFuncInput.derivative, range_);
 351
 352         minQuotient = std::min(minQuotient, thisMinQuotient);
 353     }
 354
 355     double spacing = 0.5 * std::cbrt(72.0 * std::sqrt(3.0) * tolerance * minQuotient);
 356
 357     tableScale_ = 1.0 / spacing;
 358
 359     if (range_.second * tableScale_ > 2e6)
 360     {
 361         GMX_THROW(
 362                 ToleranceError("Over a million points would be required for table; decrease range "
 363                                "or increase tolerance"));
 364     }
 365
 366     // Loop over all tables again.
 367     // Here we create the actual table for each function, and then
 368     // combine them into a multiplexed table function.
 369     std::size_t funcIndex = 0;
 370
 371     for (const auto& thisFuncInput : analyticalInputList)
 372     {
 373         try
 374         {
 375             std::vector<real> tmpYfghTableData;
 376
 377             fillSingleCubicSplineTableData(
 378                     thisFuncInput.function, thisFuncInput.derivative, range_, spacing, &tmpYfghTableData);
 379
 380             internal::fillMultiplexedTableData(
 381                     tmpYfghTableData, &yfghMultiTableData_, 4, numFuncInTable_, funcIndex);
 382
 383             funcIndex++;
 384         }
 385         catch (gmx::GromacsException& ex)
 386         {
 387             ex.prependContext("Error generating cubic spline table for function '"
 388                               + thisFuncInput.desc + "'");
 389             throw;
 390         }
 391     }
 392 }
 393
 394
 395 CubicSplineTable::CubicSplineTable(std::initializer_list<NumericalSplineTableInput> numericalInputList,
 396                                    const std::pair<real, real>&                     range,
 397                                    real                                             tolerance) :
 398     numFuncInTable_(numericalInputList.size()),
 399     range_(range)
 400 {
 401     // Sanity check on input values
 402     if (range.first < 0.0 || (range.second - range.first) < 0.001)
 403     {
 404         GMX_THROW(InvalidInputError(
 405                 "Range to tabulate cannot include negative values and must span at least 0.001"));
 406     }
 407
 408     if (tolerance < GMX_REAL_EPS)
 409     {
 410         GMX_THROW(ToleranceError("Table tolerance cannot be smaller than GMX_REAL_EPS"));
 411     }
 412
 413     double minQuotient = GMX_REAL_MAX;
 414
 415     // loop over all functions to find smallest spacing
 416     for (auto thisFuncInput : numericalInputList)
 417     {
 418         try
 419         {
 420             // We do not yet know what the margin is, but we need to test that we at least cover
 421             // the requested range before starting to calculate derivatives
 422             if (thisFuncInput.function.size() < range_.second / thisFuncInput.spacing + 1)
 423             {
 424                 GMX_THROW(
 425                         InconsistentInputError("Table input vectors must cover requested range, "
 426                                                "and a margin beyond the upper endpoint"));
 427             }
 428
 429             if (thisFuncInput.function.size() != thisFuncInput.derivative.size())
 430             {
 431                 GMX_THROW(InconsistentInputError(
 432                         "Function and derivative vectors have different lengths"));
 433             }
 434
 435             internal::throwUnlessDerivativeIsConsistentWithFunction(
 436                     thisFuncInput.function, thisFuncInput.derivative, thisFuncInput.spacing, range_);
 437         }
 438         catch (gmx::GromacsException& ex)
 439         {
 440             ex.prependContext("Error generating cubic spline table for function '"
 441                               + thisFuncInput.desc + "'");
 442             throw;
 443         }
 444         // Calculate the required table spacing h. The error we make with linear interpolation
 445         // of the derivative will be described by the third-derivative correction term.
 446         // This means we can compute the required spacing as h = sqrt(12*tolerance*min(f'/f''')),
 447         // where f'/f''' is the first and third derivative of the function, respectively.
 448
 449         double thisMinQuotient = internal::findSmallestQuotientOfFunctionAndThirdDerivative(
 450                 thisFuncInput.derivative, thisFuncInput.spacing, range_);
 451
 452         minQuotient = std::min(minQuotient, thisMinQuotient);
 453     }
 454
 455     double spacing = std::cbrt(72.0 * std::sqrt(3.0) * tolerance * minQuotient);
 456
 457     tableScale_ = 1.0 / spacing;
 458
 459     if (range_.second * tableScale_ > 1e6)
 460     {
 461         GMX_THROW(
 462                 ToleranceError("Requested tolerance would require over a million points in table"));
 463     }
 464
 465     // Loop over all tables again.
 466     // Here we create the actual table for each function, and then
 467     // combine them into a multiplexed table function.
 468     std::size_t funcIndex = 0;
 469
 470     for (auto thisFuncInput : numericalInputList)
 471     {
 472         try
 473         {
 474             if (spacing < thisFuncInput.spacing)
 475             {
 476                 GMX_THROW(
 477                         ToleranceError("Input vector spacing cannot achieve tolerance requested"));
 478             }
 479
 480             std::vector<real> tmpYfghTableData;
 481
 482             fillSingleCubicSplineTableData(thisFuncInput.function,
 483                                            thisFuncInput.derivative,
 484                                            thisFuncInput.spacing,
 485                                            range,
 486                                            spacing,
 487                                            &tmpYfghTableData);
 488
 489             internal::fillMultiplexedTableData(
 490                     tmpYfghTableData, &yfghMultiTableData_, 4, numFuncInTable_, funcIndex);
 491
 492             funcIndex++;
 493         }
 494         catch (gmx::GromacsException& ex)
 495         {
 496             ex.prependContext("Error generating cubic spline table for function '"
 497                               + thisFuncInput.desc + "'");
 498             throw;
 499         }
 500     }
 501 }
 502
 503 } // namespace gmx