.net

Development Platform:
Java

svm_c.cpp：Code Content
							#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();

  exit_optimizer();
  examples->set_initialised_alpha();
  parameters->working_set_size = param_wss;
  return the_result;
};


void svm_c::shrink(){
  // move shrinked examples to back
  if(to_shrink>examples_total/10){
    SVMINT i;
    SVMINT last_pos=examples_total;
    if(last_pos > working_set_size){
	for(i=0;i<last_pos;i++){
	    if(at_bound[i] >= shrink_const){
		// shrinxk i
		sum_alpha += all_alphas[i];
		last_pos--;
		examples->swap(i,last_pos);
		kernel->overwrite(i,last_pos);
		sum[i] = sum[last_pos];
		at_bound[i] = at_bound[last_pos];
		if(last_pos <= working_set_size){
		    break;
		};
	    };
	};
	to_shrink=0;
	examples_total = last_pos;
	kernel->set_examples_size(examples_total);
    };
    if(parameters->verbosity>=4){
      cout<<"shrinked to "<<examples_total<<" variables"<<endl;
    };
  };
};


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

  if(biased){
    // calc lambda_eq
    for(i=0;i<examples_total;i++){
      alpha = all_alphas[i];
      alpha_sum += alpha;
      if((alpha>is_zero) && (alpha < Cneg)){ //-Cneg < -is_zero)){
	// alpha^* = - nabla
	the_lambda_eq += -nabla(i); //all_ys[i]-epsilon_neg-sum[i];
	total++;
      }
      else if((alpha<-is_zero) && (alpha > -Cpos)){ //+Cpos > is_zero)){
	// alpha = nabla
	the_lambda_eq += nabla(i); //all_ys[i]+epsilon_pos-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; //(lambda_eq+4*lambda_WS)/5;
      if(parameters->verbosity>= 4){
	cout<<"*** no SVs in convergence(), lambda_eq = "<<lambda_eq<<"."<<endl;
      };
    };
    
    if(target_count>2){
      // estimate lambda from WS
      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];
	  if(is_alpha_neg(i) > 0){
	    lambda_eq = -nabla(i);
	  }
	  else{
	    lambda_eq = nabla(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) > convergence_epsilon){
      // equality constraint violated
      project_to_constraint();
      if(parameters->verbosity>= 4){
	cout<<"No convergence: equality constraint violated: |"<<(alpha_sum+sum_alpha)<<"| >> 0"<<endl;
      };
      result = 0;  
    };
  }
  else{
    // not biased
    lambda_eq = 0.0;
  };

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

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


void svm_c::minheap_heapify(SVMINT start, SVMINT size){
  // build heap of array working_set[start:start+size-1]
  // (i.e. "size" elements starting at "start"th element)

  // minheap = 1 <=> maximal element at root 
  // (i.e. we build the heap of minimal elements)

  // v_a[i] = w_s_v[start-1+i], count beginning at 1
  SVMFLOAT* value_array = working_set_values+start-1; 
  SVMINT* pos_array = working_set+start-1;

  int running = 1;
  SVMINT pos = 1;
  SVMINT left, right, largest;
  SVMFLOAT dummyf;
  SVMINT dummyi;
  while(running){
      left = 2*pos;
      right = left+1;
      if((left<=size) && 
	 (value_array[left] > value_array[pos]))
	largest = left;
      else{
	largest = pos;
      };
      if((right<=size) && 
	 (value_array[right] > value_array[largest])){
	largest = right;
      };
      if(largest == pos){
	running = 0;
      }
      else{
	//cout<<"switching "<<pos<<" and "<<largest<<endl;
	dummyf = value_array[pos];
	dummyi = pos_array[pos];
	value_array[pos] = value_array[largest];
	pos_array[pos] = pos_array[largest];
	value_array[largest] = dummyf;
	pos_array[largest] = dummyi;
	pos = largest;
      };
  };
};


void svm_c::maxheap_heapify(SVMINT start, SVMINT size){
  // build heap of array working_set[start:start+size-1]
  // (i.e. "size" elements starting at "start"th element)

  // minheap = 1 <=> maximal element at root 
  // (i.e. we build the heap of minimal elements)

  // v_a[i] = w_s_v[start-1+i], count beginning at 1
  SVMFLOAT* value_array = working_set_values+start-1; 
  SVMINT* pos_array = working_set+start-1;

  int running = 1;
  SVMINT pos = 1;
  SVMINT left, right, largest;
  SVMFLOAT dummyf;
  SVMINT dummyi;
  while(running){
      left = 2*pos;
      right = left+1;
      if((left<=size) && 
	 (value_array[left] < value_array[pos])){
	largest = left;
      }
      else{
	largest = pos;
      };
      if((right<=size) && 
	 (value_array[right] < value_array[largest])){
	  largest = right;
      };
      if(largest == pos){
	running = 0;
      }
      else{
	dummyf = value_array[pos];
	dummyi = pos_array[pos];
	value_array[pos] = value_array[largest];
	pos_array[pos] = pos_array[largest];
	value_array[largest] = dummyf;
	pos_array[largest] = dummyi;
	pos = largest;
      };
  };
};


SVMINT svm_c::maxheap_add(SVMINT size, const SVMINT element, const SVMFLOAT value){
    if(size < (working_set_size/2+working_set_size%2)){
	// add to max_heap
	working_set_values[working_set_size/2+size] = value;
	working_set[working_set_size/2+size] = element;
	size++;
	if(size == working_set_size/2+working_set_size%2){
	    // build heap
	    SVMINT j;
	    for(j=size;j>0;j--){
		maxheap_heapify(working_set_size/2+j-1,size+1-j);
	    };
	};
    }
    else if(value >= working_set_values[working_set_size/2]){
	// replace min of max_heap
	working_set_values[working_set_size/2] = value;
	working_set[working_set_size/2] = element;
	maxheap_heapify(working_set_size/2,size);
    };
    return size;
};


SVMINT svm_c::minheap_add(SVMINT size, const SVMINT element, const SVMFLOAT value){
    if(size<working_set_size/2){
	// add to min_heap
	working_set_values[size] = value;
	working_set[size] = element;
	size++;
	if(size == working_set_size/2){
	    // build heap
	    SVMINT j;
	    for(j=size;j>0;j--){
		minheap_heapify(j-1,size+1-j);
	    };
	};
    }
    else if(value < working_set_values[0]){
	// replace max of min_heap
	working_set_values[0] = value;
	working_set[0] = element;
	minheap_heapify(0,size);
    };
    return size;
};


void svm_c::calculate_working_set(){
  /**
   *
   * Find top and bottom (w.r.t. in_alpha_neg*nabla) feasible 
   * variables for working_set
   *
   */

  long time_start = get_time();

  // reset WSS
  if(working_set_size < parameters->working_set_size){
    working_set_size=parameters->working_set_size;
    if(working_set_size>examples_total) working_set_size = examples_total;
  };

  SVMINT heap_min=0;
  SVMINT heap_max=0;

  SVMINT i=0;
  SVMFLOAT sort_value;

  working_set_values[0] = infinity;
  working_set_values[working_set_size/2] = -infinity;

  SVMFLOAT the_lambda;
  SVMFLOAT the_nabla;
  int is_feasible;
  int atbound;
  SVMINT j;

  while(i<examples_total) {
    is_feasible = feasible(i,&the_nabla,&the_lambda,&atbound);
    if(0 != is_feasible){
      if(is_alpha_neg(i) > 0){
	sort_value = -the_nabla;  // - : maximum inconsistency approach
      }
      else{
	sort_value = the_nabla;
      };
      // add to heaps
      heap_min = minheap_add(heap_min,i,sort_value);
      heap_max = maxheap_add(heap_max,i,sort_value);
    };
    i++;
  };

  if(working_set_values[0] >= working_set_values[working_set_size/2]){
    if((heap_min>0) &&
       (heap_max>0)){
      // there could be the same values in the min- and maxheap,
      // sort them out (this is very unlikely)
      j=0;
      i=0;
      while(i<heap_min){
	// working_set[i] also in max-heap?
        j=working_set_size/2;
	while((j<working_set_size/2+heap_max) && 
	      (working_set[j] != working_set[i])){
	  j++;
	};
	if(j<working_set_size/2+heap_max){
          // working_set[i] equals working_set[j]
          if(heap_min<heap_max){
            // remove j from WS
	    working_set[j] = working_set[working_set_size/2-1+heap_max];
	    heap_max--;
	  }
	  else{
	    working_set[i] = working_set[heap_min-1];
	    heap_min--;
	  };
	}
	else{
          i++;
	};
      };
    };
  };

  if(heap_min+heap_max < working_set_size) {
      // condense WS
      for(i=0;i<heap_max;i++){
	  working_set[heap_min+i] = working_set[working_set_size/2+i];
      };
      working_set_size = heap_min+heap_max;
  };

  //  if((working_set_size<examples_total) && (working_set_size>0)){
  if(target_count>0){
    // convergence error on last iteration?
    // some more tests on WS
    // unlikely to happen, so speed isn't so important

    // are all variables at the bound?
    SVMINT pos_abs;
    int bounded_pos=1;
    int bounded_neg=1;
    SVMINT pos=0;
    while((pos<working_set_size) && ((1 == bounded_pos) || (1 == bounded_neg))){
      pos_abs = working_set[pos];
      if(is_alpha_neg(pos_abs) > 0){
	if(all_alphas[pos_abs]-Cneg < -is_zero){
	  bounded_pos = 0;
	};
	if(all_alphas[pos_abs] > is_zero){
	  bounded_neg = 0;
	};
      }
      else{
	if(all_alphas[pos_abs]+Cneg > is_zero){
	  bounded_neg = 0;
	};
	if(all_alphas[pos_abs] < -is_zero){
	  bounded_pos = 0;
	};
      };
      pos++;
    };
    if(0 != bounded_pos){
      // all alphas are at upper bound
      // need alpha that can be moved upward
      // use alpha with smallest lambda
      SVMFLOAT max_lambda = infinity;
      SVMINT max_pos=examples_total;
      for(pos_abs=0;pos_abs<examples_total;pos_abs++){
	if(is_alpha_neg(pos_abs) > 0){
	  if(all_alphas[pos_abs]-Cneg < -is_zero){
	    if(lambda(pos_abs) < max_lambda){
	      max_lambda = lambda(pos_abs);
	      max_pos = pos_abs;
	    };
	  };
	}
	else{
	  if(all_alphas[pos_abs] < -is_zero){
	    if(lambda(pos_abs) < max_lambda){
	      max_lambda = lambda(pos_abs);
	      max_pos = pos_abs;
	    };
	  };
	};
      };
      if(max_pos<examples_total){
	if(working_set_size<parameters->working_set_size){
	  working_set_size++;
	};
	working_set[working_set_size-1] = max_pos;
      };
    }
    else if(0 != bounded_neg){
      // all alphas are at lower bound
      // need alpha that can be moved downward
      // use alpha with smallest lambda
      SVMFLOAT max_lambda = infinity;
      SVMINT max_pos=examples_total;
      for(pos_abs=0;pos_abs<examples_total;pos_abs++){
	if(is_alpha_neg(pos_abs) > 0){
	  if(all_alphas[pos_abs] > is_zero){
	    if(lambda(pos_abs) < max_lambda){
	      max_lambda = lambda(pos_abs);
	      max_pos = pos_abs;
	    };
	  };
	}
	else{
	  if(all_alphas[pos_abs]+Cneg > is_zero){
	    if(lambda(pos_abs) < max_lambda){
	      max_lambda = lambda(pos_abs);
	      max_pos = pos_abs;
	    };
	  };
	};
      };
      if(max_pos<examples_total){
	if(working_set_size<parameters->working_set_size){
	  working_set_size++;
	};
	working_set[working_set_size-1] = max_pos;
      };
    };
  };

  if((working_set_size<parameters->working_set_size) &&
     (working_set_size<examples_total)){
    // use full working set
    SVMINT pos = (SVMINT)((SVMFLOAT)examples_total*rand()/(RAND_MAX+1.0));
    int ok;
    while((working_set_size<parameters->working_set_size) &&
	  (working_set_size<examples_total)){
      // add pos into WS if it isn't already
      ok = 1;
      for(i=0;i<working_set_size;i++){
	if(working_set[i] == pos){
	  ok=0;
	  i = working_set_size;
	};
      };
      if(1 == ok){
	working_set[working_set_size] = pos;
	working_set_size++;
      };
      pos = (pos+1)%examples_total;
    };
  };

  SVMINT ipos;
  for(ipos=0;ipos<working_set_size;ipos++){
    which_alpha[ipos] = is_alpha_neg(working_set[ipos]);
  };

  time_calc += get_time() - time_start;
  return;
};


void svm_c::project_to_constraint(){
  // project alphas to match the constraint
  if(biased){
    SVMFLOAT alpha_sum = sum_alpha;
    SVMINT SVcount=0;
    SVMFLOAT alpha;
    SVMINT i;
    for(i=0;i<examples_total;i++){
      alpha = all_alphas[i];
      alpha_sum += alpha;
      if(((alpha>is_zero) && (alpha-Cneg < -is_zero)) ||
	 ((alpha<-is_zero) && (alpha+Cpos > is_zero))){
	SVcount++;
      };
    };
    if(SVcount > 0){
      // project
      alpha_sum /= (SVMFLOAT)SVcount;
      for(i=0;i<examples_total;i++){
	alpha = all_alphas[i];
	if(((alpha>is_zero) && (alpha-Cneg < -is_zero)) ||
	   ((alpha<-is_zero) && (alpha+Cpos > is_zero))){
	  all_alphas[i] -= alpha_sum;
	};
      };
    };
  };
};


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

  project_to_constraint();
  // check bounds!
  if(examples->initialised_alpha()){
    if(parameters->verbosity >= 2){
      cout<<"Initialising variables, this may take some time."<<endl;
    };
    for(i=0; i<examples_total;i++){
      sum[i] = 0;
      at_bound[i] = 0;
      for(j=0; j<examples_total;j++){
	sum[i] += all_alphas[j]*kernel->calculate_K(i,j);
      };
    };
  }
  else{
    // skip kernel calculation as all alphas = 0
    for(i=0; i<examples_total;i++){
      sum[i] = 0;
      at_bound[i] = 0;
    };    
  };

  if(examples->initialised_alpha()){
    calculate_working_set();
  }
  else{
    // first working set is random
    j=0;
    i=0;
    while((i<working_set_size) && (j < examples_total)){
      working_set[i] = j;
      if(is_alpha_neg(j) > 0){
	which_alpha[i] = 1;
      }
      else{
	which_alpha[i] = -1;
      };
      i++;
      j++;
    };
  };   
  update_working_set();
};


