svm_c.cpp

支持向量机(SVM)的VC源代码

#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++){