Prev: About Function Handle
Next: Rudimentary network emulator(RUNE)
From: Bruno Luong on 10 Jan 2010 17:06
There are few bugs in my previous indexing sorting code. Here is a new one. The default quicksort is out: none of the pivoting strategy I tested is completely satisfied.

Bruno

/**************************************************************************
* MATLAB mex function sort_idx.c
* Calling syntax:
* >> I = sort_idx(A)
* returns I such that A(I) is sorted
*
* Work in still in progress...
*
* Compile on 32-bit platform
* mex -O -v sort_idx.c
* On 64-bit platform
* mex -v -O -largeArrayDims sort_idx.c
*
* Author Bruno Luong <brunoluong@?????.com>
* Date: 10-Jan-2010
*************************************************************************/

#include "mex.h"
#include "matrix.h"
#include "string.h"
#include "math.h"

/* Define correct type depending on platform */
#if defined(_MSC_VER) || defined(__BORLANDC__)
typedef unsigned __int32 uint32;
typedef signed __int32 int32;
typedef unsigned __int64 ulong64;
typedef signed __int64 long64;
#define INTMIN64 0x8000000000000000
#else
typedef unsigned int uint32;
typedef signed int int32;
typedef unsigned long long ulong64;
typedef signed long long long64;
#define INTMIN64 0x8000000000000000ULL
#endif

/* Different options for Pivot selection */
/* Different options for Pivot selection */
#define MIDPOINT 0
#define MEDIAN3 1
#define MEDIANMEDIANS 2

/* Pivot Strategy, use one of the above */
#define PIVOT MEDIAN3

#define INSERTSORT 0
#define SHELLSORT 1
#define QUICKSORT 2
#define HEAPSORT 3
#define INTROSORT 4

/* Method for sorting small arrays */
/* INSERTSORT | SHELLSORT */
#define BASESORT INSERTSORT

/* Method dor sorting large arrays */
/* QUICKSORT | HEAPSORT | INTROSORT */
#define MAINSORT INTROSORT

/* Increment for shellshort */
#define SizeIncs 3
mwIndex incs[SizeIncs] = { 7, 3, 1 };

/* A limit below the recusion stops, we which we switch to InsertionSort */
/* This must be smaller than */
#define SWITCH_ALGO_THRESHOLD 12

/* TYPE */
/* Type of elements to be sorted */
#define ARRAYTYPE double

/* Global variables, used to avoid stacking them during recusive call since
they do not change */
mwIndex *pos;
ARRAYTYPE *list;

/*************************************************************************/
/* Insertion sort, this is stable */
void InsertionSort(mwIndex left, mwIndex right)
{
mwIndex pival;
mwIndex *pi, *pj, *pleft, *pright;
ARRAYTYPE Value;

pleft = pos+left;
pright = pos+right;
for (pi=pleft; pi<pright; pi++) {
Value = list[pival=(*(pi+1))];
pj = pi;

while ((pj>=pleft) && (list[*pj]>Value)) { /* Comparison */
*(pj+1) = *pj; /* Permute */
pj--;
}
*(pj+1) = pival; /* Insert */
} /* for */

return;
} /* InsertionSort */

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

/* Shell sort */
void ShellSort(mwIndex left, mwIndex right)
{
mwIndex pival, pjval, h, k;
mwIndex *pi, *pj, *ph;
mwIndex *pleft, *pright;
ARRAYTYPE Value;

if (right > left) {

pleft = pos+left;
pright = pos+right;

/* Loop on gap size */
for ( k = 0; k < SizeIncs; k++) {
h = incs[k];
ph = pos + h;
/* Insertion loop */
for (pi=pleft+h; pi<=pright; pi++)
{
Value = list[pival=*pi];
pj = pi;
while ((pj>=ph) && (list[pjval=*(pj-h)]>Value)) { /* Comparison */
*pj = pjval; /* Permute */
pj -= h;
}
*pj = pival; /* Insert */
}
}
}

return;
} /* ShellSort */

/*************************************************************************/
/*Find the index of the Median of the elements
of array that occur at every "shift" positions.*/
mwIndex findMedianIndex(mwIndex left, mwIndex right, mwIndex shift)
{
mwIndex tmp, groups, n;
ARRAYTYPE minValue;
mwIndex *pi, *pj, *pk, *pright, *pminIndex;

groups = (right-left)/shift + 1;
pk = pos + (n = left + (groups/2)*shift);
pright = pos + right;
for (pi=pos+left; pi<=pk; pi+= shift)
{
pminIndex = pi;
minValue = list[*pminIndex];

for (pj=pi; pj<=pright; pj+=shift)
if (list[*pj]<minValue) /* Comparison */
minValue = list[*(pminIndex=pj)];
/* Swap pos[i] with pos[minIndex] */
tmp = *pi;
*pi = *pminIndex;
*pminIndex = tmp;
}

return n;

} /* findMedianIndex */

