Protein.h

/* Copyright 2004
Stanford University

This file is part of the DSR PDB Library.

The DSR PDB Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.

The DSR PDB Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public License
along with the DSR PDB Library; see the file COPYING.LIB.  If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */

#ifndef DSR_PDB_PROTEIN_H
#define DSR_PDB_PROTEIN_H

#include <dsrpdb/Point.h>
#include <dsrpdb/Atom.h>
#include <dsrpdb/Residue.h>
#include <iostream>
#include <vector>
#include <iterator>
#include <cassert>


namespace dsrpdb {

  class Model;
  

  class Protein {
    friend class Model;
  public:


    typedef std::pair<Atom::Index,Atom::Index> Bond;

    Protein(std::istream &in, bool print_errors=false);

    Protein();

    char chain() const;

    void set_chain(char c);
    
    typedef std::vector<Residue>::iterator Residues_iterator;
    Residues_iterator residues_begin() {
      return residues_.begin();
    }
    Residues_iterator residues_end() {
      return residues_.end();
    }
    typedef std::vector<Residue>::const_iterator Const_residues_iterator;
    Const_residues_iterator residues_begin() const {
      return residues_.begin();
    }
    Const_residues_iterator residues_end() const {
      return residues_.end();
    }
    unsigned int number_of_residues() const {
      return residues_.size();
    }

    void new_residue(const Residue &res);
    
    class Atoms_iterator;
    Atoms_iterator atoms_begin();
    Atoms_iterator atoms_end();

    class Const_atoms_iterator;
    Const_atoms_iterator atoms_begin() const;
    Const_atoms_iterator atoms_end() const;
    unsigned int number_of_atoms() const;

    class Bonds_iterator;
    Bonds_iterator bonds_begin() const;
    Bonds_iterator bonds_end() const;
    unsigned int number_of_bonds() const;


    const Residue& residue_containing_atom(Atom::Index atom_index) const;
    Residue& residue_containing_atom(Atom::Index atom_index);

    const Atom& atom(Atom::Index atom_index) const;

    void set_atom(Atom::Index atom_index, const Atom &a);


    std::vector<Residue::Type> sequence() const;

    void write(std::ostream &out) const;

    void write_pdb(std::ostream &out) const;
    void dump(std::ostream &out) const;



    const Residue& residue(Residue::Index i) const;

    bool has_residue(Residue::Index i) const;

    /*const Atom &atom(unsigned int i) const {
      for (Const_residues_iterator it = residues_begin(); it != residues_.end(); ++it){
        
      }
      }*/


    bool has_bonds() const {
      return residues_[0].has_bonds();
    }

    void set_has_bonds(bool tf) {
      for (unsigned int i=0; i< residues_.size(); ++i){
        residues_[i].set_has_bonds(tf);
      }
    }


#if 0


    class Backbone_coordinates_iterator;
    Backbone_coordinates_iterator backbone_coordinates_begin() const;
    Backbone_coordinates_iterator backbone_coordinates_end() const;
#endif
    
  protected:
    void process_line(const char *line);
    
    unsigned int residue_offset_of_atom_index(Atom::Index i) const;

    unsigned int residue_offset(Residue::Index i) const;
    
    std::vector<Residue> residues_;
    std::vector<std::string> header_;
    //static Residue dummy_residue_;
    char chain_;
  };







  class Protein::Atoms_iterator {
    friend class Protein;
    friend class Protein::Const_atoms_iterator;
  public:
    typedef Atom value_type;
    typedef std::forward_iterator_tag iterator_category;
    typedef std::size_t difference_type;
    typedef Residue::Atoms_iterator::reference reference;
    typedef Residue::Atoms_iterator::pointer pointer;

    reference operator*() {
      return *ait_;
    }
    pointer operator->() {
      return ait_.operator->();
    }
    Atoms_iterator operator++() {
      ++ait_;
      if (ait_== aend_) {
        ++rit_;
        if (rit_!= rend_){
          ait_=rit_->atoms_begin();
          aend_= rit_->atoms_end();
        }
      }
      return *this;
    }
    
    bool operator==(const Atoms_iterator& o) const {
      if (rit_ == rend_) return rit_==o.rit_;
      else return rit_== o.rit_ && ait_ == o.ait_;
    }
    bool operator!=(const Atoms_iterator& o) const {
      if (rit_== rend_) return rit_!= o.rit_;
      else return rit_!= o.rit_ || ait_ != o.ait_;
    }

