/************************************************************************
 *                                                                      *
 *  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 "nonlin.h"

extern universe *U;
extern bool verbose;


static str80 buf;

struct hypersphere_kernel_ *hypersphere_kernel_param = NULL;

#define PI              3.1415926535897931160
#define LN2		0.69314718055994530942	/* log(2) */
#define LOGPI_2		0.5723649429247000871	/* 1/2* log(PI)  */

#define ETA_DEFAULT (0.5)
#define MIN_ETA (0.0)
#define MAX_ETA (1.0)
static float eta;

/* "Numerical Recipes in C: The Art of Scientific Computing"
  published by Cambridge University Press (1988) */
static double gammln_NRC(xx)
double xx;
{
        double x,tmp,ser;
        static double cof[6]={76.18009173,-86.50532033,24.01409822,
                -1.231739516,0.120858003e-2,-0.536382e-5};
        int j;

        x=xx-1.0;
        tmp=x+5.5;
        tmp -= (x+0.5)*log(tmp);
        ser=1.0;
        for (j=0;j<=5;j++) {
                x += 1.0;
                ser += cof[j]/x;
        }
        return( -tmp+log(2.50662827465*ser) );
}


static void set_eta()
{
 float e;

 eta = ETA_DEFAULT;
 printf( "Parameter for the calculus of the radius: eta=%2.1f\b\b\b", eta );
 gets( buf );
 if( buf[0] != '\0' )
 {
   sscanf( buf, "%f", &e );
   if( e <= MIN_ETA || e >= MAX_ETA )
   {
     printf("Value %f for eta is invalid. Setting to default %f...",
          e, ETA_DEFAULT );
     e = ETA_DEFAULT;
     gets( buf );
   }
   eta = e;
 }
}


void init_nonlinear_dist()
{
 int i, dim, ni;
 double ln_aux, dd, ln_rho, aux, rho, Vol, arg;

 set_eta();
 printf("\nUsing Parzen estimator with hyperspheric kernel with eta=%f!\n",eta);
 printf("Initializing ... ");

 /* Inialize the hyperspheric kernel:
	Put for all classes and all possible dimensions the
	threshold parameter and inverse volume of the hypersphere
	into a lookup table
 */
 FREE( hypersphere_kernel_param );
 hypersphere_kernel_param = (struct hypersphere_kernel_ *)
        malloc( U->nrFeat * U->nrClass * sizeof( struct hypersphere_kernel_ ) );

 CHKPTR( hypersphere_kernel_param );

 for( dim = 1; dim < U->nrFeat; dim++ )
 {
   dd = (double)dim;
   ln_aux = LN2 + dd * LOGPI_2 - log( dd ) - gammln_NRC( dd / 2.0 );
   for( i = 0; i < U->nrClass; i++ )
   {
     /* calculate radius pp. 429, formula 16 */
     ni = U->C[i].numSampl;
     aux = - (double)eta * log( (double)ni );
     rho = exp( aux / dd );
     arg = ln_aux + aux;
     Vol = exp( arg );
     hypersphere_kernel_param[ dim * U->nrClass + i ].thresh = rho;
     hypersphere_kernel_param[ dim * U->nrClass + i ].inverseVol = 1.0 / Vol;
     /* printf("dim = %d  i = %d  rho=%lf  Vol=%lf  invVol=%lf\n",
	dim, i+1, rho, Vol, 1.0 / Vol );  /**/
   }
 }
 printf("done.\n");
}

void close_nonlinear_dist()
{
 FREE( hypersphere_kernel_param );
}


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

 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 ) return( invVol );
 return( 0.0 );
}