JAlgorithm.hh

Go to the documentation of this file.

#ifndef __JTRIGGER__JALGORITHM__
#define __JTRIGGER__JALGORITHM__

#include <vector>
#include <algorithm>
#include <iterator>


/**
 * \file
 *
 * Algorithms for hit clustering and sorting.
 * \author mdejong
 */
namespace JTRIGGER {}
namespace JPP { using namespace JTRIGGER; }

namespace JTRIGGER {

  /**
   * Anonymous structure for clustering of hits.
   */
  static const struct {

    /**
     * Partition data according given binary match operator.
     *
     * The underlying algorithm is known in literature as 'clique'.
     * The result is (likely to be) the maximal sub-set with all elements matched to each other.
     * The complexity is quadratic, i.e.\ at most (number of elements x number of elements) operations.
     * The algorithm will sort the data such that all clusterized hits are at the front. 
     * The return value points the first non clusterized hit.
     *
     * The hit iterator refers to a data structure which should conform with the match operator.
     *
     * \param  __begin        begin of data
     * \param  __end          end   of data
     * \param  match          binary match operator
     * \param  Nmin           minimum cluster size
     * \return                end   of data
     */
    template<class JHitIterator_t, class JMatch_t>
    inline JHitIterator_t operator()(JHitIterator_t  __begin,
                                     JHitIterator_t  __end,
                                     const JMatch_t& match,
                                     const int       Nmin = 1) const
    {
      return (*this)(__begin, std::distance(__begin,__end), match, Nmin, typename std::iterator_traits<JHitIterator_t>::iterator_category());
    }


    /**
     * Select best root hit according given comparator and partition remaining data
     * according given binary match operator and this root hit.
     * The complexity is linear, i.e.\ at most number of elements operations.
     * The algorithm will sort the data such that all selected hits are at the front. 
     *
     * \param  __begin        begin of data
     * \param  __end          end   of data
     * \param  comparator     comparator
     * \param  match          binary match operator
     * \return                end   of data
     */
    template<class JHitIterator_t, class JComparator_t, class JMatch_t>
    inline JHitIterator_t operator()(JHitIterator_t       __begin,
                                     JHitIterator_t       __end,
                                     const JComparator_t& comparator,
                                     const JMatch_t&      match) const
    {
      if (__begin != __end) {

        std::sort(__begin, __end, comparator);

        JHitIterator_t root = __begin;

        return std::partition(++__begin, __end, JBind2nd(match, *root));

      } else {

        return __end;
      }
    }

  private:
    /**
     * Implementation of method clusterize for random access iterators.
     *
     * \param  buffer         pointer to data
     * \param  N              number of hits
     * \param  match          binary match operator
     * \param  Nmin           minimum cluster size
     * \param  tag            iterator tag
     * \return                pointer to end of data
     */
    template<class JHitIterator_t, class JMatch_t>
    inline JHitIterator_t operator()(JHitIterator_t  buffer,
                                     const int       N,
                                     const JMatch_t& match,
                                     const int       Nmin,
                                     std::random_access_iterator_tag tag) const
    {
      using namespace std;

      if (N >= Nmin) {

        int i, j;
        int n = N;

        // Determine number of associated hits for each hit.

        count.resize(N);

        for (vector<int>::iterator __i = count.begin(); __i != count.end(); ++__i) {
          *__i = 1;                        // Assume always a match with itself.
        }

        for (i = 0; i != N; ++i) {
          for (j = i; ++j != N; ) {
            if (match(buffer[i], buffer[j])) {
              ++count[i];
              ++count[j];
            }
          }
        }

        // Remove hit with the smallest number of associated hits.
        // This procedure stops when:
        // 1) number of associated hits is equal to the number of (remaining) hits or
        // 2) maximal number of associated hits is less than the specified minimum.

        for ( ; ; ) {
      
          int M = (count[0] >= Nmin ? 1 : 0);
      
          j = 0;

          for (i = 1; i != n; ++i) {
            if (count[i] < count[j]) {
              j = i;
            }
            if (count[i] >= Nmin) {
              ++M;
            }
          }
      
          // Ready?

          if (count[j] == n) {
            return buffer + n;
          }

          if (M < Nmin) {
            return buffer;
          }
      
          // Reduce effective size.

          --n;
      
          // Swap the selected hit to end.

          swap(buffer[j], buffer[n]); 
          swap(count [j], count [n]);

          // Decrease number of associated hits for each associated hit.

          for (i = 0; i != n && count[n] != 1; ++i) {
            if (match(buffer[i], buffer[n])) {
              --count[i];
              --count[n];
            }
          }
        }
      }

      return buffer;
    }


    mutable std::vector<int> count;

  } clusterize;


