/*
  File name: xlearn.c
  Created by: Ljubomir Buturovic
  Created: 01/03/2002
  Purpose: cross-validation for supervised learning 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.
*/

/*
  The main functions in this file are p_xlearn() and xlearn().

  p_xlearn() performs multiple cross-validated learning experiments. A
  single experiment is performed by xlearn().

  xlearn() calls xpart() to partition input data set into N subsets;
  trains a classifier on N-1 subsets and tests it on the remaining
  subset. This is repeated for each of the subsets.
*/

static char rcsid[] = "$Id: xlearn.c,v 1.211 2006/05/24 06:05:02 ljubomir Exp $";

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <errno.h>
#include <ctype.h>
#include <sys/types.h>
#include <sys/stat.h>
#include "xlearn.h"
#include "mlp.h"
#include "pcp.h"
#include "bagging.h"
#include "lau.h"
#include "lmat.h"
#include "pau.h"
#include "hash_util.h"
#include "pcl_svm.h"
#include "adaboost.h"
#include "lin.h"
#include "lind.h"
#include "knn.h"
#include "cda.h"
#include "emap.h"
#include "parametric.h"
#include "fselect.h"

/*
  Utility function for partitioning set of integers in [low..high]
  interval into 'nxval' disjoint subsets. The subsets are returned in
  'sx' and 'lx' integer arrays. Local variable 'tsx' has elements of
  the subsets, local variable 'tlx[i]' has cardinality of 'i'-th
  subset. tlx[i] is either [range/nxval] or [range/nxval]+1. So,
  tsx[0..lx[0]-1] are elements of the first subset;
  tsx[lx[0]..lx[0]+lx[1]-1] are elements of the second subset,
  etc. Addresses of 'tsx' and 'tlx' are returned in '*sx' and '*lx'.

  The partitioning is performed using pseudo-random number generator
  rand().

  Return -1 in case of malloc() error, 0 otherwise. If 'nxval' is
  greater than 'range', or less than 1, set 'sx' and 'lx' to NULL.

  NOTE: notation [low..high] designates set of integers between 'low'
  and 'high', including both 'low' and 'high'.
*/
static int xss(int low, int high, int nxval, int **sx, int **lx)
{
  int status;
  int range;
  int i;
  int j;
  int ivx;
  int jx;
  int tcx;
  int idx;
  int lsmall;
  int lbig;
  int *tsx;
  int *tlx;
  int *xc;
  float rcc;

  tsx = (int *) 0;
  tlx = (int *) 0;
  status = 0;
  range = high-low+1;
  if ((nxval >= 1) && (nxval <= range))
    {
      rcc = range/nxval;
      idx = 0;
      /*
	xc has a list of remaining integers from the 0..range-1.
	Initially xc contains all integers in the range.
      */
      xc = malloc(range*sizeof(int));
      for (j = 0; j < range; j++)
	xc[j] = j;
      
      tsx = malloc(range*sizeof(int));
      tlx = malloc(nxval*sizeof(int));
      if ((xc == (int *) 0) || (tsx == (int *) 0) || (tlx == (int *) 0))
	status = -1;
      else
	{
	  /*
	    Among the 'nxval' subsets, lsmall will have cardinality
	    [range/nxval]. The remaining lbig = nxval-lsmall subsets
	    will have one more element.
	  */
	  lbig = rcc;
	  lbig = range-lbig*nxval;
	  lsmall = nxval-lbig;
	  for (j = 0; j < lsmall; j++)
	    tlx[j] = rcc;
	  for (j = lsmall; j < nxval; j++)
	    tlx[j] = rcc+1;
	  tcx = range;
	  i = 0;
	  for (j = 0; j < range; j++)
	    {
	      /*
		Choose an integer, pseudo-randomly, and insert it into
		tsx.
	      */
	      ivx = rand_int(0, tcx-1);
	      tsx[i] = xc[ivx]+low;
	      i++;
	      for (jx = ivx; jx < tcx-1; jx++)
		xc[jx] = xc[jx+1];
	      tcx--;
	    }
	  free(xc);
	}
    }
  else
    status = -1;
  *sx = tsx;
  *lx = tlx;
  return status;
}

/*
  Partition each of the 'c' classes with 'nd[i]' vectors per class,
  into 'nxval' disjoint subsets. The resulting partition is specified
  in 'sxc' and 'lxc' arrays.
  
  Local arrays 'tsxc' and 'tlxc' have list of vectors in each subset,
  and subset cardinalities, respectively. 'tsxc[i]' and 'tlxc[i]' are
  pointers to 'sx' and 'lx' arrays, respectively, for class 'i', as
  described in comments for xss(). The addresses of 'tsxc' and 'tlxc'
  are returned in '*sxc' and '*lxc'.

  The function may be used for cross-validation experiments, to return
  all subsets. Then the subsets may be analyzed one at a time.

  The function returns -1 and sets errno in case of malloc() failure,
  0 otherwise. If 'nxval' is not in a correct range for any class, set
  'sxc' and 'lxc' to NULL. 
*/
int xpart(int c, int *nd, int nxval, int ***sxc, int ***lxc)
{
  int i;
  int j;
  int done;
  int status;
  int ccc;
  int *sx;
  int *lx;
  int **tsxc;
  int **tlxc;

  status = 0;
  done = 0;
  tsxc = malloc(c*sizeof(int *));
  tlxc = malloc(c*sizeof(int *));
  if ((tsxc == (int **) 0) || (tlxc == (int **) 0))
    status = -1;
  else
    for (i = 0; (i < c) && (done == 0); i++)
      {
	ccc = ivec_sum(nd, i);
	status = xss(ccc, ccc+nd[i]-1, nxval, &sx, &lx);
	if (status == 0)
	  {
	    if ((sx != (int *) 0) && (lx != (int *) 0))
	      {
		tsxc[i] = ivec_clone(sx, nd[i]);
		tlxc[i] = ivec_clone(lx, nxval);
		free(sx);
		free(lx);
	      }
	    else
	      {
		for (j = 0; j < i; j++)
		  {
		    free(tsxc[j]);
		    free(tlxc[j]);
		  }
		free(tsxc);
		free(tlxc);
		tsxc = (int **) 0;
		tlxc = (int **) 0;
		done = 1;
	      }
	  }
      }
  *sxc = tsxc;
  *lxc = tlxc;
  return status;
}

/*
  Create learning and validation subsets of `dset' as specified by
  `lxc', `sxc'. See xpart() for details.
*/
int xset(struct dataset *dset, int **lxc, int **sxc, int idx, 
	 struct dataset **learning_dset, struct dataset **validation_dset,
	 FILE *fdbg)
{
  int status;
  int jx;
  int nvec;
  int nval;
  int c;
  int *validation_nd;
  int *learning_nd;
  float **learning_x; /* learning vectors */
  float **validation_x; /* validation vectors */
  struct dataset *vset;
  struct dataset *lset;

  status = 0;
  learning_x = (float **) 0;
  validation_x = (float **) 0;
  vset = (struct dataset *) 0;
  c = dset->c;
  validation_nd = malloc(c*sizeof(int));
  learning_nd = malloc(c*sizeof(int));
  for (jx = 0; jx < c; jx++)
    {
      validation_nd[jx] = lxc[jx][idx];
      learning_nd[jx] = dset->nd[jx]-validation_nd[jx];
    }
  extract_sets(dset->x, idx, c, dset->nd, dset->d, sxc, lxc, 
	       &learning_x, &validation_x, fdbg);
  if (validation_dset)
    {
      nvec = ivec_sum(validation_nd, c);
      vset = dataset_lt(dset->d, c, validation_nd, 
			nvec, (char **) 0, validation_x);
    }
  else
    {
      nval = 0;
      for (jx = 0; jx < c; jx++)
	nval += lxc[jx][idx];
      mx_free((void **) validation_x, nval);
    }
  if (vset || !validation_dset)
    {
      nvec = ivec_sum(learning_nd, c);
      lset = dataset_lt(dset->d, c, learning_nd, nvec,
			(char **) 0, learning_x);
      if (lset)
	{
	  if (validation_dset)
	    *validation_dset = vset;
	  *learning_dset = lset;
	}
      else
	status = -1;
    }
  else
    status = -1;
  return status; 
}

