/*---------------------------------------------------------------------------
 bsvm2_mex.c: MEX-file for multi-class B-SVM with L2-soft margin.

 Compile: 
  mex new_bsvm2_mex.c qpcsolver.c kernel_fun.c

 Synopsis:
  [Alpha,bias,exitflag,kercnt,access,trnerr,t,NA,UB,LB,History] = 
     new_bsvm2_mex(data,labels,ker,arg,C,solver,tmax,tolabs,tolrel,cache,verb)

 Input:
  data [dim x num_data] Training vectors.
  labels [1 x num_data] Labels.
  ker [string] Kernel identifier.
  arg [1 x nargs] Kernel argument.
  C [1x1] Regularization constant.
  solver [string] Solver; options are 'mdm','imdm','iimdm'
     'keerthi','kowalczyk'.
  tmax [1x1] Maximal number of iterations.
  tolabs [1x1] Absolute tolerance stopping condition.
  tolrel [1x1] Relaitve tolerance stopping condition.
  cache [1x1] Number of columns of kernel matrix to be cached.
    It takes cache*num_data*size(double) bytes of memory.
  verb [1x1] If 1 then some info about the training is printed.

 Output:
  Alpha [nclass x num_data] Weights.
  bias [1x1] Bias.
  exitflag [1x1] Indicates which stopping condition was used:
    UB-LB <= tolabs           ->  exit_flag = 1   Abs. tolerance.
    (UB-LB)/(LB+1) <= tolrel  ->  exit_flag = 2   Relative tolerance.
    t >= tmax                 ->  exit_flag = 0   Number of iterations.
  kercnt [1x1] Number of kernel evaluations.
  access [1x1] Access to the virtual kernel matrix.
  trnerr [1x1] Training error.
  t [1x1] Number of iterations.
  NA [1x1] Number of non-zero alphas returned by the QP solver.
  UB [1x1] Upper bound on the optimal solution.
  LB [1x1] Lower bound on the optimal solution.
  History [2x(t+1)] UB and LB with respect to number of iterations.

  About: Statistical Pattern Recognition Toolbox
  (C) 1999-2004, Written by Vojtech Franc and Vaclav Hlavac
  <a href="http://www.cvut.cz">Czech Technical University Prague</a>
  <a href="http://www.feld.cvut.cz">Faculty of Electrical Engineering</a>
  <a href="http://cmp.felk.cvut.cz">Center for Machine Perception</a>

 Modifications:
 28-nov-2004, VF
 26-nov-2004, VF
 24-nov-2004, VF
 20-nov-2004, VF
 31-may-2004, VF
 25-jan-2003, VF
 24-jan-2003, VF
 23-jan-2003, VF
-------------------------------------------------------------------- */

#include "mex.h"
#include "matrix.h"
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <limits.h>

#include "kernel_fun.h"

#define INDEX(ROW,COL,DIM) ((COL*DIM)+ROW)

#define MINUS_INF INT_MIN
#define PLUS_INF  INT_MAX

#define KDELTA(A,B) (A==B)
#define KDELTA4(A1,A2,A3,A4) ((A1==A2)||(A1==A3)||(A1==A4)||(A2==A3)||(A2==A4)||(A3==A4))


/* Declaration of global variables */

unsigned long access_cnt;
long num_classes;      
long num_virt_data;  /* number of virtual "single-class" examples */
long num_data;       /* number of input training examples */
double kernel_diag;  /* regularization constant */
double *labels;      /* Pointer to labels */

long Cache_Size;     /* number of cached columns (min 1) */

/* cache (FIFO) for columns of the kernel matrix */
long *cache_index;                  /* indices cached of kernel columns */
long first_kernel_inx;              /* index of first inserted column */
double **kernel_columns;            /* pointers at cached columns */

/* cache for three columns of the virtual kernel matrix */
int first_virt_inx;                 /* index of first used column */
double *virt_columns[3];            /* cache for three columns*/

/* ------------------------------------------------------------
  Returns pointer at a-th column of the kernel matrix.
  This function maintains FIFO cache of kernel columns.
------------------------------------------------------------ */
void *get_kernel_col( long a ) 
{
  double *col_ptr;
  long i;
  long inx;

  inx = -1;
  for( i=0; i < Cache_Size; i++ ) {
    if( cache_index[i] == a ) { inx = i; break; }
  }
    
  if( inx != -1 ) {
    col_ptr = kernel_columns[inx];
    return( col_ptr );
  }
   
  col_ptr = kernel_columns[first_kernel_inx];
  cache_index[first_kernel_inx] = a;

  first_kernel_inx++;
  if( first_kernel_inx >= Cache_Size ) first_kernel_inx = 0;

  for( i=0; i < num_data; i++ ) {
    col_ptr[i] = kernel(i,a);
  }

  return( col_ptr );
}

/* ------------------------------------------------------------
  Computes index of input example and its class label from 
  index of virtual "single-class" example.
------------------------------------------------------------ */
void get_indices2( long *index, long *class, long i )
{
   *index = i / (num_classes-1);
 
   *class = (i % (num_classes-1))+1;
   if( *class >= labels[ *index ]) (*class)++;

   return;
}

/* ------------------------------------------------------------
  Retures (a,b)-th element of the virtual kernel matrix 
  of size [num_virt_data x num_virt_data]. 
------------------------------------------------------------ */
double kernel_fce( long a, long b )
{
  double value;
  long i1,c1,i2,c2;

  get_indices2( &i1, &c1, a );
  get_indices2( &i2, &c2, b );

  if( KDELTA4(labels[i1],labels[i2],c1,c2) ) {
    value = (+KDELTA(labels[i1],labels[i2]) 
             -KDELTA(labels[i1],c2)
             -KDELTA(labels[i2],c1)
             +KDELTA(c1,c2)
            )*(kernel( i1, i2 )+1);
  }
  else
  {
    value = 0;
  }

  if(a==b) value += kernel_diag; 

  return( value );
}

/* ------------------------------------------------------------
  Returns pointer at the a-th column of the virtual K matrix.

  (note: the b-th column must be preserved in the cache during 
   updating but b is from (a(t-2), a(t-1)) where a=a(t) and
   thus FIFO with three columns does not have to take care od b.)
------------------------------------------------------------ */
void *get_col( long a, long b )
{
  long i;
  long inx;
  long min_usage; 
  double *col_ptr;
  double *ker_ptr;
  double value;
  long i1,c1,i2,c2;

  access_cnt = access_cnt + num_virt_data;

  col_ptr = virt_columns[first_virt_inx++];
  if( first_virt_inx >= 3 ) first_virt_inx = 0;

  get_indices2( &i1, &c1, a );
  ker_ptr = (double*) get_kernel_col( i1 );

  for( i=0; i < num_virt_data; i++ ) {
    get_indices2( &i2, &c2, i );

    if( KDELTA4(labels[i1],labels[i2],c1,c2) ) {
      value = (+KDELTA(labels[i1],labels[i2]) 
               -KDELTA(labels[i1],c2)
               -KDELTA(labels[i2],c1)
               +KDELTA(c1,c2)
              )*(ker_ptr[i2]+1);
    }
    else
    {
      value = 0;
    }

    if(a==i) value += kernel_diag; 

    col_ptr[i] = value;
  }
  
  return( col_ptr );
}


/* -------------------------------------------------------------------
 Main MEX function - interface to Matlab.
-------------------------------------------------------------------- */
void mexFunction( int nlhs, mxArray *plhs[],int nrhs, const mxArray*prhs[] )
{
  char solver[20];   /* solver identifier */
  int exitflag;      /* output arg */
  int *err_bit;      /* axiliary cache for computation of trn errors*/
  int buf_len;       /* real length of the solver identifier */