Introduce SimulatorBuilder
[gromacs.git] / src / gromacs / tables / splineutil.h
blob8a58486c9de071181c1cecb92d5f11b27bac61d9
1 /*
2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2016,2017,2018, 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.
36 /*! \internal \file
37 * \brief
38 * Declares internal utility functions for spline tables
40 * \author Erik Lindahl <erik.lindahl@gmail.com>
41 * \ingroup module_tables
43 #ifndef GMX_TABLES_SPLINEUTIL_H
44 #define GMX_TABLES_SPLINEUTIL_H
46 #include <functional>
47 #include <utility>
48 #include <vector>
50 #include "gromacs/utility/alignedallocator.h"
51 #include "gromacs/utility/arrayref.h"
52 #include "gromacs/utility/gmxassert.h"
53 #include "gromacs/utility/real.h"
55 namespace gmx
58 namespace internal
61 /*! \brief Ensure analytical derivative is the derivative of analytical function.
63 * This routine evaluates the numerical derivative of the function for
64 * a few (1000) points in the interval and checks that the relative difference
65 * between numerical and analytical derivative is within the expected error
66 * for the numerical derivative approximation we use.
68 * The main point of this routine is to make sure the user has not made a
69 * mistake or sign error when defining the functions.
71 * \param function Analytical function to differentiate
72 * \param derivative Analytical derivative to compare with
73 * \param range Range to test
75 * \throws If the provided derivative does not seem to match the function.
77 * \note The function/derivative are always double-valued to avoid accuracy loss.
79 void
80 throwUnlessDerivativeIsConsistentWithFunction(const std::function<double(double)> &function,
81 const std::function<double(double)> &derivative,
82 const std::pair<real, real> &range);
84 /*! \brief Ensure vector of derivative values is the derivative of function vector.
86 * This routine differentiates a vector of numerical values and checks
87 * that the relative difference to a provided vector of numerical derivatives
88 * is smaller than the expected error from the numerical differentiation.
90 * The main point of this routine is to make sure the user has not made a
91 * mistake or sign error when defining the functions.
93 * To avoid problems if the vectors change from zero to finite values at the
94 * start/end of the interval, we only check inside the range requested.
96 * \param function Numerical function value vector to differentiate
97 * \param derivative Numerical derivative vector to compare with
98 * \param inputSpacing Distance between input points
99 * \param range Range to test
101 * \throws If the provided derivative does not seem to match the function.
103 * \note The function/derivative vectors and spacing are always double-valued
104 * to avoid accuracy loss.
106 void
107 throwUnlessDerivativeIsConsistentWithFunction(ArrayRef<const double> function,
108 ArrayRef<const double> derivative,
109 double inputSpacing,
110 const std::pair<real, real> &range);
115 /*! \brief Find smallest quotient between analytical function and its 2nd derivative
117 * Used to calculate spacing for quadratic spline tables. This function divides the
118 * function value by the second derivative (or a very small number when that is zero),
119 * and returns the smallest such quotient found in the range.
121 * Our quadratic tables corresponds to linear interpolation of the derivative,
122 * which means the derivative will typically have larger error than the value
123 * when interpolating. The spacing required to reach a particular relative
124 * tolerance in the derivative depends on the quotient between the first
125 * derivative and the third derivative of the function itself.
127 * You should call this routine with the analytical derivative as the "function"
128 * parameter, and the quotient between "function and second derivative" will
129 * then correspond to the quotient bewteen the derivative and the third derivative
130 * of the actual function we want to tabulate.
132 * Since all functions that can be tabulated efficiently are reasonably smooth,
133 * we simply check 1,000 points in the interval rather than bother about
134 * implementing any complicated global optimization scheme.
136 * \param f Analytical function
137 * \param range Interval
139 * \return Smallest quotient found in range.
141 * \note The function is always double-valued to avoid accuracy loss.
143 real
144 findSmallestQuotientOfFunctionAndSecondDerivative(const std::function<double(double)> &f,
145 const std::pair<real, real> &range);
151 /*! \brief Find smallest quotient between vector of values and its 2nd derivative
153 * Used to calculate spacing for quadratic spline tables. This function divides the
154 * function value by the second derivative (or a very small number when that is zero),
155 * and returns the smallest such quotient found in the range.
157 * Our quadratic tables corresponds to linear interpolation of the derivative,
158 * which means the derivative will typically have larger error than the value
159 * when interpolating. The spacing required to reach a particular relative
160 * tolerance in the derivative depends on the quotient between the first
161 * derivative and the third derivative of the function itself.
163 * You should call this routine with the analytical derivative as the "function"
164 * parameter, and the quotient between "function and second derivative" will
165 * then correspond to the quotient bewteen the derivative and the third derivative
166 * of the actual function we want to tabulate.
168 * \param function Vector with function values
169 * \param inputSpacing Spacing between function values
170 * \param range Interval to check
172 * \return Smallest quotient found in range.
174 * \note The function vector and input spacing are always double-valued to
175 * avoid accuracy loss.
177 real
178 findSmallestQuotientOfFunctionAndSecondDerivative(ArrayRef<const double> function,
179 double inputSpacing,
180 const std::pair<real, real> &range);
186 /*! \brief Find smallest quotient between analytical function and its 3rd derivative
188 * Used to calculate table spacing. This function divides the function value
189 * by the second derivative (or a very small number when that is zero), and
190 * returns the smallest such quotient found in the range.
192 * Our quadratic tables corresponds to linear interpolation of the derivative,
193 * which means the derivative will typically have larger error than the value
194 * when interpolating. The spacing required to reach a particular relative
195 * tolerance in the derivative depends on the quotient between the first
196 * derivative and the third derivative of the function itself.
198 * You should call this routine with the analytical derivative as the "function"
199 * parameter, and the quotient between "function and second derivative" will
200 * then correspond to the quotient bewteen the derivative and the third derivative
201 * of the actual function we want to tabulate.
203 * Since all functions that can be tabulated efficiently are reasonably smooth,
204 * we simply check 1,000 points in the interval rather than bother about
205 * implementing any complicated global optimization scheme.
207 * \param f Analytical function
208 * \param range Interval
210 * \return Smallest quotient found in range.
212 * \note The function is always double-valued to avoid accuracy loss.
214 real
215 findSmallestQuotientOfFunctionAndThirdDerivative(const std::function<double(double)> &f,
216 const std::pair<real, real> &range);
221 /*! \brief Find smallest quotient between function and 2nd derivative (vectors)
223 * Used to calculate table spacing. This function divides the function value
224 * by the second derivative (or a very small number when that is zero), and
225 * returns the smallest such quotient found in the range.
227 * Our quadratic tables corresponds to linear interpolation of the derivative,
228 * which means the derivative will typically have larger error than the value
229 * when interpolating. The spacing required to reach a particular relative
230 * tolerance in the derivative depends on the quotient between the first
231 * derivative and the third derivative of the function itself.
233 * You should call this routine with the analytical derivative as the "function"
234 * parameter, and the quotient between "function and second derivative" will
235 * then correspond to the quotient bewteen the derivative and the third derivative
236 * of the actual function we want to tabulate.
238 * \param function Vector with function values
239 * \param inputSpacing Spacing between function values
240 * \param range Interval to check
242 * \return Smallest quotient found in range.
244 * \note The function vector and input spacing are always double-valued to
245 * avoid accuracy loss.
247 real
248 findSmallestQuotientOfFunctionAndThirdDerivative(ArrayRef<const double> function,
249 double inputSpacing,
250 const std::pair<real, real> &range);
253 /*! \brief Calculate second derivative of vector and return vector of same length
255 * 5-point approximations are used, with endpoints using non-center interpolation.
257 * \param f Vector (function) for which to calculate second derivative
258 * \param spacing Spacing of input data.
260 * \throws If the input vector has fewer than five data points.
262 * \note This function always uses double precision arguments since it is meant
263 * to be used on raw user input data for tables, where we want to avoid
264 * accuracy loss (since differentiation can be numerically fragile).
266 std::vector<double>
267 vectorSecondDerivative(ArrayRef<const double> f,
268 double spacing);
271 /*! \brief Copy (temporary) table data into aligned multiplexed vector
273 * This routine takes the temporary data generated for a single table
274 * and writes multiplexed output into a multiple-table-data vector.
275 * If the output vector is empty we will resize it to fit the data, and
276 * otherwise we assert the size is correct to add out input data.
278 * \tparam T Type of container for input data
279 * \tparam U Type of container for output data
281 * \param[in] inputData Input data for single table
282 * \param[inout] multiplexedOutputData Multiplexed output vector, many tables.
283 * \param[in] valuesPerTablePoint Number of real values for each table point,
284 * for instance 4 in DDFZ tables.
285 * \param[in] numTables Number of tables mixed into multiplexed output
286 * \param[in] thisTableIndex Index of this table in multiplexed output
288 * \note The output container type can be different from the input since the latter
289 * sometimes uses an aligned allocator so the data can be loaded efficiently
290 * in the GROMACS nonbonded kernels.
292 template<class T, class U>
293 void
294 fillMultiplexedTableData(const T inputData,
295 U * multiplexedOutputData,
296 std::size_t valuesPerTablePoint,
297 std::size_t numTables,
298 std::size_t thisTableIndex)
300 if (multiplexedOutputData->empty())
302 multiplexedOutputData->resize( inputData.size() * numTables );
304 else
306 GMX_ASSERT(inputData.size() * numTables == multiplexedOutputData->size(),
307 "Size mismatch when multiplexing table data");
310 GMX_ASSERT(inputData.size() % valuesPerTablePoint == 0,
311 "Input table size must be a multiple of valuesPerTablePoint");
313 std::size_t points = inputData.size() / valuesPerTablePoint;
315 for (std::size_t i = 0; i < points; i++)
317 std::size_t inputOffset = valuesPerTablePoint * i;
318 std::size_t outputOffset = valuesPerTablePoint * ( numTables * i + thisTableIndex );
320 for (std::size_t j = 0; j < valuesPerTablePoint; j++)
322 (*multiplexedOutputData)[outputOffset + j] = inputData[inputOffset + j];
328 } // namespace internal
330 } // namespace gmx
332 #endif // GMX_TABLES_SPLINEUTIL_H