#ifndef AminoAcidScore_hpp
#define AminoAcidScore_hpp

#include "AminoAcid.h"
#include "Atom.h"
#include "Match.h"


/*
CLASS
  AminoAcidScore
   
  This class serves for a score calculation of sequence of amino acid pairs.
  The score can be based on the PAM250 or BLOSUM62 matrix or any other user defined score matrix.

KEYWORD
    Atom, Match, ScoreMatrix

AUTHORS
    Yuval Inbar (inbaryuv@post.tau.ac.il)
    Maxim Shatsky (maxshats@post.tau.ac.il)
    
    copyright: SAMBA Group , Tel-Aviv Univ. Israel, 1997, 2003.

CHANGES LOG


USAGE
 Before calculating the score of two matched amino acid sequences,the user
 should define a score matrix.This is a 20x20 (more exact MAX_AA x MAX_AA, to allow GAP, OTHER) 
 matrix M when M[i,j] which 
 defines a value for AminoAcid_i against AminoAcid_j (0<=i,j<=19). 
 The amino acids are indexed in the following way (see AminoAcid class for enum RESIDUE_INDEX): 
  ARG,LYS,GLU,GLN,ASP,ASN,PRO,GLY,SER,
  THR,HIS,ALA,ILE,LEU,VAL,PHE,CYS,MET,TYR,TRP. (ARG=0,...,TRP=19)
 It means that M[3,0]=M[GLN,ARG] (to map amino acid to proper index see AminoAcid class) .
 There is a predefined PAM250 score matrix AminoAcidScore::PAM250 ,someone
 can use it as a template for creating another score matrix. 

 EXAMPLE
 
 vector<Atom> vModel;
 vector<Atom> vTarget;
 //...
 //fill atom information into vectors vModel,vTarget  
 //...
 vector<Match*> vMatches; 
 //..
 //find some matches 
 //...

 AminoAcidScore::PAM250 pam250; //define score matrix
 AminoAcidScore seqScore;

 //calculate and output PAM250 score for each found match. 
 for(int i=0;i < vMatches.size();i++){
  float score=seqScore.calcScore(vModel,
                          vTarget,
                          *vMatches[i],
                          pam250);
  score=score/vMatches[i].size(); //normolize the score
  cout<<score<<endl;
 }

 END

 Function calcScore is a template that is defined by 2 classes: 
 <class TMolecule= vector<Atom>,class TMatch>.
For example: TMolecule=vector<Atom>, TMatch=Match.

 Class TMolecule should support the following:
 EXAMPLE
 Molecule<Atom> pdbRec;
 pdbRec[int index].residueType(); //this function returns one-letter abbreviation of amino acid residue.
 END
 
 Class TMatch should support the following:
 EXAMPLE
 Match match;
 match[i].model; //returns index in  pdbRec
 match[i].scene; //returns index in  pdbRec
 END

 If one wishes to create its own score matrix just see how PAM250 class
 is defined.
*/



class AminoAcidScore
{

public:
  ////ScoreMatrix class defines the score for each pair of amino acid residues
  class ScoreMatrix{
  protected:
    float m_fScore[MAX_AA][MAX_AA];
    float m_fMax;
    float m_fMin;

    void updateMinMax(){
      m_fMin=99999;
      m_fMax=-99999;
      
      for(int i=0;i<MAX_AA;++i)
	for(int j=0;j<MAX_AA;++j){
	  if(m_fScore[i][j]> m_fMax)
	    m_fMax=m_fScore[i][j];
	  if(m_fScore[i][j]< m_fMin)
	    m_fMin=m_fScore[i][j];
	}
    }
    
  public:
    ////Returns the score
    float score(int i,int j)const{return m_fScore[i][j];}
    
    ////Returns min score value
    float minValue()const{return m_fMin;};
    
    ////Return max score value
    float maxValue()const{return m_fMax;};
    
    ////sets all diagonal elements to sum(a_ii)/20
    //usefull for computing conservation score.
    void normDiagonal(){
	float sum=0;
	for(int i=0;i<20;++i)
		sum+=m_fScore[i][i];

	sum/=20;
	for(int i=0;i<20;++i)
		m_fScore[i][i]=sum;

	updateMinMax();
    }
    
    ////reverses the value of diagonal elements, i.e. the higest value becomes 
    //the lowest value.
    //usefull for computing conservation score.
    void reverseValueDiagonal(){
	float fmax=m_fScore[0][0],fmin=m_fScore[0][0];
	for(int i=1;i<20;++i){
	  if(m_fScore[i][i]>fmax)
		fmax= m_fScore[i][i];
	  if(m_fScore[i][i]<fmin)
	  	fmin= m_fScore[i][i];
	}
	for(int i=0;i<20;++i)
		m_fScore[i][i]=fmax - m_fScore[i][i] + fmin;

	updateMinMax();
    }
    
    	
  };

  ////Extends ScoreMatrix. The score is defined in the following way:
  //Hydrophobic matched with Hydrophobic receives 1.
  //Hydrophilic matched with Hydrophilic receives 1.
  //Otherwise -1.
  class Hydrophobic:public ScoreMatrix{
  public:
    ////Constructor. Complexity O(400).
    Hydrophobic();
  };

  ////Extends ScoreMatrix. The score is defined as in PAM250 matrix.
  class PAM250:public ScoreMatrix{
  public:
    ////Constructor. Complexity O(400).
    PAM250();
  };
  
  
  
  ////Extends ScoreMatrix. The score is defined as in BLOSUM62 matrix.
  class BLOSUM62:public ScoreMatrix{
  public:
    ////Constructor. Complexity O(400).
    BLOSUM62();
  };

  
  
  ////Returns the score of two matched sequences.Complexity O(match.size()).
  template< class TMolecule,class TMatch>
  float calcScore( const TMolecule& pdbModel,
                   const TMolecule& pdbTarget,
		   const Match& match,
                   const ScoreMatrix& m_scoreMtrx)const{
    float fSum=0;

    for(int i=0;i<match.size();i++){
      RESIDUE_INDEX model=AminoAcid::getResidueTypeIndx(pdbModel[match[i].model].residueType());
      RESIDUE_INDEX target=AminoAcid::getResidueTypeIndx(pdbTarget[match[i].scene].residueType());
      
      fSum+= m_scoreMtrx.score( model , target );
    }
    
    return fSum;
  }

  
};


#endif