/*
  File name: p_svm.c
  Created by: Ljubomir Buturovic
  Created: 08/23/2002
  Purpose: SVM model selection, learning and prediction, using the
  LIBSVM library.
*/

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

#include <stdlib.h>
#include <errno.h>
#include <unistd.h>
#include <string.h>
#include <float.h>
#include <math.h>
#include "simplex.h"
#include "lau.h"
#include "lmat.h"
#include "pcl_svm.h"
#include "pau.h"
#include "svm.h"
#include "bagging.h"
#include "xlearn.h"
#include "xpar.h"

/*
  This function collects SVM learning parameters from the user, and
  then calls svm_train() which performs the SVM training.

  In case of successful completion, set 'errc' to 0. In case of error,
  set errc to 'errno'. If error is file access error, set 'file'.
*/
void p_svm_learn(int *errc, char **xname, int dbg)
{
  int    status;
  char   *fname;
  FILE   *outdev;
  FILE   *fdbg = (FILE *) 0;
  struct svm_problem *problem;
  struct svm_parameter *parameters;
  struct svm_model *model;

  *errc = 0;
  fname = (char *) 0; /* avoid compiler warning */
  outdev = stdout;
  fflush(outdev);
  if (dbg)
    fdbg = fopen(PCP_DBG, "a");
  problem = create_problem(tds);
  if (problem)
    {
      parameters = calloc(1, sizeof(struct svm_parameter));
      if (parameters)
	{
	  status = input_svm(outdev, tds->c, tds->nv, problem, parameters,
			     (char **) 0, (int *) 0, 0, 1, PCP_SVM_K_NONE,
			     (int *) 0);
	  parameters->probability = 1;
	  if (!status)
	    {
	      fname = input_filename(PCP_UMSG_SVM_FNAME, PCP_SVM, outdev);
	      inverse_video();
	      model = svm_train(problem, parameters);
	      if (model == (struct svm_model *) 0)
		status = -1;
	      else
		status = svm_save_model(fname, model);
	      if (status == -1)
		{
		  *errc = errno;
		  *xname = strdup(fname);
		}
	    }
	  else
	    {
	      *errc = errno;
	      *xname = strdup(fname);
	    }
	  free(parameters);
	}
      else
	*errc = errno;
    }
  else
    *errc = errno;
  reset_video();
}

#define PCP_NU_SVC   1
#define PCP_C_SVC    2

/*
  Kernel types are internally coded starting from 1 (because 0 is
  reserved for the default kernel type). This function converts the
  internal kernel types to LIBSVM kernel types.
*/
int convert_kernel_type(int type)
{
  int kernel_type;

  kernel_type = PCP_SVM_K_RBF;
  if (type == PCP_SVM_K_LINEAR)
    kernel_type = LINEAR;
  else if (type == PCP_SVM_K_POLY)
    kernel_type = POLY;
  else if (type == PCP_SVM_K_RBF)
    kernel_type = RBF;
  else if (type == PCP_SVM_K_SIGMOID)
    kernel_type = SIGMOID;
  return kernel_type;
}

int input_svm_type(char *msg, FILE *outdev)
{
  int min_range;
  int max_range;
  int type;

  min_range = PCP_NU_SVC;
  max_range = PCP_C_SVC;
  sprintf(msg, PCP_UMSG_SVM_TYPE, PCP_NU_SVC, PCP_C_SVC, PCP_C_SVC);
  type = input_integer(stdin, outdev, msg, PCP_QLEN, &max_range, &min_range, &max_range);
  if (type == PCP_NU_SVC)
    type = NU_SVC;
  else if (type == PCP_C_SVC)
    type = C_SVC;
  return type;
}

/*
  Accept kernel type from stdin.
*/
int input_kernel_type(char *msg, FILE *outdev)
{
  int kernel_type;
  int min_range;
  int max_range;
  int dflt;
  
  min_range = PCP_SVM_K_LINEAR;
  max_range = PCP_SVM_K_SIGMOID;
  dflt = PCP_SVM_K_RBF;
  sprintf(msg, PCP_UMSG_KERNEL_TYPE, PCP_SVM_K_LINEAR, PCP_SVM_K_POLY, PCP_SVM_K_RBF,
	  PCP_SVM_K_SIGMOID, dflt);
  kernel_type = input_integer(stdin, outdev, msg, PCP_QLEN, &dflt, &min_range, &max_range);
  return kernel_type;
}

/*
  Accept kernel type from a limited selection (only the types
  supported by the model selection functions pcp_svm_grid() and
  pcp_svm_simplex()).
*/
int input_kernel(char *msg, FILE *outdev)
{
  int kernel_type;
  int dflt;
  int choices[2];
  int cntr;
  
  cntr = 0;
  choices[cntr++] = PCP_SVM_K_LINEAR;
  choices[cntr++] = PCP_SVM_K_RBF;
  dflt = PCP_SVM_K_RBF;
  sprintf(msg, PCP_UMSG_KERNEL, PCP_SVM_K_LINEAR, PCP_SVM_K_RBF, dflt);
  kernel_type = input_choice(stdin, outdev, msg, &dflt, choices, cntr);
  return kernel_type;
}

