src/gromacs/tables/splineutil.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
   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
  80 throwUnlessDerivativeIsConsistentWithFunction(const std::function<double(double)>  &function,
  81                                               const std::function<double(double)>  &derivative,
  82                                               const std::pair<real, real>          &range);
  83
  84 /*! \brief Ensure vector of derivative values is the derivative of function vector.
  85  *
  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.
  89  *
  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.
  92  *
  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.
  95  *
  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
 100  *
 101  *  \throws If the provided derivative does not seem to match the function.
 102  *
 103  *  \note The function/derivative vectors and spacing are always double-valued
 104  *        to avoid accuracy loss.
 105  */
 106 void
 107 throwUnlessDerivativeIsConsistentWithFunction(ArrayRef<const double>        function,
 108                                               ArrayRef<const double>        derivative,
 109                                               double                        inputSpacing,
 110                                               const std::pair<real, real>  &range);
 111
 112
 113
 114
 115 /*! \brief Find smallest quotient between analytical function and its 2nd derivative
 116  *
 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.
 120  *
 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.
 126  *
 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.
 131  *
 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.
 135  *
 136  *  \param f          Analytical function
 137  *  \param range      Interval
 138  *
 139  *  \return Smallest quotient found in range.
 140  *
 141  *  \note The function is always double-valued to avoid accuracy loss.
 142  */
 143 real
 144 findSmallestQuotientOfFunctionAndSecondDerivative(const std::function<double(double)>  &f,
 145                                                   const std::pair<real, real>          &range);
 146
 147
 148
 149
 150
 151 /*! \brief Find smallest quotient between vector of values and its 2nd derivative
 152  *
 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.
 156  *
 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.
 162  *
 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.
 167  *
 168  *  \param function     Vector with function values
 169  *  \param inputSpacing Spacing between function values
 170  *  \param range        Interval to check
 171  *
 172  *  \return Smallest quotient found in range.
 173  *
 174  *  \note The function vector and input spacing are always double-valued to
 175  *        avoid accuracy loss.
 176  */
 177 real
 178 findSmallestQuotientOfFunctionAndSecondDerivative(ArrayRef<const double>        function,
 179                                                   double                        inputSpacing,
 180                                                   const std::pair<real, real>  &range);
 181
 182
 183
 184
 185
 186 /*! \brief Find smallest quotient between analytical function and its 3rd derivative
 187  *
 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.
 191  *
 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.
 197  *
 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.
 202  *
 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.
 206  *
 207  *  \param f          Analytical function
 208  *  \param range      Interval
 209  *
 210  *  \return Smallest quotient found in range.
 211  *
 212  *  \note The function is always double-valued to avoid accuracy loss.
 213  */
 214 real
 215 findSmallestQuotientOfFunctionAndThirdDerivative(const std::function<double(double)>  &f,
 216                                                  const std::pair<real, real>          &range);
 217
 218
 219
 220
 221 /*! \brief Find smallest quotient between function and 2nd derivative (vectors)
 222  *
 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.
 226  *
 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.
 232  *
 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.
 237  *
 238  *  \param function     Vector with function values
 239  *  \param inputSpacing Spacing between function values
 240  *  \param range        Interval to check
 241  *
 242  *  \return Smallest quotient found in range.
 243  *
 244  *  \note The function vector and input spacing are always double-valued to
 245  *        avoid accuracy loss.
 246  */
 247 real
 248 findSmallestQuotientOfFunctionAndThirdDerivative(ArrayRef<const double>        function,
 249                                                  double                        inputSpacing,
 250                                                  const std::pair<real, real>  &range);
 251
 252
 253 /*! \brief Calculate second derivative of vector and return vector of same length
 254  *
 255  *  5-point approximations are used, with endpoints using non-center interpolation.
 256  *
 257  *  \param f       Vector (function) for which to calculate second derivative
 258  *  \param spacing Spacing of input data.
 259  *
 260  *  \throws If the input vector has fewer than five data points.
 261  *
 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).
 265  */
 266 std::vector<double>
 267 vectorSecondDerivative(ArrayRef<const double>       f,
 268                        double                       spacing);
 269
 270
 271 /*! \brief Copy (temporary) table data into aligned multiplexed vector
 272  *
 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.
 277  *
 278  *  \tparam T     Type of container for input data
 279  *  \tparam U     Type of container for output data
 280  *
 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
 287  *
 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.
 291  */
 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)
 299 {
 300     if (multiplexedOutputData->empty())
 301     {
 302         multiplexedOutputData->resize( inputData.size() * numTables );
 303     }
 304     else
 305     {
 306         GMX_ASSERT(inputData.size() * numTables == multiplexedOutputData->size(),
 307                    "Size mismatch when multiplexing table data");
 308     }
 309
 310     GMX_ASSERT(inputData.size() % valuesPerTablePoint == 0,
 311                "Input table size must be a multiple of valuesPerTablePoint");
 312
 313     std::size_t points = inputData.size() / valuesPerTablePoint;
 314
 315     for (std::size_t i = 0; i < points; i++)
 316     {
 317         std::size_t inputOffset  = valuesPerTablePoint * i;
 318         std::size_t outputOffset = valuesPerTablePoint * ( numTables * i + thisTableIndex );
 319
 320         for (std::size_t j = 0; j < valuesPerTablePoint; j++)
 321         {
 322             (*multiplexedOutputData)[outputOffset + j] = inputData[inputOffset + j];
 323         }
 324     }
 325 }
 326
 327
 328 }      // namespace internal
 329
 330 }      // namespace gmx
 331
 332 #endif // GMX_TABLES_SPLINEUTIL_H