ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/math/CubicSpline.cpp
(Generate patch)

Comparing:
branches/development/src/math/CubicSpline.cpp (file contents), Revision 1479 by gezelter, Mon Jul 26 19:00:48 2010 UTC vs.
trunk/src/math/CubicSpline.cpp (file contents), Revision 2069 by gezelter, Thu Mar 5 16:30:23 2015 UTC

# Line 35 | Line 35
35   *                                                                      
36   * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37   * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 < * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 < * [4]  Vardeman & Gezelter, in progress (2009).                        
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39 > * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 > * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43   #include "math/CubicSpline.hpp"
43 – #include "utils/simError.h"
44   #include <cmath>
45 + #include <cassert>
46 + #include <cstdio>
47   #include <algorithm>
48 < #include <iostream>
48 > #include <numeric>
49  
50   using namespace OpenMD;
51   using namespace std;
52  
53 < CubicSpline::CubicSpline() : generated(false), isUniform(true) {}
53 > CubicSpline::CubicSpline() : isUniform(true), generated(false) {
54 >  x_.clear();
55 >  y_.clear();
56 > }
57  
58 < void CubicSpline::addPoint(RealType xp, RealType yp) {
59 <  data.push_back(make_pair(xp, yp));
58 > void CubicSpline::addPoint(const RealType xp, const RealType yp) {
59 >  x_.push_back(xp);
60 >  y_.push_back(yp);
61   }
62  
63   void CubicSpline::addPoints(const vector<RealType>& xps,
64                              const vector<RealType>& yps) {
65    
66 <  if (xps.size() != yps.size()) {
67 <    printf( painCave.errMsg,
68 <            "CubicSpline::addPoints was passed vectors of different length!\n");
69 <    painCave.severity = OPENMD_ERROR;
70 <    painCave.isFatal = 1;
65 <    simError();    
66 >  assert(xps.size() == yps.size());
67 >  
68 >  for (unsigned int i = 0; i < xps.size(); i++){
69 >    x_.push_back(xps[i]);
70 >    y_.push_back(yps[i]);
71    }
67 –
68 –  for (int i = 0; i < xps.size(); i++)
69 –    data.push_back(make_pair(xps[i], yps[i]));
72   }
73  
74   void CubicSpline::generate() {
75    // Calculate coefficients defining a smooth cubic interpolatory spline.
76    //
77    // class values constructed:
78 <  //   n   = number of data points.
79 <  //   x   = vector of independent variable values
80 <  //   y   = vector of dependent variable values
81 <  //   b   = vector of S'(x[i]) values.
82 <  //   c   = vector of S"(x[i])/2 values.
83 <  //   d   = vector of S'''(x[i]+)/6 values (i < n).
78 >  //   n   = number of data_ points.
79 >  //   x_  = vector of independent variable values
80 >  //   y_  = vector of dependent variable values
81 >  //   b   = vector of S'(x_[i]) values.
82 >  //   c   = vector of S"(x_[i])/2 values.
83 >  //   d   = vector of S'''(x_[i]+)/6 values (i < n).
84    // Local variables:  
85 <  
85 >
86    RealType fp1, fpn, h, p;
87    
88    // make sure the sizes match
89    
90 <  n = data.size();  
90 >  n = x_.size();  
91    b.resize(n);
92    c.resize(n);
93    d.resize(n);
# Line 95 | Line 97 | void CubicSpline::generate() {
97    bool sorted = true;
98    
99    for (int i = 1; i < n; i++) {
100 <    if ( (data[i].first - data[i-1].first ) <= 0.0 ) sorted = false;
100 >    if ( (x_[i] - x_[i-1] ) <= 0.0 ) sorted = false;
101    }
102    
103    // sort if necessary
104    
105 <  if (!sorted) sort(data.begin(), data.end());  
105 >  if (!sorted) {
106 >    vector<int> p = sort_permutation(x_);
107 >    x_ = apply_permutation(x_, p);
108 >    y_ = apply_permutation(y_, p);
109 >  }
110    
111 +  
112    // Calculate coefficients for the tridiagonal system: store
113    // sub-diagonal in B, diagonal in D, difference quotient in C.  
114    
115 <  b[0] = data[1].first - data[0].first;
116 <  c[0] = (data[1].second - data[0].second) / b[0];
115 >  b[0] = x_[1] - x_[0];
116 >  c[0] = (y_[1] - y_[0]) / b[0];
117    
118    if (n == 2) {
119  
120      // Assume the derivatives at both endpoints are zero. Another
121      // assumption could be made to have a linear interpolant between
122      // the two points.  In that case, the b coefficients below would be
123 <    // (data[1].second - data[0].second) / (data[1].first - data[0].first)
123 >    // (y_[1] - y_[0]) / (x_[1] - x_[0])
124      // and the c and d coefficients would both be zero.
125      b[0] = 0.0;
126 <    c[0] = -3.0 * pow((data[1].second - data[0].second) /
127 <                      (data[1].first-data[0].first), 2);
121 <    d[0] = -2.0 * pow((data[1].second - data[0].second) /
122 <                      (data[1].first-data[0].first), 3);
126 >    c[0] = -3.0 * pow((y_[1] - y_[0]) / (x_[1] - x_[0]), 2);
127 >    d[0] = -2.0 * pow((y_[1] - y_[0]) / (x_[1] - x_[0]), 3);
128      b[1] = b[0];
129      c[1] = 0.0;
130      d[1] = 0.0;
131 <    dx = 1.0 / (data[1].first - data[0].first);
131 >    dx = 1.0 / (x_[1] - x_[0]);
132      isUniform = true;
133      generated = true;
134      return;
# Line 132 | Line 137 | void CubicSpline::generate() {
137    d[0] = 2.0 * b[0];
138    
139    for (int i = 1; i < n-1; i++) {
140 <    b[i] = data[i+1].first - data[i].first;
140 >    b[i] = x_[i+1] - x_[i];
141      if ( fabs( b[i] - b[0] ) / b[0] > 1.0e-5) isUniform = false;
142 <    c[i] = (data[i+1].second - data[i].second) / b[i];
142 >    c[i] = (y_[i+1] - y_[i]) / b[i];
143      d[i] = 2.0 * (b[i] + b[i-1]);
144    }
145    
146    d[n-1] = 2.0 * b[n-2];
147    
148    // Calculate estimates for the end slopes using polynomials
149 <  // that interpolate the data nearest the end.
149 >  // that interpolate the data_ nearest the end.
150    
151    fp1 = c[0] - b[0]*(c[1] - c[0])/(b[0] + b[1]);
152    if (n > 3) fp1 = fp1 + b[0]*((b[0] + b[1]) * (c[2] - c[1]) /
153                                 (b[1] + b[2]) -
154 <                               c[1] + c[0]) / (data[3].first - data[0].first);
154 >                               c[1] + c[0]) / (x_[3] - x_[0]);
155    
156    fpn = c[n-2] + b[n-2]*(c[n-2] - c[n-3])/(b[n-3] + b[n-2]);
157 <
157 >  
158    if (n > 3)  fpn = fpn + b[n-2] *
159 <    (c[n-2] - c[n-3] - (b[n-3] + b[n-2]) *
160 <     (c[n-3] - c[n-4])/(b[n-3] + b[n-4]))/(data[n-1].first - data[n-4].first);
159 >                (c[n-2] - c[n-3] - (b[n-3] + b[n-2]) *
160 >                 (c[n-3] - c[n-4])/(b[n-3] + b[n-4])) /
161 >                (x_[n-1] - x_[n-4]);
162    
157 –  
163    // Calculate the right hand side and store it in C.
164    
165    c[n-1] = 3.0 * (fpn - c[n-2]);
# Line 178 | Line 183 | void CubicSpline::generate() {
183    // Calculate the coefficients defining the spline.
184    
185    for (int i = 0; i < n-1; i++) {
186 <    h = data[i+1].first - data[i].first;
186 >    h = x_[i+1] - x_[i];
187      d[i] = (c[i+1] - c[i]) / (3.0 * h);
188 <    b[i] = (data[i+1].second - data[i].second)/h - h * (c[i] + h * d[i]);
188 >    b[i] = (y_[i+1] - y_[i])/h - h * (c[i] + h * d[i]);
189    }
190    
191    b[n-1] = b[n-2] + h * (2.0 * c[n-2] + h * 3.0 * d[n-2]);
192    
193 <  if (isUniform) dx = 1.0 / (data[1].first - data[0].first);
193 >  if (isUniform) dx = 1.0 / (x_[1] - x_[0]);
194    
195    generated = true;
196    return;
197   }
198  
199 < RealType CubicSpline::getValueAt(RealType t) {
199 > RealType CubicSpline::getValueAt(const RealType& t) {
200    // Evaluate the spline at t using coefficients
201    //
202    // Input parameters
# Line 200 | Line 205 | RealType CubicSpline::getValueAt(RealType t) {
205    //   value of spline at t.
206    
207    if (!generated) generate();
203 –  RealType dt;
208    
209 <  if ( t < data[0].first || t > data[n-1].first ) {    
210 <    sprintf( painCave.errMsg,
207 <             "CubicSpline::getValueAt was passed a value outside the range of the spline!\n");
208 <    painCave.severity = OPENMD_ERROR;
209 <    painCave.isFatal = 1;
210 <    simError();    
211 <  }
209 >  assert(t >= x_.front());
210 >  assert(t <= x_.back());
211  
212    //  Find the interval ( x[j], x[j+1] ) that contains or is nearest
213    //  to t.
214  
215 <  int j;
215 >  if (isUniform) {    
216 >    
217 >    j = max(0, min(n-1, int((t - x_[0]) * dx)));  
218  
219 +  } else {
220 +
221 +    j = n-1;
222 +    
223 +    for (int i = 0; i < n; i++) {
224 +      if ( t < x_[i] ) {
225 +        j = i-1;
226 +        break;
227 +      }      
228 +    }
229 +  }
230 +  
231 +  //  Evaluate the cubic polynomial.
232 +  
233 +  dt = t - x_[j];
234 +  return y_[j] + dt*(b[j] + dt*(c[j] + dt*d[j]));  
235 + }
236 +
237 +
238 + void CubicSpline::getValueAt(const RealType& t, RealType& v) {
239 +  // Evaluate the spline at t using coefficients
240 +  //
241 +  // Input parameters
242 +  //   t = point where spline is to be evaluated.
243 +  // Output:
244 +  //   value of spline at t.
245 +  
246 +  if (!generated) generate();
247 +  
248 +  assert(t >= x_.front());
249 +  assert(t <= x_.back());
250 +
251 +  //  Find the interval ( x[j], x[j+1] ) that contains or is nearest
252 +  //  to t.
253 +
254    if (isUniform) {    
255      
256 <    j = max(0, min(n-1, int((t - data[0].first) * dx)));  
256 >    j = max(0, min(n-1, int((t - x_[0]) * dx)));  
257  
258    } else {
259  
260      j = n-1;
261      
262      for (int i = 0; i < n; i++) {
263 <      if ( t < data[i].first ) {
263 >      if ( t < x_[i] ) {
264          j = i-1;
265          break;
266        }      
# Line 233 | Line 269 | RealType CubicSpline::getValueAt(RealType t) {
269    
270    //  Evaluate the cubic polynomial.
271    
272 <  dt = t - data[j].first;
273 <  return data[j].second + dt*(b[j] + dt*(c[j] + dt*d[j]));
238 <  
272 >  dt = t - x_[j];
273 >  v = y_[j] + dt*(b[j] + dt*(c[j] + dt*d[j]));  
274   }
275  
276 + pair<RealType, RealType> CubicSpline::getLimits(){
277 +  if (!generated) generate();
278 +  return make_pair( x_.front(), x_.back() );
279 + }
280  
281 < pair<RealType, RealType> CubicSpline::getValueAndDerivativeAt(RealType t) {
281 > void CubicSpline::getValueAndDerivativeAt(const RealType& t, RealType& v,
282 >                                          RealType &dv) {
283    // Evaluate the spline and first derivative at t using coefficients
284    //
285    // Input parameters
286    //   t = point where spline is to be evaluated.
247 –  // Output:
248 –  //   pair containing value of spline at t and first derivative at t
287  
288    if (!generated) generate();
251 –  RealType dt;
289    
290 <  if ( t < data.front().first || t > data.back().first ) {    
291 <    sprintf( painCave.errMsg,
255 <             "CubicSpline::getValueAndDerivativeAt was passed a value outside the range of the spline!\n");
256 <    painCave.severity = OPENMD_ERROR;
257 <    painCave.isFatal = 1;
258 <    simError();    
259 <  }
290 >  assert(t >= x_.front());
291 >  assert(t <= x_.back());
292  
293    //  Find the interval ( x[j], x[j+1] ) that contains or is nearest
294    //  to t.
295  
264 –  int j;
265 –
296    if (isUniform) {    
297      
298 <    j = max(0, min(n-1, int((t - data[0].first) * dx)));  
298 >    j = max(0, min(n-1, int((t - x_[0]) * dx)));  
299  
300    } else {
301  
302      j = n-1;
303      
304      for (int i = 0; i < n; i++) {
305 <      if ( t < data[i].first ) {
305 >      if ( t < x_[i] ) {
306          j = i-1;
307          break;
308        }      
# Line 281 | Line 311 | pair<RealType, RealType> CubicSpline::getValueAndDeriv
311    
312    //  Evaluate the cubic polynomial.
313    
314 <  dt = t - data[j].first;
314 >  dt = t - x_[j];
315  
316 <  RealType yval = data[j].second + dt*(b[j] + dt*(c[j] + dt*d[j]));
317 <  RealType dydx = b[j] + dt*(2.0 * c[j] + 3.0 * dt * d[j]);
316 >  v = y_[j] + dt*(b[j] + dt*(c[j] + dt*d[j]));
317 >  dv = b[j] + dt*(2.0 * c[j] + 3.0 * dt * d[j]);
318 > }
319 >
320 > std::vector<int> CubicSpline::sort_permutation(std::vector<RealType>& v) {
321 >  std::vector<int> p(v.size());
322 >
323 >  // 6 lines to replace std::iota(p.begin(), p.end(), 0);
324 >  int value = 0;
325 >  std::vector<int>::iterator i;
326 >  for (i = p.begin(); i != p.end(); ++i) {
327 >    (*i) = value;
328 >    ++value;
329 >  }
330    
331 <  return make_pair(yval, dydx);
331 >  std::sort(p.begin(), p.end(), OpenMD::Comparator(v) );
332 >  return p;
333   }
334 +
335 + std::vector<RealType> CubicSpline::apply_permutation(std::vector<RealType> const& v,
336 +                                                     std::vector<int> const& p) {
337 +  int n = p.size();
338 +  std::vector<double> sorted_vec(n);
339 +  for (int i = 0; i < n; ++i) {
340 +    sorted_vec[i] = v[p[i]];
341 +  }
342 +  return sorted_vec;
343 + }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines