svm_nu.cpp

这是我找到的C语言编写的支持向量机程序包

    //      target_count=0;
    if(parameters->descend < old_target - new_target){
      target_count=0;
    }
    else{
      convError=1;
    };
    if(parameters->verbosity>=5){
      cout<<"descend = "<<old_target-new_target<<"  ("<<old_target<<" --> "<<new_target<<")"<<endl;
    };
  }
  else{
    // nothing we can do
    KKTerror=0;
    convError=1;
    if(parameters->verbosity>=5){
      cout<<"WARNING: no descend ("<<old_target<<" -> "<<new_target<<"), stopping"<<endl;
    };
  };

  if(convError){
    target_count++;
    if(new_target>=old_target){
      for(i=0;i<working_set_size;i++){
	primal[i] = qp.A[i]*all_alphas[working_set[i]];
      };  
    };                                                                        
    if(parameters->verbosity>=5){	
      cout<<"WARNING: Convergence error, setting sigfig = "<<sigfig_max<<endl;
    };
  };

  if(target_count>50){
    // non-recoverable numerical error
    feasible_epsilon=1;
    convergence_epsilon*=2;
    if(parameters->verbosity>=1)
      cout<<"WARNING: reducing KKT precision to "<<convergence_epsilon<<endl;
    target_count=0;
  };

  if(parameters->verbosity>=5){	
    cout<<"Resulting values:"<<endl;
    for(i=0;i<working_set_size;i++){
      cout<<i<<": "<<primal[i]<<endl;
    };
  };

  time_optimize += get_time() - time_start;
};


void svm_nu_regression_c::print_special_statistics(){
  // calculate tube size epsilon
  SVMFLOAT b = examples->get_b();
  SVMFLOAT epsilon_pos = 0;
  SVMFLOAT epsilon_neg = 0;
  SVMINT pos_count = 0;
  SVMINT neg_count = 0;
  SVMINT i;
  for(i=0;i<examples_total;i++){
    if((all_alphas[i] > is_zero) && (all_alphas[i]-Cpos<-is_zero)){
      epsilon_neg += all_ys[i]-sum[i]-b;
      neg_count++;
    }
    else if((all_alphas[i] <- is_zero) && (all_alphas[i]+Cneg>+is_zero)){
      epsilon_pos += -all_ys[i]+sum[i]+b;
      pos_count++;
    };
  };
  if((parameters->Lpos == parameters->Lneg) ||
     (pos_count == 0) ||
     (neg_count == 0)){
    // symmetrical
    epsilon_pos += epsilon_neg;
    pos_count += neg_count;
    if(pos_count>0){
      epsilon_pos /= (SVMINT)pos_count;
      cout<<"epsilon = "<<epsilon_pos<<endl;
    }
    else{
      cout<<"ERROR: could not calculate epsilon."<<endl;
      cout<<pos_count<<"\t"<<neg_count<<endl; // @@@@@@@
    };
  }
  else{
    // asymmetrical
    epsilon_pos /= (SVMINT)pos_count;
    cout<<"epsilon+ = "<<epsilon_pos<<endl;
    epsilon_neg /= (SVMINT)neg_count;
    cout<<"epsilon- = "<<epsilon_pos<<endl;
  };
};


/**
 *
 * svm_nu_pattern_c
 *
 **/

SVMFLOAT svm_nu_pattern_c::nabla(const SVMINT i){
  if(all_ys[i] > 0){
    return( sum[i]);
  }
  else{
    return(-sum[i]);
  };
};


void svm_nu_pattern_c::init(kernel_c* new_kernel, parameters_c* new_parameters){
  new_parameters->realC = 1;
  svm_nu_regression_c::init(new_kernel,new_parameters);
};


void svm_nu_pattern_c::init_optimizer(){
  // Cs are dived by examples_total in init_optimizer
  svm_nu_regression_c::init_optimizer();
  SVMINT i;
  for(i=0;i<working_set_size;i++){
    qp.l[i] = 0;
  };
};


