poseClustering.h

Go to the documentation of this file.

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

#ifndef BALL_DOCKING_POSECLUSTERING_H
#define BALL_DOCKING_POSECLUSTERING_H

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

#ifndef BALL_DOCKING_COMMON_CONFORMATIONSET_H
# include <BALL/DOCKING/COMMON/conformationSet.h>
#endif

#ifndef BALL_MOLMEC_COMMON_SNAPSHOT_H
# include <BALL/MOLMEC/COMMON/snapShot.h>
#endif

#ifndef BALL_STRUCTURE_ATOMBIJECTION_H
# include <BALL/STRUCTURE/atomBijection.h>
#endif

#ifndef BALL_KERNEL_SYSTEM_H
# include <BALL/KERNEL/system.h>
#endif

#ifndef BALL_DATATYPE_STRING_H
# include <BALL/DATATYPE/string.h>
#endif

#ifndef BALL_MATHS_VECTOR3_H
#       include <BALL/MATHS/vector3.h>
#endif

#ifndef BALL_MATHS_MATRIX44_H
# include <BALL/MATHS/matrix44.h>
#endif

#include <Eigen/Core>

#include <boost/shared_ptr.hpp>
#include <boost/variant.hpp>
#include <boost/graph/adjacency_list.hpp>

#include <set>
#include <iostream>

#ifdef BALL_HAS_TBB
# include <tbb/parallel_reduce.h>
# include <tbb/blocked_range.h>
#endif

//#define POSECLUSTERING_DEBUG 1
#undef POSECLUSTERING_DEBUG

namespace BALL
{
  class BALL_EXPORT PoseClustering
  {
    public:
      struct BALL_EXPORT Option
      {
        static const String CLUSTER_METHOD;

        static const String DISTANCE_THRESHOLD;

        static const String RMSD_LEVEL_OF_DETAIL;

        static const String RMSD_TYPE;

        static const String RUN_PARALLEL;

      };

      struct BALL_EXPORT Default
      {
        static const Index CLUSTER_METHOD;
        static const float DISTANCE_THRESHOLD;
        static const Index RMSD_LEVEL_OF_DETAIL;
        static const Index RMSD_TYPE;
        static const bool  RUN_PARALLEL;
        static const bool USE_CENTER_OF_MASS_PRECLINK;
      };

      enum BALL_EXPORT RMSDType
      {
        SNAPSHOT_RMSD,
        RIGID_RMSD,
        CENTER_OF_MASS_DISTANCE
      };

      enum BALL_EXPORT RMSDLevelOfDetail
      {
        C_ALPHA, //=0
        HEAVY_ATOMS,
        BACKBONE,
        ALL_ATOMS,
        PROPERTY_BASED_ATOM_BIJECTION
      };

      enum BALL_EXPORT ClusterMethod
      {
        TRIVIAL_COMPLETE_LINKAGE,
        SLINK_SIBSON,
        CLINK_DEFAYS,
        NEAREST_NEIGHBOR_CHAIN_WARD,
        CLINK_ALTHAUS
      };

      class BALL_EXPORT RigidTransformation
      {
        public:
          RigidTransformation() {};

          RigidTransformation(Eigen::Vector3f const& new_trans, Eigen::Matrix3f const& new_rot)
            : translation(new_trans),
              rotation(new_rot)
          {}

          Eigen::Vector3f translation;
          Eigen::Matrix3f rotation;
      };

      class BALL_EXPORT PosePointer
      {
        public:
          PosePointer(RigidTransformation const* new_trafo, SnapShot const* new_snap = 0)
            : trafo(new_trafo),
              snap(new_snap)
          { }

          PosePointer(SnapShot const* new_snap)
            : trafo(0),
              snap(new_snap)
          { }

          RigidTransformation const* trafo;
          SnapShot            const* snap;
      };

      class BALL_EXPORT ClusterProperties
      {
        public:
          template <class Archive>
          void serialize(Archive& ar, const unsigned int version);

          std::set<Index> poses;

          Size size;

          boost::variant<Eigen::VectorXf, RigidTransformation> center;

          float merged_at;
#ifdef POSECLUSTERING_DEBUG 

          float current_cluster_id;
#endif
      };

      typedef boost::adjacency_list<boost::vecS,
                                    boost::vecS,
                                    boost::directedS,
                                    ClusterProperties,
                                    boost::no_property,
                                    unsigned int> ClusterTree;

      typedef ClusterTree::vertex_descriptor ClusterTreeNode;

      BALL_CREATE(PoseClustering);


      PoseClustering();

      PoseClustering(ConformationSet* poses, float rmsd);

      PoseClustering(System const& base_system, String transformation_file_name);

      virtual ~PoseClustering();



      bool compute();



