/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2016,2017,2018, by the GROMACS development team, led by
+ * Copyright (c) 2016,2017,2018,2019, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
* \param[out] functionTableData Output table with function data
* \param[out] derivativeTableData OUtput table with (adjusted) derivative data
*/
-void
-fillSingleQuadraticSplineTableData(const std::function<double(double)> &function,
- const std::function<double(double)> &derivative,
- const std::pair<real, real> &range,
- double spacing,
- std::vector<real> *functionTableData,
- std::vector<real> *derivativeTableData)
+void fillSingleQuadraticSplineTableData(const std::function<double(double)>& function,
+ const std::function<double(double)>& derivative,
+ const std::pair<real, real>& range,
+ double spacing,
+ std::vector<real>* functionTableData,
+ std::vector<real>* derivativeTableData)
{
- std::size_t endIndex = static_cast<std::size_t>(range.second / spacing + 2);
+ std::size_t endIndex = static_cast<std::size_t>(range.second / spacing + 2);
functionTableData->resize(endIndex);
derivativeTableData->resize(endIndex);
- double maxMagnitude = 0.0001*GMX_REAL_MAX;
- bool functionIsInRange = true;
- std::size_t lastIndexInRange = endIndex - 1;
+ double maxMagnitude = 0.0001 * GMX_REAL_MAX;
+ bool functionIsInRange = true;
+ std::size_t lastIndexInRange = endIndex - 1;
for (int i = endIndex - 1; i >= 0; i--)
{
- double x = i * spacing;
+ double x = i * spacing;
double tmpFunctionValue;
double tmpDerivativeValue;
// Calculate third derivative term (2nd derivative of the derivative)
// Make sure we stay in range. In practice this means we use one-sided
// interpolation at the interval endpoints (indentical to an offset for 3-point formula)
- const double h = std::pow( GMX_DOUBLE_EPS, 0.25 );
- double y = std::min( std::max(x, range.first + h), range.second - h);
- double thirdDerivativeValue = ( derivative(y+h) - 2.0 * derivative(y) + derivative(y-h) ) / ( h * h );
+ const double h = std::pow(GMX_DOUBLE_EPS, 0.25);
+ double y = std::min(std::max(x, range.first + h), range.second - h);
+ double thirdDerivativeValue =
+ (derivative(y + h) - 2.0 * derivative(y) + derivative(y - h)) / (h * h);
- tmpDerivativeValue = derivative(x) - spacing * spacing * thirdDerivativeValue / 12.0;
+ tmpDerivativeValue = derivative(x) - spacing * spacing * thirdDerivativeValue / 12.0;
if (std::abs(tmpFunctionValue) > maxMagnitude || std::abs(tmpDerivativeValue) > maxMagnitude)
{
// we simply make a linear function from the last in-range value of the derivative.
double lastIndexFunction = (*functionTableData)[lastIndexInRange];
double lastIndexDerivative = (*derivativeTableData)[lastIndexInRange];
- (*functionTableData)[i] = lastIndexFunction + lastIndexDerivative * (i - lastIndexInRange) * spacing;
- (*derivativeTableData)[i] = lastIndexDerivative;
+ (*functionTableData)[i] =
+ lastIndexFunction + lastIndexDerivative * (i - lastIndexInRange) * spacing;
+ (*derivativeTableData)[i] = lastIndexDerivative;
}
}
}
* \param[out] functionTableData Output table with function data
* \param[out] derivativeTableData OUtput table with (adjusted) derivative data
*/
-void
-fillSingleQuadraticSplineTableData(ArrayRef<const double> function,
- ArrayRef<const double> derivative,
- double inputSpacing,
- const std::pair<real, real> &range,
- double spacing,
- std::vector<real> *functionTableData,
- std::vector<real> *derivativeTableData)
+void fillSingleQuadraticSplineTableData(ArrayRef<const double> function,
+ ArrayRef<const double> derivative,
+ double inputSpacing,
+ const std::pair<real, real>& range,
+ double spacing,
+ std::vector<real>* functionTableData,
+ std::vector<real>* derivativeTableData)
{
- std::size_t endIndex = static_cast<std::size_t>(range.second / spacing + 2);
+ std::size_t endIndex = static_cast<std::size_t>(range.second / spacing + 2);
functionTableData->resize(endIndex);
derivativeTableData->resize(endIndex);
- std::vector<double> thirdDerivative(internal::vectorSecondDerivative(derivative, inputSpacing));
+ std::vector<double> thirdDerivative(internal::vectorSecondDerivative(derivative, inputSpacing));
- double maxMagnitude = 0.0001*GMX_REAL_MAX;
- bool functionIsInRange = true;
- std::size_t lastIndexInRange = endIndex - 1;
+ double maxMagnitude = 0.0001 * GMX_REAL_MAX;
+ bool functionIsInRange = true;
+ std::size_t lastIndexInRange = endIndex - 1;
for (int i = endIndex - 1; i >= 0; i--)
{
- double x = i * spacing;
+ double x = i * spacing;
double tmpFunctionValue;
double tmpDerivativeValue;
double inputEps = inputXTab - inputIndex;
// Linear interpolation of input derivative and third derivative
- double thirdDerivativeValue = (1.0 - inputEps) * thirdDerivative[inputIndex] + inputEps * thirdDerivative[inputIndex+1];
- double derivativeValue = (1.0 - inputEps) * derivative[inputIndex] + inputEps * derivative[inputIndex+1];
+ double thirdDerivativeValue = (1.0 - inputEps) * thirdDerivative[inputIndex]
+ + inputEps * thirdDerivative[inputIndex + 1];
+ double derivativeValue =
+ (1.0 - inputEps) * derivative[inputIndex] + inputEps * derivative[inputIndex + 1];
// Quadratic interpolation for function value
- tmpFunctionValue = function[inputIndex] + 0.5 * (derivative[inputIndex] + derivativeValue) * inputEps * inputSpacing;
- tmpDerivativeValue = derivativeValue - spacing * spacing * thirdDerivativeValue / 12.0;
+ tmpFunctionValue = function[inputIndex]
+ + 0.5 * (derivative[inputIndex] + derivativeValue) * inputEps * inputSpacing;
+ tmpDerivativeValue = derivativeValue - spacing * spacing * thirdDerivativeValue / 12.0;
if (std::abs(tmpFunctionValue) > maxMagnitude || std::abs(tmpDerivativeValue) > maxMagnitude)
{
// we simply make a linear function from the last in-range value of the derivative.
double lastIndexFunction = (*functionTableData)[lastIndexInRange];
double lastIndexDerivative = (*derivativeTableData)[lastIndexInRange];
- (*functionTableData)[i] = lastIndexFunction + lastIndexDerivative * (i - lastIndexInRange) * spacing;
- (*derivativeTableData)[i] = lastIndexDerivative;
+ (*functionTableData)[i] =
+ lastIndexFunction + lastIndexDerivative * (i - lastIndexInRange) * spacing;
+ (*derivativeTableData)[i] = lastIndexDerivative;
}
}
}
*
* \throws If the vector lengths do not match.
*/
-void
-fillDdfzTableData(const std::vector<real> &functionTableData,
- const std::vector<real> &derivativeTableData,
- std::vector<real> *ddfzTableData)
+void fillDdfzTableData(const std::vector<real>& functionTableData,
+ const std::vector<real>& derivativeTableData,
+ std::vector<real>* ddfzTableData)
{
- GMX_ASSERT(functionTableData.size() == derivativeTableData.size(), "Mismatching vector lengths");
+ GMX_ASSERT(functionTableData.size() == derivativeTableData.size(),
+ "Mismatching vector lengths");
std::size_t points = functionTableData.size();
for (std::size_t i = 0; i < points; i++)
{
- (*ddfzTableData)[4*i] = derivativeTableData[i];
+ (*ddfzTableData)[4 * i] = derivativeTableData[i];
- double nextDerivative = ( i < functionTableData.size() - 1 ) ? derivativeTableData[i+1] : 0.0;
+ double nextDerivative = (i < functionTableData.size() - 1) ? derivativeTableData[i + 1] : 0.0;
- (*ddfzTableData)[4*i + 1] = nextDerivative - derivativeTableData[i];
- (*ddfzTableData)[4*i + 2] = functionTableData[i];
- (*ddfzTableData)[4*i + 3] = 0.0;
+ (*ddfzTableData)[4 * i + 1] = nextDerivative - derivativeTableData[i];
+ (*ddfzTableData)[4 * i + 2] = functionTableData[i];
+ (*ddfzTableData)[4 * i + 3] = 0.0;
}
}
-} // namespace
+} // namespace
+const real QuadraticSplineTable::defaultTolerance = 10.0 * GMX_FLOAT_EPS;
-const real
-QuadraticSplineTable::defaultTolerance = 10.0 * GMX_FLOAT_EPS;
-
-QuadraticSplineTable::QuadraticSplineTable(std::initializer_list<AnalyticalSplineTableInput> analyticalInputList,
- const std::pair<real, real> &range,
- real tolerance)
- : numFuncInTable_(analyticalInputList.size()), range_(range)
+QuadraticSplineTable::QuadraticSplineTable(std::initializer_list<AnalyticalSplineTableInput> analyticalInputList,
+ const std::pair<real, real>& range,
+ real tolerance) :
+ numFuncInTable_(analyticalInputList.size()),
+ range_(range)
{
// Sanity check on input values
- if (range_.first < 0.0 || (range_.second-range_.first) < 0.001)
+ if (range_.first < 0.0 || (range_.second - range_.first) < 0.001)
{
- GMX_THROW(InvalidInputError("Range to tabulate cannot include negative values and must span at least 0.001"));
+ GMX_THROW(InvalidInputError(
+ "Range to tabulate cannot include negative values and must span at least 0.001"));
}
if (tolerance < GMX_REAL_EPS)
double minQuotient = GMX_REAL_MAX;
// loop over all functions to find smallest spacing
- for (const auto &thisFuncInput : analyticalInputList)
+ for (const auto& thisFuncInput : analyticalInputList)
{
try
{
- internal::throwUnlessDerivativeIsConsistentWithFunction(thisFuncInput.function, thisFuncInput.derivative, range_);
+ internal::throwUnlessDerivativeIsConsistentWithFunction(
+ thisFuncInput.function, thisFuncInput.derivative, range_);
}
- catch (gmx::GromacsException &ex)
+ catch (gmx::GromacsException& ex)
{
- ex.prependContext("Error generating quadratic spline table for function '" + thisFuncInput.desc + "'");
+ ex.prependContext("Error generating quadratic spline table for function '"
+ + thisFuncInput.desc + "'");
throw;
}
// Calculate the required table spacing h. The error we make with linear interpolation
// This means we can compute the required spacing as h = sqrt(12*tolerance*min(f'/f''')),
// where f'/f''' is the first and third derivative of the function, respectively.
- double thisMinQuotient = internal::findSmallestQuotientOfFunctionAndSecondDerivative(thisFuncInput.derivative, range_);
+ double thisMinQuotient = internal::findSmallestQuotientOfFunctionAndSecondDerivative(
+ thisFuncInput.derivative, range_);
minQuotient = std::min(minQuotient, thisMinQuotient);
}
if (range_.second * tableScale_ > 1e6)
{
- GMX_THROW(ToleranceError("Over a million points would be required for table; decrease range or increase tolerance"));
+ GMX_THROW(
+ ToleranceError("Over a million points would be required for table; decrease range "
+ "or increase tolerance"));
}
// Loop over all tables again.
// combine them into a multiplexed table function.
std::size_t funcIndex = 0;
- for (const auto &thisFuncInput : analyticalInputList)
+ for (const auto& thisFuncInput : analyticalInputList)
{
try
{
std::vector<real> tmpDerTableData;
std::vector<real> tmpDdfzTableData;
- fillSingleQuadraticSplineTableData(thisFuncInput.function,
- thisFuncInput.derivative,
- range_,
- spacing,
- &tmpFuncTableData,
- &tmpDerTableData);
+ fillSingleQuadraticSplineTableData(thisFuncInput.function, thisFuncInput.derivative,
+ range_, spacing, &tmpFuncTableData, &tmpDerTableData);
fillDdfzTableData(tmpFuncTableData, tmpDerTableData, &tmpDdfzTableData);
- internal::fillMultiplexedTableData(tmpDerTableData, &derivativeMultiTableData_,
- 1, numFuncInTable_, funcIndex);
+ internal::fillMultiplexedTableData(tmpDerTableData, &derivativeMultiTableData_, 1,
+ numFuncInTable_, funcIndex);
- internal::fillMultiplexedTableData(tmpDdfzTableData, &ddfzMultiTableData_,
- 4, numFuncInTable_, funcIndex);
+ internal::fillMultiplexedTableData(tmpDdfzTableData, &ddfzMultiTableData_, 4,
+ numFuncInTable_, funcIndex);
funcIndex++;
}
- catch (gmx::GromacsException &ex)
+ catch (gmx::GromacsException& ex)
{
- ex.prependContext("Error generating quadratic spline table for function '" + thisFuncInput.desc + "'");
+ ex.prependContext("Error generating quadratic spline table for function '"
+ + thisFuncInput.desc + "'");
throw;
}
}
}
-QuadraticSplineTable::QuadraticSplineTable(std::initializer_list<NumericalSplineTableInput> numericalInputList,
- const std::pair<real, real> &range,
- real tolerance)
- : numFuncInTable_(numericalInputList.size()), range_(range)
+QuadraticSplineTable::QuadraticSplineTable(std::initializer_list<NumericalSplineTableInput> numericalInputList,
+ const std::pair<real, real>& range,
+ real tolerance) :
+ numFuncInTable_(numericalInputList.size()),
+ range_(range)
{
// Sanity check on input values
- if (range.first < 0.0 || (range.second-range.first) < 0.001)
+ if (range.first < 0.0 || (range.second - range.first) < 0.001)
{
- GMX_THROW(InvalidInputError("Range to tabulate cannot include negative values and must span at least 0.001"));
+ GMX_THROW(InvalidInputError(
+ "Range to tabulate cannot include negative values and must span at least 0.001"));
}
if (tolerance < GMX_REAL_EPS)
// the requested range before starting to calculate derivatives
if (thisFuncInput.function.size() < range_.second / thisFuncInput.spacing + 1)
{
- GMX_THROW(InconsistentInputError("Table input vectors must cover requested range, and a margin beyond the upper endpoint"));
+ GMX_THROW(
+ InconsistentInputError("Table input vectors must cover requested range, "
+ "and a margin beyond the upper endpoint"));
}
if (thisFuncInput.function.size() != thisFuncInput.derivative.size())
{
- GMX_THROW(InconsistentInputError("Function and derivative vectors have different lengths"));
+ GMX_THROW(InconsistentInputError(
+ "Function and derivative vectors have different lengths"));
}
- internal::throwUnlessDerivativeIsConsistentWithFunction(thisFuncInput.function, thisFuncInput.derivative, thisFuncInput.spacing, range_);
+ internal::throwUnlessDerivativeIsConsistentWithFunction(
+ thisFuncInput.function, thisFuncInput.derivative, thisFuncInput.spacing, range_);
}
- catch (gmx::GromacsException &ex)
+ catch (gmx::GromacsException& ex)
{
- ex.prependContext("Error generating quadratic spline table for function '" + thisFuncInput.desc + "'");
+ ex.prependContext("Error generating quadratic spline table for function '"
+ + thisFuncInput.desc + "'");
throw;
}
// Calculate the required table spacing h. The error we make with linear interpolation
// errors by calculating the quotient of the function and second derivative of the
// input-derivative-analytical function instead.
- double thisMinQuotient = internal::findSmallestQuotientOfFunctionAndSecondDerivative(thisFuncInput.derivative, thisFuncInput.spacing, range_);
+ double thisMinQuotient = internal::findSmallestQuotientOfFunctionAndSecondDerivative(
+ thisFuncInput.derivative, thisFuncInput.spacing, range_);
minQuotient = std::min(minQuotient, thisMinQuotient);
}
- double spacing = std::sqrt(12.0 * tolerance * minQuotient);
+ double spacing = std::sqrt(12.0 * tolerance * minQuotient);
halfSpacing_ = 0.5 * spacing;
tableScale_ = 1.0 / spacing;
if (range_.second * tableScale_ > 1e6)
{
- GMX_THROW(ToleranceError("Requested tolerance would require over a million points in table"));
+ GMX_THROW(
+ ToleranceError("Requested tolerance would require over a million points in table"));
}
// Loop over all tables again.
{
if (spacing < thisFuncInput.spacing)
{
- GMX_THROW(ToleranceError("Input vector spacing cannot achieve tolerance requested"));
+ GMX_THROW(
+ ToleranceError("Input vector spacing cannot achieve tolerance requested"));
}
std::vector<real> tmpFuncTableData;
std::vector<real> tmpDerTableData;
std::vector<real> tmpDdfzTableData;
- fillSingleQuadraticSplineTableData(thisFuncInput.function,
- thisFuncInput.derivative,
- thisFuncInput.spacing,
- range,
- spacing,
- &tmpFuncTableData,
- &tmpDerTableData);
+ fillSingleQuadraticSplineTableData(thisFuncInput.function, thisFuncInput.derivative,
+ thisFuncInput.spacing, range, spacing,
+ &tmpFuncTableData, &tmpDerTableData);
fillDdfzTableData(tmpFuncTableData, tmpDerTableData, &tmpDdfzTableData);
- internal::fillMultiplexedTableData(tmpDerTableData, &derivativeMultiTableData_,
- 1, numFuncInTable_, funcIndex);
+ internal::fillMultiplexedTableData(tmpDerTableData, &derivativeMultiTableData_, 1,
+ numFuncInTable_, funcIndex);
- internal::fillMultiplexedTableData(tmpDdfzTableData, &ddfzMultiTableData_,
- 4, numFuncInTable_, funcIndex);
+ internal::fillMultiplexedTableData(tmpDdfzTableData, &ddfzMultiTableData_, 4,
+ numFuncInTable_, funcIndex);
funcIndex++;
}
- catch (gmx::GromacsException &ex)
+ catch (gmx::GromacsException& ex)
{
- ex.prependContext("Error generating quadratic spline table for function '" + thisFuncInput.desc + "'");
+ ex.prependContext("Error generating quadratic spline table for function '"
+ + thisFuncInput.desc + "'");
throw;
}
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff