svmpattern.java

一个java程序编写的svm支持向量机小程序

package edu.udo.cs.mySVMdb.SVM;


public class SVMpattern extends SVM
{
    /**
     * Class for pattern recognition SVM
     * @author Stefan R黳ing
     * @version 1.0
     */
  
    public SVMpattern()
    {
    };


    /**
     * Calls the optimizer
     */
    protected void optimize()
	throws Exception
    {
	// optimizer-specific call

	//	qp.n = working_set_size;
	
	int i;
	int j;
	
	double[] my_primal = primal;

	// equality constraint
	qp.b[0]=0;
	for(i=0;i<working_set_size;i++){
	    qp.b[0] += the_container.get_alpha(working_set[i]);
	};
	
	// set initial optimization parameters
	double new_target=0;
	double old_target=0;
	double target_tmp;
	for(i=0;i<working_set_size;i++){
	    target_tmp = my_primal[i]*qp.H[i*working_set_size+i]/2;
	    for(j=0;j<i;j++){
		target_tmp+=my_primal[j]*qp.H[j*working_set_size+i];
	    };
	    target_tmp+=qp.c[i];
	    old_target+=target_tmp*my_primal[i];
	};
	
	double new_constraint_sum=0;
	double my_is_zero = is_zero;
	int sv_count=working_set_size;
	
	// optimize
	boolean KKTerror=true; // KKT not yet satisfied
	boolean convError=false; // KKT can still be satisfied
	
	qp.max_allowed_error = convergence_epsilon;

        qp.x = primal;
        qp.solve();
        primal = qp.x;
        lambda_WS = qp.lambda_eq;
	my_primal = primal;

	// loop while some KKT condition is not valid (alpha=0)
	int it=3;
	while(KKTerror){
	    // clip
	    sv_count=working_set_size;
	    new_constraint_sum=qp.b[0];
	    for(i=0;i<working_set_size;i++){
		// check if at bound
		if(my_primal[i] <= my_is_zero){
		    // at lower bound
		    my_primal[i] = qp.l[i];
		    sv_count--;
		}
		else if(qp.u[i]-my_primal[i] <= my_is_zero){
		    // at upper bound
		    my_primal[i] = qp.u[i];
		    sv_count--;
		};
		new_constraint_sum -= qp.A[i]*my_primal[i];
	    };
	    
	    // enforce equality constraint
	    if(sv_count>0){
		new_constraint_sum /= (double)sv_count;
		logln(5,"adjusting "+sv_count+" alphas by "+new_constraint_sum);
		for(i=0;i<working_set_size;i++){
		    if((my_primal[i] > qp.l[i]) && 
		       (my_primal[i] < qp.u[i])){
			// real sv
			my_primal[i] += qp.A[i]*new_constraint_sum;
		    };
		};
	    }
	    else if(Math.abs(new_constraint_sum)>(double)working_set_size*is_zero){
		// error, can't get feasible point
		logln(5,"WARNING: No SVs, constraint_sum = "+new_constraint_sum);
		old_target = -Double.MIN_VALUE; 
		convError=true;
	    };
	    // test descend
	    new_target=0;
	    for(i=0;i<working_set_size;i++){
		// attention: optimizer changes one triangle of H!
		target_tmp = my_primal[i]*qp.H[i*working_set_size+i]/2.0;
		for(j=0;j<i;j++){
		    target_tmp+=my_primal[j]*qp.H[j*working_set_size+i];
		};
		target_tmp+=qp.c[i];
		new_target+=target_tmp*my_primal[i];
	    };
	    
	    if(new_target < old_target){
		KKTerror = false;
		if(descend < old_target - new_target){
		    target_count=0;
		}
		else{
		    convError=true;
		};
		logln(5,"descend = "+(old_target-new_target));
	    }
	    else if(sv_count > 0){
		// less SVs
		// set my_is_zero to min_i(primal[i]-qp.l[i], qp.u[i]-primal[i])
		my_is_zero = Double.MAX_VALUE;
		for(i=0;i<working_set_size;i++){
		    if((my_primal[i] > qp.l[i]) && (my_primal[i] < qp.u[i])){
			if(my_primal[i] - qp.l[i] < my_is_zero){
			    my_is_zero = my_primal[i]-qp.l[i];
			};
			if(qp.u[i]  - my_primal[i]  < my_is_zero){
			    my_is_zero = qp.u[i] - my_primal[i];
			};
		    };
		};
		if(target_count == 0){
		    my_is_zero *= 2;
		};
		logln(5,"WARNING: no descend ("+(old_target-new_target)
		      +" <= "+descend
		      +"), adjusting is_zero to "+my_is_zero);
		logln(5,"new_target = "+new_target);
	    }
	    else{
		// nothing we can do
		logln(5,"WARNING: no descend ("+(old_target-new_target)
		      +" <= "+descend+"), stopping.");
		KKTerror=false;
		convError=true;
	    };
	};
	
	if(convError){
	    target_count++;
	    if(old_target < new_target){
		for(i=0;i<working_set_size;i++){
		    my_primal[i] = qp.A[i]*the_container.get_alpha(working_set[i]);
		};                              
		logln(5,"WARNING: Convergence error, restoring old primals");
	    };                                          
	};
	
	if(target_count>50){
	    // non-recoverable numerical error
	    convergence_epsilon*=2;
	    feasible_epsilon = convergence_epsilon;
	    logln(1,"WARNING: reducing KKT precision to "+convergence_epsilon);
	    target_count=0;
	};
    };
    


