Algorithm for finding Median: Got one or can you dissect this?

[Paul E. Johnson]

I was using an old algorithm to find the median of an array of doubles. 
It was a bit slow and I saw in the Numerical Recipes 2ed that they had
some new ways. The fastest is called "select(k,N,arr[])".  They have
issued patches to fix big mistakes and I've found them and installed
them. 

Still, I get some systematic mistakes.  I enclose code below that can
compile a test program. This compares the findings of qsort and the
select() code in a function I call findMedian.  Sometimes the results
agree on the first run, but hardly ever after that. Note my recent
mastery of xmalloc and free!

If you see what's wrong, tell me.

If you have a nice, quick algorithm to find medians, please give it to
me, because I've tried pretty hard with this one.

The qsort gets awful slow when there are hundreds or thousands of
elements in the array, otherwise I would just stick with that.


Save the below in "main.m".
If you are using RedHat linux, you should be setup to compile under
emacs. Emacs (or xemacs) will grab the environment for the compiler from
the bottom of the file, a cool feature Marcus showed me last week. If
you don't use emacs, you can use the example compiler command to tailor
for your system.  

#import <simtools.h>
#import <objectbase/SwarmObject.h>
#import <misc.h>
#import <random.h>

#define NOFSPACES 3
unsigned Dimensionality[] = {3,4,5};

  /*Numerical Recipes in C, p. 341*/
#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;

int
compare_doubles (const double *a, const double *b)
{
  //return (int) (*a - *b);
  return  *a  < *b ? -1 : *a != *b;
}

double
findMedian( int k, int n, double testarr[])
{
  unsigned long i, ir, j, l, mid; 
  float a,temp;
  double arr[n];
  l=1;
  ir=n;
  

  if(n==0)return 999;

  /*This creates a copy of the array for use in the Numerical Recipes
    algorithm. */

  for(i=0;i<n;i++){printf("IN FindMedian array is index: %d %f
\n",i,testarr[i]);arr[i]=testarr[i];}

  /*Now do qsort to find the median, the k-1th element of the array*/
  qsort(testarr,n,sizeof(double),compare_doubles);
  for(i=0;i<n;i++)printf("IN FindMedian sorted array is index: %d %f
\n",i,testarr[i]);
  printf("qsort says the %d'th array element (index k-1=%d) is: %f 
\n",k,k-1, testarr[k-1]);
  /*End of Test Array */


  /*This code is from Numerical Recipes for finding the k'th ranked item
in an array.
    It has a strange effect if called repeatedly. */
  for (;;) {
    if (ir <= l+1){
      if (ir == l+1 && arr[ir] < arr[l]){
        SWAP(arr[l],arr[ir])
          }
      printf("in findmedian,  NR's algo reaches answer. The %d'th ranked
is  %f \n", k, arr[k]);
     
      return arr[k];
    }
    else
      { 
        mid=(l+ir) >> 1;
        SWAP(arr[mid],arr[l+1])
          if (arr[l] > arr[ir]){
            SWAP(arr[l],arr[ir])
              }
        if (arr[l+1] > arr[ir]) {
          SWAP(arr[l+1],arr[ir])
            }
        if (arr[l] > arr[l+1]) {
          SWAP(arr[l],arr[l+1])
            }
        i=l+1;
        j=ir;
        a=arr[l+1];
        for (;;){
          do i++; while (arr[i] < a);
          do j--; while (arr[j] > a);
          if (j < i) break;
          SWAP(arr[i],arr[j])
            }
        arr[l+1]=arr[j];
        arr[j]=a;
        if (j >= k) ir=j-1;
        if (j <= k) l=i;
      }
  }
}


int
main (int argc, const char **argv)
{
  unsigned i, s, k;
  id world;
  int numofcases;
  double * someArray;

  initSwarm (argc, argv);

   for(s=0; s< 11; s++)
     {
      numofcases=[ uniformIntRand getIntegerWithMin: 5 withMax: 25];  
       // numofcases= 10*s;
      someArray= xmalloc(numofcases * sizeof(double));

      for(i=0; i< numofcases;i++)
      someArray[i]=[uniformDblRand getDoubleWithMin: 0.0 withMax:
100.0];

      if(numofcases%2 == 0) k=(numofcases/2);
      else k=((numofcases+1)/2); 

      findMedian(k,numofcases,someArray);

      free(someArray);  
     }
}

/*
Local Variables:
compile-command: "egcs -o pj1 main.m -g -fno-inline -L/usr/lib/swarm 
-L/usr/lib -L/usr/X11R6/lib  -I/usr/include/swarm -lanalysis -lsimtools
-lsimtoolsgui -lactivity -ltkobjc -lrandom -lobjectbase  -ldefobj 
-lspace -lanalysis -lsimtools -lsimtoolsgui -ltkobjc -ltclobjc
-lactivity -lrandom -lobjectbase -lcollections -ldefobj -lmisc -lBLT
-ltk8.0 -ltcl8.0 -lXpm -lpng -lz -lffi -lX11 -lm -lobjc -lpthread -ldl
-Wl,--rpath -Wl,/usr/lib/swarm -Wl,--rpath -Wl,/usr/lib -Wl,--rpath
-Wl,/usr/X11R6/lib -Wl,--rpath -Wl,/usr/lib -Wl,--rpath -Wl,/usr/lib  "
End:
*/

-- 
Paul E. Johnson                         email: address@hidden
Dept. of Political Science              http://lark.cc.ukans.edu/~pauljohn
University of Kansas                    Office: (785) 864-9086
Lawrence, Kansas 66045                  FAX: (785) 864-5700

                  ==================================
   Swarm-Support is for discussion of the technical details of the day
   to day usage of Swarm.  For list administration needs (esp.
   [un]subscribing), please send a message to <address@hidden>
   with "help" in the body of the message.