// 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_DOUBLESolver_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)