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, 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 <cmath>
#include <cassert>
#include <cstdio>
#include <algorithm>
#include <numeric>

using namespace OpenMD;
using namespace std;

CubicSpline::CubicSpline() : isUniform(true), generated(false) {
  x_.clear();
  y_.clear();
}

void CubicSpline::addPoint(const RealType xp, const RealType yp) {
  x_.push_back(xp);
  y_.push_back(yp);
}

void CubicSpline::addPoints(const vector<RealType>& xps, 
                            const vector<RealType>& yps) {
  
  assert(xps.size() == yps.size());
  
  for (unsigned int i = 0; i < xps.size(); i++){
    x_.push_back(xps[i]);
    y_.push_back(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, p;
  RealType h(0.0);
  
  // make sure the sizes match
  
  n = x_.size();  
  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 ( (x_[i] - x_[i-1] ) <= 0.0 ) sorted = false;
  }
  
  // sort if necessary
  
  if (!sorted) {
    vector<int> 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] = 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
    // (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((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 / (x_[1] - x_[0]);
    isUniform = true;
    generated = true;
    return;
  }
  
  d[0] = 2.0 * b[0];
  
  for (int i = 1; i < n-1; i++) {
    b[i] = x_[i+1] - x_[i];
    if ( fabs( b[i] - b[0] ) / b[0] > 1.0e-5) isUniform = false;
    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.
  
  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]) / (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])) /
                (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]);
  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 = x_[i+1] - x_[i];
    d[i] = (c[i+1] - c[i]) / (3.0 * h);
    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 / (x_[1] - x_[0]); 
  
  generated = true;
  return;
}

RealType CubicSpline::getValueAt(const 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();
  
  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];
  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<RealType, RealType> CubicSpline::getLimits(){
  if (!generated) generate();
  return make_pair( x_.front(), x_.back() );
}

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.

  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]));
  dv = b[j] + dt*(2.0 * c[j] + 3.0 * dt * d[j]); 
}

std::vector<int> CubicSpline::sort_permutation(std::vector<RealType>& v) {
  std::vector<int> p(v.size());

  // 6 lines to replace std::iota(p.begin(), p.end(), 0);
  int value = 0;
  std::vector<int>::iterator i;
  for (i = p.begin(); i != p.end(); ++i) {
    (*i) = value;
    ++value;
  }
  
  std::sort(p.begin(), p.end(), OpenMD::Comparator(v) );
  return p;
}