static void preamble(FILE *fptr, int method, int nmodels, int nxval, int nexp, 
		     int normalize, int dr_method, int fscrit, int d, int c)
{
  if (method == PALG_SVM)
    fprintf(fptr, "# classifier:                           SVM\n");
  else if (method == PALG_BAG_SVM)
    fprintf(fptr, "# classifier:                           bagging SVM\n");
  else if (method == PALG_ADABOOST_SVM)
    fprintf(fptr, "# classifier:                           Adaboost SVM\n");
  else if ((method == PALG_KNN) || (method == PALG_BAG_KNN))
    fprintf(fptr, "# classifier:                           k-nearest neighbor\n");
  else if ((method == PALG_MLP) || (method == PALG_BAG_MLP))
    fprintf(fptr, "# classifier:                           multi-layer perceptron\n");
  else if ((method == PALG_LIN) || (method == PALG_BAG_LIN))
    fprintf(fptr, "# classifier:                           linear discriminant\n");
  else if ((method == PALG_PLC) || (method == PALG_BAG_PLC))
    fprintf(fptr, "# classifier:                           parametric linear classifier\n");
  else if ((method == PALG_PQC) || (method == PALG_BAG_PQC))
    fprintf(fptr, "# classifier:                           parametric quadratic classifier\n");
  if (nmodels > 0)
    fprintf(fptr, "# number of models:                     %-d\n", nmodels);
  fprintf(fptr, "# number of experiments:                %-d\n", nexp);
  fprintf(fptr, "# number of cross-validation subsets:   %-d\n", nxval);
  fprintf(fptr, "# number of classes:                    %-d\n", c);
  if (normalize == 1)
    fprintf(fptr, "# normalization:                        yes\n");
  else
    fprintf(fptr, "# normalization:                        no\n");
  if (dr_method == PDR_NONE)
    fprintf(fptr, "# dimensionality reduction method:      %-s\n", PDS_NONE);
  else if (dr_method == PDR_SVD)
    fprintf(fptr, "# dimensionality reduction method:      %-s\n", PDS_SVD);
  else if (dr_method == PDR_FISHER)
    fprintf(fptr, "# dimensionality reduction method:      %-s\n", PDS_FISHER);
  else if (dr_method == PDR_PCA)
    fprintf(fptr, "# dimensionality reduction method:      %-s\n", PDS_PCA);
  else if (dr_method == PDR_RANKING)
    fprintf(fptr, "# dimensionality reduction method:      %-s\n", PDS_RANKING);
  else if (dr_method == PDR_EMAP)
    fprintf(fptr, "# dimensionality reduction method:      %-s\n", PDS_EMAP);
  else if (dr_method == PDR_FORWARD)
    fprintf(fptr, "# dimensionality reduction method:      %-s\n", PDS_FORWARD);
  else if (dr_method == PDR_BACKWARD)
    fprintf(fptr, "# dimensionality reduction method:      %-s\n", PDS_BACKWARD);
  else if (dr_method == PDR_FLOAT_FORWARD)
    fprintf(fptr, "# dimensionality reduction method:      %-s\n", PDS_FLOAT_FORWARD);
  else if (dr_method == PDR_FLOAT_BACKWARD)
    fprintf(fptr, "# dimensionality reduction method:      %-s\n", PDS_FLOAT_BACKWARD);
  if (dr_method != PDR_NONE)
    {
      if (fscrit == PCP_FSEL_EUCLIDEAN)
	fprintf(fptr, "# feature evaluation criterion:         %-s\n", PCP_SCRIT_EUCLIDEAN);
      else if (fscrit == PCP_FSEL_PEARSON)
	fprintf(fptr, "# feature evaluation criterion:         %-s\n", PCP_SCRIT_PEARSON);
      else if (fscrit == PCP_FSEL_GOLUB)
	fprintf(fptr, "# feature evaluation criterion:         %-s\n", PCP_SCRIT_GOLUB);
      else if ((fscrit == PCP_FSEL_KNN) || (fscrit == PCP_FSUB_KNN))
	fprintf(fptr, "# feature evaluation criterion:         %-s\n", PCP_SCRIT_KNN);
      else if ((fscrit == PCP_FSEL_BAYES) || (fscrit == PCP_FSUB_BAYES))
	fprintf(fptr, "# feature evaluation criterion:         %-s\n", PCP_SCRIT_BAYES);
      else if (fscrit == PCP_FSUB_IN_IN)
	fprintf(fptr, "# feature evaluation criterion:         %-s\n", PCP_SCRIT_IN_IN);
    }
  fprintf(fptr, "# number of features:                   %-d\n", d);
}

/*
  Write SVM results file preamble. Return 0 in case of success, errno
  error code in case of failure.
*/
static int svm_save_preamble(int method, int nmodels, int nxval, int nexp, int c, 
			     int normalize, int dr_method, int fscrit, int d, 
			     struct svm_parameter *parameters, unsigned int seed, 
			     char *fname)
{
  int   jx; 
  int   status = 0;
  FILE  *fptr;

  fptr = fopen(fname, "w");
  if (fptr)
    {
      preamble(fptr, method, nmodels, nxval, nexp, normalize, dr_method, fscrit, d, c);
      if ((parameters->nr_weight > 0) && (parameters->svm_type == C_SVC))
	for (jx = 0; jx < c; jx++)
	  fprintf(fptr, "# class %5d cost:                     %-12.5g\n", jx+1, parameters->weight[jx]);
      fprintf(fptr, "# seed:                                 %-d\n", seed);
      if (parameters->svm_type == C_SVC)
	fprintf(fptr, "# SVM type:                             C-SVM\n");
      else
	fprintf(fptr, "# SVM type:                             nu-SVM\n");
      if (parameters->kernel_type == LINEAR)
	fprintf(fptr, "# kernel type:                          linear\n");
      else if (parameters->kernel_type == POLY)
	fprintf(fptr, "# kernel type:                          polynomial\n");
      else if (parameters->kernel_type == RBF)
	fprintf(fptr, "# kernel type:                          RBF\n");
      else if (parameters->kernel_type == SIGMOID)
	fprintf(fptr, "# kernel type:                          sigmoid\n");
      if (parameters->kernel_type == POLY)
	fprintf(fptr, "# kernel degree:                        %-5.2f\n", parameters->degree);
      if ((parameters->kernel_type == POLY) ||
	  (parameters->kernel_type == RBF) ||
	  (parameters->kernel_type == SIGMOID))
	fprintf(fptr, "# gamma:                                %-14.7g\n", parameters->gamma);
      if ((parameters->kernel_type == POLY) ||
	  (parameters->kernel_type == SIGMOID))
	fprintf(fptr, "# coef0:                                %-5.2f\n", parameters->coef0);
      /*fprintf(fptr, "# cache size (MB):                      %-12.2g\n", parameters->cache_size);*/
      fprintf(fptr, "# stopping criterion:                   %-12.5g\n", parameters->eps);
      if (parameters->svm_type == C_SVC)
	fprintf(fptr, "# cost (C):                             %-12.5g\n", parameters->C);
      if (parameters->svm_type == NU_SVC)
	fprintf(fptr, "# nu:                                   %-12.5g\n", parameters->nu);
      status = fclose(fptr);
      if (status != 0)
	status = errno;
    }
  else
    status = errno;
  return status;
}

static void xlearn_free(void **models, int nmodels, int c, int method, void **problems)
{
  int i;

  if (models)
    {
      for (i = 0; i < nmodels; i++)
	{
	  if ((method == PALG_SVM) || (method == PALG_ADABOOST_SVM) || 
	      (method == PALG_BAG_SVM))
	    {
	      svm_destroy_model(models[i]);
	      if (problems)
		free_problem((struct svm_problem *) problems[i]);
	    }
	  else if ((method == PALG_LIN) || (method == PALG_BAG_LIN) ||
		   (method == PALG_PLC) || (method == PALG_BAG_PLC))
	    mx_free(models[i], c);
	  else if ((method == PALG_PQC) || (method == PALG_BAG_PQC))
	    qmodel_free(models[i], c);
	  else if ((method == PALG_MLP) || (method == PALG_BAG_MLP))
	    mlp_free(models[i]);
	}
      free(models);
      free(problems);
    }
}

/*
  Normalize attributes in 'learning_dset' and 'validation_dset'. The
  normalization attributes are computed on 'learning_dset', then
  applied to both datasets.

  In case of success, return 0. Otherwise return -1 and set errno. The
  possible failures are memory allocation errors.
*/
int normalize_attributes(struct dataset *learning_dset, struct dataset *validation_dset)
{
  int   status;
  int   nv;
  int   d;
  float *std;
  float *xmean;
  float **x;

  status = 0;
  if (learning_dset)
    {
      nv = learning_dset->nv;
      d = learning_dset->d;
      x = learning_dset->x;
      xmean = fmx_mean(x, nv, d);
      if (xmean)
	{
	  std = fmx_std(x, nv, d);
	  if (std)
	    {
	      fmx_prenorm(x, nv, d, xmean, std);
	      if (validation_dset)
		fmx_prenorm(validation_dset->x, validation_dset->nv, d, xmean, std);
	      free(std);
	    }
	  else
	    status = -1;
	  vx_free(xmean);
	}
      else
	status = -1;
    }
  return status;
}

/*
  Utility function for reduce_d().
*/
static int clone_sets(struct dataset *learning_dset, struct dataset *validation_dset,
		      struct dataset **l_dset, struct dataset **v_dset, int *errc)
{
  int status;

