research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Move definition of polyfit to polyfit.cpp

Browse source
This commit is contained in:
Ray Speth 2015-08-01 20:57:36 -04:00
parent acdf9cf0ed
commit 9cdfd12172
2 changed files with 52 additions and 97 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

97
src/numerics/funcs.cpp
View file

@ -8,27 +8,9 @@
*/
#include "cantera/numerics/funcs.h"
#include "cantera/base/ctexceptions.h"
#include "cantera/base/stringUtils.h"
using namespace std;
#ifndef FTN_TRAILING_UNDERSCORE
#define _DPOLFT_ dpolft
#define _DPCOEF_ dpcoef
#else
#define _DPOLFT_ dpolft_
#define _DPCOEF_ dpcoef_
#endif
extern "C" {
int _DPOLFT_(integer* n, doublereal* x, doublereal* y, doublereal* w,
integer* maxdeg, integer* ndeg, doublereal* eps, doublereal* r,
integer* ierr, doublereal* a);
int _DPCOEF_(integer* l, doublereal* c, doublereal* tc, doublereal* a);
}
namespace Cantera
{
@ -50,83 +32,4 @@ doublereal linearInterp(doublereal x, const vector_fp& xpts,
return ff;
}
//! Fits a polynomial function to a set of data points
/*!
* Given a collection of points X(I) and a set of values Y(I) which
* correspond to some function or measurement at each of the X(I),
* subroutine DPOLFT computes the weighted least-squares polynomial
* fits of all degrees up to some degree either specified by the user
* or determined by the routine. The fits thus obtained are in
* orthogonal polynomial form. Subroutine DP1VLU may then be
* called to evaluate the fitted polynomials and any of their
* derivatives at any point. The subroutine DPCOEF may be used to
* express the polynomial fits as powers of (X-C) for any specified
* point C.
*
* @param n The number of data points.
* @param x A set of grid points on which the data is specified.
* The array of values of the independent variable. These
* values may appear in any order and need not all be
* distinct. There are n of them.
* @param y array of corresponding function values. There are n of them
* @param w array of positive values to be used as weights. If
* W[0] is negative, DPOLFT will set all the weights
* to 1.0, which means unweighted least squares error
* will be minimized. To minimize relative error, the
* user should set the weights to: W(I) = 1.0/Y(I)**2,
* I = 1,...,N .
* @param maxdeg maximum degree to be allowed for polynomial fit.
* MAXDEG may be any non-negative integer less than N.
* Note -- MAXDEG cannot be equal to N-1 when a
* statistical test is to be used for degree selection,
* i.e., when input value of EPS is negative.
* @param ndeg output degree of the fit computed.
* @param eps Specifies the criterion to be used in determining
* the degree of fit to be computed.
* (1) If EPS is input negative, DPOLFT chooses the
* degree based on a statistical F test of
* significance. One of three possible
* significance levels will be used: .01, .05 or
* .10. If EPS=-1.0 , the routine will
* automatically select one of these levels based
* on the number of data points and the maximum
* degree to be considered. If EPS is input as
* -.01, -.05, or -.10, a significance level of
* .01, .05, or .10, respectively, will be used.
* (2) If EPS is set to 0., DPOLFT computes the
* polynomials of degrees 0 through MAXDEG .
* (3) If EPS is input positive, EPS is the RMS
* error tolerance which must be satisfied by the
* fitted polynomial. DPOLFT will increase the
* degree of fit until this criterion is met or
* until the maximum degree is reached.
* @param r Output vector containing the first LL+1 Taylor coefficients
* where LL=ABS(ndeg).
* P(X) = r[0] + r[1]*(X-C) + ... + r[ndeg] * (X-C)**ndeg
* ( here C = 0.0)
* @return returns the RMS error of the polynomial of degree ndeg .
*/
doublereal polyfit(int n, doublereal* x, doublereal* y, doublereal* w, int maxdeg, int& ndeg, doublereal eps, doublereal* r)
{
integer nn = n;
integer mdeg = maxdeg;
integer ndg = ndeg;
doublereal epss = eps;
integer ierr;
int worksize = 3*n + 3*maxdeg + 3;
vector_fp awork(worksize,0.0);
vector_fp coeffs(n+1, 0.0);
doublereal zer = 0.0;
_DPOLFT_(&nn, x, y, w, &mdeg, &ndg, &epss, &coeffs[0],
&ierr, &awork[0]);
if (ierr != 1) throw CanteraError("polyfit",
"DPOLFT returned error code IERR = " + int2str(ierr) +
"while attempting to fit " + int2str(n) + " data points "
+ "to a polynomial of degree " + int2str(maxdeg));
ndeg = ndg;
_DPCOEF_(&ndg, &zer, r, &awork[0]);
return epss;
}
}

52
src/numerics/polyfit.cpp Normal file
View file

@ -0,0 +1,52 @@
//! @file polyfit.cpp
#include "cantera/numerics/polyfit.h"
#include "cantera/base/ctexceptions.h"
#include "cantera/base/stringUtils.h"
#ifndef FTN_TRAILING_UNDERSCORE
#define _DPOLFT_ dpolft
#define _DPCOEF_ dpcoef
#else
#define _DPOLFT_ dpolft_
#define _DPCOEF_ dpcoef_
#endif
extern "C" {
int _DPOLFT_(integer* n, doublereal* x, doublereal* y, doublereal* w,
integer* maxdeg, integer* ndeg, doublereal* eps, doublereal* r,
integer* ierr, doublereal* a);
int _DPCOEF_(integer* l, doublereal* c, doublereal* tc, doublereal* a);
}
namespace Cantera
{
doublereal polyfit(int n, doublereal* x, doublereal* y, doublereal* w,
int maxdeg, int& ndeg, doublereal eps, doublereal* r)
{
integer nn = n;
integer mdeg = maxdeg;
integer ndg = ndeg;
doublereal epss = eps;
integer ierr;
int worksize = 3*n + 3*maxdeg + 3;
vector_fp awork(worksize,0.0);
vector_fp coeffs(n+1, 0.0);
doublereal zer = 0.0;
_DPOLFT_(&nn, x, y, w, &mdeg, &ndg, &epss, &coeffs[0],
&ierr, &awork[0]);
if (ierr != 1) {
throw CanteraError("polyfit",
"DPOLFT returned error code IERR = " + int2str(ierr) +
"while attempting to fit " + int2str(n) + " data points "
+ "to a polynomial of degree " + int2str(maxdeg));
}
ndeg = ndg;
_DPCOEF_(&ndg, &zer, r, &awork[0]);
return epss;
}
}