atom_classifier.h

Go to the documentation of this file.

// Copyright (c) 2009-2010  INRIA Sophia-Antipolis (France).
// All rights reserved.
//
//This file is part of ESBTL.
//
//ESBTL 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 3 of the License, or
//(at your option) any later version.
//
//ESBTL 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 General Public License for more details.
//
//You should have received a copy of the GNU General Public License
//along with ESBTL.  If not, see <http://www.gnu.org/licenses/>.
//
//
//Additional permission under GNU GPL version 3 section 7
//
//If you modify this Library, or any covered work, by linking or
//combining it with CGAL (or a modified version of that library), the
//licensors of this Library grant you additional permission to convey
//the resulting work. Corresponding Source for a non-source form of
//such a combination shall include the source code for the parts of CGAL
//used as well as that of the covered work. 
//
//
//
// Author(s)     :  Sébastien Loriot


#ifndef ESBTL_ATOM_CLASSIFIER_H
#define ESBTL_ATOM_CLASSIFIER_H

#include <boost/unordered_map.hpp>
#include<boost/lexical_cast.hpp>
#include <boost/algorithm/string.hpp>
#include <iostream>
#include <sstream>
#include <fstream>


namespace ESBTL{

template <class NT,class Atom>
class Radius_of_atom{
private:
  typedef Radius_of_atom<NT,Atom>  Self;
  NT radius_;
  unsigned index_;

public:
  typedef Atom                  Query_type;
  typedef std::string           Key_type;
  typedef NT                    Value_type;
  

  //TODO check why index is useful
  Radius_of_atom(const NT& radius,const unsigned& index):radius_(radius),index_(index){}
  NT value() const {return radius_;}

  //TODO check why index is useful  
  Radius_of_atom(std::stringstream& ss,unsigned index):index_(index){
    ss >> radius_; 
  }
  
  template <class Dictionary,class Vector_properties>
  static unsigned default_loader(Dictionary& dict,Vector_properties& vect)
  {
    #include <ESBTL/properties/Tsai_jmb_99_radii.h>
  }
  
  static std::string make_key(const Atom& atom)
  {
    return atom.residue_name()+atom.atom_name();
  }

  template <class Dictionary>
  static 
  unsigned add_classification(std::stringstream& ss,Dictionary& dict){
    std::string resname,atmname;
    ss >> resname;
    ss >> atmname;
    unsigned prop;
    ss >> prop;
    dict[resname+atmname]=prop;
    return prop;
  }
  
  static void handle_extra(std::stringstream&){}
    
  static
  int& index_of_default(){
    static int index_of_default=-1;
    return index_of_default;
  }
  
};

template <class NT,class Atom>
NT get_radius(const Radius_of_atom<NT,Atom>& property){
  return property.value();
}


template <class NT,class Atom>
class Name_and_radius_of_atom{
private:
  typedef Name_and_radius_of_atom<NT,Atom>  Self;
  // this value must be updated according to file default_atom_properties.h
  static 
  unsigned& water_index(){
    static unsigned index=9;
    return index;
  }

  static
  std::string water_name(){
    static std::string water_name="Owat";
    return water_name;
  }

public:
  
  typedef Atom                  Query_type;
  typedef std::string           Key_type;
  

  Name_and_radius_of_atom(const std::string& name_,const NT& radius_,const unsigned& index_):name(name_),index(index_),radius(radius_)
  {}
  
  std::string name;
  unsigned index;
  NT radius;
  
  bool is_water() const {
    assert(index!=water_index() || name==water_name());
    return index==water_index();
  }
    
  template <class Dictionary,class Vector_properties>
  static unsigned default_loader(Dictionary& dict,Vector_properties& vect)
  {
    #include <ESBTL/properties/default_atom_properties.h>
  }
  
  static std::string make_key(const Atom& atom)
  {
    return atom.residue_name()+atom.atom_name();
  }
  
  static
  int& index_of_default(){
    static int index_of_default=-1;
    return index_of_default;
  }
  
};

class Name_of_pair{
private:
  typedef Name_of_pair Self;


public:
  

  static
  int& index_of_default(){
    static int index_of_default=-1;
    return index_of_default;
  }
  
  
  //~ typedef std::pair<unsigned,unsigned> Key_type;
  typedef std::string Key_type;
  typedef std::pair<unsigned,unsigned> Query_type;

  Name_of_pair(const std::string& name_,const unsigned& index_):name(name_),index(index_){}
  std::string name;
  unsigned index;
    
  template <class Dictionary,class Vector_properties>
  static unsigned default_loader(Dictionary& dict,Vector_properties& vect)
  {
    #include <ESBTL/properties/default_pair_properties.h>
  }
  
  static Key_type make_key(const Query_type& indices)
  {
    unsigned i0=indices.first;
    unsigned i1=indices.second;
    if (i0<i1)
      return boost::lexical_cast<std::string>(i0)+
             boost::lexical_cast<std::string>(i1);
      //~ return indices;
    //~ return std::make_pair(i1,i0);
    return boost::lexical_cast<std::string>(i1)+
             boost::lexical_cast<std::string>(i0);
  }
};



template<class Properties_>
struct Generic_classifier{
  typedef Properties_ Properties;
  typedef typename Properties_::Key_type                      Key_type;
  typedef typename Properties_::Query_type                    Query_type;
  typedef typename std::vector<Properties_>::iterator         Properties_iterator;
  typedef typename std::vector<Properties_>::const_iterator   Properties_const_iterator;
private:
  typedef boost::unordered_map<Key_type, unsigned>      Internal_map_type;
  typedef std::vector<Properties_>                      Internal_vector_type;
  
