00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef __MATHMATICS_H_
00014 #define __MATHMATICS_H_
00015
00016 #include "lib/common.h"
00017 #include "lib/io.h"
00018 #include "lib/lapack.h"
00019 #include "base/SGObject.h"
00020 #include "base/Parallel.h"
00021
00022 #include <math.h>
00023 #include <stdio.h>
00024 #include <float.h>
00025 #include <pthread.h>
00026 #include <unistd.h>
00027 #include <sys/types.h>
00028 #include <sys/time.h>
00029 #include <time.h>
00030
00031
00032 #ifdef _WIN32
00033 #include <float.h>
00034 #ifndef finite
00035 #define finite _finite
00036 #endif
00037
00038 #ifndef isnan
00039 #define isnan _isnan
00040 #endif
00041 #endif
00042
00043
00044 #ifdef CYGWIN
00045 #ifndef finite
00046 #include <ieeefp.h>
00047 #endif
00048 #endif
00049
00050 #ifndef NAN
00051 #include <stdlib.h>
00052 #define NAN (strtod("NAN",NULL))
00053 #endif
00054
00055 #ifdef SUNOS
00056 extern "C" int finite(double);
00057 #endif
00058
00059
00060 #define RADIX_STACK_SIZE 512
00061
00062
00063 #define radix_push(a, n, i) sp->sa = a, sp->sn = n, (sp++)->si = i
00064 #define radix_pop(a, n, i) a = (--sp)->sa, n = sp->sn, i = sp->si
00065
00067 template <class T> struct radix_stack_t
00068 {
00070 T *sa;
00072 size_t sn;
00074 unsigned short si;
00075 };
00076
00078 struct pair
00079 {
00081 int idx1;
00083 int idx2;
00084 };
00085
00087 struct thread_qsort
00088 {
00090 DREAL* output;
00092 INT* index;
00094 INT size;
00095
00097 INT* qsort_threads;
00099 INT sort_limit;
00101 INT num_threads;
00102 };
00103
00107 class CMath : public CSGObject
00108 {
00109 public:
00113
00114 CMath();
00115
00117 virtual ~CMath();
00119
00123
00125
00126 template <class T>
00127 static inline T min(T a, T b)
00128 {
00129 return ((a)<=(b))?(a):(b);
00130 }
00131
00133 template <class T>
00134 static inline T max(T a, T b)
00135 {
00136 return ((a)>=(b))?(a):(b);
00137 }
00138
00140 template <class T>
00141 static inline T clamp(T value, T lb, T ub)
00142 {
00143 if (value<=lb)
00144 return lb;
00145 else if (value>=ub)
00146 return ub;
00147 else
00148 return value;
00149 }
00150
00152 template <class T>
00153 static inline T abs(T a)
00154 {
00155
00156
00157 if (a==0)
00158 return 0;
00159 else if (a>0)
00160 return a;
00161 else
00162 return -a;
00163 }
00165
00168
00170 static UINT crc32(BYTE *data, INT len);
00171
00172 static inline DREAL round(DREAL d)
00173 {
00174 return ::floor(d+0.5);
00175 }
00176
00177 static inline DREAL floor(DREAL d)
00178 {
00179 return ::floor(d);
00180 }
00181
00182 static inline DREAL ceil(DREAL d)
00183 {
00184 return ::ceil(d);
00185 }
00186
00188 template <class T>
00189 static inline T sign(T a)
00190 {
00191 if (a==0)
00192 return 0;
00193 else return (a<0) ? (-1) : (+1);
00194 }
00195
00197 template <class T>
00198 static inline void swap(T & a,T &b)
00199 {
00200 T c=a;
00201 a=b;
00202 b=c;
00203 }
00204
00206 template <class T>
00207 static inline T sq(T x)
00208 {
00209 return x*x;
00210 }
00211
00213 static inline SHORTREAL sqrt(SHORTREAL x)
00214 {
00215 return ::sqrtf(x);
00216 }
00217
00219 static inline DREAL sqrt(DREAL x)
00220 {
00221 return ::sqrt(x);
00222 }
00223
00225 static inline LONGREAL sqrt(LONGREAL x)
00226 {
00227
00228
00229 #ifdef HAVE_SQRTL
00230 return ::sqrtl(x);
00231 #else
00232 return ::sqrt(x);
00233 #endif
00234 }
00235
00236
00238 static inline long double powl(long double x, long double n)
00239 {
00240
00241
00242 #ifdef HAVE_POWL
00243 return ::powl(x, n);
00244 #else
00245 return ::pow(x, n);
00246 #endif
00247 }
00248
00249 static inline INT pow(INT x, INT n)
00250 {
00251 ASSERT(n>=0);
00252 INT result=1;
00253 while (n--)
00254 result*=x;
00255
00256 return result;
00257 }
00258
00259 static inline DREAL pow(DREAL x, DREAL n)
00260 {
00261 return ::pow(x, n);
00262 }
00263
00264 static inline DREAL log10(DREAL v)
00265 {
00266 return ::log(v)/::log(10.0);
00267 }
00268
00269 #ifdef HAVE_LOG2
00270 static inline DREAL log2(DREAL v)
00271 {
00272 return ::log2(v);
00273 }
00274 #else
00275 static inline DREAL log2(DREAL v)
00276 {
00277 return ::log(v)/::log(2.0);
00278 }
00279 #endif //HAVE_LOG2
00280
00281 static inline DREAL log(DREAL v)
00282 {
00283 return ::log(v);
00284 }
00285
00286 template <class T> static void transpose_matrix(T*& matrix, INT& num_feat, INT& num_vec)
00287 {
00288 T* transposed=new T[num_vec*num_feat];
00289 for (INT i=0; i<num_vec; i++)
00290 {
00291 for (INT j=0; j<num_feat; j++)
00292 transposed[j+i*num_feat]=matrix[i+j*num_vec];
00293 }
00294
00295 delete[] matrix;
00296 matrix=transposed;
00297
00298 INT tmp=num_feat;
00299 num_feat=num_vec;
00300 num_vec=tmp;
00301 }
00302
00303 #ifdef HAVE_LAPACK
00306 static DREAL* pinv(DREAL* matrix, INT rows, INT cols, DREAL* target=NULL);
00307
00308
00309
00310
00311 static inline void dgemm(double alpha, const double* A, int rows, int cols, CBLAS_TRANSPOSE transposeA,
00312 double *B, int cols_B, CBLAS_TRANSPOSE transposeB,
00313 double beta, double *C)
00314 {
00315 cblas_dgemm(CblasColMajor, transposeA, transposeB, rows, cols, cols_B,
00316 alpha, A, cols, B, cols_B, beta, C, cols);
00317 }
00318
00319
00320 static inline void dgemv(double alpha, const double *A, int rows, int cols, const CBLAS_TRANSPOSE transposeA,
00321 const double* X, double beta, double* Y)
00322 {
00323 cblas_dgemv(CblasColMajor, transposeA,
00324 rows, cols, alpha, A, cols,
00325 X, 1, beta, Y, 1);
00326 }
00327 #endif
00328
00329 static inline LONG factorial(INT n)
00330 {
00331 LONG res=1 ;
00332 for (int i=2; i<=n; i++)
00333 res*=i ;
00334 return res ;
00335 }
00336
00337 static void init_random(UINT initseed=0)
00338 {
00339 if (initseed==0)
00340 {
00341 struct timeval tv;
00342 gettimeofday(&tv, NULL);
00343 seed=(UINT) (4223517*getpid()*tv.tv_sec*tv.tv_usec);
00344 }
00345 else
00346 seed=initseed;
00347 #if !defined(CYGWIN) && !defined(__INTERIX)
00348
00349
00350 initstate(seed, CMath::rand_state, sizeof(CMath::rand_state));
00351 #endif
00352 }
00353
00354 static inline LONG random()
00355 {
00356 #if defined(CYGWIN) || defined(__INTERIX)
00357 return rand();
00358 #else
00359 return ::random();
00360 #endif
00361 }
00362
00363 static inline INT random(INT min_value, INT max_value)
00364 {
00365 INT ret = min_value + (INT) ((max_value-min_value+1) * (random() / (RAND_MAX+1.0)));
00366 ASSERT(ret>=min_value && ret<=max_value);
00367 return ret ;
00368 }
00369
00370 static inline SHORTREAL random(SHORTREAL min_value, SHORTREAL max_value)
00371 {
00372 SHORTREAL ret = min_value + ((max_value-min_value) * (random() / (1.0*RAND_MAX)));
00373
00374 if (ret<min_value || ret>max_value)
00375 SG_SPRINT("min_value:%10.10f value: %10.10f max_value:%10.10f", min_value, ret, max_value);
00376 ASSERT(ret>=min_value && ret<=max_value);
00377 return ret;
00378 }
00379
00380 static inline DREAL random(DREAL min_value, DREAL max_value)
00381 {
00382 DREAL ret = min_value + ((max_value-min_value) * (random() / (1.0*RAND_MAX)));
00383
00384 if (ret<min_value || ret>max_value)
00385 SG_SPRINT("min_value:%10.10f value: %10.10f max_value:%10.10f", min_value, ret, max_value);
00386 ASSERT(ret>=min_value && ret<=max_value);
00387 return ret;
00388 }
00389
00390 template <class T>
00391 static void fill_vector(T* vec, INT len, T value)
00392 {
00393 for (INT i=0; i<len; i++)
00394 vec[i]=value;
00395 }
00396 template <class T>
00397 static void range_fill_vector(T* vec, INT len, T start=0)
00398 {
00399 for (INT i=0; i<len; i++)
00400 vec[i]=i+start;
00401 }
00402
00403 template <class T>
00404 static void random_vector(T* vec, INT len, T min_value, T max_value)
00405 {
00406 for (INT i=0; i<len; i++)
00407 vec[i]=CMath::random(min_value, max_value);
00408 }
00409
00410 static inline INT* randperm(INT n)
00411 {
00412 INT* perm = new INT[n];
00413
00414 if (!perm)
00415 return NULL;
00416 for (INT i = 0; i < n; i++)
00417 perm[i] = i;
00418 for (INT i = 0; i < n; i++)
00419 swap(perm[random(0, n - 1)], perm[i]);
00420 return perm;
00421 }
00422
00423 static inline LONG nchoosek(INT n, INT k)
00424 {
00425 long res=1 ;
00426
00427 for (INT i=n-k+1; i<=n; i++)
00428 res*=i ;
00429
00430 return res/factorial(k) ;
00431 }
00432
00434 template <class T>
00435 static inline void vec1_plus_scalar_times_vec2(T* vec1,
00436 T scalar, const T* vec2, INT n)
00437 {
00438 for (INT i=0; i<n; i++)
00439 vec1[i]+=scalar*vec2[i];
00440 }
00441
00443 static inline DREAL dot(const DREAL* v1, const DREAL* v2, INT n)
00444 {
00445 DREAL r=0;
00446 #ifdef HAVE_LAPACK
00447 INT skip=1;
00448 r = cblas_ddot(n, v1, skip, v2, skip);
00449 #else
00450 for (INT i=0; i<n; i++)
00451 r+=v1[i]*v2[i];
00452 #endif
00453 return r;
00454 }
00455
00457 static inline SHORTREAL dot(const SHORTREAL* v1, const SHORTREAL* v2, INT n)
00458 {
00459 DREAL r=0;
00460 #ifdef HAVE_LAPACK
00461 INT skip=1;
00462 r = cblas_sdot(n, v1, skip, v2, skip);
00463 #else
00464 for (INT i=0; i<n; i++)
00465 r+=v1[i]*v2[i];
00466 #endif
00467 return r;
00468 }
00469
00471 template <class T>
00472 static inline void add(T* target, T alpha, const T* v1, T beta, const T* v2, INT len)
00473 {
00474 for (INT i=0; i<len; i++)
00475 target[i]=alpha*v1[i]+beta*v2[i];
00476 }
00477
00479 template <class T>
00480 static inline void add_scalar(T alpha, T* vec, INT len)
00481 {
00482 for (INT i=0; i<len; i++)
00483 vec[i]+=alpha;
00484 }
00485
00487 template <class T>
00488 static inline void scale_vector(T alpha, T* vec, INT len)
00489 {
00490 for (INT i=0; i<len; i++)
00491 vec[i]*=alpha;
00492 }
00493
00495 template <class T>
00496 static inline T sum(T* vec, INT len)
00497 {
00498 T result=0;
00499 for (INT i=0; i<len; i++)
00500 result+=vec[i];
00501
00502 return result;
00503 }
00504
00506 template <class T>
00507 static inline T max(T* vec, INT len)
00508 {
00509 ASSERT(len>0);
00510 T maxv=vec[0];
00511
00512 for (INT i=1; i<len; i++)
00513 maxv=CMath::max(vec[i], maxv);
00514
00515 return maxv;
00516 }
00517
00519 template <class T>
00520 static inline T sum_abs(T* vec, INT len)
00521 {
00522 T result=0;
00523 for (INT i=0; i<len; i++)
00524 result+=CMath::abs(vec[i]);
00525
00526 return result;
00527 }
00528
00529 static inline DREAL mean(DREAL* vec, INT len)
00530 {
00531 ASSERT(vec);
00532 ASSERT(len>0);
00533
00534 DREAL mean=0;
00535 for (INT i=0; i<len; i++)
00536 mean+=vec[i];
00537 return mean/len;
00538 }
00539
00540 static inline DREAL trace(DREAL* mat, INT cols, INT rows)
00541 {
00542 DREAL trace=0;
00543 for (INT i=0; i<rows; i++)
00544 trace+=mat[i*cols+i];
00545 return trace;
00546 }
00547
00551 static void sort(INT *a, INT cols, INT sort_col=0) ;
00552 static void sort(DREAL *a, INT*idx, INT N) ;
00553
00554
00555
00556
00557
00558
00559
00561 template <class T>
00562 inline static void radix_sort(T* array, size_t size)
00563 {
00564 radix_sort_helper(array,size,0);
00565 }
00566
00567 template <class T>
00568 static inline BYTE byte(T word, unsigned short p)
00569 {
00570 return (word >> (sizeof(T)-p-1) * 8) & 0xff;
00571 }
00572
00573 template <class T>
00574 static void radix_sort_helper(T* array, size_t size, unsigned short i)
00575 {
00576 static size_t count[256], nc, cmin;
00577 T *ak;
00578 BYTE c=0;
00579 radix_stack_t<T> s[RADIX_STACK_SIZE], *sp, *olds, *bigs;
00580 T *an, *aj, *pile[256];
00581 size_t *cp, cmax;
00582
00583
00584 sp = s;
00585 radix_push(array, size, i);
00586
00587
00588 while (sp>s) {
00589 radix_pop(array, size, i);
00590 an = array + size;
00591
00592
00593 if (nc == 0) {
00594 cmin = 0xff;
00595 for (ak = array; ak < an; ak++) {
00596 c = byte(*ak, i);
00597 count[c]++;
00598 if (count[c] == 1) {
00599
00600 if (c < cmin)
00601 cmin = c;
00602 nc++;
00603 }
00604 }
00605
00606
00607 if (sp + nc > s + RADIX_STACK_SIZE) {
00608 radix_sort_helper(array, size, i);
00609 continue;
00610 }
00611 }
00612
00613
00614
00615
00616
00617
00618
00619 olds = bigs = sp;
00620 cmax = 2;
00621 ak = array;
00622 pile[0xff] = an;
00623 for (cp = count + cmin; nc > 0; cp++) {
00624
00625 while (*cp == 0)
00626 cp++;
00627
00628 if (*cp > 1) {
00629
00630 if (*cp > cmax) {
00631 cmax = *cp;
00632 bigs = sp;
00633 }
00634
00635 if (i < sizeof(T)-1)
00636 radix_push(ak, *cp, (unsigned short) (i + 1));
00637 }
00638 pile[cp - count] = ak += *cp;
00639 nc--;
00640 }
00641
00642
00643 swap(*olds, *bigs);
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657 aj = array;
00658 while (aj<an)
00659 {
00660 T r;
00661
00662 for (r = *aj; aj < (ak = --pile[c = byte(r, i)]);)
00663 swap(*ak, r);
00664
00665 *aj = r;
00666 aj += count[c];
00667 count[c] = 0;
00668 }
00669 }
00670 }
00673 template <class T>
00674 static void insertion_sort(T* output, INT size)
00675 {
00676 for (INT i=0; i<size-1; i++)
00677 {
00678 INT j=i-1;
00679 T value=output[i];
00680 while (j >= 0 && output[j] > value)
00681 {
00682 output[j+1] = output[j];
00683 j--;
00684 }
00685 output[j+1]=value;
00686 }
00687 }
00688
00689
00692
00693
00694 template <class T>
00695 static void qsort(T* output, INT size)
00696 {
00697 if (size==2)
00698 {
00699 if (output[0] > output [1])
00700 swap(output[0],output[1]);
00701 return;
00702 }
00703 T split=output[random(0,size-1)];
00704
00705 INT left=0;
00706 INT right=size-1;
00707
00708 while (left<=right)
00709 {
00710 while (output[left] < split)
00711 left++;
00712 while (output[right] > split)
00713 right--;
00714
00715 if (left<=right)
00716 {
00717 swap(output[left],output[right]);
00718 left++;
00719 right--;
00720 }
00721 }
00722
00723 if (right+1> 1)
00724 qsort(output,right+1);
00725
00726 if (size-left> 1)
00727 qsort(&output[left],size-left);
00728 }
00729
00731 template <class T> static void display_vector(T* vector, INT n, const char* name="vector");
00732
00734 template <class T> static void display_matrix(T* matrix, INT rows, INT cols, const char* name="matrix");
00735
00741 template <class T1,class T2>
00742 static void qsort_index(T1* output, T2* index, UINT size);
00743
00744 template <class T1,class T2>
00745 static void* parallel_qsort_index(void* p);
00746
00752 template <class T1,class T2>
00753 static void qsort_backward_index(T1* output, T2* index, INT size);
00754
00755
00756
00757 template <class T>
00758 static void min(DREAL* output, T* index, INT size);
00759
00760
00761
00762 template <class T>
00763 static void nmin(DREAL* output, T* index, INT size, INT n);
00764
00765
00766
00767 template <class T>
00768 static INT unique(T* output, INT size)
00769 {
00770 qsort(output, size);
00771 INT i,j=0 ;
00772 for (i=0; i<size; i++)
00773 if (i==0 || output[i]!=output[i-1])
00774 output[j++]=output[i];
00775 return j ;
00776 }
00777
00778
00779
00780 template <class T>
00781 static INT binary_search_helper(T* output, INT size, T elem)
00782 {
00783 INT start=0;
00784 INT end=size-1;
00785
00786 if (size<1)
00787 return -1;
00788
00789 while (start<end)
00790 {
00791 INT middle=(start+end)/2;
00792
00793 if (output[middle]>elem)
00794 end=middle-1;
00795 else if (output[middle]<elem)
00796 start=middle+1;
00797 else
00798 return middle;
00799 }
00800
00801 if (output[start]==elem)
00802 return start;
00803 else
00804 return -1;
00805 }
00806
00807
00808
00809 template <class T>
00810 static inline INT binary_search(T* output, INT size, T elem)
00811 {
00812 INT ind = binary_search_helper(output, size, elem);
00813 if (ind >= 0 && output[ind] == elem)
00814 return ind;
00815 return -1;
00816 }
00817
00818
00819
00820
00821
00822 template <class T>
00823 static INT binary_search_max_lower_equal(T* output, INT size, T elem)
00824 {
00825 INT ind = binary_search_helper(output, size, elem);
00826
00827 if (output[ind]<=elem)
00828 return ind;
00829 if (ind>0 && output[ind-1] <= elem)
00830 return ind-1;
00831 return -1;
00832 }
00833
00836 static DREAL Align(CHAR * seq1, CHAR* seq2, INT l1, INT l2, DREAL gapCost);
00837
00842 static INT calcroc(DREAL* fp, DREAL* tp, DREAL* output, INT* label, INT& size, INT& possize, INT& negsize, DREAL& tresh, FILE* rocfile);
00844
00847 static double mutual_info(DREAL* p1, DREAL* p2, INT len);
00848
00851 static double relative_entropy(DREAL* p, DREAL* q, INT len);
00852
00854 static double entropy(DREAL* p, INT len);
00855
00857 inline static UINT get_seed()
00858 {
00859 return CMath::seed;
00860 }
00861
00862
00873 #ifdef USE_LOGCACHE
00874 static inline DREAL logarithmic_sum(DREAL p, DREAL q)
00875 {
00876 DREAL diff;
00877
00878 if (!finite(p))
00879 return q;
00880 if (!finite(q))
00881 {
00882 SG_SWARNING(("INVALID second operand to logsum(%f,%f) expect undefined results\n", p, q);
00883 return NAN;
00884 }
00885 diff = p - q;
00886 if (diff > 0)
00887 return diff > LOGRANGE? p : p + logtable[(int)(diff * LOGACCURACY)];
00888 return -diff > LOGRANGE? q : q + logtable[(int)(-diff * LOGACCURACY)];
00889 }
00890
00892 static void init_log_table();
00893
00895 static INT determine_logrange();
00896
00898 static INT determine_logaccuracy(INT range);
00899 #else
00900 static inline DREAL logarithmic_sum(DREAL p, DREAL q)
00901 {
00902 DREAL diff;
00903
00904 if (!finite(p))
00905 return q;
00906 if (!finite(q))
00907 return p;
00908 diff = p - q;
00909 if (diff > 0)
00910 return diff > LOGRANGE? p : p + log(1 + exp(-diff));
00911 return -diff > LOGRANGE? q : q + log(1 + exp(diff));
00912 }
00913 #endif
00914 #ifdef LOG_SUM_ARRAY
00915
00920 static inline DREAL logarithmic_sum_array(DREAL *p, INT len)
00921 {
00922 if (len<=2)
00923 {
00924 if (len==2)
00925 return logarithmic_sum(p[0],p[1]) ;
00926 if (len==1)
00927 return p[0];
00928 return -INFTY ;
00929 }
00930 else
00931 {
00932 register DREAL *pp=p ;
00933 if (len%2==1) pp++ ;
00934 for (register INT j=0; j < len>>1; j++)
00935 pp[j]=logarithmic_sum(pp[j<<1], pp[1+(j<<1)]) ;
00936 }
00937 return logarithmic_sum_array(p,len%2+len>>1) ;
00938 }
00939 #endif
00940
00941 public:
00944
00945 static const DREAL INFTY;
00946 static const DREAL ALMOST_INFTY;
00947
00949 static const DREAL ALMOST_NEG_INFTY;
00950
00952 static INT LOGRANGE;
00953
00955 static UINT seed;
00956
00957 #ifdef USE_LOGCACHE
00958
00960 static INT LOGACCURACY;
00962 protected:
00964 static DREAL* logtable;
00965 #endif
00966 static CHAR* rand_state;
00967 };
00968
00969 template <class T1,class T2>
00970 void* CMath::parallel_qsort_index(void* p)
00971 {
00972 struct thread_qsort* ps=(thread_qsort*) p;
00973 T1* output=ps->output;
00974 T2* index=ps->index;
00975 INT size=ps->size;
00976 INT* qsort_threads=ps->qsort_threads;
00977 INT sort_limit=ps->sort_limit;
00978 INT num_threads=ps->num_threads;
00979
00980 if (size==2)
00981 {
00982 if (output[0] > output [1])
00983 {
00984 swap(output[0], output[1]);
00985 swap(index[0], index[1]);
00986 }
00987 return NULL;
00988 }
00989
00990 double split=output[size/2];
00991
00992 INT left=0;
00993 INT right=size-1;
00994
00995 while (left<=right)
00996 {
00997 while (output[left] < split)
00998 left++;
00999 while (output[right] > split)
01000 right--;
01001
01002 if (left<=right)
01003 {
01004 swap(output[left], output[right]);
01005 swap(index[left], index[right]);
01006 left++;
01007 right--;
01008 }
01009 }
01010 bool lthread_start=false;
01011 bool rthread_start=false;
01012 pthread_t lthread;
01013 pthread_t rthread;
01014 struct thread_qsort t1;
01015 struct thread_qsort t2;
01016
01017 if (right+1> 1 && (right+1< sort_limit || *qsort_threads >= num_threads-1))
01018 qsort_index(output,index,right+1);
01019 else if (right+1> 1)
01020 {
01021 *qsort_threads++;
01022 lthread_start=true;
01023 t1.output=output;
01024 t1.index=index;
01025 t1.size=right+1;
01026 if (pthread_create(<hread, NULL, parallel_qsort_index<T1,T2>, (void*) &t1) != 0)
01027 {
01028 lthread_start=false;
01029 *qsort_threads--;
01030 qsort_index(output,index,right+1);
01031 }
01032 }
01033
01034
01035 if (size-left> 1 && (size-left< sort_limit || *qsort_threads >= num_threads-1))
01036 qsort_index(&output[left],&index[left], size-left);
01037 else if (size-left> 1)
01038 {
01039 *qsort_threads++;
01040 rthread_start=true;
01041 t2.output=&output[left];
01042 t2.index=&index[left];
01043 t2.size=size-left;
01044 if (pthread_create(&rthread, NULL, parallel_qsort_index<T1,T2>, (void*)&t2) != 0)
01045 {
01046 rthread_start=false;
01047 *qsort_threads--;
01048 qsort_index(&output[left],&index[left], size-left);
01049 }
01050 }
01051
01052 if (lthread_start)
01053 {
01054 pthread_join(lthread, NULL);
01055 *qsort_threads--;
01056 }
01057
01058 if (rthread_start)
01059 {
01060 pthread_join(rthread, NULL);
01061 *qsort_threads--;
01062 }
01063
01064 return NULL;
01065 }
01066
01067
01068
01069 template <class T1,class T2>
01070 void CMath::qsort_index(T1* output, T2* index, UINT size)
01071 {
01072 if (size==2)
01073 {
01074 if (output[0] > output [1])
01075 {
01076 swap(output[0],output[1]);
01077 swap(index[0],index[1]);
01078 }
01079 return;
01080 }
01081 T1 split=output[random(0,size-1)];
01082
01083 INT left=0;
01084 INT right=size-1;
01085
01086 while (left<=right)
01087 {
01088 while (output[left] < split)
01089 left++;
01090 while (output[right] > split)
01091 right--;
01092
01093 if (left<=right)
01094 {
01095 swap(output[left],output[right]);
01096 swap(index[left],index[right]);
01097 left++;
01098 right--;
01099 }
01100 }
01101
01102 if (right+1> 1)
01103 qsort_index(output,index,right+1);
01104
01105 if (size-left> 1)
01106 qsort_index(&output[left],&index[left], size-left);
01107 }
01108
01109 template <class T1,class T2>
01110 void CMath::qsort_backward_index(T1* output, T2* index, INT size)
01111 {
01112 if (size==2)
01113 {
01114 if (output[0] < output [1])
01115 {
01116 swap(output[0],output[1]);
01117 swap(index[0],index[1]);
01118 }
01119 return;
01120 }
01121
01122 T1 split=output[random(0,size-1)];
01123
01124 INT left=0;
01125 INT right=size-1;
01126
01127 while (left<=right)
01128 {
01129 while (output[left] > split)
01130 left++;
01131 while (output[right] < split)
01132 right--;
01133
01134 if (left<=right)
01135 {
01136 swap(output[left],output[right]);
01137 swap(index[left],index[right]);
01138 left++;
01139 right--;
01140 }
01141 }
01142
01143 if (right+1> 1)
01144 qsort_backward_index(output,index,right+1);
01145
01146 if (size-left> 1)
01147 qsort_backward_index(&output[left],&index[left], size-left);
01148 }
01149
01150 template <class T>
01151 void CMath::nmin(DREAL* output, T* index, INT size, INT n)
01152 {
01153 if (6*n*size<13*size*CMath::log(size))
01154 for (INT i=0; i<n; i++)
01155 min(&output[i], &index[i], size-i) ;
01156 else
01157 qsort_index(output, index, size) ;
01158 }
01159
01160
01161 template <class T>
01162 void CMath::min(DREAL* output, T* index, INT size)
01163 {
01164 if (size<=1)
01165 return;
01166 DREAL min_elem=output[0];
01167 INT min_index=0;
01168 for (INT i=1; i<size; i++)
01169 {
01170 if (output[i]<min_elem)
01171 {
01172 min_index=i;
01173 min_elem=output[i];
01174 }
01175 }
01176 swap(output[0], output[min_index]);
01177 swap(index[0], index[min_index]);
01178 }
01179 #endif