JAlgorithm.hh

Go to the documentation of this file.

#ifndef __JTRIGGER__JALGORITHM__
#define __JTRIGGER__JALGORITHM__

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

#include "JTrigger/JMatch.hh"


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

namespace JTRIGGER {


  /**
   * 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 template argument <tt>JHit_t</tt> refers to a data structure which should have the following member methods:
   *   - double %getX();       // [m]
   *   - double %getY();       // [m]
   *   - double %getZ();       // [m]
   *   - double %getT();       // [ns]
   *
   * \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 JHit_t> 
  inline JHitIterator_t clusterize(JHitIterator_t        __begin, 
                                   JHitIterator_t        __end, 
                                   const JMatch<JHit_t>& match,
                                   const int             Nmin = 1)
  {
    return clusterize(__begin, __end, match, Nmin, typename std::iterator_traits<JHitIterator_t>::iterator_category());
  }


  /**
   * Implementation of method clusterize for random access iterators.
   */
  template<class JHitIterator_t, class JHit_t> 
  inline JHitIterator_t clusterize(JHitIterator_t        __begin, 
                                   JHitIterator_t        __end, 
                                   const JMatch<JHit_t>& match,
                                   const int             Nmin,
                                   std::random_access_iterator_tag)
  {
    const int N = std::distance(__begin, __end);

    if (N >= Nmin) {

      int i, j;
      int n = N;

      JHitIterator_t buffer = __begin;

      // Determine number of associated hits for each hit.
    
      std::vector<int> count(N, 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];

          } else if (n == Nmin) {

            return __begin;
          }
        }

        if (count[i] < Nmin) {

          --n;

          if (n < Nmin) {
            return __begin;
          }
        }
      }

      n = N;

      // Remove hit with the smallest number of associated hits.
      // This procedure stops when the number of associated hits 
      // is equal to the number of (remaining) hits.
    
      for ( ; n >= Nmin; ) {
      
        j = 0;
      
        for (i = 1; i != n; ++i) {
          if (count[i] < count[j]) {
            j = i;
          }
        }
      
        // Ready?

        if (count[j] == n) {
          break;
        }
      
        // Reduce effective size.

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

        std::swap(buffer[j], buffer[n]); 
        std::swap(count [j], count [n]);
      
        // Decrease number of associated hits for each associated hit.

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

      if (n >= Nmin)
        return __begin + n;
      else
        return __begin;
    }
        
    return __begin;
  }


  /**
   * 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 template argument <tt>JHit_t</tt> refers to a data structure which should have the following member methods:
   *   - double %getX();       // [m]
   *   - double %getY();       // [m]
   *   - double %getZ();       // [m]
   *   - double %getT();       // [ns]
   *
   * \param  __begin        begin of data
   * \param  __end          end   of data
   * \param  match          binary match operator
   * \return                end   of data
   */
  template<class JHitIterator_t, class JHit_t> 
  inline JHitIterator_t reverse_clusterize(JHitIterator_t        __begin, 
                                           JHitIterator_t        __end, 
                                           const JMatch<JHit_t>& match)
  {
    return reverse_clusterize(__begin, __end, match, typename std::iterator_traits<JHitIterator_t>::iterator_category());
  }


  /**
   * Implementation of method reverse_clusterize for random access iterators.
   */
  template<class JHitIterator_t, class JHit_t> 
  inline JHitIterator_t reverse_clusterize(JHitIterator_t        __begin, 
                                           JHitIterator_t        __end, 
                                           const JMatch<JHit_t>& match,
                                           std::random_access_iterator_tag)
  {
    const int N = std::distance(__begin, __end);

    if (N != 0) {

      int i, j;

      JHitIterator_t buffer = __begin;

      // Determine number of associated hits for each hit.
    
      std::vector<int> count(N, 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 k = 0; ; ) {

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

        if (count[j] == 1) {
          return __begin + k;
        }

        // Swap the selected hit to begin.

        std::swap(buffer[j], buffer[k]);
        std::swap(count [j], count [k]);
      
        // Decrease number of associated hits for each associated hit.

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

        // Increase effective size.

        ++k;
      }
    }

    return __begin;
  }


  /**
   * 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 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. 
   *
   * \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 JHit_t>
  inline JHitIterator_t clusterize(JHitIterator_t        __begin,
                                   JHitIterator_t        __end,
                                   JComparator_t         comparator,
                                   const JMatch<JHit_t>& match)
  {
    if (__begin != __end) {

      std::sort(__begin, __end, comparator);
      
      const JHit_t& root = *__begin;

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

    } else {

      return __end;
    }
  }


  /**
   * 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 template argument <tt>JHit_t</tt> refers to a data structure which should have the following member methods:
   *   - double %getX();       // [m]
   *   - double %getY();       // [m]
   *   - double %getZ();       // [m]
   *   - double %getT();       // [ns]
   *   - double %getW();       // [au]
   *
   * \param  __begin        begin of data
   * \param  __end          end   of data
   * \param  match          binary match operator
   * \return                end   of data
   */
  template<class JHitIterator_t, class JHit_t>
  inline JHitIterator_t clusterizeWeight(JHitIterator_t        __begin,
                                         JHitIterator_t        __end,
                                         const JMatch<JHit_t>& match)
  {
    return clusterizeWeight(__begin, __end, match, typename std::iterator_traits<JHitIterator_t>::iterator_category());
  }


  /**
   * Implementation of method clusterizeWeight for random access iterators.
   */
  template<class JHitIterator_t, class JHit_t>
  inline JHitIterator_t clusterizeWeight(JHitIterator_t        __begin,
                                         JHitIterator_t        __end,
                                         const JMatch<JHit_t>& match,
                                         std::random_access_iterator_tag)
  {
    if (__begin != __end) {

      const int N = std::distance(__begin, __end);

      int i, j;

      // Determine weight of each hit.

      std::vector<double> weight(N, 0.0);
    
      JHitIterator_t buffer = __begin;

      for (i = 0; i != N; ++i) {
      
        weight[i] += buffer[i].getW();
      
        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;
          }

          if (weight[i] > W) {
            W = weight[i];
          }
        }
      
        // Ready?
      
        if (weight[j] == W) {
          break;
        }
      
        // Reduce effective size.
      
        --n;
      
        // Swap the selected hit to end.
      
        std::swap(buffer[j], buffer[n]);
        std::swap(weight[j], weight[n]);
      
        // Decrease weight of associated hits for each associated hit.
      
        for (i = 0; i != n; ++i) {
          if (match(buffer[i], buffer[n])) {
            weight[i] -= buffer[n].getW();
          }
        }
      }

      return __begin + n;

    } else {

      return __end;
    }
  }


  /**
   * 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 %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