HierarchicalClustering.h

Go to the documentation of this file.

// --------------------------------------------------------------------------
//                   OpenMS -- Open-Source Mass Spectrometry
// --------------------------------------------------------------------------
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
// ETH Zurich, and Freie Universitaet Berlin 2002-2012.
//
// This software is released under a three-clause BSD license:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name of any author or any participating institution
//    may be used to endorse or promote products derived from this software
//    without specific prior written permission.
// For a full list of authors, refer to the file AUTHORS.
// --------------------------------------------------------------------------
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
// INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// --------------------------------------------------------------------------
// $Maintainer: Lars Nilse $
// $Authors: Bastian Blank $
// --------------------------------------------------------------------------

#include <cmath>
#include <limits>
#include <map>
#include <queue>
#include <boost/unordered/unordered_set.hpp>

#include <OpenMS/COMPARISON/CLUSTERING/HashGrid.h>
#include <OpenMS/CONCEPT/Types.h>

#ifndef OPENMS_COMPARISON_CLUSTERING_HIERARCHICALCLUSTERING_H
#define OPENMS_COMPARISON_CLUSTERING_HIERARCHICALCLUSTERING_H

namespace OpenMS
{
  template <typename PointRef>
  class HierarchicalClustering
  {
public:
    typedef DPosition<2, DoubleReal> PointCoordinate;

    class BoundingBox :
      public std::pair<PointCoordinate, PointCoordinate>
    {
public:
      BoundingBox(const PointCoordinate & p) :
        std::pair<PointCoordinate, PointCoordinate>(std::make_pair(p, p))
      {}

      BoundingBox(const BoundingBox & b) :
        std::pair<PointCoordinate, PointCoordinate>(b)
      {}

      PointCoordinate size() const
      {
        return this->second - this->first;
      }

      BoundingBox & operator|=(const BoundingBox & rhs)
      {
        typename PointCoordinate::iterator lit;
        typename PointCoordinate::const_iterator rit;

        // Calculate lower bound
        lit = this->first.begin(); rit = rhs.first.begin();
        for (; lit != this->first.end(); ++lit, ++rit) *lit = std::min(*lit, *rit);

        // Calculate upper bound
        lit = this->second.begin(); rit = rhs.second.begin();
        for (; lit != this->second.end(); ++lit, ++rit) *lit = std::max(*lit, *rit);

        return *this;
      }

      BoundingBox operator|(const BoundingBox & rhs) const
      {
        BoundingBox ret(*this);
        ret |= rhs;
        return ret;
      }

      operator PointCoordinate() const
      {
        // (first + second) / 2
        return coordScalarDiv_(this->first + this->second, 2);
      }
    };

    class Cluster :
      public boost::unordered_multimap<PointCoordinate, PointRef>
    {
public:
      BoundingBox bbox;

      Cluster(const BoundingBox & bbox) :
        bbox(bbox)
      {}
    };

    typedef HashGrid<Cluster> Grid;

    Grid grid;

protected:
    class TreeNode
    {
public:
      const PointCoordinate coord;
      const BoundingBox bbox;
      TreeNode * left, * right;
      UInt points;
      const bool center;
      const PointRef ref;

      TreeNode(const PointCoordinate & coord, const PointRef & ref, bool center) :
        coord(coord), bbox(coord), left(0), right(0), points(1), center(center), ref(ref)
      {}

      TreeNode(const PointCoordinate & coord, const BoundingBox & bbox, TreeNode * left, TreeNode * right) :
        coord(coord), bbox(bbox),
        left(left), right(right),
        points(left->points + right->points),
        center(left->center && right->center),
        ref(PointRef())
      {}
    };

    typedef std::map<typename Grid::CellIndex, std::pair<typename Grid::CellContent *, bool> > ClusterCells;
    typedef boost::unordered_set<TreeNode *> ClusterTrees;