  /**
   * Anonymous structure for reverse clustering of hits.
   */
  static const struct {

    /**
     * Partition data according given binary match operator.
     *
     * The underlying algorithm is known in literature as reverse 'clique'.
     * The result is (likely to be) the maximal sub-set with all elements matched to each other.
     * The complexity is quadratic, i.e.\ at most (number of elements x number of elements) operations.
     * The algorithm will sort the data such that all clusterized hits are at the front. 
     * The return value points the first non clusterized hit.
     *
     * The hit iterator refers to a data structure which should conform with the match operator.
     *
     * \param  __begin        begin of data
     * \param  __end          end   of data
     * \param  match          binary match operator
     * \return                end   of data
     */
    template<class JHitIterator_t, class JMatch_t>
    inline JHitIterator_t operator()(JHitIterator_t  __begin,
                                     JHitIterator_t  __end,
                                     const JMatch_t& match) const
    {
      return (*this)(__begin,  std::distance(__begin,__end), match, typename std::iterator_traits<JHitIterator_t>::iterator_category());
    }

  private:
    /**
     * Implementation of method reverse_clusterize for random access iterators.
     *
     * \param  buffer         pointer to data
     * \param  N              number of hits
     * \param  match          binary match operator
     * \param  tag            iterator tag
     * \return                pointer to end of data
     */
    template<class JHitIterator_t, class JMatch_t>
    inline JHitIterator_t operator()(JHitIterator_t  buffer,
                                     const int       N,
                                     const JMatch_t& match,
                                     std::random_access_iterator_tag tag) const
    {
      using namespace std;

      if (N != 0) {

        int i, j;

        // Determine number of associated hits for each hit.
    
        count.resize(N);

        for (vector<int>::iterator __i = count.begin(); __i != count.end(); ++__i) {
          *__i = 1;                        // Assume always a match with itself.
        }

        for (i = 0; i != N; ++i) {
          for (j = 0; j != i; ++j) {
            if (match(buffer[i], buffer[j])) {
              ++count[i];
              ++count[j];
            }
          }
        }

        // Keep hit with the largest number of associated hits.
        // This procedure stops when the number of associated hits is equal to one.
 
        for (int n = 0; ; ++n) {

          j = n;
      
          for (i = j; ++i != N; ) {
            if (count[i] > count[j]) {
              j = i;
            }
          }
      
          // Ready?

          if (count[j] == 1) {
            return buffer + n;
          }

          // Swap the selected hit to begin.

          swap(buffer[j], buffer[n]);
          swap(count [j], count [n]);

          // Decrease number of associated hits for each associated hit.

          for (i = n; ++i != N; ) {
            if (match(buffer[i], buffer[n])) {
              --count[i];
              --count[n];
            }
          }
        }
      }

      return buffer;
    }


    mutable std::vector<int> count;

  } reverse_clusterize;


