Cantera  2.3.0
polyfit.cpp
Go to the documentation of this file.
1 //! @file polyfit.cpp
2
3 // This file is part of Cantera. See License.txt in the top-level directory or
5
7 #include "cantera/numerics/eigen_dense.h"
8 #include "cantera/base/global.h"
10
11 namespace Cantera
12 {
13
14 double polyfit(int n, double* xp, double* yp, double* wp,
15  int deg, int& ndeg, double eps, double* rp)
16 {
17  warn_deprecated("polyfit(n, x, y, w, maxdeg, ndeg, eps, r",
18  "The ndeg and eps arguments to polyfit are deprecated and unused. Use "
19  "the form of polyfit with signature polyfit(n, deg, x, y, w, p). To be "
20  "removed after Cantera 2.3.");
21  ndeg = deg;
22  return polyfit(n, deg, xp, yp, wp, rp);
23 }
24
25 double polyfit(size_t n, size_t deg, const double* xp, const double* yp,
26  const double* wp, double* pp)
27 {
28  ConstMappedVector x(xp, n);
29  Eigen::VectorXd y = ConstMappedVector(yp, n);
30  MappedVector p(pp, deg+1);
31
32  if (deg >= n) {
33  throw CanteraError("polyfit", "Polynomial degree ({}) must be less "
34  "than number of input data points ({})", deg, n);
35  }
36
37  // Construct A such that each row i of A has the elements
38  // 1, x[i], x[i]^2, x[i]^3 ... + x[i]^deg
39  Eigen::MatrixXd A(n, deg+1);
40  A.col(0).setConstant(1.0);
41
42  if (deg > 0) {
43  A.col(1) = x;
44  }
45  for (size_t i = 1; i < deg; i++) {
46  A.col(i+1) = A.col(i).array() * x.array();
47  }
48
49  if (wp != nullptr && wp[0] > 0) {
50  // For compatibility with old Fortran dpolft, input weights are the
51  // squares of the weight vector used in this algorithm
52  Eigen::VectorXd w = ConstMappedVector(wp, n).cwiseSqrt().eval();
53
54  // Multiply by the weights on both sides
55  A = w.asDiagonal() * A;
56  y.array() *= w.array();
57  }
58
59  // Solve W*A*p = W*y to find the polynomial coefficients
60  p = A.colPivHouseholderQr().solve(y);
61
62  // Evaluate the computed polynomial at the input x coordinates to compute
63  // the RMS error as the return value
64  return (A*p - y).eval().norm() / sqrt(n);
65 }
66
67 }
void warn_deprecated(const std::string &method, const std::string &extra)
Print a warning indicating that method is deprecated.
Definition: global.cpp:54
This file contains definitions for utility functions and text for modules, inputfiles, logs, textlogs, (see Input File Handling, Diagnostic Output, and Writing messages to the screen).
double polyfit(int n, double *xp, double *yp, double *wp, int deg, int &ndeg, double eps, double *rp)
Fits a polynomial function to a set of data points.
Definition: polyfit.cpp:14
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
Namespace for the Cantera kernel.
Definition: application.cpp:29
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...