sort.h

Go to the documentation of this file.

#ifndef SORT_H
#define SORT_H

#include <itpp/base/vec.h>
#include <itpp/base/converters.h>
#include <itpp/base/math/log_exp.h>


namespace itpp
{

enum SORTING_METHOD { INTROSORT = 0, QUICKSORT = 1, HEAPSORT = 2,
                      INSERTSORT = 3
                    };

template<class T>
class Sort
{
public:
  Sort(SORTING_METHOD method = INTROSORT): sort_method(method) {}

  void set_method(SORTING_METHOD method) { sort_method = method; }

  SORTING_METHOD get_method() const { return sort_method; }

  void sort(int low, int high, Vec<T> &data);

  ivec sort_index(int low, int high, const Vec<T> &data);

  void intro_sort(int low, int high, int max_depth, Vec<T> &data);

  ivec intro_sort_index(int low, int high, int max_depth,
                        const Vec<T> &data);

private:
  SORTING_METHOD sort_method;

  void IntroSort(int low, int high, int max_depth, T data[]);
  void IntroSort_Index(int low, int high, int max_depth, int indexlist[],
                       const T data[]);

  void QuickSort(int low, int high, T data[]);
  void QuickSort_Index(int low, int high, int indexlist[], const T data[]);

  void HeapSort(int low, int high, T data[]);
  void HeapSort_Index(int low, int high, int indexlist[], const T data[]);

