2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
38 * Implements classes for cubic spline table functions
40 * \author Erik Lindahl <erik.lindahl@gmail.com>
41 * \ingroup module_tables
45 #include "cubicsplinetable.h"
51 #include <initializer_list>
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"
61 #include "splineutil.h"
69 /*! \brief Calculate table elements from function/derivative data
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
81 void calculateCubicSplineCoefficients(double functionValue0,
82 double functionValue1,
83 double derivativeValue0,
84 double derivativeValue1,
92 *F = spacing * derivativeValue0;
93 *G = 3.0 * (functionValue1 - functionValue0) - spacing * (derivativeValue1 + 2.0 * derivativeValue0);
94 *H = -2.0 * (functionValue1 - functionValue0) + spacing * (derivativeValue1 + derivativeValue0);
97 /*! \brief Perform cubic spline interpolation in interval from function/derivative
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
108 void cubicSplineInterpolationFromFunctionAndDerivative(double functionValue0,
109 double functionValue1,
110 double derivativeValue0,
111 double derivativeValue1,
114 double* interpolatedFunctionValue,
115 double* interpolatedDerivativeValue)
119 calculateCubicSplineCoefficients(
120 functionValue0, functionValue1, derivativeValue0, derivativeValue1, spacing, &Y, &F, &G, &H);
122 double Fp = fma(fma(H, eps, G), eps, F);
124 *interpolatedFunctionValue = fma(Fp, eps, Y);
125 *interpolatedDerivativeValue = fma(eps, fma(2.0 * eps, H, G), Fp) / spacing;
129 /*! \brief Construct the data for a single cubic table from analytical functions
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
137 void fillSingleCubicSplineTableData(const std::function<double(double)>& function,
138 const std::function<double(double)>& derivative,
139 const std::pair<real, real>& range,
141 std::vector<real>* yfghTableData)
143 int endIndex = static_cast<int>(range.second / spacing + 2);
145 yfghTableData->resize(4 * endIndex);
147 double maxMagnitude = 0.0001 * GMX_REAL_MAX;
148 bool functionIsInRange = true;
149 std::size_t lastIndexInRange = endIndex - 1;
151 for (int i = endIndex - 1; i >= 0; i--)
153 double x = i * spacing;
154 double tmpFunctionValue;
155 double tmpDerivativeValue;
156 double nextHigherFunction;
157 double nextHigherDerivative;
160 if (range.first > 0 && i == 0)
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;
167 if (functionIsInRange)
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;
174 if (std::abs(tmpFunctionValue) > maxMagnitude || std::abs(tmpDerivativeValue) > maxMagnitude)
176 functionIsInRange = false; // Once this happens, it never resets to true again
180 if (functionIsInRange)
182 calculateCubicSplineCoefficients(tmpFunctionValue,
185 nextHigherDerivative,
195 double lastIndexY = (*yfghTableData)[4 * lastIndexInRange];
196 double lastIndexF = (*yfghTableData)[4 * lastIndexInRange + 1];
198 Y = lastIndexY + lastIndexF * (i - lastIndexInRange);
204 (*yfghTableData)[4 * i] = Y;
205 (*yfghTableData)[4 * i + 1] = F;
206 (*yfghTableData)[4 * i + 2] = G;
207 (*yfghTableData)[4 * i + 3] = H;
212 /*! \brief Construct the data for a single cubic table from vector data
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
221 void fillSingleCubicSplineTableData(ArrayRef<const double> function,
222 ArrayRef<const double> derivative,
224 const std::pair<real, real>& range,
226 std::vector<real>* yfghTableData)
228 int endIndex = static_cast<int>(range.second / spacing + 2);
230 std::vector<double> tmpFunction(endIndex);
231 std::vector<double> tmpDerivative(endIndex);
233 double maxMagnitude = 0.0001 * GMX_REAL_MAX;
234 bool functionIsInRange = true;
235 std::size_t lastIndexInRange = endIndex - 1;
237 // Interpolate function and derivative values in positions needed for output
238 for (int i = endIndex - 1; i >= 0; i--)
240 double x = i * spacing;
241 double xtab = x / inputSpacing;
242 int index = static_cast<int>(xtab);
243 double eps = xtab - index;
245 if (range.first > 0 && i == 0)
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;
252 if (functionIsInRange
253 && (std::abs(function[index]) > maxMagnitude || std::abs(derivative[index]) > maxMagnitude))
255 functionIsInRange = false; // Once this happens, it never resets to true again
258 if (functionIsInRange)
260 cubicSplineInterpolationFromFunctionAndDerivative(function[index],
263 derivative[index + 1],
267 &(tmpDerivative[i]));
272 double lastIndexFunction = tmpFunction[lastIndexInRange];
273 double lastIndexDerivative = tmpDerivative[lastIndexInRange];
274 tmpFunction[i] = lastIndexFunction + lastIndexDerivative * (i - lastIndexInRange) * spacing;
275 tmpDerivative[i] = lastIndexDerivative;
279 yfghTableData->resize(4 * endIndex);
281 for (int i = 0; i < endIndex; i++)
285 double nextFunction = ((i + 1) < endIndex) ? tmpFunction[i + 1] : 0.0;
286 double nextDerivative = ((i + 1) < endIndex) ? tmpDerivative[i + 1] : 0.0;
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;
301 const real CubicSplineTable::defaultTolerance = 1e-10;
303 const real CubicSplineTable::defaultTolerance = 10.0 * GMX_FLOAT_EPS;
306 CubicSplineTable::CubicSplineTable(std::initializer_list<AnalyticalSplineTableInput> analyticalInputList,
307 const std::pair<real, real>& range,
309 numFuncInTable_(analyticalInputList.size()),
312 // Sanity check on input values
313 if (range_.first < 0.0 || (range_.second - range_.first) < 0.001)
315 GMX_THROW(InvalidInputError(
316 "Range to tabulate cannot include negative values and must span at least 0.001"));
319 if (tolerance < GMX_REAL_EPS)
321 GMX_THROW(ToleranceError("Table tolerance cannot be smaller than GMX_REAL_EPS"));
324 double minQuotient = GMX_REAL_MAX;
326 // loop over all functions to find smallest spacing
327 for (const auto& thisFuncInput : analyticalInputList)
331 internal::throwUnlessDerivativeIsConsistentWithFunction(
332 thisFuncInput.function, thisFuncInput.derivative, range_);
334 catch (gmx::GromacsException& ex)
336 ex.prependContext("Error generating cubic spline table for function '"
337 + thisFuncInput.desc + "'");
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.
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.
349 double thisMinQuotient = internal::findSmallestQuotientOfFunctionAndThirdDerivative(
350 thisFuncInput.derivative, range_);
352 minQuotient = std::min(minQuotient, thisMinQuotient);
355 double spacing = 0.5 * std::cbrt(72.0 * std::sqrt(3.0) * tolerance * minQuotient);
357 tableScale_ = 1.0 / spacing;
359 if (range_.second * tableScale_ > 2e6)
362 ToleranceError("Over a million points would be required for table; decrease range "
363 "or increase tolerance"));
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;
371 for (const auto& thisFuncInput : analyticalInputList)
375 std::vector<real> tmpYfghTableData;
377 fillSingleCubicSplineTableData(
378 thisFuncInput.function, thisFuncInput.derivative, range_, spacing, &tmpYfghTableData);
380 internal::fillMultiplexedTableData(
381 tmpYfghTableData, &yfghMultiTableData_, 4, numFuncInTable_, funcIndex);
385 catch (gmx::GromacsException& ex)
387 ex.prependContext("Error generating cubic spline table for function '"
388 + thisFuncInput.desc + "'");
395 CubicSplineTable::CubicSplineTable(std::initializer_list<NumericalSplineTableInput> numericalInputList,
396 const std::pair<real, real>& range,
398 numFuncInTable_(numericalInputList.size()),
401 // Sanity check on input values
402 if (range.first < 0.0 || (range.second - range.first) < 0.001)
404 GMX_THROW(InvalidInputError(
405 "Range to tabulate cannot include negative values and must span at least 0.001"));
408 if (tolerance < GMX_REAL_EPS)
410 GMX_THROW(ToleranceError("Table tolerance cannot be smaller than GMX_REAL_EPS"));
413 double minQuotient = GMX_REAL_MAX;
415 // loop over all functions to find smallest spacing
416 for (auto thisFuncInput : numericalInputList)
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)
425 InconsistentInputError("Table input vectors must cover requested range, "
426 "and a margin beyond the upper endpoint"));
429 if (thisFuncInput.function.size() != thisFuncInput.derivative.size())
431 GMX_THROW(InconsistentInputError(
432 "Function and derivative vectors have different lengths"));
435 internal::throwUnlessDerivativeIsConsistentWithFunction(
436 thisFuncInput.function, thisFuncInput.derivative, thisFuncInput.spacing, range_);
438 catch (gmx::GromacsException& ex)
440 ex.prependContext("Error generating cubic spline table for function '"
441 + thisFuncInput.desc + "'");
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.
449 double thisMinQuotient = internal::findSmallestQuotientOfFunctionAndThirdDerivative(
450 thisFuncInput.derivative, thisFuncInput.spacing, range_);
452 minQuotient = std::min(minQuotient, thisMinQuotient);
455 double spacing = std::cbrt(72.0 * std::sqrt(3.0) * tolerance * minQuotient);
457 tableScale_ = 1.0 / spacing;
459 if (range_.second * tableScale_ > 1e6)
462 ToleranceError("Requested tolerance would require over a million points in table"));
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;
470 for (auto thisFuncInput : numericalInputList)
474 if (spacing < thisFuncInput.spacing)
477 ToleranceError("Input vector spacing cannot achieve tolerance requested"));
480 std::vector<real> tmpYfghTableData;
482 fillSingleCubicSplineTableData(thisFuncInput.function,
483 thisFuncInput.derivative,
484 thisFuncInput.spacing,
489 internal::fillMultiplexedTableData(
490 tmpYfghTableData, &yfghMultiTableData_, 4, numFuncInTable_, funcIndex);
494 catch (gmx::GromacsException& ex)
496 ex.prependContext("Error generating cubic spline table for function '"
497 + thisFuncInput.desc + "'");