/*
  File name: adaboost.c
  Created by: Ljubomir Buturovic
  Created: 09/18/2002
  Purpose: implementation of Adaboost variants of principal pattern
  classification algorithms.
*/

/*
  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: adaboost.c,v 1.46 2006/05/24 06:04:12 ljubomir Exp $";

#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <string.h>
#include "lau.h"
#include "lmat.h"
#include "adaboost.h"
#include "pau.h"
#include "mlp.h"
#include "pcl_svm.h"
#include "svm.h"
#include "heap.h"

#ifndef INIT_NO_MODELS
#define INIT_NO_MODELS 10
#endif

static int boost_dataset_malloc(struct dataset **adaset, struct dataset *dset,
				float ***x, double **vector, int **vidx)
{
  int   status = 0;
  int   *vx = (int *) 0;
  int   *nd;
  double *v = (double *) 0;
  float **matrix = (float **) 0;
  struct dataset *set = (struct dataset *) 0;

  nd = ivec_clone(dset->nd, dset->c);
  if (nd)
    {
      set = dataset_new(dset->d, dset->c, nd, dset->nv, (char **) 0, (float **) 0);
      if (set != (struct dataset *) 0)
	{
	  matrix = fmx_alloc(dset->nv, dset->d);
	  if (matrix != (float **) 0)
	    {
	      v = malloc((dset->nv+1)*sizeof(double));
	      if (v)
		{
		  vx = malloc((dset->nv+1)*sizeof(int)); /* +1 for intsort() (i.e., sedgesort()) */
		  if (!vx)
		    status = -1;
		}
	      else
		status = -1;
	    }
	  else
	    status = -1;
	}
      else
	status = -1;
      if (status == 0)
	{
	  *adaset = set;
	  *x = matrix;
	  *vector = v;
	  *vidx = vx;
	}
      if (status == -1)
	{
	  mx_free((void **) matrix, dset->nv);
	  dataset_free(set);
	}
    }
  else
    status = -1;
  return status;
}

#if 0
static int boost_dataset_realloc(struct dataset **adaset, struct dataset *dset,
				 float ***x, double **vector, int **vidx)
{
  int status = -1;
  

  return status;
}
#endif

/*
  Create boost version of 'dset'. Each vector in the boosted dataset
  is chosen from 'dset' with probability 'p[i]'.
*/
static struct dataset *boost_dataset(struct dataset *dset, double *p)
{
  int    i;
  int    nx;
  int    status;
  int    current_class;
  int    offset;
  int    cnd;
  float  r;
  double csum;
  int    *vidx;
  double *dist_func;
  float  **x;
  struct dataset *adaset = (struct dataset *) 0;

  status = boost_dataset_malloc(&adaset, dset, &x, &dist_func, &vidx);
  if (status == 0)
    {
      dvec_copy(&dist_func[1], p, dset->nv);
      dist_func[0] = 0.0;
      status = dheap(dist_func, dset->nv+1, (int **) 0);
      if (status == 0)
	{
	  /*
	    Convert dist_func to distribution function.
	  */
	  csum = 0.0;
	  for (i = 0; i < dset->nv; i++)
	    {
	      csum += dist_func[i+1];
	      dist_func[i+1] = csum;
	    }
	  for (i = 0; i < dset->nv; i++)
	    {
	      /*
		Generate pseudo-random number between 0 and 1. If
		it lands between dist_func[j] and dist_func[j+1],
		add dset vector j to adaset.  
		
		But: since we have to order these vectors according to
		their classes - first we have to store these indices,
		determine their classes, and then copy the
		corresponding vectors.
	      */
	      r = float_rand();
	      vidx[i] = dloc((double) r, dist_func, dset->nv+1);
	      r += 1; /* for breakpoints */
	    }
	  intsort(vidx, dset->nv);
	  current_class = 0;
	  nx = 0;
	  cnd = dset->nd[0];
	  for (i = 0; i < dset->nv; i++)
	    {
	      nx++;
	      offset = vidx[i];
	      if (offset >= cnd)
		{
		  adaset->nd[current_class] = nx;
		  current_class++;
		  cnd += dset->nd[current_class];
		  nx = 0;
		}
	      x[i] = fvec_clone(dset->x[offset], dset->d);
	    }
	  adaset->nd[current_class] = nx;
	}
      free(dist_func);
      free(vidx);
      adaset->x = x;
    }
  return adaset;
}

