#include "libscl.h"

/*
extern void mvnquad (INTEGER n, INTEGER M, 
                      const scl::realmat& mu, const scl::realmat& sig, 
		      scl::realmat& xvals, scl::realmat& probs);
*/

using namespace scl;
using namespace std;

int main(int argc, char** argp, char** envp)
{

  INTEGER M = 3;
  INTEGER n = 3;

  realmat mu(M,1,0.0);
  realmat sig(M,M,0.0);
  for (INTEGER i=1; i<=M; ++i) sig(i,i) = 1.0;

  realmat xvals;
  realmat probs;

  mvnquad(n,M,mu,sig,xvals,probs);

  INTEGER len = probs.size();

  realmat muhat(M,1,0.0);
  realmat sighat(M,M,0.0);
  realmat x(M,1);

  for (INTEGER i=1; i<=len; ++i) {
    for (INTEGER j=1; j<=M; ++j) x[j] = xvals(i,j);
    muhat += probs[i]*x;
    sighat += probs[i]*(x*T(x));
  }

  cout << mu << muhat;
  cout << sig << sighat - mu*T(mu);

  for (INTEGER i=1; i<=M; ++i) {
    mu[i] = i;
    for (INTEGER j=1; j<=M; ++j) {
      if (i == j) {
        sig(i,j) = i;
      }
      else {
        sig(i,j) = (i + j) / 10.0; 
      }
    }
  }

  mvnquad(n,M,mu,sig,xvals,probs);

  len = probs.size();

  for (INTEGER i=1; i<=len; ++i) {
    for (INTEGER j=1; j<=M; ++j) x[j] = xvals(i,j);
    muhat += probs[i]*x;
    sighat += probs[i]*(x*T(x));
  }

  cout << mu << muhat;
  cout << sig << sighat - mu*T(mu);

  M = 2;
  n = 4;

  mu.resize(M,1);
  sig.resize(M,M,0.0);

  mu[1] = 0.95;
  mu[2] = 12.5; 

  sig(1,1) = pow(0.01,2);
  sig(2,2) = pow(2.0,2);

  mvnquad(n,M,mu,sig,xvals,probs);

  len = probs.size();
  REAL sum = 0.0;
  for (INTEGER i=1; i<=len; ++i) sum += probs[i];

  cout << xvals << probs << "\n\t\t\t\t  sum = " << sum << '\n';

  intvec support;
  sum = 0.0;
  for (INTEGER i=1; i<=len; ++i) {
    REAL beta = xvals(i,1);
    REAL gamma = xvals(i,2);
    bool in_support = true;
    if ( (beta  <= 0.8) || (0.99 <= beta  ) ) in_support = false;
    if ( (gamma <= 0.0) || (100.0 <= gamma) ) in_support = false;
    if (in_support) {
      support.push_back(i); 
      sum += probs[i];
    }
  }
 
  xvals = xvals(support,"");
  probs = probs(support,"")/sum;

  len = probs.size();
  sum = 0.0;
  for (INTEGER i=1; i<=len; ++i) sum += probs[i];

  cout << xvals << probs << "\n\t\t\t\t  sum = " << sum << '\n';
 
  writetable("quadpoints.txt",cbind(xvals,probs),24,16);

  string filename;
  ofstream fout;
  ifstream fin;

  vector<string> parmfile;

  filename = "cbs.parm.000";
  fin.open(filename.c_str());
  if (!fin) error("Error, cannot open " + filename);
  string line;
  while(getline(fin,line)) parmfile.push_back(line);
  fin.clear(); fin.close();

  for (INTEGER i=1; i<=len; ++i) {
    line = fmt('f',20,16,xvals(i,1))() + "     Beta = crra_parms[1], REAL"; 
    parmfile[35] = line;
    line = fmt('f',20,16,xvals(i,2))() + "     Gamma = crra_parms[2], REAL";
    parmfile[36] = line;
    parmfile[37] 
      = "         0.0             Log base 10 penalty, log10lam, float";
    filename = "cbs.parm_" + fmt('d',2,i)('0');
    fout.open(filename.c_str());
    if (!fout) error("Error, cannot open " + filename);
    vector<string>::const_iterator iter;
    for (iter=parmfile.begin(); iter!=parmfile.end(); ++iter) {
      fout << *iter << '\n';
    }
    fout.clear(); fout.close();
  }

  return 0;
}
