/*
* All questions regarding the software should be addressed to
* 
*       Prof. Charles A. Bouman
*       Purdue University
*       School of Electrical and Computer Engineering
*       1285 Electrical Engineering Building
*       West Lafayette, IN 47907-1285
*       USA
*       +1 765 494 0340
*       +1 765 494 3358 (fax)
*       email:  bouman@ecn.purdue.edu
*       http://www.ece.purdue.edu/~bouman
* 
* Copyright (c) 1995 The Board of Trustees of Purdue University.
*
* Permission to use, copy, modify, and distribute this software and its
* documentation for any purpose, without fee, and without written agreement is
* hereby granted, provided that the above copyright notice and the following
* two paragraphs appear in all copies of this software.
*
* IN NO EVENT SHALL PURDUE UNIVERSITY BE LIABLE TO ANY PARTY FOR DIRECT,
* INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
* USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF PURDUE UNIVERSITY HAS
* BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
* PURDUE UNIVERSITY SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
* PARTICULAR PURPOSE.  THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" BASIS,
* AND PURDUE UNIVERSITY HAS NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT,
* UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
*/


#define PI 3.141592654

#include<stdio.h>
#include<stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>

#include "clust_defs.h"

#include "clust_io.h"
#include "clust_util.h"
#include "alloc_util.h"
#include "clustermessage.h"

#include "Main_def.h"

/* Subroutine in eigen.c */
void eigen(
  double **M,     /* Input matrix */
  double *lambda, /* Output eigenvalues */
  int n           /* Input matrix dimension */
);


static void LogLikelihood_init(struct SigSet *S);
int G_ludcmp(double **a, int n, int *indx, double *d);

int invert(
  double **a, /* input/output matrix */
  int      n  /* dimension */
);

void LogLikelihood(
  double *vector, 
  double *ll,          /* log likelihood, ll[class] */
  struct SigSet *S,    /* class signatures */
  double text_cost,
  double non_text_cost
);


int classify(
double  **vector,        /* i : feature vector */
int      num,            /* i : # of the vector */
double  **loglikelihood, /* o : loglikelihood */
int     *clus,           /* o : ML estimation */
int      dim,            /* i : vector dimension */
char    *filename,       /* i : parameter file name */
double   text_cost,      /* i : cost toward text */
double   non_text_cost
) 
{
    FILE *fp;
    int i,j;
    struct SigSet S;
    char *paramfname;
    int NDataVectors,NRead;
    double **data,tmp,*ll;
    int maxindex;
    double maxval;

    /* Read File name : Please modify the path */
    paramfname = filename;

    /* Read SigSet from parameter file */
    if((fp=fopen(paramfname,"r"))==NULL) {
      fprintf(stderr, "\nError: Can't open parameter file %s\n",paramfname); 
      exit(-1);
    }
    I_ReadSigSet (fp,&S);
    fclose(fp);

    /* Parameter check */
    if (S.nbands != dim ) {
      fprintf(stderr, "\nError: Dimension of Input feature vector is wrong.\n"); 
      exit(-1);
    }

    /* Single vector input*/
    NDataVectors = num;


    /* Read data */
    data = G_alloc_matrix(NDataVectors,S.nbands);
    for(i=0; i<NDataVectors; i++) {
      for(j=0; j<S.nbands; j++) 
        data[i][j]=vector[i][j];
    }

    /* Initialize constants for Log likelihood calculations */
    LogLikelihood_init(&S);

    /* Compute Log likelihood for each class*/
    ll = G_alloc_vector(S.nclasses);

    for(i=0; i<NDataVectors; i++) {
      LogLikelihood(data[i],ll,&S,text_cost, non_text_cost); 

      maxval = ll[0];
      maxindex = 0;
      for (j=0; j<S.nclasses; j++) {
        loglikelihood[i][j] = ll[j];
        if( ll[j] > maxval ) {
          maxval = ll[j];
          maxindex = j;
        }
      }
      /* Final output */
      clus[i] = maxindex;
    }

    G_free_vector(ll);
    G_free_matrix(data);
    return(0);
}




