svm_c.cpp

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

  if(parameters->loo_estim){
    r_delta = 0;
    SVMFLOAT r_current;
    for(SVMINT j=0;j<examples_total;j++){
      norm_x = kernel->calculate_K(j,j);
      for(i=0;i<examples_total;i++){
	r_current = norm_x-kernel->calculate_K(i,j);
	if(r_current > r_delta){
	  r_delta = r_current;
	};
      };
    };
  };

  for(i=0;i<examples_total;i++){
    alpha=all_alphas[i];
    if(alpha<min_alpha) min_alpha = alpha;
    if(alpha>max_alpha) max_alpha = alpha;
    prediction = predict(i);
    y = examples->unscale_y(all_ys[i]);
    actloss=loss(prediction,y);
    theloss+=actloss;
    MAE += abs(prediction-y);
    MSE += (prediction-y)*(prediction-y);
    if(y < prediction-parameters->epsilon_pos){
      theloss_pos += actloss;
      countpos++;
    }
    else if(y > prediction+parameters->epsilon_neg){
      theloss_neg += actloss;
      countneg++;
    };
    if(parameters->loo_estim){
      if(abs(alpha)>is_zero){ 
	if(is_alpha_neg(i)>=0){
	  loo_loss_estim += loss(prediction-(abs(alpha)*(2*kernel->calculate_K(i,i)+r_delta)+2*epsilon_neg),y);
	}
	else{
	  loo_loss_estim += loss(prediction+(abs(alpha)*(2*kernel->calculate_K(i,i)+r_delta)+2*epsilon_pos),y);
	};
      };
    }
    else{
      // loss doesn't change if non-SV is omitted
      loo_loss_estim += actloss;
    };

    if(abs(alpha)>is_zero){ 
     // a support vector
      svs++; 
      if((alpha-Cneg >= -is_zero) || (alpha+Cpos <= is_zero)){ 
	bsv++; 
      };
    };

    if(is_pattern){
      if(y>0){
	if(prediction>0){
	  correct_pos++;
	};
	if(prediction>1){
	  xi=0;
	}
	else{
	  xi=1-prediction;
	};
	if(2*alpha*r_delta+xi >= 1){
	  estim_pos++;
	};
	total_pos++;
      }
      else{
	if(prediction<=0){
	  correct_neg++;
	};
	if(prediction<-1){
	  xi=0;
	}
	else{
	  xi=1+prediction;
	};
	if(2*(-alpha)*r_delta+xi >= 1){
	  estim_neg++;
	};
	total_neg++;
      };
    };    
  };
  if(countpos != 0){
    theloss_pos /= (SVMFLOAT)countpos;
  };
  if(countneg != 0){
    theloss_neg /= (SVMFLOAT)countneg;
  };

  the_result.MAE = MAE / (SVMFLOAT)examples_total;
  the_result.MSE = MSE / (SVMFLOAT)examples_total;
  the_result.VCdim = 1+norm_w*max_norm_x;
  the_result.loss = theloss/((SVMFLOAT)examples_total);
  if(parameters->loo_estim){
    the_result.pred_loss = loo_loss_estim/((SVMFLOAT)examples_total);
  }
  else{
    the_result.pred_loss = the_result.loss;
  };
  the_result.loss_pos = theloss_pos;
  the_result.loss_neg = theloss_neg;
  the_result.number_svs = svs;
  the_result.number_bsv = bsv;
  if(is_pattern){
    the_result.accuracy = ((SVMFLOAT)(correct_pos+correct_neg))/((SVMFLOAT)(total_pos+total_neg));
    the_result.precision = ((SVMFLOAT)correct_pos/((SVMFLOAT)(correct_pos+total_neg-correct_neg)));
    the_result.recall = ((SVMFLOAT)correct_pos/(SVMFLOAT)total_pos);
    if(parameters->loo_estim){
      the_result.pred_accuracy = (1-((SVMFLOAT)(estim_pos+estim_neg))/((SVMFLOAT)(total_pos+total_neg)));
      the_result.pred_precision = ((SVMFLOAT)(total_pos-estim_pos))/((SVMFLOAT)(total_pos-estim_pos+estim_neg));
      the_result.pred_recall = (1-(SVMFLOAT)estim_pos/((SVMFLOAT)total_pos));
    }
    else{
      the_result.pred_accuracy = the_result.accuracy;
      the_result.pred_precision = the_result.precision;
      the_result.pred_recall = the_result.recall;
    };
  }
  else{
    the_result.accuracy = -1;
    the_result.precision = -1;
    the_result.recall = -1;
    the_result.pred_accuracy = -1;
    the_result.pred_precision = -1;
    the_result.pred_recall = -1;
  };


  if(convergence_epsilon > parameters->convergence_epsilon){
    cout<<"WARNING: The results were obtained using a relaxed epsilon of "<<convergence_epsilon<<" on the KKT conditions!"<<endl;
  }
  else if(parameters->verbosity>=2){
    cout<<"The results are valid with an epsilon of "<<convergence_epsilon<<" on the KKT conditions."<<endl;
  };
  if(parameters->verbosity >= 2){
    cout << "Average loss  : "<<the_result.loss<<" (loo-estim: "<< the_result.pred_loss<<")"<<endl;
    cout << "Avg. loss pos : "<<theloss_pos<<"\t ("<<countpos<<" occurences)"<<endl;
    cout << "Avg. loss neg : "<<theloss_neg<<"\t ("<<countneg<<" occurences)"<<endl;
    cout << "Mean absolute error : "<<the_result.MAE<<endl;
    cout << "Mean squared error  : "<<the_result.MSE<<endl;
    cout << "Support Vectors : "<<svs<<endl;
    cout << "Bounded SVs     : "<<bsv<<endl;
    cout<<"min SV: "<<min_alpha<<endl
	<<"max SV: "<<max_alpha<<endl;
    cout<<"|w| = "<<sqrt(norm_w)<<endl;
    cout<<"max |x| = "<<sqrt(max_norm_x)<<endl;
    cout<<"VCdim <= "<<the_result.VCdim<<endl;

    print_special_statistics();

    if((is_pattern) && (! parameters->is_distribution)){
      // output precision, recall and accuracy
      if(parameters->loo_estim){
	cout<<"performance (+estimators):"<<endl;
	cout<<"Accuracy  : "<<the_result.accuracy<<" ("<<the_result.pred_accuracy<<")"<<endl;
	cout<<"Precision : "<<the_result.precision<<" ("<<the_result.pred_precision<<")"<<endl;
	cout<<"Recall    : "<<the_result.recall<<" ("<<the_result.pred_recall<<")"<<endl;
      }
      else{
	cout<<"performance :"<<endl;
	cout<<"Accuracy  : "<<the_result.accuracy<<endl;
	cout<<"Precision : "<<the_result.precision<<endl;
	cout<<"Recall    : "<<the_result.recall<<endl;
      };
      if(parameters->verbosity>= 2){
	// nice printout ;-)
	int rows = (int)(1+log10((SVMFLOAT)(total_pos+total_neg)));
	int now_digits = rows+2;
	int i,j;
	cout<<endl;
	cout<<"Predicted values:"<<endl;
	cout<<"   |";
	for(i=0;i<rows;i++){ cout<<" "; };
	cout<<"+  |";
	for(j=0;j<rows;j++){ cout<<" "; };
	cout<<"-"<<endl;

	cout<<"---+";
	for(i=0;i<now_digits;i++){ cout<<"-"; };
	cout<<"-+-";
	for(i=0;i<now_digits;i++){ cout<<"-"; };
	cout<<endl;

	cout<<" + |  ";
	now_digits=rows-(int)(1+log10((SVMFLOAT)correct_pos))-1;
	for(i=0;i<now_digits;i++){ cout<<" "; };
	cout<<correct_pos<<"  |  ";
	now_digits=rows-(int)(1+log10((SVMFLOAT)(total_pos-correct_pos)))-1;
	for(i=0;i<now_digits;i++){ cout<<" "; };
	cout<<total_pos-correct_pos<<"    (true pos)"<<endl;

	cout<<" - |  ";
	now_digits=rows-(int)(1+log10((SVMFLOAT)(total_neg-correct_neg)))-1;
	for(i=0;i<now_digits;i++){ cout<<" "; };
	cout<<(total_neg-correct_neg)<<"  |  ";
	now_digits=rows-(int)(1+log10((SVMFLOAT)correct_neg))-1;
	for(i=0;i<now_digits;i++){ cout<<" "; };
	cout<<correct_neg<<"    (true neg)"<<endl;
	cout<<endl;
      };
    };
  };

  SVMINT dim = examples->get_dim();
  if(((parameters->print_w == 1) && (parameters->is_linear != 0)) ||
     ((dim<100) && 
      (parameters->verbosity>= 2) && 
      (parameters->is_linear != 0))
     ){
    // print hyperplane
    SVMINT j;
    svm_example example;
    SVMFLOAT* w = new SVMFLOAT[dim];
    SVMFLOAT b = examples->get_b();
    for(j=0;j<dim;j++) w[j] = 0;
    for(i=0;i<examples_total;i++){
      example = examples->get_example(i);
      alpha = examples->get_alpha(i);
      for(j=0;j<example.length;j++){
	w[((example.example)[j]).index] += alpha*((example.example)[j]).att;
      };
    };
    if(examples->initialised_scale()){
      SVMFLOAT* exp = examples->get_exp();
      SVMFLOAT* var = examples->get_var();
      for(j=0;j<dim;j++){
	if(var[j] != 0){
	  w[j] /= var[j];
	};
	if(0 != var[dim]){
	  w[j] *= var[dim];
	};
	b -= w[j]*exp[j];
      };
      b += exp[dim];
    };
    for(j=0;j<dim;j++){
      cout << "w["<<j<<"] = " << w[j] << endl;
    };
    cout << "b = "<<b<<endl;
    if(dim==1){
      cout<<"y = "<<w[0]<<"*x+"<<b<<endl;
    };
    if((dim==2) && (is_pattern)){
      cout<<"x1 = "<<-w[0]/w[1]<<"*x0+"<<-b/w[1]<<endl;
    };
    delete []w;
  };

  if(parameters->verbosity>= 2){
    cout<<"Time for learning:"<<endl
	<<"init        : "<<(time_init/100)<<"s"<<endl
	<<"optimizer   : "<<(time_optimize/100)<<"s"<<endl
	<<"convergence : "<<(time_convergence/100)<<"s"<<endl
	<<"update ws   : "<<(time_update/100)<<"s"<<endl
	<<"calc ws     : "<<(time_calc/100)<<"s"<<endl
	<<"============="<<endl
	<<"all         : "<<(time_all/100)<<"s"<<endl;
  }
  else if(parameters->verbosity>=2){
    cout<<"Time for learning: "<<(time_all/100)<<"s"<<endl;
  };

  return the_result;
};