void input_class_costs(struct svm_parameter *parameters, int c)
{
  int   min_range;
  int   max_range;
  int   dflt;
  int   ixd;
  int   i;
  float fmin;
  float fdflt;
  float fxd;
  char  *msg;

  min_range = 0;
  max_range = 1;
  dflt = 0;
  ixd = input_integer(stdin, stdout, PCP_UMSG_CCOSTS, PCP_QLEN, &dflt, &min_range, &max_range);
  if (ixd == 1)
    {
      parameters->nr_weight = c;
      parameters->weight_label = malloc(c*sizeof(int));
      parameters->weight = malloc(c*sizeof(double));
      msg = malloc((PCP_QLEN+1)*sizeof(char));
      for (i = 0; i < c; i++)
	{
	  fmin = 0.0;
	  fdflt = 1.0;
	  sprintf(msg, PCP_UMSG_CLASS_COST, i+1, fdflt);
	  fxd = input_float(stdin, stdout, msg, PCP_QLEN, &fdflt, &fmin, (float *) 0);
	  parameters->weight_label[i] = i+1;
	  parameters->weight[i] = fxd;
	}
      free(msg);
    }
}

/*
  Set the SVM parameters we decided to keep as defaults.
*/
void svm_param_defaults(void *svm_parameters)
{
  float fdflt;
  struct svm_parameter *svm_params;

  svm_params = (struct svm_parameter *) svm_parameters;
  svm_params->probability = PCP_SVM_DFLT_PROB;
  svm_params->cache_size = PCP_SVM_DFLT_CACHE;
  if (svm_params->svm_type == NU_SVC) /* default values per LIBSVM README file */
    fdflt = PCP_SVM_DFLT_EPS_NU;
  else
    fdflt = PCP_SVM_DFLT_EPS_C;
  svm_params->eps = fdflt; /* decided to always use default values. ljb, 03/06/2005 */
  svm_params->shrinking = 1; /* Decided to always use the shrinking heuristics */
}

/*
  Compute max. feasible value of nu for NU-SVM. Use Proposition 3 in
  Section 6 of
  http://www.csie.ntu.edu.tw/~cjlin/papers/nusvmtutorial.pdf.
*/
float get_svm_nu_max(int *nd, int c)
{
  int   i;
  int   j;
  int   n1;
  int   n2;
  int   nx;
  float nu;
  float numax;

  numax = 1.0;
  if (nd)
    {
      for (i = 0; i < c-1; i++)
	for (j = i+1; j < c; j++)
	  {
	    n1 = nd[i];
	    n2 = nd[j];
	    nx = n1;
	    if (n2 < n1)
	      nx = n2;
	    nu = 2.0*nx/(n1+n2);
	    if (nu < numax)
	      numax = nu;
	  }
    }
  return numax;
}

/*
  Compute max. feasible value of nu for `nxval' NU-SVM
  cross-validation experiments.
*/
float get_svm_nu_xmax(int *nd, int c, int nxval)
{
  int   i;
  int   j;
  int   n1;
  int   n2;
  float nu;
  float numax;
  int   ndx[2];

  numax = 1.0;
  if (nd)
    {
      for (i = 0; i < c-1; i++)
	for (j = i+1; j < c; j++)
	  {
	    n1 = nd[i];
	    n2 = nd[j];
	    
	    /*
	      These steps go through all four possibilities of
	      learning subset sizes for classes i and j, and pick the
	      lowest numax of the four possibilities.
	    */
	    ndx[0] = n1-n1/nxval;
	    ndx[1] = n2-n2/nxval;
	    nu = get_svm_nu_max(ndx, 2);
	    if (nu < numax)
	      numax = nu;

	    ndx[0] -= 1;
	    nu = get_svm_nu_max(ndx, 2);
	    if (nu < numax)
	      numax = nu;
	    
	    ndx[0] += 1;
	    ndx[1] -= 1;
	    nu = get_svm_nu_max(ndx, 2);
	    if (nu < numax)
	      numax = nu;
	    
	    ndx[0] -= 1;
	    nu = get_svm_nu_max(ndx, 2);
	    if (nu < numax)
	      numax = nu;
	  }
    }
  return numax;
}

