/***************************************************************************
                          cheshire_mod.cpp  -  description
                             -------------------
    begin                : Mon Dec 17 11:23:57 2007
    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/tostring.h>
#include <zeta/module.h>
#include <zeta/class.h>
#include <zeta/zextclass.h>
#include <zeta/zobj.hpp>
#include <mdb.h>
#include <pdb.h>
#include <molecules/molecules.h>
#include <coor.h>
#include <secstruct.h>
#include <iterator>


using namespace Almost;
using namespace std;
extern Coor<double> mols_coor(const Molecules & molecules);
void cheshire_bootstrap_pdb(string name,
			    string shx_file, 
			    string pdb_file, 
			    string mdb_file){
  //Make fasta
  MDB mdb(mdb_file.c_str());
  PDB pdb(pdb_file.c_str());
  Protein p(name);
  p.build_missing(pdb[0][0],mdb,"DEFAULT","DEFAULT");
  vector<string> seq = p.sequence(Protein::SHORT);
  string fasta_seq;
  for(int i=0;i<seq.size();i++)
    fasta_seq += three2one(seq[i]);
  aout<<">> Sequence "<<fasta_seq<<endl;
  {
    ofstream fasta;
    fasta.open((name+".fasta").c_str());
    if(fasta){
      fasta<<">"<<name<<"\n";
      fasta<<fasta_seq<<"\n";
    }
    else {
      throw ">> FATAL ERROR: can't open "+name+".fasta";
    }
  }

  string sec;
  {
    Molecules m;
    m.add_protein(p);
    //Secondary structure
    AlmSecStruct ss(m);
    sec = ss.secstruct(mols_coor(m));
    aout<<">> Secondary structure: "<<sec<<endl;
  }

  vector<double> ha,h,n,ca,cb,c;
  {
    ifstream in;
    in.open(shx_file.c_str());
    if(!in) throw ">> FATAL ERROR: can't open "+shx_file;

    istream_iterator<string> iter(in),end;
    while(iter!=end){
      ++iter;++iter;
      ha.push_back(atof(iter->c_str()));++iter;
      h.push_back(atof(iter->c_str()));++iter;
      n.push_back(atof(iter->c_str()));++iter;
      ca.push_back(atof(iter->c_str()));++iter;
      cb.push_back(atof(iter->c_str()));++iter;
      c.push_back(atof(iter->c_str()));++iter;
    }
  }

  vector<double> phi = p.phi();
  vector<double> psi = p.psi();
  vector<double> omega = p.omega();

  ofstream mrg;
  mrg.open((name+".mrg").c_str());
  for(int i=0;i<seq.size();i++){
    mrg<<i+1<<" "<<fasta_seq[i]<<" "
       <<phi[i]*180.0/M_PI<<" "<<psi[i]*180.0/M_PI<<" "
       <<omega[i]*180.0/M_PI<<" "
       <<ha[i]<<" "
       <<h[i]<<" "
       <<n[i]<<" "
       <<ca[i]<<" "
       <<cb[i]<<" "
       <<c[i]<<" "
       <<sec[i]<<endl;
  }
}
 

void cheshire_bootstrap(string name,
			string shx_file, 
			string fasta_file){
  //Make fasta
  string fasta_seq;
  {
    ifstream f;
    f.open(fasta_file.c_str());
    if(!f) throw ">> FATAL ERROR: can't open "+fasta_file;
    istream_iterator<string> iter(f),end;
    while(iter!=end){
      if((*iter)[0]!='>') fasta_seq += *iter;
      ++iter;
    }
  }

  aout<<">> Sequence "<<fasta_seq<<endl;


  vector<double> ha,h,n,ca,cb,c;
  {
    ifstream in;
    in.open(shx_file.c_str());
    if(!in) throw ">> FATAL ERROR: can't open "+shx_file;

    istream_iterator<string> iter(in),end;
    while(iter!=end){
      ++iter;++iter;
      ha.push_back(atof(iter->c_str()));++iter;
      h.push_back(atof(iter->c_str()));++iter;
      n.push_back(atof(iter->c_str()));++iter;
      ca.push_back(atof(iter->c_str()));++iter;
      cb.push_back(atof(iter->c_str()));++iter;
      c.push_back(atof(iter->c_str()));++iter;
    }
  }


  ofstream mrg;
  mrg.open((name+".mrg").c_str());
  for(int i=0;i<fasta_seq.size();i++){
    mrg<<i+1<<" "<<fasta_seq[i]<<" "
       <<999<<" "<<999<<" "<<999<<" "
       <<ha[i]<<" "
       <<h[i]<<" "
       <<n[i]<<" "
       <<ca[i]<<" "
       <<cb[i]<<" "
       <<c[i]<<" "
       <<"X"<<endl;
  }
}
 

string che_random_path(map<string, double> & m, int nbest){
  map<double,string> mm;
  map<string,double>::iterator iter= m.begin(), end = m.end();
  while(iter!=end){
    mm[iter->second] = iter->first;
    ++iter;
  }
  int pos = lrand48()%(min(nbest,(int)mm.size()));
  {
    map<double,string>::iterator iter= mm.begin(), end = mm.end();
    int count = 0;
    while(iter!=end){
      if(count==pos) return iter->second;
      ++iter; ++count;
    }
    //?? Should never get here
    return mm.begin()->second;
  }
}

extern "C" {
  void init_cheshire(){
    //declarations here

    Module mod = Module("cheshire","Chemocal shift restraints");
    mod.def_function("bootstrap_pdb","",cheshire_bootstrap_pdb);
    mod.def_function("bootstrap","",cheshire_bootstrap);
    mod.def_function("best_random","",che_random_path);
  }
}