#include "libscl.h"
#include "linear_interpolator.h"
#include "spline_interpolator.h"
using namespace scl;
using namespace std;



REAL h(REAL x) { return exp(x) + cos(x); }
REAL dh(REAL x) { return exp(x) + -sin(x); }

int main(int argc, char** argp, char** envp)
{
  INTEGER m = 6;
  INTEGER n = m - 1;
  REAL scale = 1.0/REAL(n);

  realmat x0(m,1);
  realmat y0(m,1);
  realmat dy0(m,1);

  for (INTEGER i=1; i<=m; ++i) {
    x0[i] = scale*REAL(i-1);
    y0[i] = h(x0[i]);
    dy0[i] = dh(x0[i]);
  }

  vector<string> header; 
  header.push_back("x"); 
  header.push_back("y");

  linear_interpolator f(x0,y0);
  spline_interpolator g(x0,y0);

  INTEGER N = 15*m;
  scale = scale*REAL(n)/REAL(N);
  realmat x(N,1);
  realmat y(N,1);
  realmat dy(N,1);

  for (INTEGER i=1; i<=N; ++i) {
    x[i] = scale*REAL(i);
  }

  for (INTEGER i=1; i<=N; ++i) {
    y[i] = f(x[i]);
  }

  cout << '\n';
  writetable(cout,cbind(x,y),header,20,8);

  for (INTEGER i=1; i<=N; ++i) {
    y[i] = f(x[i])/h(x[i]);
  }

  cout << '\n';
  writetable(cout,cbind(x,y),header,20,8);

  header[1] = "f(x)/h(x)";

  for (INTEGER i=1; i<=N; ++i) {
    y[i] = g(x[i]);
    dy[i] = g.derivative(x[i]);
  }

  header[1] = "y";

  header.push_back("dy");

  realmat xy = cbind(x,y);

  cout << '\n';
  writetable(cout,cbind(xy,dy),header,20,8);

  for (INTEGER i=1; i<=N; ++i) {
    y[i] = y[i]/h(x[i]);
    dy[i] = dy[i]/dh(x[i]);
  }

  xy = cbind(x,y);

  header[0] = "x";
  header[1] = "g(x)/h(x)";
  header[2] = "dg(x)/dh(x)";

  cout << '\n';
  writetable(cout,cbind(xy,dy),header,20,8);

  xy = cbind(x0,y0);

  header[0] = "x0";
  header[1] = "y0";
  header[2] = "dy0";

  cout << '\n';
  writetable(cout,cbind(xy,dy0),header,20,8);

  cout << '\n';
  cout << "knot intercept slope origin" << '\n';
  for (INTEGER i=0; i<f.size(); ++i) {
    cout << "\t" << i 
         << fmt('f',8,3,f[i].intercept())
         << fmt('f',8,3,f[i].slope())
         << fmt('f',8,3,f[i].origin())
         << '\n';
  }

  cout << '\n';
  cout << "knot intercept linear quaratic cubic origin" << '\n';
  for (INTEGER i=0; i<g.size(); ++i) {
    cout << "\t" << i 
         << fmt('f',8,3,g[i].intercept())
         << fmt('f',8,3,g[i].linear())
         << fmt('f',8,3,g[i].quadratic())
         << fmt('f',8,3,g[i].cubic())
         << fmt('f',8,3,g[i].origin())
         << '\n';
  }

  return 0;
}