/*
  Read SVM parameters from stdin. Check consistency with 'problem'.

  If `full_par' is not 0, read C/nu, gamma, degree and coef0
  parameters, depending on the SVM type. Otherwise, if `ktype' is
  PCP_SVM_K_NONE, read only the kernel type. Otherwise, we don't need
  these parameter values (presumably because we are called from a
  function which will optimize them).
  
  If `method' is not NULL, read the ensemble variant of PALG_SVM
  (i.e., PALG_BAG_SVM or PALG_ADABOOST_SVM).
*/
int input_svm(FILE *outdev, int c, int nvec, struct svm_problem *problem,
	      struct svm_parameter *parameters, char **fname, int *nxval, 
	      int max_nxval, int full_par, int ktype, int *ensemble_method)
{
  int   i;
  int   status = 0;
  int   ixd;
  int   kernel_type;
  int   min_range;
  int   max_range;
  int   done;
  float fxd;
  float fmin;
  float fmax;
  float fdflt;
  char  *pcheck;
  char  *msg;
  FILE  *fptr;
  
  clear_screen();
  cursor_on();
  if (ensemble_method)
    *ensemble_method = input_ensemble_method(PALG_SVM);
  msg = malloc((PCP_QLEN+1)*sizeof(char));
  done = 0;
  while (done == 0)
    {
      parameters->svm_type = input_svm_type(msg, outdev);

      if (full_par)
	{
	  kernel_type = input_kernel_type(msg, outdev);
	  ixd = convert_kernel_type(kernel_type);
	  parameters->kernel_type = ixd;
	}
      else if (ktype != PCP_SVM_K_NONE)
	{
	  kernel_type = input_kernel(msg, outdev);
	  parameters->kernel_type = convert_kernel_type(kernel_type);
	}
      else
	{
	  kernel_type = PCP_SVM_K_RBF;
	  parameters->kernel_type = convert_kernel_type(kernel_type);
	}

      if (c > 2)
	{
	  min_range = 0;
	  max_range = 1;
	  sprintf(msg, PCP_UMSG_INDET, max_range, min_range, min_range);
	  /* TBD-indet: disable indeterminate-region multi-class problem algorithm for now */
	  /*parameters->indet = input_integer(stdin, outdev, msg, PCP_QLEN, &dflt, &min_range, &max_range);*/
	}

      if (parameters->svm_type == NU_SVC && full_par)
	{
	  fmin = 0.0;
	  fmax = get_svm_nu_max(tds->nd, tds->c);
	  if (PCP_SVM_DFLT_NU < fmax)
	    fdflt = PCP_SVM_DFLT_NU;
	  else
	    fdflt = 0.5*fdflt;
	  sprintf(msg, "%s (%4.2f..%4.2f) [%4.2f]:", PCP_UMSG_NU, fmin, fmax, fdflt);
	  fxd = input_float(stdin, outdev, msg, PCP_QLEN, &fdflt, &fmin, &fmax);
	  parameters->nu = fxd;
	}
      if (parameters->svm_type == C_SVC && full_par)
	{
	  fdflt = PCP_SVM_DFLT_C;
	  fmin = 0.0;
	  sprintf(msg, PCP_UMSG_COST, fdflt);
	  fxd = input_float(stdin, outdev, msg, PCP_QLEN, &fdflt, &fmin, (float *) 0);
	  parameters->C = fxd;
	}
      if (parameters->svm_type == C_SVC)
	input_class_costs(parameters, c);
      if (full_par)
	{
	  if (parameters->kernel_type == POLY)
	    {
	      fdflt = PCP_SVM_DFLT_DEGREE;
	      fmin = 1.0;
	      sprintf(msg, PCP_UMSG_DEGREE, fdflt);
	      fxd = input_float(stdin, outdev, msg, PCP_QLEN, &fdflt, &fmin, (float *) 0);
	      parameters->degree = fxd;
	    }
	  if ((parameters->kernel_type == POLY) || (parameters->kernel_type == RBF) || 
	      (parameters->kernel_type == SIGMOID))
	    {
	      /*fdflt = 1.0/(tds->d);*/
	      fdflt = PCP_SVM_DFLT_GAMMA;
	      sprintf(msg, PCP_UMSG_GAMMA, fdflt);
	      fxd = input_float(stdin, outdev, msg, PCP_QLEN, &fdflt, (float *) 0, (float *) 0);
	      parameters->gamma = fxd;
	    }
	  if ((parameters->kernel_type == POLY) || (parameters->kernel_type == SIGMOID))
	    {
	      fdflt = PCP_SVM_DFLT_COEF0;
	      sprintf(msg, PCP_UMSG_COEF0, fdflt);
	      fxd = input_float(stdin, outdev, msg, PCP_QLEN, &fdflt, (float *) 0, (float *) 0);
	      parameters->coef0 = fxd;
	    }
	}
      svm_param_defaults(parameters);
      pcheck = (char *) 0;
      if (problem)
	pcheck = (char *) svm_check_parameter(problem, parameters);
      if (pcheck)
	{
	  sprintf(msg, PCP_UMSG_BAD_SVM_PARAMS, pcheck);
	  print_line(outdev, msg, PCP_QLEN);
	  fprintf(outdev, "\n");
	}
      else
	done = 1;
    }
  if (nxval)
    *nxval = input_nxval(stdin, outdev, max_nxval);
  if (fname)
    {
      if (ensemble_method)
	*fname = input_filename(PCP_UMSG_SVM_FNAME, PCP_SVM, outdev);
      else
	*fname = input_filename(PCP_UMSG_XSV, PCP_XSV, outdev);
    }
  if (fname)
    {
      fptr = fopen(*fname, "w");
      if (fptr)
	{
	  fclose(fptr);
	  unlink(*fname);
	}
      else
	{
	  i = errno; /* for debugging */
	  status = -1;
	}
    }
  vx_free(msg);
  return status;
}

