[Numerics] Modify polyfit to use Eigen instead of fortran/f2c code
This commit is contained in:
Ray Speth
parent
230688d31d
commit
efb068e319
6 changed files with 92 additions and 110 deletions
|
|
@ -8,4 +8,5 @@
|
|||
namespace Cantera {
|
||||
typedef Eigen::Map<Eigen::MatrixXd> MappedMatrix;
|
||||
typedef Eigen::Map<Eigen::VectorXd> MappedVector;
|
||||
typedef Eigen::Map<const Eigen::VectorXd> ConstMappedVector;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,6 +1,4 @@
|
|||
/**
|
||||
* @file polyfit.h C interface for Fortran DPOLFT subroutine
|
||||
*/
|
||||
//! @file polyfit.h
|
||||
/*
|
||||
* Copyright 2001-2003 California Institute of Technology
|
||||
* See file License.txt for licensing information
|
||||
|
|
@ -16,55 +14,32 @@ namespace Cantera
|
|||
|
||||
//! 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.
|
||||
* Given a collection of *n* points *x* and a set of values *y* of some function
|
||||
* evaluated at those points, this function computes the weighted least-squares
|
||||
* polynomial fit of degree *deg*:
|
||||
*
|
||||
* @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.
|
||||
* \f[ f(x) = p[0] + p[1]*x + p[2]*x^2 + \cdots + p[deg]*x^deg \f]
|
||||
*
|
||||
* @param r Output vector containing the first ndeg+1 Taylor coefficients
|
||||
*
|
||||
* P(X) = r[0] + r[1]*(X-C) + ... + r[ndeg] * (X-C)**ndeg
|
||||
* ( here C = 0.0)
|
||||
* @returns value of the rms of the interpolated function at x.
|
||||
* @param n The number of points at which the function is evaluated
|
||||
* @param deg The degree of the polynomial fit to be computed. deg <= n - 1.
|
||||
* @param x Array of points at which the function is evaluated. Length *n*.
|
||||
* @param y Array of function values at the points in *x*. Length *n*.
|
||||
* @param w Array of weights. If w == nullptr or w[0] < 0, then all the
|
||||
* weights will be set to 1.0.
|
||||
* @param[out] p Array of polynomial coefficients, starting with the constant
|
||||
* term. Length *deg+1*.
|
||||
* @returns the root mean squared error of the fit at the input points.
|
||||
*/
|
||||
double polyfit(size_t n, size_t deg, const double* x, const double* y,
|
||||
const double* w, double* p);
|
||||
|
||||
//! Fits a polynomial function to a set of data points
|
||||
/*!
|
||||
* @deprecated The ndeg and eps arguments to polyfit are deprecated and unused.
|
||||
* Use the form of polyfit with signature polyfit(n, deg, x, y, w, p). To be
|
||||
* removed after Cantera 2.3.
|
||||
*/
|
||||
doublereal polyfit(int n, doublereal* x, doublereal* y, doublereal* w,
|
||||
int maxdeg, int& ndeg, doublereal eps, doublereal* r);
|
||||
|
||||
}
|
||||
#endif
|
||||
|
|
|
|||
|
|
@ -1,51 +1,64 @@
|
|||
//! @file polyfit.cpp
|
||||
|
||||
#include "cantera/numerics/polyfit.h"
|
||||
#include "cantera/numerics/eigen_dense.h"
|
||||
#include "cantera/base/global.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)
|
||||
double polyfit(int n, double* xp, double* yp, double* wp,
|
||||
int deg, int& ndeg, double eps, double* rp)
|
||||
{
|
||||
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;
|
||||
warn_deprecated("polyfit(n, x, y, w, maxdeg, ndeg, eps, r",
|
||||
"The ndeg and eps arguments to polyfit are deprecated and unused. Use "
|
||||
"the form of polyfit with signature polyfit(n, deg, x, y, w, p). To be "
|
||||
"removed after Cantera 2.3.");
|
||||
ndeg = deg;
|
||||
return polyfit(n, deg, xp, yp, wp, rp);
|
||||
}
|
||||
|
||||
_DPOLFT_(&nn, x, y, w, &mdeg, &ndg, &epss, &coeffs[0],
|
||||
&ierr, &awork[0]);
|
||||
if (ierr != 1) {
|
||||
throw CanteraError("polyfit",
|
||||
"DPOLFT returned error code IERR = {} while attempting to fit {}"
|
||||
" data points to a polynomial of degree {}", ierr, n, maxdeg);
|
||||
double polyfit(size_t n, size_t deg, const double* xp, const double* yp,
|
||||
const double* wp, double* pp)
|
||||
{
|
||||
ConstMappedVector x(xp, n);
|
||||
Eigen::VectorXd y = ConstMappedVector(yp, n);
|
||||
MappedVector p(pp, deg+1);
|
||||
|
||||
if (deg >= n) {
|
||||
throw CanteraError("polyfit", "Polynomial degree ({}) must be less "
|
||||
"than number of input data points ({})", deg, n);
|
||||
}
|
||||
ndeg = ndg;
|
||||
_DPCOEF_(&ndg, &zer, r, &awork[0]);
|
||||
return epss;
|
||||
|
||||
// Construct A such that each row i of A has the elements
|
||||
// 1, x[i], x[i]^2, x[i]^3 ... + x[i]^deg
|
||||
Eigen::MatrixXd A(n, deg+1);
|
||||
A.col(0).setConstant(1.0);
|
||||
|
||||
if (deg > 0) {
|
||||
A.col(1) = x;
|
||||
}
|
||||
for (size_t i = 1; i < deg; i++) {
|
||||
A.col(i+1) = A.col(i).array() * x.array();
|
||||
}
|
||||
|
||||
if (wp != nullptr && wp[0] > 0) {
|
||||
// For compatibility with old Fortran dpolft, input weights are the
|
||||
// squares of the weight vector used in this algorithm
|
||||
Eigen::VectorXd w = ConstMappedVector(wp, n).cwiseSqrt().eval();
|
||||
|
||||
// Multiply by the weights on both sides
|
||||
A = w.asDiagonal() * A;
|
||||
y.array() *= w.array();
|
||||
}
|
||||
|
||||
// Solve W*A*p = W*y to find the polynomial coefficients
|
||||
p = A.colPivHouseholderQr().solve(y);
|
||||
|
||||
// Evaluate the computed polynomial at the input x coordinates to compute
|
||||
// the RMS error as the return value
|
||||
return (A*p - y).eval().norm() / sqrt(n);
|
||||
}
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -557,7 +557,6 @@ void GasTransport::fitCollisionIntegrals(MMCollisionInt& integrals)
|
|||
|
||||
void GasTransport::fitProperties(MMCollisionInt& integrals)
|
||||
{
|
||||
int ndeg = 0;
|
||||
// number of points to use in generating fit data
|
||||
const size_t np = 50;
|
||||
int degree = (m_mode == CK_Mode ? 3 : 4);
|
||||
|
|
@ -651,10 +650,8 @@ void GasTransport::fitProperties(MMCollisionInt& integrals)
|
|||
w2[n] = 1.0/(spcond[n]*spcond[n]);
|
||||
}
|
||||
}
|
||||
polyfit(np, tlog.data(), spvisc.data(),
|
||||
w.data(), degree, ndeg, 0.0, c.data());
|
||||
polyfit(np, tlog.data(), spcond.data(),
|
||||
w.data(), degree, ndeg, 0.0, c2.data());
|
||||
polyfit(np, degree, tlog.data(), spvisc.data(), w.data(), c.data());
|
||||
polyfit(np, degree, tlog.data(), spcond.data(), w.data(), c2.data());
|
||||
|
||||
// evaluate max fit errors for viscosity
|
||||
for (size_t n = 0; n < np; n++) {
|
||||
|
|
@ -752,8 +749,7 @@ void GasTransport::fitProperties(MMCollisionInt& integrals)
|
|||
w[n] = 1.0/(diff[n]*diff[n]);
|
||||
}
|
||||
}
|
||||
polyfit(np, tlog.data(), diff.data(),
|
||||
w.data(), degree, ndeg, 0.0, c.data());
|
||||
polyfit(np, degree, tlog.data(), diff.data(), w.data(), c.data());
|
||||
|
||||
for (size_t n = 0; n < np; n++) {
|
||||
double val, fit;
|
||||
|
|
|
|||
|
|
@ -296,7 +296,6 @@ doublereal MMCollisionInt::fitDelta(int table, int ntstar, int degree, doublerea
|
|||
{
|
||||
vector_fp w(8);
|
||||
doublereal* begin = 0;
|
||||
int ndeg=0;
|
||||
switch (table) {
|
||||
case 0:
|
||||
begin = omega22_table + 8*ntstar;
|
||||
|
|
@ -314,7 +313,7 @@ doublereal MMCollisionInt::fitDelta(int table, int ntstar, int degree, doublerea
|
|||
return 0.0;
|
||||
}
|
||||
w[0] = -1.0;
|
||||
return polyfit(8, delta, begin, w.data(), degree, ndeg, 0.0, c);
|
||||
return polyfit(8, degree, delta, begin, w.data(), c);
|
||||
}
|
||||
|
||||
doublereal MMCollisionInt::omega22(double ts, double deltastar)
|
||||
|
|
@ -417,7 +416,6 @@ void MMCollisionInt::fit_omega22(int degree, doublereal deltastar,
|
|||
doublereal* o22)
|
||||
{
|
||||
int i, n = m_nmax - m_nmin + 1;
|
||||
int ndeg=0;
|
||||
vector_fp values(n);
|
||||
doublereal rmserr;
|
||||
vector_fp w(n);
|
||||
|
|
@ -430,7 +428,7 @@ void MMCollisionInt::fit_omega22(int degree, doublereal deltastar,
|
|||
}
|
||||
}
|
||||
w[0]= -1.0;
|
||||
rmserr = polyfit(n, logT, values.data(), w.data(), degree, ndeg, 0.0, o22);
|
||||
rmserr = polyfit(n, degree, logT, values.data(), w.data(), o22);
|
||||
if (m_loglevel > 0 && rmserr > 0.01) {
|
||||
writelogf("Warning: RMS error = %12.6g in omega_22 fit"
|
||||
"with delta* = %12.6g\n", rmserr, deltastar);
|
||||
|
|
@ -441,7 +439,6 @@ void MMCollisionInt::fit(int degree, doublereal deltastar,
|
|||
doublereal* a, doublereal* b, doublereal* c)
|
||||
{
|
||||
int i, n = m_nmax - m_nmin + 1;
|
||||
int ndeg=0;
|
||||
vector_fp values(n);
|
||||
doublereal rmserr;
|
||||
vector_fp w(n);
|
||||
|
|
@ -454,7 +451,7 @@ void MMCollisionInt::fit(int degree, doublereal deltastar,
|
|||
}
|
||||
}
|
||||
w[0]= -1.0;
|
||||
rmserr = polyfit(n, logT, values.data(), w.data(), degree, ndeg, 0.0, a);
|
||||
rmserr = polyfit(n, degree, logT, values.data(), w.data(), a);
|
||||
|
||||
for (i = 0; i < n; i++) {
|
||||
if (deltastar == 0.0) {
|
||||
|
|
@ -464,7 +461,7 @@ void MMCollisionInt::fit(int degree, doublereal deltastar,
|
|||
}
|
||||
}
|
||||
w[0]= -1.0;
|
||||
rmserr = polyfit(n, logT, values.data(), w.data(), degree, ndeg, 0.0, b);
|
||||
rmserr = polyfit(n, degree, logT, values.data(), w.data(), b);
|
||||
|
||||
for (i = 0; i < n; i++) {
|
||||
if (deltastar == 0.0) {
|
||||
|
|
@ -474,7 +471,7 @@ void MMCollisionInt::fit(int degree, doublereal deltastar,
|
|||
}
|
||||
}
|
||||
w[0]= -1.0;
|
||||
rmserr = polyfit(n, logT, values.data(), w.data(), degree, ndeg, 0.0, c);
|
||||
rmserr = polyfit(n, degree, logT, values.data(), w.data(), c);
|
||||
if (m_loglevel > 2) {
|
||||
writelogf("\nT* fit at delta* = %.6g\n", deltastar);
|
||||
|
||||
|
|
|
|||
|
|
@ -17,11 +17,10 @@ TEST(Polyfit, exact_fit)
|
|||
{
|
||||
vector_fp x{0, 0.3, 1.0, 1.5, 2.0, 2.5};
|
||||
vector_fp p(6);
|
||||
vector_fp w(6, 1.0);
|
||||
vector_fp w(6, -1.0);
|
||||
for (int i = 0; i < 20; i++) {
|
||||
vector_fp y{-1.1*i, cos(i), pow(-1,i), 3.2/(i+1), 0.1*i*i, sin(i)};
|
||||
int ndeg;
|
||||
polyfit(6, x.data(), y.data(), w.data(), 5, ndeg, 0, p.data());
|
||||
polyfit(6, 5, x.data(), y.data(), w.data(), p.data());
|
||||
for (size_t j = 0; j < 6; j++) {
|
||||
EXPECT_NEAR(polyval(p, x[j]), y[j], 1e-12);
|
||||
}
|
||||
|
|
@ -47,13 +46,13 @@ TEST(Polyfit, sequential)
|
|||
0.011452361452361514, 0.10963690963690906, -0.022222222222222105}
|
||||
};
|
||||
|
||||
vector_fp w(7, 1.0);
|
||||
int ndeg;
|
||||
double rms_prev = 1e10;
|
||||
for (size_t i = 0; i < PP.size(); i++) {
|
||||
size_t N = i + 1;
|
||||
vector_fp p(N);
|
||||
polyfit(7, x.data(), y.data(), w.data(), i, ndeg, 0, p.data());
|
||||
ASSERT_EQ(ndeg, (int) i);
|
||||
double rms = polyfit(7, i, x.data(), y.data(), nullptr, p.data());
|
||||
EXPECT_LT(rms, rms_prev);
|
||||
rms_prev = rms;
|
||||
for (size_t j = 0; j < N; j++) {
|
||||
EXPECT_NEAR(PP[i][j], p[j], 1e-14);
|
||||
}
|
||||
|
|
@ -80,12 +79,13 @@ TEST(Polyfit, weighted)
|
|||
0.011482911646053995, 0.10962944760868476, -0.022222284629403764}
|
||||
};
|
||||
|
||||
int ndeg;
|
||||
double rms_prev = 1e10;
|
||||
for (size_t i = 0; i < PP.size(); i++) {
|
||||
size_t N = i + 1;
|
||||
vector_fp p(N);
|
||||
polyfit(7, x.data(), y.data(), w.data(), i, ndeg, 0, p.data());
|
||||
ASSERT_EQ(ndeg, (int) i);
|
||||
double rms = polyfit(7, i, x.data(), y.data(), w.data(), p.data());
|
||||
EXPECT_LT(rms, rms_prev);
|
||||
rms_prev = rms;
|
||||
for (size_t j = 0; j < N; j++) {
|
||||
EXPECT_NEAR(PP[i][j], p[j], 1e-14);
|
||||
}
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue