/* mysmo.h: interface for the mysmo class.*/

/*******************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include <time.h>
#include <math.h>

	int	  m_end_support_i; /* Support vectors are from 0...end_support_i */
	float m_C=1000; /* Penalty factor (HIGH->high penalty for wrong side) */
	float m_tolerance=0.001f; /* Tolerance in KKT condition */
	float m_eps=0.001f; /* Break limit */
	float m_b; /* Threshold (begin at 0) */
	float * m_x; /* The training points	*/
	float * m_y; /* Classification of training points */
	float * m_a; /*agrange multipliers (LM) */
	float * m_w; /* Weight vector */
	float * m_self_dot_product; /* vector of dot prod x*x for all x */
	float * m_e; /* vector of prediction error Ei */
	float   m_sigma;
	float	m_two_sigma_squared;
	float * m_testIns; /* The training points */
	float * m_errInsIndex;/* testing err instance's index */
	float * m_max;
	float * m_min;	
	float precision;
	unsigned int m_testNum;	/* number of testing points */
	unsigned int m_dimension; /* dimension of points */	
	unsigned int m_xNum; /* Number of training points */

	int		indicator_function(float *v);/* Returns 1 or -1 as class. for v */
	float	decision_function(float * v);/* Calculates the decision func. at v */
	int		testing(char * fileName,unsigned int number);
	int		LMToFile(char * fileName);
	int		pointToFile(char * FileName,float * instance ,float * y,unsigned int number);
	int		skipToFile(char * FileName,float * x ,float * y,unsigned int number);
	int		skip(char * FileName,unsigned int number,unsigned int dimession);
	unsigned char	getFileName(char * dest,const char * source);
	int		readPoint(char * FileName,unsigned int number,unsigned int dimession);
	void	initCsmo();
	int		allocMemory(unsigned int number,unsigned int dimession);
	float	drand48();
	float	dot_product(int i1, int i2);/* Dot product x*y */
	float	kernel(int i1, int i2);/* The kernel function */
	float	learned_func(int k);/* The learned function w*x - b= sum(ai*yi*ki) -b */
	void	RenormalizeData(float * instance,unsigned int number);
	int		optimize( int i1, int i2 );
	int		tryLM(int i1);
	int		createSmo();
	int		scaleToFile(char * FileName,float * x ,float * y,unsigned int number);
	int		scale(char * FileName,unsigned int number,unsigned int dimession);
	void	freeCsmo();


void initCsmo()
{
	m_a=NULL; /* Lagrange multipliers (LM) */
	m_w=NULL; /* Weight vector */
	m_y=NULL; /* Classification of training points */
	m_self_dot_product=NULL; /* vector of dot prod x*x for all x */
	m_e=NULL; /* vector of prediction error Ei */
	m_x=NULL; /* The training points */
	m_testIns=NULL; /* The training points */
	m_errInsIndex=NULL;/* testing err instance's index */
	m_max=NULL;
	m_min = NULL;
	m_xNum=0; /* Number of training points */
	m_testNum=0;	/* number of testing points */
	m_dimension=0; /* dimension of points */
	m_end_support_i=0; /* Support vectors are from 0...end_support_i	*/
	
	m_b = 0; /* Threshold (begin at 0) */
	m_two_sigma_squared = 2 * m_sigma * m_sigma;
}

int allocMemory(unsigned int number,unsigned int dimession)
{
	m_x = (float*)malloc(sizeof(float)*number * dimession);
	m_y = (float*)malloc(sizeof(float)*number);
	m_a = (float*)malloc(sizeof(float)*number); /* Lagrange multipliers (LM) */
	memset(m_a,0,number*sizeof(float));
	m_e  = (float*)malloc(sizeof(float)*number); /* vector of prediction error Ei */
	memset(m_e,0,number*sizeof(float));
	m_self_dot_product = (float*)malloc(sizeof(float)*number); /* vector of dot prod x*x for all x */
	m_max = (float*)malloc(sizeof(float)*dimession);
	m_min = (float*)malloc(sizeof(float)*dimession);
	return 0;
}

void freeCsmo()
{
	if(m_a) free(m_a);
	if(m_w) free(m_w);
	if(m_y) free( m_y);
	if(m_self_dot_product) free(m_self_dot_product);
	if(m_e) free(m_e);
	if(m_x) free(m_x);
	if(m_testIns) free(m_testIns);
	if(m_errInsIndex) free(m_errInsIndex);
	if(m_max) free(m_max);
	if(m_min) free(m_min);
}

