/***************************************************************************
                          blast_mod.cpp  -  description
                             -------------------
    begin                : Tue Jun 27 08:12:04 2006
    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 <blast.h>
#include <tinyxml/tinyxml.h>
#include <stack>
#include <psipred/psipred.h>
using namespace Almost;
template<>
inline string to_string(const Blast &){
  return "<blast>";
}

template<>
inline string to_string(const Blast::Hsp &){
  return "<blast::hsp>";
}

template<>
inline string to_string(const Psipred &){
  return "<psipred>";
}


template< typename Type>
string str_cast(Type t){
  stringstream s;
  s<<t;
  return s.str();
}

void Blast::parse(string out_file){
  TiXmlBase::SetCondenseWhiteSpace(false);
  TiXmlDocument doc(out_file);
  int status = doc.LoadFile();
  cout<<"LOAD STATUS "<<status<<endl;
  const TiXmlElement * root = doc.RootElement();
  if(root){
    cout<<root->Value()<<endl;
    const TiXmlElement * child;
    {
      child = root->FirstChildElement("BlastOutput_program");
      if(child){
	cout<<child->Value()<<" "<<child->GetText()<<endl;
      }
    }
    {
      child = root->FirstChildElement("BlastOutput_version");
      if(child){
	cout<<child->Value()<<" "<<child->GetText()<<endl;
      }
    }
    {
      child = root->FirstChildElement("BlastOutput_reference");
      if(child){
	cout<<child->Value()<<" "<<child->GetText()<<endl;
      }
    }
    {
      child = root->FirstChildElement("BlastOutput_db");
      if(child){
	cout<<child->Value()<<" "<<child->GetText()<<endl;
      }
    }
    {
      child = root->FirstChildElement("BlastOutput_query-ID");
      if(child){
	query_.id = child->GetText();
      }
    }
    {
      child = root->FirstChildElement("BlastOutput_query-def");
      if(child){
	query_.def = child->GetText();
      }
    }
    {
      child = root->FirstChildElement("BlastOutput_query-len");
      if(child){
	query_.len = atoi(child->GetText());
      }
    }
    {
      child = root->FirstChildElement("BlastOutput_param");
      if(child){
	parse_blast_param(child->FirstChildElement());
      }
    }
    {
      child = root->FirstChildElement("BlastOutput_iterations");
      if(child){
	parse_blast_iterations(child);
      }
    }
  }  
}


void Blast::parse_blast_iterations(const TiXmlElement * iters){
  const TiXmlElement * iter;
  for(iter = iters->FirstChildElement("Iteration");iter;
      iter = iter->NextSiblingElement ("Iteration")){
    //     cerr<<"DEBUG process iteration "; 
    hits.resize(hits.size()+1);
    //    cerr<<hits.size()<<endl;
    parse_blast_iteration(iter);
  }
}

void Blast::parse_blast_iteration(const TiXmlElement * iter){
  const TiXmlElement * iter_hits = 
    iter->FirstChildElement("Iteration_hits");
  if(!iter_hits) { cout<<"NO HITS "<<endl; return;}
  
  const TiXmlElement * hit;
  for(hit = iter_hits->FirstChildElement("Hit");hit;
      hit = hit->NextSiblingElement ("Hit")){
    //    cout<<"DEBUG process Hit\n";
    Hit h;
    h.iteration = hits.size();
    {
      const TiXmlElement * res = hit->FirstChildElement("Hit_num");
      if(res){
	h.num = atoi(res->GetText());
	//	cout<<"DEBUG num "<<h.num<<endl;
      }
      else { cout<<"ERRROR num missing ";}
    }
    {
      const TiXmlElement * res = hit->FirstChildElement("Hit_id");
      if(res){
	h.id = res->GetText();
	//	cout<<"DEBUG id "<<h.id<<endl;
      }
      else { cout<<"ERRROR id missing ";}
    }
    {
      const TiXmlElement * res = hit->FirstChildElement("Hit_def");
      if(res){
	h.def = res->GetText();
	//	cout<<"DEBUG def "<<h.def<<endl;
      }
      else { cout<<"ERRROR def missing ";}
    }
    {
      const TiXmlElement * res = hit->FirstChildElement("Hit_accession");
      if(res){
	h.accession = res->GetText();
	//	cout<<"DEBUG accession "<<h.accession<<endl;
      }
      else { cout<<"ERRROR accession missing ";}
    }
    {
      const TiXmlElement * res = hit->FirstChildElement("Hit_len");
      if(res){
	h.len = atoi(res->GetText());
	//	cout<<"DEBUG len "<<h.len<<endl;
      }
      else { cout<<"ERRROR len missing ";}
    }

    //HSPS
    const TiXmlElement * hsp, * hsps;
    hsps = hit->FirstChildElement("Hit_hsps");
    for(hsp = hsps->FirstChildElement("Hsp"); hsp;
	hsp = hsp->NextSiblingElement ("Hsp")){
      Hsp hsp_;
      parse_blast_hsp(hsp,hsp_);
      h.hsp.push_back(hsp_);
    }
    hits.back().push_back(h);
  }

}
void Blast::parse_blast_hsp(const TiXmlElement * hsp,Hsp & hsp_){
  {
    const TiXmlElement * res = hsp->FirstChildElement("Hsp_num");
    if(res){
      hsp_.num = atoi(res->GetText());
      //      cout<<"DEBUG len "<<hsp_.num<<endl;
    }
    else { cout<<"ERRROR len missing ";}
  }
  {
    const TiXmlElement * res = hsp->FirstChildElement("Hsp_bit-score");
    if(res){
      hsp_.bit_score = atof(res->GetText());
      //      cout<<"DEBUG len "<<hsp_.bit_score<<endl;
    }
    else { cout<<"ERRROR len missing ";}
  }
  {
    const TiXmlElement * res = hsp->FirstChildElement("Hsp_score");
    if(res){
      hsp_.score = atof(res->GetText());
      //      cout<<"DEBUG len "<<hsp_.score<<endl;
    }
    else { cout<<"ERRROR len missing ";}
  }
  {
    const TiXmlElement * res = hsp->FirstChildElement("Hsp_evalue");
    if(res){
      hsp_.evalue = atof(res->GetText());
      //      cout<<"DEBUG len "<<hsp_.evalue<<endl;
    }
    else { cout<<"ERRROR len missing ";}
  }
  {
    const TiXmlElement * res = hsp->FirstChildElement("Hsp_query-from");
    if(res){
      hsp_.query_from = atoi(res->GetText());
      //      cout<<"DEBUG len "<<hsp_.query_from<<endl;
    }
    else { cout<<"ERRROR len missing ";}
  }
  {
    const TiXmlElement * res = hsp->FirstChildElement("Hsp_query-to");
    if(res){
      hsp_.query_to = atoi(res->GetText());
      //      cout<<"DEBUG len "<<hsp_.query_to<<endl;
    }
    else { cout<<"ERRROR len missing ";}
  }
  {
    const TiXmlElement * res = hsp->FirstChildElement("Hsp_hit-from");
    if(res){
      hsp_.hit_from = atoi(res->GetText());
      //      cout<<"DEBUG len "<<hsp_.hit_from<<endl;
    }
    else { cout<<"ERRROR len missing ";}
  }
  {
    const TiXmlElement * res = hsp->FirstChildElement("Hsp_hit-to");
    if(res){
      hsp_.hit_to = atoi(res->GetText());
      //      cout<<"DEBUG len "<<hsp_.hit_to<<endl;
    }
    else { cout<<"ERRROR len missing ";}
  }
  {
    const TiXmlElement * res = hsp->FirstChildElement("Hsp_query-frame");
    if(res){
      hsp_.query_frame = atoi(res->GetText());
      //      cout<<"DEBUG len "<<hsp_.query_frame<<endl;
    }
    else { cout<<"ERRROR len missing ";}
  }
  {
    const TiXmlElement * res = hsp->FirstChildElement("Hsp_hit-frame");
    if(res){
      hsp_.hit_frame = atoi(res->GetText());
      //      cout<<"DEBUG len "<<hsp_.hit_frame<<endl;
    }
    else { cout<<"ERRROR len missing ";}
  }
  {
    const TiXmlElement * res = hsp->FirstChildElement("Hsp_identity");
    if(res){
      hsp_.identity = atoi(res->GetText());
      //      cout<<"DEBUG len "<<hsp_.identity<<endl;
    }
    else { cout<<"ERRROR len missing ";}
  }
  {
    const TiXmlElement * res = hsp->FirstChildElement("Hsp_positive");
    if(res){
      hsp_.positive = atoi(res->GetText());
      //      cout<<"DEBUG len "<<hsp_.positive<<endl;
    }
    else { cout<<"ERRROR len missing ";}
  }
  {
    const TiXmlElement * res = hsp->FirstChildElement("Hsp_align-len");
    if(res){
      hsp_.align_len = atoi(res->GetText());
      //      cout<<"DEBUG len "<<hsp_.align_len<<endl;
    }
    else { cout<<"ERRROR len missing ";}
  }
  {
    const TiXmlElement * res = hsp->FirstChildElement("Hsp_qseq");
    if(res){
      hsp_.qseq = res->GetText();
      //      cout<<"DEBUG len "<<hsp_.qseq<<endl;
    }
    else { cout<<"ERRROR len missing ";}
  }
  {
    const TiXmlElement * res = hsp->FirstChildElement("Hsp_hseq");
    if(res){
      hsp_.hseq = res->GetText();
      //      cout<<"DEBUG len "<<hsp_.hseq<<endl;
    }
    else { cout<<"ERRROR len missing ";}
  }
  {
    const TiXmlElement * res = hsp->FirstChildElement("Hsp_midline");
    if(res){
      hsp_.midline = res->GetText();
      //      cout<<"DEBUG len "<<hsp_.midline<<endl;
    }
    else { cout<<"ERRROR len missing ";}
  }

}

class AlmShaper {
  struct BondRestraint {
    int atom1;
    int atom2;
    int sep;
    double d;
  };

  vector<BondRestraint> restraints;
  set<int> ignore;
public:
  AlmShaper(const Molecules & trg, const Protein & ref, const Blast::Hsp & hsp){
    mkhomorest(ref,hsp,trg);
  }

  int max_sep(){
    int s=0;
    for(int i=0;i<restraints.size();i++){
      if(restraints[i].sep>s) s = restraints[i].sep;
    }
    return s;
  }

  BondConst<double> mkconst(const Molecules & mols, int m_sep, double K){
    BondConst<double> cst(mols);
    for(int i=0;i<restraints.size();i++){
      if(restraints[i].sep>m_sep) continue;
      int a1 = restraints[i].atom1;
      int a2 = restraints[i].atom2;
      int f1 = mols.find_fragment(a1);
      int f2 = mols.find_fragment(a2);
      if((ignore.find(f1)!=ignore.end())||
	 (ignore.find(f2)!=ignore.end())){
	aout<<">> Ignoring rest: "
	    <<" "<<f1+1<<" "<<mols.atom_name(a1,Molecules::BASE)<<" "
	    <<" "<<f2+1<<" "<<mols.atom_name(a2,Molecules::BASE)<<endl;
	continue;
      }

      cst.add(restraints[i].atom1,
	      restraints[i].atom2,
	      restraints[i].d,
	      K);
    }
    return cst;
  }

  void add_ignore(int i){
    ignore.insert(i-1);
  }
private:
  void mkhomorest( const Protein & ref, const Blast::Hsp & hsp, const Molecules & trg){
    int h_from = hsp.hit_from; //h_from = h_from -1;
    int h_to = hsp.hit_to;

    int q_from = hsp.query_from; //q_from = q_from -1;
    int q_to = hsp.query_to;

    
    int igap = 0;
    int iigap = 0;
    for(int i=0;i<hsp.qseq.size();i++){
      int cont = false;
      if(hsp.hseq[i]=='-'){ ++igap; cont = true;}
      if(hsp.qseq[i]=='-'){ ++iigap; cont = true;}
      if(cont) continue;
      int ca_1;
      
      stringstream str;
      str<<one2three(hsp.hseq[i])<<"_"<<i+h_from-igap<<"/CA";
      ca_1 =  ref.find_atom(str.str());
      if(ca_1==-1){
	cerr<< string("Unable to find atom: "+str.str())<<endl;
	continue;
      }
      int jgap = igap;
      int jjgap = iigap;
      for(int j=i+1;j<hsp.qseq.size();j++){
	bool cont2 = false;
	if(hsp.hseq[j]=='-'){ ++jgap; cont2 = true;}
	if(hsp.qseq[j]=='-'){ ++jjgap; cont2 = true;}
	if(cont2) continue;
	
	int ca_2;
    
	stringstream str2;
	str2<<one2three(hsp.hseq[j])<<"_"<<j+h_from-jgap<<"/CA";
	ca_2 =  ref.find_atom(str2.str());
	if(ca_2==-1){      
	  cerr<< string("Unable to find atom: "+str2.str())<<endl;
	  continue;
	}

	{
	  stringstream atm1;
	  atm1<<one2three(hsp.qseq[i+q_from-iigap-1])<<"_"<<i+q_from-iigap<<"/CA";
	  stringstream atm2;
	  atm2<<one2three(hsp.qseq[j+q_from-jjgap-1])<<"_"<<j+q_from-jjgap<<"/CA";
	  double dist = ref.distance(ca_1,ca_2);
	  int sep = abs(i+q_from-iigap-1-(j+q_from-jjgap-1));

	  if(dist<10&&sep>2)
	    addRest(atm1.str(),
		    atm2.str(),
		    ref.distance(ca_1,ca_2),
		    abs(i+q_from-iigap-1-(j+q_from-jjgap-1)),
		    trg);
	}
      }
    }
  }

  void addRest(string a1, string a2, double d, int sep, const Molecules & mols){
    //    cout<<"DEB addinmg "<<a1<<" "<<a2<<endl;
    NameSelection s1(mols,"//"+a1);
    NameSelection s2(mols,"//"+a2);
    
    if(s1.result().size()==1&&s2.result().size()==1){
      BondRestraint r;
      r.atom1 = *s1.result().begin();
      r.atom2 = *s2.result().begin();
      r.d = d;
      r.sep = sep;
      restraints.push_back(r);
    }
  }
};


void mkhomorest(Protein & ref, const Blast::Hsp & hsp, string file){
  int h_from = hsp.hit_from; //h_from = h_from -1;
  int h_to = hsp.hit_to;

  int q_from = hsp.query_from; //q_from = q_from -1;
  int q_to = hsp.query_to;

  ofstream out;
  out.open(file.c_str());

  int igap = 0;
  int iigap = 0;
  for(int i=0;i<hsp.qseq.size();i++){
    int cont = false;
    if(hsp.hseq[i]=='-'){ ++igap; cont = true;}
    if(hsp.qseq[i]=='-'){ ++iigap; cont = true;}
    if(cont) continue;
    int ca_1;

    stringstream str;
    str<<one2three(hsp.hseq[i])<<"_"<<i+h_from-igap<<"/CA";
    ca_1 =  ref.find_atom(str.str());
    if(ca_1==-1){
      cerr<< string("Unable to find atom: "+str.str())<<endl;
      continue;
    }
    int jgap = igap;
    int jjgap = iigap;
    for(int j=i+1;j<hsp.qseq.size();j++){
      bool cont2 = false;
      if(hsp.hseq[j]=='-'){ ++jgap; cont2 = true;}
      if(hsp.qseq[j]=='-'){ ++jjgap; cont2 = true;}
      if(cont2) continue;

      int ca_2;
    
      stringstream str2;
      str2<<one2three(hsp.hseq[j])<<"_"<<j+h_from-jgap<<"/CA";
      ca_2 =  ref.find_atom(str2.str());
      if(ca_2==-1){      
	cerr<< string("Unable to find atom: "+str2.str())<<endl;
	continue;
      }
      out<<one2three(hsp.qseq[i+q_from-iigap-1])<<"_"<<i+q_from-iigap<<" "<<"CA "<<" "
	 <<one2three(hsp.qseq[j+q_from-jjgap-1])<<"_"<<j+q_from-jjgap<<" "<<"CA "<<" "
	 <<ref.distance(ca_1,ca_2)<<" "<<j+q_from-jjgap-(i+q_from-iigap)<<"\n";
      
    }
  }
  
}

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

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

template<>
inline string to_string(const BondConst<double> &){
  return "<BondConst>";
}

template<>
inline string to_string(const AlmShaper &){
  return "<blast::shaper>";
}

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

    Module mod = Module("blast","Interface to ncbi blast");
    Class<Blast::Hsp>(mod.self(),"hsp");

    Class<Blast>(mod.self(),"blast")
      .def_method("set_exedir",&Blast::set_exedir)
      .def_method("exedir",&Blast::exedir)
      .def_method("set_dbdir",&Blast::set_dbdir)
      .def_method("dbdir",&Blast::dbdir)
      .def_method("set_db",&Blast::set_db)
      .def_method("db",&Blast::db)
      .def_method("set_pdbdir",&Blast::set_pdbdir)
      .def_method("pdbdir",&Blast::pdbdir)
      .def_method("blastpgp",&Blast::blastpgp)
      .def_method("parse",&Blast::parse)
      .def_method("extract_sequences",&Blast::extract_sequences)
      .def_method("hsp",&Blast::hsp)
      .def_attribute("A",&Blast::A,"Undocumented")
      .def_attribute("f",&Blast::f,"Undocumented")
      .def_attribute("e",&Blast::e,"Undocumented")
      .def_attribute("m",&Blast::m,"Undocumented")
      .def_attribute("o",&Blast::o,"Undocumented")
      .def_attribute("y",&Blast::y,"Undocumented")
      .def_attribute("P",&Blast::P,"Undocumented")
      .def_attribute("F",&Blast::F,"Undocumented")
      .def_attribute("G",&Blast::G,"Undocumented")
      .def_attribute("E",&Blast::E,"Undocumented")
      .def_attribute("X",&Blast::X,"Undocumented")
      .def_attribute("N",&Blast::N,"Undocumented")
      .def_attribute("S",&Blast::S,"Undocumented")
      .def_attribute("H",&Blast::H,"Undocumented")
      .def_attribute("a",&Blast::a,"Undocumented")
      .def_attribute("I",&Blast::I,"Undocumented")
      .def_attribute("h",&Blast::h,"Undocumented")
      .def_attribute("c",&Blast::c,"Undocumented")
      .def_attribute("j",&Blast::j,"Undocumented")
      .def_attribute("J",&Blast::J,"Undocumented")
      .def_attribute("Z",&Blast::Z,"Undocumented")
      .def_attribute("O",&Blast::O,"Undocumented")
      .def_attribute("M",&Blast::M,"Undocumented")
      .def_attribute("v",&Blast::v,"Undocumented")
      .def_attribute("b",&Blast::b,"Undocumented")
      .def_attribute("C",&Blast::C,"Undocumented")
      .def_attribute("R",&Blast::R,"Undocumented")
      .def_attribute("W",&Blast::W,"Undocumented")
      .def_attribute("z",&Blast::z,"Undocumented")
      .def_attribute("K",&Blast::K,"Undocumented")
      .def_attribute("s",&Blast::s,"Undocumented")
      .def_attribute("Y",&Blast::Y,"Undocumented")
      .def_attribute("p",&Blast::p,"Undocumented")
      .def_attribute("k",&Blast::k,"Undocumented")
      .def_attribute("T",&Blast::T,"Undocumented")
      .def_attribute("Q",&Blast::Q,"Undocumented")
      .def_attribute("B",&Blast::B,"Undocumented")
      .def_attribute("l",&Blast::l,"Undocumented")
      .def_attribute("U",&Blast::U,"Undocumented")
      ;

    Class<AlmShaper,
      ArgList3< const Molecules & , const Protein &  ,const Blast::Hsp & > >(mod.self(),"shaper")
      .def_method("max_sep",&AlmShaper::max_sep)
      .def_method("mkconst",&AlmShaper::mkconst)
      .def_method("add_ignore",&AlmShaper::add_ignore)
      ;

    Class<Psipred,ArgList1<string> >(mod.self(),"psipred")
      .def_method("run_pass1",&Psipred::run_pass1)
      .def_method("run_pass2",&Psipred::run_pass2)
      .def_method("run",&Psipred::run)
      ;

    mod.def_function("mkhomorest","",mkhomorest);
  }
}