    class TreeDistance
    {
public:
      DoubleReal distance;
      TreeNode * left, * right;

      TreeDistance(const DoubleReal & distance, TreeNode * left, TreeNode * right) :
        distance(distance), left(left), right(right)
      {}

      bool operator>(const TreeDistance & rhs) const
      {
        return distance > rhs.distance;
      }

    };

    typedef std::priority_queue<TreeDistance, std::vector<TreeDistance>, std::greater<TreeDistance> > TreeDistanceQueue;

public:
    HierarchicalClustering(const PointCoordinate & cluster_dimension) :
      grid(cluster_dimension)
    {}

    typename Grid::cell_iterator insertPoint(const PointCoordinate & d, const PointRef & ref)
    {
      typename Grid::cell_iterator it = insertCluster_(d);
      it->second.insert(std::make_pair(d, ref));
      return it;
    }

    void cluster()
    {
      // Collect coordinates of all active cells
      std::vector<typename Grid::CellIndex> cells;
      for (typename Grid::const_grid_iterator it = grid.grid_begin(); it != grid.grid_end(); ++it)
        cells.push_back(it->first);
      // Cluster each available cell
      for (typename std::vector<typename Grid::CellIndex>::const_iterator it = cells.begin(); it != cells.end(); ++it)
        clusterIndex_(*it);
    }

protected:
    template <class P>
    typename Grid::cell_iterator insertCluster_(const P & p)
    {
      return grid.insert(std::make_pair(p, Cluster(p)));
    }

    void clusterIndex_(const typename Grid::CellIndex & p);

    void gridCells5x5_(typename Grid::CellIndex cur, ClusterCells & cells);

    void gridCell_(const typename Grid::CellIndex & cur, ClusterCells & cells, bool center = false, bool ignore_missing = true)
    {
      try
      {
        cells.insert(std::make_pair(cur, std::make_pair(&grid.grid_at(cur), center)));
      }
      catch (std::out_of_range &)
      {
        if (!ignore_missing) throw;
      }
    }

    void addTreeDistance_(TreeNode * tree, ClusterTrees & trees, TreeDistanceQueue & dists)
    {
      // Infinity: no valid distance
      DoubleReal dist_min = std::numeric_limits<DoubleReal>::infinity();
      typename ClusterTrees::const_iterator dist_it = trees.end();

      // Generate minimal distance to existing trees
      for (typename ClusterTrees::const_iterator it = trees.begin(); it != trees.end(); ++it)
      {
        if (tree == *it) continue;
        DoubleReal dist = treeDistance_(tree, *it);
        if (dist < dist_min)
        {
          dist_min = dist;
          dist_it = it;
        }
      }

      // Insert distance if valid one found.
      if (dist_it != trees.end()) dists.push(TreeDistance(dist_min, tree, *dist_it));

      // Insert tree.
      trees.insert(tree);
    }

    DoubleReal treeDistance_(TreeNode * left, TreeNode * right)
    {
      const BoundingBox bbox = left->bbox | right->bbox;
      if (coordElemGreater_(bbox.size(), grid.cell_dimension))
      {
        return std::numeric_limits<DoubleReal>::infinity();
      }

      const PointCoordinate left_scaled = coordElemDiv_(left->coord, grid.cell_dimension);
      const PointCoordinate right_scaled = coordElemDiv_(right->coord, grid.cell_dimension);
      return coordDist_(left_scaled, right_scaled);
    }

    void tree2Cluster_(const TreeNode * tree, Cluster & cluster)
    {
      if (tree->left && tree->right)
      {
        tree2Cluster_(tree->left, cluster);
        tree2Cluster_(tree->right, cluster);
      }
      else
      {
        cluster.insert(std::make_pair(tree->bbox.first, tree->ref));
      }
      delete tree->left;
      delete tree->right;
    }

    void tree2Points_(const TreeNode * tree)
    {
      if (tree->left && tree->right)
      {
        tree2Points_(tree->left);
        tree2Points_(tree->right);
      }
      else
      {
        insertPoint(tree->bbox.first, tree->ref);
      }
      delete tree->left;
      delete tree->right;
    }

    static PointCoordinate coordScalarDiv_(const PointCoordinate & lhs, const DoubleReal & rhs)
    {
      PointCoordinate ret;
      typename PointCoordinate::iterator it = ret.begin();
      typename PointCoordinate::const_iterator lit = lhs.begin();
      for (; it != ret.end(); ++it, ++lit) *it = *lit / rhs;
      return ret;
    }

    static PointCoordinate coordElemDiv_(const PointCoordinate & lhs, const PointCoordinate & rhs)
    {
      PointCoordinate ret;
      typename PointCoordinate::iterator it = ret.begin();
      typename PointCoordinate::const_iterator lit = lhs.begin(), rit = rhs.begin();
      for (; it != ret.end(); ++it, ++lit, ++rit) *it = *lit / *rit;
      return ret;
    }

    static bool coordElemGreater_(const PointCoordinate & lhs, const PointCoordinate & rhs)
    {
      typename PointCoordinate::const_iterator lit = lhs.begin(), rit = rhs.begin();
      for (; lit != lhs.end(); ++lit, ++rit)
      {
        if (*lit > *rit) return true;
      }
      return false;
    }

    static DoubleReal coordDist_(const PointCoordinate & lhs, const PointCoordinate & rhs)
    {
      DoubleReal ret = 0;
      PointCoordinate p = lhs - rhs;
      typename PointCoordinate::const_iterator it = p.begin();
      for (; it != p.end(); ++it) ret += std::pow(*it, 2.);
      return std::sqrt(ret);
    }

  };

  template <typename I>
  void HierarchicalClustering<I>::clusterIndex_(const typename Grid::CellIndex & cur)
  {
    ClusterCells cells;
    ClusterTrees trees;
    TreeDistanceQueue dists;

    // Collect all cells we need
    try
    {
      gridCells5x5_(cur, cells);
    }
    catch (std::out_of_range &)
    {
      return;
    }

    // Collect and remove existing points from cells
    for (typename ClusterCells::iterator cell_it = cells.begin(); cell_it != cells.end(); ++cell_it)
    {
      typename Grid::CellContent & cell_cur = *cell_it->second.first;
      const bool & cell_center = cell_it->second.second;

      // Iterate per cluster
      typename Grid::cell_iterator cluster_tmp_it = cell_cur.begin();
      while (cluster_tmp_it != cell_cur.end())
      {
        typename Grid::cell_iterator cluster_it = cluster_tmp_it;
        ++cluster_tmp_it;

        // Check if it is not yet a cluster, aka have only one point
        if (cluster_it->second.size() == 1)
        {
          // Add each point to hash grid
          for (typename Cluster::const_iterator point_it = cluster_it->second.begin(); point_it != cluster_it->second.end(); ++point_it)
          {
            const PointCoordinate & coord = point_it->first;
            TreeNode * tree(new TreeNode(coord, point_it->second, cell_center));
            addTreeDistance_(tree, trees, dists);
          }

          // Remove point from hash grid cell
          cell_cur.erase(cluster_it);
        }
      }
    }

    // Try to join two subsets with minimum distance
    while (!dists.empty())
    {
      const typename TreeDistanceQueue::value_type cur_dist = dists.top();
      TreeNode * tree_left(cur_dist.left), *tree_right(cur_dist.right);
      dists.pop();

      // Check if both trees are not yet used with a smaller distance
      Size count_left = trees.count(tree_left), count_right = trees.count(tree_right);
      if (count_left && count_right)
      {
        trees.erase(tree_left);
        trees.erase(tree_right);

        const BoundingBox bbox = tree_left->bbox | tree_right->bbox;

        // Arithmethic mean: (left * left.points + right * right.points) / (left.points + right.points)
        const PointCoordinate & left = tree_left->coord, & right = tree_right->coord;
        const UInt & left_points = tree_left->points, & right_points = tree_right->points;
        const PointCoordinate coord = coordScalarDiv_(left * left_points + right * right_points, left_points + right_points);

        TreeNode * tree(new TreeNode(coord, bbox, tree_left, tree_right));

        addTreeDistance_(tree, trees, dists);
      }
      // Re-add a distance for the tree not yet used.
      // Otherwise this subset is lost even if it is not yet maximal.
      else if (count_left)
        addTreeDistance_(tree_left, trees, dists);
      else if (count_right)
        addTreeDistance_(tree_right, trees, dists);
    }

    // Add data back to grid
    for (typename ClusterTrees::iterator tree_it = trees.begin(); tree_it != trees.end(); ++tree_it)
    {
      // We got a finished tree with all points in the center, add cluster at centroid
      if ((**tree_it).center)
      {
        Cluster & cluster = insertCluster_((**tree_it).bbox)->second;
        tree2Cluster_(*tree_it, cluster);
      }
      // We got a finished tree but not all points in the center, readd as single points
      else
      {
        tree2Points_(*tree_it);
      }
      delete *tree_it;
    }
  }

