assignBondOrderProcessor.h

Go to the documentation of this file.

// -*- Mode: C++; tab-width: 2; -*-
// vi: set ts=2:
//

#ifndef BALL_STRUCTURE_ASSIGNBONDORDERPROCESSOR_H
#define BALL_STRUCTURE_ASSIGNBONDORDERPROCESSOR_H

#ifndef BALL_CONCEPT_PROCESSOR_H
  #include <BALL/CONCEPT/processor.h>
#endif

#ifndef BALL_KERNEL_ATOMCONTAINER_H
  #include <BALL/KERNEL/atomContainer.h>
#endif

#ifndef BALL_DATATYPE_HASHMAP_H
  #include <BALL/DATATYPE/hashMap.h>
#endif

#ifndef BALL_DATATYPE_HASHSET_H
  #include <BALL/DATATYPE/hashSet.h>
#endif

#ifndef BALL_KERNEL_BOND_H
  #include <BALL/KERNEL/bond.h>
#endif

#ifndef BALL_DATATYPE_OPTIONS_H
# include <BALL/DATATYPE/options.h>
#endif

#ifndef BALL_COMMON_EXCEPTION_H
# include <BALL/COMMON/exception.h>
#endif

#ifndef BALL_SYSTEM_TIMER_H
# include <BALL/SYSTEM/timer.h>
#endif

#ifndef BALL_STRUCTURE_BONDORDERS_BONDORDERASSIGNMENTSTRATEGY_H
# include <BALL/STRUCTURE/BONDORDERS/bondOrderAssignmentStrategy.h>
#endif

#ifndef BALL_STRUCTURE_BONDORDERS_BONDORDERASSIGNMENT_H
# include <BALL/STRUCTURE/BONDORDERS/bondOrderAssignment.h>
#endif

#ifndef BALL_STRUCTURE_BONDORDERS_PARTIALBONDORDERASSIGNMENT_H
# include <BALL/STRUCTURE/BONDORDERS/partialBondOrderAssignment.h>
#endif

#include <map>
#include <vector>

