#include "libscl.h"
#include "spline_interpolator.h"

using namespace std;
using namespace scl;

const bool debug = true;

int main(int argc, char** argp, char** envp)
{
  const INTEGER n_obs = 50;

  vector<string> header;
  realmat Z;

  string filename = "Z_" + fmt('d',4,n_obs)('0') + ".txt";

  readtable(filename.c_str(),Z,header);

  realmat z1 = Z("",1);

  z1.sort();

  INTEGER n = z1.size();

  realmat p(n,1);

  for (INTEGER i=1; i<=n; i++) p[i] = REAL(i)/REAL(n);

  ofstream fout;
  filename = "Fz.txt";
  fout.open(filename.c_str());
  if (!fout) error("Error, cannot open " + filename);
  
  for (INTEGER i=1; i<=n; i++) {
    fout << fmt('f',20,16,z1[i]) << fmt('f',20,16,p[i]) << '\n';
  }

  fout.clear(); fout.close();

  filename = "Fz_knots.txt";
  fout.open(filename.c_str());
  if (!fout) error("Error, cannot open " + filename);

  fout << fmt('f',20,16,z1[1]) << fmt('f',20,16,p[1]) << '\n';
  fout << fmt('f',20,16,z1[5]) << fmt('f',20,16,p[5]) << '\n';
  INTEGER lag = 5;
  for (INTEGER i=6; i<=135; i++) {
    if (p[i] > (p[lag] + 0.015) ) {
      fout << fmt('f',20,16,z1[i]) << fmt('f',20,16,p[i]) << '\n';
      lag = i;
    }
  }
  for (INTEGER i=136; i<=n-115; i++) {
    if (p[i] > (p[lag] + 0.05) ) {
      fout << fmt('f',20,16,z1[i]) << fmt('f',20,16,p[i]) << '\n';
      lag = i;
    }
  }
  for (INTEGER i=n-116; i<=n-15; i++) {
    if (p[i] > (p[lag] + 0.01) ) {
      fout << fmt('f',20,16,z1[i]) << fmt('f',20,16,p[i]) << '\n';
      lag = i;
    }
  }
  fout << fmt('f',20,16,z1[n-14]) << fmt('f',20,16,p[n-14]) << '\n';
  fout << fmt('f',20,16,z1[n-8]) << fmt('f',20,16,p[n-8]) << '\n';
  fout << fmt('f',20,16,z1[n-4]) << fmt('f',20,16,p[n-4]) << '\n';
  fout << fmt('f',20,16,z1[n-3]) << fmt('f',20,16,p[n-3]) << '\n';
  fout << fmt('f',20,16,z1[n-2]) << fmt('f',20,16,p[n-2]) << '\n';
  fout << fmt('f',20,16,z1[n-1]) << fmt('f',20,16,p[n-1]) << '\n';
  fout << fmt('f',20,16,z1[n]) << fmt('f',20,16,p[n]) << '\n';

  fout.clear(); fout.close();

  realmat Fz_knots;
  readtable(filename.c_str(),Fz_knots);

  filename = "Z_" + fmt('d',4,n_obs)('0') + "_knots.txt";
  writetable(filename.c_str(),Fz_knots,25,16);

  realmat z1_knots = Fz_knots("",1);
  realmat p_knots = Fz_knots("",2);

  spline_interpolator Fz_hat(z1_knots,p_knots);

  filename = "Fz_spline.txt";
  fout.open(filename.c_str());
  if (!fout) error("Error, cannot open " + filename);

  for (INTEGER i=1; i<=n; ++i) {
    fout << fmt('f',20,16,z1[i]) << fmt('f',20,16,Fz_hat(z1[i])) << '\n';
  }

  fout.clear(); fout.close();

  filename = "normal_knots.txt";
  realmat normal_knots;
  readtable(filename.c_str(),normal_knots);

  realmat xn_knots = normal_knots("",1);
  realmat pn_knots = normal_knots("",2);

  spline_interpolator qnorm(pn_knots,xn_knots);

  filename = "Z_" + fmt('d',4,n_obs)('0') + "_trf.txt";

  realmat Z_trf = Z;
  for (INTEGER i=1; i<=n; ++i) Z_trf(i,1) = qnorm(Fz_hat(Z(i,1)));
 
  writetable(filename.c_str(),Z_trf,header,20,8);

  return 0;
}

