/***************************************************************************
                          analysis.cpp  -  description
                             -------------------
    begin                : Fri Jul 11 00:49:53 2003
    copyright            : (C) 2002 by Cavalli Andrea
    email                : cavalli@bioc.unizh.ch
**************************************************************************/
/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#include <analysis.h>


namespace Almost {


      
  void Wham::add_file(string file, double tmp, int col){
    info_[file]= pair<double,int>(tmp,col);
  }

  void Wham::compute(double min, double max, int bins, double err){
    //read data
    cout<<"\t>> Reading data ..."<<"\n";
    {    
      map<string,pair<double,int> >::iterator iter = info_.begin(),
	end = info_.end();

      while(iter!=end){
	string file = iter->first;
	double temp = iter->second.first;
	int    col  = iter->second.second;
	cout<<"\t>> Reading "<<file<<"\n";
	DataFile df(file.c_str());
      
	if(data_.find(temp)==data_.end()){
	  data_[temp] = df.col(col);
	}
	else {
	  for(int i=0;i<df.col(col).size();++i)
	    data_[temp].push_back(df.col(col)[i]);
	}
	++iter;
      }
    }

    if(max<=min) throw "max >= min";
    min_ = min;
    max_ = max;
    dx_ = (max_ - min_)/(double)bins;
    
    //compute histo  
    { 
      map<double, vector<double> >::iterator iter = data_.begin(),
	end = data_.end();
      while(iter!=end){
	for(int i = 0;i<iter->second.size();i++){
	  int h = (int)ceil(iter->second[i]/dx_);
	  e_hist_[h] += 1;
	}
	++iter;
      }
    }

    
    //lets wham ....
    map<double,double> logZ1_;
    map<double,double> logZ2_;
    
    {
      map<double, vector<double> >::iterator  iter = data_.begin(),
	end = data_.end();
      while(iter!=end){
	double temp = iter->first;
	logZ1_[temp] = 100;
	logZ2_[temp] = 100;
	N[temp]   = iter->second.size();
	++iter;
      }
    }


    while(1){
      map<double,double>::iterator iter = logZ2_.begin(), end = logZ2_.end();

      while(iter!=end){
	// compute den
	double t_k = iter->first;
	iter->second = 0;
	//loop on all energies
	map<int,int>::iterator e_iter = e_hist_.begin(),
	  e_end = e_hist_.end();
	while(e_iter!=e_end){
	  double E = dx_*e_iter->first;
	  int N_E  = e_iter->second;
	  //compute den
	  map<double,double>::iterator iter1 = logZ1_.begin(), 
	    end1 = logZ1_.end();
	  double x = log(N[iter1->first])-iter1->second
	    + (1/t_k-1/iter1->first)*E;
	  ++iter1;
	  while(iter1!=end1){
	    double y = log(N[iter1->first])-iter1->second
	      + (1/t_k-1/iter1->first)*E - x;
	    y = exp(y);
	    x += log1p(y);
	    ++iter1;
	  }
	  double temp = log((double)N_E)-x-iter->second;
	  temp = exp(temp);
	  iter->second += log1p(temp);
	  //	  iter->second +=  exp(-x)*N_E;	  
	  ++e_iter;
	}
	++iter;
      }

      double error = 0;
      map<double,double>::iterator zi = logZ1_.begin(), 
	ze = logZ1_.end();
      while(zi!=ze){
	double t = zi->first;
	//error += ((Z2_[t]-Z1_[t])/Z2_[t])*((Z2_[t]-Z1_[t])/Z2_[t]);
	error += (logZ2_[t]-logZ1_[t])/logZ2_[t]*(logZ2_[t]-logZ1_[t])/logZ2_[t];
	++zi;
      }


      logZ1_ = logZ2_;


      if(error<10e-10){
	cout<<"\n\t>> Current Error "<<error<<"\n";
	cout<<"\t>> Exiting loop"<<endl;
	break;
      }
      else{
 	cout<<"\n\t>> Current Error "<<error;
      }
    }

    logZ_ = logZ1_;
    logZ_buff_ = logZ1_;
  }


