src/gromacs/tables/cubicsplinetable.cpp

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