  template <typename I>
  void HierarchicalClustering<I>::gridCells5x5_(typename Grid::CellIndex base, ClusterCells & cells)
  {
    // (0, 0)
    gridCell_(base, cells, true, false);

    typename Grid::CellIndex cur = base;
    cur[0] -= 2;
    // (-2, -2)
    cur[1] -= 2; gridCell_(cur, cells);
    // (-2, -1)
    cur[1] += 1; gridCell_(cur, cells);
    // (-2, 0)
    cur[1] += 1; gridCell_(cur, cells);
    // (-2, 1)
    cur[1] += 1; gridCell_(cur, cells);
    // (-2, 2)
    cur[1] += 1; gridCell_(cur, cells);

    cur = base; cur[0] -= 1;
    // (-1, -2)
    cur[1] -= 2; gridCell_(cur, cells);
    // (-1, -1)
    cur[1] += 1; gridCell_(cur, cells, true);
    // (-1, 0)
    cur[1] += 1; gridCell_(cur, cells, true);
    // (-1, 1)
    cur[1] += 1; gridCell_(cur, cells, true);
    // (-1, 2)
    cur[1] += 1; gridCell_(cur, cells);

    cur = base;
    // (0, -2)
    cur[1] -= 2; gridCell_(cur, cells);
    // (0, -1)
    cur[1] += 1; gridCell_(cur, cells, true);
    // (0, 0)
    cur[1] += 1;
    // (0, 1)
    cur[1] += 1; gridCell_(cur, cells, true);
    // (0, 2)
    cur[1] += 1; gridCell_(cur, cells);

    cur = base; cur[0] += 1;
    // (1, -2)
    cur[1] -= 2; gridCell_(cur, cells);
    // (1, -1)
    cur[1] += 1; gridCell_(cur, cells, true);
    // (1, 0)
    cur[1] += 1; gridCell_(cur, cells, true);
    // (1, 1)
    cur[1] += 1; gridCell_(cur, cells, true);
    // (1, 2)
    cur[1] += 1; gridCell_(cur, cells);

    cur = base; cur[0] += 2;
    // (2, -2)
    cur[1] -= 2; gridCell_(cur, cells);
    // (2, -1)
    cur[1] += 1; gridCell_(cur, cells);
    // (2, 0)
    cur[1] += 1; gridCell_(cur, cells);
    // (2, 1)
    cur[1] += 1; gridCell_(cur, cells);
    // (2, 2)
    cur[1] += 1; gridCell_(cur, cells);
  }

}

#endif /* OPENMS_COMPARISON_CLUSTERING_HIERARCHICALCLUSTERING_H */