/*
  File name: parametric.c
  Created by: Ljubomir Buturovic
  Created: 08/03/2004
  Purpose: implementation of parametric classifiers (linear and quadratic).
*/

/*
  Copyright 2004 Ljubomir J. Buturovic

  Permission is hereby granted, free of charge, to any person
  obtaining a copy of this software and associated documentation files
  (the "Software"), to deal in the Software without restriction,
  including without limitation the rights to use, copy, modify, merge,
  publish, distribute, sublicense, and/or sell copies of the Software,
  and to permit persons to whom the Software is furnished to do so,
  subject to the following conditions:

  The above copyright notice and this permission notice shall be
  included in all copies or substantial portions of the Software.

  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
  BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
  ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
  CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  SOFTWARE.
*/

static char rcsid[] = "$Id: parametric.c,v 1.11 2005/02/23 05:48:00 ljubomir Exp $";

#include <stdlib.h>
#include <errno.h>
#include <unistd.h>
#include <string.h>
#include <math.h>
#include <float.h>
#include "parametric.h"
#include "dataset.h"
#include "lin.h"
#include "lau.h"
#include "lmat.h"
#include "lerr.h"

/*
  Return weighted class-conditional estimate of the covariance matrix
  for 'dset'. The estimate is the sum of class-conditional covariance
  matrices of individual classes, weighted by the class probabilities.

  Return NULL in case of error and set errno.
*/
static float **weighted_est(struct dataset *dset)
{
  int   i;
  int   j;
  int   k;
  int   status;
  int   d;
  int   offset;
  float factor;
  float **sigma;
  float **csigma; /* class-conditional covariance matrix */

  sigma = (float **) 0;
  if (dset)
    {
      d = dset->d;
      sigma = fmx_alloc(d, d);
      if (sigma)
	{
	  fmx_set(sigma, d, d, 0.0);
	  status = 0;
	  offset = 0;
	  for (i = 0; (i < dset->c) && !status; i++)
	    {
	      csigma = cest(&dset->x[offset], d, dset->nd[i], COVARIANCE_MATRIX);
	      if (csigma)
		{
		  factor = dset->nd[i];
		  factor = factor/dset->nv;
		  for (j = 0; j < dset->d; j++)
		    for (k = 0; k < dset->d; k++)
		      sigma[j][k] += factor*csigma[j][k];
		  mx_free((void **) csigma, d);
		  offset += dset->nd[i];
		}
	      else
		{
		  status = -1;
		  sigma = (float **) mx_free((void **) sigma, d);
		}
	    }
	}
    }
  return sigma;
}

/*
  Return parametric linear classifier trained on 'dset'. The
  classifier is a matrix c by d+1, where 'c' is number of classes and
  'd' is dimension of input vectors.

  In case of error, return NULL and set errc to the error code. Error
  codes are EINVAL for bad input arguments, malloc() error codes and
  LERR_SINGCOV for singular covariance matrix.

  Uses equations (60-62) from R. O. Duda, P. E. Hart and D. G. Stork,
  Pattern Classification, Second Edition, John Wiley & Sons, Inc.,
  2001.
*/
float **lin_learn(int mode, struct dataset *dset, int *errc)
{
  int    i;
  int    j;
  int    k;
  int    d;
  int    c;
  int    offset;
  float  fnv;
  float  *mean;
  float  **wmx;
  float  **sigma; /* estimate of cumulative covariance matrix */
  float  **sig;   /* inverse of sigma */
  double s;

  wmx = (float **) 0;
  if (dset && errc)
    {
      *errc = 0;
      d = dset->d;
      c = dset->c;
      fnv = dset->nv;
      if (mode == CUMULATIVE_COV)
	sigma = cest(dset->x, d, dset->nv, COVARIANCE_MATRIX);
      else
	sigma = weighted_est(dset);
      if (sigma)
	{
	  sig = fmx_inv(sigma, d, (float *) 0, (float *) 0, errc);
	  mx_free((void **) sigma, d);
	  if (sig)
	    {
	      wmx = fmx_alloc(c, d+1);
	      if (wmx)
		{
		  offset = 0;
		  for (i = 0; (i < c) && !*errc; i++)
		    {
		      mean = mest(&dset->x[offset], d, dset->nd[i]);
		      if (mean)
			{
			  offset += dset->nd[i];
			  for (j = 0; j < d; j++)
			    {
			      s = 0;
			      for (k = 0; k < d; k++)
				s += sig[j][k]*mean[k];
			      wmx[i][j] = s;
			    }
			  s = 0;
			  for (j = 0; j < d; j++)
			    s += mean[j]*wmx[i][j];
			  wmx[i][d] = -0.5*s+log(dset->nd[i]/fnv);
			  free(mean);
			}
		      else
			{
			  *errc = errno;
			  wmx = (float **) mx_free((void **) wmx, c);
			}
		    }
		}
	      else 
		*errc = errno;
	      mx_free((void **) sig, d);
	    }
	}
      else 
	*errc = errno;
    }
  else if (errc)
    *errc = EINVAL;
  return wmx;
}

