/***************************************************************************
                          dcd.cpp  -  description
                             -------------------
    begin                : Thu Jul  3 00:23: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 <zeta/module.h>
#include <zeta/class.h>
#include <zeta/zextclass.h>
#include <zeta/zobj.hpp>
#include <dcd.h>
//#include <dynatypes.h>
#include <select.h>
#include <algorithm/maxsub.h>

using namespace Almost;


class MaxSubClusterDistance {
  Almost::MaxSub ms;
  float val;
public:
  MaxSubClusterDistance(const Molecules & m,float cutoff,float v):
    ms(m,cutoff), val(v){ ms.set_verbose(false);}
  bool in_cluster(double v){
    if(v>=val) return true;
    return false;
  }
  string description(){ return "MAXSUB";}
  float distance(const Coor<float> & c1_,
		 const Coor<float> & c2_){
    int atoms = c1_.size();
    Coor<double> c1;
    c1.resize(atoms);
    Coor<double> c2;
    c2.resize(atoms);
    for(int a=0;a<atoms;a++){
      int apos = a*Coor<double>::DIM;
      c1.coor[apos] = c1_.coor[apos];
      c1.coor[apos+1] = c1_.coor[apos+1];
      c1.coor[apos+2] = c1_.coor[apos+2];
    }
    for(int a=0;a<atoms;a++){
      int apos = a*Coor<double>::DIM;
      c2.coor[apos] = c2_.coor[apos];
      c2.coor[apos+1] = c2_.coor[apos+1];
      c2.coor[apos+2] = c2_.coor[apos+2];
    }
    ms.apply(c1,c2);      
    return ms.score();
  }
    
};


template<>
inline string to_string(const DCDTrj & dcd){
  return dcd.get_comment();
}

template<>
inline string to_string(const MaxSubClusterDistance & dcd){
  return "MaxSubDistance";
}

ostream & operator<<(ostream & out,const Coor<float> & ){
  out<<"<Coor>";
  return out;
}


template<typename float_type>
inline string to_string(const Coor<float_type> & data){
  return "Coor";
}


template<> 
inline string to_string(const ContactGroup &gr){
  char buff[256];
  sprintf(buff," %i",gr.size());
  return gr.name()+buff;
}

template<> 
inline string to_string(const Contact &c){
  return c.name();
}


template<>
inline string to_string(const NameSelection & ns){
  return to_string(ns.size());
}

template<>
inline string to_string(const DCDMonitorSet &){
  return "monitor_set";
}

template<typename ConstType>
void add_monitor( DCDMonitorSet & c, ConstType & cg){
  c.add_monitor(cg);
}

// void go_energy(const Molecules molecules, const DCDTrj & dcd,
// 	       string file,
// 	       const DefGoOptions & options, const MDB & mdb){
//   MOD_SWITCH = ENERGY;
//   set_prompt();

//     DefGoForceField *ff;

//   bool native_map = (options.native_map!="none");
//   bool const_file = (options.const_file!="none");


//   if(native_map&&const_file)
//     ff = new DefGoForceField(options.native_map,
// 				     molecules,
// 				     options,
// 				     options.const_file,
// 				     mdb);
//   else if(native_map&&(!const_file))
//     ff = new DefGoForceField(options.native_map,
// 				     molecules,
// 				     options,
// 				     mdb);
//   else if((!native_map)&&const_file)
//     ff = new DefGoForceField(molecules,
// 				     options,
// 				     options.const_file,
// 				     mdb);
//   else if((!native_map)&&(!const_file))
//     ff = new DefGoForceField(molecules,
// 				     options,
// 				     mdb);   

//   int natoms = dcd.get_natoms();
//   if(natoms!=molecules.atom_size()){
//     delete ff;
//     throw "Molecules and DCD do not contain the same number of atoms";
//   }

//   //coordinates
//   vector<double *> atomCoor;
//   for(int i=0; i<molecules.atom_size();i++){
//     double * nCoor = new double[3];
//     for(int j=0;j<3;j++)
//       nCoor[j]=molecules[i][j];
//     atomCoor.push_back(nCoor);
//   }


//   string res;
//   char buff[4048];
//   DefGoData data;
//   ostream * out;
//   if(file!="stdout"){
//     ofstream * fout = new ofstream();
//     fout->open(file.c_str());
//     out = fout;
//   }
//   else {
//     out = &cout;
//   }

//   (*out)<<"# frame energy native non-native Q const\n";
//   for(int i=0;i<dcd.get_nframes();i++){
//     for(int j=0;j<natoms;j++){
//       atomCoor[j][0]=dcd.X(i)[j];
//       atomCoor[j][1]=dcd.Y(i)[j];
//       atomCoor[j][2]=dcd.Z(i)[j];
//     }

//     ff->energy(atomCoor);
//     data = ff->data();
    

//     sprintf(buff,"%-10i  %-10.4f  %-10i  %-10i  %-10.4f %-10.4f\n",
// 	    i+1,
// 	    data.energy,
// 	    data.n_native,
// 	    data.n_non_native,
// 	    ((double)data.n_native)/((double)data.n_cont_max),
// 	    data.const_energy
// 	    );
//     (*out)<<buff;
//   }

//   //delete coor
//   for(int i=0; i<atomCoor.size();i++){
//     delete [] atomCoor[i];
//   }
//   delete ff;
//   out->flush();
//   if(file!="stdout"){delete out;}
// }


// void vdw_go_energy(const Molecules molecules, 
// 		   const DCDTrj & dcd,
// 		   string file,
// 		   const DefIGoOptions & options, const MDB & mdb){
//   MOD_SWITCH = ENERGY;
//   set_prompt();

//   DefIGoForceField *ff;

//   bool native_map = (options.native_map!="none");
//   bool const_file = (options.const_file!="none");


//   if(native_map&&const_file)
//     ff = new DefIGoForceField(options.native_map,
// 				     molecules,
// 				     options,
// 				     options.const_file,
// 				     mdb);
//   else if(native_map&&(!const_file))
//     ff = new DefIGoForceField(options.native_map,
// 				     molecules,
// 				     options,
// 				     mdb);
//   else if((!native_map)&&const_file)
//     ff = new DefIGoForceField(molecules,
// 				     options,
// 				     options.const_file,
// 				     mdb);
//   else if((!native_map)&&(!const_file))
//     ff = new DefIGoForceField(molecules,
// 				     options,
// 				     mdb);   

//   int natoms = dcd.get_natoms();
//   if(natoms!=molecules.atom_size()){
//     delete ff;
//     throw "Molecules and DCD do not contain the same number of atoms";
//   }

//   //coordinates
//   vector<double *> atomCoor;
//   for(int i=0; i<molecules.atom_size();i++){
//     double * nCoor = new double[3];
//     for(int j=0;j<3;j++)
//       nCoor[j]=molecules[i][j];
//     atomCoor.push_back(nCoor);
//   }


//   string res;
//   char buff[4048];
//   DefIGoData data;
//   ostream * out;
//   if(file!="stdout"){
//     ofstream * fout = new ofstream();
//     fout->open(file.c_str());
//     out = fout;
//   }
//   else {
//     out = &cout;
//   }

//   (*out)<<"# frame energy native non-native Q const\n";
//   for(int i=0;i<dcd.get_nframes();i++){
//     for(int j=0;j<natoms;j++){
//       atomCoor[j][0]=dcd.X(i)[j];
//       atomCoor[j][1]=dcd.Y(i)[j];
//       atomCoor[j][2]=dcd.Z(i)[j];
//     }

//     ff->energy(atomCoor);
//     data = ff->data();
    

//     sprintf(buff,"%-10i  %-10.4f  %-10i  %-10i  %-10.4f %-10.4f\n",
// 	    i+1,
// 	    data.energy,
// 	    data.n_native,
// 	    data.n_non_native,
// 	    ((double)data.n_native)/((double)data.n_cont_max),
// 	    data.const_energy
// 	    );
//     (*out)<<buff;
//   }

//   //delete coor
//   for(int i=0; i<atomCoor.size();i++){
//     delete [] atomCoor[i];
//   }
//   delete ff;
//   out->flush();
//   if(file!="stdout"){delete out;}
// }


// void c19_energy(const Molecules molecules, 
// 		const DCDTrj & dcd,
// 		string file,
// 		const DefC19Options & options, const MDB & mdb){
//   MOD_SWITCH = ENERGY;
//   set_prompt();

//   DefC19ForceField ff(molecules,options,mdb);;


//   int natoms = dcd.get_natoms();
//   if(natoms!=molecules.atom_size()){
//     throw "Molecules and DCD do not contain the same number of atoms";
//   }

//   //coordinates
//   vector<double *> atomCoor;
//   for(int i=0; i<molecules.atom_size();i++){
//     double * nCoor = new double[3];
//     for(int j=0;j<3;j++)
//       nCoor[j]=molecules[i][j];
//     atomCoor.push_back(nCoor);
//   }


//   string res;
//   char buff[4048];
//   DefC19Data data;
//   ostream * out;
//   if(file!="stdout"){
//     ofstream * fout = new ofstream();
//     fout->open(file.c_str());
//     out = fout;
//   }
//   else {
//     out = &cout;
//   }

//   (*out)<<"# frame energy bond angle dihe imph elec vdw sasa const pairs\n";
//   for(int i=0;i<dcd.get_nframes();i++){
//     for(int j=0;j<natoms;j++){
//       atomCoor[j][0]=dcd.X(i)[j];
//       atomCoor[j][1]=dcd.Y(i)[j];
//       atomCoor[j][2]=dcd.Z(i)[j];
//     }

//     ff.energy(atomCoor);
//     data = ff.data();
    

//     sprintf(buff,"%-10i  %-10.4f  %-10.4f  %-10.4f  %-10.4f  %-10.4f  %-10.4f  %-10.4f  %-10.4f  %-10.4f  %-10i\n",
// 	    i+1,
// 	    data.energy,
// 	    data.bond_energy,
// 	    data.angle_energy,
// 	    data.dihe_energy,
// 	    data.imph_energy,
// 	    data.elec_energy,
// 	    data.vdw_energy,
// 	    data.sasa_energy,  
// 	    data.const_energy,
// 	    data.pairs
// 	    );
//     (*out)<<buff;
//   }

//   //delete coor
//   for(int i=0; i<atomCoor.size();i++){
//     delete [] atomCoor[i];
//   }
//   out->flush();
//   if(file!="stdout"){delete out;}
// }


// void c22_energy(const Molecules molecules, 
// 		const DCDTrj & dcd,
// 		string file,
// 		const DefC22Options & options, const MDB & mdb){
//   MOD_SWITCH = ENERGY;
//   set_prompt();

//   DefC22ForceField ff(molecules,options,mdb);;


//   int natoms = dcd.get_natoms();
//   if(natoms!=molecules.atom_size()){
//     throw "Molecules and DCD do not contain the same number of atoms";
//   }

//   //coordinates
//   vector<double *> atomCoor;
//   for(int i=0; i<molecules.atom_size();i++){
//     double * nCoor = new double[3];
//     for(int j=0;j<3;j++)
//       nCoor[j]=molecules[i][j];
//     atomCoor.push_back(nCoor);
//   }


//   string res;
//   char buff[4048];
//   DefC22Data data;
//   ostream * out;
//   if(file!="stdout"){
//     ofstream * fout = new ofstream();
//     fout->open(file.c_str());
//     out = fout;
//   }
//   else {
//     out = &cout;
//   }

//   (*out)<<"# frame energy bond angle dihe imph elec vdw sasa const pairs\n";
//   for(int i=0;i<dcd.get_nframes();i++){
//     for(int j=0;j<natoms;j++){
//       atomCoor[j][0]=dcd.X(i)[j];
//       atomCoor[j][1]=dcd.Y(i)[j];
//       atomCoor[j][2]=dcd.Z(i)[j];
//     }

//     ff.energy(atomCoor);
//     data = ff.data();
    

//     sprintf(buff,"%-10i  %-10.4f  %-10.4f  %-10.4f  %-10.4f  %-10.4f  %-10.4f  %-10.4f  %-10.4f  %-10.4f  %-10i\n",
// 	    i+1,
// 	    data.energy,
// 	    data.bond_energy,
// 	    data.angle_energy,
// 	    data.dihe_energy,
// 	    data.imph_energy,
// 	    data.elec_energy,
// 	    data.vdw_energy,
// 	    data.ub_energy,  
// 	    data.const_energy,
// 	    data.pairs
// 	    );
//     (*out)<<buff;
//   }

//   //delete coor
//   for(int i=0; i<atomCoor.size();i++){
//     delete [] atomCoor[i];
//   }
//   out->flush();
//   if(file!="stdout"){delete out;}
// }

// void pbcc19_energy(const Molecules molecules, 
// 		   const DCDTrj & dcd,
// 		   string file,
// 		   const PBCC19Options & options, const MDB & mdb){
//   MOD_SWITCH = ENERGY;
//   set_prompt();

//   PBCC19ForceField ff(molecules,options,mdb);;


//   int natoms = dcd.get_natoms();
//   if(natoms!=molecules.atom_size()){
//     throw "Molecules and DCD do not contain the same number of atoms";
//   }

//   //coordinates
//   vector<double *> atomCoor;
//   for(int i=0; i<molecules.atom_size();i++){
//     double * nCoor = new double[3];
//     for(int j=0;j<3;j++)
//       nCoor[j]=molecules[i][j];
//     atomCoor.push_back(nCoor);
//   }


//   string res;
//   char buff[4048];
//   PBCC19Data data;
//   ostream * out;
//   if(file!="stdout"){
//     ofstream * fout = new ofstream();
//     fout->open(file.c_str());
//     out = fout;
//   }
//   else {
//     out = &cout;
//   }

//   (*out)<<"# frame energy bond angle dihe imph elec vdw sasa const pairs\n";
//   for(int i=0;i<dcd.get_nframes();i++){
//     for(int j=0;j<natoms;j++){
//       atomCoor[j][0]=dcd.X(i)[j];
//       atomCoor[j][1]=dcd.Y(i)[j];
//       atomCoor[j][2]=dcd.Z(i)[j];
//     }

//     ff.energy(atomCoor);
//     data = ff.data();
    

//     sprintf(buff,"%-10i  %-10.4f  %-10.4f  %-10.4f  %-10.4f  %-10.4f  %-10.4f  %-10.4f  %-10.4f  %-10.4f  %-10i\n",
// 	    i+1,
// 	    data.energy,
// 	    data.bond_energy,
// 	    data.angle_energy,
// 	    data.dihe_energy,
// 	    data.imph_energy,
// 	    data.elec_energy,
// 	    data.vdw_energy,
// 	    data.sasa_energy,  
// 	    data.const_energy,
// 	    data.pairs
// 	    );
//     (*out)<<buff;
//   }

//   //delete coor
//   for(int i=0; i<atomCoor.size();i++){
//     delete [] atomCoor[i];
//   }
//   out->flush();
//   if(file!="stdout"){delete out;}
// }


// void pbcc22_energy(const Molecules molecules, 
// 		   const DCDTrj & dcd,
// 		   string file,
// 		   const PBCC22Options & options, const MDB & mdb){
//   MOD_SWITCH = ENERGY;
//   set_prompt();

//   PBCC22ForceField ff(molecules,options,mdb);;


//   int natoms = dcd.get_natoms();
//   if(natoms!=molecules.atom_size()){
//     throw "Molecules and DCD do not contain the same number of atoms";
//   }

//   //coordinates
//   vector<double *> atomCoor;
//   for(int i=0; i<molecules.atom_size();i++){
//     double * nCoor = new double[3];
//     for(int j=0;j<3;j++)
//       nCoor[j]=molecules[i][j];
//     atomCoor.push_back(nCoor);
//   }


//   string res;
//   char buff[4048];
//   PBCC22Data data;
//   ostream * out;
//   if(file!="stdout"){
//     ofstream * fout = new ofstream();
//     fout->open(file.c_str());
//     out = fout;
//   }
//   else {
//     out = &cout;
//   }

//   (*out)<<"# frame energy bond angle dihe imph elec vdw sasa const pairs\n";
//   for(int i=0;i<dcd.get_nframes();i++){
//     for(int j=0;j<natoms;j++){
//       atomCoor[j][0]=dcd.X(i)[j];
//       atomCoor[j][1]=dcd.Y(i)[j];
//       atomCoor[j][2]=dcd.Z(i)[j];
//     }

//     ff.energy(atomCoor);
//     data = ff.data();
    

//     sprintf(buff,"%-10i  %-10.4f  %-10.4f  %-10.4f  %-10.4f  %-10.4f  %-10.4f  %-10.4f  %-10.4f  %-10.4f  %-10i\n",
// 	    i+1,
// 	    data.energy,
// 	    data.bond_energy,
// 	    data.angle_energy,
// 	    data.dihe_energy,
// 	    data.imph_energy,
// 	    data.elec_energy,
// 	    data.vdw_energy,
// 	    data.ub_energy,  
// 	    data.const_energy,
// 	    data.pairs
// 	    );
//     (*out)<<buff;
//   }

//   //delete coor
//   for(int i=0; i<atomCoor.size();i++){
//     delete [] atomCoor[i];
//   }
//   out->flush();
//   if(file!="stdout"){delete out;}
// }

template<>
inline string to_string(const Molecules & molecules){
  return "<Molecules>";
}

template<>
inline string to_string(const MDB & mdb){
  return "<MDB>";
}

template<>
inline string to_string(const RMSD & rmsd){
  return "<RMSD>";
}

template<>
inline string to_string(const Fit & rmsd){
  return "<Fit>";
}

template<>
inline string to_string(const Energy & c19_e){
  return "<Energy>";
}


template<>
inline string to_string(const RGYR & rmsd){
  return "<RGYR>";
}

template<>
inline string to_string(const Bond & rmsd){
  return "<Bond>";
}

template<>
inline string to_string(const Angle & rmsd){
  return "<Angle>";
}

template<>
inline string to_string(const Dihedral & rmsd){
  return "<Dihedral>";
}

template<>
inline string to_string(const Cluster & rmsd){
  return "<cluster>";
}

template<>
inline string to_string(const HBondFinderOpts & rmsd){
  return "<HBondfinderOpts>";
}

extern "C" {
  void init_dcd(){
    Module mod = Module("dcd","Trajectory Tools")
      .def_const("begin",1)
      .def_const("eof",-1)
      .def_function("dcd2pdb","Write coor to pdb bundle.",&dcd2pdb)
      .def_function("hbond_finder","Find HBonds between proteins and wather.",&hbond_finder)
      //       .def_function("go_energy",go_energy)
      //       .def_function("vdw_go_energy",vdw_go_energy)
      //       .def_function("c19_energy",c19_energy)
      //       .def_function("c22_energy",c22_energy)
      //       .def_function("pbc_c19_energy",pbcc19_energy)
      //       .def_function("pbc_c22_energy",pbcc22_energy);
      ; 

    Class< HBondFinderOpts>(mod.self(),"hbond_finder_options")
      .def_attribute("X",&HBondFinderOpts::X,"Undocumented")
      .def_attribute("Y",&HBondFinderOpts::Y,"Undocumented")
      .def_attribute("Z",&HBondFinderOpts::Z,"Undocumented")
      .def_attribute("dist",&HBondFinderOpts::d,"Undocumented")
      .def_attribute("angle",&HBondFinderOpts::angle,"Undocumented")
      ;


    Class<DCDTrj>(mod.self(),"dcd")
      .def_method("type",&DCDTrj::get_type)
      .def_method("add",&DCDTrj::add)
      .def_method("natoms",&DCDTrj::get_natoms)
      .def_method("nfix",&DCDTrj::get_nfix)
      .def_method("nfree",&DCDTrj::get_nfree)
      .def_method("version",&DCDTrj::get_version)
      .def_method("nframes",&DCDTrj::get_nframes)
      .def_method("first",&DCDTrj::get_first)
      .def_method("set_first",&DCDTrj::set_first)
      .def_method("nsteps",&DCDTrj::get_nsteps)
      .def_method("set_nsteps",&DCDTrj::set_nsteps)
      .def_method("dt",&DCDTrj::get_dt)
      .def_method("sex",&DCDTrj::get_sex)
      .def_method("comment",&DCDTrj::get_comment)
      .def_method("set_comment",&DCDTrj::set_comment)
      .def_method("from",&DCDTrj::get_from)
      .def_method("to",&DCDTrj::get_to)
      .def_method("step",&DCDTrj::get_step)
      .def_method("nsaved",&DCDTrj::get_nsaved)
      .def_method("read",&DCDTrj::read)
      .def_method("write",&DCDTrj::write)
      .def_method("append",&DCDTrj::append)
      .def_method("print_header",&DCDTrj::print_header)
      .def_method("print_save_info",&DCDTrj::print_save_info)
      .def_method("coor",&DCDTrj::coor)
      .def_method("X",&DCDTrj::X)
      .def_method("Y",&DCDTrj::Y)
      .def_method("Z",&DCDTrj::Z)
      ;

    Class<Contact,
      ArgList4<
    Molecules,string,string,double> 
      >(mod.self(),"contact")
      .def_method("eval",&Contact::eval);

    Class<ContactGroup,
      ArgList1<string> 
      >(mod.self(),"contact_group")
      .def_method("add",&ContactGroup::add)
      .def_method("eval",&ContactGroup::eval);

    Class<RMSD,
      ArgList2<Molecules,NameSelection> >(mod.self(),"RMSD")
      .def_method("eval",&RMSD::eval)
      .def_method("eval2",&RMSD::eval2)
      .def_method("eval_two",&RMSD::eval_two)
      ;

    Class<Energy,ArgList3<const Molecules & ,
      const EnergyOptions<double> &,
      const MDB &> >(mod.self(),"energy")
      .def_method("eval",&Energy::evalf)
      ;

    Class<Fit,
      ArgList2<Molecules,NameSelection> >(mod.self(),"Fit")
      .def_method("eval",&Fit::eval);

    Class<RGYR,
      ArgList2<Molecules,NameSelection> >(mod.self(),"RGYR")
      .def_method("eval",&RGYR::eval);

    Class<Bond,
      ArgList3<Molecules,string,string> >(mod.self(),"bond")
      .def_method("eval",&Bond::eval);

    Class<Angle,
      ArgList4<Molecules,string,string,string> >(mod.self(),"angle")
      .def_method("eval",&Angle::eval);

    Class<Dihedral,
      ArgList5<Molecules,string,string,string,string> >(mod.self(),"dihedral")
      .def_method("eval",&Dihedral::eval);


    Class<DCDMonitorSet>(mod.self(),"monitor_set")
      .def_method("clear_monitors",&DCDMonitorSet::clear_monitors)
      .def_method("size",&DCDMonitorSet::size)
      .def_method("apply",&DCDMonitorSet::apply)
      .def_method("export",&DCDMonitorSet::export_ds)
      ;


    mod.def_function("add_rmsd","Adds an RMSD monitor to the monitor_set",add_monitor<RMSD>);
    mod.def_function("add_rgyr","Adds a RGYR monitor to the monitor_set",add_monitor<RGYR>);
    mod.def_function("add_bond","Adds a Bond monitor to the monitor_set",add_monitor<Bond>);
    mod.def_function("add_angle","Adds a Angle monitor to the monitor_set",add_monitor<Angle>);
    mod.def_function("add_dihedral","Adds a Dihedral monitor to the monitor_set",add_monitor<Dihedral>);
    //     Class<Energy,ArgList3<Molecules,EnergyOptions<double>,MDB> >(mod.self(),
    // 								 "c19_energy")
    //       .def_method("eval",&Energy::eval);

    Class<Cluster,
      ArgList2<Molecules,NameSelection> >(mod.self(),"cluster")
      .def_method("leader",&Cluster::leader)
      .def_method("leader_maxsub",&Cluster::leader_with_distance<MaxSubClusterDistance>)
      .def_method("graph",&Cluster::graph)
      .def_method("neighbors",&Cluster::neighbors);

    Class<MaxSubClusterDistance,
      ArgList3<Molecules,float,float> >(mod.self(),"maxsub_dist");

  }

}