/***************************************************************************
 sandy_mod.cpp  -  description
 -------------------
 begin                : Sun Jan 25 22:22:42 2009
 copyright            : (C) 1999-2008 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 <almost.h>
#include <zeta/module.h>
#include <zeta/class.h>
#include <zeta/zextclass.h>
#include <zeta/zobj.hpp>
#include <pdb.h>
#include <mdb/mdb.h>
#include <molecules/molecules.h>
#include <moveset.h>
#include <iterator>
#include <sstream>
#include <nmr/chemicalshift.h>


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


namespace Almost {
  namespace NMR {
    typedef pair<int,int> int_pair;

    struct Peak2DFlat {
      double w1;
      double w2;
      double Vol;
      vector<int> f;
      vector<int> s;
      vector<double> W;
      vector<double> Dinv6;
      void move(){
	double ww = 0;
	for(int i=0;i<f.size();i++){
	  W[i] = max(0.05,W[i]+0.005*(drand48()-0.5));
	  ww +=W[i];
	}
	for(int i=0;i<f.size();i++)
	  W[i] /= ww;
      }

      void p_optimize(double K){
	//       for(int i=0;i<W.size();i++){
	// 	W[i] = 1.0/W.size();
	//       }

	//       Peak2DFlat tp =*this;
	//       Peak2DFlat tp_buff = tp;
	//       double E = tp.ene(K);      
	//       for(int i=0;i<10000;i++){
	// 	tp.move();
	// 	double E_b = tp.ene(K);
	// 	//	cout<<"ENE "<<E<<" "<<E_b<<endl;
	// 	if(E_b<E){
	// 	  E = E_b;
	// 	  tp_buff = tp;
	// 	  // cout<<"ENE "<<E<<" "<<E_b<<endl;
	// 	}
	// 	else {
	// 	  tp = tp_buff;
	// 	}
	
	
	//       }

	//       *this = tp_buff;

	double D = 0;
	for(int i=0;i<f.size();i++){
	  D += Dinv6[i];
	}
	for(int i=0;i<f.size();i++){
	  W[i] = Dinv6[i]/D;
	}
      }

      double ene(double K){
	double d = 0;
	for(int i=0;i<f.size();i++)
	  d += K*Dinv6[i]*W[i];

	return (1.0/Vol)*(1.0/Vol)*(Vol-d)*(Vol -d);
      }
    };

    struct Peak2D {
      double Ulm;
      double w1;
      double w2;
      double Vol;
      double Vol2;
      vector<int> n1;
      vector<double> W1;

      vector<int> n2;
      vector<double> W2;

      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;
      vector<vector<double> > D;
      vector<vector<bool> > S;
      vector<vector<double> > N;
      vector<vector<double> > Ntmp;
		
      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());
			
			
	D.clear();
	D.resize(n1.size());
	for(int i=0;i<n1.size();i++) D[i].resize(n2.size());

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

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

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

	for(int i=0;i<N.size();i++)
	  for(int j=0;j<N[i].size();j++){
	    Ntmp[i][j] =  N[i][j] = 1;
	  }
			
	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());
      }
		
		
      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 has_12(){ return _12.size()>0;}
      bool has_13(){ return _13.size()>0;}
      bool has_14(){ return _14.size()>0;}
		
      int n_12(){ return _12.size();}
      int n_13(){ return _13.size();}
      int n_14(){ return _14.size();}

      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;
      }

      int find_first(int a){
	vector<int>::iterator iter = find(n1.begin(),n1.end(),a);
	if(iter==n1.end())
	  throw "Atom not found";
	return iter-n1.begin();
      }

      int find_second(int a){
	vector<int>::iterator iter = find(n2.begin(),n2.end(),a);
	if(iter==n2.end())
	  throw "Atom not found";
	return iter-n2.begin();
      }
		
      bool is_symmetric(int a1,int a2){
	vector<int>::iterator p1 = find(n1.begin(),n1.end(),a1);
	if(p1==n1.end()) throw "NA: atom 1 not found";
	vector<int>::iterator p2 = find(n2.begin(),n2.end(),a2);
	if(p2==n2.end()) throw "NA: atom 2 not found";

	return S[p1-n1.begin()][p2-n2.begin()];
      }
		
      void init_D(){
	for(int i=0;i<D.size();i++){
	  for(int j=0;j<D[i].size();j++){
	    D[i][j] = Upper[i][j];
	    //	  if(D[i][i]>8) D[i][j] = 0;
	  }
	}
      }


      double NA(int a1, int a2){
	vector<int>::iterator p1 = find(n1.begin(),n1.end(),a1);
	if(p1==n1.end()) throw "NA: atom 1 not found";
	vector<int>::iterator p2 = find(n2.begin(),n2.end(),a2);
	if(p2==n2.end()) throw "NA: atom 2 not found";

	return N[p1-n1.begin()][p2-n2.begin()];
      }


      double upper(int a1, int a2){
	vector<int>::iterator p1 = find(n1.begin(),n1.end(),a1);
	if(p1==n1.end()) throw "NA: atom 1 not found";
	vector<int>::iterator p2 = find(n2.begin(),n2.end(),a2);
	if(p2==n2.end()) throw "NA: atom 2 not found";

	double d =  Upper[p1-n1.begin()][p2-n2.begin()];
	if(d<1) return 10000;
	return d;
      }

      double lower(int a1, int a2){
	vector<int>::iterator p1 = find(n1.begin(),n1.end(),a1);
	if(p1==n1.end()) throw "NA: atom 1 not found";
	vector<int>::iterator p2 = find(n2.begin(),n2.end(),a2);
	if(p2==n2.end()) throw "NA: atom 2 not found";

	double d =  Lower[p1-n1.begin()][p2-n2.begin()];
	if(d<0.1) return -1;
	return d;
      }


      void normalize_na(){
	for(int i=0;i<Ntmp.size();i++){
	  for(int j=0;j<Ntmp[i].size();j++){
	    if(is_12(n1[i],n2[j])||is_13(n1[i],n2[j])||is_14(n1[i],n2[j])
	       ||upper(n1[i],n2[j])<Ulm){
	      Ntmp[i][j] *= 10;
	    }

	    if(upper(n1[i],n2[j])<Ulm)
	      Ntmp[i][j] *=10;



	    Ntmp[i][j] *= exp(-0.5/25*(w1-W1[i])*(w1-W1[i])-0.5/25*(w2-W2[j])*(w2-W2[j]));

	    //  	 if(lower(n1[i],n2[j])>Ulm)
	    //  	   Ntmp[i][j] =0.0001;

	  }
	}

	double q= 0;
	for(int i=0;i<Ntmp.size();i++)
	  for(int j=0;j<Ntmp[i].size();j++)
	    q += Ntmp[i][j];

	if(q==0){ cout<<"Q is zero"<<endl; q=1;}

	for(int i=0;i<Ntmp.size();i++)
	  for(int j=0;j<Ntmp[i].size();j++){
	    Ntmp[i][j] /= q;
	    N[i][j] = Ntmp[i][j];
	    Ntmp[i][j] = 0;
	  }

      
      }
    };
	
    //Input
	
	
    double dist(const Peak2D &p1, const Peak2D & p2){
      double d1 = abs(p1.w1-p2.w1);
      double d2 = abs(p1.w2-p2.w2);
      return max(d1,d2);
    }

    double dist(double cs1, double cs2, const Peak2D & p2){
      double d1 = abs(cs1-p2.w1);
      double d2 = abs(cs2-p2.w2);
      return max(d1,d2);
    }
	
	
	
	
    class Sandy {
    public:
      vector<Peak2D> peak;
      vector<string>         seq;
      map<int,double>        cs;  
      vector<vector<int> > mat_12;
      vector<vector<int> > mat_13;
      vector<vector<int> > mat_14;  
		
      Protein p;
		
		
      map<int, vector<pair<double,double> > >  topos_rest;
		
		
		
		
      //Dist rest
      map<pair<int,int>, double>  Upper;
      map<pair<int,int>, double>  Lower;
		
    public:
      Sandy():p("BLT"){}
		
		
		
      void read_peak(const char * file, double U){
	ifstream in(file);
			
	if(!in){
	  string error = string("\t>> FATAL ERROR: unable to open file: ")+file
	    +string("\n");
	  cout<<error;
	  throw error;
	}
			
	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;
	  p.Ulm = U;
	  ss>>id;
	  ss>> p.w1;
	  ss>> p.w2;
	  ss>>id;
	  ss>> c;
	  ss>> p.Vol;
	  ss>> p.Vol2;
	  peak.push_back(p);
	}
      }
		
      void read_seq(const char * file){
	ifstream in(file);
			
	if(!in){
	  string error = string("\t>> FATAL ERROR: unable to open file: ")+file
	    +string("\n");
	  cout<<error;
	  throw error;
	}
			
			
	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";
	  }
	}
      }
		
		
      void read_shifts(const 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> lseq = p.sequence(Protein::LONG);
			
	while(in.good()){
	  char line[256];
	  in.getline(line,256);
	  if(line[0]=='#'||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 "<<lseq[nres]<<endl;
	  int natm = p.find_atom(lseq[nres]+"/"+atm);
	  {
	    //Hack for CD1 ILE...
	    if(natm==-1&&atm=="CD1"){
	      atm = "CD";
	      natm = p.find_atom(lseq[nres]+"/"+atm);
	    }
					
	  }
	  if(natm==-1){
	    cout<<"\t>> Atom not in protein\n";
	  }
				
	  if(natm!=-1){
	    cs[natm] = w;
	    cout<<"\t>> CS assignment "<<lseq[nres]+"/"+atm<<" "<<w<<"\n";
	  }
				
	}
      }
		
      void mkprotein(const char * file){
	Almost::MDB mdb("all19_eef1.mdb");
	Almost::PDB pdb(file);
	p.build_missing(pdb[0][0],mdb,"DEFAULT","DEFAULT");
			
			
	//    mat_12.resize(pp.atom_size());
	for(int i=0;i<p.atom_size();i++){
	  mat_12.push_back(p.bonds(i));
	  //Upper lower
	}
			
	mat_13.resize(p.atom_size());
	mat_14.resize(p.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++){
	  for(int j=0;j<mat_12[i].size();j++){
	    int_pair pp(i,mat_12[i][j]);
	    int_pair pinv(mat_12[i][j],i);
	    double d = p.distance(i,mat_12[i][j]);
	    Upper[pp] = d;
	    Lower[pp] = d;
	    Upper[pinv] = d;
	    Lower[pinv] = d;
	  }
	}
			
	for(int i=0;i<mat_13.size();i++){
	  for(int j=0;j<mat_13[i].size();j++){
	    int_pair pp(i,mat_13[i][j]);
	    int_pair pinv(mat_13[i][j],i);
	    double d = p.distance(i,mat_13[i][j]);
	    Upper[pp] = d;
	    Lower[pp] = d;
	    Upper[pinv] = d;
	    Lower[pinv] = d;
	  }
	}
      }
		
		
      void assign_peaks(double tol = 0.25){
			
	map<int,double>::const_iterator iter = cs.begin(), end = cs.end();
	while(iter!=end){
	  for(int i=0;i<peak.size();i++){
	    if(abs(peak[i].w1-iter->second)<tol){
	      peak[i].n1.push_back(iter->first);
	      peak[i].W1.push_back(iter->second);
	    }
	    if(abs(peak[i].w2-iter->second)<tol){
	      peak[i].n2.push_back(iter->first);
	      peak[i].W2.push_back(iter->second);
	    }

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



	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<peak.size();p++){
	      if(peak[p].is_assignment(a1,a2,false)){
		peak[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<peak.size();p++){
	      if(peak[p].is_assignment(a1,a2,false)){
		peak[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<peak.size();p++){
	      if(peak[p].is_assignment(a1,a2,false)){
		peak[p]._14.push_back(pair<int,int>(a1,a2));
	      }
	    }
	  }
	}
      
	//       {
	// 	vector<Peak2D>::iterator iter = peak.begin();
	// 	while(iter!=peak.end()){
	// 	  if(iter->n1.size()*iter->n2.size()>16)
	// 	    iter = peak.erase(iter);
	// 	  else
	// 	    ++iter;
	// 	}
	//       }
      }
		
		
      void read_topos(const char * file){
	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);
	  topos_rest[n].push_back(p);
	}
      }
		
      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, 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_rest[i][lrand48()%topos_rest[i].size()].first*M_PI/180.0 +(drand48()-0.5)*M_PI/180.0*25.0;
	    double psi_a = topos_rest[i][lrand48()%topos_rest[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(){
	BondRotation<double> bond_rotation(p);
	int count = 0;
	for(int i=0;i<peak.size();i++){
	  cout<<">> 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_rest)){ continue;}
		double D_Min = 0, D_Max = 0;
		//	    if(n1!="CA"&&n1!="CB"&&n1!="C") continue;
		//	    if(n2!="CA"&&n2!="CB"&&n2!="C") continue;
		int_pair ab(peak[i].n1[a],peak[i].n2[b]);
		int_pair ba(peak[i].n2[b],peak[i].n1[a]);
		if(Upper.find(ab)==Upper.end()||Upper[ab]<0.1){
		  check_dist(bond_rotation,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;
								
		  Upper[ab] = D_Max;
		  Upper[ba] = D_Max;
								
		  Lower[ab] = D_Min;
		  Lower[ba] = D_Min;
		}
		else { 
		  cout<<"DEBUG dist in db"<<endl;
		  peak[i].Lower[a][b] = Lower[ab];
		  peak[i].Upper[a][b] = Upper[ab];
		  cout<<"DEBUG dist in db " << peak[i].Upper[a][b]<<" "<<Upper[ab]<<"\n";;
		}
	      }
	    }
	    ++count;
	  }
	  else { continue;}
	  //       if(count == 10) {
	  // 	cout<<"WARNING DEBUG BREAK "<<endl;
	  // 	return;
	  //       }
	}
      }
		
      void normalize_na(){
	for(int i=0;i<peak.size();i++)
	  peak[i].normalize_na();
      }

      void dump(const char * dump_peaks){
	ofstream out;
	out.open(dump_peaks);
	for(int i=0;i<peak.size();i++){
	  // 	if(peak[i].has_12()||
	  // 	   peak[i].has_13()||
	  // 	   peak[i].has_14()) continue;
	  //    if(peak[i].n1.size()>2||peak[i].n2.size()>2) continue;
	  //    if(peak[i].gmax>10||peak[i].gmin<2) continue;
	  out<<">> Peak "<<i+1<<" "//<<P_W[i]<<" "
	     <<peak[i].n1.size()<<" ";
	  out<<peak[i].n2.size()<<" ";
	  out<< "{ ";
	  for(int a=0;a<peak[i].n1.size();a++) out<<peak[i].n1[a]<<" ";
	  out<< "} ";    
	  out<< "{ ";
	  for(int a=0;a<peak[i].n2.size();a++) out<<peak[i].n2[a]<<" ";
	  out<< "} ";    
	  out<<peak[i].Vol<<" ";
	  out<<peak[i].w1<<" "<<peak[i].w2<<" ";
	  out<<peak[i].gmin<<" "<<peak[i].gmax<<"\n";
	  for(int a1=0;a1<peak[i].n1.size();a1++){
	    for(int a2=0;a2<peak[i].n2.size();a2++){
	      out<<"\t"<<p.atom_name(peak[i].n1[a1])<<" "
		 <<"\t"<<p.atom_name(peak[i].n2[a2])<<" "
		//	     <<"\t"<<cs[peak[i].n1[a1]]<<" "<<cs[peak[i].n2[a2]]<<" "
		 <<p.distance(peak[i].n1[a1],peak[i].n2[a2])<<" "<<peak[i].n1[a1]<<" "<<peak[i].n2[a2]<<" ";
	      //	  if(peak[i].is_12(peak[i].n1[a1],peak[i].n2[a2]))
	      //	    out<<pow(K/(P_W[i]*peak[i].W[a1][a2]*peak[i].Vol),1.0/6.0)<<" ";
	      //	  else
	      //	    out<<pow(K13/(P_W[i]*peak[i].W[a1][a2]*peak[i].Vol),1.0/6.0)<<" ";
	      //	  out<<peak[i].W[a1][a2]<<" "<<peak[i].is_12(peak[i].n1[a1],peak[i].n2[a2])<<" "
	      //	     <<peak[i].is_13(peak[i].n1[a1],peak[i].n2[a2])<<" "
	      //	     <<peak[i].is_14(peak[i].n1[a1],peak[i].n2[a2])
	      //	     <<"\n";
	      out<<peak[i].Lower[a1][a2]<<" "
		 <<peak[i].Upper[a1][a2]<<" "
		 <<peak[i].D[a1][a2]<<" " <<peak[i].W[a1][a2]<<" "<<peak[i].Ntmp[a1][a2]
		 <<"\n";
	    }
	  }
	  out<<"\n";
	  for(int a=0;a<peak[i]._12.size();a++){
	    out<<"\t bond "<<peak[i]._12[a].first<<" "<<peak[i]._12[a].second<<"\n";
	  }      
				
	  for(int a=0;a<peak[i]._13.size();a++){
	    out<<"\t angle "<<peak[i]._13[a].first<<" "
	       <<peak[i]._13[a].second<<"\n";
	  }      
				
	  for(int a=0;a<peak[i]._14.size();a++){
	    out<<"\t dihe  "<<peak[i]._14[a].first<<" "
	       <<peak[i]._14[a].second<<"\n";
	  }      
				
	}
			
	//Dump dist
	{
	  map<pair<int,int>, double>::iterator  iter = Upper.begin(),end = Upper.end();
	  ofstream out;
	  out.open("upper.dist");
	  while(iter!=end){
	    if(iter->second!=0)
	      out<<iter->first.first<<" "<<iter->first.second<<" "<<iter->second<<"\n";
	    ++iter;
	  }
	}
			
	{
	  map<pair<int,int>, double>::iterator  iter = Lower.begin(),end = Lower.end();
	  ofstream out;
	  out.open("lower.dist");
	  while(iter!=end){
	    if(iter->second!=0)
	      out<<iter->first.first<<" "<<iter->first.second<<" "<<iter->second<<"\n";
	    ++iter;
	  }
	}
			
      }
		
		
		
      void dump2(vector<double> & P_W, double K){
	ofstream out;
	out.open("dump.peaks");
	for(int i=0;i<peak.size();i++){
	  //    if(peak[i].n1.size()>2||peak[i].n2.size()>2) continue;
	  //    if(peak[i].gmax>10||peak[i].gmin<2) continue;
	  out<<">> Peak "<<i+1<<" "//<<P_W[i]<<" "
	     <<peak[i].n1.size()<<" ";
	  out<<peak[i].n2.size()<<" ";
	  out<< "{ ";
	  for(int a=0;a<peak[i].n1.size();a++) out<<peak[i].n1[a]<<" ";
	  out<< "} ";    
	  out<< "{ ";
	  for(int a=0;a<peak[i].n2.size();a++) out<<peak[i].n2[a]<<" ";
	  out<< "} ";    
	  out<<peak[i].Vol<<" ";
	  out<<peak[i].w1<<" "<<peak[i].w2<<" ";
	  out<<peak[i].gmin<<" "<<peak[i].gmax<<"\n";
	  double Ip= 0;
	  for(int a1=0;a1<peak[i].n1.size();a1++){
	    for(int a2=0;a2<peak[i].n2.size();a2++){
	      out<<"\t"<<p.atom_name(peak[i].n1[a1])<<" "
		 <<"\t"<<p.atom_name(peak[i].n2[a2])<<" "
		//	     <<"\t"<<cs[peak[i].n1[a1]]<<" "<<cs[peak[i].n2[a2]]<<" "
		 <<p.distance(peak[i].n1[a1],peak[i].n2[a2])<<" "<<peak[i].n1[a1]<<" "<<peak[i].n2[a2]<<" ";
	      //	  if(peak[i].is_12(peak[i].n1[a1],peak[i].n2[a2]))
	      //	    out<<pow(K/(P_W[i]*peak[i].W[a1][a2]*peak[i].Vol),1.0/6.0)<<" ";
	      //	  else
	      //	    out<<pow(K13/(P_W[i]*peak[i].W[a1][a2]*peak[i].Vol),1.0/6.0)<<" ";
	      //	  out<<peak[i].W[a1][a2]<<" "<<peak[i].is_12(peak[i].n1[a1],peak[i].n2[a2])<<" "
	      //	     <<peak[i].is_13(peak[i].n1[a1],peak[i].n2[a2])<<" "
	      //	     <<peak[i].is_14(peak[i].n1[a1],peak[i].n2[a2])
	      //	     <<"\n";
	      out<<peak[i].Lower[a1][a2]<<" "
		 <<peak[i].Upper[a1][a2]<<" "
		 <<peak[i].D[a1][a2]<<" "
		 <<"\n";
	      if(peak[i].D[a1][a2]>1)
		Ip += pow(peak[i].D[a1][a2],-6.0)*K;
						
	    }
	  }
	  out<<"\t "<<peak[i].Vol*P_W[i]<<" "<<Ip<<" "<<P_W[i]<<"\n";
	  out<<"\n";
	}
      }
		
      void read_dist(){
	{
	  ifstream in;
	  in.open("upper.dist");
	  istream_iterator<string> iter(in),end;
				
	  while(iter != end){
	    int_pair p;
	    double d;
					
	    p.first = atoi(iter->c_str()); ++iter; 
	    p.second = atoi(iter->c_str()); ++iter; 
	    d = atof(iter->c_str()); ++iter;
	    int_pair q; q.first = p.second; q.second = p.first;
	    // 	  if(p.first==279&&p.second==255){
	    // 	    cout<<"DEDDD "<<d<<endl;
	    // 	    exit(-1);
	    // 	  }
	    if(d>0){
	      Upper[p] = d;
	      Upper[q] = d;
	    }
	  }
	}
			
			
	{
	  ifstream in;
	  in.open("lower.dist");
	  istream_iterator<string> iter(in),end;
				
	  while(iter != end){
	    int_pair p;
	    double d;
					
	    p.first = atoi(iter->c_str()); ++iter; 
	    p.second = atoi(iter->c_str()); ++iter; 
	    d = atof(iter->c_str()); ++iter;
	    int_pair q; q.first = p.second; q.second = p.first;
	    if(d>0){
	      Lower[p] = d;
	      Lower[q] = d;
	    }
	  }
	}
			
			
	//Ass dist to peaks;
	for(int i=0;i<peak.size();i++){
	  for(int a = 0; a<peak[i].n1.size();a++){
	    for(int b=0;b<peak[i].n2.size();b++){
	      peak[i].Upper[a][b] = Upper[int_pair(peak[i].n1[a],
						   peak[i].n2[b])];
						
	      peak[i].Lower[a][b] = Lower[int_pair(peak[i].n1[a],
						   peak[i].n2[b])];
	    }
	  }
	}
			
      }
				
      void dump_symmetric(){

	map<int_pair,set<int> > ab;

	for(int i=0;i<peak.size();i++){
	  for(int a=0;a<peak[i].n1.size();a++){
	    for(int b=0;b<peak[i].n2.size();b++){
	      int_pair p(peak[i].n1[a],peak[i].n2[b]);
	      ab[p].insert(i);
	    }
	  }
	}
	for(int i=0;i<peak.size();i++){
	  for(int a=0;a<peak[i].n1.size();a++){
	    for(int b=0;b<peak[i].n2.size();b++){
	      int_pair pp(peak[i].n1[a],peak[i].n2[b]);
	      int_pair qp(peak[i].n2[b],peak[i].n1[a]);
	      //	      cout<<">> Peak "<<i+1<<"\n";
	      if(ab.find(pp)!=ab.end()&&ab.find(qp)!=ab.end()){
		cout<<">>Peak "<<i+1<<" Symmetric "<<p.atom_name(pp.first)<<" "<<p.atom_name(pp.second)<<" "
		    <<p.distance(qp.first,qp.second)<<" "
		    <<Lower[pp]<<" "<<Upper[pp]
		    <<endl;
		peak[i].S[a][b]=1;
	      }
	    }
	  }
	}
      }

      void na_bond(){
      
	map<int_pair,set<int> > ab;

	for(int i=0;i<peak.size();i++){
	  for(int a=0;a<peak[i].n1.size();a++){
	    for(int b=0;b<peak[i].n2.size();b++){
	      int_pair p(peak[i].n1[a],peak[i].n2[b]);
	      ab[p].insert(i);
	    }
	  }
	}

	for(int i=0;i<peak.size();i++){
	  for(int a=0;a<peak[i].n1.size();a++){
	    for(int b=0;b<peak[i].n2.size();b++){
	      //	    if(peak[i].is_12(peak[i].n1[a],peak[i].n2[b])) continue;
	      //	    if(peak[i].is_13(peak[i].n1[a],peak[i].n2[b])) continue;
	    
	      int a1 = peak[i].n1[a];
	      for(int q=0;q<mat_12[a1].size();q++){
		if(peak[i].n2[b]==mat_12[a1][q]) continue;
		int_pair pp(peak[i].n2[b],mat_12[a1][q]);
		int_pair qp(mat_12[a1][q],peak[i].n2[b]);
		if(ab.find(pp)!=ab.end()||ab.find(qp)!=ab.end()){
		  // 		cout<<">> Triangle "<<p.atom_name(peak[i].n1[a])<<" "
		  // 		    <<p.atom_name(peak[i].n2[b])<<" "<<p.atom_name(mat_12[a1][q])<<" "
		  // 		    <<p.distance(peak[i].n1[a],peak[i].n2[b])<<" "
		  // 		    <<p.distance(peak[i].n2[b],mat_12[a1][q])<<" "
		  // 		    <<is_symmentric(peak[i].n1[a],peak[i].n2[a],ab)<<" "
		  // 		    <<is_symmentric(peak[i].n2[b],mat_12[a1][q],ab)<<" "
		  // 		    <<"\n";

		  //Contributions
		  {
		  
		    //N of bond a-b
		    int b1 = mat_12[a1][q];
		    double Vab = 0;
		    int_pair bp(a1,b1);
		    int_pair bb(b1,a1);
		    if(ab.find(bp)!=ab.end()){
		      set<int>::iterator iter = ab[bp].begin(),end = ab[bp].end();
		      Vab += peak[*iter].NA(bp.first,bp.second);
		      ++iter;
		    }

		    if(ab.find(bb)!=ab.end()){
		      set<int>::iterator iter = ab[bb].begin(),end = ab[bb].end();
		      Vab += peak[*iter].NA(bb.first,bb.second);
		      ++iter;
		    }

		    Vab = max(1.0,Vab);

		    if(ab.find(pp)!=ab.end()){

		      set<int>::iterator iter = ab[pp].begin(),end = ab[pp].end();
		      while(iter!=end){
			peak[i].Ntmp[a][b] += (Vab)*(peak[*iter].NA(peak[i].n2[b],mat_12[a1][q]));
			++iter;
		      }
		    }
		    if(ab.find(qp)!=ab.end()){
		      set<int>::iterator iter = ab[qp].begin(),end = ab[qp].end();
		      while(iter!=end){
			peak[i].Ntmp[a][b] += (Vab)*(peak[*iter].NA(mat_12[a1][q],peak[i].n2[b]));
			++iter;
		      }
		    }
		  }

		}
	      }

	      int a2 = peak[i].n2[b];
	      for(int q=0;q<mat_12[a2].size();q++){
		if(peak[i].n1[a]==mat_12[a2][q]) continue;
		int_pair pp(peak[i].n1[a],mat_12[a2][q]);
		int_pair qp(mat_12[a2][q],peak[i].n1[a]);
		if(ab.find(pp)!=ab.end()||ab.find(qp)!=ab.end()){
		  // 		cout<<">> Triangle "<<p.atom_name(peak[i].n2[b])<<" "
		  // 		    <<p.atom_name(peak[i].n1[a])<<" "<<p.atom_name(mat_12[a2][q])<<" "
		  // 		    <<p.distance(peak[i].n1[a],peak[i].n2[b])<<" "
		  // 		    <<p.distance(peak[i].n1[a],mat_12[a2][q])<<" "
		  // 		    <<is_symmentric(peak[i].n2[b],peak[i].n1[a],ab)<<" "
		  // 		    <<is_symmentric(peak[i].n1[a],mat_12[a2][q],ab)<<" "
		  // 		    <<"\n";

		  //Contributions
		  {

		  
		    //N of bond a-b
		    int b1 = mat_12[a2][q];
		    double Vab = 0;
		    int_pair bp(a1,b1);
		    int_pair bb(b1,a1);
		    if(ab.find(bp)!=ab.end()){
		      set<int>::iterator iter = ab[bp].begin(),end = ab[bp].end();
		      Vab += peak[*iter].NA(bp.first,bp.second);
		      ++iter;
		    }

		    if(ab.find(bb)!=ab.end()){
		      set<int>::iterator iter = ab[bb].begin(),end = ab[bb].end();
		      Vab += peak[*iter].NA(bb.first,bb.second);
		      ++iter;
		    }

		    Vab = max(1.0,Vab);
		  
		    if(ab.find(pp)!=ab.end()){
		      set<int>::iterator iter = ab[pp].begin(),end = ab[pp].end();
		      while(iter!=end){
			peak[i].Ntmp[a][b] += (Vab)*(peak[*iter].NA(pp.first,pp.second));
			++iter;
		      }
		    }
		    if(ab.find(qp)!=ab.end()){
		      set<int>::iterator iter = ab[qp].begin(),end = ab[qp].end();
		      while(iter!=end){
			peak[i].Ntmp[a][b] += (Vab)*(peak[*iter].NA(qp.first,qp.second));
			++iter;
		      }
		    }
		  }

		}
	      }	    

	    }
	  }
	}
      
      }


      void na_angle(){
      
	map<int_pair,set<int> > ab;

	for(int i=0;i<peak.size();i++){
	  for(int a=0;a<peak[i].n1.size();a++){
	    for(int b=0;b<peak[i].n2.size();b++){
	      int_pair p(peak[i].n1[a],peak[i].n2[b]);
	      ab[p].insert(i);
	    }
	  }
	}

	for(int i=0;i<peak.size();i++){
	  for(int a=0;a<peak[i].n1.size();a++){
	    for(int b=0;b<peak[i].n2.size();b++){
	      //	    if(peak[i].is_12(peak[i].n1[a],peak[i].n2[b])) continue;
	      //	    if(peak[i].is_13(peak[i].n1[a],peak[i].n2[b])) continue;
	    
	      int a1 = peak[i].n1[a];
	      for(int q=0;q<mat_13[a1].size();q++){
		int_pair pp(peak[i].n2[b],mat_13[a1][q]);
		int_pair qp(mat_13[a1][q],peak[i].n2[b]);
		if(ab.find(pp)!=ab.end()||ab.find(qp)!=ab.end()){
		  // 		cout<<">> TriAngle "<<p.atom_name(peak[i].n1[a])<<" "
		  // 		    <<p.atom_name(peak[i].n2[b])<<" "<<p.atom_name(mat_13[a1][q])<<" "
		  // 		    <<p.distance(peak[i].n1[a],peak[i].n2[b])<<" "
		  // 		    <<p.distance(peak[i].n2[b],mat_13[a1][q])<<" "
		  // 		    <<is_symmentric(peak[i].n1[a],peak[i].n2[a],ab)<<" "
		  // 		    <<is_symmentric(peak[i].n2[b],mat_13[a1][q],ab)<<" "
		  // 		    <<"\n";

		  //Contributions
		  {
		  
		    //N of bond a-b
		    int b1 = mat_13[a1][q];
		    double Vab = 0;
		    int_pair bp(a1,b1);
		    int_pair bb(b1,a1);
		    if(ab.find(bp)!=ab.end()){
		      set<int>::iterator iter = ab[bp].begin(),end = ab[bp].end();
		      Vab += peak[*iter].NA(bp.first,bp.second);
		      ++iter;
		    }

		    if(ab.find(bb)!=ab.end()){
		      set<int>::iterator iter = ab[bb].begin(),end = ab[bb].end();
		      Vab += peak[*iter].NA(bb.first,bb.second);
		      ++iter;
		    }

		    Vab = max(1.0,Vab);
		  
		    if(ab.find(pp)!=ab.end()){
		      set<int>::iterator iter = ab[pp].begin(),end = ab[pp].end();
		      while(iter!=end){
			peak[i].Ntmp[a][b] += (Vab)*(peak[*iter].NA(pp.first,pp.second));
			++iter;
		      }
		    }
		    if(ab.find(qp)!=ab.end()){
		      set<int>::iterator iter = ab[qp].begin(),end = ab[qp].end();
		      while(iter!=end){
			peak[i].Ntmp[a][b] += (Vab)*(peak[*iter].NA(qp.first,qp.second));
			++iter;
		      }
		    }
		  }

		}
	      }

	      int a2 = peak[i].n2[b];
	      for(int q=0;q<mat_13[a2].size();q++){
		int_pair pp(peak[i].n1[a],mat_13[a2][q]);
		int_pair qp(mat_13[a2][q],peak[i].n1[a]);
		if(ab.find(pp)!=ab.end()||ab.find(qp)!=ab.end()){
		  // 		cout<<">> TriAngle "<<p.atom_name(peak[i].n2[b])<<" "
		  // 		    <<p.atom_name(peak[i].n1[a])<<" "<<p.atom_name(mat_13[a2][q])<<" "
		  // 		    <<p.distance(peak[i].n1[a],peak[i].n2[b])<<" "
		  // 		    <<p.distance(peak[i].n1[a],mat_13[a2][q])<<" "
		  // 		    <<is_symmentric(peak[i].n2[b],peak[i].n1[a],ab)<<" "
		  // 		    <<is_symmentric(peak[i].n1[a],mat_13[a2][q],ab)<<" "
		  // 		    <<"\n";

		  //Contributions
		  {

		  
		    //N of bond a-b
		    int b1 = mat_13[a2][q];
		    double Vab = 0;
		    int_pair bp(a1,b1);
		    int_pair bb(b1,a1);
		    if(ab.find(bp)!=ab.end()){
		      set<int>::iterator iter = ab[bp].begin(),end = ab[bp].end();
		      Vab += peak[*iter].NA(bp.first,bp.second);
		      ++iter;
		    }

		    if(ab.find(bb)!=ab.end()){
		      set<int>::iterator iter = ab[bb].begin(),end = ab[bb].end();
		      Vab += peak[*iter].NA(bb.first,bb.second);
		      ++iter;
		    }

		    Vab = max(1.0,Vab);
		  
		    if(ab.find(pp)!=ab.end()){
		      set<int>::iterator iter = ab[pp].begin(),end = ab[pp].end();
		      while(iter!=end){
			peak[i].Ntmp[a][b] += (Vab)*(peak[*iter].NA(pp.first,pp.second));
			++iter;
		      }
		    }
		    if(ab.find(qp)!=ab.end()){
		      set<int>::iterator iter = ab[qp].begin(),end = ab[qp].end();
		      while(iter!=end){
			peak[i].Ntmp[a][b] += (Vab)*(peak[*iter].NA(qp.first,qp.second));
			++iter;
		      }
		    }
		  }

		}
	      }	    

	    }
	  }
	}
      
      }


      void na_dihe(){
      
	map<int_pair,set<int> > ab;

	for(int i=0;i<peak.size();i++){
	  for(int a=0;a<peak[i].n1.size();a++){
	    for(int b=0;b<peak[i].n2.size();b++){
	      int_pair p(peak[i].n1[a],peak[i].n2[b]);
	      ab[p].insert(i);
	    }
	  }
	}

	for(int i=0;i<peak.size();i++){
	  for(int a=0;a<peak[i].n1.size();a++){
	    for(int b=0;b<peak[i].n2.size();b++){
	      //	    if(peak[i].is_12(peak[i].n1[a],peak[i].n2[b])) continue;
	      //	    if(peak[i].is_13(peak[i].n1[a],peak[i].n2[b])) continue;
	    
	      int a1 = peak[i].n1[a];
	      for(int q=0;q<mat_14[a1].size();q++){
		int_pair pp(peak[i].n2[b],mat_14[a1][q]);
		int_pair qp(mat_14[a1][q],peak[i].n2[b]);
		if(ab.find(pp)!=ab.end()||ab.find(qp)!=ab.end()){
		  // 		cout<<">> TriAngle "<<p.atom_name(peak[i].n1[a])<<" "
		  // 		    <<p.atom_name(peak[i].n2[b])<<" "<<p.atom_name(mat_14[a1][q])<<" "
		  // 		    <<p.distance(peak[i].n1[a],peak[i].n2[b])<<" "
		  // 		    <<p.distance(peak[i].n2[b],mat_14[a1][q])<<" "
		  // 		    <<is_symmentric(peak[i].n1[a],peak[i].n2[a],ab)<<" "
		  // 		    <<is_symmentric(peak[i].n2[b],mat_14[a1][q],ab)<<" "
		  // 		    <<"\n";

		  //Contributions
		  {

		  
		    //N of bond a-b
		    int b1 = mat_14[a1][q];
		    double Vab = 0;
		    int_pair bp(a1,b1);
		    int_pair bb(b1,a1);
		    if(ab.find(bp)!=ab.end()){
		      set<int>::iterator iter = ab[bp].begin(),end = ab[bp].end();
		      Vab += peak[*iter].NA(bp.first,bp.second);
		      ++iter;
		    }

		    if(ab.find(bb)!=ab.end()){
		      set<int>::iterator iter = ab[bb].begin(),end = ab[bb].end();
		      Vab += peak[*iter].NA(bb.first,bb.second);
		      ++iter;
		    }

		    Vab = max(1.0,Vab);
		  
		    if(ab.find(pp)!=ab.end()){
		      set<int>::iterator iter = ab[pp].begin(),end = ab[pp].end();
		      while(iter!=end){
			peak[i].Ntmp[a][b] += (Vab)*(peak[*iter].NA(pp.first,pp.second));
			++iter;
		      }
		    }
		    if(ab.find(qp)!=ab.end()){
		      set<int>::iterator iter = ab[qp].begin(),end = ab[qp].end();
		      while(iter!=end){
			peak[i].Ntmp[a][b] += (Vab)*(peak[*iter].NA(qp.first,qp.second));
			++iter;
		      }
		    }
		  }

		}
	      }

	      int a2 = peak[i].n2[b];
	      for(int q=0;q<mat_14[a2].size();q++){
		int_pair pp(peak[i].n1[a],mat_14[a2][q]);
		int_pair qp(mat_14[a2][q],peak[i].n1[a]);
		if(ab.find(pp)!=ab.end()||ab.find(qp)!=ab.end()){
		  // 		cout<<">> TriAngle "<<p.atom_name(peak[i].n2[b])<<" "
		  // 		    <<p.atom_name(peak[i].n1[a])<<" "<<p.atom_name(mat_14[a2][q])<<" "
		  // 		    <<p.distance(peak[i].n1[a],peak[i].n2[b])<<" "
		  // 		    <<p.distance(peak[i].n1[a],mat_14[a2][q])<<" "
		  // 		    <<is_symmentric(peak[i].n2[b],peak[i].n1[a],ab)<<" "
		  // 		    <<is_symmentric(peak[i].n1[a],mat_14[a2][q],ab)<<" "
		  // 		    <<"\n";

		  //Contributions
		  {
		  
		    //N of bond a-b
		    int b1 = mat_14[a2][q];
		    double Vab = 0;
		    int_pair bp(a1,b1);
		    int_pair bb(b1,a1);
		    if(ab.find(bp)!=ab.end()){
		      set<int>::iterator iter = ab[bp].begin(),end = ab[bp].end();
		      Vab += peak[*iter].NA(bp.first,bp.second);
		      ++iter;
		    }

		    if(ab.find(bb)!=ab.end()){
		      set<int>::iterator iter = ab[bb].begin(),end = ab[bb].end();
		      Vab += peak[*iter].NA(bb.first,bb.second);
		      ++iter;
		    }

		    Vab = max(1.0,Vab);
		  
		    if(ab.find(pp)!=ab.end()){
		      set<int>::iterator iter = ab[pp].begin(),end = ab[pp].end();
		      while(iter!=end){
			peak[i].Ntmp[a][b] += (Vab)*(peak[*iter].NA(pp.first,pp.second));
			++iter;
		      }
		    }
		    if(ab.find(qp)!=ab.end()){
		      set<int>::iterator iter = ab[qp].begin(),end = ab[qp].end();
		      while(iter!=end){
			peak[i].Ntmp[a][b] += (Vab)*(peak[*iter].NA(qp.first,qp.second));
			++iter;
		      }
		    }
		  }

		}
	      }	    

	    }
	  }
	}
      
      }

      double Vval(int a1, int a2, map<int_pair,set<int> > &ab){
	double Res = 0;
	{
	  int_pair pp(a1,a2);
	  set<int>::iterator iter = ab[pp].begin(), end = ab[pp].end();
	  while(iter!=end){
	    Res += peak[*iter].NA(pp.first,pp.second);
	    ++iter;
	  }
	}

	{
	  int_pair pp(a2,a1);
	  set<int>::iterator iter = ab[pp].begin(), end = ab[pp].end();
	  while(iter!=end){
	    Res += peak[*iter].NA(pp.first,pp.second);
	    ++iter;
	  }
	}
      

	if(is_12(a1,a2)||is_13(a1,a2)||is_14(a1,a2))
	  Res = max(1.0,Res);
	return Res;
      }


      vector<int> g_atoms(int a1, int a2, int fa1, int fa2,  map<int_pair,set<int> > &ab){
	vector<int> g;

	for(int aa = max(0,fa1-1);aa<min(fa1+2,p.fragment_size());aa++){
	  vector<int> atm = p.fragment_atoms(aa);
	  for(int q=0;q<atm.size();q++){
	    if(p[atm[q]].name()[0]!='C') continue;
	    if(atm[q]==a1||atm[q]==a2) continue;
	    //is connected to a
	    int ga = atm[q];
	    bool a1_conn = false;
	    bool a2_conn = false;
	    if(is_12(a1,ga)||is_13(a1,ga)||is_14(a1,ga)||
	       ab.find(int_pair(a1,ga))!=ab.end()||
	       ab.find(int_pair(ga,a1))!=ab.end()){
	  
	      a1_conn = true;
	    }
	    if(!a1_conn) continue;

	    if(is_12(a2,ga)||is_13(a2,ga)||is_14(a2,ga)||
	       ab.find(int_pair(a2,ga))!=ab.end()||
	       ab.find(int_pair(ga,a2))!=ab.end()){
	  
	      a2_conn = true;
	    }
	    if(!a2_conn) continue;	  

	    //is connecting
	    g.push_back(ga);
	  }
	}

	for(int aa = max(0,fa2-1);aa<min(fa2+2,p.fragment_size());aa++){
	  vector<int> atm = p.fragment_atoms(aa);
	  for(int q=0;q<atm.size();q++){
	    if(p[atm[q]].name()[0]!='C') continue;
	    if(atm[q]==a1||atm[q]==a2) continue;
	    //is connected to a
	    int ga = atm[q];
	    bool a1_conn = false;
	    bool a2_conn = false;
	    if(is_12(a1,ga)||is_13(a1,ga)||is_14(a1,ga)||
	       ab.find(int_pair(a1,ga))!=ab.end()||
	       ab.find(int_pair(ga,a1))!=ab.end()){
	  
	      a1_conn = true;
	    }
	    if(!a1_conn) continue;

	    if(is_12(a2,ga)||is_13(a2,ga)||is_14(a2,ga)||
	       ab.find(int_pair(a2,ga))!=ab.end()||
	       ab.find(int_pair(ga,a2))!=ab.end()){
	  
	      a2_conn = true;
	    }
	    if(!a2_conn) continue;	  

	    //is connecting
	    g.push_back(ga);
	  }
	}
	return g;

      }

      bool is_12(int a1, int a2){
	return 
	  (find(mat_12[a1].begin(),mat_12[a1].end(),a2)!= mat_12[a1].end())||
	  (find(mat_12[a2].begin(),mat_12[a2].end(),a1)!= mat_12[a2].end());
      }

      bool is_13(int a1, int a2){
	return 
	  (find(mat_13[a1].begin(),mat_13[a1].end(),a2)!= mat_13[a1].end())||
	  (find(mat_13[a2].begin(),mat_13[a2].end(),a1)!= mat_13[a2].end());
      }

      bool is_14(int a1, int a2){
	return 
	  (find(mat_14[a1].begin(),mat_14[a1].end(),a2)!= mat_14[a1].end())||
	  (find(mat_14[a2].begin(),mat_14[a2].end(),a1)!= mat_14[a2].end());
      }



      void na(){

	map<int_pair,set<int> > ab;

	for(int i=0;i<peak.size();i++){
	  for(int a=0;a<peak[i].n1.size();a++){
	    for(int b=0;b<peak[i].n2.size();b++){
	      int_pair p(peak[i].n1[a],peak[i].n2[b]);
	      ab[p].insert(i);
	    }
	  }
	}

	//Compute NA for pair a1 a2
	map<int_pair,set<int> >::iterator iter = ab.begin(),
	  end = ab.end();

	while(iter!=end){
	  //	cout<<"Doing...."<<endl;
	  int a  = iter->first.first;
	  int fa = p.find_fragment(a);
	  int b = iter->first.second;
	  int fb = p.find_fragment(b);
	  //search g connected to a and b
	  //g has to be in [fa-1, fa +1] or [fb-1,fb+1]
	  vector<int> g  = g_atoms(a,b,fa,fb,ab);
	  vector<int> g2 = g_atoms(b,a,fb,fa,ab);
	  //	for(int gg=0;gg<g2.size();gg++) g.push_back(g2[gg]);
	  double Nab = 0;
	  for(int q=0;q<g.size();q++){
	    double Vag = Vval(a,g[q],ab);
	    double Vgb = Vval(g[q],b,ab);
	    if(Vag<0||Vgb<0) throw "Vval error";
	    Nab += Vag*Vgb; 
	    //	    sqrt(Vag)*sqrt(Vgb);
	  }
	  //	cout<<"Nab "<<Nab<<" "<<p.atom_name(a)<<" "<<p.atom_name(b)<<" "<<g.size()<<endl;
	  //assign
	  {
	    set<int>::iterator p_iter = iter->second.begin(),
	      p_end = iter->second.end();
	    while(p_iter!=p_end){
	      int p_a = peak[*p_iter].find_first(a);
	      int p_b = peak[*p_iter].find_second(b);
	      peak[*p_iter].Ntmp[p_a][p_b] = Nab;
	      ++p_iter;
	    }
	  }
	  ++iter; // done
	}
      }


      void na_residue(){
	ofstream out;
	out.open("na.res");
	map<int_pair,set<int> > ab;
      
	for(int i=0;i<peak.size();i++){
	  for(int a=0;a<peak[i].n1.size();a++){
	    for(int b=0;b<peak[i].n2.size();b++){
	      int_pair p(peak[i].n1[a],peak[i].n2[b]);
	      ab[p].insert(i);
	    }
	  }
	}
      
	vector<int> CA;
	for(int aa = 0; aa<p.atom_size();aa++)
	  if(p[aa].name()=="CA") CA.push_back(aa);

	for(int aa = 0; aa<p.fragment_size();aa++){
	  vector<int> atm1 = p.fragment_atoms(aa);
	  for(int bb = 0; bb<p.fragment_size();bb++){
	    if(abs(aa-bb)<5) continue;
	    double NA = 0;
	    int counter = 0;
	    vector<int> atm2 = p.fragment_atoms(bb);
	    for(int a1 = 0; a1< atm1.size();a1++){
	      if(p[atm1[a1]].name()[0]!='C') continue;
	      for(int a2 = 0;a2<atm2.size();a2++){
		if(p[atm2[a2]].name()[0]!='C') continue;
		int_pair pp(atm1[a1],atm2[a2]);
		set<int>::iterator p_iter = ab[pp].begin(),
		  p_end = ab[pp].end();
		while(p_iter!=p_end){
		  int p_a = peak[*p_iter].find_first(atm1[a1]);
		  int p_b = peak[*p_iter].find_second(atm2[a2]);
		  NA +=  peak[*p_iter].Ntmp[p_a][p_b];
		  ++p_iter;
		  ++counter;
		}
		//	      D += p.distance(atm1[a1],atm2[a2]);
	      }
	    }
	    out<<p.fragment_name(aa)<<" "
	       <<p.fragment_name(bb)<<" "
	       <<NA<<" "<<NA/counter<<" "
	       <<p.distance(CA[aa],CA[bb])<<endl;
	  }		
	}

      }


      bool is_symmentric(int a1,int a2,map<int_pair,set<int> > &ab){
	int_pair pp(a1,a2);
	int_pair qp(a2,a1);
	if(ab.find(pp)!=ab.end()&&ab.find(qp)!=ab.end()) return true;
	return false;
      }

      double bond_tol(){
	double T = 0;
	for(int i = 0; i< mat_12.size(); i++){
	  int a1 = i;
	  if(cs.find(a1) == cs.end()) continue;
	  for(int j=0;j<mat_12[i].size();j++){
	    int a2 = mat_12[i][j];
	    if(a2==a1) continue;
	    if(cs.find(a2) == cs.end()) continue;

	    double D = 100000;
	    //Doit 
	    for(int p=0; p<peak.size();p++){
	      double d = dist(cs[a1],cs[a2],peak[p]);
	      d = min(d,dist(cs[a1],cs[a2],peak[p]));
	      if(d<D) D =d;
	    }

	    if(D>T&&D<0.5) T=D;
	    cout<<"Tol Pair "<<a1<<" "<<a2<<" "<<D<<" "<<T<<endl;
	  }
	}
	return T;
      }

      double angle_tol(){
	double T = 0;
	for(int i = 0; i< mat_13.size(); i++){
	  int a1 = i;
	  if(cs.find(a1) == cs.end()) continue;
	  for(int j=0;j<mat_13[i].size();j++){
	    int a2 = mat_13[i][j];
	    if(a2==a1) continue;
	    if(cs.find(a2) == cs.end()) continue;

	    double D = 100000;
	    //Doit 
	    for(int p=0; p<peak.size();p++){
	      double d = dist(cs[a1],cs[a2],peak[p]);
	      d = min(d,dist(cs[a1],cs[a2],peak[p]));
	      if(d<D) D =d;
	    }

	    if(D>T&&D<0.5) T=D;
	    cout<<"Tol Pair "<<a1<<" "<<a2<<" "<<D<<" "<<T<<endl;
	  }
	}
	return T;
      }


      void assign_str(Protein & pro, double udl){
	vector<Peak2DFlat> flat;
	for(int p = 0; p<peak.size();p++){
	  if(peak[p].has_12()) continue;
	  if(peak[p].has_13()) continue;
	  //An peak
	  Peak2DFlat pp;
	  for(int i=0;i<peak[p].n1.size();i++){
	    for(int j=0;j<peak[p].n2.size();j++){
	      if(peak[p].n1[i]<pro.atom_size()&&
		 peak[p].n2[j]<pro.atom_size()){
		if(pro.distance(peak[p].n1[i],peak[p].n2[j])<udl){
		  //&&(!peak[p].is_12(peak[p].n1[i],peak[p].n2[j]))){
		  //is 12

		  pp.f.push_back(peak[p].n1[i]);
		  pp.s.push_back(peak[p].n2[j]);
		  pp.Dinv6.push_back(pow(pro.distance(peak[p].n1[i],peak[p].n2[j]),-6.0));
		  pp.W.push_back(1);
		  //		cout<<"YES"<<endl;
		}
	      }
	    }
	  }
	  if(pp.f.size()>0&&pp.s.size()>0){
	    pp.w1 = peak[p].w1;
	    pp.w2 = peak[p].w2;
	    pp.Vol = peak[p].Vol;
	    flat.push_back(pp);
	  
	  }
	
	}
	cout<<"Keeping "<<flat.size()<<" discarding "
	    <<peak.size()-flat.size()<<endl;

	double K;
	for(int t=0;t<5;t++){
	  //Get best K.
	  cout<<">> Calibrating....\n";
	  double Im = 0;
	  double Rm = 0;
	  double Dm = 0;
	  int count = 0;
	  for(int p=0;p<flat.size();p++){
	    flat[p].p_optimize(0);
	    double R6inv = 0;
	    for(int a=0;a<flat[p].f.size();a++){
	      double d = pro.distance(flat[p].f[a],flat[p].s[a]);
	      //	  if(d>7.5) continue;
	      R6inv += pow(d,-6.0);///flat[p].f.size();
		Dm += d;
		++count;
	    }
	    Im += flat[p].Vol*R6inv;
	    Rm += R6inv*R6inv;
	  
	  }
	
	  K = Im/Rm;
	  Dm = Dm/count;
	  Dm /=2;
	  cout<<">> K: "<<K<<" "<<Dm<<"\n";
	

	  for(int p=0;p<flat.size();p++){
	    if(flat[p].f.size()>1){
	      //	    cout<<"Opt peak "<<p+1<<endl;
	      //	    flat[p].p_optimize(K);
	    }
	  
	    double R6inv = 0;
	    for(int a=0;a<flat[p].f.size();a++){
	      R6inv += pow(pro.distance(flat[p].f[a],flat[p].s[a]),-6.0);
	    }
	    double E = (flat[p].Vol - K*R6inv);

	    cout<<">> Peak "<<p+1<<" "<<flat[p].Vol<<" "<<K*R6inv<<"\n";
	    for(int a=0;a<flat[p].f.size();a++){
	      double d = pro.distance(flat[p].f[a],flat[p].s[a]);
	      cout<<"\t"<<pro.atom_name(flat[p].f[a])<<" "
		  <<pro.atom_name(flat[p].s[a])<<" "<<d<<" "<< pow(K/(flat[p].Vol*flat[p].W[a]),1.0/6.0)<<" "
		  <<K*flat[p].W[a]*pow(d,-6.0)
		  <<" "
		  <<flat[p].W[a]<<" "
		  <<"\n";
	    }
	  }


	  {
	    vector<Peak2DFlat>::iterator iter = flat.begin();
	    while(iter!=flat.end()){
	      double R6inv = 0;
	      for(int a=0;a<iter->f.size();a++){
		R6inv += pow(pro.distance(iter->f[a],iter->s[a]),-6.0);
	      }
	      double E = abs(iter->Vol - K*R6inv)/min(K*R6inv,iter->Vol);
	      if(E>10){
		//	      cout<<"ERASING!!!"<<endl;
		iter = flat.erase(iter);
	      }
	      else{
		//	      cout<<"KEEP "<<E<<" "<<iter->Vol<<" "<<K*R6inv<<" "<<(iter->Vol - K*R6inv)/(K*R6inv)<<endl;
		++iter;
	      }
	    }
	  
	  }
	}

	ofstream out; out.open("opt.dat");
	for(int p=0;p<flat.size();p++){
	  // 	if(flat[p].f.size()>1){
	  // 	  out<<"Opt peak "<<p+1<<endl;
	  // 	  flat[p].p_optimize(K);
	  // 	}
	
	  double R6inv = 0;
	  for(int a=0;a<flat[p].f.size();a++){
	    R6inv += pow(pro.distance(flat[p].f[a],flat[p].s[a]),-6.0);
	  }
	  double E = (flat[p].Vol - K*R6inv);
	  out<<">> Peak "<<p+1<<" "<<flat[p].Vol<<" "<<K*R6inv<<" "<<abs(flat[p].Vol-K*R6inv)/min(flat[p].Vol,K*R6inv)<<"\n";
	  for(int a=0;a<flat[p].f.size();a++){
	    double d = pro.distance(flat[p].f[a],flat[p].s[a]);
	    out<<"\t"<<pro.atom_name(flat[p].f[a])<<" "
	       <<pro.atom_name(flat[p].s[a])<<" "<<d<<" "<< pow(K/(flat[p].Vol*flat[p].W[a]),1.0/6.0)<<" "
	       <<K*flat[p].W[a]*pow(d,-6.0)
	       <<" "
	       <<flat[p].W[a]<<" "
	       <<"\n";
	  }
	}
      }
    };


    void cluster(Sandy &s1, Sandy &s2, double tol, double tol2){
      vector<set<int> > s1_s;
      s1_s.resize(s1.peak.size());

      for(int p = 0; p< s1.peak.size();++p){
	for(int q = 0; q< s1.peak.size();q++){
	  if(p==q) continue;
	  if(dist(s1.peak[p],s1.peak[q])<tol)
	    s1_s[p].insert(q);
	}
      }


      vector<set<int> > s1_n;
      s1_n.resize(s1.peak.size());

      vector<set<int> > s1_n2;
      s1_n2.resize(s1.peak.size());

      for(int p = 0; p< s1.peak.size();++p){
	for(int q = 0; q< s2.peak.size();q++){
	  if(dist(s1.peak[p],s2.peak[q])<tol)
	    s1_n[p].insert(q);
	  if(dist(s1.peak[p],s2.peak[q])<tol2)
	    s1_n2[p].insert(q);
	}
      }

      for(int i=0;i<s1_n.size();i++){
	if(s1_n[i].size()&&s1_n[i].size()==s1_n2[i].size()){
	  cout<<">> Neighbors "<<i+1<<" "<<s1_s[i].size()<<" "
	      <<s1.peak[i].n_12()<<" "
	      <<s1.peak[i].n_13()<<" "
	      <<s1.peak[i].n_14()<<" "
	      <<s1.peak[i].Vol<<" "
	      <<" : ";
	  for(set<int>::iterator iter = s1_n[i].begin(); iter!= s1_n[i].end();
	      ++iter){
	    cout<<(*iter)+1<<" ";
	    int j = *iter;
	    cout<<s2.peak[j].n_12()<<" "
		<<s2.peak[j].n_13()<<" "
		<<s2.peak[j].n_14()<<" "
		<<s2.peak[j].Vol<<" ";
	  }
	  cout<<"\n";
	}
      }
    }
  }
}

using namespace Almost;

ALMTOSTRING(NMR::Sandy);
ALMTOSTRING(NMR::ChemicalShift);


extern "C" {
  void init_sandy(){
    //declarations here
    
    Module mod = Module("sandy","Solid state NMR");
		
    Class<NMR::Sandy>(mod.self(),"sandy")
      .def_method("read_peak",&NMR::Sandy::read_peak)
      .def_method("read_seq",&NMR::Sandy::read_seq)
      .def_method("read_shifts",&NMR::Sandy::read_shifts)
      .def_method("mkprotein",&NMR::Sandy::mkprotein)
      .def_method("assign_peaks",&NMR::Sandy::assign_peaks)
      .def_method("read_topos",&NMR::Sandy::read_topos)
      .def_method("challange_peaks",&NMR::Sandy::challange_peaks)
      .def_method("dump",&NMR::Sandy::dump)
      .def_method("read_dist",&NMR::Sandy::read_dist)
      .def_method("dump_symmetric",&NMR::Sandy::dump_symmetric)
      .def_method("na_bond",&NMR::Sandy::na_bond)
      .def_method("na_angle",&NMR::Sandy::na_angle)
      .def_method("normalize_na",&NMR::Sandy::normalize_na)
      .def_method("na_dihe",&NMR::Sandy::na_dihe)
      .def_method("na",&NMR::Sandy::na)
      .def_method("na_residue",&NMR::Sandy::na_residue)
      .def_method("bond_tol",&NMR::Sandy::bond_tol)
      .def_method("angle_tol",&NMR::Sandy::angle_tol)
      .def_method("assign_str",&NMR::Sandy::assign_str)
      ;


    Class<NMR::ChemicalShift>(mod.self(),"chemicalshift")
      .def_method("clear",&NMR::ChemicalShift::clear)
      .def_method("read",&NMR::ChemicalShift::read)
      .def_method("size",&NMR::ChemicalShift::size)
      ;

    mod
      .def_const("CHESHIRE",(int)NMR::ChemicalShift::CHESHIRE)	       
      .def_const("XEASY",(int)NMR::ChemicalShift::XEASY)
      ;

	      
    mod
      .def_function("cluster","Clusters peaks from different spectra",NMR::cluster);
  }
}