void qmodel_free(struct qmodel *model, int c)
{
  if (model)
    {
      mx_free((void **) model->wmx, c);
      mx_free((void **) model->sigma, c);
      vx_free(model);
    }
}

/*
  Return partial parametric quadratic classifier trained on
  'dset'. The partial classifier is a matrix c by d+1, where 'c' is
  number of classes and 'd' is dimension of input vectors.

  The function uses equations (66-69) from R. O. Duda, P. E. Hart and
  D. G. Stork, Pattern Classification, Second Edition, John Wiley &
  Sons, Inc., 2001, Section 2.6.3. According to these equations, the
  quadratic classifier consists of a quadratic and linear parts. The
  quadratic part is trivial; it is simply -0.5 times the inverted
  covariance matrix for the corresponding class. Therefore the
  function does not return it, but instead merely computes the
  class-conditional inverted covariance matrices and stores them in
  dset->sigma. The returned array is the linear part of the
  classifier.

  In case of error, return NULL and set errc to the error code. Error
  codes are EINVAL for bad input arguments, malloc() error codes and
  LERR_SINGCOV for singular covariance matrix.
*/
struct qmodel *pqc_learn(struct dataset *dset, int *errc)
{
  int    i;
  int    d;
  int    c;
  int    offset;
  int    status;
  int    length;
  float  fnv;
  float  det;
  double s;
  float  *mean;
  float  **wmx;
  struct qmodel *pqc;

  wmx = (float **) 0;
  c = 0;
  length = 0;
  pqc = (struct qmodel *) 0;
  if (dset && errc)
    {
      *errc = 0;
      d = dset->d;
      c = dset->c;
      length = d*(d+1)/2;
      pqc = malloc(sizeof(struct qmodel));
      if (pqc)
	{
	  /*
	    Compute class-conditional inverse covariance matrices.
	  */
	  status = dataset_sigma(dset, errc);
	  if (!status)
	    {
	      wmx = fmx_alloc(c, d+1);
	      if (wmx)
		{
		  offset = 0;
		  for (i = 0; (i < c) && !*errc; i++)
		    {
		      fnv = dset->nd[i];
		      mean = mest(&dset->x[offset], d, dset->nd[i]);
		      if (mean)
			{
			  fvec_smx(mean, dset->sigma[i], d, wmx[i]);
			  /*
			    Bias (constant) term.
			  */
			  s = fvec_dot(wmx[i], mean, d, errc);
			  det = log(dset->det[2*i])+dset->det[2*i+1];
			  wmx[i][d] = -0.5*(s+det)+log(dset->nd[i]/fnv);
			  free(mean);
			  offset += dset->nd[i];
			}
		      else
			{
			  *errc = errno;
			  wmx = (float **) mx_free((void **) wmx, c);
			  free(pqc);
			  pqc = (struct qmodel *) 0;
			}
		    }
		}
	      else 
		{
		  *errc = errno;
		  free(pqc);
		  pqc = (struct qmodel *) 0;
		}
	    }
	  else
	    {
	      free(pqc);
	      pqc = (struct qmodel *) 0;
	    }
	}
      else
	*errc = errno;
    }
  else if (errc)
    *errc = EINVAL;
  if (errc && !*errc)
    {
      pqc->wmx = wmx;
      pqc->sigma = fmx_clone(dset->sigma, c, length);
      if (!pqc->sigma)
	{
	  if (errc)
	    *errc = errno;
	  mx_free((void **) wmx, c);
	  free(pqc);
	  pqc = (struct qmodel *) 0;
	}
    }
  return pqc;
}

