#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 clusterize {

    /**
     * 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 reverse_clusterize {

    /**
     * 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 clusterizeWeight {

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