Marsyas
0.6.0-alpha
|
00001 /* 00002 ** Copyright (C) 1998-2006 George Tzanetakis <gtzan@cs.uvic.ca> 00003 ** 00004 ** This program is free software; you can redistribute it and/or modify 00005 ** it under the terms of the GNU General Public License as published by 00006 ** the Free Software Foundation; either version 2 of the License, or 00007 ** (at your option) any later version. 00008 ** 00009 ** This program is distributed in the hope that it will be useful, 00010 ** but WITHOUT ANY WARRANTY; without even the implied warranty of 00011 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00012 ** GNU General Public License for more details. 00013 ** 00014 ** You should have received a copy of the GNU General Public License 00015 ** along with this program; if not, write to the Free Software 00016 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. 00017 */ 00018 00019 #if !defined(__NumericLib_h) 00020 #define __NumericLib_h 00021 00022 /******************************************************************/ 00023 /* Numerical Routines */ 00024 /* */ 00025 /* author: B. Frenzel */ 00026 /* */ 00027 /* date: Jan 7 1993 */ 00028 /* */ 00029 /* description: */ 00030 /* test.c test environment for testing null() */ 00031 /* null.c main rootfinding routine */ 00032 /* muller.c rootfinding routine using muller method */ 00033 /* newton.c rootfinding routine using Newton method */ 00034 /* */ 00035 /* version: 1.2 (edited by LGM - INESC Porto - 06.2006) */ 00036 /* */ 00037 /* further necessary files: */ 00038 /* complex.c contains complex operations */ 00039 /* tools.c contains tools: */ 00040 /* fdvalue() computes P(x) and optional P'(x) */ 00041 /* poldef() deflates polynomial */ 00042 /* */ 00043 /* Copyright: */ 00044 /* Lehrstuhl fuer Nachrichtentechnik Erlangen */ 00045 /* Cauerstr. 7, 8520 Erlangen, FRG, 1993 */ 00046 /* e-mail: mrs_natural@nt.e-technik.uni-erlangen.de */ 00047 /* */ 00048 /******************************************************************/ 00049 00050 #include <cmath> 00051 #include <cstdio> 00052 #include <vector> 00053 #include <limits> 00054 #include <marsyas/common_header.h> 00055 #include <marsyas/realvec.h> 00056 00057 namespace Marsyas 00058 { 00069 class marsyas_EXPORT NumericLib 00070 { 00071 private: 00072 00073 /******************************************************************/ 00074 /* */ 00075 /* file: null.cpp */ 00076 /* */ 00077 /* main function: null() */ 00078 /* */ 00079 /* version: 1.3 (edited by LGM - INESC Porto - 29.04.02) */ 00080 /* */ 00081 /* author: B. Frenzel */ 00082 /* */ 00083 /* date: Jan 19 1993 */ 00084 /* */ 00085 /* input: p[] coefficient vector of the original */ 00086 /* polynomial */ 00087 /* pred[] coefficient vector of the deflated */ 00088 /* polynomial */ 00089 /* *n the highest exponent of the original */ 00090 /* polynomial */ 00091 /* flag flag = TRUE => complex coefficients */ 00092 /* flag = FALSE => real coefficients */ 00093 /* */ 00094 /* output: root[] vector of determined roots */ 00095 /* *maxerr estimation of max. error of all */ 00096 /* determined roots */ 00097 /* */ 00098 /* subroutines: poly_check(),quadratic(),lin_or_quad(),monic(), */ 00099 /* muller(),newton() */ 00100 /* */ 00101 /* description: */ 00102 /* main rootfinding routine for polynomials with complex or real */ 00103 /* coefficients using Muller's method combined with Newton's m. */ 00104 /* */ 00105 /* Copyright: */ 00106 /* Lehrstuhl fuer Nachrichtentechnik Erlangen */ 00107 /* Cauerstr. 7, 8520 Erlangen, FRG, 1993 */ 00108 /* e-mail: mrs_natural@nt.e-technik.uni-erlangen.de */ 00109 /* */ 00110 /******************************************************************/ 00111 unsigned char null(mrs_complex *p,mrs_complex *pred,mrs_natural *n,mrs_complex *root, 00112 mrs_real *maxerr,unsigned char flag); 00113 00114 unsigned char poly_check(mrs_complex *pred,mrs_natural *nred,mrs_natural *n,mrs_complex *root); 00115 void quadratic(mrs_complex *pred,mrs_complex *root); 00116 unsigned char lin_or_quad(mrs_complex *pred,mrs_natural nred,mrs_complex *root); 00117 void monic(mrs_complex *p,mrs_natural *n); 00118 00119 // MULLER 00120 /******************************************************************/ 00121 /* */ 00122 /* file: muller.c */ 00123 /* */ 00124 /* main function: muller() */ 00125 /* */ 00126 /* version: 1.2 */ 00127 /* */ 00128 /* author: B. Frenzel */ 00129 /* */ 00130 /* date: Jan 7 1993 */ 00131 /* */ 00132 /* input: pred[] coefficient vector of the deflated */ 00133 /* polynomial */ 00134 /* nred the highest exponent of the deflated */ 00135 /* polynomial */ 00136 /* */ 00137 /* return: xb determined root */ 00138 /* */ 00139 /* subroutines: initialize(),root_of_parabola(), */ 00140 /* iteration_equation(), compute_function(), */ 00141 /* check_x_value(), root_check() */ 00142 /* */ 00143 /* description: */ 00144 /* muller() determines the roots of a polynomial with complex */ 00145 /* coefficients with the Muller method; these roots are the */ 00146 /* initial estimations for the following Newton method */ 00147 /* */ 00148 /* Copyright: */ 00149 /* Lehrstuhl fuer Nachrichtentechnik Erlangen */ 00150 /* Cauerstr. 7, 91054 Erlangen, FRG, 1993 */ 00151 /* e-mail: mrs_natural@nt.e-technik.uni-erlangen.de */ 00152 /* */ 00153 /******************************************************************/ 00154 mrs_complex muller(mrs_complex *pred,mrs_natural nred); 00155 00156 mrs_complex x0MULLER_,x1MULLER_,x2MULLER_, /* common pomrs_naturals [x0,f(x0)=P(x0)], ... [x2,f(x2)] */ 00157 f0MULLER_,f1MULLER_,f2MULLER_, /* of parabola and polynomial */ 00158 h1MULLER_,h2MULLER_, /* distance between x2 and x1 */ 00159 q2MULLER_; /* smaller root of parabola */ 00160 mrs_natural iterMULLER_; /* iteration counter */ 00161 00162 void initialize(mrs_complex *pred,mrs_complex *xb,mrs_real *epsilon); 00163 void root_of_parabola(void); 00164 void iteration_equation(mrs_real *h2abs); 00165 void suppress_overflow(mrs_natural nred); 00166 void too_big_functionvalues(mrs_real *f2absq); 00167 void convergence_check(mrs_natural *overflow,mrs_real f1absq,mrs_real f2absq, 00168 mrs_real epsilon); 00169 void compute_function(mrs_complex *pred,mrs_natural nred,mrs_real f1absq, 00170 mrs_real *f2absq,mrs_real epsilon); 00171 void check_x_value(mrs_complex *xb,mrs_real *f2absqb,mrs_natural *rootd, 00172 mrs_real f1absq,mrs_real f2absq,mrs_real epsilon, 00173 mrs_natural *noise); 00174 void root_check(mrs_complex *pred,mrs_natural nred,mrs_real f2absqb,mrs_natural *seconditer, 00175 mrs_natural *rootd,mrs_natural *noise,mrs_complex xb); 00176 00177 //NEWTON 00178 /******************************************************************/ 00179 /* */ 00180 /* file: newton.c */ 00181 /* */ 00182 /* main function: newton() */ 00183 /* */ 00184 /* version: 1.2 */ 00185 /* */ 00186 /* author: B. Frenzel */ 00187 /* */ 00188 /* date: Jan 7 1993 */ 00189 /* */ 00190 /* input: p[] coefficient vector of the original */ 00191 /* polynomial */ 00192 /* n the highest exponent of the original */ 00193 /* polynomial */ 00194 /* ns determined root with Muller method */ 00195 /* */ 00196 /* output: *dxabs error of determined root */ 00197 /* */ 00198 /* return: xmin determined root */ 00199 /* subroutines: fdvalue() */ 00200 /* */ 00201 /* description: */ 00202 /* newton() determines the root of a polynomial with complex */ 00203 /* coefficients; the initial estimation is the root determined */ 00204 /* by Muller's method */ 00205 /* */ 00206 /* Copyright: */ 00207 /* Lehrstuhl fuer Nachrichtentechnik Erlangen */ 00208 /* Cauerstr. 7, 8520 Erlangen, FRG, 1993 */ 00209 /* e-mail: mrs_natural@nt.e-technik.uni-erlangen.de */ 00210 /* */ 00211 /******************************************************************/ 00212 mrs_complex newton(mrs_complex *p,mrs_natural n,mrs_complex ns,mrs_real *dxabs, 00213 unsigned char flag); 00214 00215 /******************************************************************/ 00216 /* */ 00217 /* file: tools.c */ 00218 /* */ 00219 /* version: 1.2 */ 00220 /* */ 00221 /* author: B. Frenzel */ 00222 /* */ 00223 /* date: Jan 7 1993 */ 00224 /* */ 00225 /* function: description: */ 00226 /* */ 00227 /* hornc() hornc() deflates the polynomial with coeff. */ 00228 /* stored in pred[0],...,pred[n] by Horner's */ 00229 /* method (one root) */ 00230 /* */ 00231 /* horncd() horncd() deflates the polynomial with coeff. */ 00232 /* stored in pred[0],...,pred[n] by Horner's */ 00233 /* method (two roots) */ 00234 /* */ 00235 /* poldef() decides whether to call hornc() or horncd() */ 00236 /* */ 00237 /* fdvalue() fdvalue() computes f=P(x0) by Horner's method; */ 00238 /* if flag=TRUE, it additionally computes the */ 00239 /* derivative df=P'(x0) */ 00240 /* */ 00241 /* Copyright: */ 00242 /* Lehrstuhl fuer Nachrichtentechnik Erlangen */ 00243 /* Cauerstr. 7, 8520 Erlangen, FRG, 1993 */ 00244 /* e-mail: mrs_natural@nt.e-technik.uni-erlangen.de */ 00245 /* */ 00246 /******************************************************************/ 00247 void hornc(mrs_complex *pred,mrs_natural nred,mrs_complex x0,unsigned char flag); 00248 void horncd(mrs_complex *pred,mrs_natural nred,mrs_real a,mrs_real b); 00249 mrs_natural poldef(mrs_complex *pred,mrs_natural nred,mrs_complex *root,unsigned char flag); 00250 void fdvalue(mrs_complex *p,mrs_natural n,mrs_complex *f,mrs_complex *df,mrs_complex x0, 00251 unsigned char flag) ; 00252 00253 static mrs_real add(mrs_real *a, mrs_real *b); 00254 static mrs_real pow_di(mrs_real *ap, mrs_natural *bp); 00255 00257 //hungarianAssignemt private methods 00259 static mrs_real mxGetInf(); 00260 static bool mxIsInf(mrs_real s); 00261 static void mxFree( void * s ); 00262 static void mexErrMsgTxt(const char * s); 00263 static void assignmentoptimal(mrs_natural *assignment, mrs_real *cost, mrs_real *distMatrix, mrs_natural nOfRows, mrs_natural nOfColumns); 00264 static void buildassignmentvector(mrs_natural *assignment, bool *starMatrix, mrs_natural nOfRows, mrs_natural nOfColumns); 00265 static void computeassignmentcost(mrs_natural *assignment, mrs_real *cost, mrs_real *distMatrix, mrs_natural nOfRows, mrs_natural nOfColumns); 00266 static void step2a(mrs_natural *assignment, mrs_real *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, mrs_natural nOfRows, mrs_natural nOfColumns, mrs_natural minDim); 00267 static void step2b(mrs_natural *assignment, mrs_real *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, mrs_natural nOfRows, mrs_natural nOfColumns, mrs_natural minDim); 00268 static void step3 (mrs_natural *assignment, mrs_real *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, mrs_natural nOfRows, mrs_natural nOfColumns, mrs_natural minDim); 00269 static void step4 (mrs_natural *assignment, mrs_real *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, mrs_natural nOfRows, mrs_natural nOfColumns, mrs_natural minDim, mrs_natural row, mrs_natural col); 00270 static void step5 (mrs_natural *assignment, mrs_real *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, mrs_natural nOfRows, mrs_natural nOfColumns, mrs_natural minDim); 00271 00272 public: 00273 00274 NumericLib(); 00275 ~NumericLib(); 00276 00277 bool polyRoots(std::vector<mrs_complex> coefs, bool complexCoefs, mrs_natural order, std::vector<mrs_complex> &roots); 00278 mrs_real determinant(realvec matrix); 00279 00280 static mrs_real gaussian(mrs_real x, mrs_real var, mrs_real mean) 00281 { 00282 return exp(-(x-mean)*(x-mean)/(2*var*var))/(var*sqrt(2*PI)); 00283 } 00284 00286 // metrics 00288 static mrs_real euclideanDistance(const realvec& Vi, const realvec& Vj, const realvec& covMatrix); 00289 static mrs_real mahalanobisDistance(const realvec& Vi, const realvec& Vj, const realvec& covMatrix); 00290 static mrs_real cosineDistance(const realvec& Vi, const realvec& Vj, const realvec& dummy = realvec()); 00291 static mrs_real cityblockDistance(const realvec& Vi, const realvec& Vj, const realvec& dummy = realvec()); 00292 static mrs_real correlationDistance(const realvec& Vi, const realvec& Vj, const realvec& dummy = realvec()); 00293 static mrs_real divergenceShape(const realvec& Ci, const realvec& Cj, const realvec& dummy = realvec()); 00294 static mrs_real bhattacharyyaShape(const realvec& Ci, const realvec& Cj, const realvec& dummy = realvec()); 00296 00297 // A is m x n 00298 // U is m x m 00299 // V is n x n 00300 // s is max(m,n)+1 x 1 00301 static void svd(mrs_natural m, mrs_natural n, realvec &A, realvec &U, realvec &V, realvec &s); 00302 00303 /* 00304 Householder reduction of matrix a to tridiagonal form. 00305 Algorithm: Martin et al., Num. Math. 11, 181-195, 1968. 00306 Ref: Smith et al., Matrix Eigensystem Routines -- EISPACK Guide 00307 Springer-Verlag, 1976, pp. 489-494. 00308 W H Press et al., Numerical Recipes in C, Cambridge U P, 00309 1988, pp. 373-374. 00310 00311 Source code adapted from F. Murtagh, Munich, 6 June 1989 00312 http://astro.u-strasbg.fr/~fmurtagh/mda-sw/pca.c 00313 */ 00314 //void tred2(realvec &**a, mrs_natural n, realvec &*d, realvec &*e); 00315 static void tred2(realvec &a, mrs_natural n, realvec &d, realvec &e); 00316 00317 /* 00318 Tridiagonal QL algorithm -- Implicit 00319 00320 Source code adapted from F. Murtagh, Munich, 6 June 1989 00321 http://astro.u-strasbg.fr/~fmurtagh/mda-sw/pca.c 00322 */ 00323 //void tqli(double d[], double e[], mrs_natural n, realvec &*z); 00324 static void tqli(realvec &d, realvec &e, mrs_natural n, realvec &z); 00325 00326 // Machine parameters 00327 // cmach : 00328 // 'B' | 'b' --> base 00329 // 'M' | 'm' --> digits in the mantissa 00330 // 'R' | 'r' --> approximation method : 1=rounding 0=chopping 00331 // 'E' | 'e' --> eps 00332 static mrs_real machp(const char *cmach); 00333 00335 //Hungarian assignment routine 00337 static mrs_real hungarianAssignment(realvec& matrixdist, realvec& matrixAssign); 00338 00339 }; 00340 00341 }//namespace Marsyas 00342 00343 #endif //__NumericLib_h 00344 00345