// writeDv610.cpp write *.dv files for a protein
// used for potential function CALSP

#include <iostream>
#include <string>
#include <fstream>
#include <iomanip>
#include <set>
#include <map>
#include <list>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <cctype>
#include <cstdio>
#include <math.h>

using namespace std;

#define NUM_ATOM_TYPE 20   // number of atom types
#define NUM_SIMP_GM 5      // number of simplified alphabet for geometric potential
#define MIN_DIS_AA 0.789   // minimum distance of Ca - Ca 3A
#define MIN_DIS_AB 0.658   // minimum distance of Ca - Cb 2.5A
#define MIN_DIS_BB 0.526   // minimum distance of Cb - Cb 2A
#define G_TYPE 5          // Gly 
#define SIM_G_TYPE 3       // simplified glycine type
int Num_Simp_Gm=5;
int Num_Simp_Ct=20;

typedef set<long int> ISET;  // define a set  
typedef map <string, int > CIMAP;
typedef map <string, string > SSMAP;
typedef map <int, int> IIMAP;
typedef multimap < float, string > FSMMAP;
typedef map <int, float> IFMAP;
typedef multimap <int, string> ISMMAP;
typedef map <string, float> SFMAP;

CIMAP AIMap;              // character to integer map
CIMAP::iterator aiItr;    // pointer of the CIMAP
SSMAP AAMap;              // string to string map
SSMAP::iterator aaItr;    // pointer for SSMAP
IIMAP PtMap;              // map a feature to a integer
IIMAP::iterator iiItr;
CIMAP SimGmMap;           // simplified geometric aa to int map
CIMAP::iterator sgItr;
CIMAP SimCtMap;           // simplified contact aa to int map
CIMAP::iterator scItr;
CIMAP ProtLenMap;         // string and int map for protein length
CIMAP::iterator plItr;
SSMAP ProtSeqMap;         // string and string map for protein sequence
SSMAP::iterator psItr;
SSMAP ProtSSMap;          // ssmap for protein and state sequence
SSMAP::iterator pssItr;
SFMAP ProtEnMap; 
SFMAP::iterator sfItr;
//IFMAP DVMap;		  // map feature vector to float	
IFMAP::iterator ifItr;
IFMAP NatDvMap;       // native descriptor vector
int Num_State=4;
int Num_Disbin=1;
int Simp_Ct_Type=0;    // if 0 type is Li&Liang if 1 type is Murphy&Levy
ifstream fileList;
ofstream outcst; 
float cont_Dis[NUM_ATOM_TYPE][NUM_ATOM_TYPE];
int Num_DV=610;          // number of descriptors
int isD=0;    // if =1 write DV for decoy structure
string natProtStr;
inline int toSG(int state, int sgm) { return state*Num_Simp_Gm+sgm;}
inline int gType() { return Num_Simp_Ct==10? 3 : 5;}
void error(char* err) {{cout<<err<<endl; exit(0);}}
void Init(int, char**);
void Usage();
void GetDistance();
void editMap();

