Mercurial > forge
annotate extra/NaN/src/kth_element.cpp @ 12720:52ca082757c2 octave-forge tip
Update copyright notices.
author | i7tiol |
---|---|
date | Sat, 27 Feb 2016 11:21:29 +0000 |
parents | 13815b367946 |
children |
rev | line source |
---|---|
7232 | 1 //------------------------------------------------------------------- |
2 // C-MEX implementation of kth element - this function is part of the NaN-toolbox. | |
3 // | |
7301
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
4 // usage: x = kth_element(X,k [,flag]) |
7232 | 5 // returns sort(X)(k) |
6 // | |
7301
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
7 // References: |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
8 // [1] https://secure.wikimedia.org/wikipedia/en/wiki/Selection_algorithm |
7232 | 9 // |
10 // This program is free software; you can redistribute it and/or modify | |
11 // it under the terms of the GNU General Public License as published by | |
12 // the Free Software Foundation; either version 3 of the License, or | |
13 // (at your option) any later version. | |
14 // | |
15 // This program is distributed in the hope that it will be useful, | |
16 // but WITHOUT ANY WARRANTY; without even the implied warranty of | |
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
18 // GNU General Public License for more details. | |
19 // | |
20 // You should have received a copy of the GNU General Public License | |
21 // along with this program; if not, see <http://www.gnu.org/licenses/>. | |
22 // | |
23 // | |
24 // Input: | |
7301
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
25 // X data vector, must be double/real |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
26 // k which element should be selected |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
27 // flag [optional]: |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
28 // 0: data in X might be reorded (partially sorted) in-place and |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
29 // is slightly faster because no local copy is generated |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
30 // data with NaN is not correctly handled. |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
31 // 1: data in X is never modified in-place, but a local copy is used. |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
32 // data with NaN is not correctly handled. |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
33 // 2: copies data and excludes all NaN's, the copying might be slower |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
34 // than 1, but it enables a faster selection algorithm. |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
35 // This is the save but slowest option |
7232 | 36 // |
37 // Output: | |
38 // x = sort(X)(k) | |
39 // | |
40 // $Id$ | |
8037 | 41 // Copyright (C) 2010,2011 Alois Schloegl <alois.schloegl@gmail.com> |
7232 | 42 // This function is part of the NaN-toolbox |
7889 | 43 // http://pub.ist.ac.at/~schloegl/matlab/NaN/ |
7232 | 44 // |
45 //------------------------------------------------------------------- | |
46 | |
47 | |
48 #include <math.h> | |
49 #include <stdlib.h> | |
7888 | 50 #include <stdint.h> |
7232 | 51 #include <string.h> |
52 #include "mex.h" | |
53 | |
54 | |
55 #ifdef tmwtypes_h | |
56 #if (MX_API_VER<=0x07020000) | |
57 typedef int mwSize; | |
7271
f8f0a8e59cc4
prepare next release; fix docu; xcovf: improve compatibility; quantile: empty set returns NaN
schloegl
parents:
7232
diff
changeset
|
58 typedef int mwIndex; |
7232 | 59 #endif |
60 #endif | |
61 | |
12692
13815b367946
use size_t instead of mwSize/mwIndex whenever possible
schloegl
parents:
12686
diff
changeset
|
62 /* |
13815b367946
use size_t instead of mwSize/mwIndex whenever possible
schloegl
parents:
12686
diff
changeset
|
63 math.h has isnan() defined for all sizes of floating point numbers, |
13815b367946
use size_t instead of mwSize/mwIndex whenever possible
schloegl
parents:
12686
diff
changeset
|
64 but c++ assumes isnan(double), causing possible conversions for float and long double |
13815b367946
use size_t instead of mwSize/mwIndex whenever possible
schloegl
parents:
12686
diff
changeset
|
65 */ |
13815b367946
use size_t instead of mwSize/mwIndex whenever possible
schloegl
parents:
12686
diff
changeset
|
66 #define ISNAN(a) (a!=a) |
13815b367946
use size_t instead of mwSize/mwIndex whenever possible
schloegl
parents:
12686
diff
changeset
|
67 |
7232 | 68 |
7272 | 69 #define SWAP(a,b) {temp = a; a=b; b=temp;} |
7232 | 70 |
7301
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
71 static void findFirstK(double *array, size_t left, size_t right, size_t k) |
7232 | 72 { |
7272 | 73 while (right > left) { |
12692
13815b367946
use size_t instead of mwSize/mwIndex whenever possible
schloegl
parents:
12686
diff
changeset
|
74 size_t pivotIndex = (left + right) / 2; |
7272 | 75 |
76 /* partition */ | |
77 double temp; | |
78 double pivotValue = array[pivotIndex]; | |
79 SWAP(array[pivotIndex], array[right]); | |
80 pivotIndex = left; | |
12692
13815b367946
use size_t instead of mwSize/mwIndex whenever possible
schloegl
parents:
12686
diff
changeset
|
81 for (size_t i = left; i <= right - 1; ++i ) { |
7301
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
82 // if (array[i] <= pivotValue || isnan(pivotValue)) // needed if data contains NaN's |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
83 if (array[i] <= pivotValue) |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
84 { |
7272 | 85 SWAP(array[i], array[pivotIndex]); |
86 ++pivotIndex; | |
87 } | |
88 } | |
89 SWAP(array[pivotIndex], array[right]); | |
90 | |
91 if (pivotIndex > k) | |
92 right = pivotIndex - 1; | |
93 else if (pivotIndex < k) | |
94 left = pivotIndex + 1; | |
95 else break; | |
7232 | 96 } |
97 } | |
98 | |
99 | |
100 void mexFunction(int POutputCount, mxArray* POutput[], int PInputCount, const mxArray *PInputs[]) | |
101 { | |
12692
13815b367946
use size_t instead of mwSize/mwIndex whenever possible
schloegl
parents:
12686
diff
changeset
|
102 size_t k, n; // running indices |
13815b367946
use size_t instead of mwSize/mwIndex whenever possible
schloegl
parents:
12686
diff
changeset
|
103 size_t szK, szX; |
7301
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
104 double *T,*X,*Y,*K; |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
105 char flag = 0; // default value |
7232 | 106 |
107 // check for proper number of input and output arguments | |
7301
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
108 if ( PInputCount < 2 || PInputCount > 3 ) { |
7272 | 109 mexPrintf("KTH_ELEMENT returns the K-th smallest element of vector X\n"); |
110 mexPrintf("\nusage:\tx = kth_element(X,k)\n"); | |
7301
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
111 mexPrintf("\nusage:\tx = kth_element(X,k,flag)\n"); |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
112 mexPrintf("\nflag=0: the elements in X can be modified in-place, and data with NaN's is not correctly handled. This can be useful for performance reasons, but it might modify data in-place and is not save for data with NaN's. You are warned.\n"); |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
113 mexPrintf("flag=1: prevents in-place modification of X using a local copy of the data, but does not handle data with NaN in the correct way.\n"); |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
114 mexPrintf("flag=2: prevents in-place modification of X using a local copy of the data and handles NaN's correctly. This is the save but slowest option.\n"); |
7232 | 115 |
116 mexPrintf("\nsee also: median, quantile\n\n"); | |
7301
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
117 mexErrMsgTxt("KTH_ELEMENT requires two or three input arguments\n"); |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
118 } |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
119 else if (PInputCount == 3) { |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
120 // check value of flag |
12692
13815b367946
use size_t instead of mwSize/mwIndex whenever possible
schloegl
parents:
12686
diff
changeset
|
121 size_t N = mxGetNumberOfElements(PInputs[2]); |
7301
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
122 if (N>1) |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
123 mexErrMsgTxt("KTH_ELEMENT: flag argument must be scalar\n"); |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
124 else if (N==1) { |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
125 switch (mxGetClassID(PInputs[2])) { |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
126 case mxLOGICAL_CLASS: |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
127 case mxCHAR_CLASS: |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
128 case mxINT8_CLASS: |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
129 case mxUINT8_CLASS: |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
130 flag = (char)*(uint8_t*)mxGetData(PInputs[2]); |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
131 break; |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
132 case mxDOUBLE_CLASS: |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
133 flag = (char)*(double*)mxGetData(PInputs[2]); |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
134 break; |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
135 case mxSINGLE_CLASS: |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
136 flag = (char)*(float*)mxGetData(PInputs[2]); |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
137 break; |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
138 case mxINT16_CLASS: |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
139 case mxUINT16_CLASS: |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
140 flag = (char)*(uint16_t*)mxGetData(PInputs[2]); |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
141 break; |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
142 case mxINT32_CLASS: |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
143 case mxUINT32_CLASS: |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
144 flag = (char)*(uint32_t*)mxGetData(PInputs[2]); |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
145 break; |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
146 case mxINT64_CLASS: |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
147 case mxUINT64_CLASS: |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
148 flag = (char)*(uint64_t*)mxGetData(PInputs[2]); |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
149 break; |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
150 case mxFUNCTION_CLASS: |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
151 case mxUNKNOWN_CLASS: |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
152 case mxCELL_CLASS: |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
153 case mxSTRUCT_CLASS: |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
154 default: |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
155 mexErrMsgTxt("KTH_ELEMENT: Type of 3rd input argument not supported."); |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
156 } |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
157 } |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
158 // else flag = default value |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
159 } |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
160 // else flag = default value |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
161 |
7232 | 162 if (POutputCount > 2) |
7301
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
163 mexErrMsgTxt("KTH_ELEMENT has only one output arguments."); |
7232 | 164 |
165 // get 1st argument | |
7301
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
166 if (mxIsComplex(PInputs[0]) || mxIsComplex(PInputs[1])) |
7232 | 167 mexErrMsgTxt("complex argument not supported (yet). "); |
168 if (!mxIsDouble(PInputs[0]) || !mxIsDouble(PInputs[1])) | |
169 mexErrMsgTxt("input arguments must be of type double . "); | |
170 // TODO: support of complex, and integer data | |
171 | |
172 | |
173 szK = mxGetNumberOfElements(PInputs[1]); | |
174 K = (double*)mxGetData(PInputs[1]); | |
175 | |
176 szX = mxGetNumberOfElements(PInputs[0]); | |
177 X = (double*)mxGetData(PInputs[0]); | |
178 | |
7301
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
179 if (flag==0) |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
180 T = X; |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
181 else { |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
182 //***** create temporary copy for avoiding unintended side effects (in-place sort of input data) */ |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
183 T = (double*)mxMalloc(szX*sizeof(double)); |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
184 if (flag==1) |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
185 memcpy(T,X,szX*sizeof(double)); |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
186 else { |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
187 /* do not copy NaN's */ |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
188 for (k=0,n=0; k < szX; k++) { |
12692
13815b367946
use size_t instead of mwSize/mwIndex whenever possible
schloegl
parents:
12686
diff
changeset
|
189 if (!ISNAN(X[k])) T[n++]=X[k]; |
7301
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
190 } |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
191 szX = n; |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
192 } |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
193 } |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
194 |
7232 | 195 /*********** create output arguments *****************/ |
196 POutput[0] = mxCreateDoubleMatrix(mxGetM(PInputs[1]),mxGetN(PInputs[1]),mxREAL); | |
197 Y = (double*) mxGetData(POutput[0]); | |
198 for (k=0; k < szK; k++) { | |
12686 | 199 if (K[k] > szX || K[k] < 1) |
7272 | 200 Y[k] = 0.0/0.0; // NaN: result undefined |
7232 | 201 else { |
12686 | 202 n = (size_t)(K[k]-1); // convert to zero-based indexing, round towards 0 |
7301
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
203 findFirstK(T, 0, szX-1, n); |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
204 Y[k] = T[n]; |
7232 | 205 } |
206 } | |
207 | |
7301
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
208 if (flag) mxFree(T); |
485d1594d155
addresses undesired side-effect off in-place sorting of data
schloegl
parents:
7272
diff
changeset
|
209 |
7232 | 210 return; |
211 } | |
212 |