#include "libscl.h"
#include "emitted_code.h"

using namespace scl;
using namespace std;

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

  string program(argp[0]);
  
  if (argc != 5) {
    cerr << '\n';
    cerr << "Usage:   " << program << " prefix start stop l_cmf_parms" << '\n';
    cerr << "Example: " << program << " sdf.parms_chain 0 9" << '\n';
    cerr << "Comment: sdf.parms_chain.000.dat would be first file read"<<'\n';
    cerr << "         sdf.parms_chain.009.dat would be last for exmple."<<'\n';
    cerr << '\n';
    return 1;
  }

  string rootname(argp[1]);
  INTEGER start(atoi(argp[2]));
  INTEGER stop(atoi(argp[3]));
  INTEGER l_cmf_parms(atoi(argp[4]));

  cout << '\n';
  cout << program << " settings;" << '\n';
  cout << "  rootname = " << rootname << '\n';
  cout << "  start = " << start << '\n';
  cout << "  stop = " << stop << '\n';
  cout << "  l_cmf_parms  = " << l_cmf_parms << '\n';
  cout << '\n';

  string::iterator sitr = rootname.end();

  if (*sitr != '.') rootname.push_back('.');

  string filename = rootname + fmt('d',3,start)('0') + ".dat";
  realmat parms_k;
  vecread(filename.c_str(),parms_k); 
  cout << "Read " << filename << '\n';

  INTEGER rows = parms_k.nrow();
  INTEGER cols = parms_k.ncol();

  INTEGER K = stop - start + 1;
  INTEGER max_cols = K*cols;

  realmat parms(rows,max_cols);

  INTEGER loc = 0;

  for (INTEGER j=1; j<=cols; ++j) {
    if (j > max_cols) break;
    for (INTEGER i=1; i<=rows; ++i) {
      parms(i,loc+j) = parms_k(i,j);
    }
  }

  loc += cols;

  for (INTEGER k=start+1; k<=stop; ++k) {
    string filename = rootname + fmt('d',3,k)('0') + ".dat";
    vecread(filename.c_str(),parms_k); 
    cout << "Read " << filename << '\n';

    if (rows != parms_k.nrow()) error("Error, number of rows differ");

    INTEGER cols = parms_k.ncol();
    
    for (INTEGER j=1; j<=cols; ++j) {
      if (j > max_cols) break;
      for (INTEGER i=1; i<=rows; ++i) {
        parms(i,loc+j) = parms_k(i,j);
      }
    }

    loc += cols;
  }

  cout << '\n';
  
  INTEGER da = rows;
  INTEGER n = parms.ncol();
  INTEGER l_snp_parms = da - l_cmf_parms;

  cout << "da = " << da << '\n';
  cout << "n  = " << n << '\n';
  cout << "l_snp_parms = " << l_snp_parms << '\n';
  cout << "l_cmf_parms = " << l_cmf_parms << '\n';
  cout << '\n';

  realmat x(da,1);
  realmat snp_parms(l_snp_parms,1);
  realmat cmf_parms(l_cmf_parms,1);

  realmat mean(da,1, 0.0);

  for (INTEGER t=1; t<=n; ++t) {
    for (INTEGER i=1; i<=da; ++i) x[i] = parms(i,t);
    for (INTEGER i=1; i<=l_snp_parms; ++i) snp_parms[i] = x[i];
    cmf_parms = theta(snp_parms,l_cmf_parms);
    for (INTEGER i=1; i<=l_cmf_parms; ++i) x[l_snp_parms + i] = cmf_parms[i];
    mean += x;
  }
  mean = mean/n;
  realmat V(da,da,0.0);
  for (INTEGER t=1; t<=n; ++t) {
    for (INTEGER i=1; i<=da; ++i) x[i] = parms(i,t);
    for (INTEGER i=1; i<=l_snp_parms; ++i) snp_parms[i] = x[i];
    cmf_parms = theta(snp_parms,l_cmf_parms);
    for (INTEGER i=1; i<=l_cmf_parms; ++i) x[l_snp_parms + i] = cmf_parms[i];
    V += (x-mean)*T(x-mean);
  }
  V = V/n;

  filename = rootname + "V.dat";
  vecwrite(filename.c_str(),V);

  cout << '\n';
  cout << "Wrote " << filename << '\n';

  realmat sdev(da,1);
  for (INTEGER i=1; i<=da; ++i) sdev[i] = sqrt(V(i,i));

  filename = rootname + "sdev.dat";
  vecwrite(filename.c_str(),sdev);

  cout << "Wrote " << filename << '\n';

  filename = rootname + "sdev.txt";
  ofstream fout;
  fout.open(filename.c_str());
  if (!fout) error("Error, cannot open " + filename);
  for (INTEGER i=1; i<=da; ++i) fout << fmt('f',8,5,sqrt(V(i,i))) << '\n';
  fout.clear(); fout.close();

  cout << "Wrote " << filename << '\n';
  cout << '\n';
  cout << "Which contained" << '\n';
  cout << '\n';
  for (INTEGER i=1; i<=da; ++i) cout << fmt('f',8,5,sqrt(V(i,i))) << '\n';
  cout << '\n';

  return 0;
}
