/************************************************************************
 *                                                                      *
 *  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 <stdio.h>
#include <math.h>
#include "def.h"
#include "featslct.h"
#include "nonlin.h"

#define SOON {printf("Soon at the theater near to you...type <Ctrl> c");DBG;}

extern universe *U;
extern int dist;
extern float minkowski_order;
extern int euclid_dist;
extern bool mutual_all;

extern float selectINTER_CLASS_DISTANCE_Euclid();

static str80 buf;
static FeatVector zeta_ik = NULL, zeta_jl = NULL;


static float dist_CITY_BLOCK( vec1, vec2, len )
FeatVector  vec1, vec2;
int len;
{
 int k;
 float d = 0.0;

 for( k = 0; k < len; k++ )
   d += fabs( vec1[k] - vec2[k] );
 return( d );
}


float dist_EUCLIDEAN( vec1, vec2, len )
FeatVector  vec1, vec2;
int len;
{
 int k;
 float dSqr = 0.0, d;

 for( k = 0; k < len; k++ )
   dSqr += (vec1[k] - vec2[k])*(vec1[k] - vec2[k]);
 d = (float)sqrt( (double)dSqr );
 /* printf("dsqr=%f Euclidean distance = %f between vectors:\n", dSqr, d );
 showFV( len, vec1 ); showFV( len, vec2 );	/**/
 return( d );
}


extern struct hypersphere_kernel_ *hypersphere_kernel_param;

static float dist_NONLINEAR( vec1, vec2, len, class_i )
FeatVector  vec1, vec2;
int len, class_i;
{
 float dE, rho, invVol, Kernel, dist;

 dE =  dist_EUCLIDEAN( vec1, vec2, len );
 rho = (float)hypersphere_kernel_param[len*U->nrClass + class_i].thresh; 
 invVol = (float)hypersphere_kernel_param[len*U->nrClass + class_i].inverseVol;

 if( dE <= rho )
   Kernel = invVol;
 else
   Kernel = 0.0;
 return( Kernel );
}


static float dist_CHEBYCHEV( vec1, vec2, len )
FeatVector  vec1, vec2;
int len;
{
 int k;
 float d = 0.0, maxD = -INFINITY;

 for( k = 0; k < len; k++ )
 {
   d = (float)fabs( vec1[k] - vec2[k] );
   if( d > maxD )
     maxD = d;
 }
 if( maxD == -INFINITY )
     { fprintf(stderr,"Problems in dist_CHEBYCHEV -exit...\n"); exit(1); }
 return( maxD );
}



static float dist_S( vec1, vec2, len, s )
FeatVector  vec1, vec2;
int len;
float s;
{
 int k;
 float d = 0.0, arg;

 for( k = 0; k < len; k++ )
 {
   arg = vec1[k] - vec2[k];
   d += (float)pow( fabs((double)arg), (double)s );
 }
 d = (float)pow( (double)d, 1.0/(double)s );  
 return( d );
}


static float selectINTER_CLASS_DISTANCE( FSV, len, i, j )
FeatSelectVector FSV;
int len, i, j;
{
 int k, l, a;
 float d, dij = 0.0;

 /* calculate all mutual distances between all samples of class i and j */
 for( k = 0; k < U->C[i].numSampl; k++ )
 {
   for( l = 0; l < U->C[j].numSampl; l++ )
   {
     /* prepare the two vectors */
     for( a = 0; a < len; a++ )
     {
       zeta_ik[ a ] = U->C[i].S[ k * U->nrFeat + FSV[a].rank ];
       zeta_jl[ a ] = U->C[j].S[ l * U->nrFeat + FSV[a].rank ];
     }
     /* printf("\nd( zeta(%d,%d) , zeta(%d,%d) )\n", i,k,j,l );
     showFV( len, zeta_ik ); showFV( len, zeta_jl );	/**/
     /* calculate the distance between the two vectors */
     /* dist is global and contains the distance metrik */
     switch( dist )
     {
       case MINKOWSKI : d = dist_S( zeta_ik, zeta_jl, len, minkowski_order );
	 break;
       case CITY_BLOCK : d = dist_CITY_BLOCK( zeta_ik, zeta_jl, len );break;
       case EUCLIDEAN : d = dist_EUCLIDEAN( zeta_ik, zeta_jl, len ); break;
       case CHEBYCHEV_DIST : d = dist_CHEBYCHEV( zeta_ik, zeta_jl, len ); break;
       case NONLINEAR : d = dist_NONLINEAR( zeta_ik, zeta_jl, len, i ); break;
       default: fprintf(stderr,"Unknown metrik - exit\n");exit(1);
     }
     dij += d;
   }
 }
 /* printf("d(%d,%d)=%f", i, j, dij ); DBG;	/**/
 return( dij );
}


