#include <iostream>
#include <sstream>
#include <fstream>
#include <vector>
#include <math.h>
#include <cstdlib>
using namespace std;

void read_in_twopt(vector<vector<double> > &rtwopt,
		   const string filename){

  string line;
  //  cout << "filename"<<filename<<endl;
  ifstream input(filename.c_str());
  if (input.is_open())
  {
    //first line gives number of configs and lattice size
    getline (input,line);
    //    cout << line<<endl;
    int nc;
    int tsize;
    int xsize;
    int dum1, dum2;
    sscanf(line.c_str(),"%d %d %d %d %d",&nc,&tsize,&dum1,&xsize,&dum2);
    //    cout <<"nc tsize "<<nc<<" "<<tsize<<endl;
    rtwopt.resize(tsize);
    for(int t = 0; t < tsize; t++){
      rtwopt[t].resize(nc);
    }
    int c = 0;
    int t = 0;
    while (! input.eof() && c < nc )
    {
      t = 0;
      while(t < tsize){
	//	cout << "t "<<t<<" c "<<c<<endl;
	getline (input,line);
	//cout << line << endl;
	int tdum;
	double real_p;
	double imag_p;
	sscanf(line.c_str(),"%d %lf %lf",&tdum,&real_p,&imag_p);
	rtwopt[tdum][c]=real_p;
	//	cout <<t<<" "<<c<<" "<<rtwopt[t][c]<<endl;
	t++;
      }
      c++;
    }
    if(c != nc || t != tsize){
      cout << "Incomplete reading in "<<filename<<" "<<c <<" "<<t<<endl;
    }
    input.close();
  }
  else {
    cout << "Unable to open file " << filename; 
    exit(1);
  }

}
void jack_avg_err2(const vector<double > &data,const double avg, double &var){

  const int conf = data.size();
  var = 0;
  for(int c = 0; c < conf; c++){
    var+=(avg - data[c])*(avg - data[c]);
  }
  var *= double(conf-1)/conf;  
  var = sqrt(var);
}

void avg_var_std(const vector<double> &data, double &avg, double &var){

  const int conf = data.size();
  avg = 0;
  var = 0;
  for(int c = 0; c < conf; c++){
    avg+=double(data[c]);
  }
  avg /= double(conf);
  for(int c = 0; c < conf; c++){
    var+=(avg - double(data[c]))*(avg - double(data[c]));
  }
  var /= double(conf*(conf-1));
  var = sqrt(var);
}
void average_data(const vector<vector<double > > & twopt,vector<double> &avg, vector<double> &var){

  avg.resize(twopt.size());
  var.resize(twopt.size());
  for(int t = 0; t < twopt.size(); t++){
    avg_var_std(twopt[t],avg[t],var[t]);
  }

}  
void avg_std(const vector<double> &data, double &avg){

  const int conf = data.size();
  avg = 0;
  for(int c = 0; c < conf; c++){
    avg+=data[c];
  }
  avg /= double(conf);
}
void average_data(const vector<vector<double > > & twopt,const string filename){

  ostringstream output;
  vector<double> avg;
  vector<double> var;
  average_data(twopt,avg,var);
  for(int t=0; t < twopt.size(); t++){
    output << t <<" "<<avg[t]<<" "<<var[t]<<endl;
  }
  ofstream out_file(filename.c_str());
  out_file << output.str();
  out_file.close();
}  
void calc_effmass_exp(  vector<vector<double> > & twopt, string filename){

  const int lt = twopt.size();
  const int conf = twopt[0].size();
  ostringstream fl_name;
  fl_name<<"effmass_"<<filename;
  ofstream out_file(fl_name.str().c_str());
  
  vector<double> avg(lt);
  for(int t = 0; t < lt; t++){
    avg_std(twopt[t],avg[t]);
  }

  for(int t = 0; t < lt-1; t++){

    double rat = avg[t+1]/avg[t];
    if(t > lt/2) rat = 1./rat;
    if(rat > 0 ){

      double effmass=-log(rat);
      
      vector<double> jack_avg(conf);
      for(int jack=0; jack < conf; jack++){
	jack_avg[jack]=0.;
      }  
      bool neg = false;
      for(int jack=0; jack < conf; jack++){
	double propt = (avg[t]*conf - twopt[t][jack])/double(conf-1);
	double proptp1 = (avg[t+1]*conf - twopt[t+1][jack])/double(conf-1);
	if(proptp1/propt < 0) 
	  neg = true;
	else{
	  jack_avg[jack] = -log(proptp1/propt);
	  if(t>lt/2)jack_avg[jack] = -log(propt/proptp1);
	}
      }
      if(!neg){
	double err;
	jack_avg_err2(jack_avg,effmass,err);
	out_file << t+0.5<<" "<<effmass<<" "<<err<<endl;
      }
      else {
	out_file <<"# ratio of prop goes negative "<<t<<endl;
      }
    }
    else{
      out_file <<"# ratio of prop goes negative "<<t<<endl;
    }
  }
}


int main(int argc, char **argv)
{
  cout << "Calculates the effective mass of a correlator using -log(C(t+1)/C(t))"<<endl;
  if(argc < 2){
    cout <<"usage: effmass <filename>"<<endl;
    exit(1);
  }
  vector<vector<double> > twopt;

  cout << "Reading in 2pt function " << argv[1]<<endl;
  read_in_twopt(twopt,argv[1]);

  //two point functions
  cout<< "Calculating averages of twopt functions "<<endl;
  ostringstream filename;
  filename<<"avg_"<<argv[1];
  average_data(twopt,filename.str());
  cout<< "Calculating the effective mass "<<endl;
  calc_effmass_exp(twopt,argv[1]);
}