void svm_nu_pattern_c::update_working_set(){
  svm_c::update_working_set();
  SVMINT i;
  for(i=0;i<working_set_size;i++){
    if(qp.A[i]>0){
      qp.c[i] += all_ys[working_set[i]];
    }
    else{
      qp.c[i] -= all_ys[working_set[i]];
    };
  };
};


void svm_nu_pattern_c::init_working_set(){
  // calculate nu-sum 

  if(examples->initialised_alpha()){
    project_to_constraint();
  };

  sum_alpha_nu=0;
  SVMFLOAT the_nu_sum = 0;
  SVMFLOAT the_sum=0;
  SVMINT pos_count=0;
  SVMINT neg_count=0;
  SVMINT ni;
  for(ni=0;ni<examples_total;ni++){
    the_sum += all_alphas[ni];
    the_nu_sum += abs(all_alphas[ni]);
    if(is_alpha_neg(ni)> 0){
      neg_count++;
    }
    else{
      pos_count++;
    };
  };

  if((abs(the_sum) > is_zero) || (abs(the_nu_sum-nu) > is_zero)){
    // set initial feasible point
    // neg alpha: -nu/2n
    // pos alpha:  nu/2p

    if((nu*(SVMFLOAT)examples_total>2*(SVMFLOAT)pos_count) ||
       (nu*(SVMFLOAT)examples_total>2*(SVMFLOAT)neg_count)){
      nu = 2*((SVMFLOAT)pos_count)/((SVMFLOAT)examples_total);
      if(nu > 2*((SVMFLOAT)neg_count)/((SVMFLOAT)examples_total)){
	nu = 2*((SVMFLOAT)neg_count)/((SVMFLOAT)examples_total);
      };
      nu -= is_zero; // just to make sure
      cout<<"ERROR: nu too large, setting nu = "<<nu<<endl;
    };

    for(ni=0;ni<examples_total;ni++){
      if(is_alpha_neg(ni)> 0){
	examples->put_alpha(ni,nu/(2*(SVMFLOAT)neg_count));
      }
      else{
	examples->put_alpha(ni,-nu/(2*(SVMFLOAT)pos_count));
      };
    };
    examples->set_initialised_alpha();
  };

  svm_c::init_working_set();
};


void svm_nu_pattern_c::print_special_statistics(){
  // calculate margin rho
  SVMFLOAT b = examples->get_b();
  SVMFLOAT rho_pos = 0;
  SVMFLOAT rho_neg = 0;
  SVMINT pos_count = 0;
  SVMINT neg_count = 0;
  SVMINT i;
  for(i=0;i<examples_total;i++){
    if((all_alphas[i] > is_zero) && (all_alphas[i]-Cpos<-is_zero)){
      rho_neg += sum[i]+b;
      neg_count++;
    }
    else if((all_alphas[i] <- is_zero) && (all_alphas[i]+Cneg>+is_zero)){
      rho_pos += -sum[i]-b;
      pos_count++;
    };
  };
  if((parameters->Lpos == parameters->Lneg) ||
     (pos_count == 0) ||
     (neg_count == 0)){
    // symmetrical
    rho_pos += rho_neg;
    pos_count += neg_count;
    if(pos_count>0){
      rho_pos /= (SVMINT)pos_count;
      cout<<"margin = "<<rho_pos<<endl;
    }
    else{
      cout<<"ERROR: could not calculate margin."<<endl;
    };
  }
  else{
    // asymmetrical
    rho_pos /= (SVMINT)pos_count;    cout<<"margin+ = "<<rho_pos<<endl;
    rho_neg /= (SVMINT)neg_count;
    cout<<"margin- = "<<rho_pos<<endl;
  };
};


/**
 *
 * svm_distribution_c
 *
 **/

int svm_distribution_c::is_alpha_neg(const SVMINT i){
  // variable i is alpha*
  return 1;
};


SVMFLOAT svm_distribution_c::nabla(const SVMINT i){
  return( sum[i]);
};