  status = 0;
  *v_dset = dataset_clone(validation_dset);
  if (*v_dset)
    {
      *l_dset = dataset_clone(learning_dset); 
      if (!*l_dset)
	{
	  status = -1;
	  *errc = errno;
	}
    }
  else
    {
      status = -1;
      *errc = errno;
    }
  return status;
}

/*
  Utility function for reduce_d().
*/
static int map_sets(struct dataset *learning_dset, struct dataset *validation_dset,
		    struct dataset **l_dset, struct dataset **v_dset, 
		    int idr, float **map, int ifd, int *errc)
{
  int status;

  status = 0;
  *l_dset = dataset_map(learning_dset, idr, map); 
  if (!*l_dset)
    {
      status = -1;
      *errc = errno;
    }
  else
    {
      *v_dset = dataset_map(validation_dset, idr, map);
      if (!*v_dset)
	{
	  status = -1;
	  *errc = errno;
	}
    }
  mx_free((void **) map, ifd);
  return status;
}

/*
  Reduce dimensionality of 'learning_dset' and 'validation_dset'.  The
  processed data sets are returned in 'l_dset' and 'v_dset',
  respectively. The dimensionality reduction is performed using
  'dr_method'. The number of attributes retained is 'idr'. 'fscrit' is
  the feature subset evaluation criterion used for feature subset
  selection methods.

  Note that the processing is fully cross-validated: the
  dimensionality reduction parameters are computed using the training
  set, and applied to both sets.

  In case of success, return 0. In case of failure, return -1 and set
  errc.
  
  TBD: where should this function go?
*/
int reduce_d(struct dataset *learning_dset, struct dataset *validation_dset,
	     struct dataset **l_dset, struct dataset **v_dset, 
	     int dr_method, int idr, int fscrit, int *errc)
{
  int   status;
  int   d;
  int   idx;
  int   nv;
  int   *findex;
  float **map;
  float **x;

  findex = (int *) 0;
  map = (float **) 0;
  if (!errc)
    status = -1;
  else if (!learning_dset || !validation_dset || !l_dset || !v_dset)
    {
      status = -1;
      *errc = EINVAL;
    }
  else
    {
      *errc = 0;
      status = 0;
      d = learning_dset->d;
      nv = learning_dset->nv;
      x = learning_dset->x;
      if (dr_method == PDR_NONE)
	status = clone_sets(learning_dset, validation_dset, l_dset, v_dset, errc);
      else 
	{
	  /*
	    Optionally reduce dimension of training and validation data sets.
	  */
	  if ((dr_method == PDR_RANKING) || (dr_method == PDR_FORWARD) ||
	      (dr_method == PDR_BACKWARD) || (dr_method == PDR_PLUS_L_MINUS_R) ||
	      (dr_method == PDR_BB) || (dr_method == PDR_FLOAT_FORWARD))
	    {
	      *l_dset = select_subset(learning_dset, idr, dr_method, fscrit, &findex, 1, errc);
	      if (*l_dset)
		{
		  *v_dset = dataset_select(validation_dset, findex, idr); 
		  if (!*v_dset)
		    {
		      *errc = errno;
		      status = -1;
		    }
		}
	      else
		status = -1;
	      free(findex);
	    }
	  else if (dr_method == PDR_SVD)
	    {
	      /*
		Calculate mapping matrix, then map data sets.
	      */
	      map = svd_transform(x, nv, d, errc);
	      /*map = svd_transform(tds->x, tds->nv, tds->d, errc);*/ /* TESTING - use to check 'map' */
	      if (map)
		status = map_sets(learning_dset, validation_dset, l_dset, v_dset, 
				  idr, map, nv, errc);
	      else
		status = -1;
	    }
	  else if ((dr_method == PDR_PCA) || (dr_method == PDR_FISHER))
	    {
	      if (dr_method == PDR_PCA)
		map = pca(learning_dset, errc);
	      else if (dr_method == PDR_FISHER)
		map = fld(learning_dset, errc);
	      if (map)
		status = map_sets(learning_dset, validation_dset, l_dset, v_dset, 
				  idr, map, idr, errc);
	      else
		status = -1;
	    }
	  else if (dr_method == PDR_EMAP)
	    {
	      /*
		Calculate mapping matrix, then map data sets.
	      */
	      idx = ivec_min(learning_dset->nd, learning_dset->c)-1;
	      map = emap(learning_dset, idr, 1, idx, errc);
	      if (map)
		status = map_sets(learning_dset, validation_dset, l_dset, v_dset, 
				  idr, map, idr, errc);
	      else
		status = -1;
	    }
	  else
	    {
	      /*
		Unrecognized mapping method. 
	      */
	      status = -1;
	      *errc = EINVAL;
	    }
	}
    }
  return status;
}

/*
  Write intermediate cumulative cross-validation results in
  'fname'. 

  The function saves summary results for subsets 1..ixval. Thus, for
  each cross-validation experiment, the function gets called ixval
  times: first to save intermediate cumulative results for subset 1,
  then for subsets 1..2, etc. until final invocation to save
  cumulative results for subsets 1..nxval.

  'ec' is class-conditional error count. In other words, ec[i] is the
  cumulative number of misclassified class i samples in
  cross-validation subsets 1..ixval. 'lxc[i]' is ixval-long vector of
  cardinalities of class i test subsets 1..ixval.

  Return 0 for success. In case of error, return -1 and set errno.
*/
static int xlearn_save(int ex, int ixval, int d, unsigned int seed, float amce, 
		       int c, int *ec, int indet, int **lxc, char *fname)
{
  int   icc;
  int   jx; 
  int   icntr;
  int   nl;
  int   llen;
  int   status;
  int   fd;
  int   errno_save;
  int   n_P;
  int   n_N;
  int   F_P;
  int   T_P;
  int   F_N;
  int   T_N;
  float e_P;
  float e_N;
  float xmce;
  float sensitivity;
  float specificity;
  float ppv;
  float npv;
  char  *template;
  char  *line;
  char  **tokens;
  FILE  *fptr;
  FILE  *tptr;

  status = 0;
  n_P = -1;
  e_P = -1.0;
  nl = file_info(fname, &llen, (int *) 0, '\0');
  if (nl >= 0)
    { 
      line = malloc((llen+2)*sizeof(char));
      icntr = 0;
      fptr = fopen(fname, "r");
      if (fptr)
	{
	  template = malloc(20*sizeof(char));
	  sprintf(template, "XXXXXX");
	  fd = mkstemp(template);
	  if (fd != -1)
	    {
	      status = fchmod(fd, S_IRUSR | S_IWUSR | S_IRGRP | S_IROTH);
	      if (status != -1)
		{
		  tptr = fdopen(fd, "w");
		  if (tptr)
		    {
		      while (fgets(line, llen+2, fptr))
			{
			  if (*line == '#')
			    fprintf(tptr, "%s", line);
			  else
			    {
			      tokens = str_tokenize(line, "\t");
			      if (atoi(tokens[0]) < ex)
				fprintf(tptr, "%s", line);
			      str_free(tokens);
			    }
			}
		      free(line);
		      fclose(fptr);
		      icc = 0;
		      for (jx = 0; jx < c; jx++)
			icc += ivec_sum(lxc[jx], ixval); /* cumulative cardinality of test subsets 1..ixval */
		      xmce = (100.0*amce)/icc; /* cumulative average mce */
		      fprintf(tptr, "%d\t%5.2f\t(%d/%d)\t", ex, xmce, (int) amce, icc);
		      /*
			Class-conditional error rate for class i is
			calculated as number of class i samples
			assigned to other classes, divided by
			cardinality of class i, expressed in
			percentages. Cumulative error rate is the
			number of misclassified samples divided by the
			total number of samples, expressed in
			percentages.
			
			For two classes, additional accuracy measures
			(sensitivity, specificity, positive and
			negative predictive value) are computed. These
			measures are frequently used in clinical
			studies and assume that one of the classes is
			labeled 'positive' (corresponding to patients
			suffering from a disease) and the other is
			labeled 'negative' (healthy patients). In PCP,
			class 1 is by convention considered
			'positive', and class 2 'negative'. All four
			measures are expressed as percentages.
			
			See PCP User Guide for explanation of
			sensitivity, specificity, ppv and npv
			calculations.
		      */
		      for (jx = 0; jx < c; jx++)
			{
			  icc = ivec_sum(lxc[jx], ixval); /* cumulative card. of test subsets 1..ixval for class jx */
			  fprintf(tptr, "%5.2f\t(%d/%d)\t", (100.0*ec[jx])/icc, ec[jx], icc);
			  if (c == 2)
			    {
			      if (jx == 0)
				{
				  n_P = icc;
				  e_P = ((float) ec[jx])/icc;
				}
			      else
				{
				  n_N = icc;
				  e_N = ((float) ec[jx])/icc;
				  F_P = ec[1];
				  T_P = n_P-ec[0];
				  F_N = ec[0];
				  T_N = n_N-ec[1];
				  sensitivity = 100.0*(1.0-e_P);
				  specificity = 100.0*(1.0-e_N);
				  if (T_P+F_P > 0)
				    ppv = 100.0*T_P/(T_P+F_P);
				  else
				    /*
				      If no samples were assigned to P
				      class (that is, T_P+F_P is 0),
				      ppv is undefined. We set it to
				      0.
				    */
				    ppv = 0.0;
				  if (T_N+F_N > 0)
				    npv = 100.0*T_N/(T_N+F_N);
				  else
				    npv = 0.0;
				  fprintf(tptr, "%5.2f\t%5.2f\t%5.2f\t%5.2f\t", sensitivity, specificity, ppv, npv);
				}
			    }
			}
		      fprintf(tptr, "\n");
		      status = fflush(tptr);
		      if (status == 0)
			{
			  status = close(fd);
			  if (!status)
			    {
			      status = rename(template, fname);
			      if (!status)
				free(template);
			    }
			}
		      else
			{
			  errno_save = errno;
			  close(fd);
			  free(template);
			  errno = errno_save;
			}
		    }
		  else
		    status = -1;
		}
	    }
	  else
	    status = -1;
	}
      else
	status = -1;
    }
  else
    status = -1;
  return status;
}

