Marsyas  0.6.0-alpha
/usr/src/RPM/BUILD/marsyas-0.6.0/src/marsyas/NumericLib.h
Go to the documentation of this file.
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