前言
入职第二周的任务是将导师的Python代码C化,发现Python中存在Numpy包直接调用np.polyfit就好了,但是C++不存在需要造轮子。
#include <iostream>
#include <cmath>
#include <vector>
#include <Eigen/QR>
#include "xtensor/xarray.hpp"
void polyfit( const std::vector<double> &t,const std::vector<double> &v,std::vector<double> &coeff,int order)
{// Create Matrix Placeholder of size n x k, n= number of datapoints, k = order of polynomial, for exame k = 3 for cubic polynomialEigen::MatrixXd T(t.size(), order + 1);Eigen::VectorXd V = Eigen::VectorXd::Map(&v.front(), v.size());//std::cout<<"ceshi"<<std::endl;//std::cout<<V<<std::endl;Eigen::VectorXd result;// check to make sure inputs are correctassert(t.size() == v.size());assert(t.size() >= order + 1);// Populate the matrixfor(size_t i = 0 ; i < t.size(); ++i){for(size_t j = 0; j < order + 1; ++j){T(i, j) = pow(t.at(i), j);}}std::cout<<T<<std::endl;// Solve for linear least square fitresult = T.householderQr().solve(V);coeff.resize(order+1);for (int k = 0; k < order+1; k++){coeff[k] = result[k];}}int main()
{// time valuestd::vector<double> time = {-2, 4, 6, 7, 9};std::vector<double> velocity = {5, 17, 37, 49, 82};// placeholder for storing polynomial coefficientstd::vector<double> coeff ;polyfit(time, velocity, coeff, 2);xt::xarray<double> c = xt::zeros<double>({3});for(int i = 0; i < coeff.size(); i++){c[i] = coeff[i];}std::vector<double> fitted_velocity;std::cout<< "Printing fitted values" << std::endl;for(int p = 0; p < time.size(); ++ p){double vfitted = coeff[0] + coeff[1]*time.at(p) + coeff[2]*(pow(time.at(p), 2)) ;std::cout<< vfitted<<", ";fitted_velocity.push_back(vfitted);}std::cout<<std::endl;for(int i = 0; i < c.size(); i++){std::cout<<c[i]<<std::endl;}std::cout<<std::endl;return 0;
}