int main(int argc, char* argv[])
{
  int i,j,natLength;   // length of template and native protein
  char temp[5000],protName[200];  // temporary storage and protein name
  string natName,natStr;  // natStr is the native protein, which is threaded on the template
  int pos1,pos2,index,disbin,cont_type,dvIndex;
  char aaStr1[20],aaStr2[20],atom1[20],atom2[20];
  int natSS[1500],natSGS[1500];  // integer of template and native simplified geometric sequence
  int natSCS[1500];  // integer of template and native simplified contact sequence
  ifstream infile;
  Init(argc,argv);
cout<<"Init done"<<endl;
  while(!fileList.eof()){
    fileList.getline(temp,200);
    if(strlen(temp)<5)
      break;
    if(isD==1){
//	cout<<"isD is 1"<<endl; 
     sscanf(temp,"%s",protName);
    }
    else sscanf(temp,"%s %d",protName,&natLength);
    natName=protName;
    if(isD==1){
      natName=natName.substr(0,natName.length()-4);
      strcpy(protName,natName.c_str());
    }
    NatDvMap.clear();
    cout<<natName<<endl;
    // get the protein's simplified geometric and contact aa sequence,
    if(isD==1)
      natStr=natProtStr;
    else if((psItr=ProtSeqMap.find((string)protName))!=ProtSeqMap.end())
      natStr=psItr->second;
    else error("ProtSeqMap error!");
    natLength=natStr.length();
    for(i=1;i<=natStr.length();++i){
      if((sgItr=SimGmMap.find(natStr.substr(i-1,1)))!=SimGmMap.end())
	natSGS[i]=sgItr->second;
      else error("SimGmMap error!");
      if((scItr=SimCtMap.find(natStr.substr(i-1,1)))!=SimCtMap.end())
	natSCS[i]=scItr->second;
      else error("SimCtMap error!");
      //      cout<<tempSGS[i]<<" "; 
    }
    // get the protein's state sequence
    if((pssItr=ProtSSMap.find((string)protName))!=ProtSSMap.end())
      natStr=pssItr->second;
    else { cout<<ProtSSMap.size()<<" "<<ProtSSMap.begin()->first<<" "<<natStr<<endl; error("ProtSSMap error!");}
    for(i=0;i<natStr.length();++i){
      natSS[i+3]=atoi(natStr.substr(i,1).c_str());
      //      cout<<tempSS[i+3]<<" ";
    }
    //    cout<<endl;
    // set "seudo state" for the first two residues for convenience
    natSS[1]=natSS[3];
    natSS[2]=natSS[3];
    //    natSS[natStr.length()]=natSS[natStr.length()-1];
    strcpy(temp,protName);
    strcat(temp,".610.dv");
    outcst.open(temp);
    strcpy(temp,protName);
    strcat(temp,".abc");
    infile.open(temp);
    if(!infile.is_open()){
      strcat(protName,".sce.abc");
      infile.open(protName);
      if(!infile.is_open()) {
	strcpy(temp,"/home/jinfeng/SharedData/pc30_res2.0_SC_E_Cont/");
	strcat(temp,protName);
	infile.open(temp);
	if(!infile.is_open()) {
	  strcpy(temp,"/cluster/home/jinfeng/SharedData/pc30_res2.0_SC_E_Cont/");
	  strcat(temp,protName);
	  infile.open(temp);
	  if(!infile.is_open()) {
	    //error("cannot open abc file!");
	    cout<<"cannot open abc file "<<temp<<endl;
	    exit(0);
	  }
	}
      }
    }
    while(!infile.eof()){
      infile.getline(temp,2000);
      sscanf(temp,"%d %s %s %d %s %s %*f %d",&pos1,aaStr1,atom1,&pos2,aaStr2,atom2,&disbin);
      if(strlen(temp)<4) {
	break;
      }
      if(Num_Disbin==3)
	disbin=disbin-1;
      else disbin=0;
      if(strcmp(atom1,"CA")==0&&strcmp(atom2,"CA")==0){
	cont_type=0;
	if((iiItr=PtMap.find(100000+gType()*1000+gType()*10+disbin))!=PtMap.end())
	  dvIndex=iiItr->second;
	else error("PtMap error!");
      }
      else if(strcmp(atom1,"CA")==0&&strcmp(atom2,"SC")==0){
	cont_type=1;
	if((iiItr=PtMap.find(100000+gType()*1000+natSCS[pos2]*10+disbin))!=PtMap.end())
	  dvIndex=iiItr->second;
	else error("PtMap error!");
      }
      else if(strcmp(atom1,"SC")==0&&strcmp(atom2,"CA")==0){
	cont_type=2;
	if((iiItr=PtMap.find(100000+natSCS[pos1]*1000+gType()*10+disbin))!=PtMap.end())
	  dvIndex=iiItr->second;
	else error("PtMap error!");
      }
      else if(strcmp(atom1,"SC")==0&&strcmp(atom2,"SC")==0){
	cont_type=3;
	if((iiItr=PtMap.find(100000+natSCS[pos1]*1000+natSCS[pos2]*10+disbin))!=PtMap.end())
	  dvIndex=iiItr->second;
	else error("PtMap error!");
      }
      else { cout<<temp<<" "<<atom1<<" "<<atom2<<" atom type wrong!"<<endl; exit(0);}
      // calculate contact part of DVMap
      //cout<<dvIndex<<" ";
      if((ifItr=NatDvMap.find(dvIndex))==NatDvMap.end())
	NatDvMap.insert(IFMAP::value_type(dvIndex,1));
      else
	ifItr->second++;
    }
    
    // calculate geometric part of DVMap
    for(i=3;i<(natLength-1);++i){
      if((iiItr=PtMap.find(toSG(natSS[i],natSGS[i])*100+toSG(natSS[i+1],natSGS[i+1])))!=PtMap.end())
	dvIndex=iiItr->second;
      else error("PtMap error!");
      if((ifItr=NatDvMap.find(dvIndex))==NatDvMap.end())
	NatDvMap.insert(IFMAP::value_type(dvIndex,1));
      else
	ifItr->second++;
    }
    outcst<<natName<<":"<<natLength<<" ";
    for(ifItr=NatDvMap.begin();ifItr!=NatDvMap.end();++ifItr)
      outcst<<ifItr->first<<":"<<ifItr->second<<" ";
    outcst<<endl;
    infile.close();
    outcst.close();
  }  
}


