/*-----------------------------------------------------------------------
 Mitchell-Demyanov-Malozemov (MDM) algorithm for separable single-class 
 SVM problem. 

 int single_mdm(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:
 31-may-2004, VF
 23-Jan-2004, 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 MDM algorithm.
============================================================== */

int single_mdm(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, tmp1, tmp2;
  double k11, k12, k22;
  double kernel_diag;
  double *ProjX;
  double *K_Diag;
  double lambda;
  double *tmp_ptr;
  long min_inx;
  long max_inx;
  long new_min_inx;
  long new_max_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;
    }
  }

  max_inx = 0;
  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);

  /* Stopping conditions */
  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 */
    k11 = K_Diag[max_inx];
    k22 = K_Diag[min_inx];
    k12 = kernel_fce(max_inx,min_inx);

    lambda = (ProjX[max_inx]-ProjX[min_inx])/(Alpha[max_inx]*(k11 - 2*k12 + k22));
    if( lambda < 0 ) lambda = 0; else if (lambda > 1) lambda = 1;

    UB2 = UB2 + 2*lambda*Alpha[max_inx]*(ProjX[min_inx]-ProjX[max_inx]) + 
          lambda*lambda*Alpha[max_inx]*Alpha[max_inx]*(k11 - 2*k12 + k22);

    tmp1 = Alpha[max_inx];
    Alpha[min_inx]=Alpha[min_inx]+lambda*Alpha[max_inx];
    Alpha[max_inx]=Alpha[max_inx]*(1-lambda);

    LB = PLUS_INF;
    tmp2 = MINUS_INF;
    for( i = 0; i < num_data; i++ ) 
    {

       ProjX[i] = ProjX[i] + lambda*tmp1*
          (kernel_fce(i,min_inx) - kernel_fce(i,max_inx));

       if( Alpha[i] !=0 && tmp2 < ProjX[i]) 
       {
         new_max_inx = i;
         tmp2 = ProjX[i];
       }

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

    LB = LB/sqrt(UB2);
    min_inx = new_min_inx;    
    max_inx = new_max_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 selected values */
    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 ); 
}