svm_c.cpp

支持向量机(4)mySVM

#include "svm_c.h"

/**
*
* svm_c
*
*/

svm_c::svm_c(){
	// initialise arrays
	sum =0;
	primal=0;
	which_alpha = 0;
	lambda_eq=0;
	sum_alpha = 0;
	at_bound=0;
	all_alphas=0;
	all_ys=0;
	working_set=0;
	working_set_values=0;
	time_init=0;
	time_optimize=0;
	time_convergence=0;
	time_update=0;
	time_calc=0;
	qp.m = 1;
};


void svm_c::init(kernel_c* new_kernel, parameters_c* new_parameters){
	sum =0;
	primal=0;
	which_alpha = 0;
	lambda_eq=0;
	at_bound=0;
	target_count=0;
	kernel = new_kernel;
	parameters = new_parameters;
	Cpos = parameters->get_Cpos();
	Cneg = parameters->get_Cneg();
	is_zero = parameters->is_zero; 
	is_pattern = parameters->is_pattern;
	epsilon_pos = parameters->epsilon_pos;
	epsilon_neg = parameters->epsilon_neg;
	working_set_size = parameters->working_set_size;
	sum_alpha = 0;
	convergence_epsilon = parameters->convergence_epsilon;
	lambda_threshold = 100;
	lambda_factor=1;
	feasible_epsilon = 1;
	shrink_const = parameters->shrink_const;
	time_init=0;
	time_optimize=0;
	time_convergence=0;
	time_update=0;
	time_calc=0;
	qp.m = 1;
};

void svm_c::init_optimizer(){
	SVMINT i;
	if(0 != sum) delete []sum;
	if(0 != which_alpha) delete []which_alpha;
	if(0 != at_bound) delete []at_bound;
	sum = new SVMFLOAT[examples_total];
	at_bound = new SVMINT[examples_total];
	
	// init variables
	if(working_set_size>examples_total) working_set_size = examples_total;
	
	qp.n = working_set_size;
	qp.c = new SVMFLOAT[qp.n];
	qp.H = new SVMFLOAT[qp.n*qp.n];
	qp.A = new SVMFLOAT[qp.n];
	qp.b = new SVMFLOAT[qp.m];
	qp.l = new SVMFLOAT[qp.n];
	qp.u = new SVMFLOAT[qp.n];
	
	which_alpha = new SVMINT[working_set_size];
	primal = new SVMFLOAT[qp.n];
	// reserve workspace for calculate_working_set
	working_set = new SVMINT[working_set_size];
	working_set_values = new SVMFLOAT[working_set_size];
	
	for(i=qp.n;i<qp.n;i++){
		qp.A[i] = 1;
	};
	//      qp.l[i] = 0 done in svm_pattern_c::
	for(i=0;i<working_set_size;i++){
		qp.l[i] = -is_zero; 
	};
	
	Cpos /= (SVMFLOAT)examples_total;
	Cneg /= (SVMFLOAT)examples_total;
	
	if((parameters->quadraticLossPos) && (parameters->quadraticLossNeg)){
		Cpos = infinity;
		Cneg = infinity;
	}
	else if(parameters->quadraticLossPos){
		Cpos = infinity;
	}
	else if(parameters->quadraticLossNeg){
		Cneg = infinity;
	};
	//  sigfig_max = -log10(is_zero);
	lambda_WS = 0;
	to_shrink=0;
	
	smo.init(parameters->is_zero,parameters->convergence_epsilon,working_set_size*working_set_size);
};


void svm_c::exit_optimizer(){
	delete [](qp.c);
	delete [](qp.H);
	delete [](qp.A);
	delete [](qp.b);
	delete [](qp.l);
	delete [](qp.u);
	
	delete []primal;
	delete []working_set;
	delete []working_set_values;
	delete []sum;
	delete []at_bound;
	delete []which_alpha;
	
	primal = 0;
	working_set = 0;
	working_set_values = 0;
	sum=0;
	at_bound=0;
	which_alpha=0;
};


int svm_c::is_alpha_neg(const SVMINT i){
	// variable i is alpha*
	
	// take a look at svm_pattern_c::is_alpha_neg 
	// and svm_regression_c::is_alpha_neg!
	
	int result=0;
	
	if(is_pattern){
		if(all_ys[i] > 0){
			result = 1;
		}
		else{
			result = -1;
		};
	}
	else if(all_alphas[i] > 0){
		result = 1;
	}
	else if(all_alphas[i] == 0){
		result = 2*(i%2)-1;
	}
	else{
		result = -1;
	};
	return result;
};


SVMFLOAT svm_c::nabla(const SVMINT i){
	if(is_alpha_neg(i) > 0){
		return( sum[i] - all_ys[i] + epsilon_neg);
	}
	else{
		return(-sum[i] + all_ys[i] + epsilon_pos);
	};
};