  DataSet Wham::logZ(double min, double max, int points){

    //construct return map
    map<double,double> res;
    double dx = (max-min)/points;
    double t = min;

    while(t <= max){
      res[t] = logZ__(t);
      t +=dx;
    }
    return res;
  }

  double Wham::logZ__(double t){

    //look if there
    // map<double,double>::iterator iter = logZ_buff_.find(t);
    //if(iter!=logZ_buff_.end()) return iter->second;

    double z = 0;
    double t_k = t;
    //loop on all energies
    map<int,int>::iterator e_iter = e_hist_.begin(),
      e_end = e_hist_.end();

    while(e_iter!=e_end){
      double E = dx_*e_iter->first;
      int N_E  = e_iter->second;
      //compute den
      map<double,double>::iterator iter1 = logZ_.begin(), 
	end1 = logZ_.end();
      double x = log(N[iter1->first])-iter1->second
	+ (1/t_k-1/iter1->first)*E;
      ++iter1;
      while(iter1!=end1){
	double y = log(N[iter1->first])-iter1->second 
	  + (1/t_k-1/iter1->first)*E - x;
	y = exp(y);
	x += log1p(y);
	++iter1;
      }
      double temp = log((double)N_E)-x-z;
      temp = exp(temp);
      z += log1p(temp);
      ++e_iter;
    }
    return z;
  }

  
  DataSet Wham::U(int p, double min, double max, int points){
    //construct return map
    map<double,double> res;
    double dx = (max-min)/points;
    double t = min;

    while(t<=max){
      double u = 0;
      double t_k = t;
      //loop on all energies
      map<int,int>::iterator e_iter = e_hist_.begin(),
	e_end = e_hist_.end();

      while(e_iter!=e_end){
	double E = dx_*e_iter->first;
	int N_E  = e_iter->second;
	//compute den
	map<double,double>::iterator iter1 = logZ_.begin(), 
	  end1 = logZ_.end();
	double x = log(N[iter1->first])-iter1->second
	  + (1/t_k-1/iter1->first)*E;
	++iter1;
	while(iter1!=end1){
	  double y = log(N[iter1->first])-iter1->second
	    + (1/t_k-1/iter1->first)*E - x;
	  y = exp(y);
	  x += log1p(y);
	  ++iter1;
	}
	u += pow(E,p)*exp(log((double)N_E)-x-logZ__(t));
	++e_iter;
      }
      res[t] = u;
      cout<<"U "<<p<<" "<<t<<" "<<res[t]<<"\n";
    t += dx;
    }
    return res;
  }    

  DataSet Wham::cv(double min, double max, int points){
    DataSet e = U(1,min,max,points);
    DataSet e2 = U(2,min,max,points);

    DataSet res;
    DataSet::iterator iter = e.begin(),end = e.end();
    while(iter!=end){
      double t = iter++->first;
      res[t] = (e2[t]-e[t]*e[t])/(t*t);
    }

    return res;

  }


  DataSet Wham::logR(){
    //loop on all energies
    map<int,int>::iterator e_iter = e_hist_.begin(),
      e_end = e_hist_.end();
    while(e_iter!=e_end){
      double E = dx_*e_iter->first;
      int N_E  = e_iter->second;
      //compute den
      map<double,double>::iterator iter1 = logZ_.begin(), 
	end1 = logZ_.end();
      double x = log(N[iter1->first])-iter1->second
	+ (-1/iter1->first)*E;
      ++iter1;
      while(iter1!=end1){
	double y = log(N[iter1->first])-iter1->second
	  + (-1/iter1->first)*E - x;
	y = exp(y);
	x += log1p(y);
	++iter1;
      }
      logR_[E] = -x+log((double)N_E);
      ++e_iter;
    }


    //write error log
//     ofstream err;
//     err.open("wham.dat");
//     map<double,double>::iterator r_iter = logR_.begin(),
//       r_end = logR_.end();
//     while(r_iter!=r_end){
//       err<<r_iter->first<<"\t"<<1/sqrt(N[r_iter->first])<<"\n";
//       ++r_iter;
//     }
    return logR_;
  }


