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

Comparing trunk/src/math/CubicSpline.cpp (file contents):
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

# Line 45 | Line 45
45   #include <cassert>
46   #include <cstdio>
47   #include <algorithm>
48 + #include <numeric>
49  
50   using namespace OpenMD;
51   using namespace std;
52  
53 < CubicSpline::CubicSpline() : generated(false), isUniform(true) {
54 <  data_.clear();
53 > CubicSpline::CubicSpline() : isUniform(true), generated(false) {
54 >  x_.clear();
55 >  y_.clear();
56   }
57  
58   void CubicSpline::addPoint(const RealType xp, const RealType yp) {
59 <  data_.push_back(make_pair(xp, yp));
59 >  x_.push_back(xp);
60 >  y_.push_back(yp);
61   }
62  
63   void CubicSpline::addPoints(const vector<RealType>& xps,
# Line 62 | Line 65 | void CubicSpline::addPoints(const vector<RealType>& xp
65    
66    assert(xps.size() == yps.size());
67    
68 <  for (unsigned int i = 0; i < xps.size(); i++)
69 <    data_.push_back(make_pair(xps[i], yps[i]));
68 >  for (unsigned int i = 0; i < xps.size(); i++){
69 >    x_.push_back(xps[i]);
70 >    y_.push_back(yps[i]);
71 >  }
72   }
73  
74   void CubicSpline::generate() {
# Line 71 | Line 76 | void CubicSpline::generate() {
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).
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  
86 <  RealType fp1, fpn, h, p;
86 >  RealType fp1, fpn, p;
87 >  RealType h(0.0);
88    
89    // make sure the sizes match
90    
91 <  n = data_.size();  
91 >  n = x_.size();  
92    b.resize(n);
93    c.resize(n);
94    d.resize(n);
# Line 92 | Line 98 | void CubicSpline::generate() {
98    bool sorted = true;
99    
100    for (int i = 1; i < n; i++) {
101 <    if ( (data_[i].first - data_[i-1].first ) <= 0.0 ) sorted = false;
101 >    if ( (x_[i] - x_[i-1] ) <= 0.0 ) sorted = false;
102    }
103    
104    // sort if necessary
105    
106 <  if (!sorted) sort(data_.begin(), data_.end());  
106 >  if (!sorted) {
107 >    vector<int> p = sort_permutation(x_);
108 >    x_ = apply_permutation(x_, p);
109 >    y_ = apply_permutation(y_, p);
110 >  }
111    
112 +  
113    // Calculate coefficients for the tridiagonal system: store
114    // sub-diagonal in B, diagonal in D, difference quotient in C.  
115    
116 <  b[0] = data_[1].first - data_[0].first;
117 <  c[0] = (data_[1].second - data_[0].second) / b[0];
116 >  b[0] = x_[1] - x_[0];
117 >  c[0] = (y_[1] - y_[0]) / b[0];
118    
119    if (n == 2) {
120  
121      // Assume the derivatives at both endpoints are zero. Another
122      // assumption could be made to have a linear interpolant between
123      // the two points.  In that case, the b coefficients below would be
124 <    // (data_[1].second - data_[0].second) / (data_[1].first - data_[0].first)
124 >    // (y_[1] - y_[0]) / (x_[1] - x_[0])
125      // and the c and d coefficients would both be zero.
126      b[0] = 0.0;
127 <    c[0] = -3.0 * pow((data_[1].second - data_[0].second) /
128 <                      (data_[1].first-data_[0].first), 2);
118 <    d[0] = -2.0 * pow((data_[1].second - data_[0].second) /
119 <                      (data_[1].first-data_[0].first), 3);
127 >    c[0] = -3.0 * pow((y_[1] - y_[0]) / (x_[1] - x_[0]), 2);
128 >    d[0] = -2.0 * pow((y_[1] - y_[0]) / (x_[1] - x_[0]), 3);
129      b[1] = b[0];
130      c[1] = 0.0;
131      d[1] = 0.0;
132 <    dx = 1.0 / (data_[1].first - data_[0].first);
132 >    dx = 1.0 / (x_[1] - x_[0]);
133      isUniform = true;
134      generated = true;
135      return;
# Line 129 | Line 138 | void CubicSpline::generate() {
138    d[0] = 2.0 * b[0];
139    
140    for (int i = 1; i < n-1; i++) {
141 <    b[i] = data_[i+1].first - data_[i].first;
141 >    b[i] = x_[i+1] - x_[i];
142      if ( fabs( b[i] - b[0] ) / b[0] > 1.0e-5) isUniform = false;
143 <    c[i] = (data_[i+1].second - data_[i].second) / b[i];
143 >    c[i] = (y_[i+1] - y_[i]) / b[i];
144      d[i] = 2.0 * (b[i] + b[i-1]);
145    }
146    
# Line 143 | Line 152 | void CubicSpline::generate() {
152    fp1 = c[0] - b[0]*(c[1] - c[0])/(b[0] + b[1]);
153    if (n > 3) fp1 = fp1 + b[0]*((b[0] + b[1]) * (c[2] - c[1]) /
154                                 (b[1] + b[2]) -
155 <                               c[1] + c[0]) / (data_[3].first - data_[0].first);
155 >                               c[1] + c[0]) / (x_[3] - x_[0]);
156    
157    fpn = c[n-2] + b[n-2]*(c[n-2] - c[n-3])/(b[n-3] + b[n-2]);
158    
159    if (n > 3)  fpn = fpn + b[n-2] *
160                  (c[n-2] - c[n-3] - (b[n-3] + b[n-2]) *
161                   (c[n-3] - c[n-4])/(b[n-3] + b[n-4])) /
162 <                (data_[n-1].first - data_[n-4].first);
162 >                (x_[n-1] - x_[n-4]);
163    
164    // Calculate the right hand side and store it in C.
165    
# Line 175 | Line 184 | void CubicSpline::generate() {
184    // Calculate the coefficients defining the spline.
185    
186    for (int i = 0; i < n-1; i++) {
187 <    h = data_[i+1].first - data_[i].first;
187 >    h = x_[i+1] - x_[i];
188      d[i] = (c[i+1] - c[i]) / (3.0 * h);
189 <    b[i] = (data_[i+1].second - data_[i].second)/h - h * (c[i] + h * d[i]);
189 >    b[i] = (y_[i+1] - y_[i])/h - h * (c[i] + h * d[i]);
190    }
191    
192    b[n-1] = b[n-2] + h * (2.0 * c[n-2] + h * 3.0 * d[n-2]);
193    
194 <  if (isUniform) dx = 1.0 / (data_[1].first - data_[0].first);
194 >  if (isUniform) dx = 1.0 / (x_[1] - x_[0]);
195    
196    generated = true;
197    return;
# Line 198 | Line 207 | RealType CubicSpline::getValueAt(const RealType& t) {
207    
208    if (!generated) generate();
209    
210 <  assert(t > data_.front().first);
211 <  assert(t < data_.back().first);
210 >  assert(t >= x_.front());
211 >  assert(t <= x_.back());
212  
213    //  Find the interval ( x[j], x[j+1] ) that contains or is nearest
214    //  to t.
215  
216    if (isUniform) {    
217      
218 <    j = max(0, min(n-1, int((t - data_[0].first) * dx)));  
218 >    j = max(0, min(n-1, int((t - x_[0]) * dx)));  
219  
220    } else {
221  
222      j = n-1;
223      
224      for (int i = 0; i < n; i++) {
225 <      if ( t < data_[i].first ) {
225 >      if ( t < x_[i] ) {
226          j = i-1;
227          break;
228        }      
# Line 222 | Line 231 | RealType CubicSpline::getValueAt(const RealType& t) {
231    
232    //  Evaluate the cubic polynomial.
233    
234 <  dt = t - data_[j].first;
235 <  return data_[j].second + dt*(b[j] + dt*(c[j] + dt*d[j]));  
234 >  dt = t - x_[j];
235 >  return y_[j] + dt*(b[j] + dt*(c[j] + dt*d[j]));  
236   }
237  
238  
# Line 237 | Line 246 | void CubicSpline::getValueAt(const RealType& t, RealTy
246    
247    if (!generated) generate();
248    
249 <  assert(t > data_.front().first);
250 <  assert(t < data_.back().first);
249 >  assert(t >= x_.front());
250 >  assert(t <= x_.back());
251  
252    //  Find the interval ( x[j], x[j+1] ) that contains or is nearest
253    //  to t.
254  
255    if (isUniform) {    
256      
257 <    j = max(0, min(n-1, int((t - data_[0].first) * dx)));  
257 >    j = max(0, min(n-1, int((t - x_[0]) * dx)));  
258  
259    } else {
260  
261      j = n-1;
262      
263      for (int i = 0; i < n; i++) {
264 <      if ( t < data_[i].first ) {
264 >      if ( t < x_[i] ) {
265          j = i-1;
266          break;
267        }      
# Line 261 | Line 270 | void CubicSpline::getValueAt(const RealType& t, RealTy
270    
271    //  Evaluate the cubic polynomial.
272    
273 <  dt = t - data_[j].first;
274 <  v = data_[j].second + dt*(b[j] + dt*(c[j] + dt*d[j]));  
273 >  dt = t - x_[j];
274 >  v = y_[j] + dt*(b[j] + dt*(c[j] + dt*d[j]));  
275   }
276  
277 + pair<RealType, RealType> CubicSpline::getLimits(){
278 +  if (!generated) generate();
279 +  return make_pair( x_.front(), x_.back() );
280 + }
281  
282 < pair<RealType, RealType> CubicSpline::getValueAndDerivativeAt(const RealType& t){
282 > void CubicSpline::getValueAndDerivativeAt(const RealType& t, RealType& v,
283 >                                          RealType &dv) {
284    // Evaluate the spline and first derivative at t using coefficients
285    //
286    // Input parameters
287    //   t = point where spline is to be evaluated.
274 –  // Output:
275 –  //   pair containing value of spline at t and first derivative at t
288  
289    if (!generated) generate();
290    
291 <  assert(t > data_.front().first);
292 <  assert(t < data_.back().first);
291 >  assert(t >= x_.front());
292 >  assert(t <= x_.back());
293  
294    //  Find the interval ( x[j], x[j+1] ) that contains or is nearest
295    //  to t.
296  
297    if (isUniform) {    
298      
299 <    j = max(0, min(n-1, int((t - data_[0].first) * dx)));  
299 >    j = max(0, min(n-1, int((t - x_[0]) * dx)));  
300  
301    } else {
302  
303      j = n-1;
304      
305      for (int i = 0; i < n; i++) {
306 <      if ( t < data_[i].first ) {
306 >      if ( t < x_[i] ) {
307          j = i-1;
308          break;
309        }      
# Line 300 | Line 312 | pair<RealType, RealType> CubicSpline::getValueAndDeriv
312    
313    //  Evaluate the cubic polynomial.
314    
315 <  dt = t - data_[j].first;
315 >  dt = t - x_[j];
316  
317 <  yval = data_[j].second + dt*(b[j] + dt*(c[j] + dt*d[j]));
318 <  dydx = b[j] + dt*(2.0 * c[j] + 3.0 * dt * d[j]);
307 <
308 <  return make_pair(yval, dydx);
317 >  v = y_[j] + dt*(b[j] + dt*(c[j] + dt*d[j]));
318 >  dv = b[j] + dt*(2.0 * c[j] + 3.0 * dt * d[j]);
319   }
320  
321 < void CubicSpline::getValueAndDerivativeAt(const RealType& t, RealType& v,
322 <                                          RealType &dv) {
313 <  // Evaluate the spline and first derivative at t using coefficients
314 <  //
315 <  // Input parameters
316 <  //   t = point where spline is to be evaluated.
321 > std::vector<int> CubicSpline::sort_permutation(std::vector<RealType>& v) {
322 >  std::vector<int> p(v.size());
323  
324 <  if (!generated) generate();
325 <  
326 <  assert(t > data_.front().first);
327 <  assert(t < data_.back().first);
328 <
329 <  //  Find the interval ( x[j], x[j+1] ) that contains or is nearest
324 <  //  to t.
325 <
326 <  if (isUniform) {    
327 <    
328 <    j = max(0, min(n-1, int((t - data_[0].first) * dx)));  
329 <
330 <  } else {
331 <
332 <    j = n-1;
333 <    
334 <    for (int i = 0; i < n; i++) {
335 <      if ( t < data_[i].first ) {
336 <        j = i-1;
337 <        break;
338 <      }      
339 <    }
324 >  // 6 lines to replace std::iota(p.begin(), p.end(), 0);
325 >  int value = 0;
326 >  std::vector<int>::iterator i;
327 >  for (i = p.begin(); i != p.end(); ++i) {
328 >    (*i) = value;
329 >    ++value;
330    }
331    
332 <  //  Evaluate the cubic polynomial.
333 <  
334 <  dt = t - data_[j].first;
332 >  std::sort(p.begin(), p.end(), OpenMD::Comparator(v) );
333 >  return p;
334 > }
335  
336 <  v = data_[j].second + dt*(b[j] + dt*(c[j] + dt*d[j]));
337 <  dv = b[j] + dt*(2.0 * c[j] + 3.0 * dt * d[j]);
336 > std::vector<RealType> CubicSpline::apply_permutation(std::vector<RealType> const& v,
337 >                                                     std::vector<int> const& p) {
338 >  int n = p.size();
339 >  std::vector<double> sorted_vec(n);
340 >  for (int i = 0; i < n; ++i) {
341 >    sorted_vec[i] = v[p[i]];
342 >  }
343 >  return sorted_vec;
344   }

Diff Legend

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