source: trunk/src/SVM.cc @ 115

// $Id: SVM.cc 115 2004-07-19 14:38:43Z peter $

// System includes
#include <cmath>
#include <utility>
#include <vector>

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



namespace theplu {
namespace cpptools {  

  SVM::SVM(const gslapi::matrix& kernel, const gslapi::vector& target, 
           const std::vector<size_t> train_set) 
  : alpha_(target.size(),0),
    bias_(0),
    c_(0),
    kernel_(kernel),
    target_(target),
    trained_(false),
    train_set_(train_set),
    tolerance_(0.000000000001)
       
  {
    if (!train_set_.size())
      for (size_t i=0; i<target_.size(); i++)
        train_set_.push_back(i);  
  }


  std::pair<size_t,size_t> 
  SVM::choose(const theplu::gslapi::vector& E,  
              const theplu::gslapi::vector& target_train,
              const theplu::gslapi::vector& alpha_train,
              bool& stop_condition)
  {
    std::pair<size_t, size_t> index(0,0);     
    
    // First check for violation among support vectors
    double max = -10000000;
    double min = 10000000;
    // Peter, avoid looping over all alphas
    for (size_t i=0; i<E.size(); i++){ 
      if (alpha_train(i)){
        if (E(i)>max){
          max = E(i);
          index.second = i;
        }
        if (E(i)<min){
          min = E(i);
          index.first = i;
        }
      }
    }
    //index = E.minmax_index(support_vectors); 
    if (E(index.second) - E(index.first) > 2*tolerance_)
      return index;
    
    // If no violation check among non-support vectors
    std::vector<size_t> tmp;
    tmp.resize(E.size());
    for (size_t i=0; i<E.size(); i++)
      tmp[i]=i;
    my_uniform_rng a;
    random_shuffle(tmp.begin(), tmp.end(), a);
    for (size_t i=0; i<tmp.size(); i++){
      if (target_train(tmp[i])==1){
        if ( E(tmp[i]) < E(index.first) - 2*tolerance_){
          index.second = tmp[i];
          return index;
        }
      }
      else if (target_train(tmp[i])==-1){
        if ( E(tmp[i]) > E(index.second) + 2*tolerance_){
          index.first = tmp[i];
          return index;
        }
      }
    }
    // If there is no violation then we should stop
    stop_condition = true;
    return index;
  }
  
  void SVM::train() 
  {
  using namespace std;

  // initializing variables for optimization
  gslapi::matrix kernel_train(train_set_.size(),train_set_.size());
  gslapi::vector alpha_train(train_set_.size());
  gslapi::vector target_train(train_set_.size());
  for (unsigned int i=0; i<train_set_.size(); i++) {
    alpha_train(i) = alpha_(train_set_[i]);
    target_train(i) = target_(train_set_[i]);
    for (unsigned int j=0; j<train_set_.size(); j++) 
      kernel_train(i,j) = kernel_(train_set_[i],train_set_[j]); 
  }

  // Modify the diagonal of the kernel matrix (soft margin)
  double d;
  if (c_>0)
    d = 1/c_;
  else if (c_==0)
    d = 0;
  else {
    cerr << "SVM): " << "C-parameter cannot be a negative number" << endl;
    exit(1);
  }
  for (unsigned int i=0; i<kernel_train.columns(); i++)
    kernel_train(i,i) += d; 
  
  gslapi::vector  E = ( kernel_train*target_train.mul_elements(alpha_train) 
                        - target_train ) ;

  double upper_bound = pow(10., 32.); //Peter, should be static...
  u_int epochs = 0;
  double alpha_new2;
  double alpha_new1;
  double u;
  double v;
  bool stop_condition = false;

  // Training loop
  while(!stop_condition) {
    pair<size_t, size_t> index = choose(E, target_train, 
                                        alpha_train,stop_condition); 
    if (target_train[index.first]!=target_train[index.second]) {       
      if (alpha_train[index.second] > alpha_train[index.first]) {
        v = upper_bound;
        u = alpha_train[index.second] - alpha_train[index.first];
      }
      else {
        v = upper_bound - alpha_train[index.first] + alpha_train[index.second];
        u = 0;
      }
    }
    else {       
      if (alpha_train[index.second] + alpha_train[index.first] > upper_bound) {
        u = alpha_train[index.second] + alpha_train[index.first] - upper_bound;
        v = upper_bound;    
      }
      else {
        u = 0;
        v = alpha_train[index.first] + alpha_train[index.second];
      }
    }
    
    double k = ( kernel_train(index.first, index.first) + 
                 kernel_train(index.second, index.second) - 
                 2*kernel_train(index.first, index.second));
    alpha_new2 = (alpha_train[index.second] + target_train[index.second]*
                 (E[index.first]-E[index.second])/k );
    
    if (alpha_new2 > v)
      alpha_new2 = v;
    else if (alpha_new2<u)
      alpha_new2 = u;
    
    // Updating the alphas and some more
    if (alpha_new2 < tolerance_)
      alpha_new2=0;
    
    alpha_new1 = (alpha_train[index.first] + 
                  (target_train[index.first] * target_train[index.second] * 
                   (alpha_train[index.second] - alpha_new2)) );
    if (alpha_new1 < tolerance_){
      alpha_new1=0;
    }

    alpha_train[index.first] = alpha_new1;
    alpha_train[index.second] = alpha_new2;
    E = kernel_train*target_train.mul_elements(alpha_train) - target_train;
    gslapi::vector ones(alpha_train.size(),1);
    
    epochs++;  
    if (epochs>1000000){
      cerr << "WARNING: SVM: " << "1,000,000 epochs reached" << endl;
      exit(1);
    }
  }

  // Calculating the bias from support vectors. E = output - bias - target
  Averager bias;
  for (size_t i=0; i<E.size(); i++)
    if (alpha_train(i))
      bias.add(-E(i));
  bias_ = bias.mean();
  
  trained_= true;
  alpha_.set_all(0);
  for (unsigned int i=0; i<train_set_.size(); i++) 
    alpha_(train_set_[i]) = alpha_train(i);
    
  }

  
   
    

}} // of namespace cpptools and namespace theplu