src/math/CubicSpline.cpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 *
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
 * research, please cite the appropriate papers when you publish your
 * work.  Good starting points are:
 *                                                                      
 * [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).                        
 */
 
#include "math/CubicSpline.hpp"
#include "utils/simError.h"
#include <cmath>
#include <algorithm>
#include <iostream>

using namespace OpenMD;
using namespace std;

CubicSpline::CubicSpline() : generated(false), isUniform(true) {}

void CubicSpline::addPoint(RealType xp, RealType yp) {
  data.push_back(make_pair(xp, yp));
}

void CubicSpline::addPoints(const vector<RealType>& xps, 
                            const vector<RealType>& 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();    
  }

  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).
  // Local variables:   
  
  RealType fp1, fpn, h, p;
  
  // make sure the sizes match
  
  n = data.size();  
  x.resize(n);
  y.resize(n);
  b.resize(n);
  c.resize(n);
  d.resize(n);
  
  // make sure we are monotonically increasing in x:
  
  bool sorted = true;
  
  for (int i = 1; i < n; i++) {
    if ( (data[i].first - data[i-1].first ) <= 0.0 ) sorted = false;
  }
  
  // sort if necessary
  
  if (!sorted) sort(data.begin(), data.end());  
  
  // Copy spline data out to separate arrays:
  
  for (int i = 0; i < n; i++) {
    x[i] = data[i].first;
    y[i] = data[i].second;
  }
  
  // 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];
  
  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)
    // 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);
    b[1] = b[0];
    c[1] = 0.0;
    d[1] = 0.0;
    dx = 1.0 / (data[1].first - data[0].first);
    isUniform = true;
    generated = true;
    return;
  }
  
  d[0] = 2.0 * b[0];
  
  for (int i = 1; i < n-1; i++) {
    b[i] = data[i+1].first - data[i].first;
    if ( fabs( b[i] - b[0] ) / b[0] > 1.0e-5) isUniform = false;
    c[i] = (data[i+1].second - data[i].second) / 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.
  
  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);
  
  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);
  
  
  // Calculate the right hand side and store it in C.
  
  c[n-1] = 3.0 * (fpn - c[n-2]);
  for (int i = n-2; i > 0; i--) 
    c[i] = 3.0 * (c[i] - c[i-1]);  
  c[0] = 3.0 * (c[0] - fp1);
  
  // Solve the tridiagonal system.
  
  for (int k = 1; k < n; k++) {
    p = b[k-1] / d[k-1];
    d[k] = d[k] - p*b[k-1];
    c[k] = c[k] - p*c[k-1];
  }
  
  c[n-1] = c[n-1] / d[n-1];
  
  for (int k = n-2; k >= 0; k--) 
    c[k] = (c[k] - b[k] * c[k+1]) / d[k];
  
  // Calculate the coefficients defining the spline.
  
  for (int i = 0; i < n-1; i++) {
    h = data[i+1].first - data[i].first;
    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[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); 
  
  generated = true;
  return;
}

RealType CubicSpline::getValueAt(RealType t) {
  // 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();
  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();    
  }

  //  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)));   

  } else { 

    j = n-1;
    
    for (int i = 0; i < n; i++) {
      if ( t < data[i].first ) {
        j = i-1;
        break;
      }      
    }
  }
  
  //  Evaluate the cubic polynomial.
  
  dt = t - data[j].first;
  return data[j].second + dt*(b[j] + dt*(c[j] + dt*d[j]));
  
}