/*
  Save training data set in LIBSVM data format. See LIBSVM
  documentation for the format description.

  In case of success, set 'errc' to 0, otherwise set it to 'errno'. If
  error is file access error, store the name of the offending file in
  'xname'.
*/
void p_svm_save(int *errc, char **xname)
{
  int    i;
  int    j;
  int    type;
  int    min_range;
  int    max_range;
  char   *fname;
  char   *msg;
  struct svm_problem *problem;
  struct svm_node *node;
  FILE *fptr;

  *errc = 0;
  clear_screen();
  cursor_on();
  msg = malloc((PCP_QLEN+1)*sizeof(char));
  min_range = 0;
  max_range = 1;
  sprintf(msg, " Convert training (%d) or test (%d) dataset [%d]:", 
	  min_range, max_range, min_range);
  type = input_integer(stdin, stdout, msg, PCP_QLEN, &min_range, &min_range, 
		       &max_range);
  free(msg);
  fname = input_filename(PCP_UMSG_ONAME, PCP_LVM, stdout);
  if (type == 0)
    problem = create_problem(tds);
  else
    problem = create_problem(teds);
  if (problem)
    {
      fptr = fopen(fname, "w");
      if (fptr)
	{
	  for (i = 0; i < problem->l; i++)
	    {
	      fprintf(fptr, "%d ", (int) problem->y[i]);
	      node = problem->x[i];
	      j = 0;
	      while (node[j].index != -1)
		{
		  fprintf(fptr, "%d:%g ", j+1, node[j].value);
		  j++;
		}
	      fprintf(fptr, "\n");
	    }
	  fclose(fptr);
	}
      else
	{
	  *errc = errno;
	  *xname = strdup(fname);
	}
    }
  else
    *errc = errno;
}

/*
  Classify test data set using SVM classifier stored in 'svm_fname'
  (the file name provided by user), and display classification
  results.

  In case of successful completion, set 'errc' to 0. Otherwise, set
  errc to 'errno'. In case of file error, store the name of the file
  causing the error in 'xname'.
*/
void p_svm_predict(int *errc, char **xname)
{
  int   i;
  int   j;
  int   verbose;
  int   min_range;
  int   max_range;
  int   predicted_class;
  int   nmodels;
  int   prob_flag;
  float *predictions;
  char  *svm_fname;
  char  *output_fname;
  float *weights;
  struct svm_model *model;
  struct svm_node  *svm_vector;
  struct svm_model **models;
  FILE  *outdev;
  FILE  *svm_fptr;
  FILE  *output_fptr;

  *errc = 0;
  verbose = -1;
  clear_screen();
  outdev = stdout;
  cursor_on();
  svm_fname = input_filename(PCP_UMSG_SVM_FNAME, PCP_SVM, outdev);
  svm_fptr = fopen(svm_fname, "r");
  if (svm_fptr)
    {
      fclose(svm_fptr);
      output_fname = strdup(PCP_RCL);
      output_fptr = fopen(output_fname, "w");
      if (output_fptr)
	{
	  models = pcl_svm_load(svm_fname, &nmodels, &weights);
	  if (models)
	    {
	      predictions = calloc(tds->c, sizeof(float));
	      teds->prediction = malloc(teds->nv*sizeof(int));
	      teds->ps = dmx_alloc(teds->nv, teds->c);
	      for (i = 0; i < teds->nv; i++)
		{
		  svm_vector = create_svm_vector(teds->x[i], teds->d);
		  fvec_set(predictions, tds->c, 0.0);
		  for (j = 0; j < nmodels; j++)
		    {
		      model = models[j];
		      /*
			Check if the model contains probability
			information.
		      */
		      prob_flag = svm_check_probability_model(model);
		      if (prob_flag)
			predicted_class = svm_predict_probability(model, svm_vector, 
								  teds->ps[i])-1;
		      else
			predicted_class = svm_predict(model, svm_vector)-1;
		      predictions[predicted_class] += weights[j];
		    }
		  predicted_class = fvec_argmax(predictions, tds->c);
		  free(svm_vector);
		  teds->prediction[i] = predicted_class;
		  write_rcl(output_fptr, i, teds, tds);
		}
	      free(predictions);
	      for (j = 0; j < nmodels; j++)
		svm_destroy_model(models[j]);
	      free(models);
	      fclose(output_fptr);
	      if (teds->c == tds->c)
		{
		  min_range = 0;
		  max_range = 1;
		  verbose = input_integer(stdin, outdev, OUTPUT_MSG, PCP_QLEN,
					  &min_range, &min_range, &max_range);
		}
	      /*
		Display results.
	      */
	      predict_disp(teds, verbose, PALG_SVM);
	      pwait();
	    }
	  else
	    {
	      *errc = PERR_UNRECOGNIZED_SVM;
	      *xname = strdup(svm_fname);
	      fclose(output_fptr);
	      unlink(output_fname);
	    }
	}
      else
	{
	  *errc = errno;
	  *xname = strdup(output_fname);
	}
    }
  else
    {
      *errc = errno;
      *xname = strdup(svm_fname);
    }
}