/*
  Save parametric quadratic classifier defined by 'qcx' in
  'fname'. 'd' is dimension of the problem, 'c' is number of classes.

  The classifier is saved in the following format: for each class, the
  linear portion of the classifier is saved first (which is a vector
  of length d+1), followed by the inverse of the class-conditional
  covariance matrix. To perform classification, use equations (66-69)
  from Duda, Hart & Stork.

  Return 0 for success, error code for failure. 
*/
int pqc_save(struct qmodel *qcx, int d, int c, char *fname)
{
  int  i;
  int  j;
  int  length;
  int  errc;
  int  status;
  FILE *fptr;

  errc = 0;
  status = 0;
  if (qcx && fname)
    {
      length = d*(d+1)/2;
      fptr = fopen(fname, "w");
      if (fptr)
	{
	  for (i = 0; (i < c) && (status >= 0); i++)
	    {
	      for (j = 0; (j <= d) && (status >= 0); j++)
		status = fprintf(fptr, "%g\t", qcx->wmx[i][j]);
	      if (status >= 0)
		{
		  /*
		    Save the inverse covariance matrix in Symmetric
		    Storage Mode.
		  */
		  for (j = 0; (j < length) && (status >= 0); j++)
		    {
		      if (j == length-1)
			status = fprintf(fptr, "%g\n", qcx->sigma[i][j]);
		      else
			status = fprintf(fptr, "%g\t", qcx->sigma[i][j]);
		    }
		}
	    }
	  status = fclose(fptr);
	  if (!status)
	    errc = errno;
	  /*
	    Try to close the file even if there was an error.
	  */
	  if (errc)
	    fclose(fptr);
	}
      else
	errc = errno;
    }
  else
    errc = EINVAL;
  return errc;
}

/*
  Append parametric quadratic classifier defined by 'qcx' to 'fptr',
  in 'bagging' format. The first column is total number of models
  'nmodels, followd by model index 'mdx' and model weight 'weight'.

  Return 0 for success, error code for failure. 
*/
int pqc_write(FILE *fptr, int nmodels, struct qmodel *qcx, int d, int c, int mdx, 
	      float weight)
{
  int  i;
  int  j;
  int  length;
  int  errc;
  int  status;

  errc = 0;
  status = 0;
  if (qcx && fptr)
    {
      length = d*(d+1)/2;
      for (i = 0; (i < c) && (status >= 0); i++)
	{
	  fprintf(fptr, "%d\t%d\t%g\t", nmodels, mdx, weight);
	  for (j = 0; (j <= d) && (status >= 0); j++)
	    status = fprintf(fptr, "%g\t", qcx->wmx[i][j]);
	  if (status >= 0)
	    {
	      /*
		Save the quadratic term coefficients, which are in
		fact inverse class-conditional covariance matrix in
		Symmetric Storage Mode.
	      */
	      for (j = 0; (j < length) && (status >= 0); j++)
		{
		  if (j == length-1)
		    status = fprintf(fptr, "%g\n", qcx->sigma[i][j]);
		  else
		    status = fprintf(fptr, "%g\t", qcx->sigma[i][j]);
		}
	    }
	}
    }
  else
    errc = EINVAL;
  return errc;
}

