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

int main()
{
  //const REAL root_two_pi = 2.5066282746310005024;

  INTEGER dim_x = 3;
  INTEGER deg_main = 4;
  INTEGER deg_inter = 2;

  realmat x(dim_x,1);
  //realmat x(1,dim_x);

  for (INTEGER i=1; i<=dim_x; i++) x[i]=i;

  poly regular('r',x,deg_main,deg_inter);

  realmat basis;
  regular.get_basis(basis);

  std::vector<intvec> multi;
  regular.get_multi(multi);

  std::vector<std::string> smulti;
  regular.get_multi(smulti);

  cout << starbox("/x, multi-indexes, basis//");

  for (INTEGER j=1; j<=regular.get_len(); j++) {
    if (j < 10) cout << "                          Row  " << j << "  "; 
    else        cout << "                          Row " << j << "  ";
    for (INTEGER i=1; i<=dim_x; i++) cout << x[i];
    cout << "  ";
    for (INTEGER i=1; i<=dim_x; i++) cout << multi[j-1][i];
    cout << "  ";
    cout << smulti[j-1];
    cout << "  ";
    cout << basis[j];
    cout << '\n';
  }

  INTEGER len_basis = regular.get_len();
  cout << "                          "; 
  cout << "len = " << len_basis << '\n';

  realmat jacobian;
  realmat numerical(len_basis,dim_x);
  realmat basis_hi;
  realmat basis_lo;
  const REAL eps = pow(REAL_EPSILON,0.33333333);

  cout << starbox("/Jacobian, analytical, numerical, relative error//");

  REAL save = x[dim_x];
  x[dim_x] = 0;

  regular.set_x(x);
  regular.get_jacobian(jacobian); 

  for (INTEGER j=1; j<=dim_x; j++) {
    REAL tmp = x[j];
    REAL h = eps*std::fabs(tmp);
    if (h == 0) h = eps;
    REAL hi = tmp + h;
    REAL lo = tmp - h;
    x[j] = hi;
    regular.set_x(x);
    regular.get_basis(basis_hi);
    x[j] = lo;
    regular.set_x(x);
    regular.get_basis(basis_lo);
    REAL difference = hi - lo;
    x[j] = tmp;
    for (INTEGER i=1; i<=len_basis; i++) {
      numerical(i,j) = (basis_hi[i]-basis_lo[i])/difference;
    }
  }

  x[dim_x] = save;

  for (INTEGER i=1; i<=len_basis; i++) {
    for (INTEGER j=1; j<=dim_x; j++) {
      cout << "                           ";
      if (i < 10) cout << "Row  " << i << "   ";
      else        cout << "Row " << i << "   ";
      cout << "Elem " << j << "   "
           << jacobian(i,j) << "  "
	   << numerical(i,j) << "  "
	   << (numerical(i,j) - jacobian(i,j))/(fabs(jacobian(i,j))+eps)
	   << '\n';
    }
  }
  
  for (INTEGER i=1; i<=dim_x; i++) x[i]=i;

  poly hermite('h',x,deg_main,deg_inter);

  hermite.get_basis(basis);

  hermite.get_multi(multi);

  cout << starbox("/x, multi-indexes, basis//");

  for (INTEGER j=1; j<=hermite.get_len(); j++) {
    if (j < 10) cout << "                          Row  " << j << "  "; 
    else        cout << "                          Row " << j << "  ";
    for (INTEGER i=1; i<=dim_x; i++) cout << x[i];
    cout << "  ";
    for (INTEGER i=1; i<=dim_x; i++) cout << multi[j-1][i];
    cout << "  ";
    cout << basis[j];
    cout << '\n';
  }

  len_basis = hermite.get_len();
  cout << "                          "; 
  cout << "len = " << len_basis << '\n';
 
  cout << starbox("/Jacobian, analytical, numerical, relative error//");

  x[dim_x] = 0;
  hermite.set_x(x);
  hermite.get_jacobian(jacobian); 

  for (INTEGER j=1; j<=dim_x; j++) {
    REAL tmp = x[j];
    REAL h = eps*std::fabs(tmp);
    if (h == 0) h = eps;
    REAL hi = tmp + h;
    REAL lo = tmp - h;
    x[j] = hi;
    hermite.set_x(x);
    hermite.get_basis(basis_hi);
    x[j] = lo;
    hermite.set_x(x);
    hermite.get_basis(basis_lo);
    REAL difference = hi - lo;
    x[j] = tmp;
    for (INTEGER i=1; i<=len_basis; i++) {
      numerical(i,j) = (basis_hi[i]-basis_lo[i])/difference;
    }
  }

  for (INTEGER i=1; i<=len_basis; i++) {
    for (INTEGER j=1; j<=dim_x; j++) {
      cout << "                           ";
      if (i < 10) cout << "Row  " << i << "   ";
      else        cout << "Row " << i << "   ";
      cout << "Elem " << j << "   "
           << jacobian(i,j) << "  "
	   << numerical(i,j) << "  "
	   << (numerical(i,j) - jacobian(i,j))/(fabs(jacobian(i,j))+eps)
	   << '\n';
    }
  }

  cout << starbox("/Orthogonality check//");

  INT_32BIT seed = 15;
  INTEGER reps = 100000;
  realmat sum(len_basis,len_basis,0.0);

  for (INTEGER j=1; j<=reps; j++) {
    for (INTEGER i=1; i<=dim_x; i++) {
      x[i] = unsk(seed);
    }
    hermite.set_x(x);
    hermite.get_basis(basis);
    sum += basis*T(basis);
  }
  sum = sum/reps;

  cout << sum;
  
  for (INTEGER i=1; i<=dim_x; i++) x[i]=i;

  poly laguerre('l',x,deg_main,deg_inter);

  laguerre.get_basis(basis);

  laguerre.get_multi(multi);

  cout << starbox("/x, multi-indexes, basis//");

  for (INTEGER j=1; j<=laguerre.get_len(); j++) {
    if (j < 10) cout << "                          Row  " << j << "  "; 
    else        cout << "                          Row " << j << "  ";
    for (INTEGER i=1; i<=dim_x; i++) cout << x[i];
    cout << "  ";
    for (INTEGER i=1; i<=dim_x; i++) cout << multi[j-1][i];
    cout << "  ";
    cout << basis[j];
    cout << '\n';
  }

  len_basis = laguerre.get_len();
  cout << "                          "; 
  cout << "len = " << len_basis << '\n';
 
  cout << starbox("/Jacobian, analytical, numerical, relative error//");

  x[dim_x] = 0;
  laguerre.set_x(x);
  laguerre.get_jacobian(jacobian); 

  for (INTEGER j=1; j<=dim_x; j++) {
    REAL tmp = x[j];
    REAL h = eps*std::fabs(tmp);
    if (h == 0) h = eps;
    REAL hi = tmp + h;
    REAL lo = tmp - h;
    x[j] = hi;
    laguerre.set_x(x);
    laguerre.get_basis(basis_hi);
    x[j] = lo;
    laguerre.set_x(x);
    laguerre.get_basis(basis_lo);
    REAL difference = hi - lo;
    x[j] = tmp;
    for (INTEGER i=1; i<=len_basis; i++) {
      numerical(i,j) = (basis_hi[i]-basis_lo[i])/difference;
    }
  }

  for (INTEGER i=1; i<=len_basis; i++) {
    for (INTEGER j=1; j<=dim_x; j++) {
      cout << "                           ";
      if (i < 10) cout << "Row  " << i << "   ";
      else        cout << "Row " << i << "   ";
      cout << "Elem " << j << "   "
           << jacobian(i,j) << "  "
	   << numerical(i,j) << "  "
	   << (numerical(i,j) - jacobian(i,j))/(fabs(jacobian(i,j))+eps)
	   << '\n';
    }
  }

  cout << starbox("/Orthogonality check//");

  fill(sum,0.0);

  for (INTEGER j=1; j<=reps; j++) {
    for (INTEGER i=1; i<=dim_x; i++) {
      REAL u = ran(seed);
      if (u < 1.0) {
        x[i] = -log(1.0 - u);
      }
      else {
	error("Error, tstploy, this should not happen");
      }
    }
    laguerre.set_x(x);
    laguerre.get_basis(basis);
    sum += basis*T(basis);
  }
  sum = sum/reps;

  cout << sum;

  cout << starbox("/Check default constructor//") << '\n';
  poly null;
  cout << "len_basis = " << null.get_len() << '\n';
  
  cout << starbox("/Check copy constructor//") << '\n';
  poly regcopy = regular;
  regular.get_basis(basis);
  cout << basis;
  regcopy.get_basis(basis);
  cout << basis;
  poly hermcopy = hermite;
  poly lagucopy = laguerre;
  cout << '\n';
  cout << "Lengths " << regcopy.get_len() << ' ' 
       << hermcopy.get_len() << ' ' << lagucopy.get_len() << '\n';

  cout << starbox("/Check operator = //") << '\n';
  regcopy = regular;
  regular.get_basis(basis);
  cout << basis;
  regcopy.get_basis(basis);
  cout << basis;
  hermcopy = hermite;
  lagucopy = laguerre;
  cout << '\n';
  cout << "Lengths " << regcopy.get_len() << ' ' 
       << hermcopy.get_len() << ' ' << lagucopy.get_len() << '\n';

  hermite.get_multi(multi);
  cout << multi[0] << multi[hermite.get_len()-1];
  hermcopy.get_multi(multi);
  cout << multi[0] << multi[hermcopy.get_len()-1];
  

  return 0;
}