    protected final boolean is_alpha_neg(int i)
    {
	boolean result;
	if(the_container.get_y(i) > 0){
	    result = true;
	}
	else{
	    result = false;
	};
	return result;
    };


    protected final double nabla(int i)
    {
	double result;
      if(is_alpha_neg(i)){
	result = ( sum[i] - 1);
      }
      else{
	result = (-sum[i] -1);
      };
      return result;
    };


    protected void print_statistics()
	throws Exception
    {
      int dim = the_container.get_dim();
      int i,j;
      double alpha;
      double[] x;
      int svs=0;
      int bsv = 0;
      int correct_pos=0;
      int correct_neg=0;
      int total_pos=0;
      int total_neg=0;
      double y;
      double prediction;
      double min_lambda = Double.MAX_VALUE;
      double b = the_container.get_b();
      for(i=0;i<examples_total;i++){
	  if(lambda(i) < min_lambda){
	      min_lambda = lambda(i);
	  };
 	y = the_container.get_y(i);
 	prediction = sum[i]+b;
 	alpha = the_container.get_alpha(i);
	if(y>0){
	    if(prediction>0){
		correct_pos++;
	    };
	    total_pos++;
	}
	else{
	    if(prediction<=0){
		correct_neg++;
	    };
	    total_neg++;
	};
	if(alpha != 0){
	    svs++;
	    if((alpha == Cpos) || (alpha == -Cneg)){
		bsv++;
	    };
	};
      };
      min_lambda = -min_lambda;
      logln(1,"Error on KKT is "+min_lambda);
      logln(1,svs+" SVs");
      logln(1,bsv+" BSVs");
      logln(1,"Accuracy : "+((double)(correct_pos+correct_neg)/(double)(total_pos+total_neg)));
      logln(1,"Precision: "+((double)correct_pos/(double)(correct_pos+total_neg-correct_neg)));
      logln(1,"Recall   : "+((double)correct_pos/(double)total_pos));
      logln(1,"Pred:\t+\t-");
      logln(1,"\t"+correct_pos+"\t"+(total_pos-correct_pos)+"\t(true pos)");
      logln(1,"\t"+(total_neg-correct_neg)+"\t"+correct_neg+"\t(true neg)");


      if(verbosity >= 2){
	  // print hyperplane
	  double[] w = new double[dim];
	  for(j=0;j<dim;j++) w[j] = 0;
	  for(i=0;i<examples_total;i++){
	      x = the_container.get_example(i);
	      alpha = the_container.get_alpha(i);
	      for(j=0;j<dim;j++){
		  w[j] += alpha*x[j];
	      };
	  };
	  double[] Exp = the_container.Exp;
	  double[] Dev = the_container.Dev;
	  if(Exp != null){
	      for(j=0;j<dim;j++){
		  if(Dev[j] != 0){
		      w[j] /= Dev[j];
		  };
		  if(0 != Dev[dim]){
		      w[j] *= Dev[dim];
		  };
		  b -= w[j]*Exp[j];
	      };
	      b += Exp[dim];
	  };
	  logln(2," ");
	  for(j=0;j<dim;j++){
	      logln(2,"w["+j+"] = "+w[j]);
	  };
	  logln(2,"b = "+b);
	  if(dim==1){
	      logln(2,"y = "+w[0]+"*x+"+b);
	  };
      };
    };
};