Marsyas: /home/marsyas/marsyas/src/marsyas/NumericLib.h Source 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 "common_header.h"
00055 #include "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