#define PCP_GAMMA_LEFT         -8
#define PCP_GAMMA_RIGHT         0
#define PCP_GAMMA_MIDDLE       -2
#define PCP_C_LEFT             -1
#define PCP_C_RIGHT            4
#define PCP_C_MIDDLE           1.5
#define PCP_NU_LEFT            0.01
#define PCP_NU_FACTOR          0.99     /* numeric safety factor for NU-SVM feasible region */
#define PCP_NSTEP              10       /* number of steps for grid search */
#define PCP_RSTEP              5        /* number of steps for refinement of the grid search */
#define PCP_RFACTOR            0.05     /* refinement factor */

/*#define PCP_GAMMA_LEFT          -6.75*/ /* for testing */
/*#define PCP_C_LEFT             3.609375*/ /* for testing */

/*
  Transplant optimal parameters to xpar_crit_parameters->options.
*/
static void set_optimal(struct xpar_crit_parameters *xpar_parameters)
{
  if (((struct svm_parameter *) xpar_parameters->options)->svm_type == C_SVC)
    ((struct svm_parameter *) xpar_parameters->options)->C = pow(10, xpar_parameters->x1);
  else
    ((struct svm_parameter *) xpar_parameters->options)->nu = xpar_parameters->x1;
  ((struct svm_parameter *) xpar_parameters->options)->gamma = pow(10, xpar_parameters->x2);
}

static struct xpar_crit_parameters *init_xpar(struct dataset *dset, FILE *fdbg)
{
  struct xpar_crit_parameters *xpar_params;
  
  xpar_params = calloc(1, sizeof(struct xpar_crit_parameters));
  if (xpar_params)
    {
      xpar_params->options = calloc(1, sizeof(struct svm_parameter));
      if (xpar_params->options)
	{
	  xpar_params->dset = dset;
	  xpar_params->outdev = stdout;
	  xpar_params->fdbg = fdbg;
	  xpar_params->vid = 0;
	  xpar_params->eval = FLT_MAX;
	  xpar_params->method = PALG_SVM;
	  xpar_params->nmodels = 1;
	  svm_param_defaults(xpar_params->options);
	}
      else
	vx_free(xpar_params);
    }
  return xpar_params;
}

