# [ROOT] Problem with Median()

From: cstrato (cstrato@aon.at)
Date: Fri Sep 24 2004 - 18:26:11 MEST

```Dear Rooters

The attached macro "macroMedian.C" contains three implementations
of "median" with the following problems:
- Why does implementation Median (CRAN) only run when compiled
with ACLiC but not with CINT?
- Why is the output of TMath::Median() different?

root [0] .x macroMedian.C
Median (ROOT) = 0.509148
Median (C.S.) = 0.508224

*** Break *** illegal instruction

root [0] .x macroMedian.C+
Info in <TUnixSystem::ACLiC>: creating shared library
/Users/rabbitus/ROOT/rootcode/xpsx/./macroMedian_C.so
Median (ROOT) = 0.496104
Median (C.S.) = 0.495683
Median (CRAN) = 0.495683
root [1] .x macroMedian.C+
Median (ROOT) = 0.509148
Median (C.S.) = 0.508224
Median (CRAN) = 0.508224

Best regards
Christian
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
C.h.r.i.s.t.i.a.n. .S.t.r.a.t.o.w.a
V.i.e.n.n.a.         .A.u.s.t.r.i.a
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-

/******************************************************************************
* macro to test median                                                        *
*                                                                             *
* Author: Dr. Christian Stratowa, Vienna, Austria.                            *
******************************************************************************/

#include "TMath.h"
#include "TRandom.h"

#include <math.h>
#include <Riostream.h>

/**********************************************************
** Source code from BioConductor package "affy"
**
** int sort_double(const void *a1,const void *a2)
**
** a comparison function used when sorting doubles.
**
**********************************************************/

int sort_double(const double *a1,const double *a2){
if (*a1 < *a2)
return (-1);
if (*a1 > *a2)
return (1);
return 0;
}

/**************************************************************************
** Source code from BioConductor package "affy"
**
** double median(double *x, int length)
**
** double *x - vector
** int length - length of *x
**
** returns the median of *x
**
*************************************************************************/

double  median(double *x, int length)
{
int i;
int half;
double med;
double *buffer = (double *)calloc(length,sizeof(double));

for (i = 0; i < length; i++)
buffer[i] = x[i];

qsort(buffer,length,sizeof(double), (int(*)(const void*, const void*))sort_double);
half = (length + 1)/2;
if (length % 2 == 1){
med = buffer[half - 1];
} else {
med = (buffer[half] + buffer[half-1])/2.0;
}

free(buffer);
return med;
}

//______________________________________________________________________________
Double_t XMedian(Int_t n, const Double_t *arr)
{
// Calculate median (source code C. Stratowa)

if (n == 1) return arr[0];

// Create index and sort array
Int_t *index = 0;
if (!(index = new Int_t[n])) {return 1;} //return not ok
TMath::Sort(n, arr, index);

// Find median
Int_t k;
Double_t median = 0;
if ((n % 2) == 0){
k = (Int_t)TMath::Floor(n / 2.0) - 1;
median = (arr[index[k]] + arr[index[k+1]])/2;
} else {
k = (Int_t)TMath::Floor(n / 2.0);
median = arr[index[k]];
}//if

delete [] index;

return median;
}//Median

//______________________________________________________________________________
void macroMedian()
{
Int_t n = 1000;
Double_t *x = new Double_t[n];

for (Int_t i=0; i<n; i++) x[i] = gRandom->Rndm();

cout << "Median (ROOT) = " << TMath::Median(n, x) << endl;
cout << "Median (C.S.) = " << XMedian(n, x) << endl;
cout << "Median (CRAN) = " << median(x, n) << endl;

delete [] x;
}//macroMedian
```

This archive was generated by hypermail 2b29 : Sun Jan 02 2005 - 05:50:09 MET