/************************************************************************ * * * Program package T O O L D I A G * * * * Version 1.5 * * Date: Tue Feb 8 13:39:08 1994 * * * * NOTE: This program package is copyrighted in the sense that it * * may be used for scientific purposes. The package as a whole, or * * parts thereof, cannot be included or used in any commercial * * application without written permission granted by the author. * * No programs contained in this package may be copied for commercial * * distribution. * * * * All comments concerning this program package may be sent to the * * e-mail address 'tr@fct.unl.pt'. * * * ************************************************************************/ #include #include #include "def.h" #include "featslct.h" #include "matrix.h" /* #define PEEPMAT(M) printf("M=\n");print_Matrix(&M);DBG; /**/ #define PEEPMAT(M) /**/ #define PI_TIMES_2 6.2831853071795862320 extern universe *U; extern MatElem trace(); extern bool return_log; extern float chernoff_parameter; static str80 buf; static Matrix *Means = NULL; static Matrix *Covariances = NULL; static Matrix mat0, mat1, mat2, mat3, mat4, mat5, mat6; MatElem Det1, Det2, Det12; static double select_prob_distMAHALANOBIS( a, b ) int a, b; { double Jab; Subtract_Matrix( &(Means[a]), &(Means[b]), &mat1 ); /* mat1 = m1-m2 */ PEEPMAT(mat1); Transpose_Matrix( &mat1, &mat2 ); /* mat2 = (m1 - m2) T */ PEEPMAT(mat2); Add_Matrix( &(Covariances[a]), &(Covariances[b]), &mat3 ); /* mat3=2cov */ ScalarMult_Matrix( 0.5, &mat3 ); /* mat3=cov */ PEEPMAT(mat3); Copy_Matrix( &mat3, &mat4 ); Invert_Matrix( &mat4, &Det12 ); /* mat4 = cov ^-1 */ PEEPMAT(mat4); Mult_Matrix( &mat2, &mat4, &mat5 ); Mult_Matrix( &mat5, &mat1, &mat2 ); PEEPMAT(mat2); Jab = (double)mat2.Elem[0]; /* printf("MAHALANOBIS: returning for %d and %d = %lf", a, b, Jab ); DBG; /**/ return( Jab ); } static double select_prob_distCHERNOFF( a, b ) int a, b; { double Jab, arg; MatElem s; s = (MatElem)chernoff_parameter; Subtract_Matrix( &(Means[a]), &(Means[b]), &mat1 ); /* mat1 = m1-m2 */ Transpose_Matrix( &mat1, &mat2 ); /* mat2 = (m1 - m2) T */ Copy_Matrix( &(Covariances[a]), &mat3 ); ScalarMult_Matrix( 1.0-s, &mat3 ); /* mat3 = (1-s)Cov1 */ Copy_Matrix( &(Covariances[b]), &mat4 ); ScalarMult_Matrix( s, &mat4 ); /* mat4 = s*Cov2 */ Add_Matrix( &mat3, &mat4, &mat5 ); Invert_Matrix( &mat5, &Det12 ); /* mat5 = ((1-s)Cov1+s*Cov2) ^-1 */ PEEPMAT(mat5); Copy_Matrix( &(Covariances[a]), &mat6 ); Invert_Matrix( &mat6, &Det1 ); Det1 = pow( Det1, 1.0-s ); Copy_Matrix( &(Covariances[b]), &mat6 ); Invert_Matrix( &mat6, &Det2 ); Det2 = pow( Det2, s ); Mult_Matrix( &mat2, &mat5, &mat6 ); Mult_Matrix( &mat6, &mat1, &mat2 ); mat2.Elem[0] *= 0.5*s*(1.0-s); arg = Det1*Det2; PEEPMAT(mat2); /* printf("Arg=%lf Det12=%lf", arg, Det12 ); DBG; /**/ if( arg < 0.0 || Det12 < 0.0 || (arg==0.0 && Det12<0.0) || (arg<0.0 && Det12==0.0) ) { fprintf(stderr,"Invalid argument %lf %lf- Exit...\n", arg, Det12 ); exit(1); } if( arg == 0.0 && Det12 == 0.0 ) Jab = (double)mat2.Elem[0]; else Jab = (double)( 0.5*log( Det12 / arg ) + mat2.Elem[0] ); return( Jab ); } static double select_prob_distBHATTACHARYYA( crit, a, b ) int crit; int a, b; { double Jab, arg; Subtract_Matrix( &(Means[a]), &(Means[b]), &mat1 ); /* mat1 = m1-m2 */ Transpose_Matrix( &mat1, &mat2 ); /* mat2 = (m1 - m2) T */ Add_Matrix( &(Covariances[a]), &(Covariances[b]), &mat3 ); /* mat3=2cov */ ScalarMult_Matrix( 0.5, &mat3 ); Copy_Matrix( &mat3, &mat4 ); Invert_Matrix( &mat4, &Det12 ); /* mat4 = 2cov ^-1 */ Copy_Matrix( &(Covariances[a]), &mat5 ); Invert_Matrix( &mat5, &Det1 ); Copy_Matrix( &(Covariances[b]), &mat5 ); Invert_Matrix( &mat5, &Det2 ); Mult_Matrix( &mat2, &mat4, &mat5 ); Mult_Matrix( &mat5, &mat1, &mat2 ); mat2.Elem[0] *= 0.125; arg = sqrt(Det1*Det2); PEEPMAT(mat2); /* printf("Arg=%lf Det1=%f Det2=%lf Det12=%lf",arg,Det1,Det2,Det12);NL; /**/ if( arg < 0.0 || Det12 < 0.0 || (arg==0.0 && Det12<0.0) || (arg<0.0 && Det12==0.0) ) { fprintf(stderr,"Invalid argument %lf %lf- Exit...\n", arg, Det12 ); exit(1); } if( arg == 0.0 && Det12 == 0.0 ) Jab = (double)mat2.Elem[0]; else Jab = (double)( 0.5*log( Det12 / arg ) + mat2.Elem[0] ); if( crit == BHATTACHARYYA_MATUSITA ) Jab = (double)sqrt( 2.0*(1.0 - exp(-(double)Jab)) ); return( Jab ); } static double select_prob_distPATRICK_FISHER( dim, a, b ) int dim, a, b; { double Jab, arg, pi2d, aux1, aux2, aux12; Subtract_Matrix( &(Means[a]), &(Means[b]), &mat1 ); /* mat1 = m1-m2 */ Transpose_Matrix( &mat1, &mat2 ); /* mat2 = (m1 - m2) T */ Copy_Matrix( &(Covariances[a]), &mat3 ); ScalarMult_Matrix( 2.0, &mat3 ); /* mat3 = 2*Cov1 */ Invert_Matrix( &mat3, &Det1 ); /* mat3 = 2*cov1 ^-1 */ Copy_Matrix( &(Covariances[b]), &mat4 ); ScalarMult_Matrix( 2.0, &mat4 ); /* mat4 = 2*Cov2 */ Invert_Matrix( &mat4, &Det2 ); /* mat4 = 2*cov2 ^-1 */ Add_Matrix( &(Covariances[a]), &(Covariances[b]), &mat5 ); /* mat5 = (cov1 + cov2) */ Invert_Matrix( &mat5, &Det12 ); Mult_Matrix( &mat2, &mat5, &mat6 ); Mult_Matrix( &mat6, &mat1, &mat2 ); arg = -0.5 * mat2.Elem[0]; pi2d = pow( PI_TIMES_2, (double)dim ); aux1 = 1.0 / sqrt( pi2d * Det1 ); aux2 = 1.0 / sqrt( pi2d * Det2 ); aux12 = 1.0 / sqrt( pi2d * Det12 ); Jab = aux1 + aux2 - 2.0*aux12*exp( arg ); /**/ /* printf("Det1==%lf Det2==%lf Det12==%lf aux1=%lf aux2=%lf aux12=%lf\n", Det1, Det2, Det12, aux1, aux2, aux12 ); /**/ return( Jab ); } static double select_prob_distDIVERGENCE( a, b ) int a, b; { double Jab; Subtract_Matrix( &(Means[a]), &(Means[b]), &mat1 ); /* mat1 = m1-m2 */ PEEPMAT(mat1); Transpose_Matrix( &mat1, &mat2 ); /* mat2 = (m1 - m2) T */ PEEPMAT(mat2); Copy_Matrix( &(Covariances[a]), &mat3 ); Invert_Matrix( &mat3, &Det1 ); /* mat3 = cov1 ^-1 */ PEEPMAT(mat3); Copy_Matrix( &(Covariances[b]), &mat4 ); Invert_Matrix( &mat4, &Det2 ); /* mat4 = cov2 ^-1 */ PEEPMAT(mat4); Add_Matrix( &mat3, &mat4, &mat5 ); /* mat5 = (cov1 ^-1 + cov2 ^-1) */ Mult_Matrix( &mat2, &mat5, &mat6 ); Mult_Matrix( &mat6, &mat1, &mat5 ); /* mat5 = first term */ Mult_Matrix( &mat3, &(Covariances[b]), &mat1 ); Mult_Matrix( &mat4, &(Covariances[a]), &mat2 ); Add_Matrix( &mat1, &mat2, &mat3 ); Subtract_Matrix( &mat3, &mat0, &mat1 ); PEEPMAT(mat1); Jab = 0.5 * (double)( mat5.Elem[0] + trace( &mat1 ) ); return( Jab ); } static double select_prob_distAverageClass( crit, dim ) int crit, dim; { double merit = 0.0; double Jab; int a, b; init_Matrix( &mat0 ); init_Matrix( &mat1 ); init_Matrix( &mat2 ); init_Matrix( &mat3 ); init_Matrix( &mat4 ); init_Matrix( &mat5 ); init_Matrix( &mat6 ); Unit_Matrix( &mat0, Covariances[0].row ); ScalarMult_Matrix( 2.0, &mat0 ); /* calculate the mutual criterion between two classes */ for( a = 0; a < U->nrClass-1; a++ ) { for( b = a+1; b < U->nrClass; b++ ) { switch( crit ) { case CHERNOFF : Jab = select_prob_distCHERNOFF( a, b ); break; case BHATTACHARYYA : case BHATTACHARYYA_MATUSITA : Jab = select_prob_distBHATTACHARYYA( crit, a, b ); break; case DIVERGENCE : Jab = select_prob_distDIVERGENCE( a, b ); break; case PATRICK_FISHER : Jab = select_prob_distPATRICK_FISHER( dim, a, b ); break; case MAHALANOBIS : Jab = select_prob_distMAHALANOBIS( a, b ); break; default: fprintf(stderr,"What the hell is crit %d? - exit...\n", crit ); exit(1); } merit += U->C[a].a_priori_prob * U->C[b].a_priori_prob * Jab; } } Matrix_Free( &mat0 ); Matrix_Free( &mat1 ); Matrix_Free( &mat2 ); Matrix_Free( &mat3 ); Matrix_Free( &mat4 ); Matrix_Free( &mat5 ); Matrix_Free( &mat6 ); /* printf("Returning merit = %lf", merit ); DBG; /**/ return( merit ); } float select_multivariate_prob_dist( crit, FSV, len ) int crit; FeatSelectVector FSV; int len; { float mean, scalar; double merit; int i, j, k, feat, a, b; Matrix *M = NULL; /* --- MEANS AND COVARAIANCES --- */ /* allocate space for means and covariances of all classes */ Means = (Matrix*) malloc( U->nrClass * sizeof( Matrix ) ); Covariances = (Matrix*) malloc( U->nrClass * sizeof( Matrix ) ); init_Matrix( &mat1 ); mat1.col = 1; mat1.row = len; Matrix_Alloc( &mat1 ); Reset_Matrix( &mat1 ); init_Matrix( &mat2 ); init_Matrix( &mat3 ); init_Matrix( &mat4 ); for( i = 0; i < U->nrClass; i++ ) { init_Matrix( &(Means[i]) ); init_Matrix( &(Covariances[i]) ); Means[i].col = 1; Means[i].row = len; Matrix_Alloc( &(Means[i]) ); Reset_Matrix( &(Means[i]) ); Covariances[i].col = len; Covariances[i].row = len; Matrix_Alloc( &(Covariances[i]) ); Reset_Matrix( &(Covariances[i]) ); /* generate mean for class i */ for( k = 0; k < len; k++ ) { feat = FSV[k].rank; mean = U->C[i].mean[feat]; Means[i].Elem[k] = (MatElem)mean; } /* printf("Mean vector of class %d = \n", i ); print_Matrix( &(Means[i]) ); DBG; /**/ /* generate covariance matrix for class i */ for( j = 0; j < U->C[i].numSampl; j++ ) { for( k = 0; k < len; k++ ) mat1.Elem[k] = (MatElem)U->C[i].S[j*U->nrFeat+FSV[k].rank]; Subtract_Matrix( &mat1, &(Means[i]), &mat2 ); Transpose_Matrix( &mat2, &mat3 ); Mult_Matrix( &mat2, &mat3, &mat4 ); Add_Matrix( &mat4, &(Covariances[i]), &mat2 ); Copy_Matrix( &mat2, &(Covariances[i]) ); } scalar = 1.0 / (MatElem)(U->C[i].numSampl - 1); ScalarMult_Matrix( scalar, &(Covariances[i]) ); /* Check the covariances, it may not have a zero element */ M = &(Covariances[i]); for( a = 0; a < M->row; a++ ) for( b = 0; b < M->col; b++ ) if( M->Elem[a*M->col+b] == 0.0 ) { fprintf(stderr,"\n\nHad some zero element in the covariance matrix "); fprintf(stderr,"of class %d - Exitus...\n", i+1 ); exit(1); } /* printf("Covariance matrix of class %d = \n", i ); print_Matrix( &(Covariances[i]) ); DBG; /**/ } Matrix_Free( &mat1 ); Matrix_Free( &mat2 ); Matrix_Free( &mat3 ); Matrix_Free( &mat4 ); /* ------------------------------ */ switch( crit ) { case CHERNOFF : case MAHALANOBIS : case BHATTACHARYYA : case BHATTACHARYYA_MATUSITA : case DIVERGENCE : case PATRICK_FISHER : merit = select_prob_distAverageClass( crit, len ); break; default: fprintf(stderr,"What the hell is crit %d? - exit...\n", crit ); exit(1); } /* Free means and covariances */ for( i = 0; i < U->nrClass; i++ ) { Matrix_Free( &(Means[i]) ); Matrix_Free( &(Covariances[i]) ); } FREE( Means ); FREE( Covariances ); if( return_log ) merit = log( merit ); return( (float)merit ); }