Mercurial > forge

annotate extra/NaN/src/kth_element.cpp @ 12692:13815b367946 octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
use size_t instead of mwSize/mwIndex whenever possible
author schloegl
date Sat, 12 Sep 2015 14:59:20 +0000
parents 7c1bc8d8c406
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
7232
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
1 //-------------------------------------------------------------------
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
2 // C-MEX implementation of kth element - this function is part of the NaN-toolbox.
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
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
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
5 // returns sort(X)(k)
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
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
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
9 //
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
10 // This program is free software; you can redistribute it and/or modify
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
11 // it under the terms of the GNU General Public License as published by
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
12 // the Free Software Foundation; either version 3 of the License, or
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
13 // (at your option) any later version.
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
14 //
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
15 // This program is distributed in the hope that it will be useful,
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
16 // but WITHOUT ANY WARRANTY; without even the implied warranty of
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
18 // GNU General Public License for more details.
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
19 //
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
20 // You should have received a copy of the GNU General Public License
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
21 // along with this program; if not, see <http://www.gnu.org/licenses/>.
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
22 //
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
23 //
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
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
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
36 //
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
37 // Output:
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
38 // x = sort(X)(k)
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
39 //
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
40 // $Id$
8037
6a419bec96bb update contact e-mail and www address
schloegl
parents: 7889
diff changeset
41 // Copyright (C) 2010,2011 Alois Schloegl <alois.schloegl@gmail.com>
7232
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
42 // This function is part of the NaN-toolbox
7889
c101c486d80a fix web address
schloegl
parents: 7888
diff changeset
43 // http://pub.ist.ac.at/~schloegl/matlab/NaN/
7232
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
44 //
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
45 //-------------------------------------------------------------------
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
46
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
47
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
48 #include <math.h>
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
49 #include <stdlib.h>
7888
b9f35668b55e replace <inttypes.h> with <stdint.h>
schloegl
parents: 7301
diff changeset
50 #include <stdint.h>
7232
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
51 #include <string.h>
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
52 #include "mex.h"
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
53
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
54
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
55 #ifdef tmwtypes_h
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
56 #if (MX_API_VER<=0x07020000)
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
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
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
59 #endif
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
60 #endif
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
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
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
68
7272
b91d541b9286 improve speed, avoid recursion, simplify code
schloegl
parents: 7271
diff changeset
69 #define SWAP(a,b) {temp = a; a=b; b=temp;}
7232
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
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
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
72 {
7272
b91d541b9286 improve speed, avoid recursion, simplify code
schloegl
parents: 7271
diff changeset
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
b91d541b9286 improve speed, avoid recursion, simplify code
schloegl
parents: 7271
diff changeset
75
b91d541b9286 improve speed, avoid recursion, simplify code
schloegl
parents: 7271
diff changeset
76 /* partition */
b91d541b9286 improve speed, avoid recursion, simplify code
schloegl
parents: 7271
diff changeset
77 double temp;
b91d541b9286 improve speed, avoid recursion, simplify code
schloegl
parents: 7271
diff changeset
78 double pivotValue = array[pivotIndex];
b91d541b9286 improve speed, avoid recursion, simplify code
schloegl
parents: 7271
diff changeset
79 SWAP(array[pivotIndex], array[right]);
b91d541b9286 improve speed, avoid recursion, simplify code
schloegl
parents: 7271
diff changeset
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
b91d541b9286 improve speed, avoid recursion, simplify code
schloegl
parents: 7271
diff changeset
85 SWAP(array[i], array[pivotIndex]);
b91d541b9286 improve speed, avoid recursion, simplify code
schloegl
parents: 7271
diff changeset
86 ++pivotIndex;
b91d541b9286 improve speed, avoid recursion, simplify code
schloegl
parents: 7271
diff changeset
87 }
b91d541b9286 improve speed, avoid recursion, simplify code
schloegl
parents: 7271
diff changeset
88 }
b91d541b9286 improve speed, avoid recursion, simplify code
schloegl
parents: 7271
diff changeset
89 SWAP(array[pivotIndex], array[right]);
b91d541b9286 improve speed, avoid recursion, simplify code
schloegl
parents: 7271
diff changeset
90
b91d541b9286 improve speed, avoid recursion, simplify code
schloegl
parents: 7271
diff changeset
91 if (pivotIndex > k)
b91d541b9286 improve speed, avoid recursion, simplify code
schloegl
parents: 7271
diff changeset
92 right = pivotIndex - 1;
b91d541b9286 improve speed, avoid recursion, simplify code
schloegl
parents: 7271
diff changeset
93 else if (pivotIndex < k)
b91d541b9286 improve speed, avoid recursion, simplify code
schloegl
parents: 7271
diff changeset
94 left = pivotIndex + 1;
b91d541b9286 improve speed, avoid recursion, simplify code
schloegl
parents: 7271
diff changeset
95 else break;
7232
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
96 }
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
97 }
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
98
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
99
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
100 void mexFunction(int POutputCount, mxArray* POutput[], int PInputCount, const mxArray *PInputs[])
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
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
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
106
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
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
b91d541b9286 improve speed, avoid recursion, simplify code
schloegl
parents: 7271
diff changeset
109 mexPrintf("KTH_ELEMENT returns the K-th smallest element of vector X\n");
b91d541b9286 improve speed, avoid recursion, simplify code
schloegl
parents: 7271
diff changeset
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
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
115
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
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
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
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
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
164
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
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
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
167 mexErrMsgTxt("complex argument not supported (yet). ");
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
168 if (!mxIsDouble(PInputs[0]) || !mxIsDouble(PInputs[1]))
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
169 mexErrMsgTxt("input arguments must be of type double . ");
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
170 // TODO: support of complex, and integer data
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
171
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
172
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
173 szK = mxGetNumberOfElements(PInputs[1]);
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
174 K = (double*)mxGetData(PInputs[1]);
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
175
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
176 szX = mxGetNumberOfElements(PInputs[0]);
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
177 X = (double*)mxGetData(PInputs[0]);
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
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
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
195 /*********** create output arguments *****************/
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
196 POutput[0] = mxCreateDoubleMatrix(mxGetM(PInputs[1]),mxGetN(PInputs[1]),mxREAL);
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
197 Y = (double*) mxGetData(POutput[0]);
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
198 for (k=0; k < szK; k++) {
12686
7c1bc8d8c406 fix range check of k-index
schloegl
parents: 8037
diff changeset
199 if (K[k] > szX || K[k] < 1)
7272
b91d541b9286 improve speed, avoid recursion, simplify code
schloegl
parents: 7271
diff changeset
200 Y[k] = 0.0/0.0; // NaN: result undefined
7232
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
201 else {
12686
7c1bc8d8c406 fix range check of k-index
schloegl
parents: 8037
diff changeset
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
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
205 }
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
206 }
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
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
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
210 return;
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
211 }
7df85b226263 add kth_element and improve median (performance gain)
schloegl
parents:
diff changeset
212