Author: Not specified Language: cpp
Description: Not specified Timestamp: 2018-07-06 09:00:55 +0000
View raw paste Reply
  1.  
  2. #include <Eigen/QR>
  3. #include <stdio.h>
  4. #include <vector>
  5. #include <cmath>
  6.  
  7. void polyfit(const std::vector<double> &xv, const std::vector<double> &yv, std::vector<double> &coeff, int order)
  8. {
  9.     Eigen::MatrixXd A(xv.size(), order+1);
  10.     Eigen::VectorXd yv_mapped = Eigen::VectorXd::Map(&yv.front(), yv.size());
  11.     Eigen::VectorXd result;
  12.  
  13.     assert(xv.size() == yv.size());
  14.     assert(xv.size() >= order+1);
  15.  
  16.  
  17.     for (size_t i = 0; i < xv.size(); i++)
  18.         for (size_t j = 0; j < order+1; j++)
  19.             A(i, j) = pow(xv.at(i), j);
  20.  
  21.  
  22.     result = A.householderQr().solve(yv_mapped);
  23.  
  24.     coeff.resize(order+1);
  25.     for (size_t i = 0; i < order+1; i++)
  26.         coeff[i] = result[i];
  27. }
  28.  
  29. int main()
  30. {
  31.     std::vector<double> x_values, y_values, coeff;
  32.     double x, y;
  33.  
  34.     while (scanf("%lf %lf\n", &x, &y) == 2) {
  35.         x_values.push_back(x);
  36.         y_values.push_back(y);
  37.     }
  38.  
  39.     polyfit(x_values, y_values, coeff, 3);
  40.     printf("%f + %f*x + %f*x^2 + %f*x^3\n", coeff[0], coeff[1], coeff[2], coeff[3]);
  41.  
  42.     return 0;
  43. }
  44.  
View raw paste Reply