source: trunk/src/SVM.cc @ 37

// $Id: SVM.cc 37 2004-02-13 15:46:31Z peter $

#ifndef _THEP_SVM_
#define _THEP_SVM_

// System includes
#include <math.h>

// Thep C++ Tools
#include "SVM.h"
#include "matrix.h"
#include "vector.h"
#include "random_singleton.h"

using namespace thep_cpp_tools;
using namespace std;

SVM::SVM(const thep_gsl_api::matrix& kernel, 
   const thep_gsl_api::vector& target) : trained_(false), 
                 kernel_(kernel),
                 target_(target),
                 alpha_(target.size(),true,true),
                 bias_(0)
  
{
}

void SVM::train() //Should be done so one can choose to not train on all the samples
{
  thep_gsl_api::vector E = thep_gsl_api::vector(-target_);
  
  double upper_bound = pow(10., 32); 
  u_int count = 0;
  double alpha_new;
  random_singleton* rnd=random_singleton::get_instance();
  double u;
  double v;
  thep_gsl_api::vector dalpha;  
  
  // Stop criteria  
  bool stop_condition = false;
  while(!stop_condition)
    {
     count++;  
     dalpha = alpha_.mul_elements(target_);
     E = kernel_*dalpha-target_; //should be done in another way!!
       
     // Choosing a pair of variables to modify (Should be done more clever)
     u_long index1 = rnd->get_uniform_int(kernel_.cols());
     u_long index2 = rnd->get_uniform_int(kernel_.cols()-1);
     if (index2 >= index1) 
      index2++;
  
     //Updating the two variables
     if (target_[index1]!=target_[index2])
       {       
  if (alpha_[index2] > alpha_[index1] )
    {
    v = upper_bound;
    u = alpha_[index2] - alpha_[index1];
    }
  else
    {
     v = upper_bound - alpha_[index1] + alpha_[index2];
     u = 0;
    }
       }
     else
       {       
  if (alpha_[index2] + alpha_[index1] > upper_bound)
    {
    u = alpha_[index2] + alpha_[index1] - upper_bound;
    v = upper_bound;    
    }
  else
    {
     u = 0;
     v = alpha_[index1] + alpha_[index2];
    }
       }

     double k = kernel_.get(index1, index1) + kernel_.get(index2, index2) - 
       2*kernel_.get(index1, index2);
     alpha_new = alpha_[index2] + target_[index2]*
       (E[index1]-E[index2])/k;
     if (alpha_new > v)
       alpha_new = v;
     else if (alpha_new<u)
       alpha_new = u;
     
     alpha_[index1]+=target_[index1]*target_[index2]*(alpha_[index2]-alpha_new);
     alpha_[index2]=alpha_new;
     stop_condition = stop( target_, kernel_, alpha_);
     if (count>10000000){
      cerr << "SVM): " << "more than 10,000,000 epochs reached" << endl;
       exit(1);
     }
    }

  // Calculating the bias
  double min_output_positive = 10000;
  double max_output_negative = -10000;
  thep_gsl_api::vector output_unbiased = kernel_ * alpha_.mul_elements(target_);
  for (u_int i=0; i<target_.size(); i++){
   if (target_[i]==1){
     if (output_unbiased[i] < min_output_positive)
       min_output_positive = output_unbiased[i];
   }
   else if( output_unbiased[i] > max_output_negative )
     max_output_negative = output_unbiased[i];
  
  }
  cout << max_output_negative << endl;
  cout << min_output_positive << endl;
  
  bias_ = ( -min_output_positive - max_output_negative )/2;
  cout << "bias:" << bias_ << endl;

  trained_= true;
}

bool SVM::stop( const thep_gsl_api::vector& target_,
    const thep_gsl_api::matrix& kernel_,
    const thep_gsl_api::vector& alpha_)
{
  double min_output_positive = 10000;
  double max_output_negative = -10000;
  double epsilon = 0.000001; // used twice, should perhaps be two different values
  thep_gsl_api::vector output_unbiased = kernel_ * alpha_.mul_elements(target_);
  for (u_int i=0; i<target_.size(); i++){
   if (target_[i]==1){
     if (output_unbiased[i] < min_output_positive)
       min_output_positive = output_unbiased[i];
   }
   else if( output_unbiased[i] > max_output_negative )
     max_output_negative = output_unbiased[i];
   
  }
  
  if (min_output_positive - max_output_negative < 2-epsilon){
   return false;
  }
  else{
   double primal = alpha_.mul_elements(target_)
     *  ( alpha_.mul_elements(target_) );
   double gap = 2 * primal - alpha_.sum();
   if (gap/(primal+1)<epsilon)
     return true;
   else
     return false;
  }


}



#endif