view release on metacpan or  search on metacpan

src/cluster.c  view on Meta::CPAN

}

/* ************************************************************************ */

double
median(int n, double x[])
/*
Find the median of X(1), ... , X(N), using as much of the quicksort
algorithm as is needed to isolate it.
N.B. On exit, the array X is partially ordered.
Based on Alan J. Miller's median.f90 routine.
*/

{
    int i, j;
    int nr = n / 2;
    int nl = nr - 1;
    int even = 0;
    /* hi & lo are position limits encompassing the median. */
    int lo = 0;
    int hi = n-1;

    if (n == 2*nr) even = 1;
    if (n < 3) {
        if (n < 1) return 0.;
        if (n == 1) return x[0];
        return 0.5*(x[0]+x[1]);
    }

    /* Find median of 1st, middle & last values. */
    do {
        int loop;
        int mid = (lo + hi)/2;
        double result = x[mid];
        double xlo = x[lo];
        double xhi = x[hi];
        if (xhi<xlo) {
            double temp = xlo;
            xlo = xhi;
            xhi = temp;
        }
        if (result>xhi) result = xhi;
        else if (result<xlo) result = xlo;
        /* The basic quicksort algorithm to move all values <= the sort key
         * (XMED) to the left-hand end, and all higher values to the other end.
         */
        i = lo;
        j = hi;
        do {
            while (x[i]<result) i++;
            while (x[j]>result) j--;
            loop = 0;
            if (i<j) {
                double temp = x[i];
                x[i] = x[j];
                x[j] = temp;
                i++;
                j--;
                if (i<=j) loop = 1;
            }
        } while (loop); /* Decide which half the median is in. */

        if (even) {
            if (j == nl && i == nr) {
                /* Special case, n even, j = n/2 & i = j + 1, so the median is
                 * between the two halves of the series. Find max. of the first
                 * half & min. of the second half, then average.
                 */
                int k;
                double xmax = x[0];
                double xmin = x[n-1];
                for (k = lo; k <= j; k++) xmax = max(xmax, x[k]);
                for (k = i; k <= hi; k++) xmin = min(xmin, x[k]);
                return 0.5*(xmin + xmax);
            }
            if (j<nl) lo = i;
            if (i>nr) hi = j;
            if (i == j) {
                if (i == nl) lo = nl;
                if (j == nr) hi = nr;
            }
        }
        else {
            if (j<nr) lo = i;
            if (i>nr) hi = j;
            /* Test whether median has been isolated. */
            if (i == j && i == nr) return result;
        }
    }
    while (lo<hi-1);

    if (even) return (0.5*(x[nl]+x[nr]));
    if (x[lo]>x[hi]) {
        double temp = x[lo];
        x[lo] = x[hi];
        x[hi] = temp;
    }
    return x[nr];
}

/* ********************************************************************** */

static const double* sortdata = NULL; /* used in the quicksort algorithm */

/* ---------------------------------------------------------------------- */

static int
compare(const void* a, const void* b)
/* Helper function for sort. Previously, this was a nested function under
 * sort, which is not allowed under ANSI C.
 */
{
    const int i1 = *(const int*)a;
    const int i2 = *(const int*)b;
    const double term1 = sortdata[i1];
    const double term2 = sortdata[i2];

    if (term1 < term2) return -1;
    if (term1 > term2) return +1;
    return 0;
}