void svm_c::put_optimizer_values(){
  // update nabla, sum, examples.
  // sum[i] += (primal_j^*-primal_j-alpha_j^*+alpha_j)K(i,j)
  // check for |nabla| < is_zero (nabla <-> nabla*)
  //  cout<<"put_optimizer_values()"<<endl;
  SVMINT i=0; 
  SVMINT j=0;
  SVMINT pos_i;
  SVMFLOAT the_new_alpha;
  SVMFLOAT* kernel_row;
  SVMFLOAT alpha_diff;

  long time_start = get_time();
  pos_i=working_set_size;
  while(pos_i>0){
    pos_i--;
    if(which_alpha[pos_i]>0){
      the_new_alpha = primal[pos_i];
    }
    else{
      the_new_alpha = -primal[pos_i];
    };
    // next three statements: keep this order!
    i = working_set[pos_i];
    alpha_diff = the_new_alpha-all_alphas[i];
    all_alphas[i] = the_new_alpha;

    if(alpha_diff != 0){
      // update sum ( => nabla)
      kernel_row = kernel->get_row(i);
      for(j=0;j<examples_total;j++){
	sum[j] += alpha_diff*kernel_row[j];
      };
    };
  };
  time_update += get_time() - time_start;
};


void svm_c::update_working_set(){
  long time_start = get_time();
  // setup subproblem
  SVMINT i,j;
  SVMINT pos_i, pos_j;
  SVMFLOAT* kernel_row;
  SVMFLOAT sum_WS;

  for(pos_i=0;pos_i<working_set_size;pos_i++){
    i = working_set[pos_i];

    // put row sort_i in hessian 
    kernel_row = kernel->get_row(i);
    sum_WS=0;
    //    for(pos_j=0;pos_j<working_set_size;pos_j++){
    for(pos_j=0;pos_j<pos_i;pos_j++){
      j = working_set[pos_j];
      // put all elements K(i,j) in hessian, where j in WS
      if(((which_alpha[pos_j] < 0) && (which_alpha[pos_i] < 0)) ||
	 ((which_alpha[pos_j] > 0) && (which_alpha[pos_i] > 0))){
	// both i and j positive or negative
	(qp.H)[pos_i*working_set_size+pos_j] = kernel_row[j];
	(qp.H)[pos_j*working_set_size+pos_i] = kernel_row[j];
      }
      else{
	// one of i and j positive, one negative
	(qp.H)[pos_i*working_set_size+pos_j] = -kernel_row[j];
	(qp.H)[pos_j*working_set_size+pos_i] = -kernel_row[j];
      };
    };
    for(pos_j=0;pos_j<working_set_size;pos_j++){
      j = working_set[pos_j];
      sum_WS+=all_alphas[j]*kernel_row[j];
    };
    // set main diagonal 
    (qp.H)[pos_i*working_set_size+pos_i] = kernel_row[i];

    // linear and box constraints
    if(which_alpha[pos_i]<0){
      // alpha
      (qp.A)[pos_i] = -1;
      // lin(alpha) = y_i+eps-sum_{i not in WS} alpha_i K_{ij}
      //            = y_i+eps-sum_i+sum_{i in WS}
      (qp.c)[pos_i] = all_ys[i]+epsilon_pos-sum[i]+sum_WS;
      primal[pos_i] = -all_alphas[i];
      (qp.u)[pos_i] = Cpos;
    }
    else{
      // alpha^*
      (qp.A)[pos_i] = 1;
      (qp.c)[pos_i] = -all_ys[i]+epsilon_neg+sum[i]-sum_WS;
      primal[pos_i] = all_alphas[i];
      (qp.u)[pos_i] = Cneg;
    };
  };
  if(parameters->quadraticLossNeg){
    for(pos_i=0;pos_i<working_set_size;pos_i++){
      if(which_alpha[pos_i]>0){
	(qp.H)[pos_i*(working_set_size+1)] += 1/Cneg;
	(qp.u)[pos_i] = infinity;
      };
    };
  };
  if(parameters->quadraticLossPos){
    for(pos_i=0;pos_i<working_set_size;pos_i++){
      if(which_alpha[pos_i]<0){
	(qp.H)[pos_i*(working_set_size+1)] += 1/Cpos;
	(qp.u)[pos_i] = infinity;
      };
    };
  };

  time_update += get_time() - time_start; 
};


