src/gromacs/tables/splineutil.cpp

   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  * Implements internal utility functions for spline tables
  39  *
  40  * \author Erik Lindahl <erik.lindahl@gmail.com>
  41  * \ingroup module_tables
  42  */
  43 #include "gmxpre.h"
  44
  45 #include "splineutil.h"
  46
  47 #include <cmath>
  48
  49 #include <algorithm>
  50 #include <functional>
  51 #include <utility>
  52 #include <vector>
  53
  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"
  59
  60 namespace gmx
  61 {
  62
  63 namespace internal
  64 {
  65
  66 void
  67 throwUnlessDerivativeIsConsistentWithFunction(const std::function<double(double)>  &function,
  68                                               const std::function<double(double)>  &derivative,
  69                                               const std::pair<real, real>          &range)
  70 {
  71     // Since the numerical derivative will evaluate extra points
  72     // we shrink the interval slightly to avoid calling the function with values
  73     // outside the range specified.
  74     double                     h            = std::cbrt(GMX_DOUBLE_EPS); // ideal spacing
  75     std::pair<double, double>  newRange(range.first + h, range.second - h);
  76     const int                  points       = 1000;                      // arbitrary
  77     double                     dx           = (newRange.second - newRange.first) / points;
  78     bool                       isConsistent = true;
  79     double                     minFail      = newRange.second;
  80     double                     maxFail      = newRange.first;
  81
  82     // NOLINTNEXTLINE(clang-analyzer-security.FloatLoopCounter)
  83     for (double x = newRange.first; x <= newRange.second; x += dx)
  84     {
  85         double analyticalDerivative = derivative(x);
  86         double numericalDerivative  = (function(x+h)-function(x-h))/(2*h);
  87         double thirdDerivative      = (derivative(x+h)-2.0*derivative(x)+derivative(x-h))/(h*h);
  88
  89         // We make two types of errors in numerical calculation of the derivative:
  90         // - The truncation error: eps * |f| / h
  91         // - The rounding error: h * h * |f'''| / 6.0
  92         double expectedErr = GMX_DOUBLE_EPS*std::abs(function(x))/h + h*h*std::abs(thirdDerivative)/6.0;
  93
  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)
  97         {
  98             isConsistent = false;
  99             minFail      = std::min(minFail, x);
 100             maxFail      = std::max(maxFail, x);
 101         }
 102     }
 103
 104     if (!isConsistent)
 105     {
 106         GMX_THROW(InconsistentInputError(formatString("Derivative inconsistent with analytical function in range [%d,%d]", minFail, maxFail)));
 107     }
 108 }
 109
 110
 111 void
 112 throwUnlessDerivativeIsConsistentWithFunction(ArrayRef<const double>        function,
 113                                               ArrayRef<const double>        derivative,
 114                                               double                        inputSpacing,
 115                                               const std::pair<real, real>  &range)
 116 {
 117     std::size_t     firstIndex   = range.first / inputSpacing;
 118     std::size_t     lastIndex    = range.second / inputSpacing;
 119     bool            isConsistent = true;
 120     std::size_t     minFail      = lastIndex;
 121     std::size_t     maxFail      = firstIndex;
 122
 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++)
 125     {
 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])/(inputSpacing * inputSpacing);
 129
 130         // We make two types of errors in numerical calculation of the derivative:
 131         // - The truncation error: eps * |f| / h
 132         // - The rounding error: h * h * |f'''| / 6.0
 133         double expectedErr = GMX_DOUBLE_EPS*std::abs(function[i])/inputSpacing + inputSpacing*inputSpacing*std::abs(thirdDerivative)/6.0;
 134
 135         // To avoid triggering false errors because of compiler optimization or numerical issues
 136         // in the function evalulation we allow an extra factor of 10 in the expected error
 137         if (std::abs(inputDerivative-numericalDerivative) > 10.0*expectedErr)
 138         {
 139             isConsistent = false;
 140             minFail      = std::min(minFail, i);
 141             maxFail      = std::max(maxFail, i);
 142         }
 143     }
 144     if (!isConsistent)
 145     {
 146         GMX_THROW(InconsistentInputError(formatString("Derivative inconsistent with numerical vector for elements %d-%d", minFail+1, maxFail+1)));
 147     }
 148 }
 149
 150
 151 /*! \brief Calculate absolute quotient of function and its second derivative
 152  *
 153  * This is a utility function used in the functions to find the smallest quotient
 154  * in a range.
 155  *
 156  * \param[in]    previousPoint Value of function at x-h.
 157  * \param[in]    thisPoint     Value of function at x.
 158  * \param[in]    nextPoint     Value of function at x+h.
 159  * \param[in]    spacing       Value of h.
 160  *
 161  * \return The absolute value of the quotient. If either the function or second
 162  *         derivative is smaller than sqrt(GMX_REAL_MIN), they will be set to
 163  *         that value.
 164  */
 165 static double
 166 quotientOfFunctionAndSecondDerivative(double     previousPoint,
 167                                       double     thisPoint,
 168                                       double     nextPoint,
 169                                       double     spacing)
 170 {
 171     double lowerLimit        = static_cast<double>(std::sqrt(GMX_REAL_MIN));
 172     double value             = std::max(std::abs( thisPoint ), lowerLimit );
 173     double secondDerivative  = std::abs( (previousPoint - 2.0 * thisPoint + nextPoint) / (spacing * spacing ) );
 174
 175     // Make sure we do not divide by zero. This limit is arbitrary,
 176     // but it doesnt matter since this point will have a very large value,
 177     // and the whole routine is searching for the smallest value.
 178     secondDerivative = std::max(secondDerivative, lowerLimit);
 179
 180     return (value / secondDerivative);
 181 }
 182
 183
 184 real
 185 findSmallestQuotientOfFunctionAndSecondDerivative(const std::function<double(double)>   &f,
 186                                                   const std::pair<real, real>           &range)
 187 {
 188     // Since the numerical second derivative will evaluate extra points
 189     // we shrink the interval slightly to avoid calling the function with values
 190     // outside the range specified.
 191     double                     h           = std::pow( GMX_DOUBLE_EPS, 0.25 );
 192     std::pair<double, double>  newRange(range.first + h, range.second - h);
 193     const int                  points      = 500; // arbitrary
 194     double                     dx          = (newRange.second - newRange.first) / points;
 195     double                     minQuotient = GMX_REAL_MAX;
 196
 197     // NOLINTNEXTLINE(clang-analyzer-security.FloatLoopCounter)
 198     for (double x = newRange.first; x <= newRange.second; x += dx)
 199     {
 200         minQuotient = std::min(minQuotient, quotientOfFunctionAndSecondDerivative(f(x-h), f(x), f(x+h), h));
 201     }
 202
 203     return static_cast<real>(minQuotient);
 204 }
 205
 206
 207
 208
 209
 210
 211 real
 212 findSmallestQuotientOfFunctionAndSecondDerivative(ArrayRef<const double>        function,
 213                                                   double                        inputSpacing,
 214                                                   const std::pair<real, real>   &range)
 215 {
 216
 217     std::size_t  firstIndex  = range.first  / inputSpacing;
 218     std::size_t  lastIndex   = range.second / inputSpacing;
 219     double       minQuotient = GMX_REAL_MAX;
 220
 221     for (std::size_t i = firstIndex + 1; (i + 1) < lastIndex; i++)
 222     {
 223         minQuotient = std::min(minQuotient, quotientOfFunctionAndSecondDerivative(function[i-1], function[i], function[i+1], inputSpacing));
 224     }
 225     return static_cast<real>(minQuotient);
 226 }
 227
 228
 229
 230 /*! \brief Calculate absolute quotient of function and its third derivative
 231  *
 232  * This is a utility function used in the functions to find the smallest quotient
 233  * in a range.
 234  *
 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.
 241  *
 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
 244  *         that value.
 245  */
 246 static double
 247 quotientOfFunctionAndThirdDerivative(double     previousPreviousPoint,
 248                                      double     previousPoint,
 249                                      double     thisPoint,
 250                                      double     nextPoint,
 251                                      double     nextNextPoint,
 252                                      double     spacing)
 253 {
 254     double lowerLimit       = static_cast<double>(std::sqrt(GMX_REAL_MIN));
 255     double value            = std::max(std::abs( thisPoint ), lowerLimit );
 256     double thirdDerivative  = std::abs((nextNextPoint - 2 * nextPoint + 2 * previousPoint - previousPreviousPoint) / (2 * spacing * spacing * spacing));
 257
 258     // Make sure we do not divide by zero. This limit is arbitrary,
 259     // but it doesnt matter since this point will have a very large value,
 260     // and the whole routine is searching for the smallest value.
 261     thirdDerivative = std::max(thirdDerivative, lowerLimit);
 262
 263     return (value / thirdDerivative);
 264 }
 265
 266
 267 real
 268 findSmallestQuotientOfFunctionAndThirdDerivative(const std::function<double(double)>  &f,
 269                                                  const std::pair<real, real>           &range)
 270 {
 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;
 279
 280     // NOLINTNEXTLINE(clang-analyzer-security.FloatLoopCounter)
 281     for (double x = newRange.first; x <= newRange.second; x += dx)
 282     {
 283         minQuotient = std::min(minQuotient, quotientOfFunctionAndThirdDerivative(f(x-2*h), f(x-h), f(x), f(x+h), f(x+2*h), h));
 284     }
 285     return static_cast<real>(minQuotient);
 286 }
 287
 288
 289 real
 290 findSmallestQuotientOfFunctionAndThirdDerivative(ArrayRef<const double>        function,
 291                                                  double                        inputSpacing,
 292                                                  const std::pair<real, real>   &range)
 293 {
 294
 295     std::size_t  firstIndex  = range.first  / inputSpacing;
 296     std::size_t  lastIndex   = range.second / inputSpacing;
 297     double       minQuotient = GMX_REAL_MAX;
 298
 299     for (std::size_t i = firstIndex + 2; (i + 2) < lastIndex; i++)
 300     {
 301         minQuotient = std::min(minQuotient, quotientOfFunctionAndThirdDerivative(function[i-2], function[i-1], function[i], function[i+1], function[i+2], inputSpacing));
 302     }
 303     return static_cast<real>(minQuotient);
 304 }
 305
 306
 307
 308 std::vector<double>
 309 vectorSecondDerivative(ArrayRef<const double> f, double spacing)
 310 {
 311     if (f.size() < 5)
 312     {
 313         GMX_THROW(APIError("Too few points in vector for 5-point differentiation"));
 314     }
 315
 316     std::vector<double> d(f.size());
 317     std::size_t         i;
 318
 319     // 5-point formula evaluated for points 0,1
 320     i    = 0;
 321     d[i] = (11 * f[i+4] - 56 * f[i+3] + 114 * f[i+2] - 104 * f[i+1] + 35 * f[i]) / ( 12 * spacing * spacing);
 322     i    = 1;
 323     d[i] = (-f[i+3] + 4 * f[i+2] + 6 * f[i+1] - 20 * f[i] + 11 * f[i-1]) / ( 12 * spacing * spacing);
 324
 325     for (std::size_t i = 2; i < d.size() - 2; i++)
 326     {
 327         // 5-point formula evaluated for central point (2)
 328         d[i] = (-f[i+2] + 16 * f[i+1] - 30 * f[i] + 16 * f[i-1] - f[i-2]) / (12 * spacing * spacing);
 329     }
 330
 331     // 5-point formula evaluated for points 3,4
 332     i    = d.size() - 2;
 333     d[i] = (11 * f[i+1] - 20 * f[i] + 6 * f[i-1] + 4 * f[i-2] - f[i-3]) / ( 12 * spacing * spacing);
 334     i    = d.size() - 1;
 335     d[i] = (35 * f[i] - 104 * f[i-1] + 114 * f[i-2] - 56 * f[i-3] + 11 * f[i-4]) / ( 12 * spacing * spacing);
 336
 337     return d;
 338 }
 339
 340
 341
 342 }  // namespace internal
 343
 344 }  // namespace gmx