/*
  Predict class for `svm_vector' using `binary_models' for `c'
  classes.

  The models are of the type `class i vs. the rest'. For example,
  model 1 is `class 1 vs. the rest'; model 2 is `class 2 vs. the
  rest', etc. `svm_vector' is assigned to class i if binary model i
  predicts class i, and all others predict `rest'. Otherwise, the
  prediction is considered ambiguous.

  Consistent with the PCP standard, the function returns a values in
  the range [0..c-1] in case it is able to make a prediction, and
  value PCP_UNASSIGNED if the prediction is ambiguous.
*/
static int binary_predict(struct svm_model **binary_models, 
			  struct svm_node *svm_vector, int c)
{
  int i;
  int done;
  int icl;
  int prediction;

  done = 0;
  prediction = PCP_UNASSIGNED; 
  for (i = 0; i < c && !done; i++)
    {
      icl = svm_predict(binary_models[i], svm_vector)-1;
      if (icl == 0)
	{
	  /*
	    SVM voted for class i. Check if we already have a vote for
	    another class.
	  */
	  if (prediction != PCP_UNASSIGNED)
	    {
	      /*
		We already have a vote - return ambiguous class
		prediction.
	      */
	      done = 1;
	      prediction = PCP_UNASSIGNED;
	    }
	  else
	    prediction = i;
	}
    }
  return prediction;
}

