src/gromacs/tables/splineutil.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
   8  *
   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.
  13  *
  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.
  18  *
  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.
  23  *
  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.
  31  *
  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.
  34  */
  35
  36 /*! \internal \file
  37  * \brief
  38  * Declares internal utility functions for spline tables
  39  *
  40  * \author Erik Lindahl <erik.lindahl@gmail.com>
  41  * \ingroup module_tables
  42  */
  43 #ifndef GMX_TABLES_SPLINEUTIL_H
  44 #define GMX_TABLES_SPLINEUTIL_H
  45
  46 #include <functional>
  47 #include <utility>
  48 #include <vector>
  49
  50 #include "gromacs/utility/alignedallocator.h"
  51 #include "gromacs/utility/arrayref.h"
  52 #include "gromacs/utility/gmxassert.h"
  53 #include "gromacs/utility/real.h"
  54
  55 namespace gmx
  56 {
  57
  58 namespace internal
  59 {
  60
  61 /*! \brief Ensure analytical derivative is the derivative of analytical function.
  62  *
  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.
  67  *
  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.
  70  *
  71  *  \param function   Analytical function to differentiate
  72  *  \param derivative Analytical derivative to compare with
  73  *  \param range      Range to test
  74  *
  75  *  \throws If the provided derivative does not seem to match the function.
  76  *
  77  *  \note The function/derivative are always double-valued to avoid accuracy loss.
  78  */
  79 void throwUnlessDerivativeIsConsistentWithFunction(const std::function<double(double)>& function,
  80                                                    const std::function<double(double)>& derivative,
  81                                                    const std::pair<real, real>&         range);
  82
  83 /*! \brief Ensure vector of derivative values is the derivative of function vector.
  84  *
  85  *  This routine differentiates a vector of numerical values and checks
  86  *  that the relative difference to a provided vector of numerical derivatives
  87  *  is smaller than the expected error from the numerical differentiation.
  88  *
  89  *  The main point of this routine is to make sure the user has not made a
  90  *  mistake or sign error when defining the functions.
  91  *
  92  *  To avoid problems if the vectors change from zero to finite values at the
  93  *  start/end of the interval, we only check inside the range requested.
  94  *
  95  *  \param function     Numerical function value vector to differentiate
  96  *  \param derivative   Numerical derivative vector to compare with
  97  *  \param inputSpacing Distance between input points
  98  *  \param range        Range to test
  99  *
 100  *  \throws If the provided derivative does not seem to match the function.
 101  *
 102  *  \note The function/derivative vectors and spacing are always double-valued
 103  *        to avoid accuracy loss.
 104  */
 105 void throwUnlessDerivativeIsConsistentWithFunction(ArrayRef<const double>       function,
 106                                                    ArrayRef<const double>       derivative,
 107                                                    double                       inputSpacing,
 108                                                    const std::pair<real, real>& range);
 109
 110
 111 /*! \brief Find smallest quotient between analytical function and its 2nd derivative
 112  *
 113  *  Used to calculate spacing for quadratic spline tables. This function divides the
 114  *  function value by the second derivative (or a very small number when that is zero),
 115  *  and returns the smallest such quotient found in the range.
 116  *
 117  *  Our quadratic tables corresponds to linear interpolation of the derivative,
 118  *  which means the derivative will typically have larger error than the value
 119  *  when interpolating. The spacing required to reach a particular relative
 120  *  tolerance in the derivative depends on the quotient between the first
 121  *  derivative and the third derivative of the function itself.
 122  *
 123  *  You should call this routine with the analytical derivative as the "function"
 124  *  parameter, and the quotient between "function and second derivative" will
 125  *  then correspond to the quotient bewteen the derivative and the third derivative
 126  *  of the actual function we want to tabulate.
 127  *
 128  *  Since all functions that can be tabulated efficiently are reasonably smooth,
 129  *  we simply check 1,000 points in the interval rather than bother about
 130  *  implementing any complicated global optimization scheme.
 131  *
 132  *  \param f          Analytical function
 133  *  \param range      Interval
 134  *
 135  *  \return Smallest quotient found in range.
 136  *
 137  *  \note The function is always double-valued to avoid accuracy loss.
 138  */
 139 real findSmallestQuotientOfFunctionAndSecondDerivative(const std::function<double(double)>& f,
 140                                                        const std::pair<real, real>&         range);
 141
 142
 143 /*! \brief Find smallest quotient between vector of values and its 2nd derivative
 144  *
 145  *  Used to calculate spacing for quadratic spline tables. This function divides the
 146  *  function value by the second derivative (or a very small number when that is zero),
 147  *  and returns the smallest such quotient found in the range.
 148  *
 149  *  Our quadratic tables corresponds to linear interpolation of the derivative,
 150  *  which means the derivative will typically have larger error than the value
 151  *  when interpolating. The spacing required to reach a particular relative
 152  *  tolerance in the derivative depends on the quotient between the first
 153  *  derivative and the third derivative of the function itself.
 154  *
 155  *  You should call this routine with the analytical derivative as the "function"
 156  *  parameter, and the quotient between "function and second derivative" will
 157  *  then correspond to the quotient bewteen the derivative and the third derivative
 158  *  of the actual function we want to tabulate.
 159  *
 160  *  \param function     Vector with function values
 161  *  \param inputSpacing Spacing between function values
 162  *  \param range        Interval to check
 163  *
 164  *  \return Smallest quotient found in range.
 165  *
 166  *  \note The function vector and input spacing are always double-valued to
 167  *        avoid accuracy loss.
 168  */
 169 real findSmallestQuotientOfFunctionAndSecondDerivative(ArrayRef<const double>       function,
 170                                                        double                       inputSpacing,
 171                                                        const std::pair<real, real>& range);
 172
 173
 174 /*! \brief Find smallest quotient between analytical function and its 3rd derivative
 175  *
 176  *  Used to calculate table spacing. This function divides the function value
 177  *  by the second derivative (or a very small number when that is zero), and
 178  *  returns the smallest such quotient found in the range.
 179  *
 180  *  Our quadratic tables corresponds to linear interpolation of the derivative,
 181  *  which means the derivative will typically have larger error than the value
 182  *  when interpolating. The spacing required to reach a particular relative
 183  *  tolerance in the derivative depends on the quotient between the first
 184  *  derivative and the third derivative of the function itself.
 185  *
 186  *  You should call this routine with the analytical derivative as the "function"
 187  *  parameter, and the quotient between "function and second derivative" will
 188  *  then correspond to the quotient bewteen the derivative and the third derivative
 189  *  of the actual function we want to tabulate.
 190  *
 191  *  Since all functions that can be tabulated efficiently are reasonably smooth,
 192  *  we simply check 1,000 points in the interval rather than bother about
 193  *  implementing any complicated global optimization scheme.
 194  *
 195  *  \param f          Analytical function
 196  *  \param range      Interval
 197  *
 198  *  \return Smallest quotient found in range.
 199  *
 200  *  \note The function is always double-valued to avoid accuracy loss.
 201  */
 202 real findSmallestQuotientOfFunctionAndThirdDerivative(const std::function<double(double)>& f,
 203                                                       const std::pair<real, real>&         range);
 204
 205
 206 /*! \brief Find smallest quotient between function and 2nd derivative (vectors)
 207  *
 208  *  Used to calculate table spacing. This function divides the function value
 209  *  by the second derivative (or a very small number when that is zero), and
 210  *  returns the smallest such quotient found in the range.
 211  *
 212  *  Our quadratic tables corresponds to linear interpolation of the derivative,
 213  *  which means the derivative will typically have larger error than the value
 214  *  when interpolating. The spacing required to reach a particular relative
 215  *  tolerance in the derivative depends on the quotient between the first
 216  *  derivative and the third derivative of the function itself.
 217  *
 218  *  You should call this routine with the analytical derivative as the "function"
 219  *  parameter, and the quotient between "function and second derivative" will
 220  *  then correspond to the quotient bewteen the derivative and the third derivative
 221  *  of the actual function we want to tabulate.
 222  *
 223  *  \param function     Vector with function values
 224  *  \param inputSpacing Spacing between function values
 225  *  \param range        Interval to check
 226  *
 227  *  \return Smallest quotient found in range.
 228  *
 229  *  \note The function vector and input spacing are always double-valued to
 230  *        avoid accuracy loss.
 231  */
 232 real findSmallestQuotientOfFunctionAndThirdDerivative(ArrayRef<const double>       function,
 233                                                       double                       inputSpacing,
 234                                                       const std::pair<real, real>& range);
 235
 236
 237 /*! \brief Calculate second derivative of vector and return vector of same length
 238  *
 239  *  5-point approximations are used, with endpoints using non-center interpolation.
 240  *
 241  *  \param f       Vector (function) for which to calculate second derivative
 242  *  \param spacing Spacing of input data.
 243  *
 244  *  \throws If the input vector has fewer than five data points.
 245  *
 246  * \note This function always uses double precision arguments since it is meant
 247  *       to be used on raw user input data for tables, where we want to avoid
 248  *       accuracy loss (since differentiation can be numerically fragile).
 249  */
 250 std::vector<double> vectorSecondDerivative(ArrayRef<const double> f, double spacing);
 251
 252
 253 /*! \brief Copy (temporary) table data into aligned multiplexed vector
 254  *
 255  *  This routine takes the temporary data generated for a single table
 256  *  and writes multiplexed output into a multiple-table-data vector.
 257  *  If the output vector is empty we will resize it to fit the data, and
 258  *  otherwise we assert the size is correct to add out input data.
 259  *
 260  *  \tparam T     Type of container for input data
 261  *  \tparam U     Type of container for output data
 262  *
 263  *  \param[in]    inputData               Input data for single table
 264  *  \param[inout] multiplexedOutputData   Multiplexed output vector, many tables.
 265  *  \param[in]    valuesPerTablePoint     Number of real values for each table point,
 266  *                                        for instance 4 in DDFZ tables.
 267  *  \param[in]    numTables               Number of tables mixed into multiplexed output
 268  *  \param[in]    thisTableIndex          Index of this table in multiplexed output
 269  *
 270  *  \note The output container type can be different from the input since the latter
 271  *        sometimes uses an aligned allocator so the data can be loaded efficiently
 272  *        in the GROMACS nonbonded kernels.
 273  */
 274 template<class T, class U>
 275 void fillMultiplexedTableData(const T     inputData,
 276                               U*          multiplexedOutputData,
 277                               std::size_t valuesPerTablePoint,
 278                               std::size_t numTables,
 279                               std::size_t thisTableIndex)
 280 {
 281     if (multiplexedOutputData->empty())
 282     {
 283         multiplexedOutputData->resize(inputData.size() * numTables);
 284     }
 285     else
 286     {
 287         GMX_ASSERT(inputData.size() * numTables == multiplexedOutputData->size(),
 288                    "Size mismatch when multiplexing table data");
 289     }
 290
 291     GMX_ASSERT(inputData.size() % valuesPerTablePoint == 0,
 292                "Input table size must be a multiple of valuesPerTablePoint");
 293
 294     std::size_t points = inputData.size() / valuesPerTablePoint;
 295
 296     for (std::size_t i = 0; i < points; i++)
 297     {
 298         std::size_t inputOffset  = valuesPerTablePoint * i;
 299         std::size_t outputOffset = valuesPerTablePoint * (numTables * i + thisTableIndex);
 300
 301         for (std::size_t j = 0; j < valuesPerTablePoint; j++)
 302         {
 303             (*multiplexedOutputData)[outputOffset + j] = inputData[inputOffset + j];
 304         }
 305     }
 306 }
 307
 308
 309 } // namespace internal
 310
 311 } // namespace gmx
 312
 313 #endif // GMX_TABLES_SPLINEUTIL_H