Mathematics.cpp

Go to the documentation of this file.
00001 /*
00002  * This program is free software; you can redistribute it and/or modify
00003  * it under the terms of the GNU General Public License as published by
00004  * the Free Software Foundation; either version 3 of the License, or
00005  * (at your option) any later version.
00006  *
00007  * Written (W) 1999-2009 Soeren Sonnenburg
00008  * Written (W) 1999-2008 Gunnar Raetsch
00009  * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society
00010  */
00011 
00012 // Math.cpp: implementation of the CMath class.
00013 //
00015 
00016 
00017 #include "lib/common.h"
00018 #include "lib/Mathematics.h"
00019 #include "lib/lapack.h"
00020 #include "lib/io.h"
00021 
00022 #include <stdio.h>
00023 #include <stdlib.h>
00024 #include <math.h>
00025 
00026 using namespace shogun;
00027 
00028 #ifdef USE_LOGCACHE
00029 #ifdef USE_HMMDEBUG
00030 #define MAX_LOG_TABLE_SIZE 10*1024*1024
00031 #define LOG_TABLE_PRECISION 1e-6
00032 #else //USE_HMMDEBUG
00033 #define MAX_LOG_TABLE_SIZE 123*1024*1024
00034 #define LOG_TABLE_PRECISION 1e-15
00035 #endif //USE_HMMDEBUG
00036 int32_t CMath::LOGACCURACY         = 0; // 100000 steps per integer
00037 #endif // USE_LOGCACHE
00038 
00039 int32_t CMath::LOGRANGE            = 0; // range for logtable: log(1+exp(x))  -25 <= x <= 0
00040 
00041 const float64_t CMath::INFTY            =  -log(0.0);   // infinity
00042 const float64_t CMath::ALMOST_INFTY     =  +1e+20;      //a large number
00043 const float64_t CMath::ALMOST_NEG_INFTY =  -1000;   
00044 #ifdef USE_LOGCACHE
00045 float64_t* CMath::logtable = NULL;
00046 #endif
00047 char* CMath::rand_state = NULL;
00048 uint32_t CMath::seed = 0;
00049 
00050 CMath::CMath()
00051 : CSGObject()
00052 {
00053     CMath::rand_state=new char[RNG_SEED_SIZE];
00054     init_random();
00055 
00056 #ifdef USE_LOGCACHE
00057     LOGRANGE=CMath::determine_logrange();
00058     LOGACCURACY=CMath::determine_logaccuracy(LOGRANGE);
00059     CMath::logtable=new float64_t[LOGRANGE*LOGACCURACY];
00060     init_log_table();
00061 #else
00062     int32_t i=0;
00063     while ((float64_t)log(1+((float64_t)exp(-float64_t(i)))))
00064         i++;
00065 
00066     LOGRANGE=i;
00067 #endif 
00068 }
00069 
00070 CMath::~CMath()
00071 {
00072     delete[] CMath::rand_state;
00073     CMath::rand_state=NULL;
00074 #ifdef USE_LOGCACHE
00075     delete[] CMath::logtable;
00076     CMath::logtable=NULL;
00077 #endif
00078 }
00079 
00080 #ifdef USE_LOGCACHE
00081 int32_t CMath::determine_logrange()
00082 {
00083     int32_t i;
00084     float64_t acc=0;
00085     for (i=0; i<50; i++)
00086     {
00087     acc=((float64_t)log(1+((float64_t)exp(-float64_t(i)))));
00088     if (acc<=(float64_t)LOG_TABLE_PRECISION)
00089         break;
00090     }
00091 
00092     SG_SINFO( "determined range for x in table log(1+exp(-x)) is:%d (error:%G)\n",i,acc);
00093     return i;
00094 }
00095 
00096 int32_t CMath::determine_logaccuracy(int32_t range)
00097 {
00098     range=MAX_LOG_TABLE_SIZE/range/((int)sizeof(float64_t));
00099     SG_SINFO( "determined accuracy for x in table log(1+exp(-x)) is:%d (error:%G)\n",range,1.0/(double) range);
00100     return range;
00101 }
00102 
00103 //init log table of form log(1+exp(x))
00104 void CMath::init_log_table()
00105 {
00106   for (int32_t i=0; i< LOGACCURACY*LOGRANGE; i++)
00107     logtable[i]=log(1+exp(float64_t(-i)/float64_t(LOGACCURACY)));
00108 }
00109 #endif
00110 
00111 void CMath::sort(int32_t *a, int32_t cols, int32_t sort_col)
00112 {
00113   int32_t changed=1;
00114   if (a[0]==-1) return ;
00115   while (changed)
00116   {
00117       changed=0; int32_t i=0 ;
00118       while ((a[(i+1)*cols]!=-1) && (a[(i+1)*cols+1]!=-1)) // to be sure
00119       {
00120           if (a[i*cols+sort_col]>a[(i+1)*cols+sort_col])
00121           {
00122               for (int32_t j=0; j<cols; j++)
00123                   CMath::swap(a[i*cols+j],a[(i+1)*cols+j]) ;
00124               changed=1 ;
00125           } ;
00126           i++ ;
00127       } ;
00128   } ;
00129 } 
00130 
00131 void CMath::sort(float64_t *a, int32_t* idx, int32_t N) 
00132 {
00133     int32_t changed=1;
00134     while (changed)
00135     {
00136         changed=0;
00137         for (int32_t i=0; i<N-1; i++)
00138         {
00139             if (a[i]>a[i+1])
00140             {
00141                 swap(a[i],a[i+1]) ;
00142                 swap(idx[i],idx[i+1]) ;
00143                 changed=1 ;
00144             } ;
00145         } ;
00146     } ;
00147  
00148 }
00149 
00150 float64_t CMath::Align(
00151     char* seq1, char* seq2, int32_t l1, int32_t l2, float64_t gapCost)
00152 {
00153   float64_t actCost=0 ;
00154   int32_t i1, i2 ;
00155   float64_t* const gapCosts1 = new float64_t[ l1 ];
00156   float64_t* const gapCosts2 = new float64_t[ l2 ];
00157   float64_t* costs2_0 = new float64_t[ l2 + 1 ];
00158   float64_t* costs2_1 = new float64_t[ l2 + 1 ];
00159 
00160   // initialize borders
00161   for( i1 = 0; i1 < l1; ++i1 ) {
00162     gapCosts1[ i1 ] = gapCost * i1;
00163   }
00164   costs2_1[ 0 ] = 0;
00165   for( i2 = 0; i2 < l2; ++i2 ) {
00166     gapCosts2[ i2 ] = gapCost * i2;
00167     costs2_1[ i2+1 ] = costs2_1[ i2 ] + gapCosts2[ i2 ];
00168   }
00169   // compute alignment
00170   for( i1 = 0; i1 < l1; ++i1 ) {
00171     swap( costs2_0, costs2_1 );
00172     actCost = costs2_0[ 0 ] + gapCosts1[ i1 ];
00173     costs2_1[ 0 ] = actCost;
00174     for( i2 = 0; i2 < l2; ++i2 ) {
00175       const float64_t actMatch = costs2_0[ i2 ] + ( seq1[i1] == seq2[i2] );
00176       const float64_t actGap1 = costs2_0[ i2+1 ] + gapCosts1[ i1 ];
00177       const float64_t actGap2 = actCost + gapCosts2[ i2 ];
00178       const float64_t actGap = min( actGap1, actGap2 );
00179       actCost = min( actMatch, actGap );
00180       costs2_1[ i2+1 ] = actCost;
00181     }
00182   }
00183 
00184   delete [] gapCosts1;
00185   delete [] gapCosts2;
00186   delete [] costs2_0;
00187   delete [] costs2_1;
00188   
00189   // return the final cost
00190   return actCost;
00191 }
00192 
00193 //plot x- axis false positives (fp) 1-Specificity
00194 //plot y- axis true positives (tp) Sensitivity
00195 int32_t CMath::calcroc(
00196     float64_t* fp, float64_t* tp, float64_t* output, int32_t* label,
00197     int32_t& size, int32_t& possize, int32_t& negsize, float64_t& tresh,
00198     FILE* rocfile)
00199 {
00200     int32_t left=0;
00201     int32_t right=size-1;
00202     int32_t i;
00203 
00204     for (i=0; i<size; i++)
00205     {
00206         if (!(label[i]== -1 || label[i]==1))
00207             return -1;
00208     }
00209 
00210     //sort data such that -1 labels first +1 behind
00211     while (left<right)
00212     {
00213         while ((label[left] < 0) && (left<right))
00214             left++;
00215         while ((label[right] > 0) && (left<right))  //warning: label must be either -1 or +1
00216             right--;
00217 
00218         swap(output[left],output[right]);
00219         swap(label[left], label[right]);
00220     }
00221 
00222     // neg/pos sizes
00223     negsize=left;
00224     possize=size-left;
00225     float64_t* negout=output;
00226     float64_t* posout=&output[left];
00227 
00228     // sort the pos and neg outputs separately
00229     qsort(negout, negsize);
00230     qsort(posout, possize);
00231 
00232     // set minimum+maximum values for decision-treshhold
00233     float64_t minimum=min(negout[0], posout[0]);
00234     float64_t maximum=minimum;
00235     if (negsize>0)
00236         maximum=max(maximum,negout[negsize-1]);
00237     if (possize>0)
00238         maximum=max(maximum,posout[possize-1]);
00239 
00240     float64_t treshhold=minimum;
00241     float64_t old_treshhold=treshhold;
00242 
00243     //clear array.
00244     for (i=0; i<size; i++)
00245     {
00246         fp[i]=1.0;
00247         tp[i]=1.0;
00248     }
00249 
00250     //start with fp=1.0 tp=1.0 which is posidx=0, negidx=0
00251     //everything right of {pos,neg}idx is considered to belong to +1
00252     int32_t posidx=0;
00253     int32_t negidx=0;
00254     int32_t iteration=1;
00255     int32_t returnidx=-1;
00256 
00257     float64_t minerr=10;
00258 
00259     while (iteration < size && treshhold<=maximum)
00260     {
00261         old_treshhold=treshhold;
00262 
00263         while (treshhold==old_treshhold && treshhold<=maximum)
00264         {
00265             if (posidx<possize && negidx<negsize)
00266             {
00267                 if (posout[posidx]<negout[negidx])
00268                 {
00269                     if (posout[posidx]==treshhold)
00270                         posidx++;
00271                     else
00272                         treshhold=posout[posidx];
00273                 }
00274                 else
00275                 {
00276                     if (negout[negidx]==treshhold)
00277                         negidx++;
00278                     else
00279                         treshhold=negout[negidx];
00280                 }
00281             }
00282             else
00283             {
00284                 if (posidx>=possize && negidx<negsize-1)
00285                 {
00286                     negidx++;
00287                     treshhold=negout[negidx];
00288                 }
00289                 else if (negidx>=negsize && posidx<possize-1)
00290                 {
00291                     posidx++;
00292                     treshhold=posout[posidx];
00293                 }
00294                 else if (negidx<negsize && treshhold!=negout[negidx])
00295                     treshhold=negout[negidx];
00296                 else if (posidx<possize && treshhold!=posout[posidx])
00297                     treshhold=posout[posidx];
00298                 else
00299                 {
00300                     treshhold=2*(maximum+100); // force termination
00301                     posidx=possize;
00302                     negidx=negsize;
00303                     break;
00304                 }
00305             }
00306         }
00307 
00308         //calc tp,fp
00309         tp[iteration]=(possize-posidx)/(float64_t (possize));
00310         fp[iteration]=(negsize-negidx)/(float64_t (negsize));
00311 
00312         //choose point with minimal error
00313         if (minerr > negsize*fp[iteration]/size+(1-tp[iteration])*possize/size )
00314         {
00315             minerr=negsize*fp[iteration]/size+(1-tp[iteration])*possize/size;
00316             tresh=(old_treshhold+treshhold)/2;
00317             returnidx=iteration;
00318         }
00319 
00320         iteration++;
00321     }
00322 
00323     // set new size
00324     size=iteration;
00325 
00326     if (rocfile)
00327     {
00328         const char id[]="ROC";
00329         fwrite(id, sizeof(char), sizeof(id), rocfile);
00330         fwrite(fp, sizeof(float64_t), size, rocfile);
00331         fwrite(tp, sizeof(float64_t), size, rocfile);
00332     }
00333 
00334     return returnidx;
00335 }
00336 
00337 float64_t CMath::mutual_info(float64_t* p1, float64_t* p2, int32_t len)
00338 {
00339     double e=0;
00340 
00341     for (int32_t i=0; i<len; i++)
00342         for (int32_t j=0; j<len; j++)
00343             e+=exp(p2[j*len+i])*(p2[j*len+i]-p1[i]-p1[j]);
00344 
00345     return (float64_t) e;
00346 }
00347 
00348 float64_t CMath::relative_entropy(float64_t* p, float64_t* q, int32_t len)
00349 {
00350     double e=0;
00351 
00352     for (int32_t i=0; i<len; i++)
00353         e+=exp(p[i])*(p[i]-q[i]);
00354 
00355     return (float64_t) e;
00356 }
00357 
00358 float64_t CMath::entropy(float64_t* p, int32_t len)
00359 {
00360     double e=0;
00361 
00362     for (int32_t i=0; i<len; i++)
00363         e-=exp(p[i])*p[i];
00364 
00365     return (float64_t) e;
00366 }
00367 
00368 
00369 //Howto construct the pseudo inverse (from "The Matrix Cookbook")
00370 //
00371 //Assume A does not have full rank, i.e. A is n \times m and rank(A) = r < min(n;m).
00372 //
00373 //The matrix A+ known as the pseudo inverse is unique and does always exist.
00374 //
00375 //The pseudo inverse A+ can be constructed from the singular value
00376 //decomposition A = UDV^T , by  A^+ = V(D+)U^T.
00377 
00378 #ifdef HAVE_LAPACK
00379 float64_t* CMath::pinv(
00380     float64_t* matrix, int32_t rows, int32_t cols, float64_t* target)
00381 {
00382     if (!target)
00383         target=new float64_t[rows*cols];
00384 
00385     char jobu='A';
00386     char jobvt='A';
00387     int m=rows; /* for calling external lib */
00388     int n=cols; /* for calling external lib */
00389     int lda=m; /* for calling external lib */
00390     int ldu=m; /* for calling external lib */
00391     int ldvt=n; /* for calling external lib */
00392     int info=-1; /* for calling external lib */
00393     int32_t lsize=CMath::min((int32_t) m, (int32_t) n);
00394     double* s=new double[lsize];
00395     double* u=new double[m*m];
00396     double* vt=new double[n*n];
00397 
00398     wrap_dgesvd(jobu, jobvt, m, n, matrix, lda, s, u, ldu, vt, ldvt, &info);
00399     ASSERT(info==0);
00400 
00401     for (int32_t i=0; i<n; i++)
00402     {
00403         for (int32_t j=0; j<lsize; j++)
00404             vt[i*n+j]=vt[i*n+j]/s[j];
00405     }
00406 
00407     cblas_dgemm(CblasColMajor, CblasTrans, CblasTrans, m, n, m, 1.0, vt, ldvt, u, ldu, 0, target, m);
00408 
00409     delete[] u;
00410     delete[] vt;
00411     delete[] s;
00412 
00413     return target;
00414 }
00415 #endif
00416 
00417 namespace shogun
00418 {
00419 template <>
00420 void CMath::display_vector(const uint8_t* vector, int32_t n, const char* name)
00421 {
00422     ASSERT(n>=0);
00423     SG_SPRINT("%s=[", name);
00424     for (int32_t i=0; i<n; i++)
00425         SG_SPRINT("%d%s", vector[i], i==n-1? "" : ",");
00426     SG_SPRINT("]\n");
00427 }
00428 
00429 template <>
00430 void CMath::display_vector(const int32_t* vector, int32_t n, const char* name)
00431 {
00432     ASSERT(n>=0);
00433     SG_SPRINT("%s=[", name);
00434     for (int32_t i=0; i<n; i++)
00435         SG_SPRINT("%d%s", vector[i], i==n-1? "" : ",");
00436     SG_SPRINT("]\n");
00437 }
00438 
00439 template <>
00440 void CMath::display_vector(const int64_t* vector, int32_t n, const char* name)
00441 {
00442     ASSERT(n>=0);
00443     SG_SPRINT("%s=[", name);
00444     for (int32_t i=0; i<n; i++)
00445         SG_SPRINT("%lld%s", vector[i], i==n-1? "" : ",");
00446     SG_SPRINT("]\n");
00447 }
00448 
00449 template <>
00450 void CMath::display_vector(const uint64_t* vector, int32_t n, const char* name)
00451 {
00452     ASSERT(n>=0);
00453     SG_SPRINT("%s=[", name);
00454     for (int32_t i=0; i<n; i++)
00455         SG_SPRINT("%llu%s", vector[i], i==n-1? "" : ",");
00456     SG_SPRINT("]\n");
00457 }
00458 
00459 template <>
00460 void CMath::display_vector(const float32_t* vector, int32_t n, const char* name)
00461 {
00462     ASSERT(n>=0);
00463     SG_SPRINT("%s=[", name);
00464     for (int32_t i=0; i<n; i++)
00465         SG_SPRINT("%10.10f%s", vector[i], i==n-1? "" : ",");
00466     SG_SPRINT("]\n");
00467 }
00468 
00469 template <>
00470 void CMath::display_vector(const float64_t* vector, int32_t n, const char* name)
00471 {
00472     ASSERT(n>=0);
00473     SG_SPRINT("%s=[", name);
00474     for (int32_t i=0; i<n; i++)
00475         SG_SPRINT("%10.10f%s", vector[i], i==n-1? "" : ",");
00476     SG_SPRINT("]\n");
00477 }
00478 
00479 template <>
00480 void CMath::display_matrix(
00481     const int32_t* matrix, int32_t rows, int32_t cols, const char* name)
00482 {
00483     ASSERT(rows>=0 && cols>=0);
00484     SG_SPRINT("%s=[\n", name);
00485     for (int32_t i=0; i<rows; i++)
00486     {
00487         SG_SPRINT("[");
00488         for (int32_t j=0; j<cols; j++)
00489             SG_SPRINT("\t%d%s", matrix[j*rows+i],
00490                 j==cols-1? "" : ",");
00491         SG_SPRINT("]%s\n", i==rows-1? "" : ",");
00492     }
00493     SG_SPRINT("]\n");
00494 }
00495 
00496 template <>
00497 void CMath::display_matrix(
00498     const float64_t* matrix, int32_t rows, int32_t cols, const char* name)
00499 {
00500     ASSERT(rows>=0 && cols>=0);
00501     SG_SPRINT("%s=[\n", name);
00502     for (int32_t i=0; i<rows; i++)
00503     {
00504         SG_SPRINT("[");
00505         for (int32_t j=0; j<cols; j++)
00506             SG_SPRINT("\t%lf%s", (double) matrix[j*rows+i],
00507                 j==cols-1? "" : ",");
00508         SG_SPRINT("]%s\n", i==rows-1? "" : ",");
00509     }
00510     SG_SPRINT("]\n");
00511 }
00512 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation