Distributed Information System (DIS)
  • Home
  • The blog
  • Contact

Median value selection (Fixed)

12/20/2017

0 Comments

 
In 2009 I presented a heap based median selection algorithm. It was original, and was apparently very fast when compiled with the Intel compiler (icc). Since I don't have the Intel compiler anymore, I can't test it's performance. It's slower than the nthElement code given below when compiled with g++-7 -O3.

Here is the fixed code.

float fixedHeapMedian (float *a) {
  const unsigned char HEAP_LEN = 13;
  float left[HEAP_LEN], right[HEAP_LEN], *p, median;
  unsigned char nLeft, nRight;

  // pick first value as median candidate
  p = a;
  median = *p++;
  nLeft = nRight = 0;

  for (;;) {
    //dumpState(left, nLeft, median, right, nRight, p, 27 - (p-a));
    //assert(stateIsValid(left, nLeft, median, right, nRight));

    // get next value
    float val = *p++;

    // if value is smaller than median, append to left heap
    if (val <= median) {
      // move biggest value to the top of left heap
      unsigned char child = nLeft++, parent = (child - 1) / 2;
      while (child && val > left[parent]) {
        left[child] = left[parent];
        child = parent;
        parent = (parent - 1) / 2;
      }
      left[child] = val;

      // if left heap is full
      if (nLeft == HEAP_LEN) {
        //cout << "---" << endl;
        // for each remaining value
        for (unsigned char nVal = 27-(p - a); nVal; --nVal) {
          //dumpState(left, nLeft, median, right, nRight, p, nVal);
          //assert(stateIsValid(left, nLeft, median, right, nRight));
          // get next value
          val = *p++;
          // discard values falling in other heap
          if (val >= median) {
            continue;
          }
          // if val is bigger than biggest in heap, val is new median
          if (val >= left[0]) {
            median = val;
            continue;
          }
          // biggest heap value becomes new median
          median = left[0];
          // insert val in heap
          parent = 0;
          child = 2;
          while (child < HEAP_LEN) {
            if (left[child-1] > left[child]) {
              child = child-1;
            }
            if (val >= left[child]) {
               break;
            }
            left[parent] = left[child];
            parent = child;
            child = (parent + 1) * 2;
          }
          left[parent] = val;
        }
        return median;
      }
    } else {
      // move smallest value to the top of right heap
      unsigned char child = nRight++, parent = (child - 1) / 2;
      while (child && val < right[parent]) {
        right[child] = right[parent];
        child = parent;
        parent = (parent - 1) / 2;
      }
      right[child] = val;

      // if right heap is full
      if (nRight == HEAP_LEN) {
        //cout << "---" << endl;
        // for each remaining value
        for (unsigned char nVal = 27-(p - a); nVal; --nVal) {
          //dumpState(left, nLeft, median, right, nRight, p, nVal);
          //assert(stateIsValid(left, nLeft, median, right, nRight));
          // get next value
          val = *p++;
          // discard values falling in other heap
          if (val <= median) {
            continue;
          }
          // if val is smaller than smallest in heap, val is new median
          if (val <= right[0]) {
            median = val;
            continue;
          }
          // heap top value becomes new median
          median = right[0];
          // insert val in heap
          parent = 0;
          child = 2;
          while (child < HEAP_LEN) {
            if (right[child-1] < right[child]) {
              child = child-1;
            }
            if (val <= right[child]) {
              break;
            }
            right[parent] = right[child];
            parent = child;
            child = (parent + 1) * 2;
          }
          right[parent] = val;
        }
        return median;
      }
    }
  }
}
0 Comments



Leave a Reply.

    Author

    Christophe Meessen is a  computer science engineer working in France.

    Any suggestions to make DIS more useful ? Tell me by using the contact page.

    Categories

    All
    Business Model
    Database
    Dis
    Ditp
    Dvcs
    Git
    Gob
    Idr
    Misc
    Murphys Law
    Programming Language
    Progress Status
    Startup
    Suggested Reading
    Web Site

    Archives

    December 2017
    November 2015
    September 2015
    February 2013
    December 2012
    November 2012
    May 2012
    February 2012
    March 2010
    October 2009
    September 2009
    July 2009
    June 2009
    May 2009
    February 2009
    January 2009
    November 2008
    September 2008
    August 2008
    July 2008
    May 2008
    April 2008
    March 2008
    February 2008
    January 2008
    December 2007
    October 2007
    August 2007
    July 2007
    June 2007
    May 2007

    RSS Feed

    Live traffic feed
    You have no departures or arrivals yet. Wait a few minutes and check again.
    Powered by FEEDJIT
Powered by Create your own unique website with customizable templates.