/*
  Memory allocator for pqc_load_models(). Return -1 and set errno in
  case of failure.
*/
static int pqc_load_models_alloc(struct qmodel ***models, float **wts, char **line, 
				 int model_cnt, int c, int ad, int length,
				 int llen)
{
  int   i;
  int   status;
  struct qmodel **qcx;

  status = 0;
  qcx = malloc(model_cnt*sizeof(struct qmodel *));
  if (qcx)
    {
      for (i = 0; (i < model_cnt) && !status; i++)
	{
	  qcx[i] = malloc(sizeof(struct qmodel));
	  if (qcx[i])
	    {
	      qcx[i]->wmx = fmx_alloc(c, ad);
	      if (qcx[i]->wmx)
		{
		  qcx[i]->sigma = fmx_alloc(c, length);
		  if (!qcx[i]->sigma)
		    status = -1;
		}
	      else
		status = -1;
	    }
	  else
	    status = -1;
	}
      if (!status)
	{
	  *wts = malloc(model_cnt*sizeof(float));
	  if (*wts)
	    {
	      *line = malloc(llen+2);
	      if (*line)
		{
		  if (model_cnt == 1)
		    *wts[0] = 1.0;
		}
	      else
		status = -1;
	    }
	  else
	    status = -1;
	}
      else
	status = -1;
    }
  else
    status = -1;
  if (!status)
    *models = qcx;
  return status;
}

/*
  Load a committee of quadratic classifiers from 'fname'. The function
  returns an array of pointers to 'qmodel' structures, one structure
  per model. Each struct defines a parametric quadratic classifier.

  The function automatically recognizes two types of model files:
  single model (MODEL_SINGLE) or multiple model (MODEL_MULTI). Single
  model file contains one linear classifier, while multiple model
  contains multiple classifiers.
  
  Multiple model file has total number of models in the first column,
  followed by the model ID and model weight. The rest of information
  is the same as in single model file.

  The function sets the model type file in 'type'. For multiple model
  file, the number of models is returnes in 'nmodels', and model
  weights are returned in 'weights'.

  Each element's qmodel->wmx is a 'rows' by 'columns' matrix.

  The quadratic coefficients (qmodel->sigma) are simply the inverse
  class-conditional covariance matrices. Note that for the actual
  classification, the quadratic part has to be multiplied by -0.5.

  In case of failure, return NULL and set 'errc'. Possible errors are
  memory allocation and file access errors.
*/
struct qmodel **pqc_load_models(char *fname, int *type, int *nmodels, float **weights,
				int *rows, int *columns, int *errc)
{
  int    status;
  int    i;
  int    ntk;
  int    lc;
  int    llen;
  int    nmd;
  int    lcntr;
  int    fcol;
  int    model_cnt;
  int    length;
  int    offset;
  int    c;
  int    d;
  int    first_line;
  float  fad;
  char   *line;
  float  *wts;
  char   **tokens;
  struct qmodel **models;
  FILE   *fptr;
  
  if (fname && errc)
    {
      models = (struct qmodel **) 0;
      *errc = 0;
      lc = file_info(fname, &llen, &ntk, '\0');
      if ((lc != -1) && (ntk > 0))
	{
	  *type = detect_model(fname, &model_cnt);
	  if (*type != -1)
	    {
	      fptr = fopen(fname, "r");
	      if (fptr)
		{
		  c = lc/model_cnt;
		  if (rows)
		    *rows = c;
		  if (*type == MODEL_MULTI)
		    fcol = ntk-3;
		  else
		    fcol = ntk;
		  /*
		    Magic formula to calculate the number of features
		    d of the data set which was used to build this
		    quadratic classifier. See R & D, vol. V, page
		    10/09/2004. The number of columns in the linear
		    component of the classifier equals d+1.
		  */
		  fad = 8*fcol+1;
		  fad = 0.5*(sqrt(fad)-3);
		  d = fad;
		  if (columns)
		    *columns = d+1;
		  length = (d+1)*d/2;
		  status = pqc_load_models_alloc(&models, &wts, &line, model_cnt,
						 c, d+1, length, llen);
		  if (!status)
		    {
		      nmd = 0;
		      lcntr = 0;
		      first_line = 1;
		      while (fgets(line, llen+2, fptr) && !status)
			{
			  tokens = str_tokenize(line, WHITESPACE);
			  if (tokens)
			    {
			      offset = 0;
			      if (*type == MODEL_MULTI)
				{
				  wts[nmd] = atof(tokens[1]);
				  offset = 3;
				}
			      for (i = 0; i < d+1; i++)
				models[nmd]->wmx[lcntr][i] = atof(tokens[i+offset]);
			      offset += d+1;
			      for (i = 0; i < length; i++)
				models[nmd]->sigma[lcntr][i] = atof(tokens[i+offset]);
			      if ((lcntr == c-1) && (nmd < model_cnt-1))
				{
				  /*
				    Next model.
				  */
				  nmd++;
				  lcntr = 0;
				}
			      else
				lcntr++;
			      str_free(tokens);
			    }
			  else
			    status = -1;
			}
		    }
		  else
		    *errc = errno; /* alloc() failure */
		  status = fclose(fptr);
		  if (status == -1)
		    {
		      for (i = 0; i < model_cnt; i++)
			mx_free((void **) models[i], c);
		      models = vx_free(models);
		      *errc = errno;
		    }
		  else if (weights)
		    *weights = wts;
		}
	      else
		*errc = errno;
	    }
	  else
	    *errc = errno;
	}
      else
	*errc = errno;
      if (nmodels)
	*nmodels = model_cnt;
    }
  else if (errc)
    *errc = EINVAL;
  return models;
}