void LogLikelihood(
  double *vector, 
  double *ll,          /* log likelihood, ll[class] */
  struct SigSet *S,    /* class signatures */
  double  text_cost,    /* cost */
  double  non_text_cost
)
{
    int  m;               /* class index */
    int  k;               /* subclass index */
    int  b1,b2;           /* spectral index */
    int  max_nsubclasses; /* maximum number of subclasses */
    int  nbands;          /* number of spectral bands */
    double *subll;        /* log likelihood of subclasses */
    double *diff; 
    double maxlike;
    double subsum;
    struct ClassSig *C;
    struct SubSig *SubS;

    nbands = S->nbands;

    /* determine the maximum number of subclasses */
    max_nsubclasses = 0;
    for(m=0; m<S->nclasses; m++ )
      if(S->ClassSig[m].nsubclasses>max_nsubclasses)
        max_nsubclasses = S->ClassSig[m].nsubclasses;

    /* allocate memory */
    diff = (double *)G_malloc(nbands*sizeof(double));
    subll = (double *)G_malloc(max_nsubclasses*sizeof(double));

    /* Compute log likelihood for each class */

    /* for each class */
    for(m=0; m<S->nclasses; m++ ) {
      C = &(S->ClassSig[m]);

      /* compute log likelihood for each subclass */
      for(k=0; k<C->nsubclasses; k++) {
        SubS = &(C->SubSig[k]);
        subll[k] = SubS->cnst;
        for(b1=0; b1<nbands; b1++) {
          diff[b1] = vector[b1] - SubS->means[b1];
          subll[k] -= 0.5*diff[b1]*diff[b1]*SubS->Rinv[b1][b1];
        }
        for(b1=0; b1<nbands; b1++) 
        for(b2=b1+1; b2<nbands; b2++)
          subll[k] -= diff[b1]*diff[b2]*SubS->Rinv[b1][b2];
      }

      /* shortcut for one subclass */
      if(C->nsubclasses==1) {
        ll[m] = subll[0];
      }
      /* compute mixture likelihood */
      else {
        /* find the most likely subclass */
        for(k=0; k<C->nsubclasses; k++)
        {
          if(k==0) maxlike = subll[k];
          if(subll[k]>maxlike) maxlike = subll[k];
        }

        /* Sum weighted subclass likelihoods */
        subsum = 0;
        for(k=0; k<C->nsubclasses; k++)
          subsum += exp( subll[k]-maxlike )*C->SubSig[k].pi;

        ll[m] = log(subsum) + maxlike;
        /* Cost added by Eri */
        /* log p(y|0) - log p(y|1) < or >  log ( pi_1/pi_0 ) */
        /* log p(y|0) + log pi_0 - log p(y|1) - log pi_1 < or > 0 */
        ll[1] += text_cost;  /* More text */
        ll[0] += non_text_cost;  /* More non-text */
      }
    }
    free((char *)diff);
    free((char *)subll);
}




static void LogLikelihood_init(struct SigSet *S)
{
   int m; 
   int i;
   int b1,b2;
   int nbands;
   double *lambda;
   struct ClassSig *C;
   struct SubSig *SubS;

   nbands = S->nbands;
   /* allocate scratch memory */
   lambda = (double *)G_malloc(nbands*sizeof(double));

   /* invert matrix and compute constant for each subclass */

   /* for each class */
   for(m=0; m<S->nclasses; m++ ) {
     C = &(S->ClassSig[m]);

     /* for each subclass */
     for(i=0; i<C->nsubclasses; i++) {
        SubS = &(C->SubSig[i]);

        /* Test for symetric  matrix */
        for(b1=0; b1<nbands; b1++)
        for(b2=0; b2<nbands; b2++) {
          if(SubS->R[b1][b2]!=SubS->R[b2][b1]) {
            fprintf(stderr,"\nWarning: nonsymetric covariance for class %d ",m+1);
            fprintf(stderr,"Subclass %d\n",i+1);
          }
          SubS->Rinv[b1][b2] = SubS->R[b1][b2];
        }

        /* Test for positive definite matrix */
        eigen(SubS->Rinv,lambda,nbands);
        for(b1=0; b1<nbands; b1++) {
          if(lambda[b1]<=0.0) {
            fprintf(stderr,"Warning: nonpositive eigenvalues for class %d",m+1);
            fprintf(stderr,"Subclass %d\n",i+1);
          }
        }

        /* Precomputes the cnst */
        SubS->cnst = (-nbands/2.0)*log(2*PI);
        for(b1=0; b1<nbands; b1++) {
            SubS->cnst += - 0.5*log(lambda[b1]);
        }

        /* Precomputes the inverse of tex->R */
        invert(SubS->Rinv,nbands);
      }
    }
    free((char *)lambda);
}

