src/gromacs/tables/splineutil.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2016,2017,2018,2019,2020, 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 throwUnlessDerivativeIsConsistentWithFunction(const std::function<double(double)>& function,
  67                                                    const std::function<double(double)>& derivative,
  68                                                    const std::pair<real, real>&         range)
  69 {
  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;
  80
  81     // NOLINTNEXTLINE(clang-analyzer-security.FloatLoopCounter)
  82     for (double x = newRange.first; x <= newRange.second; x += dx)
  83     {
  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);
  87
  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
  91         double expectedErr =
  92                 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(
 107                 "Derivative inconsistent with analytical function in range [%f,%f]", minFail, maxFail)));
 108     }
 109 }
 110
 111
 112 void 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   = 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;
 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])
 129                                  / (inputSpacing * inputSpacing);
 130
 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;
 136
 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)
 140         {
 141             isConsistent = false;
 142             minFail      = std::min(minFail, i);
 143             maxFail      = std::max(maxFail, i);
 144         }
 145     }
 146     if (!isConsistent)
 147     {
 148         GMX_THROW(InconsistentInputError(formatString(
 149                 "Derivative inconsistent with numerical vector for elements %zu-%zu", minFail + 1, maxFail + 1)));
 150     }
 151 }
 152
 153
 154 /*! \brief Calculate absolute quotient of function and its second derivative
 155  *
 156  * This is a utility function used in the functions to find the smallest quotient
 157  * in a range.
 158  *
 159  * \param[in]    previousPoint Value of function at x-h.
 160  * \param[in]    thisPoint     Value of function at x.
 161  * \param[in]    nextPoint     Value of function at x+h.
 162  * \param[in]    spacing       Value of h.
 163  *
 164  * \return The absolute value of the quotient. If either the function or second
 165  *         derivative is smaller than sqrt(GMX_REAL_MIN), they will be set to
 166  *         that value.
 167  */
 168 static double quotientOfFunctionAndSecondDerivative(double previousPoint,
 169                                                     double thisPoint,
 170                                                     double nextPoint,
 171                                                     double spacing)
 172 {
 173     double lowerLimit = static_cast<double>(std::sqrt(GMX_REAL_MIN));
 174     double value      = std::max(std::abs(thisPoint), lowerLimit);
 175     double secondDerivative =
 176             std::abs((previousPoint - 2.0 * thisPoint + nextPoint) / (spacing * spacing));
 177
 178     // Make sure we do not divide by zero. This limit is arbitrary,
 179     // but it doesnt matter since this point will have a very large value,
 180     // and the whole routine is searching for the smallest value.
 181     secondDerivative = std::max(secondDerivative, lowerLimit);
 182
 183     return (value / secondDerivative);
 184 }
 185
 186
 187 real findSmallestQuotientOfFunctionAndSecondDerivative(const std::function<double(double)>& f,
 188                                                        const std::pair<real, real>&         range)
 189 {
 190     // Since the numerical second derivative will evaluate extra points
 191     // we shrink the interval slightly to avoid calling the function with values
 192     // outside the range specified.
 193     double                    h = std::pow(GMX_DOUBLE_EPS, 0.25);
 194     std::pair<double, double> newRange(range.first + h, range.second - h);
 195     const int                 points      = 500; // arbitrary
 196     double                    dx          = (newRange.second - newRange.first) / points;
 197     double                    minQuotient = GMX_REAL_MAX;
 198
 199     // NOLINTNEXTLINE(clang-analyzer-security.FloatLoopCounter)
 200     for (double x = newRange.first; x <= newRange.second; x += dx)
 201     {
 202         minQuotient = std::min(minQuotient,
 203                                quotientOfFunctionAndSecondDerivative(f(x - h), f(x), f(x + h), h));
 204     }
 205
 206     return static_cast<real>(minQuotient);
 207 }
 208
 209
 210 real findSmallestQuotientOfFunctionAndSecondDerivative(ArrayRef<const double>       function,
 211                                                        double                       inputSpacing,
 212                                                        const std::pair<real, real>& range)
 213 {
 214
 215     std::size_t firstIndex  = static_cast<std::size_t>(range.first / inputSpacing);
 216     std::size_t lastIndex   = static_cast<std::size_t>(range.second / inputSpacing);
 217     double      minQuotient = GMX_REAL_MAX;
 218
 219     for (std::size_t i = firstIndex + 1; (i + 1) < lastIndex; i++)
 220     {
 221         minQuotient = std::min(minQuotient,
 222                                quotientOfFunctionAndSecondDerivative(
 223                                        function[i - 1], function[i], function[i + 1], inputSpacing));
 224     }
 225     return static_cast<real>(minQuotient);
 226 }
 227
 228
 229 /*! \brief Calculate absolute quotient of function and its third derivative
 230  *
 231  * This is a utility function used in the functions to find the smallest quotient
 232  * in a range.
 233  *
 234  * \param[in]    previousPreviousPoint Value of function at x-2h.
 235  * \param[in]    previousPoint         Value of function at x-h.
 236  * \param[in]    thisPoint             Value of function at x.
 237  * \param[in]    nextPoint             Value of function at x+h.
 238  * \param[in]    nextNextPoint         Value of function at x+2h.
 239  * \param[in]    spacing               Value of h.
 240  *
 241  * \return The absolute value of the quotient. If either the function or third
 242  *         derivative is smaller than sqrt(GMX_REAL_MIN), they will be set to
 243  *         that value.
 244  */
 245 static double quotientOfFunctionAndThirdDerivative(double previousPreviousPoint,
 246                                                    double previousPoint,
 247                                                    double thisPoint,
 248                                                    double nextPoint,
 249                                                    double nextNextPoint,
 250                                                    double spacing)
 251 {
 252     double lowerLimit = static_cast<double>(std::sqrt(GMX_REAL_MIN));
 253     double value      = std::max(std::abs(thisPoint), lowerLimit);
 254     double thirdDerivative =
 255             std::abs((nextNextPoint - 2 * nextPoint + 2 * previousPoint - previousPreviousPoint)
 256                      / (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 findSmallestQuotientOfFunctionAndThirdDerivative(const std::function<double(double)>& f,
 268                                                       const std::pair<real, real>&         range)
 269 {
 270     // Since the numerical second derivative will evaluate extra points
 271     // we shrink the interval slightly to avoid calling the function with values
 272     // outside the range specified.
 273     double h = std::pow(GMX_DOUBLE_EPS, 0.2); // optimal spacing for 3rd derivative
 274     std::pair<double, double> newRange(range.first + 2 * h, range.second - 2 * h);
 275     const int                 points      = 500; // arbitrary
 276     double                    dx          = (newRange.second - newRange.first) / points;
 277     double                    minQuotient = GMX_REAL_MAX;
 278
 279     // NOLINTNEXTLINE(clang-analyzer-security.FloatLoopCounter)
 280     for (double x = newRange.first; x <= newRange.second; x += dx)
 281     {
 282         minQuotient = std::min(minQuotient,
 283                                quotientOfFunctionAndThirdDerivative(
 284                                        f(x - 2 * h), f(x - h), f(x), f(x + h), f(x + 2 * h), h));
 285     }
 286     return static_cast<real>(minQuotient);
 287 }
 288
 289
 290 real findSmallestQuotientOfFunctionAndThirdDerivative(ArrayRef<const double>       function,
 291                                                       double                       inputSpacing,
 292                                                       const std::pair<real, real>& range)
 293 {
 294
 295     std::size_t firstIndex  = static_cast<std::size_t>(range.first / inputSpacing);
 296     std::size_t lastIndex   = static_cast<std::size_t>(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(
 302                 minQuotient,
 303                 quotientOfFunctionAndThirdDerivative(
 304                         function[i - 2], function[i - 1], function[i], function[i + 1], function[i + 2], inputSpacing));
 305     }
 306     return static_cast<real>(minQuotient);
 307 }
 308
 309
 310 std::vector<double> vectorSecondDerivative(ArrayRef<const double> f, double spacing)
 311 {
 312     if (f.size() < 5)
 313     {
 314         GMX_THROW(APIError("Too few points in vector for 5-point differentiation"));
 315     }
 316
 317     std::vector<double> d(f.size());
 318     std::size_t         i;
 319
 320     // 5-point formula evaluated for points 0,1
 321     i    = 0;
 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);
 324     i = 1;
 325     d[i] = (-f[i + 3] + 4 * f[i + 2] + 6 * f[i + 1] - 20 * f[i] + 11 * f[i - 1]) / (12 * spacing * spacing);
 326
 327     for (std::size_t i = 2; i < d.size() - 2; i++)
 328     {
 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);
 331     }
 332
 333     // 5-point formula evaluated for points 3,4
 334     i = d.size() - 2;
 335     d[i] = (11 * f[i + 1] - 20 * f[i] + 6 * f[i - 1] + 4 * f[i - 2] - f[i - 3]) / (12 * spacing * spacing);
 336     i    = d.size() - 1;
 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);
 339
 340     return d;
 341 }
 342
 343
 344 } // namespace internal
 345
 346 } // namespace gmx