/*
  Calculate criterion function values in vertices of `smx'.
*/
static int calc_simplex(int type, float *fval, float **smx, int ndim, void *xpar_parameters, 
			int triangle, float *ftrv, int *errc)
{
  int i;
  int status;
  
  status = 0;
  if (triangle == 0)
    {
      /*
	First attempt, do the `base' triangle.
      */
      if (type == C_SVC)
	{
	  smx[0][0] = PCP_C_LEFT;  /* C, log10-scale */
	  smx[1][0] = PCP_C_RIGHT;
	  smx[2][0] = PCP_C_MIDDLE;
	}
      else
	{
	  smx[0][0] = PCP_NU_LEFT;  /* nu */
	  smx[1][0] = get_svm_nu_max(tds->nd, tds->c);
	  smx[2][0] = 0.5*(smx[0][0]+smx[1][0]);
	}
      smx[0][1] = PCP_GAMMA_LEFT;  /* gamma, log10-scale */
      smx[1][1] = PCP_GAMMA_LEFT;
      smx[2][1] = PCP_GAMMA_MIDDLE;
      fval[0] = xpar_func(smx[0], ndim, 0, xpar_parameters, errc);
      if (*errc == 0)
	{
	  ftrv[0] = fval[0];
	  for (i = 1; (i < ndim+1) && !status; i++)
	    {
	      ((struct xpar_crit_parameters *) xpar_parameters)->vid = i;
	      inverse_video();
	      fval[i] = xpar_func(smx[i], ndim, 0, xpar_parameters, errc);
	      if (*errc != 0)
		status = -1;
	      else
		ftrv[i] = fval[i];
	    }
	}
      else
	status = -1;
    }
  else if (triangle == 1)
    {
      /*
	Next. Simplex vertex 0 moves diagonally across, the other two
	remain.
      */
      if (type == C_SVC)
	smx[0][0] = PCP_C_RIGHT;  
      else
	smx[0][0] = get_svm_nu_max(tds->nd, tds->c);
      smx[0][1] = PCP_GAMMA_RIGHT; 
      ((struct xpar_crit_parameters *) xpar_parameters)->vid = 0;
      fval[0] = xpar_func(smx[0], ndim, 0, xpar_parameters, errc);
      if (*errc == 0)
	{
	  ftrv[0] = fval[0];
	  fval[1] = ftrv[1];
	  fval[2] = ftrv[2];
	}
      else
	status = -1;
    }
  else if (triangle == 2)
    {
      /*
	Next. Simplex vertex 1 moves diagonally across, the other two
	remain.
      */
      if (type == C_SVC)
	smx[1][0] = PCP_C_LEFT;  
      else
	smx[1][0] = PCP_NU_LEFT; 
      smx[1][1] = PCP_GAMMA_RIGHT; 
      ((struct xpar_crit_parameters *) xpar_parameters)->vid = 1;
      fval[1] = xpar_func(smx[1], ndim, 0, xpar_parameters, errc);
      if (*errc == 0)
	{
	  ftrv[1] = fval[1];
	  fval[0] = ftrv[0];
	  fval[2] = ftrv[2];
	}
      else
	status = -1;
    }
  else if (triangle == 3)
    {
      /*
	Last triangle. Simplex vertex 0 moves diagonally across, the
	other two remain.
      */
      if (type == C_SVC)
	smx[0][0] = PCP_C_LEFT;  
      else
	smx[0][0] = PCP_NU_LEFT; 
      smx[0][1] = PCP_GAMMA_LEFT; 
      ((struct xpar_crit_parameters *) xpar_parameters)->vid = 0;
      fval[0] = xpar_func(smx[0], ndim, 0, xpar_parameters, errc);
      if (*errc == 0)
	{
	  ftrv[0] = fval[0];
	  fval[1] = ftrv[1];
	  fval[2] = ftrv[2];
	}
      else
	status = -1;
    }
  return status;
}

/*
  Optimize SVM parameters using Simplex algorithm.
*/
void pcp_svm_simplex(int *errc, int dbg, char **xname)
{
  int   ndim;
  int   itmax;
  int   i;
  int   iter;
  int   status;
  int   ex;
  int   type;
  int   done;
  int   triangle;
  int   dflt;
  int   idx;
  float ftol = 1e-6;
  float val;
  char  *msg;
  float *fval;
  float **smx;
  FILE  *outdev;
  struct xpar_crit_parameters *xpar_parameters;
  struct svm_parameter *svm_parm;
  float ftrv[5];
  FILE  *fdbg = (FILE *) 0;

  iter = 0;
  if (dbg > 0)
    fdbg = fopen(PCP_DBG, "a");
  else
    fdbg = (FILE *) 0;
  ndim = 2;
  status = smplx_alloc(&fval, &smx, ndim);
  if (!status)
    {
      outdev = stdout;
      xpar_parameters = init_xpar(tds, fdbg);
      idx = ivec_min(tds->nd, tds->c);
      svm_parm = calloc(1, sizeof(struct svm_parameter));
      status = input_svm(outdev, tds->c, tds->nv, (struct svm_problem *) 0, svm_parm,
			 (char **) 0, (int *) 0, idx, 0, PCP_SVM_K_NONE, (int *) 0);
      xpar_parameters->options = svm_parm;
      input_xpar(tds, xpar_parameters);
      dflt = PCP_DFLT_NITER;
      ex = 1;
      msg = malloc((PCP_QLEN+1)*sizeof(char));
      snprintf(msg, PCP_QLEN+1, PCP_UMSG_MAXIT, dflt);
      itmax = input_integer(stdin, outdev, msg, PCP_QLEN, &dflt, &ex, (int *) 0);
      free(msg);

      type = ((struct svm_parameter *) xpar_parameters->options)->svm_type;

      /*
	Compute feature transformations/subsets here, reuse them
	later.
       */
      if (xpar_parameters->dr_method != PDR_NONE)
	status = compute_dr(xpar_parameters, xpar_parameters->dr_method, 
			    xpar_parameters->idr, xpar_parameters->fscrit, fdbg, errc);
      done = 0;
      /*
	We try four initial simplexes (triangles) until we execute
	more than one iteration. ftrv holds the function values in
	triangle vertices.
       */
      triangle = 0;
      unlink(PCP_MSL);
      while (!done && !status)
	{
	  inverse_video();
	  status = calc_simplex(type, fval, smx, ndim, (void *) xpar_parameters, 
				triangle, ftrv, errc);
	  if (!status)
	    {
	      inverse_video();
	      val = simplex(smx, ndim, fval, xpar_func, (void *) xpar_parameters,
			    itmax, ftol, &iter, errc);
	      if (*errc == 0)
		{
		  if ((iter > 1) || (triangle == 3))
		    {
		      done = 1;
		      i = fvec_argmin(fval, ndim+1);
		      printf("\n");
		      viprint_line(4, 1, "Optimization completed in %7d iterations; minimum error rate: %7.2f%%.", 
				   iter, xpar_parameters->eval);
		      if (type == C_SVC)
			viprint_line(4, 1, "Optimal parameters: C = %12.5g; gamma = %12.5g", 
				     pow(10, xpar_parameters->x1), pow(10, xpar_parameters->x2));
		      else
			viprint_line(4, 1, "Optimal parameters: nu = %12.5g; gamma = %12.5g", 
				     xpar_parameters->x1, pow(10, xpar_parameters->x2));
		      if (!teds)
			pwait();
		    }
		  else
		    /*
		      If we finished in one iteration, that means no
		      progress has been made. Try one more time, with
		      a different initial simplex (i.e., next
		      triangle).
		    */
		    triangle++;
		}
	      else
		status = -1;
	    }
	}
      /*
	We have the optimal parameters - now apply to the test set, if
	defined.
      */
      if (teds)
	{
	  printf("\n");
	  viprint_line(2, 1, "Applying the optimal classifier to the test dataset...");
	  set_optimal(xpar_parameters);
	  xtest_optimal(xpar_parameters, errc, xname);
	}
      xpar_free(xpar_parameters);
      mx_free((void **) smx, ndim+1);
      vx_free(fval);
    }
}