SVMFLOAT svm_distribution_c::lambda(const SVMINT i){
  // size lagrangian multiplier of the active constraint

  SVMFLOAT alpha;
  SVMFLOAT result = 0;

  alpha=all_alphas[i];

  if(alpha>is_zero){
    // alpha*
    if(alpha-Cneg >= - is_zero){
      // upper bound active
      result = -lambda_eq-sum[i];
    }
    else{
      result = -abs(sum[i]+lambda_eq);
    };
  }
  else{
    // lower bound active
    result = sum[i] + lambda_eq;
  };

  return result;
};


int svm_distribution_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 - sum[i];
    if(result>=-feasible_epsilon){
      return 0; 
    };
  }
  else if(alpha<=is_zero){
    // lower bound active
    result = sum[i]+lambda_eq;
    if(result>=-feasible_epsilon){
      return 0; 
    };
  }
  else{
    // not at bound
    result= abs(sum[i]+lambda_eq);
    if(result<=feasible_epsilon){
      return 0; 
    };
  };
  return 1;
};


int svm_distribution_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;

  alpha=all_alphas[i];
  *the_nabla = sum[i];

  if(alpha >= Cneg){ //alpha-Cneg >= - is_zero){
    // alpha* at upper bound
    *atbound = 1;
    *the_lambda = -lambda_eq - *the_nabla; //sum[i] + 1;
    if(*the_lambda >= 0){
      at_bound[i]++;
      if(at_bound[i] == shrink_const) to_shrink++;
    }
    else{
      at_bound[i] = 0;
    };
  }
  else if(alpha <= 0){
    // lower bound active
    *atbound = -1;
    *the_lambda = lambda_eq + *the_nabla; //sum[i] + 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;
    *the_lambda = -abs(*the_nabla+lambda_eq);
    at_bound[i] = 0;
  };
  if(*the_lambda >= feasible_epsilon){
    is_feasible = 0; 
  };

  return is_feasible;
};


void svm_distribution_c::init(kernel_c* new_kernel, parameters_c* new_parameters){
  new_parameters->realC = 1;
  nu = new_parameters->nu;
  convergence_epsilon = 1e-4;
  svm_pattern_c::init(new_kernel,new_parameters);
  //  is_pattern = 1;
};

void svm_distribution_c::init_optimizer(){
  // Cs are dived by examples_total in init_optimizer
  svm_pattern_c::init_optimizer();
};


void svm_distribution_c::project_to_constraint(){
  SVMINT total = 0;
  SVMFLOAT alpha_sum=sum_alpha-nu;
  SVMFLOAT alpha=0;
  SVMINT i;
  for(i=0;i<examples_total;i++){
    alpha = all_alphas[i];
    alpha_sum += alpha;
    if((alpha>is_zero) && (alpha-Cneg < -is_zero)){
      total++;
    };
  };
  if(total>0){
    // equality constraint violated
    alpha_sum /= (SVMFLOAT)total;
    for(i=0;i<examples_total;i++){
      if((alpha>is_zero) && (alpha-Cneg < -is_zero)){
	all_alphas[i] -= alpha_sum;
      };
    };
  };
};