  /**
   * Anonymous struct for weighed clustering of hits.
   */
  static const struct {

    /**
     * Partition data according given binary match operator.
     *
     * The underlying algorithm is known in literature as reverse 'clique'.
     * The result is (likely to be) the maximal sub-set with all elements matched to each other.
     * The complexity is quadratic, i.e.\ at most (number of elements x number of elements) operations.
     * The algorithm will sort the data such that all clusterized hits are at the front. 
     * The return value points the first non clusterized hit.
     * Note that the weight is assumed to be positive definite and the larger the weight the better.
     *
     * The hit iterator refers to a data structure which should conform with the match operator.
     * In addition, it should have the following member method:
     *   - double %getW();       // [a.u.]
     *
     * \param  __begin        begin of data
     * \param  __end          end   of data
     * \param  match          binary match operator
     * \return                end   of data
     */
    template<class JHitIterator_t, class JMatch_t>
    inline JHitIterator_t operator()(JHitIterator_t  __begin,
                                     JHitIterator_t  __end,
                                     const JMatch_t& match) const
    {
      return (*this)(__begin,  std::distance(__begin,__end), match, typename std::iterator_traits<JHitIterator_t>::iterator_category());
    }

  private:
    /**
     * Implementation of method clusterizeWeight for random access iterators.
     *
     * \param  buffer         pointer to data
     * \param  N              number of hits
     * \param  match          binary match operator
     * \param  tag            iterator tag
     * \return                pointer to end of data
     */
    template<class JHitIterator_t, class JMatch_t>
    inline JHitIterator_t operator()(JHitIterator_t  buffer,
                                     const int       N,
                                     const JMatch_t& match,
                                     std::random_access_iterator_tag tag) const
    {
      using namespace std;

      if (N != 0) {

        int i, j;

        // Determine weight of each hit.

        weight.resize(N);

        for (i = 0; i != N; ++i) {      
          weight[i] = buffer[i].getW();    // Assume always a match with itself.
        }

        for (i = 0; i != N; ++i) {
          for (j = i; ++j != N; ) {
            if (match(buffer[i], buffer[j])) {
              weight[i] += buffer[j].getW();
              weight[j] += buffer[i].getW();
            }
          }
        }

        // Remove hit with the smallest weight of associated hits.
        // This procedure stops when the weight of associated hits
        // is equal to the maximal weight of (remaining) hits.
    
        int    n = N;
        double W;
    
        for ( ; ; ) {
      
          j = 0;
          W = weight[j];

          for (i = 1; i != n; ++i) {
            if      (weight[i] < weight[j])
              j = i;
            else if (weight[i] > W)
              W = weight[i];
          }
      
          // Ready?
      
          if (weight[j] == W) {
            return buffer + n;
          }
      
          // Reduce effective size.
      
          --n;
      
          // Swap the selected hit to end.

          swap(buffer[j], buffer[n]);
          swap(weight[j], weight[n]);

          // Decrease weight of associated hits for each associated hit.

          for (i = 0; i != n && weight[n] > buffer[n].getW(); ++i) {
            if (match(buffer[i], buffer[n])) {
              weight[i] -= buffer[n].getW();
              weight[n] -= buffer[i].getW();
            }
          }
        }
      }

      return buffer;
    }


    mutable std::vector<double> weight;

  } clusterizeWeight;


  /**
   * Compare two hits by weight.
   *
   * The template argument <tt>JHit_t</tt> refers to a data structure which should have the following member method:
   *   - double %getW();       // [a.u.]
   *
   * \param  first      first  hit
   * \param  second     second hit
   * \return            true if first hit has larger weight than second; else false
   */
  template<class JHit_t>
  inline bool weightSorter(const JHit_t& first, const JHit_t& second)
  {
    return first.getW() > second.getW();
  }


  /**
   * Compare two hits by weight.
   *
   * The template argument <tt>JHit_t</tt> refers to a data structure which should have the following member methods:
   *   - double %getT();       // [a.u.]
   *
   * \param  first      first  hit
   * \param  second     second hit
   * \return            true if first hit earlier than second; else false
   */
  template<class JHit_t>
  inline bool timeSorter(const JHit_t& first, const JHit_t& second)
  {
    return first.getT() < second.getT();
  }
}

#endif