SVMFLOAT svm_c::lambda(const SVMINT i){
	// size lagrangian multiplier of the active constraint
	
	SVMFLOAT alpha;
	SVMFLOAT result = -abs(nabla(i)+is_alpha_neg(i)*lambda_eq);
    //= -infinity; // default = not at bound
	
	alpha=all_alphas[i];
	
	if(alpha>is_zero){
		// alpha*
		if(alpha-Cneg >= - is_zero){
			// upper bound active
			result = -lambda_eq-nabla(i);
		};
	}
	else if(alpha >= -is_zero){
		// lower bound active
		if(is_alpha_neg(i) > 0){
			result = nabla(i) + lambda_eq;
		}
		else{
			result = nabla(i)-lambda_eq;
		};
	}
	else if(alpha+Cpos <= is_zero){
		// upper bound active
		result = lambda_eq - nabla(i);
	};
	
	return result;
};


int svm_c::feasible(const SVMINT i, SVMFLOAT* the_nabla, SVMFLOAT* the_lambda, int* atbound){
	// is direction i feasible to minimize the target function
	// (includes which_alpha==0)
	
	int is_feasible=1;
	
	//  if(at_bound[i] >= shrink_const){ is_feasible = 0; };
	
	SVMFLOAT alpha;
	
	*the_nabla = nabla(i);
	*the_lambda = lambda(i);
	
	alpha=all_alphas[i];
	
	if(alpha-Cneg >= - is_zero){
		// alpha* at upper bound
		*atbound = 1;
		if(*the_lambda >= lambda_threshold){
			at_bound[i]++;
			if(at_bound[i] == shrink_const) to_shrink++;
		}
		else{
			at_bound[i] = 0;
		};
	}
	else if((alpha<=is_zero) && (alpha >= -is_zero)){
		// lower bound active
		*atbound = 1;
		if(*the_lambda >= lambda_threshold){
			at_bound[i]++;
			if(at_bound[i] == shrink_const) to_shrink++;
		}
		else{
			at_bound[i] = 0;
		};
	}
	else if(alpha+Cpos <= is_zero){
		// alpha at upper bound
		*atbound = 1;
		if(*the_lambda >= lambda_threshold){
			at_bound[i]++;
			if(at_bound[i] == shrink_const) to_shrink++;
		}
		else{
			at_bound[i] = 0;
		};
	}
	else{
		// not at bound
		*atbound = 0;
		at_bound[i] = 0;
	};
	if((*the_lambda >= -feasible_epsilon) || (at_bound[i] >= shrink_const)){
		is_feasible = 0; 
	};
	return is_feasible;
};


int svm_c::feasible(const SVMINT i){
	// is direction i feasible to minimize the target function
	// (includes which_alpha==0)
	
	if(at_bound[i] >= shrink_const){ return 0; };
	
	SVMFLOAT alpha;
	SVMFLOAT result;
	
	alpha=all_alphas[i];
	if(alpha-Cneg >= - is_zero){
		// alpha* at upper bound
		result = -lambda_eq-nabla(i);
		if(result>=-feasible_epsilon){
			return 0; 
		};
	}
	else if((alpha<=is_zero) && (alpha >= -is_zero)){
		// lower bound active
		if(is_alpha_neg(i) > 0){
			result = nabla(i) + lambda_eq;
		}
		else{
			result = nabla(i)-lambda_eq;
		};
		if(result>=-feasible_epsilon){
			return 0; 
		};
	}
	else if(alpha+Cpos <= is_zero){
		// alpha at upper bound
		result = lambda_eq - nabla(i);
		if(result>=-feasible_epsilon){
			return 0; 
		};
	}
	else{
		// not at bound
		result= abs(nabla(i)+(SVMFLOAT)(is_alpha_neg(i))*lambda_eq);
		if(result<=feasible_epsilon){
			return 0; 
		};
	};
	return 1;
};



void svm_c::set_svs(example_set_c* training_examples){
	// initialised already trained sv (for predicting or testing)
	
	examples = training_examples; 
	test_set = training_examples;
	examples_total = training_examples->size();
	all_alphas = examples->get_alphas();
};


void svm_c::reset_shrinked(){
	SVMINT old_ex_tot=examples_total;
	target_count=0;
	examples_total = examples->size();
	// unshrink, recalculate sum for all variables 
	SVMINT i,j;
	// reset all sums
	SVMFLOAT Kij;
	for(i=old_ex_tot;i<examples_total;i++){
		sum[i] = 0;
		at_bound[i] = 0;
	};
	for(i=0;i<examples_total;i++){
		if(abs(all_alphas[i])>is_zero){
			for(j=old_ex_tot;j<examples_total;j++){
				Kij = kernel->calculate_K(i,j);
				sum[j]+=all_alphas[i]*Kij;
			};
		};
	};
	sum_alpha=0;
};


