loqo_solver.cpp

支持向量分类算法在linux操作系统下的是实现

	{
	  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");