int svm_distribution_c::convergence(){
  long time_start = get_time();
  SVMFLOAT the_lambda_eq = 0;
  SVMINT total = 0;
  SVMFLOAT alpha_sum=0;
  SVMFLOAT alpha=0;
  SVMINT i;
  int result=1;

  // actual convergence-test
  total = 0; alpha_sum=0;
  //  cout<<Cneg<<"\t"<<nu<<"\t"<<all_alphas[0]<<endl;
  for(i=0;i<examples_total;i++){
    alpha = all_alphas[i];
    alpha_sum += alpha;
    if((alpha>is_zero) && (alpha-Cneg < -is_zero)){
      // alpha^* = - nabla
      the_lambda_eq += -sum[i];
      total++;
    };
  };

  if(parameters->verbosity>= 4){
    cout<<"lambda_eq = "<<(the_lambda_eq/total)<<endl;
  };
  if(total>0){
    lambda_eq = the_lambda_eq / total;
  }
  else{
    // keep WS lambda_eq
    lambda_eq = lambda_WS;
    if(parameters->verbosity>= 4){
      cout<<"*** no SVs in convergence(), lambda_eq = "<<lambda_eq<<"."<<endl;
    };
  };

  if(target_count>2){
    if(target_count>20){
      // desperate!
      lambda_eq = ((40-target_count)*lambda_eq + (target_count-20)*lambda_WS)/20;
      if(parameters->verbosity>=5){
	cout<<"Re-Re-calculated lambda from WS: "<<lambda_eq<<endl;
      };
      if(target_count>40){
	// really desperate, kick one example out!
	i = working_set[target_count%working_set_size];
	lambda_eq = -sum[i];
	if(parameters->verbosity>=5){
	  cout<<"set lambda_eq to nabla("<<i<<"): "<<lambda_eq<<endl;
	};
      };
    }
    else{
      lambda_eq = lambda_WS;
      if(parameters->verbosity>=5){
	cout<<"Re-calculated lambda_eq from WS: "<<lambda_eq<<endl;
      };
    };
  };

  // check linear constraint
  if(abs(alpha_sum+sum_alpha-nu) > convergence_epsilon){
    // equality constraint violated
    if(parameters->verbosity>= 4){
      cout<<"No convergence: equality constraint violated: |"<<(alpha_sum+sum_alpha)<<"| >> 0"<<endl;
    };
    project_to_constraint();
    result = 0;  
  };

  i=0;
  while((i<examples_total) && (result != 0)){
    if(lambda(i)>=-convergence_epsilon){
      i++;
    }
    else{
      result = 0;
    };
  };

  time_convergence += get_time() - time_start;
  return result;
};


void svm_distribution_c::init_working_set(){
  // calculate sum
  SVMINT i,j;

  if(nu>1){
    cout<<"ERROR: nu too large, setting nu to 1"<<endl;
    nu = 1-is_zero;
  };
  SVMFLOAT the_sum=0;
  for(i=0; i<examples_total;i++){
    the_sum += all_alphas[i];
  };
  if(abs(the_sum-nu) > is_zero){
    for(i=0; i<examples_total;i++){
      examples->put_alpha(i,nu/((SVMFLOAT)examples_total));
    };
    examples->set_initialised_alpha();
  };

  if(parameters->verbosity >= 3){
    cout<<"Initialising variables, this may take some time."<<endl;
  };
  for(i=0; i<examples_total;i++){
    all_ys[i] = 1;
    sum[i] = 0;
    at_bound[i] = 0;
    for(j=0; j<examples_total;j++){
      sum[i] += all_alphas[j]*kernel->calculate_K(i,j);
    };
  };

  calculate_working_set();
  update_working_set();
};


void svm_distribution_c::print_special_statistics(){
  // calculate margin rho
  SVMFLOAT rho = 0;
  SVMINT count = 0;
  SVMFLOAT norm_x;
  SVMFLOAT max_norm_x=-infinity;
  //  SVMFLOAT xi_i;
  //  SVMINT estim_loo=examples_total;
  //  SVMINT estim_loo2=examples_total;
  SVMINT svs=0;
  SVMINT i;
  for(i=0;i<examples_total;i++){
    if((all_alphas[i] > is_zero) && (all_alphas[i]-Cpos<-is_zero)){
      rho += sum[i];
      count++;
    };
    if(all_alphas[i] != 0){
      svs++;
      norm_x = kernel->calculate_K(i,i);
      if(norm_x>max_norm_x){
	max_norm_x = norm_x;
      };
    };
  };
  if(count == 0){
    cout<<"ERROR: could not calculate margin."<<endl;
  }
  else{
    // put -rho as b (same decision function)
    rho /= (SVMINT)count;
    examples->put_b(-rho);
    cout<<"margin = "<<rho<<endl;
  };

  cout<<"examples in distribution support : "<<count<<" ("<<((SVMINT)(10000.0*(SVMFLOAT)count/((SVMFLOAT)examples_total)))/100.0<<"%)."<<endl;


};