      void setConformationSet(ConformationSet* new_set, bool precompute_atombijection = false);

      void setBaseSystemAndPoses(System const& base_system, std::vector<PosePointer> const& poses);

      //       will be applied upon the current conformation
      void setBaseSystemAndTransformations(System const& base_system, String transformation_file_name);

      const ConformationSet* getConformationSet() const {return  current_set_;}

      ConformationSet* getConformationSet() {return  current_set_;}

      const std::vector<RigidTransformation> & getRigidTransformations() const {return transformations_;}

      std::vector<Vector3> & getCentersOfMass() {return com_;}

      std::vector<Vector3> const & getCentersOfMass() const {return com_;}

      const System& getSystem() const;

      System& getSystem();

      Size getNumberOfPoses() const {return poses_.size();}

      Size getNumberOfClusters() const {return clusters_.size();}

      const std::set<Index>& getCluster(Index i) const;

      std::set<Index>& getCluster(Index i);

      Size getClusterSize(Index i) const;

      float getClusterScore(Index i) const;

      float getScore(const System sys_a, const System sys_b, Options options) const;

      AtomBijection& getAtomBijection() {return atom_bijection_;}

      AtomBijection const& getAtomBijection() const {return atom_bijection_;}

      void applyTransformation2System(Index i, System& target_system);

      void convertTransformations2Snaphots();

      void convertSnaphots2Transformations();

      float computeCompleteLinkageRMSD(Index i, Options options, bool initialize = true);

      //float computeCompleteLinkageRMSD(boost::shared_ptr<ConformationSet> cluster, Option options) const;

      boost::shared_ptr<System> getPose(Index i) const;

      std::vector<PosePointer> const&  getPoses() const {return poses_;}

      boost::shared_ptr<System> getClusterRepresentative(Index i);

      Index findClusterRepresentative(Index i);

      boost::shared_ptr<ConformationSet> getClusterConformationSet(Index i);

      boost::shared_ptr<ConformationSet> getReducedConformationSet();

      bool refineClustering(Options const& refined_options);


      Options options;

      void setDefaultOptions();



      static float getRigidRMSD(Eigen::Vector3f const& t_ab, Eigen::Matrix3f const& M_ab, Eigen::Matrix3f const& covariance_matrix);

      static float getSquaredRigidRMSD(Eigen::Vector3f const& t_ab, Eigen::Matrix3f const& M_ab, Eigen::Matrix3f const& covariance_matrix);

      static Eigen::Matrix3f computeCovarianceMatrix(System const& system, Index rmsd_level_of_detail = C_ALPHA);



      std::vector<std::set<Index> > extractClustersForThreshold(float threshold, Size min_size = 0);

      std::vector<std::set<Index> > extractNBestClusters(Size n);

      std::vector<std::set<Index> > filterClusters(Size min_size = 1);

      void serializeWardClusterTree(std::ostream& out, bool binary = false);

      void deserializeWardClusterTree(std::istream& in, bool binary = false);

      void exportWardClusterTreeToGraphViz(std::ostream& out);



      void printClusters(std::ostream& out = std::cout) const;

      void printClusterScores(std::ostream& out = std::cout);


    protected:

#ifdef BALL_HAS_TBB

      class ComputeNearestClusterTask_
      {
        public:
          ComputeNearestClusterTask_(PoseClustering* parent,
                                     const std::vector<ClusterTreeNode>& active_clusters,
                                     Position current_cluster,
                                     Index rmsd_type);

          ComputeNearestClusterTask_(ComputeNearestClusterTask_& cnct, tbb::split);

          void join(ComputeNearestClusterTask_ const& cnct);

          void operator() (const tbb::blocked_range<size_t>& r);

          Position getMinIndex() {return my_min_index_;}

          float getMinValue() {return my_min_value_;}

        protected:
          // the PoseClustering instance that called us
          PoseClustering* parent_;

          // the array we work on
          const std::vector<ClusterTreeNode>& active_clusters_;

          // the cluster to compare to everything else
          Position current_cluster_;

          // the kind of rmsd computation desired
          Index rmsd_type_;

          // the minimum index in our own block
          Position my_min_index_;

          // the minimum value in our own block
          float my_min_value_;
      };
#endif

      class ClusterTreeWriter_
      {
        public:
          ClusterTreeWriter_(ClusterTree const* cluster_tree)
            : cluster_tree_(cluster_tree)
          { }

          void operator() (std::ostream& out, const ClusterTreeNode& v) const;

        protected:
          ClusterTree const* cluster_tree_;
      };

      class ClusterTreeNodeComparator
      {
        public:
          ClusterTreeNodeComparator(ClusterTree& cluster_tree)
            : cluster_tree_(&cluster_tree)
          {}

