/***************************************************************************
                          sandy.cpp  -  description
                             -------------------
    begin                : Tue Nov 18 01:45:35 2003
    copyright            : (C) 1999-2008 by Cavalli Andrea
    author               : : cavalli $
    date                 : : 2003/12/09 11:10:33 $
    id                   : : coor.h,v 1.1.2.2 2003/12/09 11:10:33 cavalli Exp $
    email                : cavalli@bioc.unizh.ch amc82@cam.ac.uk
**************************************************************************/
/***************************************************************************
 *                                                                         *
 *   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 <almost.h>
#include <sstream>
#include <mdb.h>
#include <pdb.h>
#include <molecules/molecules.h>
#include <iterator>
#include <moveset.h>
namespace Almost {


  struct Peak2D {
    double w1;
    double w2;
    double Vol;
    double Vol2;
    vector<int> n1;
    vector<int> n2;
    vector<pair<int,int> > _12;
    vector<pair<int,int> > _13;
    vector<pair<int,int> > _14;
    
    int gmin;
    int gmax;

    vector<vector<double> > Upper;
    vector<vector<double> > Lower;

    vector<vector<double> > W;

    void init_w(){
      W.clear();
      W.resize(n1.size());

      Upper.clear();
      Upper.resize(n1.size());
      for(int i=0;i<n1.size();i++) Upper[i].resize(n2.size());

      Lower.clear();
      Lower.resize(n1.size());
      for(int i=0;i<n1.size();i++) Lower[i].resize(n2.size());

      for(int i=0;i<n1.size();i++) W[i].resize(n2.size());

      for(int i=0;i<W.size();i++)
	for(int j=0;j<W[i].size();j++)
	  W[i][j] = 1.0/(n1.size()*n2.size());
    }

    void move(const Protein & p){
      double W_tot = 0;
      for(int i=0;i<W.size();i++){
	for(int j=0;j<W[i].size();j++){
	  if(drand48()<0.05){
	    W[i][j] = W[i][j]+0.01*(drand48()-0.5);
	    if(W[i][j]<0) W[i][j] = 0;
	    if(p.distance(n1[i],n2[j])>6) W[i][j] = 0;
	  }
	  // 	  if(_12.size()||_13.size()||_14.size()){
	  // 	    bool a = is_12(n1[i],n2[j]);
	  // 	    bool b = is_13(n1[i],n2[j]);
	  // 	    bool c = is_14(n1[i],n2[j]);
	  // 	    if(!(a||b||c))
	  // 	      W[i][j] = 0;

	  // 	  }

	  W_tot += W[i][j];
	}
      }

      W_tot = 1.0/W_tot;
      for(int i=0;i<W.size();i++){
	for(int j=0;j<W[i].size();j++){
	  W[i][j] *= W_tot;
	}
      }
    }

    bool is_intra(const vector<int> & atm){
      for(int i=0;i<atm.size();i++){
	for(int j = i+1;j<atm.size();j++){
	  if(find(n1.begin(),n1.end(),atm[i])!=n1.end() &&
	     find(n2.begin(),n2.end(),atm[j])!=n2.end()) return true;

	  if(find(n1.begin(),n1.end(),atm[j])!=n1.end() &&
	     find(n2.begin(),n2.end(),atm[i])!=n2.end()) return true;



	}
      }
      return false;
    }

    bool is_14(int a1, int a2){
      pair<int,int> q(a1,a2);
      return find(_14.begin(),_14.end(),q)!=_14.end();
    }


    bool is_13(int a1, int a2){
      pair<int,int> q(a1,a2);
      return find(_13.begin(),_13.end(),q)!=_13.end();
    }

    bool is_12(int a1, int a2){
      pair<int,int> q(a1,a2);
      return find(_12.begin(),_12.end(),q)!=_12.end();
    }



    bool is_assignment(int a1, int a2, bool both = true){
      if(find(n1.begin(),n1.end(),a1)!=n1.end() &&
	 find(n2.begin(),n2.end(),a2)!=n2.end()) return true;
      if(both)
	if(find(n1.begin(),n1.end(),a2)!=n1.end() &&
	   find(n2.begin(),n2.end(),a1)!=n2.end()) return true;

      return false;
    }
  };

  Protein mkprot(vector<string> seq){
    MDB mdb("all19_eef1.mdb");
    Protein p("BLT");
    p.build_sequ(seq,mdb,"DEFAULT","DEFAULT");
    return p;
  }

  vector<Peak2D> read_peak(char * file){
    ifstream in(file);
    
    if(!in){
      string error = string("\t>> FATAL ERROR: unable to open file: ")+file
	+string("\n");
      cout<<error;
      throw error;
    }

    vector<Peak2D> peak;
    while(in.good()){
      char line[256];
      in.getline(line,256);
      cout<<line<<"\n";
      if(line[0]=='#'||in.gcount()<50) continue;

      stringstream ss;
      ss<<line;
      
      //
      int id;
      char c;
      Peak2D p;
      ss>>id;
      ss>> p.w1;
      ss>> p.w2;
      ss>>id;
      ss>> c;
      ss>> p.Vol;
      ss>> p.Vol2;
      peak.push_back(p);
    }

    return peak;
  }


  vector<string> read_seq(char * file){
    ifstream in(file);
    
    if(!in){
      string error = string("\t>> FATAL ERROR: unable to open file: ")+file
	+string("\n");
      cout<<error;
      throw error;
    }

    vector<string> seq;

    while(in.good()){
      char line[256];
      in.getline(line,256);
      stringstream ss(line);
      string s,n;
      ss>>s;
      if(s.size()<3) continue;
      ss>>n;
      n="_"+n;
      if(s.size()==3) seq.push_back(s+n);
      else if(s.size()>3) seq.push_back(string(s.begin(),s.begin()+3)+n);
      else {
	cout<<"\t>> Error in reading sequence file: "<<s<<"\n";
	throw "\t>> Error in reading sequence";
      }
    }

    return seq;
  }


  map<int, double> read_shifts(char * file, const Protein &p){
    map<int,double>        cs;
    ifstream in(file);
    
    if(!in){
      string error = string("\t>> FATAL ERROR: unable to open file: ")+file
	+string("\n");
      cout<<error;
      throw error;
    }

    vector<string> seq = p.sequence(Protein::LONG);

    while(in.good()){
      char line[256];
      in.getline(line,256);
      if(in.gcount()<10) continue;
      stringstream ss(line);

      int id;
      double w;
      double err;
      string atm;
      int    nres;
      ss>>id;
      ss>>w;
      ss>>err;
      ss>>atm;
      ss>>nres;
      --nres;

      cout<<"\t>> Chemical shift for atom "<<atm<<" in res "<<seq[nres]<<endl;
      int natm = p.find_atom(seq[nres]+"/"+atm);
      {
	//Hack for CD1 ILE...
	if(natm==-1&&atm=="CD1"){
	  atm = "CD";
	  natm = p.find_atom(seq[nres]+"/"+atm);
	}

      }
      if(natm==-1){
	cout<<"\t>> Atom not in protein\n";
      }

      if(natm!=-1){
	cs[natm] = w;
	cout<<"\t>> CS assignment "<<seq[nres]+"/"+atm<<" "<<w<<"\n";
      }

    }

    return cs;
  }

  void assign_peaks(vector<Peak2D> &  p1, const Protein & p,
		    const map<int,double> &cs,
		    double tol = 0.25){
    
    map<int,double>::const_iterator iter = cs.begin(), end = cs.end();
    while(iter!=end){
      for(int i=0;i<p1.size();i++){
	if(abs(p1[i].w1-iter->second)<tol) p1[i].n1.push_back(iter->first);
	if(abs(p1[i].w2-iter->second)<tol) p1[i].n2.push_back(iter->first);
      }
      ++iter;
    }

    for(int i=0;i<p1.size();i++){
      p1[i].gmax = -1;
      p1[i].gmin = 1000000;
      for(int a1=0;a1<p1[i].n1.size();a1++){
	int f1 = p.find_fragment(p1[i].n1[a1]);
	for(int a2=0;a2<p1[i].n2.size();a2++){
	  int f2 = p.find_fragment(p1[i].n2[a2]);
	  int q = abs(f1-f2);
	  if(q>p1[i].gmax) p1[i].gmax = q;
	  if(q<p1[i].gmin) p1[i].gmin = q;
	}
      }
    }
    {
      vector<Peak2D>::iterator iter= p1.begin(), end = p1.end();
      while(iter!=p1.end()){
	if(iter->n1.size()==0||iter->n2.size()==0)
	  iter = p1.erase(iter);
	else
	  ++iter;
      }
    }

    for(int i=0;i<p1.size();i++) p1[i].init_w();
  }



  void dump(vector<Almost::Peak2D> &p1, 
	    const Protein & pp ){
    ofstream out;
    out.open("dump.peaks");
    for(int i=0;i<p1.size();i++){
      //    if(p1[i].n1.size()>2||p1[i].n2.size()>2) continue;
      //    if(p1[i].gmax>10||p1[i].gmin<2) continue;
      out<<">> Peak "<<i+1<<" "//<<P_W[i]<<" "
	 <<p1[i].n1.size()<<" ";
      out<<p1[i].n2.size()<<" ";
      out<< "{ ";
      for(int a=0;a<p1[i].n1.size();a++) out<<p1[i].n1[a]<<" ";
      out<< "} ";    
      out<< "{ ";
      for(int a=0;a<p1[i].n2.size();a++) out<<p1[i].n2[a]<<" ";
      out<< "} ";    
      out<<p1[i].Vol<<" ";
      out<<p1[i].w1<<" "<<p1[i].w2<<" ";
      out<<p1[i].gmin<<" "<<p1[i].gmax<<"\n";
      for(int a1=0;a1<p1[i].n1.size();a1++){
	for(int a2=0;a2<p1[i].n2.size();a2++){
	  out<<"\t"<<pp.atom_name(p1[i].n1[a1])<<" "
	     <<"\t"<<pp.atom_name(p1[i].n2[a2])<<" "
	    //	     <<"\t"<<cs[p1[i].n1[a1]]<<" "<<cs[p1[i].n2[a2]]<<" "
	     <<pp.distance(p1[i].n1[a1],p1[i].n2[a2])<<" "<<p1[i].n1[a1]<<" "<<p1[i].n2[a2]<<" ";
	  //	  if(p1[i].is_12(p1[i].n1[a1],p1[i].n2[a2]))
	  //	    out<<pow(K/(P_W[i]*p1[i].W[a1][a2]*p1[i].Vol),1.0/6.0)<<" ";
	  //	  else
	  //	    out<<pow(K13/(P_W[i]*p1[i].W[a1][a2]*p1[i].Vol),1.0/6.0)<<" ";
	  //	  out<<p1[i].W[a1][a2]<<" "<<p1[i].is_12(p1[i].n1[a1],p1[i].n2[a2])<<" "
	  //	     <<p1[i].is_13(p1[i].n1[a1],p1[i].n2[a2])<<" "
	  //	     <<p1[i].is_14(p1[i].n1[a1],p1[i].n2[a2])
	  //	     <<"\n";
	  out<<p1[i].Lower[a1][a2]<<" "
	     <<p1[i].Upper[a1][a2]<<"\n";
	}
      }
      out<<"\n";
      for(int a=0;a<p1[i]._12.size();a++){
	out<<"\t bond "<<p1[i]._12[a].first<<" "<<p1[i]._12[a].second<<"\n";
      }      
      
      for(int a=0;a<p1[i]._13.size();a++){
	out<<"\t angle "<<p1[i]._13[a].first<<" "
	   <<p1[i]._13[a].second<<"\n";
      }      

      for(int a=0;a<p1[i]._14.size();a++){
	out<<"\t dihe  "<<p1[i]._14[a].first<<" "
	   <<p1[i]._14[a].second<<"\n";
      }      
      
    }
  }

  map<int, vector<pair<double,double> >  > read_topos(char * file){
    map<int, vector<pair<double,double> > > tps;
    ifstream in;
    in.open(file);
    if(!in){
      string error = string("\t>> FATAL ERROR: unable to open file: ")+file
	+string("\n");
      cout<<error;
      throw error;
    }

    istream_iterator<string> iter(in),end;
    
    while(iter!=end){
      int n = atoi(iter->c_str())-1; ++iter;
      ++iter;
      double phi  = atof(iter->c_str()); ++iter;
      double psi  = atof(iter->c_str()); ++iter;
      pair<double,double> p(phi,psi);
      tps[n].push_back(p);
    }
    return tps;
  }

  bool have_topos(int m, int M, map<int, vector<pair<double,double> > > &topos){
    for(int i=m;i<=M;i++)
      if(topos.find(i)== topos.end()) return false;
    return true;
  }


  void check_dist(BondRotation<double> & bond_rotation, Protein & p, map<int, vector<pair<double,double> > > &topos, int m, int M, int a1, int a2, int f1, int f2, double & D_m, double &D_M){

    vector<vector<int> > phi = p.phi_atoms();
    vector<vector<int> > psi = p.psi_atoms();
    
    Coor<double> c;
    c.resize(p.atom_size());
    for(int i=0;i<p.atom_size();i++){
      c.coor[Coor<double>::DIM*i] = p[i][0];
      c.coor[Coor<double>::DIM*i+1] = p[i][1];
      c.coor[Coor<double>::DIM*i+2] = p[i][2];
    }

    double d_min = 100000;
    double d_max = -1;
    for(int q=0;q<(M-m+1)*1000;q++){
      for(int i=m;i<=M;i++){
	double phi_a = topos[i][lrand48()%topos[i].size()].first*M_PI/180.0 +(drand48()-0.5)*M_PI/180.0*25.0;
	double psi_a = topos[i][lrand48()%topos[i].size()].second*M_PI/180.0+(drand48()-0.5)*M_PI/180.0*25.0;
	if(phi[i].size()<4||psi[i].size()<4) continue;
	bond_rotation.set_dihedral(c,
				   phi[i][0],
				   phi[i][1],
				   phi[i][2],
				   phi[i][3],
				   phi_a);
	
	bond_rotation.set_dihedral(c,
				   psi[i][0],
				   psi[i][1],
				   psi[i][2],
				   psi[i][3],
				   psi_a);
	
	
      }
      double d = 0;
      for(int i=0;i<3;i++){
	double dd = c.coor[Coor<double>::DIM*a1+i] -
	  c.coor[Coor<double>::DIM*a2+i];
	d +=dd*dd;
      }
      d = sqrt(d);
      if(d<d_min) d_min = d;
      if(d>d_max) d_max = d;
    }
    //    cout<<"Distance "<<d_min<<" "<<d_max<<endl;
    D_m = d_min;
    D_M = d_max;
  }
  void challange_peaks(vector<Peak2D> & peak,
		       map<int, vector<pair<double,double> > > &topos,
		       Protein &p){
    BondRotation<double> bond_rotation(p);
    int count = 0;
    for(int i=0;i<peak.size();i++){
      cerr<<">> Challenging peak "<<i+1<<" of "<<peak.size()<<endl;
      if(peak[i].n1.size()&&peak[i].n2.size()){
	for(int a=0;a<peak[i].n1.size();a++){
	  for(int b=0;b<peak[i].n2.size();b++){
	    cout<<">> Atom "<<peak[i].n1[a]<<" "
		<<p.atom_name(peak[i].n1[a])<<" "
		<<peak[i].n2[b]<<" "
		<<p.atom_name(peak[i].n2[b])<<" "
		<<p.find_fragment(peak[i].n1[a])<<" "
		<<p.find_fragment(peak[i].n2[b])<<" "
		<<"\n";
	    int f1 = p.find_fragment(peak[i].n1[a]);
	    int f2 = p.find_fragment(peak[i].n2[b]);
	    string n1 = p.atom_name(peak[i].n1[a],Protein::SHORT);
	    string n2 = p.atom_name(peak[i].n2[b],Protein::SHORT);
	    int m = min(f1,f2);
	    int M = max(f1,f2);
	    if(M-m>12){ continue;};
	    if(!have_topos(m,M,topos)){ continue;}
	    double D_Min = 0, D_Max = 0;
	    //	    if(n1!="CA"&&n1!="CB"&&n1!="C") continue;
	    //	    if(n2!="CA"&&n2!="CB"&&n2!="C") continue;
	    check_dist(bond_rotation,p,topos,m,M,peak[i].n1[a],peak[i].n2[b],f1,f2,D_Min,D_Max);
	    peak[i].Lower[a][b] = D_Min;
	    peak[i].Upper[a][b] = D_Max;
	    cout<<"Real "<<p.distance(peak[i].n1[a],peak[i].n2[b])<<" "<<peak[i].n1[a]<<" "<<peak[i].n2[b]<<" "
		<<D_Min<<" "<<D_Max<<" "<<M-m<<" "
		<<p.atom_name(peak[i].n1[a])<<" "
		<<p.atom_name(peak[i].n2[b])<<" "
		<<endl;
	  }
	}
	++count;
      }
      else { continue;}
//       if(count == 10) {
// 	cout<<"WARNING DEBUG BREAK "<<endl;
// 	return;
//       }
    }
  }
}




int main(int argc, char ** argv){

  vector<Almost::Peak2D> p1;
  vector<string>         seq;
  map<int,double>        cs;
  map<pair<int,int>,vector<int> > peak_connection;
  map<int, vector<pair<double,double> > >  topos_rest;

  vector<vector<int> > mat_12;
  vector<vector<int> > mat_13;
  vector<vector<int> > mat_14;

  p1 = Almost::read_peak(argv[1]);
  cout<<"\t>> Read "<<p1.size()<<" peaks from "<<argv[1]<<endl;

  seq = Almost::read_seq(argv[2]);
  cout<<"\t>> Read "<<seq.size()<<" amin-acids from sequence file: "<<argv[2]<<endl;


  Almost::Protein p = Almost::mkprot(seq);

  cs = Almost::read_shifts(argv[3],p);

 Almost:assign_peaks(p1,
		     p,
		     cs,
		     0.3);


  Almost::MDB mdb("all19_eef1.mdb");
  Almost::Protein pp("XRAY");
  Almost::PDB pdb(argv[4]);
  pp.build_missing(pdb[0][0],mdb,"DEFAULT","DEFAULT");




  topos_rest = Almost::read_topos(argv[5]);


  
  //    mat_12.resize(pp.atom_size());
  for(int i=0;i<pp.atom_size();i++){
    mat_12.push_back(pp.bonds(i));      
  }
  
  mat_13.resize(pp.atom_size());
  mat_14.resize(pp.atom_size());
  for(int i=0;i<mat_12.size();i++){
    int a1 = i;
    for(int j=0;j<mat_12[i].size();j++){
      int a2 = mat_12[i][j];
      for(int k=0;k<mat_12[a2].size();k++){
	int a3 = mat_12[a2][k];
	if(a1!=a3){
	  mat_13[a1].push_back(a3);
	  for(int l=0;l<mat_12[a3].size();l++){
	    int a4 = mat_12[a3][l];
	    if(a4!=a2&&a4!=a3&&a4!=a2){
	      mat_14[a1].push_back(a4);
	    }
	  }
	}
      }
    }
  }
  
  for(int i=0;i<mat_12.size(); i++){
    int a1 = i;
    for(int j=0;j<mat_12[i].size();j++){
      int a2 = mat_12[i][j];
      for(int p=0;p<p1.size();p++){
	if(p1[p].is_assignment(a1,a2,false)){
	  p1[p]._12.push_back(pair<int,int>(a1,a2));
	}
      }
    }
  }
  
  
  for(int i=0;i<mat_13.size(); i++){
    int a1 = i;
    for(int j=0;j<mat_13[i].size();j++){
      int a2 = mat_13[i][j];
      for(int p=0;p<p1.size();p++){
	if(p1[p].is_assignment(a1,a2,false)){
	  p1[p]._13.push_back(pair<int,int>(a1,a2));
	}
      }
    }
  }


  for(int i=0;i<mat_14.size(); i++){
    int a1 = i;
    for(int j=0;j<mat_14[i].size();j++){
      int a2 = mat_14[i][j];
      for(int p=0;p<p1.size();p++){
	if(p1[p].is_assignment(a1,a2,false)){
	  p1[p]._14.push_back(pair<int,int>(a1,a2));
	}
      }
    }
  }

  cout<<">> Connecting....\n";
  for(int i=0;i<p1.size();i++){
    for(int a1 = 0; a1<p1[i].n1.size();a1++){
      for(int a2 = 0; a2<p1[i].n2.size();a2++){
	pair<int,int> p(p1[i].n1[a1],p1[i].n2[a2]);
	peak_connection[p].push_back(i);
      }
    }
  }
  {
    map<pair<int,int>,vector<int> >::iterator iter = peak_connection.begin(),
      end = peak_connection.end();
    while(iter!=end){
      cout<<">> Atoms "<<iter->first.first<<" "
	  <<iter->first.second<<" "<<" ";
      for(int i=0;i<iter->second.size();i++)
	cout<<iter->second[i]<<" ";
      cout<<"\n";
      ++iter;
    }
  }


  Almost::challange_peaks(p1,topos_rest,pp); //should use LINEAR P!!!
  Almost::dump(p1,
	       pp);
}

// Local Variables:
// mode:C++
// End: