diff --git a/include/cantera/numerics/eigen_dense.h b/include/cantera/numerics/eigen_dense.h index c925d433b..17de99157 100644 --- a/include/cantera/numerics/eigen_dense.h +++ b/include/cantera/numerics/eigen_dense.h @@ -8,4 +8,5 @@ namespace Cantera { typedef Eigen::Map MappedMatrix; typedef Eigen::Map MappedVector; + typedef Eigen::Map ConstMappedVector; } diff --git a/include/cantera/numerics/polyfit.h b/include/cantera/numerics/polyfit.h index c7cf472b2..588342362 100644 --- a/include/cantera/numerics/polyfit.h +++ b/include/cantera/numerics/polyfit.h @@ -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 diff --git a/src/numerics/polyfit.cpp b/src/numerics/polyfit.cpp index 9435a8298..c1bd0b5db 100644 --- a/src/numerics/polyfit.cpp +++ b/src/numerics/polyfit.cpp @@ -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); } } diff --git a/src/transport/GasTransport.cpp b/src/transport/GasTransport.cpp index 7aa22bb26..4ab0193ba 100644 --- a/src/transport/GasTransport.cpp +++ b/src/transport/GasTransport.cpp @@ -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; diff --git a/src/transport/MMCollisionInt.cpp b/src/transport/MMCollisionInt.cpp index b3f275c2d..ce9b79428 100644 --- a/src/transport/MMCollisionInt.cpp +++ b/src/transport/MMCollisionInt.cpp @@ -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); diff --git a/test/general/test_numerics.cpp b/test/general/test_numerics.cpp index c23713d1f..9eea78a33 100644 --- a/test/general/test_numerics.cpp +++ b/test/general/test_numerics.cpp @@ -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); }