static void pcp_grid_loop(int svm_type, int kernel_type, int *iter, float *grid_point, 
			  float step, float gleft, float gstep, int nstep1, int nstep2, 
			  int ndim, struct xpar_crit_parameters *xpar_parameters, int *errc)
{
  int   i;
  int   j;
  float fval;
  float pm; /* probability of error of the majority classifier */
  float p_error; 
  float p_error_opt; 


  pm = (float) ivec_max(xpar_parameters->dset->nd, xpar_parameters->dset->c);
  pm = 100.0*(1.0-pm/xpar_parameters->dset->nv);
  for (i = 0; (i <= nstep1) && !*errc; i++) /* C/nu loop */
    {
      grid_point[0] += step;
      grid_point[1] = gleft-gstep;
      for (j = 0; (j <= nstep2) && !*errc; j++) /* other parameter loop */
	{
	  grid_point[1] += gstep;
	  fval = xpar_func(grid_point, ndim, *iter, xpar_parameters, errc);
	  p_error = fval/pm;
	  p_error_opt = xpar_parameters->eval/pm;
	  if (!*errc)
	    {
	      if (kernel_type == PCP_SVM_K_LINEAR)
		{
		  viprint_line(5, 1, "Current: C = %10.3g; error rate: %5.1f%%; p-error: %7.3f.", 
			       pow(10, grid_point[0]), fval, p_error);
		  viprint_line(6, 1, "Optimal: C = %10.3g; error rate: %5.1f%%; p-error: %7.3f.", 
			       pow(10, xpar_parameters->x1), xpar_parameters->eval, p_error_opt);
		}
	      else if (kernel_type == PCP_SVM_K_RBF)
		{
		  if (svm_type == C_SVC)
		    {
		      viprint_line(6, 1, "Current: C = %10.3g, gamma = %10.3g; error: %5.1f%%; p-error: %7.3f.", 
				   pow(10, grid_point[0]), pow(10, grid_point[1]), fval, p_error);
		      viprint_line(6, 1, "Optimal: C = %10.3g, gamma = %10.3g; error: %5.1f%%; p-error: %7.3f.", 
				   pow(10, xpar_parameters->x1), pow(10, xpar_parameters->x2), 
				   xpar_parameters->eval, p_error_opt);
		    }
		  else
		    {
		      viprint_line(6, 1, "Current: nu = %10.3g, gamma = %10.3g; error: %5.1f%%; p-error: %7.3f.", 
				   grid_point[0], pow(10, grid_point[1]), fval, p_error);
		      viprint_line(6, 1, "Optimal: nu = %10.3g, gamma = %10.3g; error: %5.1f%%; p-error: %7.3f.", 
				   xpar_parameters->x1, pow(10, xpar_parameters->x2), xpar_parameters->eval, 
				   fval, p_error_opt);
		    }
		}
	      printf("\n");
	      (*iter)++;
	    }
	}
    }
}