void Init(int argc, char** argv)
{
  int i=1,j,len,dvIndex,dvKey;
  char potential[50],protName[50],temp[2000],seq[1500],outfile[200],natProtName[50];
  float energy;
  ifstream infile;
  string str;
  
  if(argc==1) {
    Usage();
  }
  strcpy(potential,"bmt_best.cgp");
  strcpy(outfile,"tmp.dat");

  while(i<argc) {
    if(strcmp(argv[i],"-fl")==0) { // protein/sequence name
      if(argc==i)
	Usage();
      strcpy(temp,argv[i+1]);
      fileList.open(temp);
      i=i+2;
    }
    else if(strcmp(argv[i],"-ns")==0) { // number of states
      if(argc==i)
	Usage();
      Num_State=atoi(argv[i+1]);
      i=i+2;
    }
    else if(strcmp(argv[i],"-nd")==0) { // number of descriptors
      if(argc==i)
	Usage();
      Num_DV=atoi(argv[i+1]);
      i=i+2;
    }
    else if(strcmp(argv[i],"-d")==0) { // number of descriptors
      isD=1;
      i=i+1;
    }
    else if(strcmp(argv[i],"-pn")==0) { // protein/sequence name
      if(argc==i)
	Usage();
      strcpy(natProtName,argv[i+1]);
      i=i+2;
    }
    else if(strcmp(argv[i],"-nsc")==0) { // number of simplified contact type
      if(argc==i)
	Usage();
      Num_Simp_Ct=atoi(argv[i+1]);
      i=i+2;
    }
    else if(strcmp(argv[i],"-ndb")==0) { // number of distance bin
      if(argc==i)
	Usage();
      Num_Disbin=atoi(argv[i+1]);
      i=i+2;
    }
    else if(strcmp(argv[i],"-nsg")==0) { // number of simplified geometric type
      if(argc==i)
	Usage();
      Num_Simp_Gm=atoi(argv[i+1]);
      i=i+2;
    }
    else if(strcmp(argv[i],"-sct")==0) { // type of simplified contacts
      if(argc==i)
	Usage();
      Simp_Ct_Type=atoi(argv[i+1]);
      i=i+2;
    }
  }

  // get the protein names and their sequence
  if(isD==1){
    strcpy(temp,natProtName);
    strcat(temp,".seq");
    infile.open(temp);
    if(!infile.is_open()) {
      error(".seq file error!");
    }
    infile.getline(temp,2000);
    sscanf(temp,"%s %s",protName,seq);
    natProtStr=(string)seq;
    ProtSeqMap.insert(SSMAP::value_type((string)protName,(string)seq));
  }
  else {
    infile.open("/home/jinfeng/SharedData/sequence.db");
    if(!infile.is_open()) {
      infile.open("/cluster/home/jinfeng/SharedData/sequence.db");
    }
    while(!infile.eof()){
      infile.getline(temp,2000);
      sscanf(temp,"%s %s",protName,seq);
      ProtSeqMap.insert(SSMAP::value_type((string)protName,(string)seq));
    }
  }
  infile.close();

  // get the protein name and their state sequence
  if(isD==1){
    strcpy(temp,natProtName);
    strcat(temp,".dec.stsq");
    infile.open(temp);
    if(!infile.is_open())
      error("dec.stsq error!");
    while(!infile.eof()){
      infile.getline(temp,2000);
      sscanf(temp,"%s %s",protName,seq);
      ProtSSMap.insert(SSMAP::value_type((string)protName,(string)seq));
    }
    infile.close();
  }
  else {
    infile.open("/home/jinfeng/SharedData/stateSeq.dat");
    if(!infile.is_open()) {
      infile.open("/cluster/home/jinfeng/SharedData/stateSeq.dat");
      if(!infile.is_open())
	error("stateSeq file error");
    }
    while(!infile.eof()){
      infile.getline(temp,2000);
      sscanf(temp,"%s %s",protName,seq);
      ProtSSMap.insert(SSMAP::value_type((string)protName,(string)seq));
    }
    infile.close();
  }
  //  GetDistance();
  editMap();
}