svm_result svm_c::test(example_set_c* test_examples, int verbose){
  svm_result the_result;
  test_set = test_examples;

  SVMINT i;
  SVMFLOAT MAE=0;
  SVMFLOAT MSE=0;
  SVMFLOAT actloss=0;
  SVMFLOAT theloss=0;
  SVMFLOAT theloss_pos=0;
  SVMFLOAT theloss_neg=0;
  SVMINT countpos=0;
  SVMINT countneg=0;
  // for pattern:
  SVMINT correct_pos=0;
  SVMINT correct_neg=0;
  SVMINT total_pos=0;
  SVMINT total_neg=0;

  SVMFLOAT prediction;
  SVMFLOAT y;
  svm_example example;
  for(i=0;i<test_set->size();i++){
    example = test_set->get_example(i);
    prediction = predict(example);
    y = examples->unscale_y(test_set->get_y(i));
    MAE += abs(y-prediction);
    MSE += (y-prediction)*(y-prediction);
    actloss=loss(prediction,y);
    theloss+=actloss;
    if(y < prediction-parameters->epsilon_pos){
      theloss_pos += actloss;
      countpos++;
    }
    else if(y > prediction+parameters->epsilon_neg){
      theloss_neg += actloss;
      countneg++;
    };
    // if pattern!
    if(is_pattern){
      if(y>0){
	if(prediction>0){
	  correct_pos++;
	};
	total_pos++;
      }
      else{
	if(prediction<=0){
	  correct_neg++;
	};
	total_neg++;
      };
    };    
  };
  if(countpos != 0){
    theloss_pos /= (SVMFLOAT)countpos;
  };
  if(countneg != 0){
    theloss_neg /= (SVMFLOAT)countneg;
  };

  the_result.MAE =  MAE / (SVMFLOAT)test_set->size();
  the_result.MSE =  MSE / (SVMFLOAT)test_set->size();
  the_result.loss = theloss/test_set->size();
  the_result.loss_pos = theloss_pos;
  the_result.loss_neg = theloss_neg;
  the_result.number_svs = 0;
  the_result.number_bsv = 0;
  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);
  }
  else{
    the_result.accuracy = -1;
    the_result.precision = -1;
    the_result.recall = -1;
  };

  if(verbose){
    cout << "Average loss  : "<<(theloss/test_set->size())<<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;

    if(is_pattern){
      // output precision, recall and accuracy
      cout<<"Accuracy  : "<<the_result.accuracy<<endl;
      cout<<"Precision : "<<the_result.precision<<endl;
      cout<<"Recall    : "<<the_result.recall<<endl;
      // 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;
    };
  };
  return the_result;
};


