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,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
  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()), range_(range)
 310 {
 311     // Sanity check on input values
 312     if (range_.first < 0.0 || (range_.second - range_.first) < 0.001)
 313     {
 314         GMX_THROW(InvalidInputError(
 315                 "Range to tabulate cannot include negative values and must span at least 0.001"));
 316     }
 317
 318     if (tolerance < GMX_REAL_EPS)
 319     {
 320         GMX_THROW(ToleranceError("Table tolerance cannot be smaller than GMX_REAL_EPS"));
 321     }
 322
 323     double minQuotient = GMX_REAL_MAX;
 324
 325     // loop over all functions to find smallest spacing
 326     for (const auto& thisFuncInput : analyticalInputList)
 327     {
 328         try
 329         {
 330             internal::throwUnlessDerivativeIsConsistentWithFunction(
 331                     thisFuncInput.function, thisFuncInput.derivative, range_);
 332         }
 333         catch (gmx::GromacsException& ex)
 334         {
 335             ex.prependContext("Error generating cubic spline table for function '"
 336                               + thisFuncInput.desc + "'");
 337             throw;
 338         }
 339         // Calculate the required table spacing h. The error we make with a third order polynomial
 340         // (second order for derivative) will be described by the fourth-derivative correction term.
 341         //
 342         // This means we can compute the required spacing as h =
 343         // 0.5*cbrt(72*sqrt(3)*tolerance**min(f'/f'''')), where f'/f'''' is the first and fourth
 344         // derivative of the function, respectively. Since we already have an analytical form of the
 345         // derivative, we reduce the numerical errors by calculating the quotient of the function
 346         // and third derivative of the input-derivative-analytical function instead.
 347
 348         double thisMinQuotient = internal::findSmallestQuotientOfFunctionAndThirdDerivative(
 349                 thisFuncInput.derivative, range_);
 350
 351         minQuotient = std::min(minQuotient, thisMinQuotient);
 352     }
 353
 354     double spacing = 0.5 * std::cbrt(72.0 * std::sqrt(3.0) * tolerance * minQuotient);
 355
 356     tableScale_ = 1.0 / spacing;
 357
 358     if (range_.second * tableScale_ > 2e6)
 359     {
 360         GMX_THROW(
 361                 ToleranceError("Over a million points would be required for table; decrease range "
 362                                "or increase tolerance"));
 363     }
 364
 365     // Loop over all tables again.
 366     // Here we create the actual table for each function, and then
 367     // combine them into a multiplexed table function.
 368     std::size_t funcIndex = 0;
 369
 370     for (const auto& thisFuncInput : analyticalInputList)
 371     {
 372         try
 373         {
 374             std::vector<real> tmpYfghTableData;
 375
 376             fillSingleCubicSplineTableData(
 377                     thisFuncInput.function, thisFuncInput.derivative, range_, spacing, &tmpYfghTableData);
 378
 379             internal::fillMultiplexedTableData(
 380                     tmpYfghTableData, &yfghMultiTableData_, 4, numFuncInTable_, funcIndex);
 381
 382             funcIndex++;
 383         }
 384         catch (gmx::GromacsException& ex)
 385         {
 386             ex.prependContext("Error generating cubic spline table for function '"
 387                               + thisFuncInput.desc + "'");
 388             throw;
 389         }
 390     }
 391 }
 392
 393
 394 CubicSplineTable::CubicSplineTable(std::initializer_list<NumericalSplineTableInput> numericalInputList,
 395                                    const std::pair<real, real>&                     range,
 396                                    real                                             tolerance) :
 397     numFuncInTable_(numericalInputList.size()), range_(range)
 398 {
 399     // Sanity check on input values
 400     if (range.first < 0.0 || (range.second - range.first) < 0.001)
 401     {
 402         GMX_THROW(InvalidInputError(
 403                 "Range to tabulate cannot include negative values and must span at least 0.001"));
 404     }
 405
 406     if (tolerance < GMX_REAL_EPS)
 407     {
 408         GMX_THROW(ToleranceError("Table tolerance cannot be smaller than GMX_REAL_EPS"));
 409     }
 410
 411     double minQuotient = GMX_REAL_MAX;
 412
 413     // loop over all functions to find smallest spacing
 414     for (auto thisFuncInput : numericalInputList)
 415     {
 416         try
 417         {
 418             // We do not yet know what the margin is, but we need to test that we at least cover
 419             // the requested range before starting to calculate derivatives
 420             if (thisFuncInput.function.size() < range_.second / thisFuncInput.spacing + 1)
 421             {
 422                 GMX_THROW(
 423                         InconsistentInputError("Table input vectors must cover requested range, "
 424                                                "and a margin beyond the upper endpoint"));
 425             }
 426
 427             if (thisFuncInput.function.size() != thisFuncInput.derivative.size())
 428             {
 429                 GMX_THROW(InconsistentInputError(
 430                         "Function and derivative vectors have different lengths"));
 431             }
 432
 433             internal::throwUnlessDerivativeIsConsistentWithFunction(
 434                     thisFuncInput.function, thisFuncInput.derivative, thisFuncInput.spacing, range_);
 435         }
 436         catch (gmx::GromacsException& ex)
 437         {
 438             ex.prependContext("Error generating cubic spline table for function '"
 439                               + thisFuncInput.desc + "'");
 440             throw;
 441         }
 442         // Calculate the required table spacing h. The error we make with linear interpolation
 443         // of the derivative will be described by the third-derivative correction term.
 444         // This means we can compute the required spacing as h = sqrt(12*tolerance*min(f'/f''')),
 445         // where f'/f''' is the first and third derivative of the function, respectively.
 446
 447         double thisMinQuotient = internal::findSmallestQuotientOfFunctionAndThirdDerivative(
 448                 thisFuncInput.derivative, thisFuncInput.spacing, range_);
 449
 450         minQuotient = std::min(minQuotient, thisMinQuotient);
 451     }
 452
 453     double spacing = std::cbrt(72.0 * std::sqrt(3.0) * tolerance * minQuotient);
 454
 455     tableScale_ = 1.0 / spacing;
 456
 457     if (range_.second * tableScale_ > 1e6)
 458     {
 459         GMX_THROW(
 460                 ToleranceError("Requested tolerance would require over a million points in table"));
 461     }
 462
 463     // Loop over all tables again.
 464     // Here we create the actual table for each function, and then
 465     // combine them into a multiplexed table function.
 466     std::size_t funcIndex = 0;
 467
 468     for (auto thisFuncInput : numericalInputList)
 469     {
 470         try
 471         {
 472             if (spacing < thisFuncInput.spacing)
 473             {
 474                 GMX_THROW(
 475                         ToleranceError("Input vector spacing cannot achieve tolerance requested"));
 476             }
 477
 478             std::vector<real> tmpYfghTableData;
 479
 480             fillSingleCubicSplineTableData(thisFuncInput.function,
 481                                            thisFuncInput.derivative,
 482                                            thisFuncInput.spacing,
 483                                            range,
 484                                            spacing,
 485                                            &tmpYfghTableData);
 486
 487             internal::fillMultiplexedTableData(
 488                     tmpYfghTableData, &yfghMultiTableData_, 4, numFuncInTable_, funcIndex);
 489
 490             funcIndex++;
 491         }
 492         catch (gmx::GromacsException& ex)
 493         {
 494             ex.prependContext("Error generating cubic spline table for function '"
 495                               + thisFuncInput.desc + "'");
 496             throw;
 497         }
 498     }
 499 }
 500
 501 } // namespace gmx