void GetDistance()
{
  int i,j;
  char tmp[200];
  ifstream infile;
  infile.open("ContDis.txt");
  if(!infile.is_open()) {
    infile.open("/home/jinfeng/SharedData/ContDis.txt");
    if(!infile.is_open()) {
      infile.open("/cluster/home/jinfeng/SharedData/ContDis.txt");
      if(!infile.is_open()) {
	cout<<"Erro opening file ContDis.txt correctly"<<endl;
	exit(0);
      }
    }
  }
  for(i=0;i<NUM_ATOM_TYPE;++i) {
    for(j=0;j<NUM_ATOM_TYPE;++j) {
      infile>>cont_Dis[i][j];
      //            cout<<cont_Dis[i][j]<<" ";
    }
    //  cout<<endl;
  }
infile.close();
}

void editMap()
{
  AAMap.insert(SSMAP::value_type("ALA", "A"));
  AAMap.insert(SSMAP::value_type("CYS", "C"));
  AAMap.insert(SSMAP::value_type("ASP", "D"));
  AAMap.insert(SSMAP::value_type("GLU", "E"));
  AAMap.insert(SSMAP::value_type("PHE", "F"));
  AAMap.insert(SSMAP::value_type("GLY", "G"));
  AAMap.insert(SSMAP::value_type("HIS", "H"));
  AAMap.insert(SSMAP::value_type("ILE", "I"));
  AAMap.insert(SSMAP::value_type("LYS", "K"));
  AAMap.insert(SSMAP::value_type("LEU", "L"));
  AAMap.insert(SSMAP::value_type("MET", "M"));
  AAMap.insert(SSMAP::value_type("ASN", "N"));
  AAMap.insert(SSMAP::value_type("PRO", "P"));
  AAMap.insert(SSMAP::value_type("GLN", "Q"));
  AAMap.insert(SSMAP::value_type("ARG", "R"));
  AAMap.insert(SSMAP::value_type("SER", "S"));
  AAMap.insert(SSMAP::value_type("THR", "T"));
  AAMap.insert(SSMAP::value_type("VAL", "V"));
  AAMap.insert(SSMAP::value_type("TRP", "W"));
  AAMap.insert(SSMAP::value_type("TYR", "Y"));
  AAMap.insert(SSMAP::value_type("A", "ALA"));
  AAMap.insert(SSMAP::value_type("C", "CYS"));
  AAMap.insert(SSMAP::value_type("D", "ASP"));
  AAMap.insert(SSMAP::value_type("E", "GLU"));
  AAMap.insert(SSMAP::value_type("F", "PHE"));
  AAMap.insert(SSMAP::value_type("G", "GLY"));
  AAMap.insert(SSMAP::value_type("H", "HIS"));
  AAMap.insert(SSMAP::value_type("I", "ILE"));
  AAMap.insert(SSMAP::value_type("K", "LYS"));
  AAMap.insert(SSMAP::value_type("L", "LEU"));
  AAMap.insert(SSMAP::value_type("M", "MET"));
  AAMap.insert(SSMAP::value_type("N", "ASN"));
  AAMap.insert(SSMAP::value_type("P", "PRO"));
  AAMap.insert(SSMAP::value_type("Q", "GLN"));
  AAMap.insert(SSMAP::value_type("R", "ARG"));
  AAMap.insert(SSMAP::value_type("S", "SER"));
  AAMap.insert(SSMAP::value_type("T", "THR"));
  AAMap.insert(SSMAP::value_type("V", "VAL"));
  AAMap.insert(SSMAP::value_type("W", "TRP"));
  AAMap.insert(SSMAP::value_type("Y", "TYR"));

  // initialize AIMap
  AIMap.insert(CIMAP::value_type("A",0));
  AIMap.insert(CIMAP::value_type("C",1));
  AIMap.insert(CIMAP::value_type("D",2));
  AIMap.insert(CIMAP::value_type("E",3));
  AIMap.insert(CIMAP::value_type("F",4));
  AIMap.insert(CIMAP::value_type("G",5));
  AIMap.insert(CIMAP::value_type("H",6));
  AIMap.insert(CIMAP::value_type("I",7));
  AIMap.insert(CIMAP::value_type("K",8));
  AIMap.insert(CIMAP::value_type("L",9));
  AIMap.insert(CIMAP::value_type("M",10));
  AIMap.insert(CIMAP::value_type("N",11));
  AIMap.insert(CIMAP::value_type("P",12));
  AIMap.insert(CIMAP::value_type("Q",13));
  AIMap.insert(CIMAP::value_type("R",14));
  AIMap.insert(CIMAP::value_type("S",15));
  AIMap.insert(CIMAP::value_type("T",16));
  AIMap.insert(CIMAP::value_type("V",17));
  AIMap.insert(CIMAP::value_type("W",18));
  AIMap.insert(CIMAP::value_type("Y",19));
  
  // initialize SimGmMap (following NEW aa_clus result)
  SimGmMap.insert(CIMAP::value_type("A",1));
  SimGmMap.insert(CIMAP::value_type("C",0));
  SimGmMap.insert(CIMAP::value_type("D",4));
  SimGmMap.insert(CIMAP::value_type("E",1));
  SimGmMap.insert(CIMAP::value_type("F",0));
  SimGmMap.insert(CIMAP::value_type("G",2));
  //  SimGmMap.insert(CIMAP::value_type("H",2));
  SimGmMap.insert(CIMAP::value_type("H",1));
  SimGmMap.insert(CIMAP::value_type("I",0));
  SimGmMap.insert(CIMAP::value_type("K",1));
  SimGmMap.insert(CIMAP::value_type("L",0));
  SimGmMap.insert(CIMAP::value_type("M",0));
  SimGmMap.insert(CIMAP::value_type("N",4));
  SimGmMap.insert(CIMAP::value_type("P",3));
  SimGmMap.insert(CIMAP::value_type("Q",1));
  SimGmMap.insert(CIMAP::value_type("R",1));
  //  SimGmMap.insert(CIMAP::value_type("S",2));
  //  SimGmMap.insert(CIMAP::value_type("T",2));
  SimGmMap.insert(CIMAP::value_type("S",1));
  SimGmMap.insert(CIMAP::value_type("T",1));
  SimGmMap.insert(CIMAP::value_type("V",0));
  SimGmMap.insert(CIMAP::value_type("W",0));
  SimGmMap.insert(CIMAP::value_type("Y",0));
  if(Num_Simp_Ct==10&&Simp_Ct_Type==1){ 
    // initialize SimCtMap (following results by Levy & Murphy)
    SimCtMap.insert(CIMAP::value_type("A",2));
    SimCtMap.insert(CIMAP::value_type("C",1));
    SimCtMap.insert(CIMAP::value_type("D",7));
    SimCtMap.insert(CIMAP::value_type("E",7));
    SimCtMap.insert(CIMAP::value_type("F",6));
    SimCtMap.insert(CIMAP::value_type("G",3));
    SimCtMap.insert(CIMAP::value_type("H",9));
    SimCtMap.insert(CIMAP::value_type("I",0));
    SimCtMap.insert(CIMAP::value_type("K",8));
    SimCtMap.insert(CIMAP::value_type("L",0));
    SimCtMap.insert(CIMAP::value_type("M",0));
    SimCtMap.insert(CIMAP::value_type("N",7));
    SimCtMap.insert(CIMAP::value_type("P",5));
    SimCtMap.insert(CIMAP::value_type("Q",7));
    SimCtMap.insert(CIMAP::value_type("R",8));
    SimCtMap.insert(CIMAP::value_type("S",4));
    SimCtMap.insert(CIMAP::value_type("T",4));
    SimCtMap.insert(CIMAP::value_type("V",0));
    SimCtMap.insert(CIMAP::value_type("W",6));
    SimCtMap.insert(CIMAP::value_type("Y",6));
  }
  else if(Num_Simp_Ct==10&&Simp_Ct_Type==0){
  // initialize SimCtMap, exact number is still 10 (following Li and Liang)
    SimCtMap.insert(CIMAP::value_type("A",8));
    SimCtMap.insert(CIMAP::value_type("C",0));
    SimCtMap.insert(CIMAP::value_type("D",1));
    SimCtMap.insert(CIMAP::value_type("E",1));
    SimCtMap.insert(CIMAP::value_type("F",6));
    SimCtMap.insert(CIMAP::value_type("G",9));
    SimCtMap.insert(CIMAP::value_type("H",4));
    SimCtMap.insert(CIMAP::value_type("I",5));
    SimCtMap.insert(CIMAP::value_type("K",2));
    SimCtMap.insert(CIMAP::value_type("L",5));
    SimCtMap.insert(CIMAP::value_type("M",6));
    SimCtMap.insert(CIMAP::value_type("N",3));
    SimCtMap.insert(CIMAP::value_type("P",4));
    SimCtMap.insert(CIMAP::value_type("Q",3));
    SimCtMap.insert(CIMAP::value_type("R",2));
    SimCtMap.insert(CIMAP::value_type("S",3));
    SimCtMap.insert(CIMAP::value_type("T",3));
    SimCtMap.insert(CIMAP::value_type("V",5));
    SimCtMap.insert(CIMAP::value_type("W",7));
    SimCtMap.insert(CIMAP::value_type("Y",7));	  
  }  
  else if(Num_Simp_Ct==20){
    // initialize SimCtMap with each aa as a type
    SimCtMap.insert(CIMAP::value_type("A",0));
    SimCtMap.insert(CIMAP::value_type("C",1));
    SimCtMap.insert(CIMAP::value_type("D",2));
    SimCtMap.insert(CIMAP::value_type("E",3));
    SimCtMap.insert(CIMAP::value_type("F",4));
    SimCtMap.insert(CIMAP::value_type("G",5));
    SimCtMap.insert(CIMAP::value_type("H",6));
    SimCtMap.insert(CIMAP::value_type("I",7));
    SimCtMap.insert(CIMAP::value_type("K",8));
    SimCtMap.insert(CIMAP::value_type("L",9));
    SimCtMap.insert(CIMAP::value_type("M",10));
    SimCtMap.insert(CIMAP::value_type("N",11));
    SimCtMap.insert(CIMAP::value_type("P",12));
    SimCtMap.insert(CIMAP::value_type("Q",13));
    SimCtMap.insert(CIMAP::value_type("R",14));
    SimCtMap.insert(CIMAP::value_type("S",15));
    SimCtMap.insert(CIMAP::value_type("T",16));
    SimCtMap.insert(CIMAP::value_type("V",17));
    SimCtMap.insert(CIMAP::value_type("W",18));
    SimCtMap.insert(CIMAP::value_type("Y",19));
  }
  int i,j,k,l,m;
  k=1;
  // geometric
  for(i=0;i<(Num_State*Num_Simp_Gm);++i)
    for(j=0;j<(Num_State*Num_Simp_Gm);++j) { 
      PtMap.insert(IIMAP::value_type(i*100+j,k));
      k++;	
    }	

  for(i=0;i<Num_Simp_Ct;++i)
    for(j=i;j<Num_Simp_Ct;++j)
      for(l=0;l<Num_Disbin;++l){
	PtMap.insert(IIMAP::value_type(100000+i*1000+j*10+l,k));
	PtMap.insert(IIMAP::value_type(100000+j*1000+i*10+l,k));	
	k++;	
      }
  // insert the three over-packing parameters
  PtMap.insert(IIMAP::value_type(-1,k++));   // for CA-CA
  PtMap.insert(IIMAP::value_type(-2,k++));   // for CA-SC
  PtMap.insert(IIMAP::value_type(-3,k++));   // for SC-SC
}
 
void Usage()
{
  cout<<"Usage: "<<endl;
  // to be added
}