#include "libscl.h"
using namespace scl;
using namespace std;

int main(int argc, char **argp, char **envp)
{
  if (argc != 5) {
    cerr << '\n';
    cerr << "Usage: " << argp[0] << " path temperature nodes reps\n";
    cerr << "Example: " << argp[0] << " results 5 15 5\n";
    cerr << '\n';
    return 1;
  }

  string pathname(argp[1]);
  pathname = pathname + string("/");
  REAL temp = atof(argp[2]);
  INTEGER n_snd = atoi(argp[3]);
  INTEGER n_num = atoi(argp[4]);

  cout << "pathname = " << pathname << '\n';
  cout << "temp = " << temp << '\n';
  cout << "n_snd = " << n_snd << '\n';
  cout << "n_num = " << n_num << '\n';
  cout << '\n';
  
  const realmat null;

  string fullname;
  char sender[5], number[5], filename[80];

  realmat theta;

  realmat grand_mu;
  realmat grand_S;
  realmat grand_sd;
  INTEGER grand_n=0;

  INTEGER n = 0;
  INTEGER p = 0; 

  for (INTEGER snd=1; snd<=n_snd; ++snd) {
    for (INTEGER num=0; num<=n_num; ++num) {
  
      sprintf(sender,"%03d",snd);
      sprintf(number,"%03d",num);
    
      filename[0]='\0';
      strcat(filename,"theta");
      strncat(filename,sender,3);
      strncat(filename,number,3);
      strcat(filename,".dat");
      fullname = pathname+filename;

      ifstream sim_stream(fullname.c_str());

      vecread(sim_stream, theta);

      cout << filename << '\n';

      p=theta.nrow();
      n=theta.ncol();

      if (grand_mu == null) {
        grand_mu.resize(p,1,0.0);
        grand_S.resize(p,p,0.0);
	grand_sd.resize(p,1);
	grand_n = 0;
      }

      realmat mu(p,1,0.0);
      for (INTEGER t=1; t<=n; t++) {
        mu += theta("",t);
      }
      grand_mu += mu;

      mu = mu/n;

      realmat S(p,p,0.0);
      for (INTEGER t=1; t<=n; t++) {
        for (INTEGER j=1; j<=p; j++) {
	  for (INTEGER i=1; i<=p; i++) {
	    S(i,j) += (theta(j,t)-mu[j])*(theta(i,t)-mu[i]);
          }
        }
      }
      grand_S += S;

      S = S/n;
      S = temp*S;

      realmat sd(p,1);
      for (INTEGER i=1; i<=p; i++) sd[i] = sqrt(S(i,i));

      cout << cbind(mu,sd);
      grand_n += n;
    }
  }

  grand_mu = grand_mu/grand_n;
  grand_S = grand_S/grand_n;
  grand_S = temp*grand_S;

  for (INTEGER i=1; i<=p; i++) grand_sd[i] = sqrt(grand_S(i,i));

  cout << "grand mean & std. dev." << cbind(grand_mu,grand_sd);
  cout << "grand S" << grand_S << '\n';
  cout << "grand n = " << grand_n << '\n';

  //vecwrite("theta_hat.sav", grand_mu);
  //vecwrite("V_hat.sav", grand_S);

  vecwrite("theta_mu.dat", grand_mu);
  vecwrite("theta_sig.dat", grand_S);

  cout << "\n\n";
  for (INTEGER i=1; i<=p; i++) {
    cout << "const REAL XXX_start  = " << grand_mu[i] << ";\n";
  }
  cout << "\n\n";
  for (INTEGER j=1; j<=p; j++) {
    for (INTEGER i=1; i<=p; i++) {
      cout << "def[0].Vmat("<<i<<","<<j<<") = "<<grand_S(i,j)<<";\n";
    }
  }

  return 0;
}