/*
  Train classifier `method' on `learning_dset', test it on
  `validation_dset', and report cumulative and class-conditional error
  counts. Algorithm-specific parameters of `method' are stored in
  `options'.

  Optionally reduce dimension of the training/validation data sets to
  `idr' using `dr_method'. If `dr_method' is a feature selection
  method, use `fscrit' feature subset evaluation criterion.

  If `normalize' is 1, normalize the data sets before applying the
  learning method and dimension reduction. NOTE: this will change
  'learning_dset' and 'validation_dset'. Otherwise, the sets are
  unchanged on return from xlearn().

  `nmodels' is the number of models for bagging/Adaboost methods.

  The function returns the number of misclassified samples, and sets
  class-conditional error counts in `ccer'. In case of failure,
  returns -1 and sets `errc'.

  TBD: throwaway `seed'.
*/
int xlearn(int method, int dr_method, int idr, int fscrit, int normalize, 
	   struct dataset *learning_dset, struct dataset *validation_dset,
	   unsigned int seed, int nmodels, int bag_size, void *options, 
	   int **ccer, FILE *outdev, FILE *fdbg, int *errc)
{
  int   status;
  int   i;
  int   j;
  int   c;
  int   jx;
  int   icl;
  int   tcl;
  int   tcum;
  int   mce;
  int   imdl;
  int   indet;
  char  *temp_fname;
  float *output;
  float *combined_output = (float *) 0;
  float *weights = (float *) 0;
  int   *ccr = (int *) 0;
  void  **models = (void **) 0;
  float **target;
  struct mlp *perceptron;
  struct mlp **perceptrons;
  struct dataset *l_dset; /* preprocessed learning_dset */
  struct dataset *v_dset; /* preprocessed validation_dset */
  struct svm_node *svm_vector;
  struct mlp_options *mlp_optional;
  struct knn_options *knn_options;
  struct svm_problem *problem;
  void   **problems = (void **) 0;
  struct svm_model *model;
  struct svm_model **binary_models;
  char  *func = "xlearn()";

  binary_models = (struct svm_model **) 0;
  indet = 0;
  if ((method == PALG_SVM) && ((struct svm_parameter *) options)->indet == 1)
    indet = 1;
  perceptrons = (struct mlp **) 0;
  svm_vector = (struct svm_node *) 0;
  icl = -1;
  *errc = 0;
  mce = -1;
  c = learning_dset->c;
  temp_fname = malloc((L_tmpnam+1)*sizeof(char));
  tmpnam(temp_fname);
  status = 0;
  if (normalize)
    status = normalize_attributes(learning_dset, validation_dset);
  if (!status)
    {
      status = reduce_d(learning_dset, validation_dset, &l_dset, &v_dset,
			dr_method, idr, fscrit, errc);
      if (!status)
	{
	  if (outdev)
	    inverse_video();
	  /*
	    Build models (except for k-NN methods).
	  */
	  if (method == PALG_MLP)
	    {
	      target = mlp_target(c, l_dset->nd);
	      mlp_optional = (struct mlp_options *) options;
	      perceptron = mlp_learn(mlp_optional->opt_method, l_dset->x, l_dset->nv,
				     l_dset->nd, idr, target, mlp_optional->nlayers,
				     mlp_optional->npl, mlp_optional->itmax, 
				     mlp_optional->range, mlp_optional->eta, 
				     mlp_optional->mu, outdev, 0, (char *) 0, seed,
				     errc, fdbg);
	      mx_free((void **) target, l_dset->nv);
	      if (perceptron)
		{
		  nmodels = 1;
		  perceptrons = calloc(2, sizeof((struct mlp **) 0));
		  perceptrons[0] = perceptron;
		  models = (void **) perceptrons;
		}
	      else
		status = -1;
	    }
	  else if (method == PALG_BAG_MLP)
	    {
	      mlp_optional = (struct mlp_options *) options;
	      models = bagging(l_dset, outdev, method, temp_fname, seed, nmodels,
			       bag_size, options, (void ***) 0, errc, fdbg);
	      if (models)
		perceptrons = (struct mlp **) models;
	      else
		status = -1;
	    }
	  else if ((method == PALG_ADABOOST_MLP) || (method == PALG_ADABOOST_SVM))
	    {
	      models = adaboost(l_dset, method, &nmodels, &weights, temp_fname, seed,
				options, errc, fdbg);
	      if (!models)
		status = -1;
	    }
	  else if (method == PALG_SVM)
	    {
	      if (indet)
		{
		  /*
		    Create the special structure for the case of
		    multi-class SVM using indeterminate region.
		  */
		  binary_models = malloc(l_dset->c*sizeof(struct svm_model *));
		  for (i = 0; i < l_dset->c; i++)
		    {
		      problem = create_binary_problem(l_dset, i);
		      binary_models[i] = svm_train(problem, (struct svm_parameter *) options);
		    }
		}
	      else
		{
		  problem = create_problem(l_dset);
		  problems = malloc(sizeof(struct svm_problem *));
		  problems[0] = problem;
		  model = svm_train(problem, (struct svm_parameter *) options);
		  if (model)
		    {
		      models = malloc(sizeof(void **));
		      models[0] = model;
		      weights = malloc((c+1)*sizeof(float));
		      fvec_set(weights, c+1, 1.0);
		    }
		  else
		    status = -1;
		}
	    }
	  else if (method == PALG_BAG_SVM)
	    {
	      models = bagging(l_dset, outdev, method, temp_fname, seed, nmodels,
			       bag_size, options, &problems, errc, fdbg);
	      if (models)
		{
		  weights = malloc((c+1)*sizeof(float));
		  fvec_set(weights, c+1, 1.0);
		}
	      else
		status = -1;
	    }
	  else if ((method == PALG_LIN) || (method == PALG_BAG_LIN) || 
		   (method == PALG_PLC) || (method == PALG_BAG_PLC))
	    {
	      if ((method == PALG_BAG_LIN) || (method == PALG_BAG_PLC))
		models = bagging(l_dset, outdev, method, temp_fname, seed, nmodels,
				 bag_size, options, &problems, errc, fdbg);
	      else if (method == PALG_LIN)
		{
		  models = malloc(sizeof(void **));
		  models[0] = (void **) lind_learn(l_dset, errc, fdbg);
		}
	      else if (method == PALG_PLC)
		{
		  models = malloc(sizeof(void **));
		  models[0] = (void **) lin_learn(WEIGHTED_COV, l_dset, errc);
		}
	      if (models && models[0])
		{
		  weights = malloc((c+1)*sizeof(float));
		  fvec_set(weights, c+1, 1.0);
		}
	      else
		status = -1;
	    }
	  else if ((method == PALG_PQC) || (method == PALG_BAG_PQC))
	    {
	      if (method == PALG_BAG_PQC)
		models = bagging(l_dset, outdev, method, temp_fname, seed, nmodels,
				 bag_size, options, &problems, errc, fdbg);
	      else
		{
		  models = malloc(sizeof(void **));
		  models[0] = (void **) pqc_learn(l_dset, errc);
		}
	      if (models && models[0])
		{
		  weights = malloc((c+1)*sizeof(float));
		  fvec_set(weights, c+1, 1.0);
		}
	      else
		status = -1;
	    }
	  /*
	    This fills the disk quite quickly in the debug (i.e.,
	    regression) mode. Have to remove the temporary file for
	    now. If needed for special debugging, change this
	    statement to preserve the temporary file. ljb, 07/13/2004
	  */
	  /*if (!fdbg)*/
	  unlink(temp_fname);
	  free(temp_fname);
	  if ((method != PALG_KNN) && (method != PALG_BAG_KNN))
	    l_dset = dataset_free(l_dset);
	  /*
	    Classify using the built models.
	  */
	  if (!status && (models || (method == PALG_KNN) || (method == PALG_BAG_KNN)))
	    {
	      ccr = calloc(c+1, sizeof(float));
	      mce = 0;
	      combined_output = malloc(c*sizeof(float));
	      tcl = 0; /* correct class of vector j; icl is predicted class of vector j */
	      tcum = v_dset->nd[0]; /* cumulative cardinality of classes 0..tcl */
	      for (j = 0; (j < v_dset->nv) && !status; j++)
		{
		  /*
		    For each validation set vector, generate
		    predictions for all models. combined_output[i] has
		    weighted number of votes for class i.
		  */
		  if (j == tcum)
		    {
		      tcl++;
		      tcum += v_dset->nd[tcl];
		    }
		  fvec_set(combined_output, c, 0.0);
		  if ((method == PALG_MLP) || (method == PALG_BAG_MLP))
		    {
		      for (imdl = 0; imdl < nmodels; imdl++)
			{
			  output = mlp_output(perceptrons[imdl], v_dset->x[j]);
			  for (jx = 0; jx < c; jx++)
			    combined_output[jx] += output[jx];
			  free(output);
			}
		      for (jx = 0; jx < c; jx++)
			combined_output[jx] /= nmodels;
		    }
		  else if ((method == PALG_ADABOOST_MLP) || (method == PALG_ADABOOST_SVM) || 
			   (method == PALG_SVM) || (method == PALG_BAG_SVM))
		    {
		      if ((method == PALG_ADABOOST_SVM) || (method == PALG_SVM) || 
			  (method == PALG_BAG_SVM))
			svm_vector = create_svm_vector(v_dset->x[j], idr);
		      if (indet)
			{
			  icl = binary_predict(binary_models, svm_vector, v_dset->c);
			  combined_output[icl] += weights[icl];
			}
		      else
			{
			  for (imdl = 0; imdl < nmodels; imdl++)
			    {
			      if ((method == PALG_ADABOOST_SVM) || (method == PALG_SVM) || 
				  (method == PALG_BAG_SVM))
				icl = svm_predict(models[imdl], svm_vector)-1;
			      else if (method == PALG_ADABOOST_MLP)
				{
				  ; /* TBD: implement MLP Adaboost cross-validation */
				}
			      combined_output[icl] += weights[icl];
			    }
			}
		      if ((method == PALG_ADABOOST_SVM) || (method == PALG_SVM) || 
			  (method == PALG_BAG_SVM))
			free(svm_vector);
		    }
		  else if ((method == PALG_LIN) || (method == PALG_BAG_LIN) ||
			   (method == PALG_PLC) || (method == PALG_BAG_PLC) ||
			   (method == PALG_PQC) || (method == PALG_BAG_PQC))
		    {
		      for (imdl = 0; imdl < nmodels; imdl++)
			{
			  if ((method == PALG_PQC) || (method == PALG_BAG_PQC))
			    output = pqc_output(models[imdl], c, idr+1, v_dset->x[j]);
			  else
			    output = lin_output(models[imdl], c, idr+1, v_dset->x[j]);
			  for (jx = 0; jx < c; jx++)
			    combined_output[jx] += output[jx];
			  free(output);
			}
		      for (jx = 0; jx < c; jx++)
			combined_output[jx] /= nmodels;
		    }
		  else if (method == PALG_KNN)
		    {
		      /* classify v_dset->x[j]) */
		      knn_options = (struct knn_options *) options;
		      icl = knn(v_dset->x[j], tcl, l_dset, knn_options->k,
				knn_options->dist, errc, fdbg);
		    }

		  else if (method == PALG_BAG_KNN)
		    {
		      /* classify v_dset->x[j]) using k-NN bagging */
		      knn_options = (struct knn_options *) options;
		      icl = knn_bagging(v_dset->x[j], tcl, l_dset, nmodels, bag_size,
					knn_options->k, knn_options->dist, errc, fdbg);
		      if (*errc != 0)
			status = -1;
		    }
		  if ((method != PALG_KNN) && (method != PALG_BAG_KNN))
		    icl = fvec_argmax(combined_output, c); /* assigned class for non-k-NN methods */
		  if (fdbg)
		    {
		      fprintf(fdbg, "%s; vector %d; assigned class: %d; class: %d; ",
			      func, j, icl, tcl);
		      for (jx = 0; jx < c; jx++)
			fprintf(fdbg, "combined_output[%d]: %12.5g; ", jx, combined_output[jx]);
		      fprintf(fdbg, "\n");
		    }
		  if (!status)
		    {
		      if (icl != tcl)
			{
			  if (icl != PCP_UNASSIGNED)
			    {
			      mce++;
			      ccr[tcl]++;
			    }
			  else 
			    ccr[c]++; /* ambiguous prediction */
			}
		    }
		  else
		    mce = -1;
		}
	      free(combined_output);
	      free(weights);
	      dataset_free(l_dset);
	      dataset_free(v_dset);
	    }
	  xlearn_free(models, nmodels, c, method, problems);
	  remove_file(PCP_ADA);
	  if (ccer)
	    *ccer = ccr;
	  else
	    free(ccr);
	  reset_video();
	}
    }
  return mce;
}

/*
  Write cross-validation results file preamble. Return 0 in case of
  success, errno error code in case of failure.
*/
static int save_preamble(int method, int nmodels, int nxval, int nexp, int c,
			 int normalize, int dr_method, int fscrit, int d, void *options,
			 int itmax, int nlayers, int *npl, unsigned int seed, char *fname)
{
  int   jx; 
  int   status = 0;
  int   dist;
  FILE  *fptr;
  struct knn_options *knn_optional;

  fptr = fopen(fname, "w");
  if (fptr)
    {
      preamble(fptr, method, nmodels, nxval, nexp, normalize, dr_method, fscrit, d, c);
      fprintf(fptr, "# seed:                                 %-d\n", seed);
      if ((method == PALG_MLP) || (method == PALG_BAG_MLP) || (method == PALG_ADABOOST_MLP))
	{
	  fprintf(fptr, "# number of iterations:                 %-d\n", itmax);
	  fprintf(fptr, "# number of inputs:                     %-d\n", d);
	  fprintf(fptr, "# number of layers:                     %-d\n", nlayers);
	  for (jx = 0; jx < nlayers; jx++)
	    fprintf(fptr, "# number of nodes in layer %5d:       %-d\n", jx, npl[jx]);
	}
      else if ((method == PALG_KNN) || (method == PALG_BAG_KNN))
	{
	  knn_optional = (struct knn_options *) options;
	  fprintf(fptr, "# number of neighbors:                  %-d\n",
		  knn_optional->k);
	  dist = knn_optional->dist;
	  if (dist == EUCLIDEAN_DIST)
	    fprintf(fptr, "# distance measure:                     %-d (Euclidean)\n",
		    dist);
	  else if (dist == CITY_BLOCK_DIST)
	    fprintf(fptr, "# distance measure:                     %-d (city-block)\n",
		    dist);
	  else if (dist == MAHALANOBIS_DIST)
	    fprintf(fptr, "# distance measure:                     %-d (Mahalanobis)\n",
		    dist);
	}
      status = fclose(fptr);
      if (status != 0)
	status = errno;
    }
  else
    status = errno;
  return status;
}

/*
  Calculate mean classification error for 'experiments'. 'emn' and
  'imn' are misclassification count and total count integer totals,
  currently unused.
*/
static float get_mean(char **experiments, int offset, int nexp, int *emn, int *imn)
{
  int   idx;
  int   iexp;
  int   ntok;
  float xmn;
  float ft;
  char  *str;
  char  *ptr;
  char  **tokens;

  if (emn != (int *) 0)
    {
      *emn = 0;
      *imn = 0;
    }
  iexp = 0;
  xmn = 0.0;
  for (idx = 0; idx < nexp; idx++)
    {
      tokens = str_tokenize(experiments[idx], WHITESPACE);
      ntok = str_length(tokens);
      if (tokens && (ntok > offset) && tokens[offset])
	{
	  ft = atof(tokens[offset]);
	  iexp++;
	  xmn += ft;
	  if (emn != (int *) 0)
	    {
	      str = tokens[offset+1];
	      str++;
	      ptr = strchr(str, '/');
	      *ptr = '\0';
	      *emn += atoi(str);
	      ptr++;
	      str = strchr(ptr, ')');
	      *str = '\0';
	      *imn += atoi(ptr);
	    }
	}
      str_free(tokens);
    }
  if (iexp)
    xmn = xmn/iexp;
  else
    xmn = xmn/nexp;
  return xmn;
}

static float get_token(char **tokens, int offset)
{
  int   ntok;
  float ft;

  ft = 0.0;
  ntok = str_length(tokens);
  if (tokens && (ntok > offset) && tokens[offset])
    ft = atof(tokens[offset]);
  return ft;
}

static void disp_param(char *line)
{
  char *ptr;
  
  ptr = line+2;
  *ptr = toupper(*ptr);
  printf("| ");
  print_ln(stdout, ptr, PCP_QLEN-3);
  printf("|\n");
}

void xlearn_disp(int *errc, char **xname, char *default_name)
{
  int   i;
  int   ncl;
  int   nf;
  int   nl;
  int   llen;
  int   idx;
  int   iex;
  int   nexp;
  int   offset;
  int   iexp;
  int   verbose;
  float ft;
  float xmn;
  float xstd;
  float sensitivity;
  float specificity;
  float ppv;
  float npv;
  float mean_sensitivity;
  float mean_specificity;
  float mean_ppv;
  float mean_npv;
  float std_sensitivity;
  float std_specificity;
  float std_ppv;
  float std_npv;
  char  *fname;
  char  *line;
  char  *name;
  char  **tokens;
  char  **experiments;
  FILE  *fptr;
  FILE  *fs;

  ncl = 0;
  nf = 0;
  nexp = 0;
  name = (char *) 0;
  experiments = (char **) 0;
  clear_screen();
  cursor_on();
  fname = input_filename(PCP_UMSG_XFNAME, default_name, stdout);
  nl = file_info(fname, &llen, (int *) 0, '\0');
  if (nl >= 0)
    { 
      /*
	The summary results are automatically saved in file
	PCP_XSM. TBD: for now, all summary file errors are quietly
	ignored.
      */
      fs = fopen(PCP_XSM, "w");
      idx = 0;
      iex = 1;
      verbose = input_integer(stdin, stdout, OUTPUT_MSG, PCP_QLEN, &idx, &idx, &iex);
      line = malloc((llen+2)*sizeof(char));
      fptr = fopen(fname, "r");
      idx = 0;
      iex = 0;
      inverse_video();
      printf(PCP_XLINE);
      printf(PCP_EMPTY_LINE);
      printf("|              C R O S S - V A L I D A T I O N   R E S U L T S               |\n");
      printf(PCP_EMPTY_LINE);
      printf(PCP_XLINE);
      while (fgets(line, llen+2, fptr))
	{
	  if (*line == '#')
	    {
	      if (strstr(line, "number of experiments"))
		{
		  nexp = atoi(&line[PCP_HFT]);
		  experiments = calloc(nexp, sizeof(char *));
		  disp_param(line);
		}
	      else if (strstr(line, "number of cross-validation subsets"))
		disp_param(line);
	      else if (strstr(line, "number of classes"))
		ncl = atoi(&line[PCP_HFT]);
	      else if (strstr(line, "dimensionality reduction"))
		disp_param(line);
	      else if (strstr(line, "feature evaluation"))
		disp_param(line);
	      else if (strstr(line, "number of features"))
		nf = atoi(&line[PCP_HFT]);
	      else if (strstr(line, "number of neighbors"))
		disp_param(line);
	      else if (strstr(line, "distance measure"))
		disp_param(line);
	      else if (strstr(line, "seed"))
		disp_param(line);
	      else if (strstr(line, "classifier"))
		disp_param(line);
	      else if (strstr(line, "SVM type"))
		disp_param(line);
	      else if (strstr(line, "kernel type"))
		disp_param(line);
	      else if (strstr(line, "kernel degree"))
		disp_param(line);
	      else if (strstr(line, "class") && strstr(line, "cost"))
		disp_param(line);
	      else if (strstr(line, "cost"))
		disp_param(line);
	      else if (strstr(line, "nu"))
		disp_param(line);
	      else if (strstr(line, "gamma"))
		disp_param(line);
	      else if (strstr(line, "normalization"))
		disp_param(line);
	    }
	  /*
	    Ignore empty lines.
	  */
	  else if (strcspn(line, WHITESPACE) > 0)
	    experiments[iex++] = strdup(line);
	}
      free(line);
      fclose(fptr);
      xmn = get_mean(experiments, 1, nexp, (int *) 0, (int *) 0);
      mean_sensitivity = get_mean(experiments, 7, nexp, (int *) 0, (int *) 0);
      mean_specificity = get_mean(experiments, 8, nexp, (int *) 0, (int *) 0);
      mean_ppv = get_mean(experiments, 9, nexp, (int *) 0, (int *) 0);
      mean_npv = get_mean(experiments, 10, nexp, (int *) 0, (int *) 0);
      xstd = 0.0;
      std_sensitivity = 0.0;
      std_specificity = 0.0;
      std_ppv = 0.0;
      std_npv = 0.0;
      iexp = 0;
      for (idx = 0; idx < nexp; idx++)
	{
	  if (idx == 0)
	    {
	      printf(PCP_XLINE);
	      if (ncl == 2)
		{
		  printf("| Experiment | Sensitivity | Specificity | Pos. Pred. Value | Neg. Pred Value|\n");
		  if (verbose)
		    printf(PCP_XLINE);
		}
	    }
	  tokens = str_tokenize(experiments[idx], WHITESPACE);
	  if (tokens)
	    {
	      ft = get_token(tokens, 1);
	      sensitivity = get_token(tokens, 7);
	      specificity = get_token(tokens, 8);
	      ppv = get_token(tokens, 9);
	      npv = get_token(tokens, 10);
	      iexp++;
	      if (verbose)
		{
		  if (ncl == 2)
		    printf("|  %5d     |   %7.2f   |   %7.2f   |     %7.2f      |    %7.2f     |\n",
			   idx+1, sensitivity, specificity, ppv, npv);
		}
	      xstd += (ft-xmn)*(ft-xmn);
	      std_sensitivity += (sensitivity-mean_sensitivity)*(sensitivity-mean_sensitivity);
	      std_specificity += (specificity-mean_specificity)*(specificity-mean_specificity);
	      std_ppv += (ppv-mean_ppv)*(ppv-mean_ppv);
	      std_npv += (npv-mean_npv)*(npv-mean_npv);
	    }
	  str_free(tokens);
	}
      if (iexp > 1)
	{
	  xstd = sqrt(xstd/(iexp-1));
	  std_sensitivity = sqrt(std_sensitivity/(iexp-1));
	  std_specificity = sqrt(std_specificity/(iexp-1));
	  std_ppv = sqrt(std_ppv/(iexp-1));
	  std_npv = sqrt(std_npv/(iexp-1));
	}
      printf(PCP_XLINE);
      if (ncl == 2)
	{	
	  printf("|    Mean    |    %6.2f   |    %6.2f   |      %6.2f      |     %6.2f     |\n",
		 mean_sensitivity, mean_specificity, mean_ppv, mean_npv);
	  printf("| Std.  dev. |    %6.2f   |    %6.2f   |      %6.2f      |     %6.2f     |\n",
		 std_sensitivity, std_specificity, std_ppv, std_npv);
	  if (fs)
	    fprintf(fs, "%6.2f\t%6.2f\t%6.2f\t%6.2f\t%6.2f\t%6.2f\t%6.2f\t%6.2f\t%6.2f\t%6.2f\t\n",
		    xmn, xstd, mean_sensitivity, std_sensitivity, mean_specificity,
		    std_specificity, mean_ppv, std_ppv, mean_npv, std_npv);
	  printf(PCP_XLINE);
	}
      if (verbose)
	printf("| Experiment |                     Cumulative error rate                     |\n");
      else
	printf("|            |                     Cumulative error rate                     |\n");
      printf(PCP_XLINE);
      if (verbose)
	{
	  for (idx = 0; idx < nexp; idx++)
	    {
	      tokens = str_tokenize(experiments[idx], WHITESPACE);
	      if (tokens)
		{
		  ft = get_token(tokens, 1);
		  printf("|  %5d     |                       %6.2f %10s                       |\n",
			 idx+1, ft, tokens[2]);
		  str_free(tokens);
		}
	    }
	  printf(PCP_XLINE);
	}
      printf("|    Mean    |                       %6.2f                                  |\n",
	     xmn);
      printf("| Std.  dev. |                       %6.2f                                  |\n",
	     xstd);
      if (fs)
	fprintf(fs, "%6.2f\t%6.2f\n", xmn, xstd);
      offset = 3;
      for (i = 0; i < ncl; i++)
	{
	  xmn = get_mean(experiments, offset, nexp, (int *) 0, (int *) 0);
	  xstd = 0.0;
	  iexp = 0;
	  for (idx = 0; idx < nexp; idx++)
	    {
	      iexp++;
	      if (idx == 0)
		{
		  printf(PCP_XLINE);
		  if (tds && (tds->c == ncl) && (tds->d == nf))
		    name = bname(tds->fnames[i]);
		  if (verbose)
		    {
		      if (name)
			printf("| Experiment |        Error rate for class %-30s    |\n", name);
		      else
			printf("| Experiment |        Error rate for class %5d                             |\n",
			       i+1);
		      printf(PCP_XLINE);
		    }
		  else
		    {
		      if (name)
			printf("|            |        Error rate for class %-30s    |\n", name);
		      else
			printf("|            |        Error rate for class %5d                             |\n", 
			       i+1);
		    }
		  free(name);
		}
	      tokens = str_tokenize(experiments[idx], WHITESPACE);
	      if (tokens)
		{
		  ft = atof(tokens[offset]);
		  if (verbose)
		    printf("|  %5d     |                       %6.2f %10s                       |\n",
			   idx+1, ft, tokens[offset+1]);
		  xstd += (ft-xmn)*(ft-xmn);
		}
	      str_free(tokens);
	    }
	  offset += 2;
	  if (iexp > 1)
	    xstd = sqrt(xstd/(iexp-1));
	  printf(PCP_XLINE);
	  printf("|    Mean    |                       %6.2f                                  |\n",
		 xmn);
	  printf("| Std.  dev. |                       %6.2f                                  |\n",
		 xstd);
	  if (fs)
	    fprintf(fs, "%6.2f\t%6.2f\n", xmn, xstd);
	}
      if (fs)
	fclose(fs);
      for (i = 0; i < nexp; i++)
	free(experiments[i]);
      free(experiments);
      printf(PCP_XLINE);
      cursor_off();
      reset_video();
      pwait();
    }
  else
    {
      *errc = errno;
      *xname = strdup(fname);
    }
}

/*
  Return average error rate (in percentages) of `nexp' `nxval'-fold
  cross-validation experiments for classifier `method', using `dset'.

  In case of error, return -1.0 and set `errc'.

  TBD: currently only supports PALG_KNN method.
*/
float xlearn_machine(struct dataset *dset, int method, int nexp, int nxval, int *errc)
{
  int   status;
  int   ex;
  int   idx;
  int   c;
  int   error_count;
  float ecrit;
  float amce;
  int   **sxc;
  int   **lxc;
  void  *learning_params;
  struct dataset *learning_dset;
  struct dataset *validation_dset;

  status = 0;
  *errc = 0;
  learning_params = (void *) 0;
  c = dset->c;
  if (method == PALG_KNN)
    {
      learning_params = malloc(sizeof(struct knn_options));
      ((struct knn_options *) learning_params)->dist = EUCLIDEAN_DIST;
      ((struct knn_options *) learning_params)->k = 1;
    }
  /*
    Compute cross-validation error rates for the given learning
    method and parameters.
  */
  ecrit = 0.0;
  for (ex = 0; (ex < nexp) && !status; ex++)
    {
      amce = 0.0;
      status = xpart(c, dset->nd, nxval, &sxc, &lxc);
      /*
	lxc[i][j] is number of class i vectors in cross-validation
	subset j.
	
	sxc[i] is list of class i vectors with arranged so that
	cross-validation subsets are in sequence.  For example, if
	there are 3 cross-validation subsets with 10 vectors each,
	sxc[i][0..9] is list of class i vectors in subset 0,
	sxc[i][10..19] is list of class i vectors in subset 1, and
	sxc[i][20..29] is list of class i vectors in subset 2.
      */
      if (!status)
	{
	  for (idx = 0; (idx < nxval) && !status; idx++)
	    {
	      status = xset(dset, lxc, sxc, idx, &learning_dset, &validation_dset, (FILE *) 0);
	      error_count = xlearn(method, PDR_NONE, 0, 0, 0, learning_dset, validation_dset, 
				   1, 1, learning_dset->nv, (void *) learning_params, (int **) 0, 
				   (FILE *) 0, (FILE *) 0, errc);
	      if (error_count >= 0)
		amce += error_count;
	      else
		status = -1;
	      dataset_free(validation_dset);
	      dataset_free(learning_dset);
	    }
	  mx_free((void **) sxc, c);
	  mx_free((void **) lxc, c);
	  if (!status)
	    ecrit += 100.0*amce/dset->nv; /* error rate for this experiment */
	}
      else
	*errc = errno;
    }
  if (!status)
    ecrit = ecrit/nexp;
  else
    ecrit = -1.0;
  free(learning_params);
  return ecrit;
}

/*
  Perform multiple cross-validation learning experiments.
  
  Each experiment consists in the following cross-validation test:
  train the classifier using 'method' on a pseudo-randomly chosen
  subset of the training data set (the learning data set), and test it
  on the remainder of the data set (the validation data set); then
  choose another learning/validation partition, repeat the test,
  etc. This is repeated 'nxval' times, where 'nxval' is the number of
  cross-validation subsets.
*/
void p_xlearn(int method, int *errc, char **xname, int dbg)
{
  int   idx;
  int   jx;
  int   c;
  int   nvec; /* cardinality of validation set */
  int   mce;
  int   status;
  int   nexp;
  int   ex;
  int   nxval;
  int   nmodels = -1;
  int   bag_size;
  int   dr_method; /* dimension reduction method code */
  int   idr; /* reduced dimension */
  int   fscrit; /* feature subset selection criterion */
  int   normalize;
  int   dflt;
  int   opt_method;
  int   dist;
  int   indet;
  int   mtd;
  float eta;
  float mu;
  float fmce;
  float amce;
  float xmce;
  char  *msg;
  char  *xval_fname;
  int   *learning_nd; /* cardinality of learning classes */
  int   *validation_nd; /* cardinality of validation classes */
  int   *ccer;
  int   *tccer;
  void  *options;
  float **learning_x; /* learning vectors */
  float **validation_x; /* validation vectors */
  int   **sxc;
  int   **lxc;
  FILE  *outdev;
  FILE  *fdbg = (FILE *) 0;
  unsigned int   seed;
  struct dataset *learning_dset;
  struct dataset *validation_dset;

  int   mlp_nlayers; 
  int   *mlp_npl;
  int   mlp_itmax;
  float mlp_range; 
  struct mlp_options *mlp_optional;

  int    k;
  struct knn_options *knn_optional;
  
  struct svm_parameter *parameters;
  char  *func = "p_xlearn();";

  indet = 0;
  mtd = method;
  options = (void **) 0;
  parameters = (struct svm_parameter *) 0;
/*
  Accept input data.
*/
  if (tds)
    {
      status = 0;
      c = tds->c;
      outdev = stdout;
      if (dbg > 0)
	fdbg = fopen(PCP_DBG, "a");
      tccer = calloc(c+1, sizeof(int)); /* cumulative class-conditional error rates */
      idx = ivec_min(tds->nd, tds->c);
      if ((mtd == PALG_MLP) || (mtd == PALG_BAG_MLP) || (mtd == PALG_ADABOOST_MLP))
	{
	  if (mtd == PALG_BAG_MLP)
	    input_mlp(outdev, tds->d, tds->c, tds->nv, (int *) 0, &mlp_nlayers, &mlp_npl, 
		      &mlp_itmax, &mlp_range, &opt_method, &eta, &mu, &nxval, idx, &nmodels,
		      &seed, &xval_fname);
	  else
	    input_mlp(outdev, tds->d, tds->c, tds->nv, (int *) 0, &mlp_nlayers, &mlp_npl, 
		      &mlp_itmax, &mlp_range, &opt_method, &eta, &mu, &nxval, idx, (int *) 0,
		      &seed, &xval_fname);
	  mlp_optional = calloc(1, sizeof(struct mlp_options));
	  mlp_optional->nlayers = mlp_nlayers;
	  mlp_optional->npl = mlp_npl;
	  mlp_optional->itmax = mlp_itmax;
	  mlp_optional->range = mlp_range;
	  mlp_optional->eta = eta;
	  mlp_optional->opt_method = opt_method;
	  options = mlp_optional;
	}
      else if ((mtd == PALG_ADABOOST_SVM) || (mtd == PALG_SVM) || (mtd == PALG_BAG_SVM))
	{
	  parameters = calloc(1, sizeof(struct svm_parameter));
	  if (mtd == PALG_SVM)
	    status = input_svm(outdev, tds->c, tds->nv, (struct svm_problem *) 0,
			       parameters, &xval_fname, &nxval, idx, 1, PCP_SVM_K_NONE,
			       (int *) 0);
	  else
	    status = input_svm(outdev, tds->c, tds->nv, (struct svm_problem *) 0,
			       parameters, &xval_fname, &nxval, idx, 1, PCP_SVM_K_NONE,
			       &mtd);
	  if (!status)
	    {
	      if (mtd == PALG_SVM)
		nmodels = 1;
	      else if (mtd == PALG_BAG_SVM)
		nmodels = input_nmodels(stdin, outdev);
	      else if (mtd == PALG_ADABOOST_SVM)
		nmodels = boosting_nmodels(stdin, outdev);
	      seed = input_seed(stdin, outdev);
	      parameters->verbose = 1;
	      options = parameters;
	      if ((mtd == PALG_SVM) && ((struct svm_parameter *) options)->indet == 1)
		indet = 1;
	    }
	}
      else if ((method == PALG_LIN) || (method == PALG_BAG_LIN) || 
	       (method == PALG_PLC) || (method == PALG_BAG_PLC) ||
	       (method == PALG_PQC) || (method == PALG_BAG_PQC))
	{
	  clear_screen();
	  cursor_on();
	  if ((method == PALG_LIN) || (method == PALG_BAG_LIN))
	    xval_fname = input_filename(LD_MSG, PCP_XLD, outdev);
	  else
	    xval_fname = input_filename(CLASSIFIER_MSG, PCP_XPL, outdev);
	  nxval = input_nxval(stdin, outdev, idx);
	  if ((method == PALG_LIN) || (method == PALG_PLC) || (method == PALG_PQC))
	    nmodels = 1;
	  else if ((method == PALG_BAG_LIN) || (method == PALG_BAG_PLC) || (method == PALG_BAG_PQC))
	    nmodels = input_nmodels(stdin, outdev);
	  else if (method == PALG_ADABOOST_LIN)
	    nmodels = boosting_nmodels(stdin, outdev);
	  seed = input_seed(stdin, outdev);
	}
      else if ((method == PALG_KNN) || (method == PALG_BAG_KNN))
	{
	  knn_optional = calloc(1, sizeof(struct knn_options));
	  if (method == PALG_KNN)
	    status = input_knn(outdev, &nxval, idx, (int *) 0, &k, -1, &seed, &xval_fname, &dist);
	  else
	    status = input_knn(outdev, &nxval, idx, &nmodels, &k, -1, &seed, &xval_fname, &dist);
	  knn_optional->k = k;
	  knn_optional->dist = dist;
	  options = knn_optional;
	}
      if (!status)
	{
	  if (fdbg)
	    fprintf(fdbg, "%s seed: %d\n", func, seed);
	  srand(seed);
	  dflt = PCP_DFLT_NEXP;
	  jx = 1;
	  msg = malloc(1000);
	  sprintf(msg, PCP_UMSG_NEXP, dflt);
	  nexp = input_integer(stdin, outdev, msg, PCP_QLEN, &dflt, &jx, (int *) 0);
	  free(msg);

	  jx = 0;
	  ex = 1;
	  normalize = input_integer(stdin, outdev, PCP_UMSG_RAW, PCP_QLEN, &jx,
				    &jx, &ex);
	  
	  ex = 0;
	  for (jx = 0; jx < tds->c; jx++)
	    ex += (tds->nd[jx]-1)/nxval+1;
	  dr_method = input_dr(outdev, tds->nv-ex, tds->d, tds->c, &idr, &fscrit);

	  cursor_off();
	  inverse_video();
	  if ((mtd == PALG_SVM) || (mtd == PALG_BAG_SVM) || (mtd == PALG_ADABOOST_SVM))
	    status = svm_save_preamble(mtd, nmodels, nxval, nexp, tds->c, normalize, 
				       dr_method, fscrit, idr, parameters, seed, xval_fname);
	  else
	    status = save_preamble(method, nmodels, nxval, nexp, tds->c, normalize, 
				   dr_method, fscrit, idr, options, mlp_itmax, 
				   mlp_nlayers, mlp_npl, seed, xval_fname);
	  if (!status)
	    {
	      mce = 0;
	      for (ex = 0; (ex < nexp) && !status && (mce >= 0); ex++)
		{
		  amce = 0.0;
		  status = xpart(c, tds->nd, nxval, &sxc, &lxc);
		  /*
		    lxc[i][j] is number of class i vectors in
		    cross-validation subset j.
		    
		    sxc[i] is list of class i vectors with arranged so
		    that cross-validation subsets are in sequence.  For
		    example, if there are 3 cross-validation subsets with
		    10 vectors each, sxc[i][0..9] is list of class i
		    vectors in subset 0, sxc[i][10..19] is list of class i
		    vectors in subset 1, and sxc[i][20..29] is list of
		    class i vectors in subset 2.
		  */
		  if (!status)
		    {
		      ivec_set(tccer, c, 0);
		      mce = 0;
		      for (idx = 0; (idx < nxval) && (mce >= 0) && !status; idx++)
			{
			  validation_nd = malloc(c*sizeof(int));
			  learning_nd = malloc(c*sizeof(int));
			  for (jx = 0; jx < c; jx++)
			    {
			      validation_nd[jx] = lxc[jx][idx];
			      learning_nd[jx] = tds->nd[jx]-validation_nd[jx];
			      if (fdbg)
				{
				  fprintf(fdbg, "%s learning_nd[%d]: %d\n", func, jx, learning_nd[jx]);
				  fprintf(fdbg, "%s validation_nd[%d]: %d\n", func, jx, validation_nd[jx]);
				}
			    }
			  log_ites(fdbg, idx, c, tds->nv, lxc, sxc);
			  extract_sets(tds->x, idx, c, tds->nd, tds->d, sxc, lxc, 
				       &learning_x, &validation_x, fdbg);
			  log_tt(fdbg, func, c, tds->d, learning_nd, validation_nd, 
				 learning_x, validation_x);
			  nvec = ivec_sum(validation_nd, c);
			  validation_dset = dataset_lt(tds->d, c, validation_nd, 
						       nvec, (char **) 0, validation_x);
			  nvec = ivec_sum(learning_nd, c);
			  learning_dset = dataset_lt(tds->d, c, learning_nd, nvec,
						     (char **) 0, learning_x);
			  /*
			    NOTE: in November 2002 decided to fix the
			    size of each bagging set to the
			    cardinality of the training data
			    set. Perhaps revisit.
			  */
			  bag_size = learning_dset->nv;
			  mce = xlearn(method, dr_method, idr, fscrit, normalize, 
				       learning_dset, validation_dset, seed, nmodels, 
				       bag_size, options, &ccer, outdev, fdbg, errc);
			  if (mce >= 0)
			    {
			      for (jx = 0; jx < c+1; jx++)
				tccer[jx] += ccer[jx];
			      free(ccer);
			      amce += mce;
			      nvec = validation_dset->nv;
			      fmce = 100.0*((float) mce)/nvec; /* mce for this subset */
			      inverse_video();
			      fprintf(outdev, "experiment %5d; cross-validation subset%3d; error rate: "
				      "%5.2f%% (%5d/%5d)\n", ex+1, idx+1, fmce, mce, nvec);
			      if (fdbg)
				fprintf(fdbg, "experiment %5d; cross-validation subset%3d; error rate: "
					"%5.2f%% (%5d/%5d)\n", ex+1, idx+1, fmce, mce, nvec);
			      status = xlearn_save(ex+1, idx+1, tds->d, seed, amce, c, 
						   tccer, indet, lxc, xval_fname);
			      if ((status == -1) && errc)
				{
				  *errc = errno;
				  *xname = strdup(xval_fname);
				}
			    }
			  dataset_free(validation_dset);
			  dataset_free(learning_dset);
			}
		      /*sleep(1);*/
		      mx_free((void **) sxc, c);
		      mx_free((void **) lxc, c);
		      if ((status == 0) && (mce >= 0))
			{
			  xmce = 100.0*amce/tds->nv; /* cumulative average mce */
			  fprintf(outdev, "\nAverage misclassification error rate: "
				  "%5.2f%%                                  \n", xmce);
			}
		    }
		  else
		    *errc = errno;
		}
	      free(tccer);
	    }
	  else
	    {
	      *xname = strdup(xval_fname);
	      *errc = status;
	    }
	}
      else
	{
	  *errc = errno;
	  *xname = strdup(xval_fname);
	}
      if (!*errc)
	pwait();
      fflush(stdout);
      reset_video();
    }
  else
    *errc = PERR_UNDEFINED_TDS;
}



syntax highlighted by Code2HTML, v. 0.9.1