std::vector<RealType> CubicSpline::apply_permutation(std::vector<RealType> const& v,
                                                     std::vector<int> const& p) { 
  int n = p.size();
  std::vector<double> sorted_vec(n);
  for (int i = 0; i < n; ++i) {
    sorted_vec[i] = v[p[i]];
  }
  return sorted_vec;
}
Revision:	2071
Committed:	Sat Mar 7 21:41:51 2015 UTC (10 years, 1 month ago) by gezelter
File size:	9359 byte(s)
Log Message:	Reducing the number of warnings when using g++ to compile.
#	User	Rev	Content
1	gezelter	1475	/*
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	gezelter	1879	* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
39	gezelter	1665	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40			* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	gezelter	1475	*/
42
43			#include "math/CubicSpline.hpp"
44			#include <cmath>
45	gezelter	1879	#include <cassert>
46	gezelter	1618	#include <cstdio>
47	gezelter	1475	#include <algorithm>
48	gezelter	2057	#include <numeric>
49	gezelter	1475
50			using namespace OpenMD;
51			using namespace std;
52
53	gezelter	2069	CubicSpline::CubicSpline() : isUniform(true), generated(false) {
54	gezelter	2057	x_.clear();
55			y_.clear();
56	gezelter	1536	}
57	gezelter	1475
58	gezelter	1536	void CubicSpline::addPoint(const RealType xp, const RealType yp) {
59	gezelter	2057	x_.push_back(xp);
60			y_.push_back(yp);
61	gezelter	1475	}
62
63			void CubicSpline::addPoints(const vector<RealType>& xps,
64			const vector<RealType>& yps) {
65
66	gezelter	1879	assert(xps.size() == yps.size());
67
68	gezelter	2057	for (unsigned int i = 0; i < xps.size(); i++){
69			x_.push_back(xps[i]);
70			y_.push_back(yps[i]);
71			}
72	gezelter	1475	}
73
74			void CubicSpline::generate() {
75			// Calculate coefficients defining a smooth cubic interpolatory spline.
76			//
77			// class values constructed:
78	gezelter	1536	// n = number of data_ points.
79	gezelter	2057	// 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	gezelter	1475	// Local variables:
85	gezelter	1536
86	gezelter	2071	RealType fp1, fpn, p;
87			RealType h(0.0);
88	gezelter	1475
89			// make sure the sizes match
90
91	gezelter	2057	n = x_.size();
92	gezelter	1475	b.resize(n);
93			c.resize(n);
94			d.resize(n);
95
96			// make sure we are monotonically increasing in x:
97
98			bool sorted = true;
99
100			for (int i = 1; i < n; i++) {
101	gezelter	2057	if ( (x_[i] - x_[i-1] ) <= 0.0 ) sorted = false;
102	gezelter	1475	}
103
104			// sort if necessary
105
106	gezelter	2057	if (!sorted) {
107			vector<int> p = sort_permutation(x_);
108			x_ = apply_permutation(x_, p);
109			y_ = apply_permutation(y_, p);
110			}
111	gezelter	1475
112	gezelter	2057
113	gezelter	1475	// Calculate coefficients for the tridiagonal system: store
114			// sub-diagonal in B, diagonal in D, difference quotient in C.
115
116	gezelter	2057	b[0] = x_[1] - x_[0];
117			c[0] = (y_[1] - y_[0]) / b[0];
118	gezelter	1475
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	gezelter	2057	// (y_[1] - y_[0]) / (x_[1] - x_[0])
125	gezelter	1475	// and the c and d coefficients would both be zero.
126			b[0] = 0.0;
127	gezelter	2057	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	gezelter	1475	b[1] = b[0];
130			c[1] = 0.0;
131			d[1] = 0.0;
132	gezelter	2057	dx = 1.0 / (x_[1] - x_[0]);
133	gezelter	1475	isUniform = true;
134			generated = true;
135			return;
136			}
137
138			d[0] = 2.0 * b[0];
139
140			for (int i = 1; i < n-1; i++) {
141	gezelter	2057	b[i] = x_[i+1] - x_[i];
142	gezelter	1475	if ( fabs( b[i] - b[0] ) / b[0] > 1.0e-5) isUniform = false;
143	gezelter	2057	c[i] = (y_[i+1] - y_[i]) / b[i];
144	gezelter	1475	d[i] = 2.0 * (b[i] + b[i-1]);
145			}
146
147			d[n-1] = 2.0 * b[n-2];
148
149			// Calculate estimates for the end slopes using polynomials
150	gezelter	1536	// that interpolate the data_ nearest the end.
151	gezelter	1475
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	gezelter	2057	c[1] + c[0]) / (x_[3] - x_[0]);
156	gezelter	1475
157			fpn = c[n-2] + b[n-2]*(c[n-2] - c[n-3])/(b[n-3] + b[n-2]);
158	gezelter	1879
159	gezelter	1475	if (n > 3) fpn = fpn + b[n-2] *
160	gezelter	1879	(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	gezelter	2057	(x_[n-1] - x_[n-4]);
163	gezelter	1475
164			// Calculate the right hand side and store it in C.
165
166			c[n-1] = 3.0 * (fpn - c[n-2]);
167			for (int i = n-2; i > 0; i--)
168			c[i] = 3.0 * (c[i] - c[i-1]);
169			c[0] = 3.0 * (c[0] - fp1);
170
171			// Solve the tridiagonal system.
172
173			for (int k = 1; k < n; k++) {
174			p = b[k-1] / d[k-1];
175			d[k] = d[k] - p*b[k-1];
176			c[k] = c[k] - p*c[k-1];
177			}
178
179			c[n-1] = c[n-1] / d[n-1];
180
181			for (int k = n-2; k >= 0; k--)
182			c[k] = (c[k] - b[k] * c[k+1]) / d[k];
183
184			// Calculate the coefficients defining the spline.
185
186			for (int i = 0; i < n-1; i++) {
187	gezelter	2057	h = x_[i+1] - x_[i];
188	gezelter	1475	d[i] = (c[i+1] - c[i]) / (3.0 * h);
189	gezelter	2057	b[i] = (y_[i+1] - y_[i])/h - h * (c[i] + h * d[i]);
190	gezelter	1475	}
191
192			b[n-1] = b[n-2] + h * (2.0 * c[n-2] + h * 3.0 * d[n-2]);
193
194	gezelter	2057	if (isUniform) dx = 1.0 / (x_[1] - x_[0]);
195	gezelter	1475
196			generated = true;
197			return;
198			}
199
200	gezelter	1879	RealType CubicSpline::getValueAt(const RealType& t) {
201	gezelter	1475	// Evaluate the spline at t using coefficients
202			//
203			// Input parameters
204			// t = point where spline is to be evaluated.
205			// Output:
206			// value of spline at t.
207
208			if (!generated) generate();
209
210	gezelter	2057	assert(t >= x_.front());
211			assert(t <= x_.back());
212	gezelter	1879
213			// Find the interval ( x[j], x[j+1] ) that contains or is nearest
214			// to t.
215
216			if (isUniform) {
217
218	gezelter	2057	j = max(0, min(n-1, int((t - x_[0]) * dx)));
219	gezelter	1879
220			} else {
221
222			j = n-1;
223
224			for (int i = 0; i < n; i++) {
225	gezelter	2057	if ( t < x_[i] ) {
226	gezelter	1879	j = i-1;
227			break;
228			}
229			}
230	gezelter	1475	}
231	gezelter	1879
232			// Evaluate the cubic polynomial.
233
234	gezelter	2057	dt = t - x_[j];
235			return y_[j] + dt(b[j] + dt(c[j] + dt*d[j]));
236	gezelter	1879	}
237	gezelter	1475
238	gezelter	1879
239			void CubicSpline::getValueAt(const RealType& t, RealType& v) {
240			// Evaluate the spline at t using coefficients
241			//
242			// Input parameters
243			// t = point where spline is to be evaluated.
244			// Output:
245			// value of spline at t.
246
247			if (!generated) generate();
248
249	gezelter	2057	assert(t >= x_.front());
250			assert(t <= x_.back());
251	gezelter	1879
252	gezelter	1475	// Find the interval ( x[j], x[j+1] ) that contains or is nearest
253			// to t.
254
255			if (isUniform) {
256
257	gezelter	2057	j = max(0, min(n-1, int((t - x_[0]) * dx)));
258	gezelter	1475
259			} else {
260
261			j = n-1;
262
263			for (int i = 0; i < n; i++) {
264	gezelter	2057	if ( t < x_[i] ) {
265	gezelter	1475	j = i-1;
266			break;
267			}
268			}
269			}
270
271			// Evaluate the cubic polynomial.
272
273	gezelter	2057	dt = t - x_[j];
274			v = y_[j] + dt(b[j] + dt(c[j] + dt*d[j]));
275	gezelter	1475	}
276
277	gezelter	1928	pair<RealType, RealType> CubicSpline::getLimits(){
278			if (!generated) generate();
279	gezelter	2057	return make_pair( x_.front(), x_.back() );
280	gezelter	1928	}
281
282	gezelter	1879	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.
288
289			if (!generated) generate();
290
291	gezelter	2057	assert(t >= x_.front());
292			assert(t <= x_.back());
293	gezelter	1879
294			// Find the interval ( x[j], x[j+1] ) that contains or is nearest
295			// to t.
296
297			if (isUniform) {
298
299	gezelter	2057	j = max(0, min(n-1, int((t - x_[0]) * dx)));
300	gezelter	1879
301			} else {
302
303			j = n-1;
304
305			for (int i = 0; i < n; i++) {
306	gezelter	2057	if ( t < x_[i] ) {
307	gezelter	1879	j = i-1;
308			break;
309			}
310			}
311			}
312
313			// Evaluate the cubic polynomial.
314
315	gezelter	2057	dt = t - x_[j];
316	gezelter	1879
317	gezelter	2057	v = y_[j] + dt(b[j] + dt(c[j] + dt*d[j]));
318	gezelter	1879	dv = b[j] + dt(2.0 c[j] + 3.0 * dt * d[j]);
319			}
320	gezelter	2057
321			std::vector<int> CubicSpline::sort_permutation(std::vector<RealType>& v) {
322			std::vector<int> p(v.size());
323	gezelter	2058
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	gezelter	2057	std::sort(p.begin(), p.end(), OpenMD::Comparator(v) );
333			return p;
334			}
335
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			}