--- branches/development/src/math/CubicSpline.cpp 2010/07/26 19:00:48 1479 +++ trunk/src/math/CubicSpline.cpp 2015/03/07 21:41:51 2071 @@ -35,57 +35,60 @@ * * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). - * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). - * [4] Vardeman & Gezelter, in progress (2009). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ #include "math/CubicSpline.hpp" -#include "utils/simError.h" #include +#include +#include #include -#include +#include using namespace OpenMD; using namespace std; -CubicSpline::CubicSpline() : generated(false), isUniform(true) {} +CubicSpline::CubicSpline() : isUniform(true), generated(false) { + x_.clear(); + y_.clear(); +} -void CubicSpline::addPoint(RealType xp, RealType yp) { - data.push_back(make_pair(xp, yp)); +void CubicSpline::addPoint(const RealType xp, const RealType yp) { + x_.push_back(xp); + y_.push_back(yp); } void CubicSpline::addPoints(const vector& xps, const vector& yps) { - if (xps.size() != yps.size()) { - printf( painCave.errMsg, - "CubicSpline::addPoints was passed vectors of different length!\n"); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); + assert(xps.size() == yps.size()); + + for (unsigned int i = 0; i < xps.size(); i++){ + x_.push_back(xps[i]); + y_.push_back(yps[i]); } - - for (int i = 0; i < xps.size(); i++) - data.push_back(make_pair(xps[i], yps[i])); } void CubicSpline::generate() { // Calculate coefficients defining a smooth cubic interpolatory spline. // // class values constructed: - // n = number of data points. - // x = vector of independent variable values - // y = vector of dependent variable values - // b = vector of S'(x[i]) values. - // c = vector of S"(x[i])/2 values. - // d = vector of S'''(x[i]+)/6 values (i < n). + // n = number of data_ points. + // x_ = vector of independent variable values + // y_ = vector of dependent variable values + // b = vector of S'(x_[i]) values. + // c = vector of S"(x_[i])/2 values. + // d = vector of S'''(x_[i]+)/6 values (i < n). // Local variables: + + RealType fp1, fpn, p; + RealType h(0.0); - RealType fp1, fpn, h, p; - // make sure the sizes match - n = data.size(); + n = x_.size(); b.resize(n); c.resize(n); d.resize(n); @@ -95,35 +98,38 @@ void CubicSpline::generate() { bool sorted = true; for (int i = 1; i < n; i++) { - if ( (data[i].first - data[i-1].first ) <= 0.0 ) sorted = false; + if ( (x_[i] - x_[i-1] ) <= 0.0 ) sorted = false; } // sort if necessary - if (!sorted) sort(data.begin(), data.end()); + if (!sorted) { + vector p = sort_permutation(x_); + x_ = apply_permutation(x_, p); + y_ = apply_permutation(y_, p); + } + // Calculate coefficients for the tridiagonal system: store // sub-diagonal in B, diagonal in D, difference quotient in C. - b[0] = data[1].first - data[0].first; - c[0] = (data[1].second - data[0].second) / b[0]; + b[0] = x_[1] - x_[0]; + c[0] = (y_[1] - y_[0]) / b[0]; if (n == 2) { // Assume the derivatives at both endpoints are zero. Another // assumption could be made to have a linear interpolant between // the two points. In that case, the b coefficients below would be - // (data[1].second - data[0].second) / (data[1].first - data[0].first) + // (y_[1] - y_[0]) / (x_[1] - x_[0]) // and the c and d coefficients would both be zero. b[0] = 0.0; - c[0] = -3.0 * pow((data[1].second - data[0].second) / - (data[1].first-data[0].first), 2); - d[0] = -2.0 * pow((data[1].second - data[0].second) / - (data[1].first-data[0].first), 3); + c[0] = -3.0 * pow((y_[1] - y_[0]) / (x_[1] - x_[0]), 2); + d[0] = -2.0 * pow((y_[1] - y_[0]) / (x_[1] - x_[0]), 3); b[1] = b[0]; c[1] = 0.0; d[1] = 0.0; - dx = 1.0 / (data[1].first - data[0].first); + dx = 1.0 / (x_[1] - x_[0]); isUniform = true; generated = true; return; @@ -132,29 +138,29 @@ void CubicSpline::generate() { d[0] = 2.0 * b[0]; for (int i = 1; i < n-1; i++) { - b[i] = data[i+1].first - data[i].first; + b[i] = x_[i+1] - x_[i]; if ( fabs( b[i] - b[0] ) / b[0] > 1.0e-5) isUniform = false; - c[i] = (data[i+1].second - data[i].second) / b[i]; + c[i] = (y_[i+1] - y_[i]) / b[i]; d[i] = 2.0 * (b[i] + b[i-1]); } d[n-1] = 2.0 * b[n-2]; // Calculate estimates for the end slopes using polynomials - // that interpolate the data nearest the end. + // that interpolate the data_ nearest the end. fp1 = c[0] - b[0]*(c[1] - c[0])/(b[0] + b[1]); if (n > 3) fp1 = fp1 + b[0]*((b[0] + b[1]) * (c[2] - c[1]) / (b[1] + b[2]) - - c[1] + c[0]) / (data[3].first - data[0].first); + c[1] + c[0]) / (x_[3] - x_[0]); fpn = c[n-2] + b[n-2]*(c[n-2] - c[n-3])/(b[n-3] + b[n-2]); - + if (n > 3) fpn = fpn + b[n-2] * - (c[n-2] - c[n-3] - (b[n-3] + b[n-2]) * - (c[n-3] - c[n-4])/(b[n-3] + b[n-4]))/(data[n-1].first - data[n-4].first); + (c[n-2] - c[n-3] - (b[n-3] + b[n-2]) * + (c[n-3] - c[n-4])/(b[n-3] + b[n-4])) / + (x_[n-1] - x_[n-4]); - // Calculate the right hand side and store it in C. c[n-1] = 3.0 * (fpn - c[n-2]); @@ -178,20 +184,20 @@ void CubicSpline::generate() { // Calculate the coefficients defining the spline. for (int i = 0; i < n-1; i++) { - h = data[i+1].first - data[i].first; + h = x_[i+1] - x_[i]; d[i] = (c[i+1] - c[i]) / (3.0 * h); - b[i] = (data[i+1].second - data[i].second)/h - h * (c[i] + h * d[i]); + b[i] = (y_[i+1] - y_[i])/h - h * (c[i] + h * d[i]); } b[n-1] = b[n-2] + h * (2.0 * c[n-2] + h * 3.0 * d[n-2]); - if (isUniform) dx = 1.0 / (data[1].first - data[0].first); + if (isUniform) dx = 1.0 / (x_[1] - x_[0]); generated = true; return; } -RealType CubicSpline::getValueAt(RealType t) { +RealType CubicSpline::getValueAt(const RealType& t) { // Evaluate the spline at t using coefficients // // Input parameters @@ -200,31 +206,23 @@ RealType CubicSpline::getValueAt(RealType t) { // value of spline at t. if (!generated) generate(); - RealType dt; - if ( t < data[0].first || t > data[n-1].first ) { - sprintf( painCave.errMsg, - "CubicSpline::getValueAt was passed a value outside the range of the spline!\n"); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } + assert(t >= x_.front()); + assert(t <= x_.back()); // Find the interval ( x[j], x[j+1] ) that contains or is nearest // to t. - int j; - if (isUniform) { - j = max(0, min(n-1, int((t - data[0].first) * dx))); + j = max(0, min(n-1, int((t - x_[0]) * dx))); } else { j = n-1; for (int i = 0; i < n; i++) { - if ( t < data[i].first ) { + if ( t < x_[i] ) { j = i-1; break; } @@ -233,46 +231,79 @@ RealType CubicSpline::getValueAt(RealType t) { // Evaluate the cubic polynomial. - dt = t - data[j].first; - return data[j].second + dt*(b[j] + dt*(c[j] + dt*d[j])); + dt = t - x_[j]; + return y_[j] + dt*(b[j] + dt*(c[j] + dt*d[j])); +} + + +void CubicSpline::getValueAt(const RealType& t, RealType& v) { + // Evaluate the spline at t using coefficients + // + // Input parameters + // t = point where spline is to be evaluated. + // Output: + // value of spline at t. + if (!generated) generate(); + + assert(t >= x_.front()); + assert(t <= x_.back()); + + // Find the interval ( x[j], x[j+1] ) that contains or is nearest + // to t. + + if (isUniform) { + + j = max(0, min(n-1, int((t - x_[0]) * dx))); + + } else { + + j = n-1; + + for (int i = 0; i < n; i++) { + if ( t < x_[i] ) { + j = i-1; + break; + } + } + } + + // Evaluate the cubic polynomial. + + dt = t - x_[j]; + v = y_[j] + dt*(b[j] + dt*(c[j] + dt*d[j])); } +pair CubicSpline::getLimits(){ + if (!generated) generate(); + return make_pair( x_.front(), x_.back() ); +} -pair CubicSpline::getValueAndDerivativeAt(RealType t) { +void CubicSpline::getValueAndDerivativeAt(const RealType& t, RealType& v, + RealType &dv) { // Evaluate the spline and first derivative at t using coefficients // // Input parameters // t = point where spline is to be evaluated. - // Output: - // pair containing value of spline at t and first derivative at t if (!generated) generate(); - RealType dt; - if ( t < data.front().first || t > data.back().first ) { - sprintf( painCave.errMsg, - "CubicSpline::getValueAndDerivativeAt was passed a value outside the range of the spline!\n"); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } + assert(t >= x_.front()); + assert(t <= x_.back()); // Find the interval ( x[j], x[j+1] ) that contains or is nearest // to t. - int j; - if (isUniform) { - j = max(0, min(n-1, int((t - data[0].first) * dx))); + j = max(0, min(n-1, int((t - x_[0]) * dx))); } else { j = n-1; for (int i = 0; i < n; i++) { - if ( t < data[i].first ) { + if ( t < x_[i] ) { j = i-1; break; } @@ -281,10 +312,33 @@ pair CubicSpline::getValueAndDeriv // Evaluate the cubic polynomial. - dt = t - data[j].first; + dt = t - x_[j]; - RealType yval = data[j].second + dt*(b[j] + dt*(c[j] + dt*d[j])); - RealType dydx = b[j] + dt*(2.0 * c[j] + 3.0 * dt * d[j]); + v = y_[j] + dt*(b[j] + dt*(c[j] + dt*d[j])); + dv = b[j] + dt*(2.0 * c[j] + 3.0 * dt * d[j]); +} + +std::vector CubicSpline::sort_permutation(std::vector& v) { + std::vector p(v.size()); + + // 6 lines to replace std::iota(p.begin(), p.end(), 0); + int value = 0; + std::vector::iterator i; + for (i = p.begin(); i != p.end(); ++i) { + (*i) = value; + ++value; + } - return make_pair(yval, dydx); + std::sort(p.begin(), p.end(), OpenMD::Comparator(v) ); + return p; } + +std::vector CubicSpline::apply_permutation(std::vector const& v, + std::vector const& p) { + int n = p.size(); + std::vector sorted_vec(n); + for (int i = 0; i < n; ++i) { + sorted_vec[i] = v[p[i]]; + } + return sorted_vec; +}