namespace BALL
{
  class BALL_EXPORT AssignBondOrderProcessor
    : public UnaryProcessor<AtomContainer>
  {
    protected:
      friend class PartialBondOrderAssignment;
      friend class BondOrderAssignment;

      class PQ_Entry_;
      friend class PQ_Entry_;

      friend class BondOrderAssignmentStrategy;

      friend class AStarBondOrderStrategy;
      friend class BranchAndBoundBondOrderStrategy;
      friend class FPTBondOrderStrategy;
      friend class ILPBondOrderStrategy;
      friend class KGreedyBondOrderStrategy;

    public:

      struct BALL_EXPORT Option
      {
        static const char* OVERWRITE_SINGLE_BOND_ORDERS;

        static const char* OVERWRITE_DOUBLE_BOND_ORDERS;

        static const char* OVERWRITE_TRIPLE_BOND_ORDERS;

        static const char* OVERWRITE_SELECTED_BONDS;

        static const char* ADD_HYDROGENS;             //TODO

        static const char* COMPUTE_ALSO_CONNECTIVITY; //TODO

        static const char* CONNECTIVITY_CUTOFF;       //TODO 

        static const char* USE_FINE_PENALTY;

        static const char* KEKULIZE_RINGS;

        static const char* ALGORITHM;

        static const char* INIFile;

        static const char* MAX_BOND_ORDER;

        static const char* MAX_NUMBER_OF_SOLUTIONS;

        static const char* MAX_PENALTY;

        static const char* COMPUTE_ALSO_NON_OPTIMAL_SOLUTIONS;

        static const char* BOND_LENGTH_WEIGHTING;

        static const char* APPLY_FIRST_SOLUTION;

      };

      struct BALL_EXPORT Default
      {
        static const bool OVERWRITE_SINGLE_BOND_ORDERS;
        static const bool OVERWRITE_DOUBLE_BOND_ORDERS;
        static const bool OVERWRITE_TRIPLE_BOND_ORDERS;
        static const bool OVERWRITE_SELECTED_BONDS;
        static const bool ADD_HYDROGENS;
        static const bool COMPUTE_ALSO_CONNECTIVITY;
        static const float CONNECTIVITY_CUTOFF;
        static const bool USE_FINE_PENALTY;
        static const bool KEKULIZE_RINGS;
        static const String ALGORITHM;
        static const String INIFile;
        static const int MAX_BOND_ORDER;
        static const int MAX_NUMBER_OF_SOLUTIONS;
        static const int MAX_PENALTY;
        static const bool COMPUTE_ALSO_NON_OPTIMAL_SOLUTIONS;
        static const float BOND_LENGTH_WEIGHTING;
        static const bool APPLY_FIRST_SOLUTION;
      };

      struct BALL_EXPORT Algorithm
      {
        static const String A_STAR;

        static const String ILP;

        static const String FPT;

        static const String K_GREEDY;
        static const String BRANCH_AND_BOUND;
      };



      AssignBondOrderProcessor();

      // constructor with parameter filename //TODO
      //AssignBondOrderProcessor(const String& file_name) throw(Exception::FileNotFound);

      virtual ~AssignBondOrderProcessor();


      virtual bool start();

      void clear();

      virtual Processor::Result operator ()(AtomContainer& ac);

      virtual bool finish();



      Size getNumberOfAddedHydrogens(Position i)
      {
        if (i >= solutions_.size())
        {
          Log.error() << "AssignBondOrderProcessor: No solution with index " << i << std::endl;
          return 0;
        }
        int num_hydrogens = 0;

        HashMap<Atom*, int>::Iterator it = solutions_[i].number_of_virtual_hydrogens.begin();
        for (; it != solutions_[i].number_of_virtual_hydrogens.end(); it++)
          num_hydrogens += it->second;
        return num_hydrogens;
      }

      Size getNumberOfComputedSolutions() {return solutions_.size();}


      AtomContainer* getAtomContainer() {return ac_;}

      AtomContainer const* getAtomContainer() const {return ac_;}

      const System& getSolution(Position i) throw(Exception::IndexOverflow);

      float getTotalCharge(Position i)
      {
        if (i >= solutions_.size())
        {
          Log.error() << "AssignBondOrderProcessor: No solution with index " << i << std::endl;

          return std::numeric_limits<float>::max();
        }
        else
        {
          return getTotalCharge_(solutions_[i]);
        }
      }

      float getTotalPenalty(Position i=0)
      {
        if (i >= solutions_.size())
        {
          Log.error() << "AssignBondOrderProcessor: No solution with index " << i << std::endl;

          return std::numeric_limits<float>::max();
        }
        else
          return getTotalPenalty_(solutions_[i]);
      }

      /* Returns the number of node expansions before solution i was found.
       *
       * param    i  index of the solution, whose number of node expansions should be returned.
       * return  int -   number of node expansions before solution i was found.   
       */
      int getNumberOfNodeExpansions(Position i)
      {
        if (i >= solutions_.size())
        {
          Log.error() << "AssignBondOrderProcessor: No solution with index " << i << std::endl;

          return -1;
        }
        else
          return getNumberOfNodeExpansions_(solutions_[i]);
      }

      /* Returns the number of node expansions before solution i was found.
       *
       * param    i  index of the solution, whose  queue size should be returned. 
       * return  int -  queue size when solution i was found.  
       */
      int getQueueSize(Position i)
      {
        if (i >= solutions_.size())
        {
          Log.error() << "AssignBondOrderProcessor: No solution with index " << i << std::endl;

          return -1;
        }
        else
          return getQueueSize_(solutions_[i]);
      }

      bool apply(Position i);

      void resetBondOrders();

      bool computeNextSolution(bool apply_solution = true);

      void setDefaultOptions();

      bool hasValidOptions(){return readOptions_();}

      float evaluatePenalty(AtomContainer* ac);

      Options options;


    protected:

      bool readOptions_();


      bool readAtomPenalties_() throw(Exception::FileNotFound());

      bool preassignPenaltyClasses_();

      int getPenaltyClass_(Atom* atom);


      bool precomputeBondLengthPenalties_();

      float computeVirtualHydrogens_(Atom* atom);

      //TODO: move to solution!
      bool apply_(BondOrderAssignment& solution);

      void storeOriginalConfiguration_();

      /* Returns the queue's size at the moment the given solution was found.
       *
       * param   sol  solution, whose queue size should be returned. 
       * return  int -  queue size when the given solution was found.  
       */
      int getQueueSize_(const BondOrderAssignment& sol){return sol.getQueueSize();}

      float getTotalCharge_(const BondOrderAssignment& sol)
      {
        if (sol.valid)
        {
          return sol.total_charge;
        }
        else
        {
          return 0;
        }
      }

      float getTotalPenalty_(const BondOrderAssignment& sol)
      {
        return sol.coarsePenalty();
      }

      /* Returns the number of node expansions before the given solution was found.
       *
       * param   sol  solution, whose number of node expansions should be returned. 
       * return  int -  number of node expansions before solution i was found.  
       */
      int getNumberOfNodeExpansions_(const BondOrderAssignment& sol){return sol.getNumberOfNodeExpansions();}

      bool valid_;

      bool evaluation_mode_;

      // Map for storing the bonds fixed orders
      // if a bond is free, the map returns 0
      std::map<Bond*, short> bond_fixed_;

      // all free bonds in the atom container
      std::vector<Bond*> free_bonds_;

      // Map for storing the bonds associated index (all bonds)
      HashMap<Bond*, Index> bond_to_index_;

      // Vector for mapping from variable indices onto bonds (all bonds)
      std::vector<Bond*> index_to_bond_;



      // ***************** datastructures for virtual hydrogen bonds ****************** 
      //
      //  NOTE: a single virtual bond represents ALL possible hydrogen 
      //        bonds for a given atom
      //
      // the atoms with upto n possible additional hydrogens
      HashMap<Atom*, int> number_of_virtual_hydrogens_;
      //
      // the max number of virtual hydrogens per virtual bond index
      std::vector<int> virtual_bond_index_to_number_of_virtual_hydrogens_;
      //  
      // the number of virtual bonds
      Size num_of_virtual_bonds_;
      //
      // the virtual bond index assigned to this atom!
      vector<Atom*> virtual_bond_index_to_atom_;
      HashMap<Atom*, int> atom_to_virtual_bond_index_;
      //
      //
      // a virtual dummy bond
      Bond* virtual_bond_;

      // ******************* general datastructures *********************

      // the number of bonds given (free + fixed!)
      Position total_num_of_bonds_;

      // num of free bonds without virtual bonds!
      int num_of_free_bonds_;

      // store for all atom-indices the atoms fixed valences 
      std::vector<Position> fixed_val_;

      // storing the solutions
      vector<BondOrderAssignment> solutions_;

      // the original conformation before we computed anything
      // this is a vector because we can have multiple molecules...
      vector<BondOrderAssignment> starting_configuration_;

      // the inverse of the atom type penalty normalization factor
      float atom_type_normalization_factor_;

      // the inverse of the bond length penalty normalization factor
      float bond_length_normalization_factor_;

      // denotes the index of the last applied solution
      // -1 if there was no valid solution applied
      int last_applied_solution_;

      // the AtomContainer, the processor is operating on
      AtomContainer* ac_;

      // max bond order to consider
      int max_bond_order_;

      // balance parameter between atom type and bond length penalty
      float alpha_;

      // the max number of solutions to compute 
      int max_number_of_solutions_;

      // the max penalty score 
      int max_penalty_;

      // flag to indicate, whether also non-optimal solutions should be computed 
      bool compute_also_non_optimal_solutions_;

      // flag for adding missing hydrogens
      bool add_missing_hydrogens_;

      // flag for computing also the bond connectivity
      bool compute_also_connectivity_;

      // flag for using fine penalties derived from 3d information
      bool use_fine_penalty_;

      // ////////              general stuff                      /////////


      // The penalty administration datastructures.
      //  filled by readAtomPenalties_
      //  organized in imaginarey blocks of length  
      //  block_to_length_[i], starting from 
      //  block_to_start_idx_[i] associating 
      //  block_to_start_valence_[i] to the start_idx
      vector<int> penalties_;
      vector<Position> block_to_start_idx_;
      vector<Size> block_to_length_;
      vector<int> block_to_start_valence_;
      // stores the defining element and the SMART-string of each block
      vector<std::pair<String, String> > block_definition_;


      // Stores which atom belongs to which penalty block.
      // The first vector element of each atom block denotes the penalty block 
      // assigned to the atom without any additional VIRTUAL Hydrogens,
      // the second element with one additional Hydrogen and so on. 
      vector< vector<int> > atom_to_block_;

      // Stores the possible bond lengths penalties per order.
      HashMap<Bond*, vector<float> > bond_lengths_penalties_;

      Timer timer_;

      AssignBondOrderProcessor(const AssignBondOrderProcessor& abop);
      AssignBondOrderProcessor& operator = (const AssignBondOrderProcessor& abop);

      // The strategies this class can use
      StringHashMap<boost::shared_ptr<BondOrderAssignmentStrategy> > strategies_;
    };

} // namespace BALL 


#endif // BALL_STRUCTURE_ASSIGNBONDORDERPROCESSOR_H