/************************************************************************
 *                                                                      *
 *  Program package         T O O L D I A G                             *
 *                                                                      *
 *  Version   1.5                                                       *
 *  Date:     Tue Feb  8 13:39:07 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 <stdio.h>
#include <math.h>
#include "def.h"
#include "matrix.h"

extern universe *U;
extern bool verbose;

static bool done;
static str80 buf;


#define EPSILON_WR (0.001)
#define EPSILON_WI (0.001)

static void feat_extract_discr_anal( SB, SW, W )
Matrix *SB, *SW, *W;
{
 Matrix SW1SB, SW_aux;
 MatElem Det;
 float maxLambda, *wr = NULL, *wi = NULL, *sortEigVal = NULL, *eigVec = NULL;
 int i, j, dim, maxI, nrEigVal, maxIndex;
 bool found;

 dim = SW->row;

 if( dim == 1 )
 {
   W->row = dim; W->col = dim; Matrix_Alloc( W );
   W->Elem[0] = 1.0;
   return;
 }

 init_Matrix( &SW1SB ); init_Matrix( &SW_aux );
 wr = (float*) malloc( dim * sizeof(float) ); CHKPTR(wr);
 wi = (float*) malloc( dim * sizeof(float) ); CHKPTR(wi);
 sortEigVal = (float*) malloc( dim * sizeof(float) ); CHKPTR(sortEigVal);
 eigVec = (float*) malloc( dim * sizeof(float) ); CHKPTR(eigVec);

 Copy_Matrix( SW, &SW_aux );
 printf("Inverting within classes scatter matrix..."); fflush(stdout);
 Invert_Matrix( &SW_aux, &Det );
 printf(" done.\nMultiplying matrices..."); fflush(stdout);
 Mult_Matrix( &SW_aux, SB, &SW1SB );
 /* printf("\n\nSW1SB:\n"); print_Matrix( &SW1SB ); /**/

 printf(" done.\nCalculating eigenvalues...\n\n");
 if( dim == 1 )
 {
   wr[0] = SW1SB.Elem[0]; wi[0] = 0.0;
 }
 else
 {
   balanc( &SW1SB ); /* printf("balanced:\n"); print_Matrix( &SW1SB ); /**/
   elmhes( &SW1SB ); /* printf("elmhes:\n"); print_Matrix( &SW1SB ); /**/
   hqr( &SW1SB, wr, wi );  /**/
 }
 for( i = 0; i < dim; i++ )
 {
   printf("Calculating eigenvalue nr. %d\r", i+1 ); fflush(stdout);
   /* printf("wr[%d]=%32.30f  wi[%d]=%32.30f\n", i, wr[i], i, wi[i] ); /**/
   /* Is there an imaginary part of a resonably large real eigenvalue? */
   if( fabs(wi[i]) > EPSILON_WI && fabs(wr[i]) > EPSILON_WR )
   {
     printf("Eigenvalue %d had some imaginary part: %.5f\n", i+1, wi[i] );
   }
 }

 /* Eigenvalues are determined - Calaculate corresponding eigenvectors */
 /* Sort in descending order */

 /* Analyse the eigenvalues first */
 maxIndex = MIN( dim, U->nrClass - 1 );
 /* printf("Maxindex=MIN(%d,%d) = %d", dim, U->nrClass - 1, maxIndex );DBG;/**/
 nrEigVal = 0; found = FALSE;
 while( ! found )
 {
   maxLambda = -INFINITY;
   for( i = 0; i < maxIndex; i++ )
   {
     if( wr[i] > maxLambda )
     {
       maxLambda = wr[i];
       maxI = i;
     }
   }
   if( maxLambda <= EPSILON_WR )
     found = TRUE;
   else
   {
     sortEigVal[ nrEigVal ] = maxLambda;
     nrEigVal++;
     wr[ maxI ] = -INFINITY;
     if( nrEigVal == maxIndex )
       found = TRUE;
   }
 }
 printf("\nN O N - Z E R O  E I G E N V A L U E S : ");
 if( nrEigVal == 0 )
 {
   printf("none\n Sorry, this data seems inappropriate for feature extraction");
   printf("..."); gets( buf );
 }
 else
 {
   printf("\n");
   for( i = 0; i < nrEigVal; i++ )
     printf("Eigenvalue nr. %3d of %d = %15.6f\n",i+1, nrEigVal, sortEigVal[i]);

   /* Ask user for new dimension */
   printf("\n");
   printf(" - DETERMINATION OF DIMENSION OF THE REDUCED FEATURE SPACE -\n\n");
   printf("As the minimum of the original dimension of the feature space(%d)",
	dim );
   printf("\nand the number of classes minus one(%d)\n", U->nrClass - 1 );
   printf("is %d, and since the number of non-zero eigenvalues was %d,\n",
	MIN( dim, U->nrClass - 1 ), nrEigVal );
   printf("choose a value for the reduced dimension between %d and %d: ",
	1, nrEigVal );
   get_d_range( &nrEigVal, 1, nrEigVal, LEFT_CLOSED__RIGHT_CLOSED );

   /* Create matrix with dimensions D x d */
   W->row = dim; W->col = nrEigVal;
   Matrix_Alloc( W );
   Reset_Matrix( W );

   for( i = 0; i < nrEigVal; i++ )
   {
     printf("Calculating eigenvector for eigenvalue nr. %d\r", i+1 );
     fflush(stdout);
     calcEigVec( &SW1SB, &(sortEigVal[i]), eigVec, i );
     for( j = 0; j < W->row; j++ )
       W->Elem[j*nrEigVal+i] = eigVec[j];
   }
   printf("\nNew values from the inverse iteration procedure:\n");
   for( i = 0; i < nrEigVal; i++ )
     printf("Eigenvalue nr. %3d of %d = %15.6f\n",i+1, nrEigVal, sortEigVal[i]);
 } 

 FREE(wr); FREE(wi); FREE(sortEigVal); FREE(eigVec);
 Matrix_Free( &SW1SB ); Matrix_Free( &SW_aux );
}


void featExtrDiscrimAnalysis( W )
Matrix *W;
{
 int i, j, dim;
 Matrix T, SB, SW;
 FeatSelectVector FSV = NULL;

 dim = U->nrFeat;

 FSV = (feature*) malloc( sizeof(feature) * dim ); CHKPTR(FSV);
 for( j = 0; j < dim; j++ )
 {
   FSV[ j ].rank = j;
   FSV[ j ].crit = (float)EMPTY;
 }

 init_Matrix( &T ); init_Matrix( &SB ); init_Matrix( &SW );
 printf("\n --- FEATURE EXTRACTION BASED ON DISCRIMINANT ANALYSIS ---\n\n");
 printf("Calculating %d x %d between and within classes scatter matrices...",
	dim, dim ); fflush(stdout);
 scatter_mats_kittler( FSV, dim, &T, &SB, &SW );
 printf(" done.\n");

 feat_extract_discr_anal( &SB, &SW, W );

 FREE(FSV);
 Matrix_Free( &T ); Matrix_Free( &SB ); Matrix_Free( &SW );
}