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