/***************************************************************************
                          topos.cpp  -  description
                             -------------------
    begin                : Tue Nov 18 01:45:35 2003
    copyright            : (C) 1999-2006 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 <fstream>
#include <iterator>
#include <iostream>
#include <vector>
#include <map>
#include <set>
#include <cmath>
#include <cheshire/topos.h>
#include <cheshire/topos_rc.h>
namespace Almost {

  namespace Topos {

    double T = 1.5;
    double kpos = 0;
    double ks = 1;
    //double zero[6] = { 0.51, 0.51, 4.56, 2.49, 2.01, 2.02};
    double zero[6] = { 0, 0, 0,0,0,0};

    double k[6][3] = { { 14.66, 17.54, 15.25},
		       { 15, 15, 15},
		       { 0.16, 0.18, 0.20},
		       { 0.72, 0.99, 0.72},
		       { 0.76, 0.91, 0.70},
		       { 1.15, 1.21, 1.04}
    };

    double ko[3] = { 0.74, 1.48, 0.74};

    double kk = 1;

    double k_homo[20][20] = {

      { 0, 1, 1, 1, 1, 1, 1, 2, 1, 2, 1, 1, 1, 2, 3, 1, 2, 2, 2, 2 },
      { 1, 0, 1, 1, 1, 1, 1, 2, 1, 2, 1, 0, 1, 1, 3, 1, 2, 1, 1, 2 },
      { 1, 1, 0, 0, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 3, 1, 2, 1, 1, 2 },
      { 1, 1, 0, 0, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 3, 1, 2, 1, 1, 2 },
      { 1, 1, 1, 1, 0, 1, 1, 2, 1, 2, 1, 1, 1, 1, 3, 1, 2, 1, 1, 2 },
      { 1, 1, 1, 1, 1, 0, 1, 2, 1, 2, 1, 1, 1, 1, 3, 1, 2, 1, 1, 2 },
      { 1, 1, 1, 1, 1, 1, 0, 2, 1, 2, 2, 2, 1, 1, 3, 1, 2, 1, 1, 2 },
      { 2, 2, 2, 2, 2, 2, 2, 0, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3 },
      { 1, 1, 1, 1, 1, 1, 1, 3, 0, 2, 1, 2, 2, 1, 3, 2, 2, 1, 1, 2 },
      { 2, 2, 2, 2, 2, 2, 2, 3, 2, 0, 1, 2, 2, 2, 3, 2, 1, 2, 2, 0 },
      { 1, 1, 1, 1, 1, 1, 2, 3, 1, 1, 0, 1, 1, 1, 3, 2, 2, 1, 1, 2 },
      { 1, 0, 1, 1, 1, 1, 2, 3, 2, 2, 1, 0, 1, 2, 3, 1, 2, 2, 2, 2 },
      { 1, 1, 1, 1, 1, 1, 1, 3, 2, 2, 1, 1, 0, 2, 3, 1, 2, 2, 2, 2 },
      { 2, 1, 1, 1, 1, 1, 1, 3, 1, 2, 1, 2, 2, 0, 3, 2, 2, 0, 0, 1 },
      { 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 0, 3, 3, 3, 3, 3 },
      { 1, 1, 1, 1, 1, 1, 1, 3, 2, 2, 2, 1, 1, 2, 3, 0, 1, 2, 2, 2 },
      { 2, 2, 2, 2, 2, 2, 2, 3, 2, 1, 2, 2, 2, 2, 3, 1, 0, 1, 1, 1 },
      { 2, 1, 1, 1, 1, 1, 1, 3, 1, 2, 1, 2, 2, 0, 3, 2, 1, 0, 0, 1 },
      { 2, 1, 1, 1, 1, 1, 1, 3, 1, 2, 1, 2, 2, 0, 3, 2, 1, 0, 0, 1 },
      { 2, 2, 2, 2, 2, 2, 2, 3, 2, 0, 2, 2, 2, 1, 3, 2, 1, 1, 1, 0 }

    };

    //double k_dpmf =0;

    vector<vector<double> > matx;
    vector<double> Ph,Pb,Pc;
    vector<ResidueData> res_data;
    vector<ResidueDataTriple> res_data_triple;
    vector<ResidueDataTriple> res_triple_db;

    void read_blast_mtx(ifstream & in){
      istream_iterator<string> iter(in),end;
      int size = atoi(iter->c_str());++iter;
      cout<<"Sequence "<<*iter<<"\n";++iter;
      for(int i=0;i<12;i++){
	cout<<*iter<<"\n";++iter;
      }
      matx.resize(size);
      char *ncbicodes = "XAXCDEFGHIKLMNPQRSTVWXYXXX";
      for(int i=0;i<size;i++){
	for(int p=0;p<26;p++){
	  if(ncbicodes[p]=='X'){
	    ++iter;
	  }
	  else {
	    double d = atof(iter->c_str());++iter;
	    d = -d/100;
	    matx[i].push_back(d);
	  }
	}
      }
    }

    void read_prob(istream & prob){
      istream_iterator<string> iter(prob),end;

      while(iter!=end){
	++iter;
	double w;
	w = atof(iter->c_str()); ++iter;
	if(w<0.0001) w = 10000;
	else {
	  w = -log(w);
	}
	Ph.push_back(w);

	w = atof(iter->c_str()); ++iter;
	if(w<0.0001) w = 10000;
	else {
	  w = -log(w);
	}
	Pb.push_back(w);

	w = atof(iter->c_str()); ++iter;
	if(w<0.0001) w = 10000;
	else {
	  w = -log(w);
	}
	Pc.push_back(w);

      }

    }

    int toind(char c){
      if(c=='A') return 0;
      if(c=='R') return 1;
      if(c=='N') return 2;
      if(c=='D') return 3;
      if(c=='C') return 4; 
      if(c=='Q') return 5; 
      if(c=='E') return 6; 
      if(c=='G') return 7; 
      if(c=='H') return 8; 
      if(c=='I') return 9; 
      if(c=='L') return 10; 
      if(c=='K') return 11; 
      if(c=='M') return 12; 
      if(c=='F') return 13; 
      if(c=='P') return 14; 
      if(c=='S') return 15; 
      if(c=='T') return 16; 
      if(c=='W') return 17; 
      if(c=='Y') return 18; 
      if(c=='V') return 19;

      return -1;

    }

    int toinda(char c){
      if(c=='A') return 0;
      if(c=='C') return 1;
      if(c=='D') return 2;
      if(c=='E') return 3;
      if(c=='F') return 4;
      if(c=='G') return 5;
      if(c=='H') return 6;
      if(c=='I') return 7;
      if(c=='K') return 8;
      if(c=='L') return 9;
      if(c=='M') return 10;
      if(c=='N') return 11;
      if(c=='P') return 12;
      if(c=='Q') return 13;
      if(c=='R') return 14;
      if(c=='S') return 15;
      if(c=='T') return 16;
      if(c=='V') return 17;
      if(c=='W') return 18;
      if(c=='Y') return 19;
      return -1;
    }


    void read_database(string name){
      ifstream in;
      in.open(name.c_str());
      if(!in){
	cout<<"Can't read database\n";
	exit(1);
      }

      int size;
      in.read((char*)&size,sizeof(int));
      cout<<"Reading "<<size<<" data\n";

      for(int i=0;i<size;i++){
	ResidueDataTriple dt;
	in.read((char*)&dt,sizeof(ResidueDataTriple));
	res_triple_db.push_back(dt);

      }
      cout<<"Read "<<res_triple_db.size()<<" data\n";
    }

    void read_data(ifstream & in, string name){
      istream_iterator<string> iter(in),end;

      while(iter!=end){
	ResidueData rd;
	strcpy(rd.prot, name.c_str());
	rd.num = atoi((*iter).c_str()); ++iter;
	rd.name = (*iter)[0]; ++iter;
	rd.phi = atof((*iter).c_str()); ++iter;
	rd.psi = atof((*iter).c_str()); ++iter;
	rd.ome = atof((*iter).c_str()); ++iter;
	rd.ha = atof((*iter).c_str()); ++iter;
	rd.h = atof((*iter).c_str()); ++iter;
	rd.n = atof((*iter).c_str()); ++iter;
	rd.ca = atof((*iter).c_str()); ++iter;
	rd.cb = atof((*iter).c_str()); ++iter;
	rd.c = atof((*iter).c_str()); ++iter;
	rd.ss = (*iter)[0]; ++iter;
	res_data.push_back(rd);
      }
      cout<<"Read protein "<<res_data.size()<<" data\n";
    }

    void bin2ascii(string bin, string ascii){
      ifstream in;
      in.open(bin.c_str());
      if(!in){
	cout<<"Can't read database "+bin+"\n";
	exit(1);
      }

      ofstream out;
      out.open(ascii.c_str());
      if(!out){
	cout<<"Can't read database "+ascii+"\n";
	exit(1);
      }

      int size;
      in.read((char*)&size,sizeof(int));
      out<<size<<"\n";
      cout<<"Reading "<<size<<" data\n";

      for(int i=0;i<size;i++){
	ResidueDataTriple dt;
	in.read((char*)&dt,sizeof(ResidueDataTriple));
	out<<dt.prot<<" "
	   <<dt.num[0]<<" "<<dt.num[1]<<" "<<dt.num[2]<<" "
	   <<dt.name<<" "
	   <<dt.ss[0]<<" "<<dt.ss[1]<<" "<<dt.ss[2]<<" "
	   <<dt.phi[0]<<" "<<dt.phi[1]<<" "<<dt.phi[2]<<" "
	   <<dt.psi[0]<<" "<<dt.psi[1]<<" "<<dt.psi[2]<<" "
	   <<dt.ome[0]<<" "<<dt.ome[1]<<" "<<dt.ome[2]<<" "
	   <<dt.ha[0]<<" "<<dt.ha[1]<<" "<<dt.ha[2]<<" "
	   <<dt.h[0]<<" "<<dt.h[1]<<" "<<dt.h[2]<<" "
	   <<dt.n[0]<<" "<<dt.n[1]<<" "<<dt.n[2]<<" "
	   <<dt.ca[0]<<" "<<dt.ca[1]<<" "<<dt.ca[2]<<" "
	   <<dt.cb[0]<<" "<<dt.cb[1]<<" "<<dt.cb[2]<<" "
	   <<dt.c[0]<<" "<<dt.c[1]<<" "<<dt.c[2]<<" "
	   <<"\n";
	//      res_triple_db.push_back(dt);

      }
    }

    void read_ascii_database(string name){
      ifstream in;
      in.open(name.c_str());
      if(!in){
	cout<<"Can't read database\n";
	exit(1);
      }

      istream_iterator<string> iter(in), end;
      int size = atoi(iter->c_str()); ++iter;
      cout<<"Reading "<<size<<" data\n";

      for(int i=0;i<size;i++){
	ResidueDataTriple dt;
	//      in.read((char*)&dt,sizeof(ResidueDataTriple));
	string tok;
	tok = (*iter); ++iter;

	strcpy(dt.prot,tok.c_str());

	for(int q = 0;q<3;q++){
	  dt.num[q] = atoi(iter->c_str()); ++iter;
	}

	tok = (*iter); ++iter;
	strcpy(dt.name,tok.c_str());

	for(int q = 0;q<3;q++){
	  dt.ss[q] = iter->c_str()[0]; ++iter;
	}


	for(int q = 0;q<3;q++){
	  dt.phi[q] = atof(iter->c_str()); ++iter;
	}

	for(int q = 0;q<3;q++){
	  dt.psi[q] = atof(iter->c_str()); ++iter;
	}

	for(int q = 0;q<3;q++){
	  dt.ome[q] = atof(iter->c_str()); ++iter;
	}

	for(int q = 0;q<3;q++){
	  dt.ha[q] = atof(iter->c_str()); ++iter;
	}

	for(int q = 0;q<3;q++){
	  dt.h[q] = atof(iter->c_str()); ++iter;
	}

	for(int q = 0;q<3;q++){
	  dt.n[q] = atof(iter->c_str()); ++iter;
	}

	for(int q = 0;q<3;q++){
	  dt.ca[q] = atof(iter->c_str()); ++iter;
	}

	for(int q = 0;q<3;q++){
	  dt.cb[q] = atof(iter->c_str()); ++iter;
	}

	for(int q = 0;q<3;q++){
	  dt.c[q] = atof(iter->c_str()); ++iter;
	}

	res_triple_db.push_back(dt);

      }
      cout<<"Read "<<res_triple_db.size()<<" data\n";
    }

    void read_ascii_database_with_observer(string name, ToposObserver & obs){
      ifstream in;
      in.open(name.c_str());
      if(!in){
	cout<<"Can't read database\n";
	exit(1);
      }

      istream_iterator<string> iter(in), end;
      int size = atoi(iter->c_str()); ++iter;
      //    cout<<"Reading "<<size<<" data\n";

      for(int i=0;i<size;i++){
	if(i%10000)
	  obs.on_db_read(i,size);
	ResidueDataTriple dt;
	//      in.read((char*)&dt,sizeof(ResidueDataTriple));
	string tok;
	tok = (*iter); ++iter;

	strcpy(dt.prot,tok.c_str());

	for(int q = 0;q<3;q++){
	  dt.num[q] = atoi(iter->c_str()); ++iter;
	}

	tok = (*iter); ++iter;
	strcpy(dt.name,tok.c_str());

	for(int q = 0;q<3;q++){
	  dt.ss[q] = iter->c_str()[0]; ++iter;
	}


	for(int q = 0;q<3;q++){
	  dt.phi[q] = atof(iter->c_str()); ++iter;
	}

	for(int q = 0;q<3;q++){
	  dt.psi[q] = atof(iter->c_str()); ++iter;
	}

	for(int q = 0;q<3;q++){
	  dt.ome[q] = atof(iter->c_str()); ++iter;
	}

	for(int q = 0;q<3;q++){
	  dt.ha[q] = atof(iter->c_str()); ++iter;
	}

	for(int q = 0;q<3;q++){
	  dt.h[q] = atof(iter->c_str()); ++iter;
	}

	for(int q = 0;q<3;q++){
	  dt.n[q] = atof(iter->c_str()); ++iter;
	}

	for(int q = 0;q<3;q++){
	  dt.ca[q] = atof(iter->c_str()); ++iter;
	}

	for(int q = 0;q<3;q++){
	  dt.cb[q] = atof(iter->c_str()); ++iter;
	}

	for(int q = 0;q<3;q++){
	  dt.c[q] = atof(iter->c_str()); ++iter;
	}

	res_triple_db.push_back(dt);

      }

      obs.on_db_read(size,size);
      //    cout<<"Read "<<res_triple_db.size()<<" data\n";
    }

    void mktriple(){
      int i=0;
      int size = res_data.size();
      while((i+2)<size){
	ResidueDataTriple dt;
	strcpy(dt.prot,res_data[i].prot);
	dt.num[0] = res_data[i].num;
	dt.num[1] = res_data[i+1].num;
	dt.num[2] = res_data[i+2].num;

	dt.name[0] = res_data[i].name;
	dt.name[1] = res_data[i+1].name;
	dt.name[2] = res_data[i+2].name;
	dt.name[3] = '\0';

	for(int j=0;j<3;j++)
	  dt.ss[j] = res_data[i+j].ss;

	dt.phi[0] = res_data[i].phi;
	dt.phi[1] = res_data[i+1].phi;
	dt.phi[2] = res_data[i+2].phi;

	dt.psi[0] = res_data[i].psi;
	dt.psi[1] = res_data[i+1].psi;
	dt.psi[2] = res_data[i+2].psi;

	dt.ome[0] = res_data[i].ome;
	dt.ome[1] = res_data[i+1].ome;
	dt.ome[2] = res_data[i+2].ome;

	dt.ha[0] = res_data[i].ha;
	dt.ha[1] = res_data[i+1].ha;
	dt.ha[2] = res_data[i+2].ha;

	dt.h[0] = res_data[i].h;
	dt.h[1] = res_data[i+1].h;
	dt.h[2] = res_data[i+2].h;

	dt.n[0] = res_data[i].n;
	dt.n[1] = res_data[i+1].n;
	dt.n[2] = res_data[i+2].n;

	dt.ca[0] = res_data[i].ca;
	dt.ca[1] = res_data[i+1].ca;
	dt.ca[2] = res_data[i+2].ca;

	dt.cb[0] = res_data[i].cb;
	dt.cb[1] = res_data[i+1].cb;
	dt.cb[2] = res_data[i+2].cb;

	dt.c[0] = res_data[i].c;
	dt.c[1] = res_data[i+1].c;
	dt.c[2] = res_data[i+2].c;

	res_data_triple.push_back(dt);
	i = i+1;
      }
    }

    double score(ResidueDataTriple & a,
		 ResidueDataTriple & b){


      for(int i=0;i<3;i++){
	if((toind(a.name[i])==-1)||(toind(b.name[i])==-1)){
	  cout<<a.name<<" "<<b.name<<"\n";
	  cout<<a.prot<<" "<<b.prot<<"\n";
	  cout<<"NAME ERROR"<<endl;
	}
      }
  

      double s = 0;

      for(int i=0;i<3;i++){
	double d = k_homo[toind(a.name[i])][toind(b.name[i])];
	s += ko[i]*d*d;
      }

      //     //pos score
      //     for(int i=0;i<3;i++){
      //       int apos = a.num[i]-1;
      //       int bidx = toinda(b.name[i]);
      //       s -= kpos*matx[apos][bidx];

      //     }


      //  if(strcmp(a.name,b.name)!=0) return 1000000;

      for(int i=0;i<3;i++){
	int ii,jj;
	ii = toinda(a.name[i]);
	jj = toinda(b.name[i]);
	double d1,d2;
	if(a.ha[i]==0){
	  d1 = 1.5*zero[0];
	}
	else d1 = (a.ha[i] -HARC[ii]);

	if(b.ha[i]==0){
	  d2 = 1.5*zero[0];
	}
	else d2 = (b.ha[i]-HARC[jj]);

	double d = (d1-d2);
	s += kk*k[0][i]*d*d;
      }

      //    for(int i=0;i<3;i++){
      //      int ii,jj;
      //      ii = toinda(a.name[i]);
      //      jj = toinda(b.name[i]);
      //      double d1,d2;
      //      if(a.h[i]==0){
      //        d1 = 1.5*zero[1];
      //      }
      //      else d1 = (a.h[i]-HNRC[ii]);
     
      //      if(b.h[i]==0){
      //        d2 = 1.5*zero[1];
      //      }
      //      else d2 = (b.h[i]-HNRC[jj]);
     
      //      double d = (d1-d2);
      //      s += kk*k[1][i]*d*d;
     
      //    }

      for(int i=0;i<3;i++){
	int ii,jj;
	ii = toinda(a.name[i]);
	jj = toinda(b.name[i]);
	double d1,d2;
	if(a.n[i]==0){
	  d1 = 1.5*zero[2];
	}
	else d1 = (a.n[i]-N15RC[ii]);
    
	if(b.n[i]==0){
	  d2 = 1.5*zero[2];
	}
	else d2 = (b.n[i]-N15RC[jj]);
    
	double d = (d1-d2);
	s += kk*k[2][i]*d*d;
      }

      for(int i=0;i<3;i++){
	int ii,jj;
	ii = toinda(a.name[i]);
	jj = toinda(b.name[i]);
	double d1,d2;
	if(a.ca[i]==0){
	  d1 = 1.5*zero[3];
	}
	else d1 = (a.ca[i]-CARC[ii]);
    
	if(b.ca[i]==0){
	  d2 = 1.5*zero[3];
	}
	else d2 = (b.ca[i]-CARC[jj]);
    
	double d = (d1-d2);
	s += kk*k[3][i]*d*d;

      }

      for(int i=0;i<3;i++){
	int ii,jj;
	ii = toinda(a.name[i]);
	jj = toinda(b.name[i]);
	double d1,d2;
	if((a.cb[i]==0)&&(a.name[i]!='G')&&(b.name[i]!='G')){
	  d1 = 1.5*zero[4];
	}
	else d1 = (a.cb[i]-CBRC[ii]);//best was 1.0
    
	if((b.cb[i]==0)&&(a.name[i]!='G')&&(b.name[i]!='G')){
	  d2 = 1.5*zero[4];
	}
	else d2 = (b.cb[i]-CBRC[jj]);
    
	double d = (d1-d2);
	s += kk*k[4][i]*d*d;
      }

      int pos = a.num[1]-1;
      char c = b.ss[1];
      if((c=='H')||(c=='G')||(c=='I')) s += ks*Ph[pos];
      else if((c=='B')||(c=='E')||(c=='b')) s+= ks*Pb[pos];
      else s += ks*Pc[pos];

      //    double h = (64.0/360.0);
      //    double fpsi,fphi;
      //    fphi = floor((b.phi[1]+180.0)*h);
      //    int iphi = (int)fphi;

      //    fpsi = floor((b.psi[1]+180.0)*h);
      //    int ipsi = (int)fpsi;
      //    string    aa;
      //    aa+= b.name[1];
      //   cout<<"DEB "<<k_dpmf*H[aa].G[iphi][ipsi]<<endl;
      //    s += k_dpmf*H[aa].G[iphi][ipsi];
      //    for(int i=0;i<3;i++){
      //      int ii,jj;
      //      ii = toinda(a.name[i]);
      //      jj = toinda(b.name[i]);
      //     double d1,d2;
      //     if(a.c[i]==0){
      //       d1 = 1.5*zero[5];
      //     }
      //     else d1 = (a.c[i]-CORC[ii]);
    
      //     if(b.c[i]==0){
      //       d2 = 1.5*zero[5];
      //     }
      //     else d2 = (b.c[i]-CORC[jj]);
    
      //     double d = (d1-d2);
      //     s += kk*k[5][i]*d*d;

      //    }

      return s;
    }

    double ramadist(double phi1, double psi1, 
		    double phi2, double psi2){
      double d1 = abs(phi1-phi2);
      if(d1>180) d1 = d1-180;
      double d2 = abs(psi1-psi2);
      if(d2>180) d2 = d2-180;
      //   cout<<"DEBUG "<< phi1<<" "<<psi1<<" "<<d1<<" "
      //       << phi2<<" "<<psi2<<" "<<d2<<"\n";
      return sqrt(d1*d1+d2*d2);

    }

    double extscore(ResidueDataTriple & a,
		    ResidueDataTriple & b){


      for(int i=0;i<3;i++){
	if((toind(a.name[i])==-1)||(toind(b.name[i])==-1)){
	  cout<<a.name<<" "<<b.name<<"\n";
	  cout<<a.prot<<" "<<b.prot<<"\n";
	}
      }
  

      double s = 0;

      for(int i=0;i<3;i++){
	double d = k_homo[toind(a.name[i])][toind(b.name[i])];
	s += ko[i]*d*d;
      }

      //  if(strcmp(a.name,b.name)!=0) return 1000000;

      for(int i=0;i<3;i++){
	int ii,jj;
	ii = toinda(a.name[i]);
	jj = toinda(b.name[i]);
	double d = (a.ha[i]-HARC[ii])-(b.ha[i]-HARC[jj]);
	s += kk*k[0][i]*d*d;
      }

      for(int i=0;i<3;i++){
	int ii,jj;
	ii = toinda(a.name[i]);
	jj = toinda(b.name[i]);
	double d = (a.h[i]-HNRC[ii])-(b.h[i]-HNRC[jj]);
	s += kk*k[1][i]*d*d;
      }

      for(int i=0;i<3;i++){
	int ii,jj;
	ii = toinda(a.name[i]);
	jj = toinda(b.name[i]);
	double d = (a.n[i]-N15RC[ii])-(b.n[i]-N15RC[jj]);
	s += kk*k[2][i]*d*d;
      }

      for(int i=0;i<3;i++){
	int ii,jj;
	ii = toinda(a.name[i]);
	jj = toinda(b.name[i]);
	double d = (a.ca[i]-CARC[ii])-(b.ca[i]-CARC[jj]);
	s += kk*k[3][i]*d*d;
      }

      for(int i=0;i<3;i++){
	int ii,jj;
	ii = toinda(a.name[i]);
	jj = toinda(b.name[i]);
	if((a.name[i]!='G')&&(b.name[i]!='G')){
	  double d = (a.cb[i]-CBRC[ii])-(b.cb[i]-CBRC[jj]);
	  s += kk*k[4][i]*d*d;
	}
      }

      for(int i=0;i<3;i++){
	int ii,jj;
	ii = toinda(a.name[i]);
	jj = toinda(b.name[i]);
	double d = (a.c[i]-CORC[ii])-(b.c[i]-CORC[jj]);
	s += kk*k[5][i]*d*d;
      }


      double r = ramadist(a.phi[1],a.psi[1],b.phi[1],b.psi[1]);
      s += r*r;
      return s;
    }

    pair<double,double> cluster(vector<double> & phi,
				vector<double> & psi, vector<double> &scr){
      pair<double,double> phi_psi;
      vector<vector<int> > proxy;
      proxy.resize(phi.size());
      for(int i=0;i<proxy.size();i++)
	proxy[i].resize(phi.size());

      for(int i=0;i<phi.size();i++){
	for(int j=i+1;j<phi.size();j++){
	  double r = ramadist(phi[i],psi[i],phi[j],psi[j]);
	  if(r<25){
	    ++proxy[i][i];
	    ++proxy[j][j];
	    proxy[i][j] =  proxy[j][i] = 1;
	  }
	  else {
	    proxy[i][j] =  proxy[j][i] = 0;
	  }
	}
      }
	       

      {
	double bscore=10000;
	double bphi,bpsi;
	int bsize;
    
	for(int i=0;i<phi.size();i++){
	  int p = i;
	  int n = proxy[i][i];
	  if(n>0){
	    double pp = 0,qq=0, s = 0;
	    for(int i=0;i<phi.size();i++){
	      if((i!=p)&&(proxy[p][i]==1)){
		pp += phi[i];
		qq += psi[i];
		s += scr[i];
	      }
	    }
	    cout<<"\t>> Cluster  "<<pp/n<<" "<<qq/n<<" "<<n<<" "<<s/n<<"\n";
	    if(s/n<bscore){
	      bphi = pp/n;
	      bpsi = qq/n;
	      bscore = s/n;
	      bsize = n;
	    }
	  }
	}
	cout<<"\t>> Cluster Best Score  "<<bphi<<" "<<bpsi<<" "<<bsize<<" "<<bscore<<"\n";
      }


      //biggst cluster
      int p = -1;
      int  n = -1; 
      for(int i=0;i<phi.size();i++){
	if(proxy[i][i]>n){
	  p=i;
	  n=proxy[i][i];
	}
      }
  
      double pp = 0,qq=0;
      for(int i=0;i<phi.size();i++){
	if((i!=p)&&(proxy[p][i]==1)){
	  pp += phi[i];
	  qq += psi[i];
	}
      }

      if(n>2){
	cout<<"\n\t>> CLUSTER  "<<pp/n<<" ("<<qq/n<<") "<<n<<"\n";
	phi_psi.first = pp/n;
	phi_psi.second= qq/n;
      }
      else {
	cout<<"\n\t>> Cluster  "<<pp/n<<" ("<<qq/n<<") "<<n<<"\n";
	phi_psi.first = 999;
	phi_psi.second= qq/n;
      }
 
      return phi_psi;
    }

    double cluster_weight(ResidueDataTriple &a, vector<ResidueDataTriple> &b,
			  double S,pair<double,double> &data){

      double Stot = 0;
      double phi = 0;
      double psi = 0;
      double dphi = 0;
      double dpsi = 0;
      for(int i=0;i<b.size();i++){
	double s = score(a,b[i]);    
	phi += b[i].phi[1]*exp(-s/T); 
	psi += b[i].psi[1]*exp(-s/T); 
	dphi +=b[i].phi[1]*b[i].phi[1]*exp(-s/T);
	dpsi +=b[i].psi[1]*b[i].psi[1]*exp(-s/T);
	Stot += exp(-s/T);
      }
      if(b.size()>1){
	dphi = sqrt(dphi/Stot - phi/Stot*phi/Stot);
	dpsi = sqrt(dpsi/Stot - psi/Stot*psi/Stot);
	//    cout<<phi/Stot<<" "<<psi/Stot<<" "<<dphi<<" "<<dpsi<<" "<<Stot/S;
	data.first = phi/Stot;
	data.second = psi/Stot;
      }
      else {
	//    cout<<phi/Stot<<" "<<psi/Stot<<" "<<-1<<" "<<-1<<" "<<Stot/S;
	data.first = phi/Stot;
	data.second = psi/Stot;
      }



      return Stot/S;
    } 

    void mkprediction(ofstream & log, ofstream & rest){
      for(int i=0;i<res_data_triple.size();i++){
	vector<ResidueDataTriple> top;
	vector<vector<ResidueDataTriple> > clusters;
	vector<double> phi;
	vector<double> psi;

	vector<double> scr;
	clusters.push_back(vector<ResidueDataTriple>());
	map<double,ResidueDataTriple > res;
	set<double> scores;
	for(int j=0;j<res_triple_db.size();j++){
	  double s = score(res_data_triple[i],res_triple_db[j]);
	  if(s<15){//WAS 50
	    res[s] = res_triple_db[j];
	    scores.insert(s);
	  }
	}
	cerr<<"Triple "<<i+1<<" "<<res.size()<<endl;
	set<double>::iterator iter = scores.begin(),end = scores.end();
	int count = 0;
	int n = 0;
	double PHI = 0;
	double PSI = 0;
	double DPHI = 0;
	double DPSI = 0;
	double Stot2 = 0;
	while((iter!=end)&&(count<500)&&(n<500)){

	  ++n;
	  PHI += res[*iter].phi[1]*exp(-(*iter)/T);
	  PSI += res[*iter].psi[1]*exp(-(*iter)/T);
	  DPHI += ( res[*iter].phi[1])*( res[*iter].phi[1]*exp(-(*iter)/T));
	  DPSI += ( res[*iter].psi[1])*( res[*iter].psi[1]*exp(-(*iter)/T));
	  Stot2 +=exp(-(*iter)/T);
	  top.push_back(res[*iter]);
	  double S = 10000;
	  for(int l = 0;l<phi.size();l++){
	    double r = ramadist(phi[l],psi[l],res[*iter].phi[1],
				res[*iter].psi[1]);
	    if(r<S) S=r;
	    if(r<15){ //Was 25 Wed Nov 29 18:08:44 CET 2006
	      clusters[l].push_back(res[*iter]);
	      break;
	    }
	  }
	  if(S>90){ //Was 50
	    phi.push_back(res[*iter].phi[1]);
	    psi.push_back(res[*iter].psi[1]);
	    scr.push_back(*iter);
	    clusters.push_back(vector<ResidueDataTriple>());
	    clusters.back().push_back(res[*iter]);
	    ++iter; 
	  }
	  else{ 
	    ++iter;	
	  }
	  ++count;
	}

	if(n>0){
	  PHI /= Stot2;
	  PSI /= Stot2;
	  DPHI = sqrt(DPHI/Stot2-PHI*PHI);
	  DPSI = sqrt(DPSI/Stot2-PSI*PSI);


	  map<double,int> rank;
	  vector<pair<double,double> > val;
	  for(int l =0;l<clusters.size();l++){
	    pair<double,double> tmp;
	    double w = cluster_weight(res_data_triple[i],clusters[l],Stot2,tmp);
	    if(w>0.00005){
	      val.push_back(tmp);
	      rank[w] = l;
	    }
	  }
	  log<<res_data_triple[i].num[1]<<" "<<res_data_triple[i].name[1]<<" "<<PHI<<" "<<PSI<<" "<<rank.size()<<" ";
	  log<<res_data_triple[i].phi[1]<<" "<<res_data_triple[i].psi[1]<<" "
	     <<sqrt(DPHI)<<" "<<sqrt(DPSI)<<" "<<res_data_triple[i].ss[1]
	     <<" XXX\n";
	  map<double,int>::reverse_iterator iter = rank.rbegin(),
	    end = rank.rend();
	  int count = 0;
	  while(iter!=end){
	    int l = iter->second;
	    log<<"\tCluster "<<clusters[l][0].prot<<" "
	       << clusters[l][0].name<<" "<<clusters[l][0].num[0]<<"--"
	       <<clusters[l][0].num[2]<<" "
	       << clusters[l][0].ss[0]
	       << clusters[l][0].ss[1]
	       << clusters[l][0].ss[2];
	    log<<" "<<clusters[l].size()<<" "
	       <<val[l].first<<" "<<val[l].second<<" "
	       <<iter->first;
	    ++iter;
	    log<<"\n";
	    if(count<3){
	      rest <<res_data_triple[i].num[1]<<" "
		   <<res_data_triple[i].name[1]<<" "
		   <<val[l].first<<" "<<val[l].second<<"\n";
	    }
	    ++count;	
	  }
	  log<<"\n";
	}
	else {
	  log<<res_data_triple[i].num[1]<<" "<<res_data_triple[i].name[1]<<" "<<999<<" "<<999<<" "<<psi.size()<<" ";
	  log<<360<<" "<<360<<"\n";
	}
      }
    }

    void mkprediction_with_observer(ToposObserver & obs, ofstream & log, ofstream & rest){
      for(int i=0;i<res_data_triple.size();i++){
	obs.on_sequence(i,res_data_triple.size()-1);
	vector<ResidueDataTriple> top;
	vector<vector<ResidueDataTriple> > clusters;
	vector<double> phi;
	vector<double> psi;

	vector<double> scr;
	clusters.push_back(vector<ResidueDataTriple>());
	map<double,ResidueDataTriple > res;
	set<double> scores;
	for(int j=0;j<res_triple_db.size();j++){
	  double s = score(res_data_triple[i],res_triple_db[j]);
	  if(s<15){//WAS 50
	    res[s] = res_triple_db[j];
	    scores.insert(s);
	  }
	}
	//      cerr<<"Triple "<<i+1<<" "<<res.size()<<endl;
	set<double>::iterator iter = scores.begin(),end = scores.end();
	int count = 0;
	int n = 0;
	double PHI = 0;
	double PSI = 0;
	double DPHI = 0;
	double DPSI = 0;
	double Stot2 = 0;
	while((iter!=end)&&(count<500)&&(n<500)){

	  ++n;
	  PHI += res[*iter].phi[1]*exp(-(*iter)/T);
	  PSI += res[*iter].psi[1]*exp(-(*iter)/T);
	  DPHI += ( res[*iter].phi[1])*( res[*iter].phi[1]*exp(-(*iter)/T));
	  DPSI += ( res[*iter].psi[1])*( res[*iter].psi[1]*exp(-(*iter)/T));
	  Stot2 +=exp(-(*iter)/T);
	  top.push_back(res[*iter]);
	  double S = 10000;
	  for(int l = 0;l<phi.size();l++){
	    double r = ramadist(phi[l],psi[l],res[*iter].phi[1],
				res[*iter].psi[1]);
	    if(r<S) S=r;
	    if(r<15){ //Was 25 Wed Nov 29 18:08:44 CET 2006
	      clusters[l].push_back(res[*iter]);
	      break;
	    }
	  }
	  if(S>90){ //Was 50
	    phi.push_back(res[*iter].phi[1]);
	    psi.push_back(res[*iter].psi[1]);
	    scr.push_back(*iter);
	    clusters.push_back(vector<ResidueDataTriple>());
	    clusters.back().push_back(res[*iter]);
	    ++iter; 
	  }
	  else{ 
	    ++iter;	
	  }
	  ++count;
	}

	if(n>0){
	  PHI /= Stot2;
	  PSI /= Stot2;
	  DPHI = sqrt(DPHI/Stot2-PHI*PHI);
	  DPSI = sqrt(DPSI/Stot2-PSI*PSI);


	  map<double,int> rank;
	  vector<pair<double,double> > val;
	  for(int l =0;l<clusters.size();l++){
	    pair<double,double> tmp;
	    double w = cluster_weight(res_data_triple[i],clusters[l],Stot2,tmp);
	    if(w>0.00005){
	      val.push_back(tmp);
	      rank[w] = l;
	    }
	  }
	  log<<res_data_triple[i].num[1]<<" "<<res_data_triple[i].name[1]<<" "<<PHI<<" "<<PSI<<" "<<rank.size()<<" ";
	  log<<res_data_triple[i].phi[1]<<" "<<res_data_triple[i].psi[1]<<" "
	     <<sqrt(DPHI)<<" "<<sqrt(DPSI)<<" "<<res_data_triple[i].ss[1]
	     <<" XXX\n";
	  map<double,int>::reverse_iterator iter = rank.rbegin(),
	    end = rank.rend();
	  int count = 0;
	  while(iter!=end){
	    int l = iter->second;
	    log<<"\tCluster "<<clusters[l][0].prot<<" "
	       << clusters[l][0].name<<" "<<clusters[l][0].num[0]<<"--"
	       <<clusters[l][0].num[2]<<" "
	       << clusters[l][0].ss[0]
	       << clusters[l][0].ss[1]
	       << clusters[l][0].ss[2];
	    log<<" "<<clusters[l].size()<<" "
	       <<val[l].first<<" "<<val[l].second<<" "
	       <<iter->first;
	    ++iter;
	    log<<"\n";
	    if(count<3){
	      rest <<res_data_triple[i].num[1]<<" "
		   <<res_data_triple[i].name[1]<<" "
		   <<val[l].first<<" "<<val[l].second<<"\n";
	    }
	    ++count;
	  }
	  log<<"\n";
	}
	else {
	  log<<res_data_triple[i].num[1]<<" "<<res_data_triple[i].name[1]<<" "<<999<<" "<<999<<" "<<psi.size()<<" ";
	  log<<360<<" "<<360<<"\n";
	}
      }
    }

    void topos(string mrg, string prob,  double Ks, double Kk,
	       string log_file,
	       string rest_file)  
	     
    {
      cout<<"Reading database\n";
      //    read_database();

      ifstream in;
      in.open(mrg.c_str());
      {
	ifstream in;
	in.open(prob.c_str());
	read_prob(in);
      }
      read_data(in,string(mrg.begin(),mrg.end()-4));
      ks = Ks;
      cerr<<"KS set to "<<ks<<endl;

      kk = Kk;


      mktriple();
      ofstream log;
      log.open(log_file.c_str());
      ofstream rest;
      rest.open(rest_file.c_str());
      mkprediction(log,rest);
    }

    void topos_with_observer(string mrg, string prob,  double Ks, double Kk,
			     string log_file,
			     string rest_file,
			     ToposObserver & obs)  
	     
    {
      //    cout<<"Reading database\n";
      //    read_database();

      ifstream in;
      in.open(mrg.c_str());
      {
	ifstream in;
	in.open(prob.c_str());
	read_prob(in);
      }
      read_data(in,string(mrg.begin(),mrg.end()-4));
      ks = Ks;
      //    cerr<<"KS set to "<<ks<<endl;

      kk = Kk;


      mktriple();
      ofstream log;
      log.open(log_file.c_str());
      ofstream rest;
      rest.open(rest_file.c_str());
      mkprediction_with_observer(obs,log,rest);
    }

    void maketriple(string mrg, string trp, string name){
      //Read mrg
      ifstream in;
      in.open(mrg.c_str());
      istream_iterator<string> iter(in),end;
      vector<ResidueData> res;
      vector<ResidueDataTriple> res_triple;
      while(iter!=end){
	ResidueData rd;
	strcpy(rd.prot, name.c_str());
	rd.num = atoi((*iter).c_str()); ++iter;
	rd.name = toupper((*iter)[0]); ++iter;
	rd.phi = atof((*iter).c_str()); ++iter;
	rd.psi = atof((*iter).c_str()); ++iter;
	rd.ome = atof((*iter).c_str()); ++iter;
	rd.ha = atof((*iter).c_str()); ++iter;
	rd.h = atof((*iter).c_str()); ++iter;
	rd.n = atof((*iter).c_str()); ++iter;
	rd.ca = atof((*iter).c_str()); ++iter;
	rd.cb = atof((*iter).c_str()); ++iter;
	rd.c = atof((*iter).c_str()); ++iter;
	rd.ss = (*iter)[0]; ++iter;
	res.push_back(rd);
      }
      cout<<"Read protein "<<res.size()<<" data\n";
      //Mk triple
      int i=0;
      int size = res.size();
      while((i+2)<size){
	ResidueDataTriple dt;
	strcpy(dt.prot,res[i].prot);
	dt.num[0] = res[i].num;
	dt.num[1] = res[i+1].num;
	dt.num[2] = res[i+2].num;

	dt.name[0] = res[i].name;
	dt.name[1] = res[i+1].name;
	dt.name[2] = res[i+2].name;
	dt.name[3] = '\0';

	for(int j=0;j<3;j++)
	  dt.ss[j] = res[i+j].ss;

	dt.phi[0] = res[i].phi;
	dt.phi[1] = res[i+1].phi;
	dt.phi[2] = res[i+2].phi;

	dt.psi[0] = res[i].psi;
	dt.psi[1] = res[i+1].psi;
	dt.psi[2] = res[i+2].psi;

	dt.ome[0] = res[i].ome;
	dt.ome[1] = res[i+1].ome;
	dt.ome[2] = res[i+2].ome;

	dt.ha[0] = res[i].ha;
	dt.ha[1] = res[i+1].ha;
	dt.ha[2] = res[i+2].ha;

	dt.h[0] = res[i].h;
	dt.h[1] = res[i+1].h;
	dt.h[2] = res[i+2].h;

	dt.n[0] = res[i].n;
	dt.n[1] = res[i+1].n;
	dt.n[2] = res[i+2].n;

	dt.ca[0] = res[i].ca;
	dt.ca[1] = res[i+1].ca;
	dt.ca[2] = res[i+2].ca;

	dt.cb[0] = res[i].cb;
	dt.cb[1] = res[i+1].cb;
	dt.cb[2] = res[i+2].cb;

	dt.c[0] = res[i].c;
	dt.c[1] = res[i+1].c;
	dt.c[2] = res[i+2].c;

	res_triple.push_back(dt);
	i = i+1;
      }

      size  = res_triple.size();
      ofstream out; out.open(trp.c_str());
      for(int i=0;i<size;i++){
	ResidueDataTriple &dt = res_triple[i];
	out<<dt.prot<<" "
	   <<dt.num[0]<<" "<<dt.num[1]<<" "<<dt.num[2]<<" "
	   <<dt.name<<" "
	   <<dt.ss[0]<<" "<<dt.ss[1]<<" "<<dt.ss[2]<<" "
	   <<dt.phi[0]<<" "<<dt.phi[1]<<" "<<dt.phi[2]<<" "
	   <<dt.psi[0]<<" "<<dt.psi[1]<<" "<<dt.psi[2]<<" "
	   <<dt.ome[0]<<" "<<dt.ome[1]<<" "<<dt.ome[2]<<" "
	   <<dt.ha[0]<<" "<<dt.ha[1]<<" "<<dt.ha[2]<<" "
	   <<dt.h[0]<<" "<<dt.h[1]<<" "<<dt.h[2]<<" "
	   <<dt.n[0]<<" "<<dt.n[1]<<" "<<dt.n[2]<<" "
	   <<dt.ca[0]<<" "<<dt.ca[1]<<" "<<dt.ca[2]<<" "
	   <<dt.cb[0]<<" "<<dt.cb[1]<<" "<<dt.cb[2]<<" "
	   <<dt.c[0]<<" "<<dt.c[1]<<" "<<dt.c[2]<<" "
	   <<"\n";
	//      res_triple_db.push_back(dt);

      }    
    }

    const vector<ResidueDataTriple> & database(){
      return res_triple_db;
    }

    //End topos
  }

}


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