          bool operator() (ClusterTreeNode const first, ClusterTreeNode const second) const
          {
            float first_value  = (*cluster_tree_)[ first].merged_at;
            float second_value = (*cluster_tree_)[second].merged_at;

            return first_value < second_value;
          }

        protected:
          ClusterTree* cluster_tree_;
      };

      // trivial complete linkage implementation
      // with O(n^2) space request
      bool trivialCompute_();

      // space efficient (SLINK or CLINK) clustering
      bool linearSpaceCompute_();

      //
      bool althausCompute_();

      //  implementation of a single linkage clustering as described in 
      //       R. Sibson: SLINK: an optimally efficient algorithm for the single-link cluster method. 
      //       The Computer Journal. 16, 1, British Computer Society, 1973, p. 30-34
      void slinkInner_(int current_level);

      //  implementation of a complete linkage clustering as described in 
      //          D. Defays: An efficient algorithm for a complete link method. 
      //          The Computer Journal. 20, 4, British Computer Society, 1977, p. 364-366. 
      void clinkInner_(int current_level);

      // implememtation of the nearest neighbor chain clustering algorithm 
      // as described in:
      //          Murtagh, Fionn (1983): "A survey of recent advances in hierarchical clustering algorithms", 
      //          The Computer Journal 26 (4): 354–359
      bool nearestNeighborChainCompute_();

      void initWardDistance_(Index rmsd_type);

      void updateWardDistance_(ClusterTreeNode parent, ClusterTreeNode i, ClusterTreeNode j, Index rmsd_type);

      float computeWardDistance_(ClusterTreeNode i, ClusterTreeNode j, Index rmsd_type);

      std::set<Index> collectClusterBelow_(ClusterTreeNode const& v);

      // compute the center of masses
      void computeCenterOfMasses_();

      // precompute an atom bijection for faster access
      void precomputeAtomBijection_();

      // check the given atom wrt choice of option RMSD_LEVEL_OF_DETAIL
      bool static isExcludedByLevelOfDetail_(Atom const* atom, Index rmsd_level_of_detail);

      // distance between cluster i and cluster j
      float getClusterRMSD_(Index i, Index j, Index rmsd_type);

      // reads the poses given as transformations from a file
      // Note: the previously given system will be taken
      //       as untransformed reference, e.g. all transformations 
      //       will be applied upon the current conformation
      bool readTransformationsFromFile_(String filename);

      // compute the RMSD between two "poses"  
      float getRMSD_(Index i, Index j, Index rmsd_type);

      // store pointers to the snapshots in the poses vector
      void storeSnapShotReferences_();

      //
      void printCluster_(Index i, std::ostream& out = std::cout) const;

      // 
      void printVariables_(int a, int b, double c, int d, double e, int current_level);

      //
      void clear_();

      // only used by trivial clustering
      Eigen::MatrixXd                  pairwise_scores_;

      ConformationSet*                 current_set_;

      std::vector< std::set<Index> >   clusters_;

      std::vector< Index >   cluster_representatives_;

      std::vector< float >             cluster_scores_;

      Index                            rmsd_level_of_detail_;


      // ----- unified access to the poses, independent of their type
      //       (the poses are either stored in transformations_, or in current_set_)
      std::vector<PosePointer>         poses_;

      // ----- data structure for transformation input (instead of snapshots)
      std::vector<RigidTransformation> transformations_;

      Eigen::Matrix3f                  covariance_matrix_;

      // TODO: maybe use a const - ptr instead?
      System                           base_system_;

      // the reference state
      SnapShot                         base_conformation_;

      // flag indicating the use of transformation as input
      bool                             has_rigid_transformations_;

      // do we need to delete the conformation set, that was
      // created by converting transformations to snapshots
      bool                             delete_conformation_set_;

      // ------ data structures for slink and clink 

      // stores the distance at which this indexed element has longer 
      // the largest index of its cluster
      std::vector<double>    lambda_;

      // the index of the cluster representative at merge-time 
      // (element with largest index)
      std::vector<int>       pi_;

      std::vector<double>    mu_;


      // ----- data structures for nearest neighbor chain ward
      Size                   number_of_selected_atoms_;

      // ----- data structure for CENTER_OF_GRAVITY_CLINK

      // the geometric center of mass
      std::vector<Vector3>   com_;

      // ----- general data structures

      // We cache the atom bijection for faster
      // RMSD computation; this is possible, since the system topology does
      // not change
      AtomBijection     atom_bijection_;

      // helper dummies to speed up snapshot application
      System            system_i_;
      System            system_j_;

      ClusterTree       cluster_tree_;
  }; //class PoseClustering
} //namesspace BALL

#endif // BALL_DOCKING_POSECLUSTERING_H