/*Computes the median of each group of 5 elements and stores
it as the first element of the group (left). Recursively does this
till there is only one group and hence only one Median */
mwIndex findMedianOfMedians(mwIndex left, mwIndex right)
{
mwIndex i, shift, step, tmp;
mwIndex endIndex, medianIndex;

if (left==right) return left;

shift = 1;
while (shift <= (right-left))
{
step=shift*5;
for (i=left; i<=right; i+=step)
{
if ((endIndex=i+step-1)>=right)
endIndex=right;
medianIndex = findMedianIndex(i, endIndex, shift);
/* Swap pos[i] with pos[medianIndex] */
tmp = pos[i];
pos[i] = pos[medianIndex];
pos[medianIndex] = tmp;
}
shift = step;
}
return left;
} /* findMedianOfMedians */

/*************************************************************************/
/*Computes the median of three points (left,right,and mid) */
mwIndex findMedianThree(mwIndex left, mwIndex right)
{
ARRAYTYPE vleft, vright, vmid;
mwIndex mid;

if (left==right) return left;

vleft = list[pos[left]];
vright = list[pos[right]];
vmid = list[pos[mid = (left+right+1)/2]];

if (vleft<vright)
{
if (vmid>vright)
return right;
else if (vmid<vleft)
return left;
else
return mid;

} else { /* (vleft>=vright) */

if (vmid>vleft)
return left;
else if (vmid<vright)
return right;
else
return mid;

}
} /* findMedianThree */

/*************************************************************************/
/* Partitioning the list around pivot pivotValue := l[pivotIndex];
* After runing, at exit we obtain:
l[left]...l[index-1] < pivotValue <= l[index] ... l[right]
where l[i] := list[pos[i]] for all i */
mwIndex PartitionV2(mwIndex left, mwIndex right, mwIndex pivotIndex) {

ARRAYTYPE pivotValue, li;
mwIndex *pindex, *pright, *pleft;
mwIndex *plast;
mwIndex ip, il, ir;

plast=pos+right;
pindex=pos+pivotIndex;
pivotValue = list[ip = *pindex];
/* Swap pos[pivotIndex] with pos[right] */
*pindex = *plast;
*plast = ip;

pleft = pos+left;
pright = plast-1;
while (1) {

li = list[il = *pleft];
while ((li<pivotValue) ||
(li==pivotValue && il<ip)) /* <- test for stable sort */
li = list[il = *(++pleft)];

li = list[ir = *(pright)];
while ((pright>pleft) &&
((li>pivotValue) ||
(li==pivotValue && ir>ip))) /* <- test for stable sort */
li = list[ir = *(--pright)];

if (pleft<pright) {
*(pleft++) = ir;
*(pright--) = il;
}
else {
/* swap left with the pivot */
*pleft = *plast;
*plast = il;
return (mwIndex)(pleft-pos); /* Pointer arithmetic */
}
}
} /* PartitionV2 */

/*************************************************************************/
void DownHeap(mwIndex i, mwIndex n, mwIndex lo) {
ARRAYTYPE d;
mwIndex child, posd, nhalf;
mwIndex *plo, *pchild;

plo = pos + lo;
posd = *(plo+i-1);
d = list[posd];
nhalf = n/2;
while (i<=nhalf) {
child = 2*i;
pchild = plo + child;
if ((child < n) &&
(list[*(pchild-1)] < list[*pchild]))
{
child++;
pchild++;
}
if (d >= list[*(pchild-1)]) break;
*(plo+i-1) = *(pchild-1);
i = child;
}
*(plo+i-1) = posd;

return;
} /* DownHeap */

/*http://users.encs.concordia.ca/~chvatal/notes/hsort.html*/
void HeapSort(mwIndex lo, mwIndex hi) { /* (lo,hi) included */
mwIndex n, i, tmp;
mwIndex *plo;

hi++;
n = hi-lo;

for (i=n/2; i>=1; i--)
DownHeap(i, n, lo);

plo = pos + lo;
for (i=n; i>1; i--) {

tmp = *(plo);
*(plo) = *(plo+i-1);
*(plo+i-1) = tmp;

DownHeap(1, i-1, lo);
}
return;

} /* HeapSort */

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

void introsort_loop(mwIndex lo, mwIndex hi, int depth_limit)
{
mwIndex pivotIndex;
while (hi+1 >= lo + SWITCH_ALGO_THRESHOLD) {
if (depth_limit == 0) {
HeapSort(lo, hi);
return;
}
depth_limit--;

#if (PIVOT==MEDIANMEDIANS)
pivotIndex = findMedianOfMedians(lo, hi);
#elif (PIVOT==MEDIAN3)
pivotIndex = findMedianThree(lo, hi);
#else /* MIDPOINT */
pivotIndex = (lo+hi+1)/2;
#endif
pivotIndex = PartitionV2(lo, hi, pivotIndex);
introsort_loop(pivotIndex+1, hi, depth_limit);
hi = pivotIndex-1;
}
return;
} /* introsort_loop */

