{	  count += 2;	  work_status[pidx[top]] = WORK_B;	  work_status[pidx[bot]] = WORK_B;	  --top;	  ++bot;	}    }  if(count < n)    {      // Compute subset of indices in previous working set       // which were not yet selected      int j=0;      int *work_count_subset = new int[l-count];      int *subset = new int[l-count];      int *psubset = new int[l-count];      for(int i=0; i<l; ++i)	{	  if(work_status[i] == WORK_N && work_count[i] > -1)	    {	      work_count_subset[j] = work_count[i];	      subset[j] = i;	      psubset[j] = j;	      ++j;	    }	}      quick_sort(work_count_subset, psubset, 0, j-1);      // Fill up with j \in B, 0 < alpha[j] < C      for(int i=0; i<j; ++i)	{	  if(count == n) break;	  if(is_free(subset[psubset[i]])) 	    {	      work_status[subset[psubset[i]]] = WORK_B;	      ++count;	    }	}      // Fill up with j \in B, alpha[j] = 0      for(int i=0; i<j; ++i)	{	  if(count == n) break;	  if(is_lower_bound(subset[psubset[i]]))	    {	      work_status[subset[psubset[i]]] = WORK_B;	      ++count;	    }	}      // Fill up with j \in B, alpha[j] = C      for(int i=0; i<j; ++i)	{	  if(count == n) break;	  if(is_upper_bound(subset[psubset[i]]))	    {	      work_status[subset[psubset[i]]] = WORK_B;	      ++count;	    }	}      // clean up      delete[] work_count_subset;      delete[] subset;      delete[] psubset;    } // if(count < n)  // Setup work_set and not_work_set  // update work_count  int nnew=0; int i=0; int j=0; n=0;  for(int t=0; t<l; ++t)    {      if(work_status[t] == WORK_B)	{	  if(work_count[t] == -1)	    ++nnew;	  work_set[i] = t; ++i; ++n;	  ++work_count[t];	}      else	{	  not_work_set[j] = t; ++j;	  work_count[t] = -1;	}    }  // Update q  printf("nnew = %d\n", nnew);  int kin = nnew;  nnew = nnew % 2 == 0 ? nnew : nnew-1;  int L = n/10 % 2 == 0 ? n/10 : (n/10)-1;  q = min(q, max( max( 10, L ), nnew ) );  printf("q = %d\n", q);  printf("n = %d\n",n);  if(kin == 0)    {      // 1st: Increase precision of solver.      if(opt_precision > 1e-20)	opt_precision /= 100;      else	{	  info("Error: Unable to select a suitable working set!!!\n");	  return 1;	}    }  // Clean up  delete[] yG;  delete[] pidx;  printf("done.\n");  return 0;}class Solver_LOQO_NU : public Solver_LOQO{public:  // Ensure minimum working set size =3,  // next even =4.  Solver_LOQO_NU(int n, int q, int m, double nu) : Solver_LOQO(max(4,n),q, m)  {    this->nu = nu;  };  void Solve(int l, const QMatrix& Q, const double *b, const schar *y,	     double *alpha, double Cp, double Cn, double eps,	     SolutionInfo* si, int shrinking)  {    Solver_LOQO::Solve(l,Q,b,y,alpha,Cp,Cn,eps,si,shrinking);  }private:  double nu;  void allocate_a() { a = new double[2*n]; }  void allocate_d() { d = new double[2]; }  void allocate_work_space() { work_space = new double[5*n+2]; }  void setup_a(int *work_set);  void setup_d(int *not_work_set);  double calculate_rho();  //double calculate_rho(){ si->r = work_space[3*n+1]; return work_space[3*n]; }  int check_working_set();  int select_working_set(int *work_set, int *not_work_set);  void init_working_set(int *work_set, int *not_work_set);};void Solver_LOQO_NU::init_working_set(int *work_set, int *not_work_set){  int j;  NEXT_RAND = 1;  for (int i=0; i<n; ++i)  {    do      {	j = next_rand_pos() % l;      } while (work_status[j] != WORK_N);    work_status[j] = WORK_B;  }  int k=0; j=0;  for(int i=0; i<l; ++i)    {      if(work_status[i] == WORK_B)	{	  work_set[j] = i; ++j;	  work_count[i] = 0;	}      else	{	  not_work_set[k] = i; ++k;	  work_count[i] = -1;	}    }}void Solver_LOQO_NU::setup_a(int *work_set){  for(int i=0; i<n; ++i)    {      a[i] = y[work_set[i]];      a[n+i] = 1;    }}void Solver_LOQO_NU::setup_d(int *not_work_set){  d[0]=0;  for(int i=0; i<lmn; ++i)    d[0] -= y[not_work_set[i]]*alpha[not_work_set[i]];  d[1] = nu*l;  for(int i=0; i<lmn; ++i)    d[1] -= alpha[not_work_set[i]];}int Solver_LOQO_NU::select_working_set(int *work_set, int *not_work_set){  printf("selecting working set...");  // reset work status  n = n_old;  for(int i=0; i<l; ++i)    {      work_status[i] = WORK_N;    }  double Gmin1 = INF; double Gmin2 = INF;  double Gmax1 = -INF; double Gmax2 = -INF;  int min1 = -1; int min2 = -1;  int max1 = -1; int max2 = -1;  for(int t=0; t<l; ++t)    {      if(y[t] == +1)	{	  if(!is_upper_bound(t))	    {	      if(G[t] < Gmin1)		{		  Gmin1 = G[t];		  min1 = t;		}	    }	  if(!is_lower_bound(t))	    {	      if(G[t] > Gmax1)		{		  Gmax1 = G[t];		  max1 = t;		}	    }	}      else	{	  if(!is_upper_bound(t))	    {	      if(G[t] < Gmin2)		{		  Gmin2 = G[t];		  min2 = t;		}	    }	  if(!is_lower_bound(t))	    {	      if(G[t] > Gmax2)		{		  Gmax2 = G[t];		  max2 = t;		}	    }	}    }  // check for optimality, max. violating pair.  printf("max(Gmax1-Gmin1,Gmax2-Gmin2) = %g < %g\n", 	 max(Gmax1-Gmin1,Gmax2-Gmin2),eps);  if(max(Gmax1-Gmin1,Gmax2-Gmin2) < eps)    return 1;  // Sort gradient  double *Gtmp = new double[l];  int *pidx = new int[l];  for(int i=0; i<l; ++i)    {      Gtmp[i] = G[i];      pidx[i] = i;    }  quick_sort(Gtmp, pidx, 0, l-1);//   printf("pidx = ");//   for(int i=0; i<l; ++i)//     printf(" %d", pidx[i]);//   printf("\n");  int top1=l-1; int top2=l-1;  int bot1=0; int bot2=0;  int count=0;  // Select a full set initially  int nselect = iter == 0 ? n : q;  while(count < nselect)    {      if(top1 > bot1)	{ 	  while(!( (is_free(pidx[top1]) || is_upper_bound(pidx[top1])) 		   && y[pidx[top1]] == +1))	    {	      if(top1 <= bot1) break;	      --top1;	    } 	  while(!( (is_free(pidx[bot1]) || is_lower_bound(pidx[bot1])) 		   && y[pidx[bot1]] == +1))	    {	      if(bot1 >= top1) break;	      ++bot1;	    }	}      if(top2 > bot2)	{	  while(!( (is_free(pidx[top2]) || is_upper_bound(pidx[top2]))		   && y[pidx[top2]] == -1))	    {	      if(top2 <= bot2) break;	      --top2;	    }	  while(!( (is_free(pidx[bot2]) || is_lower_bound(pidx[bot2]))		   && y[pidx[bot2]] == -1))	    {	      if(bot2 >= top2) break;	      ++bot2;	    }	}      if(top1 > bot1 && top2 > bot2)	{	  if(G[pidx[top1]]-G[pidx[bot1]] > G[pidx[top2]]-G[pidx[bot2]])	    {	      work_status[pidx[top1]] = WORK_B;	      work_status[pidx[bot1]] = WORK_B;	      --top1;	      ++bot1;	    }	  else	    {	      work_status[pidx[top2]] = WORK_B;	      work_status[pidx[bot2]] = WORK_B;	      --top2;	      ++bot2;	    }	  count += 2;	}      else if(top1 > bot1)	{	  work_status[pidx[top1]] = WORK_B;	  work_status[pidx[bot1]] = WORK_B;	  --top1;	  ++bot1;	  count += 2;	}      else if(top2 > bot2)	{	  work_status[pidx[top2]] = WORK_B;	  work_status[pidx[bot2]] = WORK_B;	  --top2;	  ++bot2;	  count += 2;	}      else	break;    } // while(count < nselect)  if(count < n)    {      // Compute subset of indices in previous working set       // which were not yet selected      int j=0;      int *work_count_subset = new int[l-count];      int *subset = new int[l-count];      int *psubset = new int[l-count];      for(int i=0; i<l; ++i)	{	  if(work_status[i] == WORK_N && work_count[i] > -1)	    {	      work_count_subset[j] = work_count[i];	      subset[j] = i;	      psubset[j] = j;	      ++j;	    }	}      quick_sort(work_count_subset, psubset, 0, j-1);      // Fill up with j \in B, 0 < alpha[j] < C      for(int i=0; i<j; ++i)	{	  if(count == n) break;	  if(is_free(subset[psubset[i]])) 	    {	      work_status[subset[psubset[i]]] = WORK_B;	      ++count;	    }	}      // Fill up with j \in B, alpha[j] = 0      for(int i=0; i<j; ++i)	{	  if(count == n) break;	  if(is_lower_bound(subset[psubset[i]]))	    {	      work_status[subset[psubset[i]]] = WORK_B;	      ++count;	    }	}      // Fill up with j \in B, alpha[j] = C      for(int i=0; i<j; ++i)	{	  if(count == n) break;	  if(is_upper_bound(subset[psubset[i]]))	    {	      work_status[subset[psubset[i]]] = WORK_B;	      ++count;	    }	}      // clean up      delete[] work_count_subset;      delete[] subset;      delete[] psubset;    } // if(count < n)  // Setup work_set and not_work_set  // update work_count  int nnew=0; int i=0; int j=0; n=0;  for(int t=0; t<l; ++t)    {      if(work_status[t] == WORK_B)	{	  if(work_count[t] == -1)	    ++nnew;	  work_set[i] = t; ++i; ++n;	  ++work_count[t];	}      else	{	  not_work_set[j] = t; ++j;	  work_count[t] = -1;	}    }  // Update q  printf("nnew = %d\n", nnew);  int kin = nnew;  nnew = nnew % 2 == 0 ? nnew : nnew-1;  int L = n/10 % 2 == 0 ? n/10 : (n/10)-1;  q = min(q, max( max( 10, L ), nnew ) );  printf("q = %d\n", q);  printf("n = %d\n",n);  if(kin == 0)    {      // 1st: Increase precision of solver.      if(opt_precision > 1e-20)	opt_precision /= 100;      else	{	  info("Error: Unable to select a suitable working set!!!\n");	  return 1;	}    }  // clean up  delete[] Gtmp;  delete[] pidx;  printf("done.\n");  return 0;  }void Solver_LOQO::setup_a(int *work_set){  info("Solver_LOQO::setting up a..."); info_flush();  for(int i=0; i<n; ++i)    {      a[i] = y[work_set[i]];    }  info("done.\n"); info_flush();}void Solver_LOQO::setup_d(int *not_work_set){  info("setting up d..."); info_flush();  d[0] = 0;  for(int i=0; i<lmn; ++i)    {      if(fabs(alpha[not_work_set[i]]) > TOL_ZERO)	d[0] -= y[not_work_set[i]]*alpha[not_work_set[i]];    }  info("done.\n"); info_flush();}void Solver_LOQO::setup_problem(int *work_set, int *not_work_set){  info("setting up Q_bb..."); info_flush();  for(int i=0; i<n; ++i)    {      const Qfloat *Q_i = Q->get_Q_subset(work_set[i],work_set,n);      //  	  const Qfloat *Q_i = Q.get_Q(work_set[i],l);      c[i] = G[work_set[i]];      for(int j=0; j<n; ++j)	{	  Q_bb[i*n+j] = (double) Q_i[work_set[j]];	  //   	      Q_bb[i*n+j] = Q_i[work_set[j]];	  if(alpha[work_set[j]] > TOL_ZERO)	    c[i] -= Q_bb[i*n+j]*alpha[work_set[j]];	}    }  setup_a(work_set);  setup_d(not_work_set);  info("done.\n"); info_flush();}void Solver_LOQO::print_problem(){  printf("G=");  for(int i=0; i<l; ++i)    {      printf(" %g", G[i]);    }  printf("\n");//   printf("Q_bb=");//   for(int i=0; i<n; ++i)//     {//       for(int j=0; j<n; ++j)// 	printf(" %g", Q_bb[i*n+j]);//       printf("\n");//     }  printf("d[0] = % g\n",d[0]);  printf("d[1] = % g\n",d[1]);  printf("a=");  for(int i=0; i<n; ++i)    printf(" %g", a[i]);  printf("\n");  printf("a2=");  for(int i=0; i<n; ++i)    printf(" %g", a[i+n]);  printf("\n");