@@ -586,6 +586,34 @@ real THTensor_(maxall)(THTensor *tensor)
586
586
return theMax ;
587
587
}
588
588
589
+ static void THTensor_ (quickselectnoidx )(real * arr , long k , long elements , long stride );
590
+
591
+ real THTensor_ (medianall )(THTensor * tensor )
592
+ {
593
+ THArgCheck (tensor -> nDimension > 0 , 1 , "tensor must have one dimension" );
594
+ THArgCheck (THTensor_ (isContiguous )(tensor ), 1 , "input is not contiguous" );
595
+
596
+ real theMedian ;
597
+ ptrdiff_t numel ;
598
+ long k ;
599
+ THTensor * temp_ ;
600
+ real * temp__data ;
601
+
602
+ numel = THTensor_ (nElement )(tensor );
603
+ k = (numel - 1 ) >> 1 ;
604
+
605
+ temp_ = THTensor_ (newClone )(tensor );
606
+ temp__data = THTensor_ (data )(temp_ );
607
+
608
+ THTensor_ (quickselectnoidx )(temp__data , k , numel , 1 );
609
+
610
+ theMedian = temp__data [k ];
611
+
612
+ THTensor_ (free )(temp_ );
613
+
614
+ return theMedian ;
615
+ }
616
+
589
617
accreal THTensor_ (sumall )(THTensor * tensor )
590
618
{
591
619
accreal sum = 0 ;
@@ -2044,6 +2072,9 @@ void THTensor_(reshape)(THTensor *r_, THTensor *t, THLongStorage *size)
2044
2072
#define LONG_SWAP (AAA , BBB ) swap = AAA; AAA = BBB; BBB = swap
2045
2073
#define REAL_SWAP (AAA , BBB ) rswap = AAA; AAA = BBB; BBB = rswap
2046
2074
2075
+ #define ARR_SWAP (III , JJJ ) \
2076
+ REAL_SWAP(ARR(III), ARR(JJJ));
2077
+
2047
2078
#define BOTH_SWAP (III , JJJ ) \
2048
2079
REAL_SWAP(ARR(III), ARR(JJJ)); \
2049
2080
LONG_SWAP(IDX(III), IDX(JJJ))
@@ -2261,6 +2292,53 @@ void THTensor_(sort)(THTensor *rt_, THLongTensor *ri_, THTensor *t, int dimensio
2261
2292
}
2262
2293
}
2263
2294
2295
+ /* Implementation of the Quickselect algorithm, based on Nicolas Devillard's
2296
+ public domain implementation at http://ndevilla.free.fr/median/median/
2297
+ Adapted similarly to the above Quicksort algorithm.
2298
+ This version does not produce indices along with values. */
2299
+ static void THTensor_ (quickselectnoidx )(real * arr , long k , long elements , long stride )
2300
+ {
2301
+ long P , L , R , i , j , swap ;
2302
+ real rswap , piv ;
2303
+ L = 0 ;
2304
+ R = elements - 1 ;
2305
+
2306
+ do {
2307
+ if (R <= L ) /* One element only */
2308
+ return ;
2309
+
2310
+ if (R == L + 1 ) { /* Two elements only */
2311
+ if (ARR (L ) > ARR (R )) {
2312
+ ARR_SWAP (L , R );
2313
+ }
2314
+ return ;
2315
+ }
2316
+
2317
+ /* Use median of three for pivot choice */
2318
+ P = (L + R )>>1 ;
2319
+ ARR_SWAP (P , L + 1 );
2320
+ if (ARR (L + 1 ) > ARR (R )) { ARR_SWAP (L + 1 , R ); }
2321
+ if (ARR (L ) > ARR (R )) { ARR_SWAP (L , R ); }
2322
+ if (ARR (L + 1 ) > ARR (L )) { ARR_SWAP (L + 1 , L ); }
2323
+
2324
+ i = L + 1 ;
2325
+ j = R ;
2326
+ piv = ARR (L );
2327
+ do {
2328
+ do i ++ ; while (ARR (i ) < piv );
2329
+ do j -- ; while (ARR (j ) > piv );
2330
+ if (j < i )
2331
+ break ;
2332
+ ARR_SWAP (i , j );
2333
+ } while (1 );
2334
+ ARR_SWAP (L , j );
2335
+
2336
+ /* Re-set active partition */
2337
+ if (j <= k ) L = i ;
2338
+ if (j >= k ) R = j - 1 ;
2339
+ } while (1 );
2340
+ }
2341
+
2264
2342
/* Implementation of the Quickselect algorithm, based on Nicolas Devillard's
2265
2343
public domain implementation at http://ndevilla.free.fr/median/median/
2266
2344
Adapted similarly to the above Quicksort algorithm. */
0 commit comments