assignBondOrderProcessor.h

Go to the documentation of this file.

#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_LIMITS_H
# include <BALL/COMMON/limits.h>
#endif

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

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

#include <map>
#include <vector>
#include <queue>

#ifdef BALL_HAS_LPSOLVE
    struct _lprec;
    typedef struct _lprec lprec;
#endif

namespace BALL
{
  class BALL_EXPORT AssignBondOrderProcessor
    : public UnaryProcessor<AtomContainer>
  {
    protected:
      class Solution_;
      friend class Solution;
      class PQ_Entry_;
      friend class PQ_Entry_;

    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* HEURISTIC;

        static const char* INIFile;

        static const char* MAX_BOND_ORDER;

        static const char* MAX_NUMBER_OF_SOLUTIONS;

        static const char* COMPUTE_ALSO_NON_OPTIMAL_SOLUTIONS;

        static const char* BOND_LENGTH_WEIGHTING;

        static const char* APPLY_FIRST_SOLUTION;

        static const char* GREEDY_K_SIZE;

        static const char* BRANCH_AND_BOUND_CUTOFF;


      };

      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 HEURISTIC;
        static const String INIFile;
        static const int MAX_BOND_ORDER;
        static const int MAX_NUMBER_OF_SOLUTIONS;
        static const bool COMPUTE_ALSO_NON_OPTIMAL_SOLUTIONS;
        static const float BOND_LENGTH_WEIGHTING;
        static const bool APPLY_FIRST_SOLUTION;
        static const int GREEDY_K_SIZE;
        static const float BRANCH_AND_BOUND_CUTOFF;
      };

      struct BALL_EXPORT Algorithm
      {
        static const String A_STAR;
        static const String ILP;
        static const String K_GREEDY;
        static const String BRANCH_AND_BOUND;
      };

      struct BALL_EXPORT Heuristic
      {
        static const String SIMPLE;
        static const String MEDIUM;
        static const String TIGHT;
      };



      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 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 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();

      float evaluatePenalty(AtomContainer* ac);


      Options options;


    protected:

      // Nested class storing the parameters of a solution to our ILP
      class Solution_
      {
        friend class AssignBondOrderProcessor;

        public:
          // Default constructor
          Solution_();

          // Detailed constructor 
          Solution_(PQ_Entry_& entry ,AssignBondOrderProcessor* abop,
                    int number_of_node_expansions, int search_queue_size);

          // Destructor
          virtual ~Solution_();

          // 
          void clear();

          // 
          int getNumberOfNodeExpansions() const {return node_expansions;}

          //
          int getQueueSize() const {return queue_size;}

          // denotes whether the problem could be solved or not  
          bool valid;

          // the result : the set of bond orders for _ALL_ original bonds
          HashMap<Bond*, int> bond_orders;

          // the result part2: the atoms with n additional hydrogens
          HashMap<Atom*, int> number_of_virtual_hydrogens;

          // the virtual atoms and bonds that should be deleted when the next 
          // solution is applied
          vector<Atom*> atoms_to_delete;
          //vector<Bond*> bonds_to_delete;

          // the values of the objective function
          float atom_type_penalty;
          float bond_length_penalty;
          float total_charge;
          int node_expansions;
          int queue_size;
          AssignBondOrderProcessor* parent;
          AtomContainer* ac;
      };

      // Nested class storing a priority queue entry for the A-STAR-Option
      class PQ_Entry_
      {
        friend class AssignBondOrderProcessor;

        public:

          // Default constructor
          PQ_Entry_();

          // Destructor
          ~PQ_Entry_();

          // 
          void clear();

          // the less operator.
          // NOTE: we want a reverse sort, hence we actually return a "greater" 
          bool operator < (const PQ_Entry_& b) const;
            // the penalty

          float coarsePenalty(float atom_type_penalty, float bond_length_penalty) const
          {
            return ( ( (    (atom_type_normalization_factor < 0.0001)
                         || (bond_length_normalization_factor < 0.0001)
                         || (alpha < 0.0001)) ?
                       atom_type_penalty :
                       ((1.-alpha) * (atom_type_penalty / atom_type_normalization_factor)
                        + (alpha* bond_length_penalty / bond_length_normalization_factor))));
          }

          // the penalty
          float coarsePenalty() const
          {
            return coarsePenalty(estimated_atom_type_penalty, estimated_bond_length_penalty);
          }

          // the bond length penalty 
          float finePenalty() const {return estimated_bond_length_penalty;}

          // the estimated atom type penalty
          float estimated_atom_type_penalty;
          // the estimated bond length penalty
          float estimated_bond_length_penalty;

          // the bond orders 
          // the i-th entry denotes the bondorder of the i-th bond
          // unset bonds get the order 0
          vector<short> bond_orders;

          // the index of the bond last considered 
          Position last_bond;

          // NOTE: these variables are here to speed up
          //       a number of decisions, but we do not
          //       want to store them in each PQ_Entry_
          //       object since they are the same in the
          //       whole queue
          //
          //       But this approach is not thread-safe!
          //       Two simultaneously running bond order threads
          //       might overwrite their alpha, ... values!
          static float alpha;
          static float atom_type_normalization_factor;
          static float bond_length_normalization_factor;
          static bool  use_fine_penalty;
        };

      bool readOptions_();


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

      bool preassignPenaltyClasses_();

      int getPenaltyClass_(Atom* atom);


      bool precomputeBondLengthPenalties_();

      float computeVirtualHydrogens_(Atom* atom);

      bool apply_(Solution_& 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 Solution_& sol){return sol.getQueueSize();}

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

      float getTotalPenalty_(const Solution_& sol)
      {
        if (   (atom_type_normalization_factor_   < 0.00001)
            || (bond_length_normalization_factor_ < 0.00001)
            || (alpha_ < 0.0001))
        {
          return sol.atom_type_penalty;
        }
        else
        {
          return (  (1.-alpha_) * (sol.atom_type_penalty/atom_type_normalization_factor_)
                  + (alpha_*sol.bond_length_penalty/bond_length_normalization_factor_));
        }
      }

      /* 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 Solution_& sol){return sol.getNumberOfNodeExpansions();}


#ifdef BALL_HAS_LPSOLVE

      bool createILP_();

      bool solveILP_();
#endif

      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_;
      //
      //

      // ********************* ILP stuff ***********************
      //
      // Vector for mapping from variable indices onto free bonds in the
      // order used by the ILP
      std::vector<Bond*> ilp_index_to_free_bond_;

      // number of bond variables in the ILP
      Position ilp_number_of_free_bonds_;

      // Constant penalty (fixed bonds)
      float ilp_const_penalty_;

      // ******************* 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<Solution_> solutions_;

      // the original conformation before we computed anything
      // this is a vector because we can have multiple molecules...
      vector<Solution_> 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_;

      // 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_;

      // this enum allows faster access to the type of the chosen heuristic than a string compare
      enum HEURISTIC_INDEX
      {
        SIMPLE,
        MEDIUM,
        TIGHT
      };

      HEURISTIC_INDEX heuristic_index_;

      // //////// ************ for Algorithm::BRANCH_AND_BOUND ************ /////////
      bool performBranchAndBound_();
      float greedy_atom_type_penalty_;
      float greedy_bond_length_penalty_;

      // //////// ************ for Algorithm::K_GREEDY ************ /////////
      vector<PQ_Entry_> performGreedy_(PQ_Entry_& entry, Size greedy_k  = 10);
      int greedy_node_expansions_;

      // //////// ************ for Algorithm::A_STAR ************ /////////
      bool  performAStarStep_();


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

      std::priority_queue<PQ_Entry_> queue_;

      bool estimatePenalty_(PQ_Entry_& entry, bool include_heuristic_term = true);

      float estimateAtomTypePenalty_(Atom* atom,
                                     Index atom_index,    // the atom index
                                     int fixed_valence,   // its so far fixed valence (incl. virtual H's)
                                     int fixed_virtual_order, // its so far fixed virtual H's
                                     int num_free_bonds,  // its number of unfixed original bonds
                                     PQ_Entry_& entry);
      //  NOTE: virtual bonds are excluded!
      float estimateBondLengthPenalty_(Index atom_index, // the atom index
                                       const vector<Bond*>& free_bonds,
                                       int fixed_virtual_order,
                                       int fixed_valence,
                                       int num_free_bonds);

      // 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_;

      // The current number of node expansions. 
      // step_ + queue_.size() gives the number of touched nodes.
      int step_;

      Timer timer_;
#ifdef BALL_HAS_LPSOLVE
      lprec* ilp_;
#endif
      AssignBondOrderProcessor(const AssignBondOrderProcessor& abop);
      AssignBondOrderProcessor& operator = (const AssignBondOrderProcessor& abop);
    };

} // namespace BALL 


#endif // BALL_STRUCTURE_ASSIGNBONDORDERPROCESSOR_H

---