float drand48()
{
	return rand()/(float)RAND_MAX;
}

float dot_product(int i1, int i2)
{
/* 	Dot product x*y */
	unsigned int i;
	float dot = 0.;
	float *p1 = m_x + i1*m_dimension;
	float *p2 = m_x + i2*m_dimension;
	for(i = 0; i < m_dimension; i++)
		dot += p1[i] * p2[i];
	return dot;
}
float kernel(int i1, int i2)/* The kernel function */
{
	float s = dot_product(i1, i2); /* Dot product x*y */
	s *= -2;
	s += m_self_dot_product[i1] + m_self_dot_product[i2];
	return exp( -s / m_two_sigma_squared ); /* Gaussian RBF */
}
float learned_func(int k)/* The learned function w*x - b= sum(ai*yi*ki) -b */
{
	int i;
	float s = 0.;
	for(i = 0; i < m_end_support_i; i++)
	{
		if(m_a[i] > 0) s += m_a[i] * m_y[i] * kernel(i,k);
	}
	s -= m_b;
	return s;
}

int	optimize( int i1, int i2 )
{
	int i;
	float y1, y2, s;
	float a1o, a2o; /* old LM values */
	float a1n, a2n; /* new LM values */
	float E1, E2, L, H, eta;
	float k11,k12,k22;
	float diff;
	float c1,c2;
	float Lobj,Hobj;
	float t1,t2;
	/* Update threshold b to reflect change in a */
	float b1, b2, bnew, deltab;
	/* Assignment of values */
	if(i1 == i2) return 0;
	a1o = m_a[i1];
	a2o = m_a[i2];
	y1 = m_y[i1];
	y2 = m_y[i2];
	if(a1o > 0 && a1o < m_C) E1 = m_e[i1];
	else E1 = learned_func(i1) - y1;
	if((a2o > 0) && (a2o < m_C)) E2 = m_e[i2];
	else E2 = learned_func(i2) - y2;
	s = y1 * y2;
	/* Computation of L and H, low and high end a2n on line segment */
	if( y1 == y2 ) {
		float sum = a1o + a2o;
		if(sum > m_C){
			L = sum - m_C;
			H = m_C;
		}
		else{
			L = 0;
			H = sum;
		}
	}
	else{
		diff = a1o - a2o;
		if(diff > 0){
			L = 0;
			H = m_C - diff;
		}
		else{
			L = -diff;
			H = m_C;
		}
	}
	if(L == H) return 0; /* Error, no line segment */
	/* Computation of eta */
	k11 = kernel(i1,i1);
	k22 = kernel(i2,i2);
	k12 = kernel(i1,i2);
	eta = 2 * k12 - k11 - k22; /* Computation of eta from kernel */
	if(eta < -0.01){
		a2n = a2o + y2 * (E2 - E1) / eta; /* New a2 on line segment */
		if( a2n < L ) a2n = L; /* Lowest a2n on line segment */
		else if( a2n > H ) a2n = H; /* Highest a2n on line segment */
	}
	else
	{ /* eta~0 */
		c1 = eta/2;
		c2 = y2 * (E1-E2) - eta * a2o;
		Lobj = c1 * L * L + c2 * L;
		Hobj = c1 * H * H + c2 * H;
		if( Lobj > Hobj + m_eps ) a2n = L;
		else if( Lobj < Hobj - m_eps ) a2n = H;
		else a2n = a2o;
	}
	if( fabs(a2n - a2o) < (m_eps * (a2n + a2o + m_eps)) )
		return 0;
	a1n = a1o - s * (a2n - a2o); /* New a1 */
	if(a1n < 0)
	{ /* Disallowed case */
		a2n += s * a1n;
		a1n = 0; /* Limit case */
	}
	else if(a1n > m_C)
	{ /* Disallowed case */
		a2n += s * (a1n - m_C);
		a1n = m_C; /* Limit case */
	}

	
	if((a1n > 0) && (a1n < m_C))
	{
		bnew = m_b + E1 + y1 * (a1n - a1o) * k11 + y2 * (a2n - a2o) * k12;
	}
	else
	{
		if((a2n > 0) && (a2n < m_C))
		{
			bnew = m_b + E2 + y1 * (a1n - a1o) * k12 + y2 * (a2n - a2o) * k22;
		}
		else
		{
			b1 = m_b + E1 + y1 * (a1n - a1o) * k11 + y2 * (a2n - a2o) * k12;
			b2 = m_b + E2 + y1 * (a1n - a1o) * k12 + y2 * (a2n - a2o) * k22;
			bnew = (b1 + b2) / 2;
		}
	}
	deltab = bnew - m_b;
	m_b = bnew;
	/* Update e vector */
	t1 = y1 * (a1n - a1o);
	t2 = y2 * (a2n - a2o);
	for(i = 0; i < m_end_support_i; i++)
	{
		if((0 < m_a[i]) && (m_a[i] < m_C))
		{
			m_e[i] += t1 * kernel(i1,i) + t2 * kernel(i2,i) - deltab;
		}
	}
	m_e[i1] = 0.0;
	m_e[i2] = 0.0;
	m_a[i1] = a1n; /*Store the optimized pair of LMs */
	m_a[i2] = a2n;
	return 1; /*Success! */
}

