2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2016,2017,2018,2019, 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 functiojn
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 std::size_t lastIndexInRange = endIndex - 1;
176 for (int i = endIndex - 1; 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 double lastIndexFunction = (*functionTableData)[lastIndexInRange];
224 double lastIndexDerivative = (*derivativeTableData)[lastIndexInRange];
225 (*functionTableData)[i] =
226 lastIndexFunction + lastIndexDerivative * (i - lastIndexInRange) * spacing;
227 (*derivativeTableData)[i] = lastIndexDerivative;
232 /*! \brief Create merged DDFZ vector from function & derivative data
234 * \param functionTableData Function values
235 * \param derivativeTableData Derivative values. We have already subtracted the
236 * small third derivative component when calling this
237 * function, but in practice it is just an arbitrary
239 * \param ddfzTableData Vector four times longer, filled with
240 * the derivative, the difference to the next derivative
241 * point, the function value, and zero.
243 * \throws If the vector lengths do not match.
245 void fillDdfzTableData(const std::vector<real>& functionTableData,
246 const std::vector<real>& derivativeTableData,
247 std::vector<real>* ddfzTableData)
249 GMX_ASSERT(functionTableData.size() == derivativeTableData.size(),
250 "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()),
280 // Sanity check on input values
281 if (range_.first < 0.0 || (range_.second - range_.first) < 0.001)
283 GMX_THROW(InvalidInputError(
284 "Range to tabulate cannot include negative values and must span at least 0.001"));
287 if (tolerance < GMX_REAL_EPS)
289 GMX_THROW(ToleranceError("Table tolerance cannot be smaller than GMX_REAL_EPS"));
292 double minQuotient = GMX_REAL_MAX;
294 // loop over all functions to find smallest spacing
295 for (const auto& thisFuncInput : analyticalInputList)
299 internal::throwUnlessDerivativeIsConsistentWithFunction(
300 thisFuncInput.function, thisFuncInput.derivative, range_);
302 catch (gmx::GromacsException& ex)
304 ex.prependContext("Error generating quadratic spline table for function '"
305 + thisFuncInput.desc + "'");
308 // Calculate the required table spacing h. The error we make with linear interpolation
309 // of the derivative will be described by the third-derivative correction term.
310 // This means we can compute the required spacing as h = sqrt(12*tolerance*min(f'/f''')),
311 // where f'/f''' is the first and third derivative of the function, respectively.
313 double thisMinQuotient = internal::findSmallestQuotientOfFunctionAndSecondDerivative(
314 thisFuncInput.derivative, range_);
316 minQuotient = std::min(minQuotient, thisMinQuotient);
319 double spacing = std::sqrt(12.0 * tolerance * minQuotient);
321 halfSpacing_ = 0.5 * spacing;
322 tableScale_ = 1.0 / spacing;
324 if (range_.second * tableScale_ > 1e6)
327 ToleranceError("Over a million points would be required for table; decrease range "
328 "or increase tolerance"));
331 // Loop over all tables again.
332 // Here we create the actual table for each function, and then
333 // combine them into a multiplexed table function.
334 std::size_t funcIndex = 0;
336 for (const auto& thisFuncInput : analyticalInputList)
340 std::vector<real> tmpFuncTableData;
341 std::vector<real> tmpDerTableData;
342 std::vector<real> tmpDdfzTableData;
344 fillSingleQuadraticSplineTableData(thisFuncInput.function, thisFuncInput.derivative,
345 range_, spacing, &tmpFuncTableData, &tmpDerTableData);
347 fillDdfzTableData(tmpFuncTableData, tmpDerTableData, &tmpDdfzTableData);
349 internal::fillMultiplexedTableData(tmpDerTableData, &derivativeMultiTableData_, 1,
350 numFuncInTable_, funcIndex);
352 internal::fillMultiplexedTableData(tmpDdfzTableData, &ddfzMultiTableData_, 4,
353 numFuncInTable_, funcIndex);
357 catch (gmx::GromacsException& ex)
359 ex.prependContext("Error generating quadratic spline table for function '"
360 + thisFuncInput.desc + "'");
367 QuadraticSplineTable::QuadraticSplineTable(std::initializer_list<NumericalSplineTableInput> numericalInputList,
368 const std::pair<real, real>& range,
370 numFuncInTable_(numericalInputList.size()),
373 // Sanity check on input values
374 if (range.first < 0.0 || (range.second - range.first) < 0.001)
376 GMX_THROW(InvalidInputError(
377 "Range to tabulate cannot include negative values and must span at least 0.001"));
380 if (tolerance < GMX_REAL_EPS)
382 GMX_THROW(ToleranceError("Table tolerance cannot be smaller than GMX_REAL_EPS"));
385 double minQuotient = GMX_REAL_MAX;
387 // loop over all functions to find smallest spacing
388 for (auto thisFuncInput : numericalInputList)
392 // We do not yet know what the margin is, but we need to test that we at least cover
393 // the requested range before starting to calculate derivatives
394 if (thisFuncInput.function.size() < range_.second / thisFuncInput.spacing + 1)
397 InconsistentInputError("Table input vectors must cover requested range, "
398 "and a margin beyond the upper endpoint"));
401 if (thisFuncInput.function.size() != thisFuncInput.derivative.size())
403 GMX_THROW(InconsistentInputError(
404 "Function and derivative vectors have different lengths"));
407 internal::throwUnlessDerivativeIsConsistentWithFunction(
408 thisFuncInput.function, thisFuncInput.derivative, thisFuncInput.spacing, range_);
410 catch (gmx::GromacsException& ex)
412 ex.prependContext("Error generating quadratic spline table for function '"
413 + thisFuncInput.desc + "'");
416 // Calculate the required table spacing h. The error we make with linear interpolation
417 // of the derivative will be described by the third-derivative correction term.
418 // This means we can compute the required spacing as h = sqrt(12*tolerance*min(f'/f''')),
419 // where f'/f''' is the first and third derivative of the function, respectively.
420 // Since we already have an analytical form of the derivative, we reduce the numerical
421 // errors by calculating the quotient of the function and second derivative of the
422 // input-derivative-analytical function instead.
424 double thisMinQuotient = internal::findSmallestQuotientOfFunctionAndSecondDerivative(
425 thisFuncInput.derivative, thisFuncInput.spacing, range_);
427 minQuotient = std::min(minQuotient, thisMinQuotient);
430 double spacing = std::sqrt(12.0 * tolerance * minQuotient);
432 halfSpacing_ = 0.5 * spacing;
433 tableScale_ = 1.0 / spacing;
435 if (range_.second * tableScale_ > 1e6)
438 ToleranceError("Requested tolerance would require over a million points in table"));
441 // Loop over all tables again.
442 // Here we create the actual table for each function, and then
443 // combine them into a multiplexed table function.
444 std::size_t funcIndex = 0;
446 for (auto thisFuncInput : numericalInputList)
450 if (spacing < thisFuncInput.spacing)
453 ToleranceError("Input vector spacing cannot achieve tolerance requested"));
456 std::vector<real> tmpFuncTableData;
457 std::vector<real> tmpDerTableData;
458 std::vector<real> tmpDdfzTableData;
460 fillSingleQuadraticSplineTableData(thisFuncInput.function, thisFuncInput.derivative,
461 thisFuncInput.spacing, range, spacing,
462 &tmpFuncTableData, &tmpDerTableData);
464 fillDdfzTableData(tmpFuncTableData, tmpDerTableData, &tmpDdfzTableData);
466 internal::fillMultiplexedTableData(tmpDerTableData, &derivativeMultiTableData_, 1,
467 numFuncInTable_, funcIndex);
469 internal::fillMultiplexedTableData(tmpDdfzTableData, &ddfzMultiTableData_, 4,
470 numFuncInTable_, funcIndex);
474 catch (gmx::GromacsException& ex)
476 ex.prependContext("Error generating quadratic spline table for function '"
477 + thisFuncInput.desc + "'");