2 * This file is part of the GROMACS molecular simulation package.
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.
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()), range_(range)
311 // Sanity check on input values
312 if (range_.first < 0.0 || (range_.second - range_.first) < 0.001)
314 GMX_THROW(InvalidInputError(
315 "Range to tabulate cannot include negative values and must span at least 0.001"));
318 if (tolerance < GMX_REAL_EPS)
320 GMX_THROW(ToleranceError("Table tolerance cannot be smaller than GMX_REAL_EPS"));
323 double minQuotient = GMX_REAL_MAX;
325 // loop over all functions to find smallest spacing
326 for (const auto& thisFuncInput : analyticalInputList)
330 internal::throwUnlessDerivativeIsConsistentWithFunction(
331 thisFuncInput.function, thisFuncInput.derivative, range_);
333 catch (gmx::GromacsException& ex)
335 ex.prependContext("Error generating cubic spline table for function '"
336 + thisFuncInput.desc + "'");
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.
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.
348 double thisMinQuotient = internal::findSmallestQuotientOfFunctionAndThirdDerivative(
349 thisFuncInput.derivative, range_);
351 minQuotient = std::min(minQuotient, thisMinQuotient);
354 double spacing = 0.5 * std::cbrt(72.0 * std::sqrt(3.0) * tolerance * minQuotient);
356 tableScale_ = 1.0 / spacing;
358 if (range_.second * tableScale_ > 2e6)
361 ToleranceError("Over a million points would be required for table; decrease range "
362 "or increase tolerance"));
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;
370 for (const auto& thisFuncInput : analyticalInputList)
374 std::vector<real> tmpYfghTableData;
376 fillSingleCubicSplineTableData(
377 thisFuncInput.function, thisFuncInput.derivative, range_, spacing, &tmpYfghTableData);
379 internal::fillMultiplexedTableData(
380 tmpYfghTableData, &yfghMultiTableData_, 4, numFuncInTable_, funcIndex);
384 catch (gmx::GromacsException& ex)
386 ex.prependContext("Error generating cubic spline table for function '"
387 + thisFuncInput.desc + "'");
394 CubicSplineTable::CubicSplineTable(std::initializer_list<NumericalSplineTableInput> numericalInputList,
395 const std::pair<real, real>& range,
397 numFuncInTable_(numericalInputList.size()), range_(range)
399 // Sanity check on input values
400 if (range.first < 0.0 || (range.second - range.first) < 0.001)
402 GMX_THROW(InvalidInputError(
403 "Range to tabulate cannot include negative values and must span at least 0.001"));
406 if (tolerance < GMX_REAL_EPS)
408 GMX_THROW(ToleranceError("Table tolerance cannot be smaller than GMX_REAL_EPS"));
411 double minQuotient = GMX_REAL_MAX;
413 // loop over all functions to find smallest spacing
414 for (auto thisFuncInput : numericalInputList)
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)
423 InconsistentInputError("Table input vectors must cover requested range, "
424 "and a margin beyond the upper endpoint"));
427 if (thisFuncInput.function.size() != thisFuncInput.derivative.size())
429 GMX_THROW(InconsistentInputError(
430 "Function and derivative vectors have different lengths"));
433 internal::throwUnlessDerivativeIsConsistentWithFunction(
434 thisFuncInput.function, thisFuncInput.derivative, thisFuncInput.spacing, range_);
436 catch (gmx::GromacsException& ex)
438 ex.prependContext("Error generating cubic spline table for function '"
439 + thisFuncInput.desc + "'");
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.
447 double thisMinQuotient = internal::findSmallestQuotientOfFunctionAndThirdDerivative(
448 thisFuncInput.derivative, thisFuncInput.spacing, range_);
450 minQuotient = std::min(minQuotient, thisMinQuotient);
453 double spacing = std::cbrt(72.0 * std::sqrt(3.0) * tolerance * minQuotient);
455 tableScale_ = 1.0 / spacing;
457 if (range_.second * tableScale_ > 1e6)
460 ToleranceError("Requested tolerance would require over a million points in table"));
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;
468 for (auto thisFuncInput : numericalInputList)
472 if (spacing < thisFuncInput.spacing)
475 ToleranceError("Input vector spacing cannot achieve tolerance requested"));
478 std::vector<real> tmpYfghTableData;
480 fillSingleCubicSplineTableData(thisFuncInput.function,
481 thisFuncInput.derivative,
482 thisFuncInput.spacing,
487 internal::fillMultiplexedTableData(
488 tmpYfghTableData, &yfghMultiTableData_, 4, numFuncInTable_, funcIndex);
492 catch (gmx::GromacsException& ex)
494 ex.prependContext("Error generating cubic spline table for function '"
495 + thisFuncInput.desc + "'");