source: trunk/src/SVM.cc @ 227

// $Id: SVM.cc 227 2005-02-01 00:52:27Z jari $

#include "SVM.h"

// System includes
#include <cmath>
#include <limits.h>
#include <utility>
#include <vector>

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



namespace theplu {
namespace cpptools {  

  SVM::SVM(const Kernel& kernel, 
           const gslapi::vector& target, 
           const std::vector<size_t>& train_set) 
  : alpha_(target.size(),0),
    bias_(0),
    c_(0),
    kernel_(kernel),
    max_epochs_(10000000),
    target_(target),
    trained_(false),
    train_set_(train_set),
    tolerance_(0.000001)
       
  {
    if (max_epochs_>ULONG_MAX)
      max_epochs_=ULONG_MAX;
    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;
  }
  
  

  bool 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++) {
      target_train(i) = target_(train_set_[i]);
      for (unsigned int j=0; j<train_set_.size(); j++) 
        kernel_train(i,j) = kernel_.get()(train_set_[i],train_set_[j]); 
    }
    // Modify the diagonal of the kernel matrix (soft margin)
    if (c_>0){
      double d = 1/c_;
      for (unsigned int i=0; i<kernel_train.columns(); i++)
        kernel_train(i,i) += d; 
    }
  
    gslapi::vector  E = ( kernel_train*target_train.mul(alpha_train) 
                          - target_train ) ;

    double upper_bound = pow(10., 32); //Peter, should be static...
    unsigned long 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); 
      // Defining upper and lower bounds for alpha
      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(alpha_train) - target_train;
      gslapi::vector ones(alpha_train.size(),1);
    
      epochs++;  
      if (epochs>max_epochs_){
        return false;
      }
    }

    // Calculating the bias from support vectors. E = output - bias - target
    statistics::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);
    
    // training performed succesfully
    return true;
    
  }

  
   
    

}} // of namespace cpptools and namespace theplu