Getting indexes of rows of matrix with more than n repetitions
in [Matlab]
|
From: Bruno Luong on 9 Jan 2010 14:07 Code for you Jan: /************************************************************************** * MATLAB mex function sort_idx.c * Calling syntax: * >> I = sort_idx(A) * returns I such that A(I) is sorted * Work in progress * Date: 09-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 QUICKSORT /* Pivot Strategy, use one of the above */ #define PIVOT MEDIAN3 /* 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; if (right > left) { pleft = pos+left; pright = pos+right; for (pi=pleft+1; pi<=pright; pi++) { Value = list[pival=(*pi)]; pj = pi-1; while ((pj>=pleft) && (list[*pj]>Value)) { /* Comparison */ *(pj+1) = *pj; /* Permute */ pj--; } *(pj+1) = pival; /* Insert */ } } 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[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 *pfirst, *plast; mwIndex tmp, ip, il, ir; plast=pos+right; pindex=pos+pivotIndex; pivotValue = list[tmp = *pindex]; /* Swap pos[pivotIndex] with pos[right] */ *pindex = *plast; ip = *plast = tmp; pfirst = pleft = pos+left; pright = plast-1; for (;;) { li = list[il = *pleft]; while ((li<pivotValue) || (li==pivotValue && il<ip)) /* <- test for stable quicksort */ li = list[il = *(++pleft)]; li = list[ir = *(pright)]; while ((pright>pleft) && ((li>pivotValue) || (li==pivotValue && ir>ip))) /* <- test for stable quicksort */ li = list[ir = *(--pright)]; if (pleft<pright) { *pleft = ir; *pright = il; ++pleft, --pright; } else { *pleft = *plast; *plast = il; return (mwIndex)(pleft-pos); /* Pointer arithmetic */ } } } /* PartitionV2 */ int CheckPartition(mwIndex left, mwIndex right, mwIndex pivotIndex) { mwIndex i; ARRAYTYPE pivotValue; pivotValue = list[pos[pivotIndex]]; for (i=left;i<=pivotIndex-1;i++) if (list[pos[i]]>pivotValue) return 0; for (i=right;i>=pivotIndex+1;i--) if (list[pos[i]]<pivotValue) return 0; return 1; } /*************************************************************************/ 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 = (left+right+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 (zero-based index) */ for (n=0; n<l; n++) pos[n]=n; } return; } /* AllocatePos */ /* Free the position */ void FreePos(void) { /* Free the array of position */ if (pos) mxFree(pos); pos = NULL; /* clean programming */ /* Free the array of values int64 */ /* if (list) mxFree(list); */ list = NULL; /* clean programming */ return; } /* FreePos */ /* Return the platform dependent maximum size of for 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(ARRAYTYPE *list, mwIndex l) { mwIndex n; /* 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) */ /* Convert zero-base indexing to Matlab one-base */ for (n=0; n<l; n++) (pos[n])++; return; } /* MainSortDbl */ /* Free a pointer */ void* MyFree(void* Ptr) { if (Ptr) mxFree(Ptr); return NULL; } /* Free */ /* 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; /* Allocate the array of index and get the data pointer */ AllocatePos(l); /* Pointers on data */ list = mxGetPr(A); /* Sort the indexes, the result is put on pos[] */ MainSortDbl(list, l); /* Return the index array */ dimI[0] = 0; dimI[1] = 0; I = mxCreateNumericArray(2, dimI, IdxClass(), mxREAL); mxSetPr(I, pos); mxSetM(I, m); mxSetN(I, n); plhs[0] = I; return; } /* sort_idx */
From: Bruno Luong on 9 Jan 2010 14:28 > > > > But simultaneously I try to find an efficient C-algorithm to create permutations: Choose K elements from a vector without repetitions and with order. I hope I do not confuse both programming threads. > > > > If you want you can compare with my implementation MIN/MAX Selection on FEX. > Sorry, I missread your post. Both problems have no common relation whatsoever. So please ignore. Bruno
From: Jan Simon on 9 Jan 2010 17:59 Dear Bruno! Nice piece of code. While I can often read some Matlab Haikus, which gain the efficiency from the compactness of Matlab, this sort squeezes the dull simplicity of C. Thanks! Jan
From: Bruno Luong on 10 Jan 2010 07:30 "Jan Simon" <matlab.THIS_YEAR(a)nMINUSsimon.de> wrote in message <hib1no$q8a$1(a)fred.mathworks.com>... > Dear Bruno! > > Nice piece of code. > While I can often read some Matlab Haikus, which gain the efficiency from the compactness of Matlab, this sort squeezes the dull simplicity of C. > Thanks! Jan Still I'm intrigued: why my code 2-3 times slower than Matlab. Is it due to double indexing because of the reference sorting as I suspected? Or Matlab library uses faster compiler (my code compiled with LCC is much slower - a factor 2 with respect to MSVC)? Algorithm optimization? Something else? I have cleaned up the code a little bit to make it slightly faster, but it get nowhere close to Matlab sorting speed. Bruno
From: Jan Simon on 10 Jan 2010 13:29
Dear Bruno! > I have cleaned up the code a little bit to make it slightly faster, but it get nowhere close to Matlab sorting speed. After I've read your code, I found this: http://www.netlib.org/napack/sort.f I cannot find the correct terms to describe my reaction - so I stay at :-) LCC3.8 from the net creates faster code. The MEXOPTS.BAT needs the surprising modification, that to RSP-file indicator like "@" can be used with this LCC. Let's ask Loren for the source of SORT - I do not know, if it is her code, but I assume, she knows the author at least. Kind regards, Jan |