svm_c.cpp

vc环境下的支持向量机源码

#include "stdafx.h"
#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;
  feasible_epsilon = convergence_epsilon;
  shrink_const = parameters->shrink_const;
  time_init=0;
  time_optimize=0;
  time_convergence=0;
  time_update=0;
  time_calc=0;
  qp.m = 1;
  biased = parameters->biased;
};

void svm_c::init_optimizer(){
  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];
  if(! biased){
    qp.m = 0;
  };

  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];

  if(parameters->do_scale_y){
    epsilon_pos /= examples->get_y_var();
    epsilon_neg /= examples->get_y_var();
  };

  SVMINT i;
  //      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){
    Cpos = infinity;
  };
  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 >= 0){
      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 >= 0){
      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 >= 0){
      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;
};


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();
  kernel->set_examples_size(examples_total);
  // 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;
};


SVMFLOAT svm_c::avg_norm2(){
  SVMFLOAT avg=0.0;
  for(SVMINT i=0;i<examples_total;i++){
      avg += kernel->calculate_K(i,i);
  };
  avg /= (SVMFLOAT)examples_total;
  return(avg);
};


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;
  };
  if(parameters->realC <= 0.0){
    parameters->realC = 1.0/avg_norm2();
    Cpos = parameters->get_Cpos();
    Cneg = parameters->get_Cneg();
    if(parameters->verbosity>= 2){
      cout<<"C set to "<<parameters->realC<<endl;
    };
  };
  if(parameters->balance_cost){
    parameters->Lpos = (SVMFLOAT)training_examples->size_pos()/(SVMFLOAT)training_examples->size();
    parameters->Lneg = (SVMFLOAT)training_examples->size_neg()/(SVMFLOAT)training_examples->size();
    Cpos = parameters->get_Cpos();
    Cneg = parameters->get_Cneg();
  };

  // 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(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();

    SVMINT i;

//     SVMFLOAT now_target=0;
//     SVMFLOAT now_target_dummy=0;
//     for(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();
    
    // old, is_zero is not touched any more
    /// @@@
    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(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)){
	    // check convergence for all alphas  
	    if(parameters->verbosity>= 2){
		cout<<"***** Checking convergence for all variables"<<endl;
	    };
	    SVMINT old_target_count=target_count; // t_c set to 0 in reset_shrinked
	    reset_shrinked();
	    conv = convergence();
	    if(0 == conv){
		kernel->set_examples_size(examples_total);
		target_count = old_target_count;
	    };
	};
	
	if((0 == is_pattern) && (conv)){
	    // why???
	    for(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(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;
	};
	
	// set variables free again
	shrink_const += 10;
	for(i=0;i<examples_total;i++){
	    at_bound[i]=0;
      };
    };

    shrink();

    calculate_working_set();
    update_working_set();

    if(parameters->verbosity >= 4){
      SVMINT shrinked=0;
      SVMINT upper_bound=0;
      SVMINT lower_bound=0;

      for(i=0;i<examples_total;i++){
	if(at_bound[i] >= shrink_const){ 
	  shrinked++;
	};
	if(abs(all_alphas[i]) < is_zero){
	  lower_bound++;
	}
	else if((all_alphas[i]-Cneg>-is_zero) || (all_alphas[i]+Cpos<is_zero)){
	  upper_bound++;
	};
      };
      cout<<examples_total<<" variables total, ";
      cout<<lower_bound<<" variables at lower bound, ";
      cout<<upper_bound<<" variables at upper bound, ";
      cout<<shrinked<<" variables shrinked"<<endl;
    };
  };

  SVMINT i;
  if(iteration >= max_iterations){
    cout<<"*** No convergence: Time up."<<endl;
    if(examples_total<examples->size()){
      // set sums for all variables for statistics
      reset_shrinked();
      SVMFLOAT new_convergence_epsilon = 0;
      SVMFLOAT the_lambda;
      for(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;
    };
  };

  time_all = get_time() - time_all;

  // calculate b
  SVMFLOAT new_b=0;
  SVMINT new_b_count=0;
  for(i=0;i<examples_total;i++){
    if((all_alphas[i]-Cneg < -is_zero) && (all_alphas[i]>is_zero)){
      new_b +=  all_ys[i] - sum[i]-epsilon_neg;
      new_b_count++;
    }
    else if((all_alphas[i]+Cpos > is_zero) && (all_alphas[i]<-is_zero)){
      new_b +=  all_ys[i] - sum[i]+epsilon_pos;
      new_b_count++;
    };
    examples->put_alpha(i,all_alphas[i]);
  };

  if(new_b_count>0){
    examples->put_b(new_b/((SVMFLOAT)new_b_count));
  }
  else{
    // unlikely
    for(i=0;i<examples_total;i++){
      if((all_alphas[i]<is_zero) && (all_alphas[i]>-is_zero)) {
	new_b +=  all_ys[i]- sum[i];
	new_b_count++;
      };
    };
    if(new_b_count>0){
      examples->put_b(new_b/((SVMFLOAT)new_b_count));
    }
    else{
      // even unlikelier
      for(i=0;i<examples_total;i++){
	new_b +=  all_ys[i]- sum[i];
	new_b_count++;
      };
      examples->put_b(new_b/((SVMFLOAT)new_b_count));
    };
  };


  if(parameters->verbosity>= 2){
    cout<<"Done training: "<<iteration<<" iterations."<<endl;
    if(parameters->verbosity>= 2){
      // cout<<"lambda_eq = "<<lambda_eq<<endl;
      SVMFLOAT now_target=0;
      SVMFLOAT now_target_dummy=0;
      for(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;
    };
  };

  the_result = print_statistics();