/*
  Internal function for quadratic classifier. Computes the value of
  the parametric quadratic classifier discriminant function.
*/
static float pqc_compute(float *model, float *sigma, int d, float *vector, float *fv)
{
  float gx;
  
  gx = fvec_dot(model, vector, d, (int *) 0);
  gx += model[d];
  /*
    Add quadratic term.
  */
  fvec_smx(vector, sigma, d, fv);
  gx += -0.5*fvec_dot(vector, fv, d, (int *) 0);
  return gx;
}

/*
  Return predicted classification of 'vector' by parametric quadratic
  classifer stored in 'model'. It is assumed that the length of
  'vector' is ad-1.

  Return -1 in case of bad arguments or malloc() error and set
  errno. Otherwise, the result is in the [0..c-1] interval.
*/
int pqc_predict(struct qmodel *model, int c, int ad, float *vector)
{
  int   i;
  int   d;
  int   predicted_class;
  float gx;
  float gmax;
  float *fv;

  predicted_class = -1;
  gmax = -FLT_MAX;
  d = ad-1;
  if (model && (c > 0) && (ad > 0) && vector)
    {
      fv = malloc(d*sizeof(float));
      if (fv)
	{
	  for (i = 0; i < c; i++)
	    {
	      gx = pqc_compute(model->wmx[i], model->sigma[i], d, vector, fv);
	      if (gx > gmax)
		{
		  predicted_class = i;
		  gmax = gx;
		}
	    }
	  free(fv);
	}
    }
  else
    errno = EINVAL;
  return predicted_class;
}

/*
  Calculate quadratic classifier output for 'vector'. The classifier
  is defined by 'model'. It is assumed that the length of 'vector' is
  ad-1. The output is a vector of length c which contains values of
  discriminant function for each class.

  Return -1 in case of bad arguments or malloc() error and set errno.
*/
float *pqc_output(struct qmodel *model, int c, int ad, float *vector)
{
  int   i;
  int   d;
  int   predicted_class;
  float gmax;
  float *fv;
  float *output;

  output = (float *) 0;
  predicted_class = -1;
  gmax = -FLT_MAX;
  d = ad-1;
  if (model && (c > 0) && (ad > 0) && vector)
    {
      output = malloc(c*sizeof(float));
      if (output)
	{
	  fv = malloc(d*sizeof(float));
	  if (fv)
	    {
	      for (i = 0; i < c; i++)
		output[i] = pqc_compute(model->wmx[i], model->sigma[i], d, vector, fv);
	      free(fv);
	    }
	}
    }
  else
    errno = EINVAL;
  return output;
}