static float selectInterClassdistAverageClass( crit, FSV, len )
int crit;
FeatSelectVector FSV;
int len;
{
 float merit = 0.0, Jij, auxJij, auxJ, auxI, H, dummy=0.0;
 int i, j, boundI, boundJ;
 
 zeta_ik = (FeatVector) malloc(len*sizeof(FeatVector*));
 zeta_jl = (FeatVector) malloc(len*sizeof(FeatVector*));

 /* calculate the mutual criterion between two classes */
 /* between *ALL* classes, inclusively between the same class */
 /* Therefor loop i = 0, j = 0 */
 auxI = 0.0;
 if( mutual_all )
   boundI = U->nrClass;
 else
   boundI = U->nrClass - 1;
 for( i = 0; i < boundI; i++ )
 {
   auxJ = 0.0;
   if( mutual_all )
     boundJ = 0;
   else
     boundJ = i+1;
   for( j = boundJ; j < U->nrClass; j++ )
   {
     switch( crit )
     {
       case INTER_CLASS_DISTANCE :
         switch( dist )
         {
           case MINKOWSKI :
           case CITY_BLOCK :
           case CHEBYCHEV_DIST :
           case EUCLIDEAN :
           case NONLINEAR :
	     Jij = selectINTER_CLASS_DISTANCE( FSV, len, i, j );
             break;
           default: fprintf(stderr,"Unknown metrik - exit\n");exit(1);
         }
         /* dummy += Jij;	/**/
         auxJij = Jij / ( (float)U->C[i].numSampl * (float)U->C[j].numSampl );
         if( dist == NONLINEAR )
         {
           H = (float)hypersphere_kernel_param[len*U->nrClass + i].inverseVol;
           auxJij = H - auxJij;
         }
         auxJ += U->C[j].a_priori_prob * auxJij;
         break;
       default: fprintf(stderr,"What the hell is crit %d? - exit...\n", crit );
                exit(1);
     }
     /* printf("auxI = %f auxJ = %f\n", auxI, auxJ ); /**/
   }
   auxI += U->C[i].a_priori_prob * auxJ; 
 }
 merit = auxI / 2.0; /**/
 /* dummy /= (2.0*(float)U->sumSampl*(float)U->sumSampl);
 printf("INSIDE: merit=%f  dummy=%f", merit, dummy ); DBG; /**/
 
 FREE( zeta_ik ); FREE( zeta_jl );
 return( merit );
}


float select_multivariate_InterClassDist( crit, FSV, len )
int crit;
FeatSelectVector FSV;
int len;
{
 float merit;

 switch( crit )
 {
   case INTER_CLASS_DISTANCE :
     switch( dist )
     {
       case MINKOWSKI :
       case CITY_BLOCK :
       case CHEBYCHEV_DIST :
       case NONLINEAR :
         merit = selectInterClassdistAverageClass( crit, FSV, len ); break;
       case EUCLIDEAN :
         if( euclid_dist == EMPTY )
           merit = selectInterClassdistAverageClass( crit, FSV, len );
         else
	   merit = selectINTER_CLASS_DISTANCE_Euclid( FSV, len );
         break;
       default: fprintf(stderr,"Unknown metrik - exit\n");exit(1);
     }
     break;
   default: fprintf(stderr,"What the hell is crit %d? - exit...\n", crit );
        exit(1);
 }
 return( merit );
}