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 internal utility functions for spline tables
40 * \author Erik Lindahl <erik.lindahl@gmail.com>
41 * \ingroup module_tables
45 #include "splineutil.h"
54 #include "gromacs/utility/alignedallocator.h"
55 #include "gromacs/utility/arrayref.h"
56 #include "gromacs/utility/exceptions.h"
57 #include "gromacs/utility/real.h"
58 #include "gromacs/utility/stringutil.h"
66 void throwUnlessDerivativeIsConsistentWithFunction(const std::function<double(double)>& function,
67 const std::function<double(double)>& derivative,
68 const std::pair<real, real>& range)
70 // Since the numerical derivative will evaluate extra points
71 // we shrink the interval slightly to avoid calling the function with values
72 // outside the range specified.
73 double h = std::cbrt(GMX_DOUBLE_EPS); // ideal spacing
74 std::pair<double, double> newRange(range.first + h, range.second - h);
75 const int points = 1000; // arbitrary
76 double dx = (newRange.second - newRange.first) / points;
77 bool isConsistent = true;
78 double minFail = newRange.second;
79 double maxFail = newRange.first;
81 // NOLINTNEXTLINE(clang-analyzer-security.FloatLoopCounter)
82 for (double x = newRange.first; x <= newRange.second; x += dx)
84 double analyticalDerivative = derivative(x);
85 double numericalDerivative = (function(x + h) - function(x - h)) / (2 * h);
86 double thirdDerivative = (derivative(x + h) - 2.0 * derivative(x) + derivative(x - h)) / (h * h);
88 // We make two types of errors in numerical calculation of the derivative:
89 // - The truncation error: eps * |f| / h
90 // - The rounding error: h * h * |f'''| / 6.0
92 GMX_DOUBLE_EPS * std::abs(function(x)) / h + h * h * std::abs(thirdDerivative) / 6.0;
94 // To avoid triggering false errors because of compiler optimization or numerical issues
95 // in the function evalulation we allow an extra factor of 10 in the expected error
96 if (std::abs(analyticalDerivative - numericalDerivative) > 10.0 * expectedErr)
99 minFail = std::min(minFail, x);
100 maxFail = std::max(maxFail, x);
106 GMX_THROW(InconsistentInputError(formatString(
107 "Derivative inconsistent with analytical function in range [%f,%f]", minFail, maxFail)));
112 void throwUnlessDerivativeIsConsistentWithFunction(ArrayRef<const double> function,
113 ArrayRef<const double> derivative,
115 const std::pair<real, real>& range)
117 std::size_t firstIndex = static_cast<std::size_t>(range.first / inputSpacing);
118 std::size_t lastIndex = static_cast<std::size_t>(range.second / inputSpacing);
119 bool isConsistent = true;
120 std::size_t minFail = lastIndex;
121 std::size_t maxFail = firstIndex;
123 // The derivative will access one extra point before/after each point, so reduce interval
124 for (std::size_t i = firstIndex + 1; (i + 1) < lastIndex; i++)
126 double inputDerivative = derivative[i];
127 double numericalDerivative = (function[i + 1] - function[i - 1]) / (2.0 * inputSpacing);
128 double thirdDerivative = (derivative[i + 1] - 2.0 * derivative[i] + derivative[i - 1])
129 / (inputSpacing * inputSpacing);
131 // We make two types of errors in numerical calculation of the derivative:
132 // - The truncation error: eps * |f| / h
133 // - The rounding error: h * h * |f'''| / 6.0
134 double expectedErr = GMX_DOUBLE_EPS * std::abs(function[i]) / inputSpacing
135 + inputSpacing * inputSpacing * std::abs(thirdDerivative) / 6.0;
137 // To avoid triggering false errors because of compiler optimization or numerical issues
138 // in the function evalulation we allow an extra factor of 10 in the expected error
139 if (std::abs(inputDerivative - numericalDerivative) > 10.0 * expectedErr)
141 isConsistent = false;
142 minFail = std::min(minFail, i);
143 maxFail = std::max(maxFail, i);
148 GMX_THROW(InconsistentInputError(
149 formatString("Derivative inconsistent with numerical vector for elements %zu-%zu",
150 minFail + 1, maxFail + 1)));
155 /*! \brief Calculate absolute quotient of function and its second derivative
157 * This is a utility function used in the functions to find the smallest quotient
160 * \param[in] previousPoint Value of function at x-h.
161 * \param[in] thisPoint Value of function at x.
162 * \param[in] nextPoint Value of function at x+h.
163 * \param[in] spacing Value of h.
165 * \return The absolute value of the quotient. If either the function or second
166 * derivative is smaller than sqrt(GMX_REAL_MIN), they will be set to
169 static double quotientOfFunctionAndSecondDerivative(double previousPoint,
174 double lowerLimit = static_cast<double>(std::sqrt(GMX_REAL_MIN));
175 double value = std::max(std::abs(thisPoint), lowerLimit);
176 double secondDerivative =
177 std::abs((previousPoint - 2.0 * thisPoint + nextPoint) / (spacing * spacing));
179 // Make sure we do not divide by zero. This limit is arbitrary,
180 // but it doesnt matter since this point will have a very large value,
181 // and the whole routine is searching for the smallest value.
182 secondDerivative = std::max(secondDerivative, lowerLimit);
184 return (value / secondDerivative);
188 real findSmallestQuotientOfFunctionAndSecondDerivative(const std::function<double(double)>& f,
189 const std::pair<real, real>& range)
191 // Since the numerical second derivative will evaluate extra points
192 // we shrink the interval slightly to avoid calling the function with values
193 // outside the range specified.
194 double h = std::pow(GMX_DOUBLE_EPS, 0.25);
195 std::pair<double, double> newRange(range.first + h, range.second - h);
196 const int points = 500; // arbitrary
197 double dx = (newRange.second - newRange.first) / points;
198 double minQuotient = GMX_REAL_MAX;
200 // NOLINTNEXTLINE(clang-analyzer-security.FloatLoopCounter)
201 for (double x = newRange.first; x <= newRange.second; x += dx)
203 minQuotient = std::min(minQuotient,
204 quotientOfFunctionAndSecondDerivative(f(x - h), f(x), f(x + h), h));
207 return static_cast<real>(minQuotient);
211 real findSmallestQuotientOfFunctionAndSecondDerivative(ArrayRef<const double> function,
213 const std::pair<real, real>& range)
216 std::size_t firstIndex = static_cast<std::size_t>(range.first / inputSpacing);
217 std::size_t lastIndex = static_cast<std::size_t>(range.second / inputSpacing);
218 double minQuotient = GMX_REAL_MAX;
220 for (std::size_t i = firstIndex + 1; (i + 1) < lastIndex; i++)
222 minQuotient = std::min(
223 minQuotient, quotientOfFunctionAndSecondDerivative(function[i - 1], function[i],
224 function[i + 1], inputSpacing));
226 return static_cast<real>(minQuotient);
230 /*! \brief Calculate absolute quotient of function and its third derivative
232 * This is a utility function used in the functions to find the smallest quotient
235 * \param[in] previousPreviousPoint Value of function at x-2h.
236 * \param[in] previousPoint Value of function at x-h.
237 * \param[in] thisPoint Value of function at x.
238 * \param[in] nextPoint Value of function at x+h.
239 * \param[in] nextNextPoint Value of function at x+2h.
240 * \param[in] spacing Value of h.
242 * \return The absolute value of the quotient. If either the function or third
243 * derivative is smaller than sqrt(GMX_REAL_MIN), they will be set to
246 static double quotientOfFunctionAndThirdDerivative(double previousPreviousPoint,
247 double previousPoint,
250 double nextNextPoint,
253 double lowerLimit = static_cast<double>(std::sqrt(GMX_REAL_MIN));
254 double value = std::max(std::abs(thisPoint), lowerLimit);
255 double thirdDerivative =
256 std::abs((nextNextPoint - 2 * nextPoint + 2 * previousPoint - previousPreviousPoint)
257 / (2 * spacing * spacing * spacing));
259 // Make sure we do not divide by zero. This limit is arbitrary,
260 // but it doesnt matter since this point will have a very large value,
261 // and the whole routine is searching for the smallest value.
262 thirdDerivative = std::max(thirdDerivative, lowerLimit);
264 return (value / thirdDerivative);
268 real findSmallestQuotientOfFunctionAndThirdDerivative(const std::function<double(double)>& f,
269 const std::pair<real, real>& range)
271 // Since the numerical second derivative will evaluate extra points
272 // we shrink the interval slightly to avoid calling the function with values
273 // outside the range specified.
274 double h = std::pow(GMX_DOUBLE_EPS, 0.2); // optimal spacing for 3rd derivative
275 std::pair<double, double> newRange(range.first + 2 * h, range.second - 2 * h);
276 const int points = 500; // arbitrary
277 double dx = (newRange.second - newRange.first) / points;
278 double minQuotient = GMX_REAL_MAX;
280 // NOLINTNEXTLINE(clang-analyzer-security.FloatLoopCounter)
281 for (double x = newRange.first; x <= newRange.second; x += dx)
283 minQuotient = std::min(minQuotient,
284 quotientOfFunctionAndThirdDerivative(f(x - 2 * h), f(x - h), f(x),
285 f(x + h), f(x + 2 * h), h));
287 return static_cast<real>(minQuotient);
291 real findSmallestQuotientOfFunctionAndThirdDerivative(ArrayRef<const double> function,
293 const std::pair<real, real>& range)
296 std::size_t firstIndex = static_cast<std::size_t>(range.first / inputSpacing);
297 std::size_t lastIndex = static_cast<std::size_t>(range.second / inputSpacing);
298 double minQuotient = GMX_REAL_MAX;
300 for (std::size_t i = firstIndex + 2; (i + 2) < lastIndex; i++)
302 minQuotient = std::min(minQuotient, quotientOfFunctionAndThirdDerivative(
303 function[i - 2], function[i - 1], function[i],
304 function[i + 1], function[i + 2], inputSpacing));
306 return static_cast<real>(minQuotient);
310 std::vector<double> vectorSecondDerivative(ArrayRef<const double> f, double spacing)
314 GMX_THROW(APIError("Too few points in vector for 5-point differentiation"));
317 std::vector<double> d(f.size());
320 // 5-point formula evaluated for points 0,1
322 d[i] = (11 * f[i + 4] - 56 * f[i + 3] + 114 * f[i + 2] - 104 * f[i + 1] + 35 * f[i])
323 / (12 * spacing * spacing);
325 d[i] = (-f[i + 3] + 4 * f[i + 2] + 6 * f[i + 1] - 20 * f[i] + 11 * f[i - 1]) / (12 * spacing * spacing);
327 for (std::size_t i = 2; i < d.size() - 2; i++)
329 // 5-point formula evaluated for central point (2)
330 d[i] = (-f[i + 2] + 16 * f[i + 1] - 30 * f[i] + 16 * f[i - 1] - f[i - 2]) / (12 * spacing * spacing);
333 // 5-point formula evaluated for points 3,4
335 d[i] = (11 * f[i + 1] - 20 * f[i] + 6 * f[i - 1] + 4 * f[i - 2] - f[i - 3]) / (12 * spacing * spacing);
337 d[i] = (35 * f[i] - 104 * f[i - 1] + 114 * f[i - 2] - 56 * f[i - 3] + 11 * f[i - 4])
338 / (12 * spacing * spacing);
344 } // namespace internal