Cantera  2.5.1
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
4 // at https://cantera.org/license.txt for license and copyright information.
5 
7 #include "cantera/numerics/eigen_dense.h"
8 #include "cantera/base/global.h"
10 
11 namespace Cantera
12 {
13 
14 double polyfit(size_t n, size_t deg, const double* xp, const double* yp,
15  const double* wp, double* pp)
16 {
17  ConstMappedVector x(xp, n);
18  Eigen::VectorXd y = ConstMappedVector(yp, n);
19  MappedVector p(pp, deg+1);
20 
21  if (deg >= n) {
22  throw CanteraError("polyfit", "Polynomial degree ({}) must be less "
23  "than number of input data points ({})", deg, n);
24  }
25 
26  // Construct A such that each row i of A has the elements
27  // 1, x[i], x[i]^2, x[i]^3 ... + x[i]^deg
28  Eigen::MatrixXd A(n, deg+1);
29  A.col(0).setConstant(1.0);
30 
31  if (deg > 0) {
32  A.col(1) = x;
33  }
34  for (size_t i = 1; i < deg; i++) {
35  A.col(i+1) = A.col(i).array() * x.array();
36  }
37 
38  if (wp != nullptr && wp[0] > 0) {
39  // For compatibility with old Fortran dpolft, input weights are the
40  // squares of the weight vector used in this algorithm
41  Eigen::VectorXd w = ConstMappedVector(wp, n).cwiseSqrt().eval();
42 
43  // Multiply by the weights on both sides
44  A = w.asDiagonal() * A;
45  y.array() *= w.array();
46  }
47 
48  // Solve W*A*p = W*y to find the polynomial coefficients
49  p = A.colPivHouseholderQr().solve(y);
50 
51  // Evaluate the computed polynomial at the input x coordinates to compute
52  // the RMS error as the return value
53  return (A*p - y).eval().norm() / sqrt(n);
54 }
55 
56 }
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:61
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...
This file contains definitions for utility functions and text for modules, inputfiles,...
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:264
double polyfit(size_t n, size_t deg, const double *xp, const double *yp, const double *wp, double *pp)
Fits a polynomial function to a set of data points.
Definition: polyfit.cpp:14