/**
 *
 * Pattern SVM
 *
 */

int svm_pattern_c::is_alpha_neg(const SVMINT i){
  // variable i is alpha*

  int result;
  if(all_ys[i] > 0){
    result = 1;
  }
  else{
    result = -1;
  };
  return result;
};


SVMFLOAT svm_pattern_c::nabla(const SVMINT i){
  // = is_alpha_neg(i) * sum[i] - 1
  if(all_ys[i]>0){
    return(sum[i]-1);
  }
  else{
    return(-sum[i]-1);
  };
};


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

  SVMFLOAT alpha;
  SVMFLOAT result=0;

  alpha=all_alphas[i];

  if(alpha == Cneg){ //-Cneg >= - is_zero){
    // upper bound active
    result = -lambda_eq - sum[i] + 1;
  }
  else if(alpha == 0){ //(alpha <= is_zero) && (alpha >= -is_zero)){
    // lower bound active
    if(all_ys[i]>0){
      result = (sum[i]+lambda_eq) - 1;
    }
    else{
      result = -(sum[i]+lambda_eq) - 1;
    };
  }
  else if(alpha == -Cpos){ //+Cpos <= is_zero){
    // upper bound active
    result = lambda_eq + sum[i] + 1;
  }
  else{
    if(all_ys[i]>0){
      result = -abs(sum[i]+lambda_eq - 1);
    }
    else{
      result = -abs(-sum[i]-lambda_eq - 1);
    };
  };

  return result;
};


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

  return is_feasible;
};


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