  Internal_map_type hmap_;
  Internal_vector_type properties_;
  unsigned max_index_;
  
public:
  
  Generic_classifier(){
    max_index_=Properties_::default_loader(hmap_,properties_);
    assert(max_index_==properties_.size());
  }
  
  Generic_classifier(std::string filename,int){
    /* may need a map from name to property position in the vector */
    /* This map can be filled by a static function of Properties_  */
  }
  
  Generic_classifier(const std::string& filename){
    enum Read_state {PROPERTIES, CLASSIFICATION, END, EXTRA,DEFAULT};
    
    std::ifstream input(filename.c_str());
    
    if (!input){
      std::cerr << "Cannot open file " << filename << std::endl;
      exit( EXIT_FAILURE);
    }
    
    std::string buffer;
    Read_state state=END;

    std::map<unsigned,Properties> properties_map;
    unsigned max_property=0;
    
    
    while(input){
      std::getline(input,buffer); //read a complete line
      
      if (!input) {
        if (input.eof()) break;
        std::cerr << "Fatal error: Error while reading file "<< filename <<"\n";
        exit( EXIT_FAILURE );
      }
      
      boost::trim(buffer); //remove white space at the beginning and at the end of the line
      if (buffer.empty() || buffer[0] == '#') continue; //handle comments and empty lines
      
      switch (state){
        case END:
          if (buffer == "PROPERTIES")     {state=PROPERTIES;     break;}
          if (buffer == "CLASSIFICATION") {state=CLASSIFICATION; break;}
          if (buffer == "EXTRA")          {state=EXTRA;          break;}
          if (buffer == "DEFAULT")          {state=DEFAULT;          break;}
          std::cerr << "Fatal error: Error while reading file "<< filename;
          std::cerr << ", unexpected <|" << buffer << "|> found.\n";
          exit( EXIT_FAILURE );
        case CLASSIFICATION:
        {
          if (buffer=="END")              {state=END;            break;}
          std::stringstream ss (std::stringstream::in | std::stringstream::out);
          ss << buffer;
          //inside a function
          unsigned index=Properties::add_classification(ss,hmap_);
          if (max_property <  index) max_property=index;
          break;
        }
        case PROPERTIES:
        {
          if (buffer=="END")              {state=END;            break;}
          std::stringstream ss (std::stringstream::in | std::stringstream::out);
          ss << buffer;
          unsigned index;
          ss >> index;
          //inside a function
          assert(properties_map.find(index)==properties_map.end());
          properties_map.insert(std::make_pair(index,Properties(ss,index)));
          break;
        }        
        case DEFAULT:
        {
          if (buffer=="END")              {state=END;            break;}
          unsigned i;
          std::stringstream ss (std::stringstream::in | std::stringstream::out);
          ss << buffer;
          ss >> i;
          if ( max_property<i ) max_property=i;
          Properties::index_of_default()=i;
          break;
        }
        case EXTRA:
        {
          if (buffer=="END")              {state=END;            break;}
          //call some function
          std::stringstream ss (std::stringstream::in | std::stringstream::out);
          ss << buffer;          
          Properties::handle_extra(ss);
          break;
        }
      };
    }
    
    
    unsigned i=0;
    for (typename std::map<unsigned,Properties>::const_iterator it=properties_map.begin();it!=properties_map.end();++it){
      if (i++ != it->first){
        std::cerr << "Fatal error: properties number " << i-1 << " could not be found in file " << filename <<"\n";
        exit( EXIT_FAILURE );
      }
      properties_.push_back(it->second);
    }  
    
    if (i <= max_property){
      std::cerr << "Fatal error: Up to " << max_property+1 << " properties are used, but only " << i << " are declared.\n";
      exit( EXIT_FAILURE );    
    }
    
    max_index_=i;
    
  }
  
  const Properties& get_properties(unsigned i) const{
    assert(i <= max_index_);
    return properties_[i];
  }
    
  const Properties& get_properties(const Query_type& query) const{
    typename Internal_map_type::const_iterator it=hmap_.find( Properties_::make_key(query) );
    if ( it==hmap_.end() ){
      if (Properties::index_of_default()!=-1)
        return properties_[Properties::index_of_default()];
      std::cerr << "Fatal error: Could not find an entry for " << Properties_::make_key(query);
      std::cerr << " and no default have been defined.\n";
      exit( EXIT_FAILURE );
    }
    return properties_[it->second];
  }
  
  unsigned number_of_properties() const {return max_index_;}
  
  Properties_iterator properties_begin(){ return properties_.begin();}
  Properties_const_iterator properties_begin() const { return properties_.begin();}
  Properties_iterator properties_end(){ return properties_.end(); }
  Properties_const_iterator properties_end() const { return properties_.end(); }
  
};

//TODO better specify this class: value and Value_type are not in the concept.
template <class Atom_classifier>
struct Weight_of_atoms{
  
  Weight_of_atoms(const Atom_classifier* atom_classifier):atom_classifier_(atom_classifier){}
    
  typename Atom_classifier::Properties::Value_type
  operator()(const typename Atom_classifier::Query_type& atom) const {
    typename Atom_classifier::Properties::Value_type value=get_radius(atom_classifier_->get_properties(atom));
    return value*value;
  }
private:
  const Atom_classifier* atom_classifier_;
};

}//namespace ESBTL

#endif //ESBTL_ATOM_CLASSIFIER_H