int floor_lg(mwIndex a) {
return (int)(floor(log((double)a)/log(2.0)));
} /* floor_lg */


void IntroSort(mwIndex left, mwIndex right)
{
int depthlimit;
depthlimit=2*floor_lg(right-left+1);

introsort_loop(left, right, depthlimit);
#if (BASESORT==INSERTSORT)
InsertionSort(left, right);
#else /* #if (BASESORT==SHELLSORT)*/
ShellSort(left, right);
#endif

return;
} /* IntroSort */

/*************************************************************************/
/* Recursive engine quicksort */
void quicksort(mwIndex left, mwIndex right) {

mwIndex pivotIndex;

/* Switch to Insertion sort for small size array */
if (right+1 < left + SWITCH_ALGO_THRESHOLD)
{
#if (BASESORT==INSERTSORT)
InsertionSort(left, right);
#else /* #if (BASESORT==SHELLSORT)*/
ShellSort(left, right);
#endif
}
else
{
#if (PIVOT==MEDIANMEDIANS)
pivotIndex = findMedianOfMedians(left, right);
#elif (PIVOT==MEDIAN3)
pivotIndex = findMedianThree(left, right);
#else /* MIDPOINT */
pivotIndex = (left+right+1)/2;
#endif

pivotIndex = PartitionV2(left, right, pivotIndex);
quicksort(left, pivotIndex-1);
quicksort(pivotIndex+1, right);
}

return;
} /* quicksort */

/* Allocate and Initialize the table of position before sorting */
void AllocatePos(mwIndex l)
{
mwIndex n;

/* clean programming */
pos=NULL;

if (l>0) {
pos = mxMalloc(sizeof(mwIndex)*l);
if (pos==NULL)
mexErrMsgTxt("Out of memory (pos).");
/* Initialize the array of position (one-based Matlab index) */
pos--;
for (n=1; n<=l; n++) pos[n]=n;
pos++;
}

return;
} /* AllocatePos */


/* Return the platform dependent the class corresponding to mwIndex */
mxClassID IdxClass(void) {
switch (sizeof(mwIndex)) {
case 4:
return mxINT32_CLASS; /* 2^31-1 */
case 8:
return mxINT64_CLASS; /* 2^48-1 */
}
} /* MaxSize */

/* MainSortDbl, sorting interface for double */
void MainSortDbl(mwIndex l) {

/* Allocate the array of index and get the data pointer */
AllocatePos(l);

/* We reduce the pointer by 1 temporary
to cope with Matlab one-base indexing */
list--;
/* Work for non-empty array */
if (l>1 && list!=NULL && pos!=NULL) {
/* Call the recursive engine */
#if (MAINSORT==QUICKSORT)
quicksort(0, l-1);
#elif (MAINSORT==INTROSORT)
IntroSort(0, l-1);
#elif (MAINSORT==HEAPSORT)
HeapSort(0, l-1);
#else
#error "Non valid MAINSORT"
#endif
} /* if (l>0) */
list++;

return;

} /* MainSortDbl */

/* Gateway routine sort_idx */
void mexFunction(int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[]) {

const mxArray *A;
mwIndex m, n, l;
mxClassID DataClass;
mxArray* I;
mwSize dimI[2];

if (nrhs!=1)
mexErrMsgTxt("One input argument is required.");

/* Check data type of input argument */

A = prhs[0];
if (mxIsSparse(A))
mexErrMsgTxt("First input argument must be a full matrix.");

DataClass = mxGetClassID(A);
if ((DataClass != mxDOUBLE_CLASS))
mexErrMsgTxt("A must be of class DOUBLE.");

if (mxGetNumberOfDimensions(A) != 2)
mexErrMsgTxt("First input argument must be a vector.");

m = (mwIndex)mxGetM(A);
n = (mwIndex)mxGetN(A);

if ((m > 1) && (n > 1))
mexErrMsgTxt("First input argument must be a vector.");

/* Get the total number elements of A */
l = m*n;

/* Pointers on data */
list = mxGetPr(A);

/* Sort the indexes, the result is put on pos[] */
MainSortDbl(l);

/* Return the index array */
dimI[0] = 0; dimI[1] = 0;
I = mxCreateNumericArray(2, dimI, IdxClass(), mxREAL);
mxSetData(I, pos);
mxSetM(I, m);
mxSetN(I, n);
plhs[0] = I;

return;

} /* sort_idx */
First  |  Prev  | 
Pages: 2 3 4 5 6 7 8 9 10 11 12
Prev: About Function Handle
Next: Rudimentary network emulator(RUNE)