void svm_c::optimize(){
  // optimizer-specific call
  // get time
  long time_start = get_time();

  qp.n = working_set_size;

  SVMINT i;
  SVMINT j;

  // equality constraint
  qp.b[0]=0;
  for(i=0;i<working_set_size;i++){
    qp.b[0] += all_alphas[working_set[i]];
  };

  // set initial optimization parameters
  SVMFLOAT new_target=0;
  SVMFLOAT old_target=0;
  SVMFLOAT target_tmp;
  for(i=0;i<working_set_size;i++){
    target_tmp = primal[i]*qp.H[i*working_set_size+i]/2;
    for(j=0;j<i;j++){
      target_tmp+=primal[j]*qp.H[j*working_set_size+i];
    };
    target_tmp+=qp.c[i];
    old_target+=target_tmp*primal[i];
  };

  SVMFLOAT new_constraint_sum=0;
  SVMFLOAT my_is_zero = is_zero;
  SVMINT sv_count=working_set_size;

  qp.n = working_set_size;
  // optimize
  int KKTerror=1;
  int convError=0;

  smo.set_max_allowed_error(convergence_epsilon);

  // loop while some KKT condition is not valid (alpha=0)

  int result;
  if(biased){
    result = smo.smo_solve(&qp,primal);
    lambda_WS = smo.get_lambda_eq();
  }
  else{
    result = smo.smo_solve_single(&qp,primal);
    lambda_WS = 0.0;
  };

  /////////// new
  SVMINT it=3;
  if(! is_pattern){
      // iterate optimization 3 times with changed sign on variables, if KKT conditions are not satisfied
    SVMFLOAT lambda_lo;
    while(KKTerror && (it>0)){
      KKTerror = 0;
      it--;
      for(i=0;i<working_set_size;i++){
	if(primal[i]<is_zero){
	  lambda_lo =  epsilon_neg + epsilon_pos - qp.c[i];
	  for(j=0;j<working_set_size;j++){
	    lambda_lo -= primal[j]*qp.H[i*working_set_size+j];
	  };
	  if(qp.A[i] > 0){
	    lambda_lo -= lambda_WS;
	  }
	  else{
	    lambda_lo += lambda_WS;
	  };

	  if(lambda_lo<-convergence_epsilon){
	    // change sign of i
	    KKTerror=1;
	    qp.A[i] = -qp.A[i];
	    which_alpha[i] = -which_alpha[i];
	    primal[i] = -primal[i];
	    qp.c[i] = epsilon_neg + epsilon_pos - qp.c[i];
	    if(qp.A[i]>0){
	      qp.u[i] = Cneg;
	    }
	    else{
	      qp.u[i] = Cpos;
	    };
	    for(j=0;j<working_set_size;j++){
	      qp.H[i*working_set_size+j] = -qp.H[i*working_set_size+j];
	      qp.H[j*working_set_size+i] = -qp.H[j*working_set_size+i];
	    };
	    if(parameters->quadraticLossNeg){
	      if(which_alpha[i]>0){
		(qp.H)[i*(working_set_size+1)] += 1/Cneg;
		(qp.u)[i] = infinity;
	      }
	      else{
		// previous was neg
		(qp.H)[i*(working_set_size+1)] -= 1/Cneg;
	      };
	    };
	    if(parameters->quadraticLossPos){
	      if(which_alpha[i]<0){
		(qp.H)[i*(working_set_size+1)] += 1/Cpos;
		(qp.u)[i] = infinity;
	      }
	      else{
		//previous was pos
		(qp.H)[i*(working_set_size+1)] -= 1/Cpos;
	      };
	    };
	  };
	};
      };
      result = smo.smo_solve(&qp,primal);
      if(biased){
	lambda_WS = smo.get_lambda_eq();
      }; // sonst 0 beibehalten
    };
  };

  KKTerror = 1;
  //////////////////////

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

  while(KKTerror){
    // clip
    sv_count=working_set_size;
    new_constraint_sum=qp.b[0];
    if(biased){
      // more checks for feasibility
      for(i=0;i<working_set_size;i++){
	// check if at bound
	if(primal[i] <= my_is_zero){
	  // at lower bound
	  primal[i] = qp.l[i];
	  sv_count--;
	}
	else if(qp.u[i]-primal[i] <= my_is_zero){
	  // at upper bound
	  primal[i] = qp.u[i];
	  sv_count--;
	};
	new_constraint_sum -= qp.A[i]*primal[i];
      };
      
      // enforce equality constraint
      if(sv_count>0){
	new_constraint_sum /= (SVMFLOAT)sv_count;
	if(parameters->verbosity>=5){
	  cout<<"adjusting "<<sv_count<<" alphas by "<<new_constraint_sum<<endl;
	};
	for(i=0;i<working_set_size;i++){
	  if((primal[i] > qp.l[i]) && 
	     (primal[i] < qp.u[i])){
	    // real sv
	    primal[i] += qp.A[i]*new_constraint_sum;
	  };
	};
      }
      else if(abs(new_constraint_sum)>(SVMFLOAT)working_set_size*is_zero){
	// error, can't get feasible point
	if(parameters->verbosity>=5){
	  cout<<"WARNING: No SVs, constraint_sum = "<<new_constraint_sum<<endl;
	};
	old_target = -infinity; 
	//is_ok=0;
	convError=1;
      };
    };
    // test descend
    new_target=0;
    for(i=0;i<working_set_size;i++){
      // attention: loqo changes one triangle of H!
      target_tmp = primal[i]*qp.H[i*working_set_size+i]/2;
      for(j=0;j<i;j++){
	target_tmp+=primal[j]*qp.H[j*working_set_size+i];
      };
      target_tmp+=qp.c[i];
      new_target+=target_tmp*primal[i];
    };

    if(new_target < old_target){
      KKTerror = 0;
      if(parameters->descend < old_target - new_target){
	target_count=0;
      }
      else{
	convError=1;
      };
      if(parameters->verbosity>=5){
	cout<<"descend = "<<old_target-new_target<<endl;
      };
    }
    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 = 1e20;
      for(i=0;i<working_set_size;i++){
	if((primal[i] > qp.l[i]) && (primal[i] < qp.u[i])){
	  if(primal[i] - qp.l[i] < my_is_zero){
	    my_is_zero = primal[i]-qp.l[i];
	  };
	  if(qp.u[i]  - primal[i]  < my_is_zero){
	    my_is_zero = qp.u[i] - primal[i];
	  };
	};
      };
      if(target_count == 0){
      	my_is_zero *= 2;
      };
      if(parameters->verbosity>=5){
	cout<<"WARNING: no descend ("<<old_target-new_target
	    <<" <= "<<parameters->descend
	  //	    <<", alpha_diff = "<<alpha_diff
	    <<"), adjusting is_zero to "<<my_is_zero<<endl;
	cout<<"new_target = "<<new_target<<endl;
      };
    }
    else{
      // nothing we can do
      if(parameters->verbosity>=5){
	cout<<"WARNING: no descend ("<<old_target-new_target
	    <<" <= "<<parameters->descend<<"), stopping."<<endl;
      };
      KKTerror=0;
      convError=1;
    };
  };

  if(1 == convError){
    target_count++;
    //    sigfig_max+=0.05;
    if(old_target < new_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, restoring old primals"<<endl; 
      };
    };                                          
  };

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

  time_optimize += get_time() - time_start;
};