svm_result svm_c::train(example_set_c* training_examples){
	svm_result the_result;
	examples = training_examples; 
	test_set = training_examples;
	if(parameters->verbosity>= 4){
		cout<<"training started"<<endl;
	};
	time_init = get_time();
	time_all = time_init;
	
	examples_total = training_examples->size();
	all_alphas = examples->get_alphas();
	all_ys = examples->get_ys();
	SVMINT param_wss = parameters->working_set_size;
	if(param_wss > examples_total){
		parameters->working_set_size = examples_total;
		working_set_size = examples_total;
	};
	
	// set up data structure for optimizer
	if(parameters->verbosity>= 4){
		cout<<"Initing optimizer"<<endl;
	};
	init_optimizer();
	
	SVMINT iteration = 0;
	SVMINT max_iterations = parameters->max_iterations;
	
	if(parameters->verbosity>= 4){
		cout<<"Initing working set"<<endl;
	};
	
	// WS-optimization
	init_working_set();
	
	time_init = get_time() - time_init;
	while(iteration < max_iterations){
		iteration++;
		if(parameters->verbosity>= 5){
			cout<<"WS is: ";
			//      SVMINT uyt;
			SVMINT ws_i;
			for(ws_i=0;ws_i<working_set_size;ws_i++){
				if(! feasible(working_set[ws_i])) cout<<"!";
				if(is_alpha_neg(working_set[ws_i])>0) cout<<"+"<<working_set[ws_i]<<" ";
				if(is_alpha_neg(working_set[ws_i])<0) cout<<"-"<<working_set[ws_i]<<" ";
			};
			cout<<endl;
			cout<<"WS alphas: ";
			for(ws_i=0;ws_i<working_set_size;ws_i++){
				cout<<all_alphas[working_set[ws_i]];
				if((all_alphas[working_set[ws_i]] > Cneg-is_zero) ||
					(all_alphas[working_set[ws_i]] < -Cpos+is_zero)){
					cout<<"^";
				};
				cout<<" ";
			};
			cout<<endl;
			cout<<"WS nablas: ";
			for(ws_i=0;ws_i<working_set_size;ws_i++){
				cout<<nabla(working_set[ws_i])<<" ";
				
			};
			cout<<endl;
			cout<<"WS lambdas: ";
			for(ws_i=0;ws_i<working_set_size;ws_i++){
				cout<<lambda(working_set[ws_i])<<" ";
				
			};
			cout<<endl;
			cout<<"WS at_bounds: ";
			for(ws_i=0;ws_i<working_set_size;ws_i++){
				cout<<at_bound[working_set[ws_i]]<<" ";
				
			};
			cout<<endl;
		};
		
		if(parameters->verbosity>= 4){
			cout<<"optimizer iteration "<<iteration<<endl;
		}
		else if(parameters->verbosity>=3){
			cout<<".";
			cout.flush();
		};
		
		optimize(); 
		put_optimizer_values();
		
		//     SVMFLOAT now_target_dummy=0;
		//     for(SVMINT i=0;i<examples_total;i++){
		//       now_target_dummy=sum[i]/2-all_ys[i];
		//       if(is_alpha_neg(i)){
		// 	now_target_dummy+= epsilon_pos;
		//       }
		//       else{
		// 	now_target_dummy-= epsilon_neg;
		//       };
		//       now_target+=all_alphas[i]*now_target_dummy;
		//     };
		//     cout<<"Target function: "<<now_target<<endl;
		
		int conv = convergence();
		
		if(conv && (is_zero > parameters->is_zero)){ 
			is_zero =  parameters->is_zero;
			conv=0;
		}
		else if(conv){ 
			if(parameters->verbosity>=3){
				// dots
				cout<<endl;
			};
			// check convergence for all parameters
			if(0 == is_pattern){
				for(SVMINT i=0;i<examples_total;i++){
					if((all_alphas[i]<=is_zero) && (all_alphas[i]>=-is_zero)){
						all_alphas[i]=0;
					};
				};
				project_to_constraint(); 
			};
			if((examples_total<examples->size()) || (target_count>0)){
				if(conv){
					// check convergence for all alphas  
					if(parameters->verbosity>= 2){
						cout<<"***** Checking convergence for all variables"<<endl;
					};
					SVMINT old_target_count=target_count; // why???
					reset_shrinked();
					conv = convergence();
					if(0 == conv){
						kernel->set_examples_size(examples_total);
						target_count = old_target_count;
					};
				};
			};
			if((0 == is_pattern) && (conv)){
				for(SVMINT i=0;i<examples_total;i++){
					if((all_alphas[i]<=is_zero) && (all_alphas[i]>=-is_zero)){
						at_bound[i]++;
					};
				};
				conv = convergence();
			};
			
			if(conv){
				if(parameters->verbosity>= 1){
					cout<<"*** Convergence"<<endl;
				};
				if((parameters->verbosity>=2) || 
					(convergence_epsilon > parameters->convergence_epsilon)){
					// had to relax the KKT conditions on the way. 
					// Maybe this isn't necessary any more
					SVMFLOAT new_convergence_epsilon = 0;
					SVMFLOAT the_lambda;
					for(SVMINT i=0;i<examples_total;i++){
						the_lambda = lambda(i);
						if(the_lambda < new_convergence_epsilon){
							new_convergence_epsilon = the_lambda;
						};
					};
					convergence_epsilon = -new_convergence_epsilon;
				};
				
				break;
			};
			//       else if(parameters->verbosity>= 2){
			// 	// shrinking did not work, be more careful
			// 	cout<<"*** convergence error due to shrinking in iteration "<<iteration<<", restarting"<<endl;
			//       };