/*
  SVM model selection using grid (i.e., exhaustive) search.
*/
void pcp_svm_grid(int *errc, int dbg, char **xname)
{
  int   iter;
  int   ndim;
  int   svm_type;
  int   status;
  int   nstep1;
  int   nstep2;
  int   kernel_type;
  float left;
  float right;
  float step;
  float gleft;
  float gright;
  float gstep;
  float numax;
  float grid_point[2];
  char  *msg;
  struct xpar_crit_parameters *xpar_parameters;
  FILE  *fdbg;

  if (dbg > 0)
    fdbg = fopen(PCP_DBG, "a");
  else
    fdbg = (FILE *) 0;
  status = 0;
  clear_screen();
  cursor_on();
  ndim = 2;
  xpar_parameters = init_xpar(tds, fdbg);
  msg = malloc((PCP_QLEN+1)*sizeof(char));
  svm_type = input_svm_type(msg, stdout);
  kernel_type = input_kernel(msg, stdout);
  free(msg);
  ((struct svm_parameter *) xpar_parameters->options)->svm_type = svm_type;
  ((struct svm_parameter *) xpar_parameters->options)->kernel_type = 
    convert_kernel_type(kernel_type);
  nstep1 = PCP_NSTEP;
  if (kernel_type == PCP_SVM_K_RBF)
    nstep2 = PCP_NSTEP;
  else if (kernel_type == PCP_SVM_K_LINEAR)
    {
      ndim = 1;
      nstep2 = 0;
    }
  if (svm_type == C_SVC)
    input_class_costs(xpar_parameters->options, tds->c);
  input_xpar(tds, xpar_parameters);
  if (svm_type == C_SVC)
    {
      left = PCP_C_LEFT;
      right = PCP_C_RIGHT;
    }
  else
    {
      left = PCP_NU_LEFT;
      /*
	Shrink the feasible interval for nu slightly, to avoid
	violating the constraint due to round-off error in the nu loop
	(the outer loop below, indexed by `i').
      */
      numax = PCP_NU_FACTOR*get_svm_nu_xmax(tds->nd, tds->c, xpar_parameters->nxval);
      right = numax;
    }
  step = (right-left)/PCP_NSTEP;
  gleft = PCP_GAMMA_LEFT;
  gright = PCP_GAMMA_RIGHT;
  gstep = (gright-gleft)/PCP_NSTEP;
  cursor_off();
  inverse_video();
  grid_point[0] = left-step;
  *errc = 0;
  unlink(PCP_MSL);
  iter = 1;
  /*
    Compute optimal feature transformations/subsets once, reuse them
    later.
  */
  if (xpar_parameters->dr_method != PDR_NONE)
    status = compute_dr(xpar_parameters, xpar_parameters->dr_method, 
			xpar_parameters->idr, xpar_parameters->fscrit, fdbg, errc);
  if (!status)
    {
      pcp_grid_loop(svm_type, kernel_type, &iter, grid_point, step, gleft, gstep, nstep1,
		    nstep2, ndim, xpar_parameters, errc);
      /*
	Narrow down the search around the best pair of parameters
	found so far. TBD: this code is tested and fully
	functional. However it is presently disabled since it doesn't
	seem to improve the predictions.
      */
      if (!*errc && 0)
	{
	  if (svm_type == NU_SVC)
	    {
	      nstep1 = PCP_RSTEP;
	      nstep2 = PCP_RSTEP;
	      left = (1-PCP_RFACTOR)*xpar_parameters->x1;
	      right = (1+PCP_RFACTOR)*xpar_parameters->x1;
	      step = (right-left)/nstep1;
	      if (right > numax)
		right = numax;
	    }
	  else
	    {
	      /*
		C-SVM is relatively insensitive to the choice of the C
		parameter.  Consequently, we don't tune it any
		more. Instead, we increase the precision of tuning for
		gamma (see `nstep2' variable).
	      */
	      left = xpar_parameters->x1;
	      step = 0.0;
	      nstep1 = 0;
	      nstep2 = 5*PCP_RSTEP;
	    }
	  grid_point[0] = left-step;
	  gleft = (1-PCP_RFACTOR)*xpar_parameters->x2;
	  gright = (1+PCP_RFACTOR)*xpar_parameters->x2;
	  gstep = (gright-gleft)/nstep2;
	  pcp_grid_loop(svm_type, kernel_type, &iter, grid_point, step, gleft, 
			gstep, nstep1, nstep2, ndim, xpar_parameters, errc);
	}
      /*
	We have the optimal parameters - now apply to the test set, if
	defined.
      */
      if (!*errc)
	{
	  if (teds)
	    {
	      viprint_line(2, 1, "Applying the optimal classifier to the test dataset...");
	      set_optimal(xpar_parameters);
	      xtest_optimal(xpar_parameters, errc, xname);
	    }
	  else
	    pwait();
	}
    }
  xpar_free(xpar_parameters);
}