void svm_c::predict(example_set_c* test_examples){
  test_set = test_examples;
  SVMINT i;
  SVMFLOAT prediction;
  svm_example example;

  for(i=0;i<test_set->size();i++){
    example = test_set->get_example(i);
    prediction = predict(example);
    test_set->put_y(i,prediction);
  };
  test_set->set_initialised_y();
  test_set->put_b(examples->get_b());
  if(parameters->verbosity>=4){
    cout<<"Prediction generated"<<endl;
  };
};


SVMFLOAT svm_c::predict(svm_example example){ 
  SVMINT i;
  svm_example sv;
  SVMFLOAT the_sum=examples->get_b();

  for(i=0;i<examples_total;i++){
    if(all_alphas[i] != 0){
      sv = examples->get_example(i);
      the_sum += all_alphas[i]*kernel->calculate_K(sv,example);
    };
  };
  the_sum = examples->unscale_y(the_sum);
  if(parameters->use_min_prediction){
    if(the_sum < parameters->min_prediction){
      the_sum = parameters->min_prediction;
    };
  };
  return the_sum;
};


SVMFLOAT svm_c::predict(SVMINT i){
  //   return (sum[i]+examples->get_b());
  // numerically more stable:
  return predict(examples->get_example(i));
};

SVMFLOAT svm_c::loss(SVMINT i){
  return loss(predict(i),examples->unscale_y(all_ys[i]));
};