  DataSet Wham::p(double t){
    map<double,double> res;
    
    map<double,double>::iterator iter = logR_.begin(), end = logR_.end();


    while(iter!=end){
      double E = iter->first;
      res[E] = iter->second-1/t*E-logZ__(t);
      res[E] = exp(res[E]);      
      ++iter;
    }

    return res;
  }
    

  double mean(const DataFile & df, int col){
    if(col<0) throw "invalid column index";
    if(col>=df.cols()) throw "invalid column index";

    if(df[col].size()==0) throw "empty column";
    
    double m = 0;
    int size = df[col].size();

    for(int i=1;i<size;++i){
      m += df[col][i];
    }
    return m/size;
  }

  double var(const DataFile & df, int col){

    if(col<0) throw "invalid column index";
    if(col>=df.cols()) throw "invalid column index";

    if(df[col].size()==0) throw "empty column";
    
    double val = df[col][0];
    double valP2 = val*val;
    double min = val;
    double max = val;

    int size = df[col].size();

    for(int i=1;i<size;++i){
      double d = df[col][i];
      val += d;
      valP2 += d*d;
      if(d>max) max = d;
      if(d<min) min = d;
    }

    val /= size;
    valP2 /= size;

    cout<<"\n\t>> Variance:   "<<valP2-val*val<<"\n";
    cout<<  "\t>> Minimum:    "<<min<<"\n";
    cout<<  "\t>> Maximum:    "<<max<<"\n";
    cout<<  "\t>> Mean:       "<<val<<"\n\n";

    return valP2-val*val;
  }


  double min(const DataFile & df, int col){
    if(col<0) throw "invalid column index";
    if(col>=df.cols()) throw "invalid column index";
    double val = df[col][0];
    double min = val;
    int size = df[col].size();

    for(int i=1;i<size;++i){
      double d = df[col][i];
      if(d<min) min = d;
    }

    return min;
  }


  double max(const DataFile & df, int col){
    if(col<0) throw "invalid column index";
    if(col>=df.cols()) throw "invalid column index";
    double val = df[col][0];
    double max = val;
    int size = df[col].size();

    for(int i=1;i<size;++i){
      double d = df[col][i];
      if(d>max) max = d;
    }

    return max;
  }

  DataSet hist(const DataFile & df, int col, int num){
    if(col<0) throw "invalid column index";
    if(col>=df.cols()) throw "invalid column index";
    if(num<=0) throw "invalid number of bins";
    DataSet ds;
    
    double min_ = min(df,col);
    double max_ = max(df,col);

    double dx = (max_ - min_)/num;

    for(int i = 0; i<df.rows();i++){
      double h = ceil(df[col][i]/dx)*dx;
      ds[h] += 1;
    }

    return ds;
  }

  DataSet hist(const DataFile & df, int col, double min, double max, int num){
    if(col<0) throw "invalid column index";
    if(col>=df.cols()) throw "invalid column index";
    if(num<=0) throw "invalid number of bins";
    DataSet ds;
    

    double dx = (max - min)/num;
    double Ninv = 1.0/df.rows();
    double w = Ninv/dx;
    double sum = 0;
    for(int i = 0; i<df.rows();i++){
      double h = ceil((df[col][i])/dx)*dx;
      ds[h] += w;
      sum += w;
    }
    return ds;
  }

  void export_ds(const DataSet & ds, string file){
    ofstream out;
    out.open(file.c_str());

    DataSet::const_iterator iter = ds.begin(),end = ds.end();
    while(iter!=end){
      out<<iter->first<<"\t"<<iter->second<<"\n";
      ++iter;
    }
  }
};