/**
 *
 * Regression SVM
 *
 */

int svm_regression_c::is_alpha_neg(const SVMINT i){
  // variable i is alpha*
  int result;
  if(all_alphas[i] > 0){
    result = 1;
  }
  else if(all_alphas[i] == 0){
    if(sum[i] - all_ys[i] + lambda_eq>0){
      result = -1;
    }
    else{
      result = 1;
    };
    //    result = 2*((i+at_bound[i])%2)-1;
  }
  else{
    result = -1;
  };
  return result;
};


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


SVMFLOAT svm_regression_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-sum[i] + all_ys[i] - epsilon_neg;
    };
  }
  else if(alpha >= -is_zero){
    // lower bound active
    if(all_alphas[i] > 0){
      result = sum[i] - all_ys[i] + epsilon_neg + lambda_eq;
    }
    else if(all_alphas[i] == 0){
      if(is_alpha_neg(i)>0){
	result = sum[i] - all_ys[i] + epsilon_neg + lambda_eq;
      }
      else{
	result = -sum[i] + all_ys[i] + epsilon_pos-lambda_eq;
      };
    }
    else{
      result = -sum[i] + all_ys[i] + epsilon_pos-lambda_eq;
    };
  }
  else if(alpha+Cpos <= is_zero){
    // upper bound active
    result = lambda_eq + sum[i] - all_ys[i] - epsilon_pos;
  };

  return result;
};


