median.c

这是用matlab编写的支持向量机的函数工具箱

/* median.c: Return the median from an array.
 *
 * The following 3 routines are used to find the median of an array "a" of "k"
 * numbers in expected linear time. The value of k must be at least 1. The
 * elements in "a" are reordered by the function. The algorithm uses selection
 * and randomized partitioning, see for example Cormen, Leiserson and Rivest,
 * "Introduction to Algorithms", MIT Press, 1993.
 *
 * Copyright (c) 1996 Carl Edward Rasmussen */

#include <stdlib.h>
#include "util.h"

static real select(real *, int, int);                   /* private functions */
static int  partition(real *, int);

real median(real *a, int k)               /* return the median from an array */
{
  if (k/2 != (k+1)/2)                                         /* if k is odd */
    return select(a, k, k/2+1);          /* then the median is in the middle */
  else                               /* otherwise the mean of the two middle */
    return 0.5*(select(a, k, k/2) + select(a, k, k/2+1));         /* numbers */
}

/* Recursive function that returns the i'th smallest element from an array, 
 * i=1..k. The elements are rearranged by the function. */

static real select(real *a, int k, int i)
{
  static int q;

  if (k == 1) return a[0];
  q = partition(a, k)+1;
  if (i <= q)
    return select(a, q, i);
  else
    return select(&a[q], k-q, i-q);
}

/* Partition an array around an element chosen at random and retun the index of
 * the partiton element. Upon returning all the elements with indexes smaller
 * than or equal to the index of the partition element have values that are
 * less than or equal to the partition element; the rest of the array have
 * values larger than or equal to the partition element. */

static int partition(real *a, int k)
{
  static real x, temp;
  static int  i, j;

  x = a[(int) (k*rand()/(RAND_MAX+1.0))]; i = -1; j = k;
  for (;;) {
    do j--; while (a[j] > x);
    do i++; while (a[i] < x);
    if (i < j) { temp = a[i]; a[i] = a[j]; a[j] = temp; } else return j;
  }
}