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 quadratic spline table functions
40 * \author Erik Lindahl <erik.lindahl@gmail.com>
41 * \ingroup module_tables
45 #include "quadraticsplinetable.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 Construct the data for a single quadratic table from analytical functions
71 * \param[in] function Analytical function
72 * \param[in] derivative Analytical derivative
73 * \param[in] range Upper/lower limit of region to tabulate
74 * \param[in] spacing Distance between table points
75 * \param[out] functionTableData Output table with function data
76 * \param[out] derivativeTableData OUtput table with (adjusted) derivative data
78 void fillSingleQuadraticSplineTableData(const std::function<double(double)>& function,
79 const std::function<double(double)>& derivative,
80 const std::pair<real, real>& range,
82 std::vector<real>* functionTableData,
83 std::vector<real>* derivativeTableData)
85 std::size_t endIndex = static_cast<std::size_t>(range.second / spacing + 2);
87 functionTableData->resize(endIndex);
88 derivativeTableData->resize(endIndex);
90 double maxMagnitude = 0.0001 * GMX_REAL_MAX;
91 bool functionIsInRange = true;
92 std::size_t lastIndexInRange = endIndex - 1;
94 for (int i = endIndex - 1; i >= 0; i--)
96 double x = i * spacing;
97 double tmpFunctionValue;
98 double tmpDerivativeValue;
100 if (range.first > 0 && i == 0)
102 // Avoid x==0 if it is not in the range, since it can lead to
103 // singularities even if the value for i==1 was within or required magnitude
104 functionIsInRange = false;
107 if (functionIsInRange)
109 tmpFunctionValue = function(x);
111 // Calculate third derivative term (2nd derivative of the derivative)
112 // Make sure we stay in range. In practice this means we use one-sided
113 // interpolation at the interval endpoints (indentical to an offset for 3-point formula)
114 const double h = std::pow(GMX_DOUBLE_EPS, 0.25);
115 double y = std::min(std::max(x, range.first + h), range.second - h);
116 double thirdDerivativeValue =
117 (derivative(y + h) - 2.0 * derivative(y) + derivative(y - h)) / (h * h);
119 tmpDerivativeValue = derivative(x) - spacing * spacing * thirdDerivativeValue / 12.0;
121 if (std::abs(tmpFunctionValue) > maxMagnitude || std::abs(tmpDerivativeValue) > maxMagnitude)
123 functionIsInRange = false; // Once this happens, it never resets to true again
127 if (functionIsInRange)
129 (*functionTableData)[i] = tmpFunctionValue;
130 (*derivativeTableData)[i] = tmpDerivativeValue;
135 // Once the function or derivative (more likely) has reached very large values,
136 // we simply make a linear function from the last in-range value of the derivative.
137 double lastIndexFunction = (*functionTableData)[lastIndexInRange];
138 double lastIndexDerivative = (*derivativeTableData)[lastIndexInRange];
139 (*functionTableData)[i] =
140 lastIndexFunction + lastIndexDerivative * (i - lastIndexInRange) * spacing;
141 (*derivativeTableData)[i] = lastIndexDerivative;
147 /*! \brief Construct the data for a single quadratic table from vector data
149 * \param[in] function Input vector with function data
150 * \param[in] derivative Input vector with derivative data
151 * \param[in] inputSpacing Distance between points in input vectors
152 * \param[in] range Upper/lower limit of region to tabulate
153 * \param[in] spacing Distance between table points
154 * \param[out] functionTableData Output table with function data
155 * \param[out] derivativeTableData OUtput table with (adjusted) derivative data
157 void fillSingleQuadraticSplineTableData(ArrayRef<const double> function,
158 ArrayRef<const double> derivative,
160 const std::pair<real, real>& range,
162 std::vector<real>* functionTableData,
163 std::vector<real>* derivativeTableData)
165 std::size_t endIndex = static_cast<std::size_t>(range.second / spacing + 2);
167 functionTableData->resize(endIndex);
168 derivativeTableData->resize(endIndex);
170 std::vector<double> thirdDerivative(internal::vectorSecondDerivative(derivative, inputSpacing));
172 double maxMagnitude = 0.0001 * GMX_REAL_MAX;
173 bool functionIsInRange = true;
174 int lastIndexInRange = static_cast<int>(endIndex) - 1;
176 for (int i = lastIndexInRange; i >= 0; i--)
178 double x = i * spacing;
179 double tmpFunctionValue;
180 double tmpDerivativeValue;
182 if (range.first > 0 && i == 0)
184 // Avoid x==0 if it is not in the range, since it can lead to
185 // singularities even if the value for i==1 was within or required magnitude
186 functionIsInRange = false;
189 if (functionIsInRange)
191 // Step 1: Interpolate the function value at x from input table.
192 double inputXTab = x / inputSpacing;
193 int inputIndex = static_cast<std::size_t>(inputXTab);
194 double inputEps = inputXTab - inputIndex;
196 // Linear interpolation of input derivative and third derivative
197 double thirdDerivativeValue = (1.0 - inputEps) * thirdDerivative[inputIndex]
198 + inputEps * thirdDerivative[inputIndex + 1];
199 double derivativeValue =
200 (1.0 - inputEps) * derivative[inputIndex] + inputEps * derivative[inputIndex + 1];
202 // Quadratic interpolation for function value
203 tmpFunctionValue = function[inputIndex]
204 + 0.5 * (derivative[inputIndex] + derivativeValue) * inputEps * inputSpacing;
205 tmpDerivativeValue = derivativeValue - spacing * spacing * thirdDerivativeValue / 12.0;
207 if (std::abs(tmpFunctionValue) > maxMagnitude || std::abs(tmpDerivativeValue) > maxMagnitude)
209 functionIsInRange = false; // Once this happens, it never resets to true again
213 if (functionIsInRange)
215 (*functionTableData)[i] = tmpFunctionValue;
216 (*derivativeTableData)[i] = tmpDerivativeValue;
221 // Once the function or derivative (more likely) has reached very large values,
222 // we simply make a linear function from the last in-range value of the derivative.
223 GMX_ASSERT(lastIndexInRange >= 0, "Array index is unexpectedly negative.");
224 double lastIndexFunction = (*functionTableData)[lastIndexInRange];
225 double lastIndexDerivative = (*derivativeTableData)[lastIndexInRange];
226 (*functionTableData)[i] =
227 lastIndexFunction + lastIndexDerivative * (i - lastIndexInRange) * spacing;
228 (*derivativeTableData)[i] = lastIndexDerivative;
233 /*! \brief Create merged DDFZ vector from function & derivative data
235 * \param functionTableData Function values
236 * \param derivativeTableData Derivative values. We have already subtracted the
237 * small third derivative component when calling this
238 * function, but in practice it is just an arbitrary
240 * \param ddfzTableData Vector four times longer, filled with
241 * the derivative, the difference to the next derivative
242 * point, the function value, and zero.
244 * \throws If the vector lengths do not match.
246 void fillDdfzTableData(const std::vector<real>& functionTableData,
247 const std::vector<real>& derivativeTableData,
248 std::vector<real>* ddfzTableData)
250 GMX_ASSERT(functionTableData.size() == derivativeTableData.size(), "Mismatching vector lengths");
252 std::size_t points = functionTableData.size();
254 ddfzTableData->resize(4 * points);
256 for (std::size_t i = 0; i < points; i++)
258 (*ddfzTableData)[4 * i] = derivativeTableData[i];
260 double nextDerivative = (i < functionTableData.size() - 1) ? derivativeTableData[i + 1] : 0.0;
262 (*ddfzTableData)[4 * i + 1] = nextDerivative - derivativeTableData[i];
263 (*ddfzTableData)[4 * i + 2] = functionTableData[i];
264 (*ddfzTableData)[4 * i + 3] = 0.0;
271 const real QuadraticSplineTable::defaultTolerance = 10.0 * GMX_FLOAT_EPS;
274 QuadraticSplineTable::QuadraticSplineTable(std::initializer_list<AnalyticalSplineTableInput> analyticalInputList,
275 const std::pair<real, real>& range,
277 numFuncInTable_(analyticalInputList.size()), range_(range)
279 // Sanity check on input values
280 if (range_.first < 0.0 || (range_.second - range_.first) < 0.001)
282 GMX_THROW(InvalidInputError(
283 "Range to tabulate cannot include negative values and must span at least 0.001"));
286 if (tolerance < GMX_REAL_EPS)
288 GMX_THROW(ToleranceError("Table tolerance cannot be smaller than GMX_REAL_EPS"));
291 double minQuotient = GMX_REAL_MAX;
293 // loop over all functions to find smallest spacing
294 for (const auto& thisFuncInput : analyticalInputList)
298 internal::throwUnlessDerivativeIsConsistentWithFunction(
299 thisFuncInput.function, thisFuncInput.derivative, range_);
301 catch (gmx::GromacsException& ex)
303 ex.prependContext("Error generating quadratic spline table for function '"
304 + thisFuncInput.desc + "'");
307 // Calculate the required table spacing h. The error we make with linear interpolation
308 // of the derivative will be described by the third-derivative correction term.
309 // This means we can compute the required spacing as h = sqrt(12*tolerance*min(f'/f''')),
310 // where f'/f''' is the first and third derivative of the function, respectively.
312 double thisMinQuotient = internal::findSmallestQuotientOfFunctionAndSecondDerivative(
313 thisFuncInput.derivative, range_);
315 minQuotient = std::min(minQuotient, thisMinQuotient);
318 double spacing = std::sqrt(12.0 * tolerance * minQuotient);
320 halfSpacing_ = 0.5 * spacing;
321 tableScale_ = 1.0 / spacing;
323 if (range_.second * tableScale_ > 1e6)
326 ToleranceError("Over a million points would be required for table; decrease range "
327 "or increase tolerance"));
330 // Loop over all tables again.
331 // Here we create the actual table for each function, and then
332 // combine them into a multiplexed table function.
333 std::size_t funcIndex = 0;
335 for (const auto& thisFuncInput : analyticalInputList)
339 std::vector<real> tmpFuncTableData;
340 std::vector<real> tmpDerTableData;
341 std::vector<real> tmpDdfzTableData;
343 fillSingleQuadraticSplineTableData(thisFuncInput.function,
344 thisFuncInput.derivative,
350 fillDdfzTableData(tmpFuncTableData, tmpDerTableData, &tmpDdfzTableData);
352 internal::fillMultiplexedTableData(
353 tmpDerTableData, &derivativeMultiTableData_, 1, numFuncInTable_, funcIndex);
355 internal::fillMultiplexedTableData(
356 tmpDdfzTableData, &ddfzMultiTableData_, 4, numFuncInTable_, funcIndex);
360 catch (gmx::GromacsException& ex)
362 ex.prependContext("Error generating quadratic spline table for function '"
363 + thisFuncInput.desc + "'");
370 QuadraticSplineTable::QuadraticSplineTable(std::initializer_list<NumericalSplineTableInput> numericalInputList,
371 const std::pair<real, real>& range,
373 numFuncInTable_(numericalInputList.size()), range_(range)
375 // Sanity check on input values
376 if (range.first < 0.0 || (range.second - range.first) < 0.001)
378 GMX_THROW(InvalidInputError(
379 "Range to tabulate cannot include negative values and must span at least 0.001"));
382 if (tolerance < GMX_REAL_EPS)
384 GMX_THROW(ToleranceError("Table tolerance cannot be smaller than GMX_REAL_EPS"));
387 double minQuotient = GMX_REAL_MAX;
389 // loop over all functions to find smallest spacing
390 for (auto thisFuncInput : numericalInputList)
394 // We do not yet know what the margin is, but we need to test that we at least cover
395 // the requested range before starting to calculate derivatives
396 if (thisFuncInput.function.size() < range_.second / thisFuncInput.spacing + 1)
399 InconsistentInputError("Table input vectors must cover requested range, "
400 "and a margin beyond the upper endpoint"));
403 if (thisFuncInput.function.size() != thisFuncInput.derivative.size())
405 GMX_THROW(InconsistentInputError(
406 "Function and derivative vectors have different lengths"));
409 internal::throwUnlessDerivativeIsConsistentWithFunction(
410 thisFuncInput.function, thisFuncInput.derivative, thisFuncInput.spacing, range_);
412 catch (gmx::GromacsException& ex)
414 ex.prependContext("Error generating quadratic spline table for function '"
415 + thisFuncInput.desc + "'");
418 // Calculate the required table spacing h. The error we make with linear interpolation
419 // of the derivative will be described by the third-derivative correction term.
420 // This means we can compute the required spacing as h = sqrt(12*tolerance*min(f'/f''')),
421 // where f'/f''' is the first and third derivative of the function, respectively.
422 // Since we already have an analytical form of the derivative, we reduce the numerical
423 // errors by calculating the quotient of the function and second derivative of the
424 // input-derivative-analytical function instead.
426 double thisMinQuotient = internal::findSmallestQuotientOfFunctionAndSecondDerivative(
427 thisFuncInput.derivative, thisFuncInput.spacing, range_);
429 minQuotient = std::min(minQuotient, thisMinQuotient);
432 double spacing = std::sqrt(12.0 * tolerance * minQuotient);
434 halfSpacing_ = 0.5 * spacing;
435 tableScale_ = 1.0 / spacing;
437 if (range_.second * tableScale_ > 1e6)
440 ToleranceError("Requested tolerance would require over a million points in table"));
443 // Loop over all tables again.
444 // Here we create the actual table for each function, and then
445 // combine them into a multiplexed table function.
446 std::size_t funcIndex = 0;
448 for (auto thisFuncInput : numericalInputList)
452 if (spacing < thisFuncInput.spacing)
455 ToleranceError("Input vector spacing cannot achieve tolerance requested"));
458 std::vector<real> tmpFuncTableData;
459 std::vector<real> tmpDerTableData;
460 std::vector<real> tmpDdfzTableData;
462 fillSingleQuadraticSplineTableData(thisFuncInput.function,
463 thisFuncInput.derivative,
464 thisFuncInput.spacing,
470 fillDdfzTableData(tmpFuncTableData, tmpDerTableData, &tmpDdfzTableData);
472 internal::fillMultiplexedTableData(
473 tmpDerTableData, &derivativeMultiTableData_, 1, numFuncInTable_, funcIndex);
475 internal::fillMultiplexedTableData(
476 tmpDdfzTableData, &ddfzMultiTableData_, 4, numFuncInTable_, funcIndex);
480 catch (gmx::GromacsException& ex)
482 ex.prependContext("Error generating quadratic spline table for function '"
483 + thisFuncInput.desc + "'");