int svm_regression_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; // <=> at_bound < shrink and lambda < -feas_eps

  if(at_bound[i] >= shrink_const){ is_feasible = 0; };

  SVMFLOAT alpha;

  alpha=all_alphas[i];

  if(alpha-Cneg >= -is_zero){
    // alpha* at upper bound
    *atbound = 1;
    *the_nabla = sum[i] - all_ys[i] + epsilon_neg;
    *the_lambda = -lambda_eq- *the_nabla;
    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(all_alphas[i] > 0){
      *the_nabla = sum[i] - all_ys[i] + epsilon_neg;
      *the_lambda = *the_nabla + lambda_eq;
    }
    else if(all_alphas[i]==0){
      if(is_alpha_neg(i)>0){
	*the_nabla = sum[i] - all_ys[i] + epsilon_neg;
	*the_lambda = *the_nabla + lambda_eq;
      }
      else{
	*the_nabla = -sum[i] + all_ys[i] + epsilon_pos;
	*the_lambda = *the_nabla - lambda_eq;
      };
    }
    else{
      *the_nabla = -sum[i] + all_ys[i] + epsilon_pos;
      *the_lambda = *the_nabla - lambda_eq;
    };
    if(*the_lambda >= 0){
      if(all_alphas[i] != 0){
	// check both constraints!
	all_alphas[i]=0;
	if(is_alpha_neg(i)>0){
	  *the_lambda = *the_nabla + lambda_eq;
	}
	else{
	  *the_lambda = *the_nabla - lambda_eq;
	};
	if(*the_lambda >= 0){
	  at_bound[i]++;
	};
      }
      else{
	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;
    *the_nabla = -sum[i] + all_ys[i] + epsilon_pos;
    *the_lambda = lambda_eq - *the_nabla;
    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;
    if(is_alpha_neg(i)>0){
      *the_nabla = sum[i] - all_ys[i] + epsilon_neg;
      *the_lambda = -abs(*the_nabla + lambda_eq);
    }
    else{
      *the_nabla = -sum[i] + all_ys[i] + epsilon_pos;
      *the_lambda = -abs(lambda_eq - *the_nabla);
    };
  };
  if(*the_lambda>=feasible_epsilon){
    is_feasible = 0; 
  };
  return is_feasible;
};