#ifndef INCLUDE_POLYFIT_H
#define INCLUDE_POLYFIT_H

#include <vector>
#include <string>
#include "common.h"

namespace SIGMOD {



/** Polynomial fitting using discrete orthogonal polynomials to diagonalize the normal equations
 *  Input data is given as a vector<DataPoints>, to be evaluated either at indeces in vector
 *  or at values determined by the user. 
 *  The recursive polynomials are represented as vector<FourierPolynomial> of size(order) and is defined
 *  with the three-step recursion:
 *  
 *  p_0(x) = 1
 *  p_1(x) = (x-a0) 
 *  p_k+1(x) = (x-ak) * p_k(x) - bk  * p_k-1(x)
 * 
 *  a_k = sum( x(i) p_k(x)^2 ) / sum(p_k(x)^2)  k >= 0
 *
 *  b_k = sum(p_k(x)^2) / sum(p_k-1(x)^2) k > 0
 *
 *  The final least mean square polynomial is then
 *
 *  pol(k) = sum( t_k * p_k(x) )
 * 
 *  where the Fourier coefficients t_k are given by
 *
 *  t_k = sum (f(x) * p_k(x)^2 / sum(p_k-1(x)^2)
 *
 */

struct FourierPolynomial {
  double alfa;
  double beta;
};

class PolyFit {
private: 
  int m_polorder;
  std::vector<FourierPolynomial> m_pol; // Orthogonal polynomial
  std::vector<double> m_coeff; // and its fourier coeff
  std::vector<DataPoint> m_data; // raw input data;
  std::vector<DataPoint> m_processedData; // contains data with polynomial fitted and removed
  std::vector<DataPoint> m_polyData; // contains polynomial fit
  std::string m_input_filename;
  bool m_remove_pol_from_data;
  double m_stdv; // stdv of signal with polynomial removed
  ReadData m_access; // read data;
public:
  PolyFit();
  PolyFit(int order);
  PolyFit(int order, std::string filename, std::vector<int> columns);
  PolyFit(int order, std::vector<double>& data);
  PolyFit(int order, std::vector<DataPoint>& data);
  void transformData(std::vector<double>&);
  void doFit(int order, std::vector<double>& data); // interface to fitting
  void doFit(int order, std::vector<DataPoint>& data); // interface to fitting with user defined indep. values
  double evalPol(int ind); // evaluate the final interpolation polynomial at index ind
  void removePolFromData();
//   bool readData(std::string filename, std::vector<int> cols); // either 1 or columns
  void writeData();

  std::vector<DataPoint> data() const {return m_data;}
  std::vector<DataPoint>& data() {return m_data;}

  std::vector<DataPoint> processedData() const {return m_processedData;}
  std::vector<DataPoint>& processedData() {return m_processedData;}

  std::vector<DataPoint> polyData() const {return m_polyData;}
  std::vector<DataPoint>& polyData() {return m_polyData;}

  double stdv() const {return m_stdv;} // standard deviation of fitted data
private: 
  double P(double val, int order); // Orthogonal Polynom P(order) in expansion t_k * P_k(x)
  ReadData access() const {return m_access;}
};


} // namespace


#endif