psmo_solver.cpp

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

// D. Brugger, december 2006
// $Id: psmo_solver.cpp 7 2006-12-16 16:45:32Z beeblbrox $
//
// Copyright (C) 2006 Dominik Brugger
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

#include "psmo_solver.h"

#define CheckError(n) if(n){printf("line %d, file %s\n",__LINE__,__FILE__);}
#define MPIfloat MPI_DOUBLE

Solver_Parallel_SMO::Solver_Parallel_SMO(int n, int q, MPI_Comm comm)
{
  // Ensure that n,q are even numbers.
  this->n = n % 2 == 0 ? n : n-1; 
  this->n_old = this->n; 
  this->q = q % 2 == 0 ? q : q-1; 
  // Ensure sane q
  this->q = this->q > this->n ? this->n : this->q;
  this->comm = comm;
  ierr = MPI_Comm_rank(comm, &this->rank); CheckError(ierr);
  ierr = MPI_Comm_size(comm, &this->size); CheckError(ierr);
  NEXT_RAND = 1;
}

unsigned int Solver_Parallel_SMO::next_rand_pos()
{
  NEXT_RAND = NEXT_RAND*1103515245L + 12345L;
  return NEXT_RAND & 0x7fffffff;
}

Solver_Parallel_SMO_NU::Solver_Parallel_SMO_NU(int n, int q, MPI_Comm comm)
  : Solver_Parallel_SMO(n,q,comm)
{}

double Solver_Parallel_SMO_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;

//   printf("alpha = ");
//   for(int i=0; i<l; ++i)
//     printf(" %g",alpha[i]);
//   printf("\n");

  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_Parallel_SMO_NU::solve_inner()
{
  Solver_NU sl;
  sl.Solve(n, SVQ_No_Cache(Q_bb, QD_b, n), c, a, alpha_b, Cp, Cn, eps,
	   si, /* shrinking */ 0);
}

int Solver_Parallel_SMO_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(eps > 1e-10)
	eps /= 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_Parallel_SMO::solve_inner()
{
  Solver s;
  s.Solve(n, SVQ_No_Cache(Q_bb, QD_b, n), c, a, alpha_b, Cp, Cn, eps,
	  si, /* shrinking */ 0);
}

int Solver_Parallel_SMO::select_working_set(int *work_set, int *not_work_set)
{
  // printf("selecting working set...");
  // reset work status
  n = n_old;
  int *old_work_set = new int[n];
  for(int i=0; i<n; ++i)
    old_work_set[i] = work_set[i];
  for(int i=0; i<l; ++i)
      work_status[i] = WORK_N;

  double Gmax1 = -INF;		// max { -y_i * grad(f)_i | i in I_up(\alpha) }
  double Gmax2 = -INF;		// max { y_i * grad(f)_i | i in I_low(\alpha) }
  for(int i=0; i<l; ++i)
    {
      if(!is_upper_bound(i))
	{
	  if(y[i] == +1)
	    {
	      if(-G[i] > Gmax1)
		Gmax1 = -G[i];
	    }
	  else
	    {
	      if(-G[i] > Gmax2)
		Gmax2 = -G[i];
	    }
	}
      if(!is_lower_bound(i))
	{
	  if(y[i] == +1)
	    {
	      if(G[i] > Gmax2)
		Gmax2 = G[i];
	    }
	  else
	    {
	      if(G[i] > Gmax1)
		Gmax1 = G[i];
	    }
	}
    }
  // check for optimality
  //  printf("Gmax1 + Gmax2 = %g < %g\n", Gmax1+Gmax2,eps);
  info(" %g\n", Gmax1+Gmax2);
  if(Gmax1 + Gmax2 < eps)
      return 1;

  // Compute yG
  double *yG = new double[l];
  int *pidx = new int[l];
  for(int i=0; i<l; ++i)
    {
      if(y[i] == +1)
	yG[i] = G[i];
      else
	yG[i] = -G[i];
      pidx[i] = i;
    }
  //  double sort_time = MPI_Wtime();
  quick_sort(yG, pidx, 0, l-1);
  // printf("sort_time = %g\n", MPI_Wtime() - sort_time);
//   printf("yG = ");
//   for(int i=0; i<l; ++i)
//     printf(" %g",yG[i]);
//   printf("\n");
//   printf("pidx = ");
//   for(int i=0; i<l; ++i)
//     printf(" %d",pidx[i]);
//   printf("\n");
  int top=l-1; int bot=0; int count=0;
  // Select a full set initially
  int nselect = iter == 0 ? n : q;
  while(top > bot && count < nselect)
    {
      while(!(is_free(pidx[top]) 
	      || (is_upper_bound(pidx[top]) && y[pidx[top]] == +1)
	      || (is_lower_bound(pidx[top]) && y[pidx[top]] == -1)
	      ))
	--top;
      while(!(is_free(pidx[bot])
	      || (is_upper_bound(pidx[bot]) && y[pidx[bot]] == -1)
	      || (is_lower_bound(pidx[bot]) && y[pidx[bot]] == +1)
	      ))
	++bot;
      if(top > bot)
	{
	  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)
	    {
	      old_idx[i] = -1;
	      ++nnew;
	    }
	  else
	    {
	      for(int tt=0; tt<n_old; ++tt)
		{
		  if(old_work_set[tt] == t)
		    {
		      old_idx[i] = tt;
		      break;
		    }
		}
	    }
	  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(eps > 1e-10)
	eps /= 100;
      else
	{
	  info("Error: Unable to select a suitable working set!!!\n");
	  return 1;
	}
    }
  // Clean up
  delete[] yG;
  delete[] pidx;
  delete[] old_work_set;
  //  printf("done.\n");
  return 0;
}

void Solver_Parallel_SMO::init_working_set(int *work_set, int *not_work_set)