/***************************************************************************
                          analysis_mod.cpp  -  description
                             -------------------
    begin                : Sat Jul 12 00:34:50 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 <zeta/module.h>
#include <zeta/class.h>
#include <zeta/zextclass.h>
#include <zeta/zobj.hpp>
#include <analysis.h>

using namespace Almost;

#define ALMTOSTRING(s) \
template<> \
inline string to_string(const s & ){ \
  return "<" #s ">"; \
}

ALMTOSTRING(Molecules);
ALMTOSTRING(Coor<float>);
ALMTOSTRING(DCDTrj);
ALMTOSTRING(NameSelection);

template<>
inline string to_string(const DataFile & sp){
  return "<DataFile>";
}

template<> string to_string(const DataSet & ds){
  return "<DataSet>";
}

template<> string to_string(const Wham & ds){
  return "<Wham>";
}

template<> string to_string(const XPlorNOE & ){
  return "<XPlorNOE>";
}

template<> string to_string(const S2 & ){
  return "<S2>";
}

template<> string to_string(const S2RMSD & ){
  return "<S2RMSD>";
}

template<> string to_string(const jcoupling & ){
  return "<jcoupling>";
}

template<> string to_string(const Velocity & ){
  return "<velocity>";
}


//go detailed balance test
void detailed_balance_test(string data_file, string out_file,
			   int col, double min, double max, int bins){
  cout<<"\t>> Reading data ....\n";
  DataFile df(data_file);
  vector<double>  ds = df.col(col-1);
  vector<double>    ene_hist;
  map<int,map<int,double> > markov;
  cout<<"\t>> Done\n";

  cout<<"\t>> Generating energy histogram ....\n";
  //do energy-histo
  ene_hist.resize(bins+1);
  double e_min = min;
  double e_max = max;
  double delta = (e_max-e_min)/(double)bins;  

  for(int i=0;i<ds.size();++i){
    if((ds[i]<e_min)||(ds[i]>e_max)) throw "Energies out of window";
    int ind = (int)ceil((ds[i]-e_min)/delta);
    ++ene_hist[ind];
  }

  for(int i=0;i<ene_hist.size();i++){
    ene_hist[i] /= ds.size()*delta;
    //    cout<<ene_hist[i]<<endl;
  }

  cout<<"\t>> Done\n";

  cout<<"\t>> Building Markov matrix ....\n";
  //prob matrix
  //for(int i=0;i<ene_hist.size();i++) markov[i][i] =1;
  for(int i=0;i<ds.size()-1;i++){
    int ind1 = (int)ceil((ds[i]-e_min)/delta);
    int ind2 = (int)ceil((ds[i+1]-e_min)/delta);
    markov[ind1][ind2] +=1;
  }
  cout<<"\t>> Done\n";
  cout<<"\t>> Normalizing ....\n";
  //normalize
  {
    map<int,map<int,double> >::iterator iter = markov.begin(),
      end = markov.end();
  
    while(iter!=end){
      double N_tot = 0;
      map<int,double>::iterator iter2 = iter->second.begin(), 
	end2 = iter->second.end();
    
      while(iter2!=end2){
	N_tot += iter2->second;
	++iter2;
      }
      if(N_tot!=0){
	iter2 = iter->second.begin(), 
	  end2 = iter->second.end(); 
	while(iter2!=end2){
	  iter2->second /=N_tot;
	  ++iter2;
	}
      }
      ++iter;
    }
  }
  cout<<"\t>> Done\n";
  //db-test
  cout<<"\t>> Writing data ....\n";
  ofstream out;
  out.open(out_file.c_str());

  for(int i=0;i<ene_hist.size();i++){
    for(int j=i;j<ene_hist.size();j++){
      //Test N_a P_AB = N_B P_BA
      // if((ene_hist[i]!=0)&&(ene_hist[j]!=0))
	{
	out<<i<<"\t"<<j<<"\t"<<ene_hist[i]*markov[i][j] - ene_hist[j]*markov[j][i]
	   <<"\t"<<ene_hist[i]*markov[i][j]<<"\t"<<ene_hist[j]*markov[j][i]
	   <<"\t"<<markov[i][j]<<"\t"<<markov[j][i]<<"\n";
	}
    }
  }
  out.close();
  cout<<"\t>> Done\n";
}
      
extern "C" {
  void init_analysis(){
    Module mod = Module("analysis","Analysis Tools")
      .def_function("mean","Computes mean value of column n in the data set",Almost::mean)
      .def_function("var","Computes variance of column n in the data set",Almost::var)
      .def_function("min","Computes min value of column n in the data set",Almost::min)
      .def_function("max","Computes max value of column n in the data set",Almost::max)
      .def_function<DataSet (*)(const DataFile & , int , double, double, int)>("hist","Creates histogram",hist)
      .def_function("export","Writes data set to file",export_ds)
      .def_function("dbt","Does detailed balance test",detailed_balance_test)
      .def_function("dihedral","Mesures dihedral of selected atoms.",dihedral)
      .def_function("correlation","Computes dataset correlation.",correlation)
      ;



    Class<DataFile,ArgList1<string> >(mod.self(),"data_file")
      .def_method("cols",&DataFile::cols)
      .def_method("rows",&DataFile::rows)
      .def_method("data_set",&DataFile::data_set)
      ;

    Class<DataSet>(mod.self(),"data_set");

//     ZObj::BinaryOperator<INDEX>::BOM[pair<TypeInfo,TypeInfo>
// 				     (typeid(ZType<DataSet &>),
// 				      typeid(ZType<doub))] =
//     index_operator<DataSet ,double,DataSet &>;


    Class<Wham>(mod.self(),"wham")
      .def_method("add_file",&Wham::add_file)
      .def_method("compute",&Wham::compute)
      .def_method("logZ",&Wham::logZ)
      .def_method("U",&Wham::U)
      .def_method("cv",&Wham::cv)
      .def_method("logR",&Wham::logR)
      .def_method("p",&Wham::p)
      ;

    Class<XPlorNOE,ArgList2<const Molecules &,int> >(mod.self(),"XPlorNOE")
      .def_method("assign",&XPlorNOE::assign)
      .def_method("eval",&XPlorNOE::eval)
      .def_method("eval_coor",&XPlorNOE::eval_coor);

    Class<S2,ArgList1<const Molecules &> >(mod.self(),"s2")
      .def_method("assign",&S2::assign)
      .def_method("eval",&S2::eval);

    Class<S2RMSD,ArgList2<const Molecules &, const NameSelection &> >(mod.self(),"s2rmsd")
      .def_method("assign",&S2RMSD::assign)
      .def_method("eval",&S2RMSD::eval);

    Class<jcoupling,ArgList1<const Molecules &> >(mod.self(),"jcoupling")
      .def_method("assign",&jcoupling::assign)
      .def_method("eval",&jcoupling::eval)
      .def_method("eval_mol",&jcoupling::eval_mol);
    ;

    Class<Velocity,ArgList2<const Molecules &, NameSelection &> >(mod.self(),"velocity")
      .def_method("temp",&Velocity::temp);
  }
}