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 {
 
  /**
   * 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 clusterize(JHitIterator_t  __begin,
                                   JHitIterator_t  __end,
                                   const JMatch_t& match,
                                   const int       Nmin = 1)
  {
    using namespace std;
 
    const int N = distance(__begin,__end);
 
    JHitIterator_t buffer = __begin;
    
    if (N >= Nmin) {
 
      int i, j;
      int n = N;
 
      // 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];
          }
        }
      }
 
      // 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;
  }
 
 
  /**
   * 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 clusterize(JHitIterator_t       __begin,
                                   JHitIterator_t       __end,
                                   const JComparator_t& comparator,
                                   const JMatch_t&      match)
  {
    if (__begin != __end) {
 
      std::sort(__begin, __end, comparator);
 
      JHitIterator_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 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 reverse_clusterize(JHitIterator_t  __begin,
                                           JHitIterator_t  __end,
                                           const JMatch_t& match)
  {
    using namespace std;
 
    const int N = distance(__begin,__end);
 
    JHitIterator_t buffer = __begin;
 
    if (N != 0) {
 
      int i, j;
 
      // 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 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;
  }
 
 
  /**
   * 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 clusterizeWeight(JHitIterator_t  __begin,
                                         JHitIterator_t  __end,
                                         const JMatch_t& match)
  {
    using namespace std;
 
    const int N = distance(__begin,__end);
 
    JHitIterator_t buffer = __begin;
 
    if (N != 0) {
 
      int i, j;
 
      // Determine weight of each hit.
 
      std::vector<double> weight(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;
  }
 
 
  /**
   * 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