#include "strutil.h"
#include "common.h"
#include "polyfit.h"
#include "math.h"
#include <fstream>
#include <iostream>


namespace SIGMOD { // namespace SIGMOD
  using namespace std;

PolyFit::PolyFit() : m_remove_pol_from_data(false) {
}

PolyFit::PolyFit(int order) : m_polorder(order), m_remove_pol_from_data(false) {}
  
PolyFit::PolyFit(int order, vector<double>& data) : 
  m_remove_pol_from_data(false),
  m_polorder(order) {
  transformData(data);
  doFit(order,m_data);
}


PolyFit::PolyFit(int order, string filename, vector<int> columns) : 
  m_input_filename(filename),
  m_remove_pol_from_data(false),
  m_polorder(order) {
//   bool ok = readData(filename,columns);
  bool ok = access().readData(filename,columns,m_data);
  m_processedData = m_data;
  doFit(order,m_data);
}



void PolyFit::writeData() {
  int i;

  string base, extension;
  Strutil::fileBaseExtension(m_input_filename,base,extension);
  string gnuplot_file = base + string("_plot.plt");
  ofstream gnu(gnuplot_file.c_str(),ios::out);

  string data_file = base + "_fit.txt";
  std::ofstream out(data_file.c_str(),ios::out);


  cerr << "Writing gnuplot file "<<gnuplot_file<<" for visualizing data"<<endl;  
  gnu << "plot '"<<data_file<<"' using 1:2 title 'data' with lines,";
  gnu << "'"<<data_file<<"' using 1:3 title 'fitted polynomial' with lines\n"<<endl;
  gnu << "pause -1"<<endl;
  gnu << "plot '"<<data_file<<"' using 1:4 title 'data with pol removed' with lines\n"<<endl;
  gnu << "pause -1"<<endl;
  gnu.close();

  cerr << "Writing data, fitted polynomial and difference "<<endl;

  int nofp = data().size();  
  for (i=0;i<nofp;i++) {
    float abscissa = data()[i].xvalue;
    float datap =  data()[i].yvalue;
    float fit = polyData()[i].yvalue;
    float pol_removed = fit - datap;
    out << abscissa  << "  " 
        << datap  << "  " 
        << fit << "  "
        << pol_removed
        << endl;
  }
  out.close();
  

}


void PolyFit::transformData(std::vector<double>& vdata) {
  int i;
  m_data.reserve(vdata.size());
  for (i=0;i<vdata.size();i++)
    m_data.push_back(DataPoint(i,vdata[i],true));
}


double PolyFit::P(double xval, int order) {
  int j;
  double p_k_2, p_k_1, p_k_0;
  double x = xval;
  double retval;

  for (j=0;j<=order;j++) {
    if (j==0)
      p_k_0 = 1;
    else if (j==1)
      p_k_1 = (xval-m_pol[j-1].alfa);
    else {
      p_k_2 = (xval-m_pol[j-1].alfa) * p_k_1 - m_pol[j-1].beta * p_k_0;
      p_k_0 = p_k_1;
      p_k_1 = p_k_2;
    }
  }
  
  if (order == 0)
    retval = p_k_0;
  else if (order == 1)
    retval = p_k_1;
  else
    retval = p_k_2;

  return retval;
}


/** Calculate coeffs a_k and b_k in orthog. polynomial and fourier coefficients t_k
 *  and add these to global vectors. 
 *  
 *
 */
void PolyFit::doFit(int order, vector<DataPoint>& data) {
  int p,j;
  double nomsum_alfa, denomsum_alfa;
  double nomsum_beta, denomsum_beta;
  double nomsum_fourier, denomsum_fourier;
  double polval, polval2, polval_1, polval_1_2;

  // calculate alfa and beta coeffs and store in m_pol data struct. 

  cerr << "Evaluating Orthogonal polynomial..."<<endl;

  for (p=0;p<=m_polorder;p++) {
    nomsum_alfa = denomsum_alfa = 0.0;
    FourierPolynomial pol;
    if (p==0) {
      for (j=0;j<data.size();j++)
        nomsum_alfa += data[j].xvalue;
      pol.alfa = nomsum_alfa / data.size();
      pol.beta = 0;
      m_pol.push_back(pol);
    } else {
      nomsum_alfa = nomsum_beta = denomsum_alfa = denomsum_beta = 0;
      for (j=0;j<data.size();j++) {
        double x =  data[j].xvalue;
        polval = P(x,p);
        polval2 = polval*polval;
        polval_1 = P(x,p-1);
        polval_1_2 = polval_1*polval_1;
        nomsum_alfa += x*polval2;
        denomsum_alfa += polval2;
        nomsum_beta += polval2;
        denomsum_beta += polval_1_2;
      }
      pol.alfa = nomsum_alfa/denomsum_alfa;
      pol.beta = nomsum_beta/denomsum_beta;
      m_pol.push_back(pol);
    }
  }

//   for (p=0;p<m_polorder;p++) {
//     cerr << "alfa["<<p<<"] = "<<m_pol[p].alfa<<endl;
//     cerr << "beta["<<p<<"] = "<<m_pol[p].beta<<endl;
//   }

  // Fourier coefficients

  cerr << "Calcualting Fourier coefficients..."<<endl;
  for (p=0;p<=m_polorder;p++) {
    nomsum_fourier = denomsum_fourier = 0;
    for (j=0;j<data.size();j++) {
      double y =  data[j].yvalue;
      double x =  data[j].xvalue;
      polval = P(x,p);
      polval2 = polval*polval;
      nomsum_fourier += y*polval;
      denomsum_fourier += polval2;
    }
    double tk = nomsum_fourier / denomsum_fourier;
    m_coeff.push_back(tk);
    cerr << "fourier coeff["<<p<<"] = "<<m_coeff[p]<<endl;
  }

  // finally evaluate the polynomial in the data range and store

  cerr << "Evaluating polynomial..."<<endl;
  for (j=0;j<data.size();j++) {  
    cerr << "\r processing data point "<<j;
    double yval = evalPol(j);
    double xval = m_data[j].xvalue;
    DataPoint dataP(xval,yval,true);
    m_polyData.push_back(dataP);
  }
  cerr << endl;
}


double PolyFit::evalPol(int ind) {
  int p;
  double sum = 0;
  for (p=0;p<=m_polorder;p++) {
    double x = m_data[ind].xvalue;
    sum += m_coeff[p] * P(x,p);
  }
  return sum;
}

void PolyFit::removePolFromData() {
  m_remove_pol_from_data = true;
  int j;
  double sum, sum2;
  sum = sum2 = 0;
  int nofp = m_data.size();
  for (j=0;j<nofp;j++) {
    DataPoint rempol(m_data[j].xvalue,evalPol(j) - m_data[j].yvalue,true);
    sum += rempol.yvalue;
    sum2 += rempol.yvalue*rempol.yvalue;
    m_processedData.push_back(rempol);
  }
  double mean = sum/nofp;
  m_stdv = sqrt(sum2/nofp - mean*mean);
}  

} // namespace