pair<RealType, RealType> CubicSpline::getValueAndDerivativeAt(RealType t) {
  // 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();    
  }

  //  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)));   

  } else { 

    j = n-1;
    
    for (int i = 0; i < n; i++) {
      if ( t < data[i].first ) {
        j = i-1;
        break;
      }      
    }
  }
  
  //  Evaluate the cubic polynomial.
  
  dt = t - data[j].first;

  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]);
  
  return make_pair(yval, dydx);
}
Revision:	1475
Committed:	Wed Jul 21 17:51:22 2010 UTC (15 years, 3 months ago) by gezelter
Original Path:	*branches/development/src/math/CubicSpline.cpp*
File size:	8856 byte(s)
Log Message:	Added CubicSpline class to replace our interpolation F90 module
#	Content
1	/*
2	* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3	*
4	* The University of Notre Dame grants you ("Licensee") a
5	* non-exclusive, royalty free, license to use, modify and
6	* redistribute this software in source and binary code form, provided
7	* that the following conditions are met:
8	*
9	* 1. Redistributions of source code must retain the above copyright
10	* notice, this list of conditions and the following disclaimer.
11	*
12	* 2. Redistributions in binary form must reproduce the above copyright
13	* notice, this list of conditions and the following disclaimer in the
14	* documentation and/or other materials provided with the
15	* distribution.
16	*
17	* This software is provided "AS IS," without a warranty of any
18	* kind. All express or implied conditions, representations and
19	* warranties, including any implied warranty of merchantability,
20	* fitness for a particular purpose or non-infringement, are hereby
21	* excluded. The University of Notre Dame and its licensors shall not
22	* be liable for any damages suffered by licensee as a result of
23	* using, modifying or distributing the software or its
24	* derivatives. In no event will the University of Notre Dame or its
25	* licensors be liable for any lost revenue, profit or data, or for
26	* direct, indirect, special, consequential, incidental or punitive
27	* damages, however caused and regardless of the theory of liability,
28	* arising out of the use of or inability to use software, even if the
29	* University of Notre Dame has been advised of the possibility of
30	* such damages.
31	*
32	* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33	* research, please cite the appropriate papers when you publish your
34	* work. Good starting points are:
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).
40	*/
41
42	#include "math/CubicSpline.hpp"
43	#include "utils/simError.h"
44	#include <cmath>
45	#include <algorithm>
46	#include <iostream>
47
48	using namespace OpenMD;
49	using namespace std;
50
51	CubicSpline::CubicSpline() : generated(false), isUniform(true) {}
52
53	void CubicSpline::addPoint(RealType xp, RealType yp) {
54	data.push_back(make_pair(xp, yp));
55	}
56
57	void CubicSpline::addPoints(const vector<RealType>& xps,
58	const vector<RealType>& yps) {
59
60	if (xps.size() != yps.size()) {
61	printf( painCave.errMsg,
62	"CubicSpline::addPoints was passed vectors of different length!\n");
63	painCave.severity = OPENMD_ERROR;
64	painCave.isFatal = 1;
65	simError();
66	}
67
68	for (int i = 0; i < xps.size(); i++)
69	data.push_back(make_pair(xps[i], yps[i]));
70	}
71
72	void CubicSpline::generate() {
73	// Calculate coefficients defining a smooth cubic interpolatory spline.
74	//
75	// class values constructed:
76	// n = number of data points.
77	// x = vector of independent variable values
78	// y = vector of dependent variable values
79	// b = vector of S'(x[i]) values.
80	// c = vector of S"(x[i])/2 values.
81	// d = vector of S'''(x[i]+)/6 values (i < n).
82	// Local variables:
83
84	RealType fp1, fpn, h, p;
85
86	// make sure the sizes match
87
88	n = data.size();
89	x.resize(n);
90	y.resize(n);
91	b.resize(n);
92	c.resize(n);
93	d.resize(n);
94
95	// make sure we are monotonically increasing in x:
96
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;
101	}
102
103	// sort if necessary
104
105	if (!sorted) sort(data.begin(), data.end());
106
107	// Copy spline data out to separate arrays:
108
109	for (int i = 0; i < n; i++) {
110	x[i] = data[i].first;
111	y[i] = data[i].second;
112	}
113
114	// Calculate coefficients for the tridiagonal system: store
115	// sub-diagonal in B, diagonal in D, difference quotient in C.
116
117	b[0] = data[1].first - data[0].first;
118	c[0] = (data[1].second - data[0].second) / b[0];
119
120	if (n == 2) {
121
122	// Assume the derivatives at both endpoints are zero. Another
123	// assumption could be made to have a linear interpolant between
124	// the two points. In that case, the b coefficients below would be
125	// (data[1].second - data[0].second) / (data[1].first - data[0].first)
126	// and the c and d coefficients would both be zero.
127	b[0] = 0.0;
128	c[0] = -3.0 * pow((data[1].second - data[0].second) /
129	(data[1].first-data[0].first), 2);
130	d[0] = -2.0 * pow((data[1].second - data[0].second) /
131	(data[1].first-data[0].first), 3);
132	b[1] = b[0];
133	c[1] = 0.0;
134	d[1] = 0.0;
135	dx = 1.0 / (data[1].first - data[0].first);
136	isUniform = true;
137	generated = true;
138	return;
139	}
140
141	d[0] = 2.0 * b[0];
142
143	for (int i = 1; i < n-1; i++) {
144	b[i] = data[i+1].first - data[i].first;
145	if ( fabs( b[i] - b[0] ) / b[0] > 1.0e-5) isUniform = false;
146	c[i] = (data[i+1].second - data[i].second) / b[i];
147	d[i] = 2.0 * (b[i] + b[i-1]);
148	}
149
150	d[n-1] = 2.0 * b[n-2];
151
152	// Calculate estimates for the end slopes using polynomials
153	// that interpolate the data nearest the end.
154
155	fp1 = c[0] - b[0]*(c[1] - c[0])/(b[0] + b[1]);
156	if (n > 3) fp1 = fp1 + b[0]((b[0] + b[1]) (c[2] - c[1]) /
157	(b[1] + b[2]) -
158	c[1] + c[0]) / (data[3].first - data[0].first);
159
160	fpn = c[n-2] + b[n-2]*(c[n-2] - c[n-3])/(b[n-3] + b[n-2]);
161
162	if (n > 3) fpn = fpn + b[n-2] *
163	(c[n-2] - c[n-3] - (b[n-3] + b[n-2]) *
164	(c[n-3] - c[n-4])/(b[n-3] + b[n-4]))/(data[n-1].first - data[n-4].first);
165
166
167	// Calculate the right hand side and store it in C.
168
169	c[n-1] = 3.0 * (fpn - c[n-2]);
170	for (int i = n-2; i > 0; i--)
171	c[i] = 3.0 * (c[i] - c[i-1]);
172	c[0] = 3.0 * (c[0] - fp1);
173
174	// Solve the tridiagonal system.
175
176	for (int k = 1; k < n; k++) {
177	p = b[k-1] / d[k-1];
178	d[k] = d[k] - p*b[k-1];
179	c[k] = c[k] - p*c[k-1];
180	}
181
182	c[n-1] = c[n-1] / d[n-1];
183
184	for (int k = n-2; k >= 0; k--)
185	c[k] = (c[k] - b[k] * c[k+1]) / d[k];
186
187	// Calculate the coefficients defining the spline.
188
189	for (int i = 0; i < n-1; i++) {
190	h = data[i+1].first - data[i].first;
191	d[i] = (c[i+1] - c[i]) / (3.0 * h);
192	b[i] = (data[i+1].second - data[i].second)/h - h * (c[i] + h * d[i]);
193	}
194
195	b[n-1] = b[n-2] + h * (2.0 * c[n-2] + h * 3.0 * d[n-2]);
196
197	if (isUniform) dx = 1.0 / (data[1].first - data[0].first);
198
199	generated = true;
200	return;
201	}
202
203	RealType CubicSpline::getValueAt(RealType t) {
204	// Evaluate the spline at t using coefficients
205	//
206	// Input parameters
207	// t = point where spline is to be evaluated.
208	// Output:
209	// value of spline at t.
210
211	if (!generated) generate();
212	RealType dt;
213
214	if ( t < data[0].first \|\| t > data[n-1].first ) {
215	sprintf( painCave.errMsg,
216	"CubicSpline::getValueAt was passed a value outside the range of the spline!\n");
217	painCave.severity = OPENMD_ERROR;
218	painCave.isFatal = 1;
219	simError();
220	}
221
222	// Find the interval ( x[j], x[j+1] ) that contains or is nearest
223	// to t.
224
225	int j;
226
227	if (isUniform) {
228
229	j = max(0, min(n-1, int((t - data[0].first) * dx)));
230
231	} else {
232
233	j = n-1;
234
235	for (int i = 0; i < n; i++) {
236	if ( t < data[i].first ) {
237	j = i-1;
238	break;
239	}
240	}
241	}
242
243	// Evaluate the cubic polynomial.
244
245	dt = t - data[j].first;
246	return data[j].second + dt(b[j] + dt(c[j] + dt*d[j]));
247
248	}
249
250
251	pair<RealType, RealType> CubicSpline::getValueAndDerivativeAt(RealType t) {
252	// Evaluate the spline and first derivative at t using coefficients
253	//
254	// Input parameters
255	// t = point where spline is to be evaluated.
256	// Output:
257	// pair containing value of spline at t and first derivative at t
258
259	if (!generated) generate();
260	RealType dt;
261
262	if ( t < data.front().first \|\| t > data.back().first ) {
263	sprintf( painCave.errMsg,
264	"CubicSpline::getValueAndDerivativeAt was passed a value outside the range of the spline!\n");
265	painCave.severity = OPENMD_ERROR;
266	painCave.isFatal = 1;
267	simError();
268	}
269
270	// Find the interval ( x[j], x[j+1] ) that contains or is nearest
271	// to t.
272
273	int j;
274
275	if (isUniform) {
276
277	j = max(0, min(n-1, int((t - data[0].first) * dx)));
278
279	} else {
280
281	j = n-1;
282
283	for (int i = 0; i < n; i++) {
284	if ( t < data[i].first ) {
285	j = i-1;
286	break;
287	}
288	}
289	}
290
291	// Evaluate the cubic polynomial.
292
293	dt = t - data[j].first;
294
295	RealType yval = data[j].second + dt(b[j] + dt(c[j] + dt*d[j]));
296	RealType dydx = b[j] + dt(2.0 c[j] + 3.0 * dt * d[j]);
297
298	return make_pair(yval, dydx);
299	}