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

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

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

  time_optimize += get_time() - time_start;
};


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


/**
 *
 * svm_nu_pattern_c
 *
 **/

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


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


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


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


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

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

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

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

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

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

  svm_c::init_working_set();
};


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


/**
 *
 * svm_distribution_c
 *
 **/

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


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


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

  SVMFLOAT alpha;
  SVMFLOAT result = 0;

  alpha=all_alphas[i];

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

  return result;
};


int svm_distribution_c::feasible(const SVMINT i){
  // is direction i feasible to minimize the target function
  // (includes which_alpha==0)

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

  SVMFLOAT alpha;
  SVMFLOAT result;

  alpha=all_alphas[i];

  if(alpha-Cneg >= - is_zero){
    // alpha* at upper bound
    result = -lambda_eq - sum[i];
    if(result>=-feasible_epsilon){
      return 0; 
    };
  }
  else if(alpha<=is_zero){
    // lower bound active
    result = sum[i]+lambda_eq;
    if(result>=-feasible_epsilon){
      return 0; 
    };
  }
  else{
    // not at bound
    result= abs(sum[i]+lambda_eq);
    if(result<=feasible_epsilon){
      return 0; 
    };
  };
  return 1;
};


int svm_distribution_c::feasible(const SVMINT i, SVMFLOAT* the_nabla, SVMFLOAT* the_lambda, int* atbound){
  // is direction i feasible to minimize the target function
  // (includes which_alpha==0)
  int is_feasible=1;

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

  SVMFLOAT alpha;

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

  return is_feasible;
};


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

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


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


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

  // actual convergence-test
  total = 0; alpha_sum=0;

  for(i=0;i<examples_total;i++){
    alpha = all_alphas[i];
    alpha_sum += alpha;
    if((alpha>is_zero) && (alpha-Cneg < -is_zero)){
      // alpha^* = - nabla
      the_lambda_eq += -sum[i];
      total++;
    };
  };

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

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

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

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

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


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

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

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

  calculate_working_set();
  update_working_set();
};


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


//   cout<<"upper bound = "<<nu<<endl;
//   cout<<"lower bound = "<<1/((SVMFLOAT)(examples_total*(examples_total-1)))<<endl;
//   cout<<"max_norm_x = "<<max_norm_x<<endl;
//   for(SVMINT i=0;i<examples_total;i++){
//     the_y = predict(i);
//     if(all_alphas[i] != 0){
//       // loo error?
//       if(the_y < 0){
// 	xi_i = -the_y;
//       }
//       else{
// 	xi_i = 0;
//       };
//       // if(xi_i+(all_alphas[i]+1/((SVMFLOAT)examples_total))*max_norm_x > rho){
//       cout<<all_alphas[i]<<"\t"<<xi_i/2<<" + "<<(all_alphas[i]+nu)*max_norm_x<<" = "<<xi_i/2+(all_alphas[i]+nu)*max_norm_x<<"\t"<<rho<<endl;
//       if(xi_i/2+(all_alphas[i]+nu)*max_norm_x >= rho){
// 	estim_loo--;
// 	estim_loo2--;
//       }
//       else if(all_alphas[i]>nu){
// 	estim_loo2--;
//       }
//       else if(all_alphas[i]*(SVMFLOAT)(examples_total*(examples_total-1))<1){
// 	estim_loo2--;
//       };
//     };
//     examples->put_y(i,the_y);
//     if(the_y>=0){
//       count++;
//     };
//   };
//   examples->set_initialised_y();
//   if(nu*((SVMFLOAT)examples_total-1)>svs-1){
//     estim_loo2 = 0;
//   };
  //cout<<"opt. loo-estim. of the support   : "<<estim_loo<<" ("<<((SVMINT)(10000.0*(SVMFLOAT)estim_loo/((SVMFLOAT)examples_total)))/100.0<<"%)."<<endl;
  //cout<<"pess. loo-estim. of the support  : "<<estim_loo2<<" ("<<((SVMINT)(10000.0*(SVMFLOAT)estim_loo2/((SVMFLOAT)examples_total)))/100.0<<"%)."<<endl;
  cout<<"examples in distribution support : "<<count<<" ("<<((SVMINT)(10000.0*(SVMFLOAT)count/((SVMFLOAT)examples_total)))/100.0<<"%)."<<endl;


};