#include "libscl.h"
#include "lprior.h"
#include "lpeqns.h"

using namespace scl;
using namespace std;

int main(int argc, char** argp, char** envp)
{
  string filename;
  ifstream fin;

  filename = "../data/DCFdataL.txt";
  fin.open(filename.c_str());
  if (!fin) error("Error, cannot open " + filename);

  const INTEGER data_row = 85;  
  const INTEGER data_col = 4;

  realmat data(data_row,data_col);
  
  for (INTEGER t=1; t<=data_row; ++t) {
    REAL year, log_mrs, log_gdp_prof, del_log_prof;
    fin >> year >> log_mrs >> log_gdp_prof >> del_log_prof;
    data(t,1) = year;
    data(t,2) = log_mrs;
    data(t,3) = log_gdp_prof;
    data(t,4) = del_log_prof;
  }

  intvec idx(2);
  idx[1] = 2;
  idx[2] = 4;

  const INTEGER lag = 1;
  const INTEGER M = 2;

  if (M != idx.size()) error("Error, bad dimensions");

  realmat b0, B, S, C;

  INTEGER rank = varcoef(data, idx, lag, b0, B, S, C);

  realmat E;
  INTEGER ier;
  REAL maxlam = eigen(B,E,ier);

  cout << '\n';
  cout << "rank = " << rank << '\n';
  cout << "ier = " << ier << '\n';
  cout << "maxlam = " << maxlam << '\n';
  cout << '\n';
  cout << "b0 = " << b0 << '\n';
  cout << "B = " << B << '\n';
  cout << "S = " << S << '\n';

  const REAL pi = 4.0*atan(1.0);
  const REAL sqrt2 = sqrt(2.0);
  
  realmat IB = -B;
  for (INTEGER i=1; i<=M; ++i) IB(i,i) += 1.0;
  realmat IBinv;
  IBinv = inv(IB);

  realmat R=S;
  if(factor(R)!=0) error("Error, factor failed");

  cout << '\n';
  cout << "R = " << R << '\n';
  cout << "IB = " << IB << '\n';
  cout << "IBinv = " << IBinv << '\n';

  realmat y0(M,1);
  for (INTEGER i=1; i<=M; ++i) y0[i] = data(1,idx[i]);
  INTEGER mrs_pos = 1;
  INTEGER cf_pos = 2;
  INTEGER years = 30;
  realmat ecf,pvcf,pv1,cecf,dcf,yld;

  dcfyld(b0,B,S,y0,mrs_pos,cf_pos,years,ecf,pvcf,pv1,cecf,dcf,yld); 

  cout << '\n';
  cout << "years = " << years << '\n';
  cout << '\n';
  cout << "y0 = " << y0 << '\n';
  cout << "ecf = " << ecf << '\n';
  cout << "pvcf = " << pvcf << '\n';
  cout << "pv1 = " << pv1 << '\n';
  cout << "cecf = " << cecf << '\n';
  cout << "dcf = " << dcf << '\n';
  cout << "yld = " << yld << '\n';

  const INTEGER quad_n = 5;

  realmat x, w;

  if (hquad(quad_n,x,w)!=0) error("Error, hquad failed");

  cout << '\n';
  cout << "hquad x = " << x << '\n';
  cout << "hquad w = " << w << '\n';

  realmat z(M,1);

  realmat dcf_mean(years,1,0.0),yld_mean(years,1,0.0),cecf_mean(years,1,0.0);
  realmat dcf_var(years,years,0.0), yld_var(years,years,0.0);

  for (INTEGER i=1; i<=quad_n; ++i) {
    for (INTEGER j=1; j<=quad_n; ++j) {

      z[1] = sqrt2*x[i];
      z[2] = sqrt2*x[j];

      y0 = IBinv*(b0 + R*z);

      REAL weight = w[i]*w[j]/pi;

      dcfyld(b0,B,S,y0,mrs_pos,cf_pos,years,ecf,pvcf,pv1,cecf,dcf,yld);

      cecf_mean += weight*dcf;
      dcf_mean += weight*dcf;
      yld_mean += weight*yld;
    }
  }

  cout << '\n';
  cout << "cecf_mean = " << cecf_mean;
  cout << "dcf_mean = " << dcf_mean;
  cout << "yld_mean = " << yld_mean;
   
  for (INTEGER i=1; i<=quad_n; ++i) {
    for (INTEGER j=1; j<=quad_n; ++j) {

      z[1] = sqrt2*x[i];
      z[2] = sqrt2*x[j];

      y0 = IBinv*(b0 + R*z);

      REAL weight = w[i]*w[j]/pi;

      dcfyld(b0,B,S,y0,mrs_pos,cf_pos,years,ecf,pvcf,pv1,cecf,dcf,yld);

      dcf_var += weight*(dcf-dcf_mean)*T(dcf-dcf_mean);
      yld_var += weight*(yld-yld_mean)*T(yld-yld_mean);
    }
  }

  cout << '\n';
  cout << "dcf_var = " << dcf_var;
  cout << "yld_var = " << yld_var;

  return 0;
}
