Skip to content

GitLab

Commit 629a9631 authored by Peter Eastman's avatar Peter Eastman
Browse files

Created SplineFitter

parent e4ca66e4
Hide whitespace changes
Inline Side-by-side
Showing with 380 additions and 0 deletions
openmmapi/include/openmm/internal/SplineFitter.h 0 → 100644
View file @ 629a9631
#ifndef OPENMM_SPLINEFITTER_H_
#define OPENMM_SPLINEFITTER_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2010 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "windowsExport.h"
#include <vector>
namespace OpenMM {
/**
* SplineFitter provides routines for performing cubic spline interpolation.
*/
class SplineFitter {
public:
/**
* Fit a natural cubic spline to a set of data points. The resulting spline interpolates all the
* data points, has a continuous second derivative everywhere, and has a second derivative of 0 at
* its end points.
*
* @param x the values of the independent variable at the data points to interpolate. They must
* be strictly increasing: x[i] > x[i-1].
* @param y the values of the dependent variable at the data points to interpolate
* @param deriv on exit, this contains the second derivative of the spline at each of the data points
*/
static void createNaturalSpline(const std::vector<double>& x, const std::vector<double>& y, std::vector<double>& deriv);
/**
* Fit a periodic cubic spline to a set of data points. The resulting spline interpolates all the
* data points, has a continuous second derivative everywhere, and has identical second derivatives
* at the end points.
*
* @param x the values of the independent variable at the data points to interpolate. They must
* be strictly increasing: x[i] > x[i-1].
* @param y the values of the dependent variable at the data points to interpolate. The first and
* last entries must be identical.
* @param deriv on exit, this contains the second derivative of the spline at each of the data points
*/
static void createPeriodicSpline(const std::vector<double>& x, const std::vector<double>& y, std::vector<double>& deriv);
/**
* Evaluate a spline generated by one of the other methods in this class.
*
* @param x the values of the independent variable at the data points to interpolate
* @param y the values of the dependent variable at the data points to interpolate
* @param deriv the vector of second derivatives that was calculated by one of the other methods
* @param t the value of the independent variable at which to evaluate the spline
* @return the value of the spline at the specified point
*/
static double evaluateSpline(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& deriv, double t);
/**
* Evaluate the derivative of a spline generated by one of the other methods in this class.
*
* @param x the values of the independent variable at the data points to interpolate
* @param y the values of the dependent variable at the data points to interpolate
* @param deriv the vector of second derivatives that was calculated by one of the other methods
* @param t the value of the independent variable at which to evaluate the spline
* @return the value of the spline's derivative at the specified point
*/
static double evaluateSplineDerivative(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& deriv, double t);
private:
static void solveTridiagonalMatrix(const std::vector<double>& a, const std::vector<double>& b, const std::vector<double>& c, const std::vector<double>& rhs, std::vector<double>& sol);
};
} // namespace OpenMM
#endif /*OPENMM_SPLINEFITTER_H_*/
openmmapi/src/SplineFitter.cpp 0 → 100644
View file @ 629a9631
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2010 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include <vector>
#include "openmm/internal/SplineFitter.h"
#include "openmm/OpenMMException.h"
using namespace OpenMM;
using namespace std;
void SplineFitter::createNaturalSpline(const vector<double>& x, const vector<double>& y, vector<double>& deriv) {
int n = x.size();
if (y.size() != n)
throw OpenMMException("createNaturalSpline: x and y vectors must have same length");
if (n < 3)
throw OpenMMException("createNaturalSpline: the length of the input array must be at least 3");
deriv.resize(n);
// Create the system of equations to solve.
vector<double> a(n), b(n), c(n), rhs(n);
a[0] = 0.0;
b[0] = 1.0;
c[0] = 0.0;
rhs[0] = 0.0;
for (int i = 1; i < n-1; i++) {
a[i] = x[i]-x[i-1];
b[i] = 2.0*(x[i+1]-x[i-1]);
c[i] = x[i+1]-x[i];
rhs[i] = 6.0*((y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]));
}
a[n-1] = 0.0;
b[n-1] = 1.0;
c[n-1] = 0.0;
rhs[n-1] = 0.0;
// Solve them.
solveTridiagonalMatrix(a, b, c, rhs, deriv);
}
void SplineFitter::createPeriodicSpline(const vector<double>& x, const vector<double>& y, vector<double>& deriv) {
int n = x.size();
if (y.size() != n)
throw OpenMMException("createPeriodicSpline: x and y vectors must have same length");
if (n < 3)
throw OpenMMException("createPeriodicSpline: the length of the input array must be at least 3");
if (y[0] != y[n-1])
throw OpenMMException("createPeriodicSpline: the first and last points must have the same value");
deriv.resize(n);
// Create the system of equations to solve.
vector<double> a(n), b(n), c(n), rhs(n);
a[0] = 0.0;
b[0] = 2.0*(x[1]-x[0]+x[n-1]-x[n-2]);
c[0] = x[1]-x[0];
rhs[0] = 6.0*((y[1]-y[0])/(x[1]-x[0]) - (y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
for (int i = 1; i < n-1; i++) {
a[i] = x[i]-x[i-1];
b[i] = 2.0*(x[i+1]-x[i-1]);
c[i] = x[i+1]-x[i];
rhs[i] = 6.0*((y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]));
}
a[n-1] = 0.0;
b[n-1] = 1.0;
c[n-1] = 0.0;
rhs[n-1] = 0.0;
double beta = x[n-1]-x[n-2];
double alpha = -1.0;
double gamma = -b[0];
// This is a cyclic tridiagonal matrix. We solve it using the Sherman-Morrison method,
// which involves solving two tridiagonal systems.
b[0] -= gamma;
b[n-1] -= alpha*beta/gamma;
solveTridiagonalMatrix(a, b, c, rhs, deriv);
vector<double> u(n, 0.0), z(n);
u[0] = gamma;
u[n-1] = alpha;
solveTridiagonalMatrix(a, b, c, u, z);
double scale = (deriv[0]+beta*deriv[n-1]/gamma)/(1.0+z[0]+beta*z[n-1]/gamma);
for (int i = 0; i < n; i++)
deriv[i] -= scale*z[i];
}
double SplineFitter::evaluateSpline(const vector<double>& x, const vector<double>& y, const vector<double>& deriv, double t) {
int n = x.size();
if (t < x[0] || t > x[n-1])
throw OpenMMException("evaluateSpline: specified point is outside the range defined by the spline");
// Perform a binary search to identify the interval containing the point to evaluate.
int lower = 0;
int upper = n-1;
while (upper-lower > 1) {
int middle = (upper+lower)/2;
if (x[middle] > t)
upper = middle;
else
lower = middle;
}
// Evaluate the spline.
double dx = x[upper]-x[lower];
double a = (x[upper]-t)/dx;
double b = 1.0-a;
return a*y[lower]+b*y[upper]+((a*a*a-a)*deriv[lower] + (b*b*b-b)*deriv[upper])*dx*dx/6.0;
}
double SplineFitter::evaluateSplineDerivative(const vector<double>& x, const vector<double>& y, const vector<double>& deriv, double t) {
int n = x.size();
if (t < x[0] || t > x[n-1])
throw OpenMMException("evaluateSplineDerivative: specified point is outside the range defined by the spline");
// Perform a binary search to identify the interval containing the point to evaluate.
int lower = 0;
int upper = n-1;
while (upper-lower > 1) {
int middle = (upper+lower)/2;
if (x[middle] > t)
upper = middle;
else
lower = middle;
}
// Evaluate the spline.
double dx = x[upper]-x[lower];
double a = (x[upper]-t)/dx;
double b = 1.0-a;
double dadx = -1.0/dx;
return dadx*y[lower]-dadx*y[upper]+((1.0-3.0*a*a)*deriv[lower] + (3.0*b*b-1.0)*deriv[upper])*dx/6.0;
}
void SplineFitter::solveTridiagonalMatrix(const vector<double>& a, const vector<double>& b, const vector<double>& c, const vector<double>& rhs, vector<double>& sol) {
int n = a.size();
vector<double> gamma(n);
// Decompose the matrix.
sol[0] = rhs[0]/b[0];
double beta = b[0];
for (int i = 1; i < n; i++) {
gamma[i] = c[i-1]/beta;
beta = b[i]-a[i]*gamma[i];
sol[i] = (rhs[i]-a[i]*sol[i-1])/beta;
}
// Perform backsubstitation.
for (int i = n-2; i >= 0; i--)
sol[i] -= gamma[i+1]*sol[i+1];
}
tests/TestSplineFitter.cpp 0 → 100644
View file @ 629a9631
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2010 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#ifdef WIN32
#define _USE_MATH_DEFINES // Needed to get M_PI
#endif
#include "AssertionUtilities.h"
#include "openmm/internal/SplineFitter.h"
#include <cmath>
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
void testNaturalSpline() {
vector<double> x(20);
vector<double> y(20);
for (int i = 0; i < x.size(); i++) {
x[i] = 0.5*i+0.1*sin(i);
y[i] = sin(x[i]);
}
vector<double> deriv;
SplineFitter::createNaturalSpline(x, y, deriv);
ASSERT_EQUAL_TOL(deriv[0], 0.0, 1e-6);
ASSERT_EQUAL_TOL(deriv[deriv.size()-1], 0.0, 1e-6);
for (int i = 0; i < x.size(); i++) {
ASSERT_EQUAL_TOL(y[i], SplineFitter::evaluateSpline(x, y, deriv, x[i]), 1e-6);
ASSERT_EQUAL_TOL(cos(x[i]), SplineFitter::evaluateSplineDerivative(x, y, deriv, x[i]), 0.05);
}
for (int i = 1; i < 9; i++) {
ASSERT_EQUAL_TOL(sin(i), SplineFitter::evaluateSpline(x, y, deriv, i), 0.05);
ASSERT_EQUAL_TOL(cos(i), SplineFitter::evaluateSplineDerivative(x, y, deriv, i), 0.05);
}
}
void testPeriodicSpline() {
vector<double> x(26);
vector<double> y(26);
for (int i = 0; i < x.size()-1; i++) {
x[i] = 0.5*i+0.1*sin(i);
y[i] = sin(x[i]);
}
x[x.size()-1] = 4*M_PI;
y[y.size()-1] = y[0];
vector<double> deriv;
SplineFitter::createPeriodicSpline(x, y, deriv);
ASSERT_EQUAL_TOL(deriv[0], deriv[deriv.size()-1], 1e-6);
for (int i = 0; i < x.size(); i++) {
ASSERT_EQUAL_TOL(y[i], SplineFitter::evaluateSpline(x, y, deriv, x[i]), 1e-6);
ASSERT_EQUAL_TOL(cos(x[i]), SplineFitter::evaluateSplineDerivative(x, y, deriv, x[i]), 0.05);
}
for (int i = 1; i < 9; i++) {
ASSERT_EQUAL_TOL(sin(i), SplineFitter::evaluateSpline(x, y, deriv, i), 0.05);
ASSERT_EQUAL_TOL(cos(i), SplineFitter::evaluateSplineDerivative(x, y, deriv, i), 0.05);
}
}
int main() {
try {
testNaturalSpline();
testPeriodicSpline();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment