/*-----------------------------------------------------------------------
 Kozinec's algorithm for separable single-class SVM problem.

 int single_kozinec(TKerFun ker,
             long num_data,
             long tmax,
             double tolabs,
             double tolrel,
             double *Alpha,
             double *UB,
             double *LB,
             long *t,
             double *History)

 tmax, tolabs, tolrel ... Define stopping conditions: 

    UB <= tolabs              ->  exit_flag = 1   Abs. tolerance.
    (UB-LB)/(LB+1) <= tolrel  ->  exit_flag = 2   Relative tolerance.
    t >= tmax                 ->  exit_flag = 0   Number of iterations.
 
 Alpha ... Lagrangians defining found decision rule. 
 UB ... Achieved upper bound on the optimal solution.
 LB ... Achieved lower bound on the optimal solution.
 t  ... Number of iterations.
 History ... Value of LB and UB wrt. number of iterations.

 Modifications:
 22-Jan-2004, VF
 14-Oct-2003, VF
-------------------------------------------------------------------- */

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

#define HISTORY_BUF 1000000

#define MINUS_INF INT_MIN
#define PLUS_INF  INT_MAX

#define ABS(A) ((A >= 0) ? A : -A)
#define MIN(A,B) ((A < B) ? A : B)
#define INDEX(ROW,COL,DIM) ((COL*DIM)+ROW)

/* ==============================================================
 Kernel Kozinec's algorithm.
============================================================== */

int single_kozinec(const double (*kernel_fce)(long, long),
            long num_data, 
            long tmax,
            double tolabs,
            double tolrel,
            double **out_Alpha,
            double *out_UB,
            double *out_LB,
            long  *out_t,
            double **out_History)
{
  double *Alpha;
  double *History;
  double LB;
  double UB2;
  double tmp;
  double kernel_diag;
  double *ProjX;
  double *K_Diag;
  double lambda;
  double *tmp_ptr;
  long min_inx;
  long new_min_inx;
  long i;
  long t;
  long History_size;
  int exitflag;

  /* allocate memory */
  Alpha = mxCalloc(num_data, sizeof(double));
  if( Alpha == NULL ) mexErrMsgTxt("Not enough memory.");

  ProjX = mxCalloc(num_data, sizeof(double));
  if( ProjX == NULL ) mexErrMsgTxt("Not enough memory.");

  K_Diag = mxCalloc(num_data, sizeof(double));
  if( K_Diag == NULL ) mexErrMsgTxt("Not enough memory.");

  History_size = (tmax < HISTORY_BUF ) ? tmax+1 : HISTORY_BUF;
  History = mxCalloc(History_size*2,sizeof(double));
  if( History == NULL ) mexErrMsgTxt("Not enough memory.");
  
  /* == inicialization == */
  for( LB =  PLUS_INF, i = 0; i < num_data; i++ ) 
  {
    Alpha[i] = 0;

    ProjX[i] = kernel_fce( 0, i );

    K_Diag[i] = kernel_fce( i, i );

    if( ProjX[i] < LB ) {
      LB = ProjX[i];
      min_inx = i;
    }
  }

  UB2 = K_Diag[0];
  LB = LB/sqrt(UB2);
  Alpha[0] = 1;
  t = 0;
  History[INDEX(0,0,2)] = LB;
  History[INDEX(1,0,2)] = sqrt(UB2);

  if( sqrt(UB2) <= tolabs ) exitflag = 1;
  else if((sqrt(UB2)-LB)/(ABS(LB)+1) <= tolrel ) exitflag = 2;
  else exitflag = -1;

  /* == Main cycle == */

  while( exitflag == -1 ) 
  {
    t++;     
    
    /* Adaptation rule and update */
    lambda = (UB2 - ProjX[min_inx])/(UB2 - 2*ProjX[min_inx] + K_Diag[min_inx]);
    lambda = MIN( lambda, 1);

    UB2 = UB2*(1-lambda)*(1-lambda) + 
          2*lambda*(1-lambda)*ProjX[min_inx]+
          lambda*lambda*K_Diag[min_inx];

    LB = PLUS_INF;
    for( i = 0; i < num_data; i++ ) 
    {
       if( Alpha[i] != 0) Alpha[i] = Alpha[i]*(1-lambda);
       ProjX[i] = ProjX[i]*(1-lambda) + lambda*kernel_fce(i,min_inx);

       if( ProjX[i] < LB ) 
       { 
         LB = ProjX[i];
         new_min_inx = i;
       }
    }    
    
    LB = LB/sqrt(UB2);
    Alpha[min_inx] = Alpha[min_inx] + lambda;
    min_inx = new_min_inx;    

    /* Stopping conditions */
    if( sqrt(UB2) <= tolabs ) exitflag = 1; 
    else if( ((sqrt(UB2)-LB)/(ABS(LB)+1)) <= tolrel ) exitflag = 2; 
    else if(t >= tmax) exitflag = 0; 

    /* Store lower and upper bounds */
    if( t < History_size ) {
      History[INDEX(0,t,2)] = LB;
      History[INDEX(1,t,2)] = sqrt(UB2);
    }
    else {
      tmp_ptr = mxCalloc((History_size+HISTORY_BUF)*2,sizeof(double));
      if( tmp_ptr == NULL ) mexErrMsgTxt("Not enough memory.");
      for( i = 0; i < History_size; i++ ) {
        tmp_ptr[INDEX(0,i,2)] = History[INDEX(0,i,2)];
        tmp_ptr[INDEX(1,i,2)] = History[INDEX(1,i,2)];
      }
      tmp_ptr[INDEX(0,t,2)] = LB;
      tmp_ptr[INDEX(1,t,2)] = sqrt(UB2);
      
      History_size += HISTORY_BUF;
      mxFree( History );
      History = tmp_ptr;
    }
  }

  /* transform Alphas to obtain canonical hyperplane representation */
  for( i = 0; i < num_data; i++ ) {
    Alpha[i] = Alpha[i] / (LB*sqrt(UB2));
  }

  /* outputs */
  (*out_Alpha) = Alpha;
  (*out_UB) = sqrt(UB2);
  (*out_LB) = LB;
  (*out_t) = t;
  (*out_History) = History;

  /**/
  mxFree( ProjX );
  mxFree( K_Diag );

  return( exitflag ); 
}