/*
  Train a classifier on 'dset', using AdaBoost algorithm and
  'method'. Optional parameters for 'method' are expected to be placed
  in 'options', which should be pointer to a structure of the
  method-specific parameters. Return array of 'nmodels' weighted
  classifiers, with weights in 'weights'. The function saves the
  resulting classifier in 'fname'.

  The function implements the logic and follows notation of: Amanda
  J. C. Sharkey (Ed.), Combining Artificial Neural Nets, Chapter 2,
  Combining Predictors, Section 2.5.1. Springer, London, 1999.

  In case of error, return NULL and set 'errc'. Possible errors are:
  malloc() errors; mlp_learn() errors (for method == PALG_MLP);
  svm_train() errors (for method == PALG_SVM).
*/
void **adaboost(struct dataset *dset, int method, int *nmodels, float **weights, 
		char *fname, unsigned int seed, void *options, int *errc, FILE *fdbg)
{
  int    i;
  int    j;
  int    limit;
  int    status;
  int    mlp_nlayers;
  int    mlp_itmax;
  int    true_class;
  int    class_count;
  int    mcv;
  int    nmd;
  int    done;
  int    capacity;
  float  epsilon; /* weighted error rate in a particular step */
  float  zt;
  float  x;
  char   *fnm;
  int    *mlp_npl; 
  float  *beta; /* weight assigned to a classifier */
  float  *bt; /* beta for incorrect classification, 1 for correct */ 
  double *di; /* weight assigned to a vector */
  float  *output;
  float  **target;
  void   **models = (void **) 0;
  void   *model;
  struct mlp_options *mlp_optional;
  struct svm_problem *problem;
  struct svm_parameter *parameters;
  struct svm_node *svm_vector;
  struct dataset *adaset;
  struct dataset *cset;
  FILE   *fptr;
  char   *func = "adaboost()";
  
  status = 0;
  mlp_nlayers = 0;
  mlp_itmax = 0;
  nmd = 0;
  mlp_npl = (int *) 0;
  di = (double *) 0;
  model = (void *) 0;
  beta = (float *) 0;
  mlp_optional = (struct mlp_options *) 0;
  parameters = (struct svm_parameter *) 0;
  if (dset)
    {
      di = malloc(dset->nv*sizeof(double)); 
      bt = malloc(dset->nv*sizeof(float)); 
      if (di && bt)
	{
	  for (i = 0; i < dset->nv; i++)
	    di[i] = 1.0/dset->nv;
	  /*
	    We don't know the final number of classifiers. But we set
	    it to INIT_NO_MODELS and reallocate later, if needed.
	  */
	  capacity = INIT_NO_MODELS;
	  beta = malloc(capacity*sizeof(float)); 
	  models = malloc(capacity*sizeof(void *));
	  if (models != (void **) 0)
	    {
	      status = 0;
	      fptr = fopen(fname, "w");
	      if (fptr != (FILE *) 0)
		{
		  if ((method == PALG_MLP) || (method == PALG_ADABOOST_MLP))
		    {
		      mlp_optional = (struct mlp_options *) options;
		      mlp_nlayers = mlp_optional->nlayers;
		      mlp_npl = mlp_optional->npl;
		      mlp_itmax = mlp_optional->itmax; 
		    }
		  else if ((method == PALG_SVM) || (method == PALG_ADABOOST_SVM))
		    parameters = (struct svm_parameter *) options;
		  cset = (struct dataset *) 0;
		  done = 0;
		  nmd = 0;
		  if (*nmodels == 0)
		    limit = INT_MAX;
		  else
		    limit = *nmodels;
		  for (i = 0; (i < limit) && (status == 0) && (done == 0); i++)
		    {
		      if (i == 0)
			{
			  adaset = dataset_clone(dset);
			  adaset->prediction = malloc(dset->nv*sizeof(int));
			}
		      else
			adaset = boost_dataset(cset, di);
		      dataset_free(cset);
		      cset = adaset;
		      if ((method == PALG_MLP) || (method == PALG_ADABOOST_MLP))
			{
			  target = mlp_target(adaset->c, adaset->nd);
			  if (method == PALG_MLP)
			    fnm = fname;
			  else
			    fnm = PCP_ADA;
			  model = mlp_learn(mlp_optional->opt_method, adaset->x, adaset->nv,
					    adaset->nd, adaset->d, target, mlp_nlayers,
					    mlp_npl, mlp_itmax, mlp_optional->range,
					    mlp_optional->eta, mlp_optional->mu, stdout,
					    0, fnm, seed, errc, fdbg);
			  mx_free((void **) target, adaset->c);
			}
		      else if ((method == PALG_SVM) || (method == PALG_ADABOOST_SVM))
			{
			  problem = create_problem(adaset); 
			  model = svm_train(problem, parameters);
			}
		      epsilon = 0.0;
		      true_class = 0;
		      fvec_set(bt, dset->nv, 1.0);
		      class_count = adaset->nd[0];
		      mcv = 0;
		      for (j = 0; j < adaset->nv; j++)
			{
			  if (j == class_count)
			    {
			      true_class++;
			      class_count += adaset->nd[true_class];
			    }
			  if ((method == PALG_SVM) || (method == PALG_ADABOOST_SVM))
			    {
			      svm_vector = create_svm_vector(adaset->x[j], adaset->d);
			      adaset->prediction[j] = svm_predict(model, svm_vector)-1;
			      free(svm_vector);
			    }
			  else if ((method == PALG_MLP) || (method == PALG_ADABOOST_MLP))
			    {
			      output = mlp_output(model, adaset->x[j]);
			      adaset->prediction[j] = fvec_valmax(output, adaset->c, (float *) 0);
			      free(output);
			    }
			  if (adaset->prediction[j] != true_class)
			    {
			      epsilon += di[j];
			      mcv++;
			    }
			}
		      beta[i] = epsilon/(1.0-epsilon);
		      true_class = 0;
		      class_count = adaset->nd[0];
		      for (j = 0; j < adaset->nv; j++)
			{
			  if (j == class_count)
			    {
			      true_class++;
			      class_count += adaset->nd[true_class];
			    }
			  if (adaset->prediction[j] != true_class)
			    bt[j] = beta[i];
			  di[j] = di[j]*bt[j];
			}
		      zt = dvec_sum(di, dset->nv);
		      zt = 1.0/zt;
		      for (j = 0; j < dset->nv; j++)
			di[j] = di[j]*zt;
		      if (fdbg != (FILE *) 0)
			{
			  if (beta[i] == 0.0)
			    x = 1.0;
			  else
			    x = -log(beta[i]);
			  fprintf(fdbg, "%s; model: %5d; epsilon: %12.5g; beta[%5d]: %12.5g; weight[%5d]: %12.5g; mcv: %5d\n", 
				  func, i, epsilon, i, beta[i], i, x, mcv);
			  fflush(fdbg);
			}
		      if (model)
			{
			  models[i] = model;
			  if ((method == PALG_SVM) || (method == PALG_ADABOOST_SVM))
			    status = save_svm(fptr, model, i+1, di[i]);
			  else if ((method == PALG_MLP) || (method == PALG_ADABOOST_MLP))
			    status = mlp_write(fptr, model, MLP_MODE_APPEND, i+1, di[i]);
			  if (status == 0)
			    {
			      nmd++;
			      if (nmd >= capacity)
				{
				  capacity += capacity;
				  beta = realloc(beta, capacity*sizeof(float)); 
				  models = realloc(models, capacity*sizeof(void *));
				}
			    }
			  if (epsilon == 0.0)
			    done = 1;
			}
		      else
			status = -1;
		    }
		  fclose(fptr);
		}
	      else
		status = -1;
	    }
	}
    }
  if (status == 0)
    {
      for (i = 0; i < nmd; i++)
	beta[i] = -log(beta[i]);
      *weights = beta;
      *nmodels = nmd;
    }
  else
    free(di);
  return models;
}

