printf("c=");  for(int i=0; i<n; ++i)    printf(" %g", c[i]);  printf("\n");  printf("low=");  for(int i=0; i<n; ++i)    printf(" %g", low[i]);  printf("\n");  printf("up=");  for(int i=0; i<n; ++i)    printf(" %g", up[i]);  printf("\n");}double Solver_LOQO_NU::calculate_rho(){  int nr_free1 = 0,nr_free2 = 0;  double ub1 = INF, ub2 = INF;  double lb1 = -INF, lb2 = -INF;  double sum_free1 = 0, sum_free2 = 0;  for(int i=0;i<active_size;i++)    {      if(y[i]==+1)	{	  if(is_lower_bound(i))	    ub1 = min(ub1,G[i]);	  else if(is_upper_bound(i))	    lb1 = max(lb1,G[i]);	  else	    {	      ++nr_free1;	      sum_free1 += G[i];	    }	}      else	{	  if(is_lower_bound(i))	    ub2 = min(ub2,G[i]);	  else if(is_upper_bound(i))	    lb2 = max(lb2,G[i]);	  else	    {	      ++nr_free2;	      sum_free2 += G[i];	    }	}    }  printf("nr_free1 = %d\n", nr_free1);  printf("sum_free1 = %g\n",sum_free1);  printf("nr_free2 = %d\n", nr_free2);  printf("sum_freee = %g\n",sum_free2);  double r1,r2;  if(nr_free1 > 0)    r1 = sum_free1/nr_free1;  else    r1 = (ub1+lb1)/2;	  if(nr_free2 > 0)    r2 = sum_free2/nr_free2;  else    r2 = (ub2+lb2)/2;	  si->r = (r1+r2)/2;  printf("(r1+r2)/2 = %g\n", (r1+r2)/2);  printf("(r1+r2)/2 = %g\n", (r1-r2)/2);  return (r1-r2)/2;}void Solver_LOQO::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;	}    }}int Solver_LOQO::solve_inner(){  Solver_NU sl;  schar *sy = new schar[n];  Qfloat *Q_bb_f = new Qfloat[n*n];  Qfloat *QD_f = new Qfloat[n];  for(int i=0; i<n; ++i)    {      sy[i] = (schar) a[i];      QD_f[i] = (Qfloat) Q_bb[i*n+i];      for(int j=0; j<n; ++j)	Q_bb_f[i*n+j] = (Qfloat) Q_bb[i*n+j];    }  sl.Solve(n, SVQ_No_Cache(Q_bb_f, QD_f, n), c, sy, work_space,	   Cp, Cn, eps, si, /* shrinking */ 0);  delete[] sy;  delete[] QD_f;  return 0;}// int Solver_LOQO::solve_inner()// {//   int result = 0;//   double sigdig;//   double epsilon_loqo = 1e-10;//   double margin;//   int iteration;//   for(margin=init_margin, iteration=init_iter; //       margin <= 0.9999999 && result != OPTIMAL_SOLUTION;)//     {//       sigdig = -log10(opt_precision);//       // run pr_loqo//       result = pr_loqo(n, m, c, Q_bb, a, d, low, up, work_space,// 		       &work_space[3*n], dist, 3, sigdig, iteration, margin, // 		       up[0]/4, 0);//       // Note: No need to check for choldc problem, as we employ//       // Manteuffel shifting.//       if(result == -1)// 	{// 	  info("NOTICE: Restarting PR_LOQO with more conservative");// 	  info("parameters.\n"); info_flush();// 	  if(init_margin < 0.80)// 	    init_margin = (4.0*margin+1.0)/5.0;// 	  margin = (margin+1.0)/2.0;// 	  opt_precision *= 10.0;// 	  info("NOTICE: Reducing precision of PR_LOQO.\n");// 	}//       else if(result != OPTIMAL_SOLUTION)// 	{// 	  // increase number of iterations// 	  iteration += 2000;// 	  init_iter += 10;// 	  // reduce precision// 	  opt_precision *= 10.0;// 	  info("NOTICE: Reducing precision of PR_LOQO!\n"); info_flush();// 	}//     }//   // Check precision of alphas//   for(int i=0; i<n; ++i)//     {//       if(work_space[i] < up[i]-epsilon_loqo && dist[i] < -eps)// 	epsilon_loqo = 2*(up[i]-work_space[i]);//       else if(work_space[i] > low[i]+epsilon_loqo && dist[i] > eps)// 	epsilon_loqo = 2*work_space[i];//     }//   info("Using epsilon_loqo= %g\n", epsilon_loqo); info_flush();//   // Clip alphas to bounds//   for(int i=0; i<n; ++i)//     {//       if(fabs(work_space[i]) < epsilon_loqo)// 	{// 	  work_space[i] = 0;// 	}//       if(fabs(work_space[i]-get_C(i)) < epsilon_loqo)// 	{// 	  work_space[i] = get_C(i);// 	}//     }//   // Compute obj after optimization//   double obj_after = 0.0;//   for(int i=0; i<n; ++i)//     {//       obj_after += work_space[i]*c[i];//       obj_after += 0.5*work_space[i]*work_space[i]*Q_bb[i*n+i];//       for(int j=0; j<i; ++j)// 	{// 	  obj_after += work_space[j]*work_space[i]*Q_bb[j*n+i];// 	}//     }//   printf("obj_after/before = %g / %g\n", obj_after, obj_before); fflush(stdout);//   printf("delta a/b = %g\n", obj_after-obj_before); fflush(stdout);//   // Check for progress//   if(obj_after >= obj_before)//     {//       // Increase precision//       opt_precision /= 100.0;//       ++precision_violations;//       info("NOTICE: Increasing Precision of PR_LOQO.\n"); info_flush();//     }//   if(precision_violations > 500)//     {//       // Relax stopping criterion//       eps *= 10.0;//       precision_violations = 0;//       info("WARNING: Relaxing epsilon on KKT-Conditions.\n"); info_flush();//     }//   return result;// }void Solver_LOQO::Solve(int l, const QMatrix& Q, const double *b_, 			const schar *y_, double *alpha_, double Cp, 			double Cn, double eps, SolutionInfo* si, 			int shrinking){  this->l = l;  this->si = si;  this->Q = &Q;  QD = Q.get_QD(); // we need diagonal for working_set_selection  clone(b,b_,l);  clone(y,y_,l);  clone(alpha,alpha_,l);  this->Cp = Cp;  this->Cn = Cn;  this->eps = eps;  this->active_size = l;  this->lmn = l - n;  int *work_set = new int[n];  int *not_work_set = new int[l];  double *delta_alpha = new double[n];  // init alpha status and work status  {    alpha_status = new char[l];    work_status = new char[l];    work_count = new int[l];    for(int i=0; i<l; ++i)      {	update_alpha_status(i);	work_status[i] = WORK_N;	work_count[i] = -1;      }  }  // init gradient  info("initializing gradient..."); info_flush();  {    G = new double[l];    for(int i=0; i<l; ++i)      {	G[i] = b[i];      }    for(int i=0; i<l; ++i)      { 	if(!is_lower_bound(i)) 	  {	    const Qfloat *Q_i = Q.get_Q(i,l); 	    double alpha_i = alpha[i];	    for(int j=0; j<l; ++j)	      {		G[j] += alpha_i * Q_i[j];	      } 	  }      }  }  info("done.\n"); info_flush();  // Allocate space for pr_loqo  Q_bb = new LOQOfloat[n*n];  c = new LOQOfloat[n];  up = new LOQOfloat[n];  low = new LOQOfloat[n];  dist = new LOQOfloat[n];  allocate_a();  allocate_d();  allocate_work_space();  iter=0;  // optimization loop  while(1)    {      // select working set and check for optimality      if(iter > 0)	{	  if(select_working_set(work_set, not_work_set) != 0)	    break;	}      else	{	  info("initializing working set..."); info_flush();   	  init_working_set(work_set, not_work_set);	  info("done.\n"); info_flush(); 	}      lmn = l - n;//       printf("working_set=");//       for(int i=0; i<n; ++i)// 	{// 	  printf("%d ",work_set[i]);// 	}//       printf("\n");//       printf("not_working_set=");//       for(int i=0; i<lmn; ++i)// 	{// 	  printf("%d ",not_work_set[i]);// 	}//       printf("\n");      ++iter;      // setup problem for pr_loqo      info("setting up problem for pr_loqo..."); info_flush();      setup_problem(work_set, not_work_set);      setup_up(work_set);      setup_low();//       print_problem();      info("done.\n");      // run inner solver      for(int i=0; i<n; ++i)	work_space[i]=alpha[work_set[i]];      // Compute obj before optimization      obj_before = 0.0;      for(int i=0; i<n; ++i)	{	  obj_before += alpha[work_set[i]]*c[i];	  obj_before += 0.5*alpha[work_set[i]]*alpha[work_set[i]]*Q_bb[i*n+i];	  for(int j=0; j<i; ++j)	    {	      obj_before += alpha[work_set[j]]*alpha[work_set[i]]*Q_bb[j*n+i];	    }	}      printf("obj_before = %g\n", obj_before); fflush(stdout);      int status = solve_inner();      printf("pr_loqo status = %d\n",status);      // Restore Q_bb, lower triangle, overwritten by pr_loqo      for(int i=0; i<n; ++i)	{	  for(int j=i+1; j<n; ++j)	    {	      Q_bb[n*j+i] = Q_bb[n*i+j];	    }	}      // update gradient, compute G_b += Q_bb*delta_alpha and       // G_n += Q_nb*delta_alpha.      int *nz = new int[n];      double sum_delta_alpha_pos=0;      double sum_delta_alpha_neg=0;      for(int i=0; i<n; ++i)	{	  delta_alpha[i] = work_space[i] - alpha[work_set[i]];	  if(y[work_set[i]]==+1)	    sum_delta_alpha_pos += delta_alpha[i];	  else	    sum_delta_alpha_neg += delta_alpha[i];   	  if(fabs(delta_alpha[i]) > TOL_ZERO)	    nz[i] = 1;  	  else  	    nz[i] = 0;	}      printf("sum_delta_alpha_pos = %g\n", sum_delta_alpha_pos);      printf("sum_delta_alpha_neg = %g\n", sum_delta_alpha_neg);      for(int i=0; i<n; ++i)	for(int j=0; j<n; ++j) // G_b	  {	    if(nz[j])	      G[work_set[i]] += Q_bb[n*i+j]*delta_alpha[j];	  }      for(int i=0; i<n; ++i) // G_n	{	  if(nz[i])	    { 	      const Qfloat *Q_i =  		Q.get_Q_subset(work_set[i],not_work_set,lmn);//  	      const Qfloat *Q_i = Q.get_Q(work_set[i],l);	      for(int j=0; j<lmn; ++j)		G[not_work_set[j]] += Q_i[not_work_set[j]] * delta_alpha[i];	    }	}      delete[] nz;      // update alpha      for(int i=0; i<n; ++i)	{ 	  alpha[work_set[i]] = work_space[i];	  update_alpha_status(work_set[i]);	}//       double sum_alpha=0;//       printf("alpha = ");//       double sum_alpha_pos=0;//       double sum_alpha_neg=0;//       for(int i=0; i<l; ++i)// 	{// 	  sum_alpha += alpha[i];// 	  if(y[i]==+1)// 	    sum_alpha_pos += alpha[i];// 	  else// 	    sum_alpha_neg += alpha[i];// 	  printf(" %g", alpha[i]);// 	}//       printf("\n");//       printf("sum_alpha = %g\n",sum_alpha);//       printf("sum_alpha_pos = %g\n",sum_alpha_pos);//       printf("sum_alpha_neg = %g\n",sum_alpha_neg);//       double delta_alpha_pos = (sum_alpha/2)-sum_alpha_pos;//       double delta_alpha_neg = (sum_alpha/2)-sum_alpha_neg;      // Restore consistency//       if(fabs(delta_alpha_pos) > 0)// 	{// 	  for(int i=0; i<l; ++i)// 	    {// 	      if(y[i] == +1)// 		{// 		  if(alpha[i]+delta_alpha_pos >= low[i] &&// 		     alpha[i]+delta_alpha_pos <= up[i])// 		    {// 		      alpha[i] += delta_alpha_pos;// 		      break;// 		    }// 		}// 	    }// 	}//       if(fabs(delta_alpha_neg) > 0)// 	{// 	  for(int i=0; i<l; ++i)// 	    {// 	      if(y[i] == -1)// 		{// 		  if(alpha[i]+delta_alpha_neg >= low[i] &&// 		     alpha[i]+delta_alpha_neg <= up[i])// 		    {// 		      alpha[i] += delta_alpha_neg;// 		      break;// 		    }// 		}// 	    }// 	}      // Recheck//       sum_alpha_pos = 0;//       sum_alpha_neg = 0;//       sum_alpha = 0;//       for(int i=0; i<l; ++i)// 	{// 	  sum_alpha += alpha[i];// 	  if(y[i] == +1)// 	    sum_alpha_pos += alpha[i];// 	  else// 	    sum_alpha_neg += alpha[i];// 	}//       printf("sum_alpha = %g\n",sum_alpha);//       printf("sum_alpha_pos = %g\n",sum_alpha_pos);//       printf("sum_alpha_neg = %g\n",sum_alpha_neg);    } // while(1)  // calculate rho  si->rho=calculate_rho();  // calculate objective value  {    double v = 0;    int i;    for(i=0;i<l;i++)      v += alpha[i] * (G[i] + b[i]);    si->obj = v/2;  }  // put back the solution  {    for(int i=0;i<l;i++)      alpha_[i] = alpha[i];  }  si->upper_bound_p = Cp;  si->upper_bound_n = Cn;  info("\noptimization finished, #iter = %d\n",iter);  // clean up  delete[] b;  delete[] y;  delete[] alpha;  delete[] alpha_status;  delete[] delta_alpha;  delete[] work_status;  delete[] work_count;  delete[] work_set;  delete[] not_work_set;  delete[] G;  delete[] Q_bb;  delete[] c;  delete[] up;  delete[] low;  delete[] dist;  delete[] a;  delete[] d;  delete[] work_space;}