int	tryLM(int i1)
{
	
	int k,k0,i2; 
	float E1,E2,y1,a1,r1,tmax,temp; /*Functional values */
	y1 = m_y[i1]; /*Classification of point i1 */
	a1 = m_a[i1]; /*LM of point i1 */
	if( (a1 > 0) && (a1 < m_C) )
	{ /*should mean r1=0 from KKT condition */
		E1 = m_e[i1];
	}
	else
	{
		E1 = learned_func(i1) - y1; /*Calculation of prediction error */
	}
	r1 = y1 * E1; /*r1 quantity of KKT condition */
	if(((r1 < -m_tolerance) && (a1 < m_C)) || ((r1 > m_tolerance) && (a1 > 0)))
	{ /*If KKT is violated */
		/*Try to find suitable example LM (index i2) */
		/*Try boundary points with the largest difference |E1 - E2| */
		for( i2 = (-1), tmax = 0, k = 0; k < m_end_support_i; k++)
		{
			if((m_a[k] > 0) && (m_a[k] < m_C))
			{ /*boundary points */
				E2 = m_e[k];
				temp = fabs(E1 - E2);
				if(temp > tmax){
					tmax = temp;
					i2 = k;
				}
			}			
		}
		if(i2 >= 0)
		{
			if(optimize( i1, i2 ) == 1) return 1;
		}
		/*Try all boundary points */
		for(k0 = (int)(drand48() * m_end_support_i), k = k0; k < m_end_support_i + k0; k++)
		{
			i2 = k % m_end_support_i;
			if((m_a[i2] > 0) && (m_a[i2] < m_C))
			{
				if(optimize(i1,i2) == 1) return 1;
			}
		}
		/*Try the entire set */
		for(k0 = (int)(drand48() * m_end_support_i), k = k0; k < m_end_support_i + k0; k++){
			i2 = k % m_end_support_i;
			if(optimize(i1,i2) == 1) return 1;
		}
	}
	return 0;
}

float decision_function(float * v)/*Calculates the decision func. at v */
{
	unsigned int i,j;
	float dec = 0;
	float s = 0; /*v - x[i] */
	float *p = m_x;
	float temp;
	for(i = 0; i < m_xNum; i++)
	{
		if(m_a[i]>0)
		{
			s = 0;
			for(j = 0; j < m_dimension; j++)
			{
				temp = v[j] - p[j];
				s = s + temp * temp;
			}
			dec += m_y[i] * m_a[i] * exp( - s / m_two_sigma_squared );
		}
			p = p + m_dimension;
	}
	dec -= m_b;
	return dec;
}

int	indicator_function(float *v)/*Returns 1 or -1 as class. for v */
{
	float dec = decision_function(v);
	if(dec > 0)
		return 1;
	else
		return -1;
}

unsigned char getFileName(char * dest,const char * source)
{
	unsigned int i,j;
	if(source==NULL)
		return 1;
	j = strlen(source) -1;
	if(j==0)
		return 1;
	for(i=j; i>0; i--)
	{
		if(source[i]=='.')
			j = i;
		else if(source[i]=='\\')
		{
			i++;
			break;
		}
	}
	strncpy(dest,source+i,j-i);
	dest[j-i]='\0';
	return 0;
}

int	testing(char * fileName,unsigned int number)
{
	unsigned int right=0,wrong=0;