#ifndef __JALGORITHM__
#define __JALGORITHM__

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

#include "JTrigger/JMatch.hh"


namespace JTRIGGER {


  /**
   * Partition data according given binary match operator.
   *
   * 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.
   * This algorithm is known in litterature as 'clique'.
   *
   * \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 = 0)
  {
    return clusterize(__begin, __end, match, Nmin, typename std::iterator_traits<JHitIterator_t>::iterator_category());
  }


  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;

      // Count number of associated hits for each hit.
    
      std::vector<int> numberOfAssociatedHits(N, 1);      // Assume always a match with itself.
      
      for (i = 0; i != N; ++i) {

	for (j = i; ++j != N; ) {

	  if (match(*(__begin + i),
		    *(__begin + j))) {

	    ++numberOfAssociatedHits[i];
	    ++numberOfAssociatedHits[j];

	  } else if (n == Nmin)

	    return __begin;
	}

	if (numberOfAssociatedHits[i] < Nmin) {

	  --n;

	  if (n < Nmin)
	    return __begin;
	}
      }

      n = N;

      // Remove hit with the smallest number of associated hits.
      // This procdure 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 (numberOfAssociatedHits[i] < numberOfAssociatedHits[j])
	    j = i;
	}
      
	// Ready?

	if (numberOfAssociatedHits[j] == n)
	  break;
      
	// Reduce effective size.

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

	std::iter_swap(__begin + j, 
		       __begin + n);
      
	// Idem dito for the two counters.

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

	for (i = 0; i != n; ++i) {
	  if (match(*(__begin + i),
		    *(__begin + n)))
	    --numberOfAssociatedHits[i];
	}
      }

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


  /**
   * Partition data according given binary match operator.
   *
   * 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.
   * This algorithm is known in lietarture as 'clique'.
   *
   * \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());
  }


  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;

      // Count number of associated hits for each hit.
    
      std::vector<int> numberOfAssociatedHits(N, 1);      // Assume always a match with itself.

      for (i = 0; i != N; ++i) {
      
	for (j = 0; j != i; ++j) {
	  if (match(*(__begin + i),
		    *(__begin + j))) {
	    ++numberOfAssociatedHits[i];
	    ++numberOfAssociatedHits[j];
	  }
	}
      }

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

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

	if (numberOfAssociatedHits[j] == 1)
	  return __begin + k;

	// Swap the selected hit to begin.

        std::iter_swap(__begin + j, __begin + k);

	// Idem dito for the two counters.

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

	for (i = k; ++i != N; ) {
	  if (match(*(__begin + i),
		    *(__begin + k)))
	    --numberOfAssociatedHits[i];
	}

	// Increase effective size.

        ++k;
      }
    }

    return __begin;
  }
}

#endif