  void InsertSort(int low, int high, T data[]);
  void InsertSort_Index(int low, int high, int indexlist[], const T data[]);
};


template<class T>
void sort(Vec<T> &data, SORTING_METHOD method = INTROSORT)
{
  Sort<T> s(method);
  s.sort(0, data.size() - 1, data);
}

template<class T>
ivec sort_index(const Vec<T> &data, SORTING_METHOD method = INTROSORT)
{
  Sort<T> s(method);
  return s.sort_index(0, data.size() - 1, data);
}


// ----------------------------------------------------------------------
// Public functions for various sorting methods
// ----------------------------------------------------------------------

template<class T>
void Sort<T>::sort(int low, int high, Vec<T> &data)
{
  int N = data.size();
  // Nothing to sort if data vector has only one or zero elements
  if (N < 2)
    return;

  it_assert((low >= 0) && (high > low) && (high < N), "Sort::sort(): "
            "low or high out of bounds");

  switch (sort_method) {
  case INTROSORT:
    IntroSort(low, high, levels2bits(N), data._data());
    break;
  case QUICKSORT:
    QuickSort(low, high, data._data());
    break;
  case HEAPSORT:
    HeapSort(low, high, data._data());
    break;
  case INSERTSORT:
    InsertSort(low, high, data._data());
    break;
  default:
    it_error("Sort<T>::sort(): Unknown sorting method");
  }
}


template<class T>
ivec Sort<T>::sort_index(int low, int high, const Vec<T> &data)
{
  int N = data.size();
  // Nothing to sort if data vector has only one or zero elements
  if (N == 1)
    return ivec("0");
  else if (N == 0)
    return ivec();

  it_assert((low >= 0) && (high > low) && (high < N), "Sort::sort(): "
            "low or high out of bounds");

  ivec indexlist(N);
  for (int i = 0; i < N; ++i) {
    indexlist(i) = i;
  }

  switch (sort_method) {
  case INTROSORT:
    IntroSort_Index(low, high, levels2bits(N), indexlist._data(),
                    data._data());
    break;
  case QUICKSORT:
    QuickSort_Index(low, high, indexlist._data(), data._data());
    break;
  case HEAPSORT:
    HeapSort_Index(low, high, indexlist._data(), data._data());
    break;
  case INSERTSORT:
    InsertSort_Index(low, high, indexlist._data(), data._data());
    break;
  default:
    it_error("Sort<T>::sort_index(): Unknown sorting method");
  }

  return indexlist;
}


// INTRO SORT
template<class T>
void Sort<T>::intro_sort(int low, int high, int max_depth, Vec<T> &data)
{
  it_assert((low >= 0) && (high > low) && (high < data.size()),
            "Sort::sort(): low or high out of bounds");
  IntroSort(low, high, max_depth, data._data());
}

// INTRO SORT INDEX
template<class T>
ivec Sort<T>::intro_sort_index(int low, int high, int max_depth,
                               const Vec<T> &data)
{
  int N = data.size();
  it_assert((low >= 0) && (high > low) && (high < N),
            "Sort::sort(): low or high out of bounds");

  ivec indexlist(N);
  for (int i = 0; i < N; ++i) {
    indexlist(i) = i;
  }

  IntroSort_Index(low, high, max_depth, indexlist._data(), data._data());

  return indexlist;
}


// ----------------------------------------------------------------------
// Private functions for sorting methods
// ----------------------------------------------------------------------

template<class T>
void Sort<T>::IntroSort(int low, int high, int max_depth, T data[])
{
  if (high - low > 16) {
    max_depth--;
    if (max_depth == 0) {
      HeapSort(low, high, data);
      return;
    }

    if (high > low) {
      T a = data[low];
      int plow = low;
      int phigh = high;
      T test = data[phigh];
      while (plow < phigh) {
        if (test < a) {
          data[plow] = test;
          plow++;
          test = data[plow];
        }
        else {
          data[phigh] = test;
          phigh--;
          test = data[phigh];
        }
      }
      data[plow] = a;
      IntroSort(low, plow - 1, max_depth, data);
      IntroSort(plow + 1, high, max_depth, data);
      return;
    }
  }
  else {
    InsertSort(low, high, data);
    return;
  }
}

template<class T>
void Sort<T>::IntroSort_Index(int low, int high, int max_depth,
                              int indexlist[], const T data[])
{
  if (high - low > 16) {
    max_depth--;
    if (max_depth == 0) {
      HeapSort_Index(low, high, indexlist, data);
      return;
    }

    if (high > low) {
      int aindex = indexlist[low];
      T a = data[aindex];
      int plow = low;
      int phigh = high;
      int testindex = indexlist[phigh];
      T test = data[testindex];
      while (plow < phigh) {
        if (test < a) {
          indexlist[plow] = testindex;
          plow++;
          testindex = indexlist[plow];
          test = data[testindex];
        }
        else {
          indexlist[phigh] = testindex;
          phigh--;
          testindex = indexlist[phigh];
          test = data[testindex];
        }
      }
      indexlist[plow] = aindex;
      IntroSort_Index(low, plow - 1, max_depth, indexlist, data);
      IntroSort_Index(plow + 1, high, max_depth, indexlist, data);
    }
  }
  else {
    InsertSort_Index(low, high, indexlist, data);
    return;
  }
}

template <class T>
void Sort<T>::QuickSort(int low, int high, T data[])
{
  if (high > low) {
    T a = data[low];
    int plow = low;
    int phigh = high;
    T test = data[phigh];
    while (plow < phigh) {
      if (test < a) {
        data[plow] = test;
        plow++;
        test = data[plow];
      }
      else {
        data[phigh] = test;
        phigh--;
        test = data[phigh];
      }
    }
    data[plow] = a;
    QuickSort(low, plow - 1, data);
    QuickSort(plow + 1, high, data);
  }
}

template<class T>
void Sort<T>::QuickSort_Index(int low, int high, int indexlist[],
                              const T data[])
{
  if (high > low) {
    int aindex = indexlist[low];
    T a = data[aindex];
    int plow = low;
    int phigh = high;
    int testindex = indexlist[phigh];
    T test = data[testindex];
    while (plow < phigh) {
      if (test < a) {
        indexlist[plow] = testindex;
        plow++;
        testindex = indexlist[plow];
        test = data[testindex];
      }
      else {
        indexlist[phigh] = testindex;
        phigh--;
        testindex = indexlist[phigh];
        test = data[testindex];
      }
    }
    indexlist[plow] = aindex;
    QuickSort_Index(low, plow - 1, indexlist, data);
    QuickSort_Index(plow + 1, high, indexlist, data);
  }
}

template<class T>
void Sort<T>::HeapSort(int low, int high, T data[])
{
  int size = (high + 1) - low;
  int i = size / 2;
  T temp;
  while (1) {
    if (i > 0)
      temp = data[--i + low];
    else {
      if (size-- == 0)
        break;
      temp = data[size + low];
      data[size + low] = data[low];
    }

    int parent = i;
    int child = i * 2 + 1;

    while (child < size) {
      if (child + 1 < size  &&  data[child + 1 + low] > data[child + low])
        child++;
      if (data[child + low] > temp) {
        data[parent + low] = data[child + low];
        parent = child;
        child = parent * 2 + 1;
      }
      else
        break;
    }
    data[parent + low] = temp;
  }
}

template<class T>
void Sort<T>::HeapSort_Index(int low, int high, int indexlist[],
                             const T data[])
{
  int size = (high + 1) - low;
  int i = size / 2;

  while (1) {
    T tempValue;
    int tempIndex;
    if (i > 0) {
      i--;
      tempValue = data[indexlist[i + low]];
      tempIndex = indexlist[i + low];
    }
    else {
      if (size-- == 0)
        break;
      tempValue = data[indexlist[size + low]];
      tempIndex = indexlist[size + low];
      indexlist[size+low] = indexlist[low];
    }

    int parent = i;
    int child = i * 2 + 1;

    while (child < size) {
      if ((child + 1 < size)
          && data[indexlist[child + 1 + low]] > data[indexlist[child + low]])
        child++;
      if (data[indexlist[child + low]] > tempValue) {
        indexlist[parent + low] = indexlist[child + low];
        parent = child;
        child = parent * 2 + 1;
      }
      else
        break;
    }
    indexlist[parent + low] = tempIndex;
  }
}

template<class T>
void Sort<T>::InsertSort(int low, int high, T data[])
{
  for (int i = low + 1; i <= high; i++) {
    T value = data[i];
    int j;
    for (j = i - 1; j >= low && data[j] > value; j--) {
      data[j + 1] = data[j];
    }
    data[j + 1] = value;
  }
}

template<class T>
void Sort<T>::InsertSort_Index(int low, int high, int indexlist[],
                               const T data[])
{
  for (int i = low + 1; i <= high; i++) {
    T value = data[indexlist[i]];
    int tempIndex = indexlist[i];
    int j;
    for (j = i - 1; j >= low && data[indexlist[j]] > value; j--) {
      indexlist[j + 1] = indexlist[j];
    }
    indexlist[j + 1] = tempIndex;
  }
}


} // namespace itpp

#endif // #ifndef SORT_H