#ifndef _NUCLEOTIDE_H
#define _NUCLEOTIDE_H

#include "NucleicAcidAtom.h"

#include <vector>
#include <iostream>
#include <fstream>
#define MIN_BACKBONE 12
#define MIN_BASE_SIZE 8

using std::vector;

/*
CLASS
  Nucleotide
  A container that stores all NucleicAcidAtoms of a nucleotide          

KEYWORDS
  Nucleic Acids, RNA, DNA

AUTHORS
  Oranit Dror (mailto: oranit@post.tau.ac.il)

  Copyright: SAMBA group, Tel-Aviv Univ. Israel, 2004.

CHANGES LOG
<UL>
<LI>
18.04.07 Alex, modified parsePDB to check the validity of the nucleotide records
added getBaseCenter

</LI>
</UL>

  OVERVIEW TEXT

USAGE
  
  EXAMPLE    
      Molecule<NucleicAcidAtom> rna;
      ifstream rnaFile("rna.pdb");
      rna.readPDBfile(rnaFile);

      vector<Nucleotide> rna_nc = Nucleotide::parsePDB(rna);
      
  END
END
*/
class Nucleotide : 
public vector<NucleicAcidAtom> {

public:
  // GROUP: Iterators

  ////
  class BaseIterator {
  public:
    ////
    BaseIterator(const Nucleotide* nucleotide_) : 
      nucleotide(nucleotide_), nucleotideIter(nucleotide_->begin()) {
      if (!nucleotideIter->isBaseAtom()) {
	operator++();
      }
    }

    ////
    // Proceede to the next base atom of the nucleotide
    void operator++() {
      const_iterator iter = nucleotideIter;
      ++iter;
      for (; iter != nucleotide->end() ; iter++) {	
	if (iter->isBaseAtom()) {
	  nucleotideIter = iter;
	  return;
	}
      }

      // No more base atoms. 
      // Thus, the iterator is past-the-end. 
      nucleotideIter = nucleotide->end();
    }


    ////
    // Returns the current base atom (Dereference)
    const NucleicAcidAtom& operator*()  { 
      return *nucleotideIter;      
    }

    ////
    // Returns true if the iterator points beyond the last base atom
    // of the nucleotide 
    bool isDone() {
      return (nucleotideIter == nucleotide->end());      
    }
  private:
    ////
    const Nucleotide* nucleotide;
    
    ////
    const_iterator nucleotideIter;
  };
  // GROUP: Query
  
  ////
  // Returns the phosphate atom of the nucleotide
  const NucleicAcidAtom& getP() const;

  ////
  // Returns the O3 atom of the nucleotide
  const NucleicAcidAtom& getO3() const;


  ////
  // Returns the chain identifier of the nucleotide
  char chainId() const;


  ////
  // Returns the residiue identifier of the nucleotide
  unsigned short residueIndex() const;

  ////
  // Returns the 3D coordinates of the phosphate atom
  Vector3 position() const{
    return getP();
  }  

  ////
  // Returns a normal to the nucleotide base, which is a planar
  // molecule. 
  Vector3 getBaseNormal() const;
  float getBaseMeanDistance(Nucleotide *nc) const;
  float getBaseBaseDistance(Nucleotide *nc) const;
  
  ////
  // Returns the base centroid 
  Vector3 getBaseCenter(void) const; 

  // GROUP: Operators

  ////
  // Returns the 3D coordinates of the nucleotide's phosphate atom 
  // To Do: To check if we need this method. currently we need it so
  // we will be able to use the Match::calculateRMSD method.
  operator Vector3()const { 
    return getP();
  }
  
  ////
  // Applies a rigid transformation on all the nucleotide's atoms.
  Nucleotide& operator*=(const RigidTrans3& transformation);

  //  static vector<Nucleotide> Nucleotide::parsePDB(const vector<NucleicAcidAtom>  i_contAtom);
  // GROUP: Output
  
  ////
  // Prints a PDB record for each atom of the nucleotide
  friend ostream& operator<<(ostream& output, const Nucleotide& nucleotide); 

  ////
  //  Return the difference between the base main directions
  static float getBaseDirectionDiff(Nucleotide *nc1, Nucleotide *nc2);

  ////
  // Parses the input molecule of NucleicAcidAtom to an array of Nucleotides
  static vector<Nucleotide> parsePDB(const vector<NucleicAcidAtom>  i_contAtom){

    vector<Nucleotide> o_vAA;

    //nucleotid index in order to differentiate between
    //different residues
    int iCurrRes;
    string sCurrRes;

    if (i_contAtom.size()==0)
        return o_vAA;

    iCurrRes=i_contAtom[0].residueIndex();     
    sCurrRes=i_contAtom[0].residueSequenceID();

    Nucleotide vRes;
    int backbone_count = 0;
    int base_count = 0;
    bool phosphate_found = false;

    for (int i=0;i<(int)i_contAtom.size();i++) {

        if (i_contAtom[i].residueIndex()!=iCurrRes || i_contAtom[i].residueSequenceID()!=sCurrRes) {
	  if ((backbone_count >= MIN_BACKBONE) && (base_count >= MIN_BASE_SIZE) && (phosphate_found)){
            o_vAA.push_back(vRes);
	  }
            vRes.clear();     
	    backbone_count = 0;
	    base_count = 0;
	    phosphate_found = false;

            iCurrRes=i_contAtom[i].residueIndex();     
            sCurrRes=i_contAtom[i].residueSequenceID();
        }
	if (i_contAtom[i].isBackbone()){
	  backbone_count++;
	}
	else if (i_contAtom[i].isBaseAtom()){
	  base_count++;
	}
	if (i_contAtom[i].isP()){
	  phosphate_found = true;
	}
        vRes.push_back(i_contAtom[i]);

    }    

    //last residue
    if ((vRes.size() > 0) && (backbone_count >= MIN_BACKBONE) && (base_count >= MIN_BASE_SIZE)  && (phosphate_found)) {
        o_vAA.push_back(vRes);        
    }


    return o_vAA; 
}

};

/*
CLASS
     NucleotideData
     A container that stores additional information for the class Nucleotid

KEYWORDS
     Nucleic Acids, RNA, DNA

AUTHORS
  Alex Shulman (mailto: shulmana@post.tau.ac.il)

  Copyright: SAMBA group, Tel-Aviv Univ. Israel, 2007.

  OVERVIEW TEXT

*/

class NucleotideData : public Nucleotide{
public:
  Vector3 p1;    
  Vector3 p2;
  Vector3 normal;
  
  ////
  // Basic constructor. Gets the vector of the main ring direction 
 
  NucleotideData(Nucleotide *nc) {
  BaseIterator baseIter1(nc);
  BaseIterator baseIter2(nc);
  
  Vector3 s1,s2;
  float max_dist = 0;
  while (!baseIter1.isDone()){
    s1=(Vector3)(*baseIter1);
    while (!baseIter2.isDone()){
      s2=(Vector3)(*baseIter2);
      float dist = s1.dist2(s2);
      if (dist > max_dist){
	p1 = s1;
	p2 = s2;
	max_dist = dist;
      }
      ++baseIter2;
    }
    ++baseIter1;
    break;
  }
  normal = (Vector3)(p2-p1);
  normal = normal.getUnitVector();
  }
};
  
#endif //_NUCLEOTIDE_H