SVMFLOAT svm_c::loss(SVMFLOAT prediction, SVMFLOAT value){
  SVMFLOAT theloss = prediction - value;
  if(is_pattern){
    if(((value > 0) && (prediction > 0)) ||
       ((value <= 0) && (prediction <= 0))){
      theloss = 0;
    }
  };
  if(theloss > parameters->epsilon_pos){ 
    if(parameters->quadraticLossPos){
      theloss = parameters->Lpos*(theloss-parameters->epsilon_pos)
	*(theloss-parameters->epsilon_pos); 
    }
    else{
      theloss =  parameters->Lpos*(theloss-parameters->epsilon_pos); 
    };
  }
  else if(theloss >= -parameters->epsilon_neg){ theloss = 0; }
  else{ 
    if(parameters->quadraticLossNeg){
      theloss = parameters->Lneg*(-theloss-parameters->epsilon_neg)
	*(-theloss-parameters->epsilon_neg);
    }
    else{
      theloss = parameters->Lneg*(-theloss-parameters->epsilon_neg);
    };
  };
  return theloss;
};


void svm_c::print_special_statistics(){
  // nothing special here!
};


svm_result svm_c::print_statistics(){
  // # SV, # BSV, pos&neg, Loss, VCdim
  // Pattern: Acc, Rec, Pred

  if(parameters->verbosity>=2){
    cout<<"----------------------------------------"<<endl;
  };

  svm_result the_result;
  if(test_set->size() <= 0){
    if(parameters->verbosity>= 0){
      cout << "No training set given" << endl;
    };
    the_result.loss = -1;
    return the_result;  // undefined
  };
  SVMINT i;
  SVMINT svs = 0;
  SVMINT bsv = 0;
  SVMFLOAT actloss = 0;
  SVMFLOAT theloss = 0;
  SVMFLOAT theloss_pos=0;
  SVMFLOAT theloss_neg=0;
  SVMINT countpos=0;
  SVMINT countneg=0;
  SVMFLOAT min_alpha=infinity;
  SVMFLOAT max_alpha=-infinity;
  SVMFLOAT norm_w=0;
  SVMFLOAT max_norm_x=0;
  SVMFLOAT min_norm_x=1e20;
  SVMFLOAT norm_x=0;
  SVMFLOAT loo_loss_estim=0;
  // for pattern:
  SVMFLOAT correct_pos=0;
  SVMINT correct_neg=0;
  SVMINT total_pos=0;
  SVMINT total_neg=0;
  SVMINT estim_pos=0;
  SVMINT estim_neg=0;
  SVMFLOAT MSE=0;
  SVMFLOAT MAE=0;
  SVMFLOAT alpha;
  SVMFLOAT prediction;
  SVMFLOAT y;
  SVMFLOAT xi;

  for(i=0;i<examples_total;i++){
    // needed before test-loop for performance estimators
    norm_w+=all_alphas[i]*sum[i];

    alpha=all_alphas[i];
    if(alpha!=0){
      norm_x = kernel->calculate_K(i,i);
      if(norm_x>max_norm_x){
	max_norm_x = norm_x;
      };
      if(norm_x<min_norm_x){
	min_norm_x = norm_x;
      };
    };
  };
  
  SVMFLOAT r_delta = max_norm_x;
  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;
};

						