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

using namespace scl;
using namespace std;

int main(int argc, char** argp, char** envp)
{
  const INTEGER first = 1;
  const INTEGER last = 16;
  const INTEGER ndraws = 1000;

  const INTEGER d = 3;

  INT_32BIT seed = 8005629;

  realmat quadpoints;
  readtable("quadpoints.txt",quadpoints);

  cout << quadpoints("",3);

  if (last != quadpoints.nrow()) error("Error, zsample, last != quadpoints");
  realmat cumprob(last,1);
  cumprob[1] = quadpoints(1,3);
  for (INTEGER i=2; i<=last; ++i) cumprob[i] = cumprob[i-1] + quadpoints(i,3);

//cerr << cumprob << '\n';

  realmat Z_save(ndraws,d);

  realmat probs(last,1,0.0);

  for (INTEGER rep=1; rep<=ndraws; ++rep) {
    REAL u = ran(seed);
    INTEGER k = last;
    for (INTEGER i=first; i<=last; ++i) {
      if (u <= cumprob[i]) {
        k = i;
	break;
      }
    }
    probs[k] += 1.0;
    string filename = "../pedist_" + fmt('d',2,k)('0') + "/Z.txt";
    vector<std::string> header;
    realmat zsim;
    readtable(filename.c_str(),zsim,header);
    INTEGER rows = zsim.nrow();
    k = iran(seed,rows-1) + 1; 
    for (INTEGER j=1; j<=d; ++j) Z_save(rep,j) = zsim(k,j);
  }

  cout << probs/ndraws;

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

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

  return 0;
}
