forge: extra/NaN/src/kth_element.cpp annotate

annotate extra/NaN/src/kth_element.cpp @ 12686:7c1bc8d8c406 octave-forge

fix range check of k-index

author	schloegl
date	Sat, 12 Sep 2015 07:33:55 +0000
parents	6a419bec96bb
children	13815b367946

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
7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	62
7272 b91d541b9286 improve speed, avoid recursion, simplify code schloegl parents: 7271 diff changeset	63 #define SWAP(a,b) {temp = a; a=b; b=temp;}
7232 7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	64
7301 485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	65 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	66 {
7272 b91d541b9286 improve speed, avoid recursion, simplify code schloegl parents: 7271 diff changeset	67 while (right > left) {
b91d541b9286 improve speed, avoid recursion, simplify code schloegl parents: 7271 diff changeset	68 mwIndex pivotIndex = (left + right) / 2;
b91d541b9286 improve speed, avoid recursion, simplify code schloegl parents: 7271 diff changeset	69
b91d541b9286 improve speed, avoid recursion, simplify code schloegl parents: 7271 diff changeset	70 /* partition */
b91d541b9286 improve speed, avoid recursion, simplify code schloegl parents: 7271 diff changeset	71 double temp;
b91d541b9286 improve speed, avoid recursion, simplify code schloegl parents: 7271 diff changeset	72 double pivotValue = array[pivotIndex];
b91d541b9286 improve speed, avoid recursion, simplify code schloegl parents: 7271 diff changeset	73 SWAP(array[pivotIndex], array[right]);
b91d541b9286 improve speed, avoid recursion, simplify code schloegl parents: 7271 diff changeset	74 pivotIndex = left;
b91d541b9286 improve speed, avoid recursion, simplify code schloegl parents: 7271 diff changeset	75 for (mwIndex i = left; i <= right - 1; ++i ) {
7301 485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	76 // 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	77 if (array[i] <= pivotValue)
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	78 {
7272 b91d541b9286 improve speed, avoid recursion, simplify code schloegl parents: 7271 diff changeset	79 SWAP(array[i], array[pivotIndex]);
b91d541b9286 improve speed, avoid recursion, simplify code schloegl parents: 7271 diff changeset	80 ++pivotIndex;
b91d541b9286 improve speed, avoid recursion, simplify code schloegl parents: 7271 diff changeset	81 }
b91d541b9286 improve speed, avoid recursion, simplify code schloegl parents: 7271 diff changeset	82 }
b91d541b9286 improve speed, avoid recursion, simplify code schloegl parents: 7271 diff changeset	83 SWAP(array[pivotIndex], array[right]);
b91d541b9286 improve speed, avoid recursion, simplify code schloegl parents: 7271 diff changeset	84
b91d541b9286 improve speed, avoid recursion, simplify code schloegl parents: 7271 diff changeset	85 if (pivotIndex > k)
b91d541b9286 improve speed, avoid recursion, simplify code schloegl parents: 7271 diff changeset	86 right = pivotIndex - 1;
b91d541b9286 improve speed, avoid recursion, simplify code schloegl parents: 7271 diff changeset	87 else if (pivotIndex < k)
b91d541b9286 improve speed, avoid recursion, simplify code schloegl parents: 7271 diff changeset	88 left = pivotIndex + 1;
b91d541b9286 improve speed, avoid recursion, simplify code schloegl parents: 7271 diff changeset	89 else break;
7232 7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	90 }
7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	91 }
7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	92
7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	93
7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	94 void mexFunction(int POutputCount, mxArray* POutput[], int PInputCount, const mxArray *PInputs[])
7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	95 {
7272 b91d541b9286 improve speed, avoid recursion, simplify code schloegl parents: 7271 diff changeset	96 mwIndex k, n; // running indices
7301 485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	97 mwSize szK, szX;
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	98 double T,X,Y,K;
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	99 char flag = 0; // default value
7232 7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	100
7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	101 // 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	102 if ( PInputCount < 2 \|\| PInputCount > 3 ) {
7272 b91d541b9286 improve speed, avoid recursion, simplify code schloegl parents: 7271 diff changeset	103 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	104 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	105 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	106 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	107 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	108 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	109
7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	110 mexPrintf("\nsee also: median, quantile\n\n");
7301 485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	111 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	112 }
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	113 else if (PInputCount == 3) {
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	114 // check value of flag
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	115 mwSize N = mxGetNumberOfElements(PInputs[2]);
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	116 if (N>1)
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	117 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	118 else if (N==1) {
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	119 switch (mxGetClassID(PInputs[2])) {
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	120 case mxLOGICAL_CLASS:
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	121 case mxCHAR_CLASS:
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	122 case mxINT8_CLASS:
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	123 case mxUINT8_CLASS:
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	124 flag = (char)(uint8_t)mxGetData(PInputs[2]);
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	125 break;
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	126 case mxDOUBLE_CLASS:
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	127 flag = (char)(double)mxGetData(PInputs[2]);
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	128 break;
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	129 case mxSINGLE_CLASS:
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	130 flag = (char)(float)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 mxINT16_CLASS:
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	133 case mxUINT16_CLASS:
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	134 flag = (char)(uint16_t)mxGetData(PInputs[2]);
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	135 break;
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	136 case mxINT32_CLASS:
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	137 case mxUINT32_CLASS:
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	138 flag = (char)(uint32_t)mxGetData(PInputs[2]);
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	139 break;
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	140 case mxINT64_CLASS:
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	141 case mxUINT64_CLASS:
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	142 flag = (char)(uint64_t)mxGetData(PInputs[2]);
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	143 break;
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	144 case mxFUNCTION_CLASS:
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	145 case mxUNKNOWN_CLASS:
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	146 case mxCELL_CLASS:
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	147 case mxSTRUCT_CLASS:
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	148 default:
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	149 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	150 }
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	151 }
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	152 // else flag = default value
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	153 }
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	154 // else flag = default value
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	155
7232 7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	156 if (POutputCount > 2)
7301 485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	157 mexErrMsgTxt("KTH_ELEMENT has only one output arguments.");
7232 7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	158
7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	159 // get 1st argument
7301 485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	160 if (mxIsComplex(PInputs[0]) \|\| mxIsComplex(PInputs[1]))
7232 7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	161 mexErrMsgTxt("complex argument not supported (yet). ");
7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	162 if (!mxIsDouble(PInputs[0]) \|\| !mxIsDouble(PInputs[1]))
7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	163 mexErrMsgTxt("input arguments must be of type double . ");
7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	164 // TODO: support of complex, and integer data
7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	165
7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	166
7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	167 szK = mxGetNumberOfElements(PInputs[1]);
7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	168 K = (double*)mxGetData(PInputs[1]);
7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	169
7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	170 szX = mxGetNumberOfElements(PInputs[0]);
7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	171 X = (double*)mxGetData(PInputs[0]);
7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	172
7301 485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	173 if (flag==0)
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	174 T = X;
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	175 else {
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	176 //***** 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	177 T = (double)mxMalloc(szXsizeof(double));
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	178 if (flag==1)
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	179 memcpy(T,X,szX*sizeof(double));
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	180 else {
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	181 /* do not copy NaN's */
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	182 for (k=0,n=0; k < szX; k++) {
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	183 if (!isnan(X[k])) T[n++]=X[k];
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	184 }
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	185 szX = n;
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	186 }
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	187 }
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	188
7232 7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	189 /********* create output arguments ***************/
7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	190 POutput[0] = mxCreateDoubleMatrix(mxGetM(PInputs[1]),mxGetN(PInputs[1]),mxREAL);
7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	191 Y = (double*) mxGetData(POutput[0]);
7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	192 for (k=0; k < szK; k++) {
12686 7c1bc8d8c406 fix range check of k-index schloegl parents: 8037 diff changeset	193 if (K[k] > szX \|\| K[k] < 1)
7272 b91d541b9286 improve speed, avoid recursion, simplify code schloegl parents: 7271 diff changeset	194 Y[k] = 0.0/0.0; // NaN: result undefined
7232 7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	195 else {
12686 7c1bc8d8c406 fix range check of k-index schloegl parents: 8037 diff changeset	196 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	197 findFirstK(T, 0, szX-1, n);
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	198 Y[k] = T[n];
7232 7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	199 }
7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	200 }
7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	201
7301 485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	202 if (flag) mxFree(T);
485d1594d155 addresses undesired side-effect off in-place sorting of data schloegl parents: 7272 diff changeset	203
7232 7df85b226263 add kth_element and improve median (performance gain) schloegl parents: diff changeset	204 return;
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

Mercurial > forge

annotate extra/NaN/src/kth_element.cpp @ 12686:7c1bc8d8c406 octave-forge