int boosting_nmodels(FILE *indev, FILE *outdev)
{
  int   nmodels;
  int   min_range;
  int   int_default;
  char  *msg;

  min_range = 0;
  int_default = 0;
  msg = malloc((PCP_QLEN+1)*sizeof(char));
  sprintf(msg, PCP_UMSG_NBOOST, int_default);
  nmodels = input_integer(indev, outdev, msg, PCP_QLEN, &int_default, &min_range, 
			 (int *) 0);
  free(msg);
  return nmodels;
}

/*
  Accept input parameters and pass them to adaboost learning function
  adaboost().

  The function provides Adaboost driver for two learning algorithms,
  MLP (*method == PALG_MLP) and SVM (*method == PALG_SVM). The resulting
  classifier is saved in user-provided file 'fname'.


  In case of error, set 'errc'. If error is file access error, set
  'xname'.
*/
void p_boost_learn(int *method, int *errc, char **xname, int *dbg)
{
  int   status;
  int   nmodels;
  int   mlp_nlayers;
  int   mlp_itmax;
  int   opt_method;
  int   *mlp_npl;
  unsigned int seed;
  float mlp_range;
  float eta;
  float mu;
  float *weights;
  char  *fname;
  void  **models;
  FILE  *outdev;
  void  *options;
  struct svm_parameter *parameters;
  FILE  *fdbg = (FILE *) 0;

  models = (void **) 0;
  outdev = stdout;
  fflush(outdev);
  if (*dbg > 0)
    fdbg = fopen(PCP_DBG, "a");
  if ((*method == PALG_MLP) || (*method == PALG_ADABOOST_MLP))
    {
      status = input_mlp(outdev, tds->d, tds->c, tds->nv, (int *) 0, &mlp_nlayers,
			 &mlp_npl, &mlp_itmax, &mlp_range, &opt_method, &eta, &mu,
			 (int *) 0, 0, (int *) 0, &seed, &fname);
      options = (struct mlp_options *) calloc(1, sizeof(struct mlp_options));
      ((struct mlp_options *) options)->npl = mlp_npl;
      ((struct mlp_options *) options)->nlayers = mlp_nlayers;
      ((struct mlp_options *) options)->itmax = mlp_itmax;
      ((struct mlp_options *) options)->range = mlp_range;
      ((struct mlp_options *) options)->opt_method = opt_method;
      ((struct mlp_options *) options)->eta = eta;
      ((struct mlp_options *) options)->mu = mu;
    }
  else
    {
      parameters = calloc(1, sizeof(struct svm_parameter));
      status = input_svm(outdev, tds->c, tds->nv, (struct svm_problem *) 0,
			 parameters, &fname, (int *) 0, 0, 1, PCP_SVM_K_NONE, 
			 (int *) 0);
      seed = input_seed(stdin, outdev);
      options = (struct svm_parameter *) parameters;
    }
  nmodels = boosting_nmodels(stdin, outdev);
  if (status == 0)
    {
      inverse_video();
      srand(seed);
      models = adaboost(tds, *method, &nmodels, &weights, fname, seed, options, errc, fdbg);
    }
  if (models == (void **) 0)
    *xname = strdup(fname);
  if (*errc == 0)
    pwait();
}