  protected:
    Atoms_iterator(std::vector<Residue>::iterator b, 
                   std::vector<Residue>::iterator e): rit_(b), rend_(e){
      if (b != e) {
        ait_= rit_->atoms_begin();
        aend_= rit_->atoms_end();
      }
    }
    std::vector<Residue>::iterator rit_, rend_;
    Residue::Atoms_iterator ait_, aend_;
  };


  class Protein::Const_atoms_iterator {
    friend class Protein;
  public:
    typedef Atom value_type;
    typedef std::forward_iterator_tag iterator_category;
    typedef std::size_t difference_type;
    typedef Residue::Const_atoms_iterator::reference reference;
    typedef Residue::Const_atoms_iterator::pointer pointer;

    reference operator*() const {
      return *ait_;
    }
    pointer operator->() const {
      return ait_.operator->();
    }
    Const_atoms_iterator operator++() {
      ++ait_;
      if (ait_== aend_) {
        ++rit_;
        if (rit_!= rend_){
          ait_=rit_->atoms_begin();
          aend_= rit_->atoms_end();
        }
      }
      return *this;
    }
    
    bool operator==(const Const_atoms_iterator& o) const {
      if (rit_ == rend_) return rit_==o.rit_;
      else return rit_== o.rit_ && ait_ == o.ait_;
    }
    bool operator!=(const Const_atoms_iterator& o) const {
      if (rit_== rend_) return rit_!= o.rit_;
      else return rit_!= o.rit_ || ait_ != o.ait_;
    }

    Const_atoms_iterator(Protein::Atoms_iterator it){
      rit_= it.rit_;
      rend_= it.rend_;
      ait_= it.ait_;
      aend_= it.aend_;
    }

  protected:
    Const_atoms_iterator(std::vector<Residue>::const_iterator b, 
                         std::vector<Residue>::const_iterator e): rit_(b), rend_(e){
      if (b != e) {
        ait_= rit_->atoms_begin();
        aend_= rit_->atoms_end();
      }
    }
    std::vector<Residue>::const_iterator rit_, rend_;
    Residue::Const_atoms_iterator ait_, aend_;
  };


  class Protein::Bonds_iterator {
    friend class Protein;
  public:
    typedef Bond value_type;
    typedef std::forward_iterator_tag iterator_category;
    typedef std::size_t difference_type;
    typedef const Bond& reference;
    typedef const Bond* pointer;
    
    reference operator*() const {
      return cur_bond();
    }
    pointer operator->() const {
      return &cur_bond();
    }
    Bonds_iterator operator++() {
      if (nc_) {
        nc_=false;
      } else {
        ++ait_;
        if (ait_== aend_) {
          std::vector<Residue>::const_iterator orit= rit_;
          ++rit_;
          if (rit_!= rend_){
            if (static_cast<unsigned int>(orit->index()) +1 == static_cast<unsigned int>(rit_->index())
                && orit->has_atom(Residue::AL_C) && rit_->has_atom(Residue::AL_N)){
              nc_=true;
            } /*else {
              std::cout << "Skipping bond between residues " << orit->index() << " and " << rit_->index() << "\n";
              std::cout << static_cast<unsigned int>(orit->index()) << " " << static_cast<unsigned int>(rit_->index()) << std::endl;
              rit_->dump(std::cout);
              orit->dump(std::cout);
              }*/
            ait_=rit_->bonds_begin();
            aend_= rit_->bonds_end();
          }
        }
      }
      return *this;
    }
     bool operator==(const Bonds_iterator& o) const {
       if (rit_ == rend_) return rit_==o.rit_;
       else return rit_== o.rit_ && ait_ == o.ait_ && nc_== o.nc_;
    }
    bool operator!=(const Bonds_iterator& o) const {
      if (rit_== rend_) return rit_!= o.rit_;
      else return rit_!= o.rit_ || ait_ != o.ait_ || nc_ != o.nc_;
    }
    
    Bonds_iterator(){}
  protected:
    Bonds_iterator(std::vector<Residue>::const_iterator b, 
                   std::vector<Residue>::const_iterator e): rit_(b), rend_(e), nc_(false){
      if (b != e) {
        ait_= rit_->bonds_begin();
        aend_= rit_->bonds_end();
      }
    }

    const Bond &cur_bond() const {
      if (nc_) {
        assert(rit_ != rend_);
        static Bond b;
        b= Bond((rit_-1)->atom(Residue::AL_C).index(),
                rit_->atom(Residue::AL_N).index());
        return b;
      } else return *ait_;
    }

    std::vector<Residue>::const_iterator rit_, rend_;
    Residue::Bonds_iterator ait_, aend_;
    bool nc_;
  };
};


#endif