#include "libscl.h"
#include "crra_mf.h"

using namespace scl;
using namespace std;

int main(int argc, char** argp, char** envp)
{
  const INTEGER rows = 2;
  const INTEGER cols = 500;

  const INTEGER mflags = 0;
  const INTEGER HAC_lags = 0;

  string filename;
  ifstream fin;
  ofstream fout;

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

  realmat theta(2,1);

  string line;
  getline(fin,line);
  string theta_1 = cutstr(line);
  getline(fin,line);
  string theta_2 = cutstr(line);
  fin.clear(); fin.close();

  if(!isREAL(theta_1,theta[1])) error("Error, zsim, bad theta_1, " + theta_1);
  if(!isREAL(theta_2,theta[2])) error("Error, zsim, bad theta_2, " + theta_2);

  cout << "theta = " << theta << '\n';

  filename = "cbs.usrvar.sim";
  fin.open(filename.c_str());
  if (!fin) error("Error, cannot open " + filename);

  realmat simdat;
  REAL simvar;
  while(fin >> simvar) simdat.push_back(simvar);
  fin.clear(); fin.close();

  cout << "simdat.size() = " << simdat.size() << '\n';

  INTEGER n_data_rows = rows;
  INTEGER n_data_cols = simdat.size()/n_data_rows;

  realmat data(n_data_rows,n_data_cols);

  for (INTEGER j=1; j<=n_data_cols; ++j) {
    for (INTEGER i=1; i<=n_data_rows; ++i) {
      INTEGER ij = n_data_rows*(j-1) + i;
      data(i,j) = simdat[ij];
    }
  }

  cout << '\n';
  cout << "first 12 data cols" << data("",seq(1,12));
  cout << '\n';
  cout << "last 12 data cols" << data("",seq(n_data_cols-11,n_data_cols));
  cout << '\n';
  cout << "stats" << simple(T(data));
  cout.flush();

  mle::crra_moment_function crramf;

  bool rv;

  INTEGER first = 1;
  INTEGER last = first - 1 + cols;

  intvec idx = seq(1,rows);
  intvec jdx = seq(first,last);

  realmat subsample = data(idx,jdx);

  rv = crramf.set_data(&subsample);
  if (!rv) error("Error, zsim, crramf.set_data failed ");

  rv = crramf.set_L(mflags);
  if (!rv) error("Error, zsim, crramf.set_L failed ");

  gmm crragmm(&crramf,mflags,&subsample,cols,HAC_lags);

  const bool correct_W_for_mean = true;
  crragmm.set_correct_W_for_mean(correct_W_for_mean);

  const bool regularize_W = true;
  const REAL ridge = 0.0;
  crragmm.set_regularize_W(regularize_W,ridge);

  const bool warning_messages = true;
  crragmm.set_warning_messages(warning_messages);

  cout << '\n';
  cout << "moment_function_lags = " << mflags << '\n';
  cout << "Lhac = " << HAC_lags << '\n';
  cout << boolalpha;
  cout << "correct_W_for_mean set to " << correct_W_for_mean << '\n';
  cout << "regularize_W set to (" << regularize_W << ", " << ridge <<')'<<'\n';
  cout << "warning_messages set to " << warning_messages << '\n';
  cout.flush();

  cout << '\n';
  cout << "subsample first = " << subsample(idx,1) << '\n';
  cout << "subsample last  = " << subsample(idx,last) << '\n';
  cout << "first = " << first << '\n';
  cout << "last = " << last << '\n';
  cout << "sample size = " << cols << '\n';
  cout.flush();

  INTEGER d = crramf.get_d();
  INTEGER n_Z_reps = n_data_cols/cols - 1;

  realmat Z_save(n_Z_reps,d);

  first = 1;
  last  = first - 1 + cols;

  realmat z;

  cout << '\n';
  cout << "compute Z" << '\n';

  for (INTEGER rep=1; rep<=n_Z_reps; ++rep) {

    cout << "first, last = " << first << ' ' << last << '\n';
    jdx = seq(first,last);
    subsample = data(idx,jdx);
    rv = crragmm.set_data(&subsample);
    if (!rv) error("Error, zsim, crragmm.set_data failed");

    denval likelihood = crragmm.likelihood(theta,z);

    INTEGER ier = crragmm.get_W_numerr();

    if (ier != 0) warn("Warning, zsim, W_numerr false");
    if (!likelihood.positive) warn("Warning, zsim, likelihood fail");

    for (INTEGER j=1; j<=d; ++j) Z_save(rep,j) = z[j];

    first += cols;
    last  += cols;
  }

  vector<string> header(d);
  header[0] = "e";
  header[1] = "lsrlag*e";
  header[2] = "lcglag*e";

  filename = "Z.txt";
  writetable(filename.c_str(),Z_save,header,25,16);

  return 0;
}
