source: trunk/lib/classifier/SVM.cc @ 544

// $Id: SVM.cc 544 2006-03-05 17:33:33Z peter $

#include <c++_tools/classifier/SVM.h>

#include <c++_tools/classifier/DataLookup2D.h>
#include <c++_tools/gslapi/matrix.h>
#include <c++_tools/gslapi/vector.h>
#include <c++_tools/statistics/Averager.h>
#include <c++_tools/statistics/Score.h>
#include <c++_tools/random/random.h>

#include <algorithm>
#include <cassert>
#include <cmath>
#include <limits>
#include <utility>
#include <vector>


namespace theplu {
namespace classifier {  

  SVM::SVM(const KernelLookup& kernel, const Target& target)
    : SupervisedClassifier(target),
      alpha_(target.size(),0),
      bias_(0),
      C_inverse_(0),
      kernel_(&kernel),
      max_epochs_(10000000),
      output_(target.size(),0),
      sample_(target.size()),
      trained_(false),
      tolerance_(0.00000001)
  {
  }

  SVM::SVM(const KernelLookup& kernel, const Target& target, 
           statistics::Score& score, const size_t nof_inputs)
    : SupervisedClassifier(target, &score, nof_inputs),
      alpha_(target.size(),0),
      bias_(0),
      C_inverse_(0),
      kernel_(&kernel),
      max_epochs_(10000000),
      output_(target.size(),0),
      sample_(target.size()),
      trained_(false),
      tolerance_(0.00000001)
  {
  }

  

  SupervisedClassifier* SVM::make_classifier(const DataLookup2D& data, 
                                             const Target& target) const
  {
    // Peter, should check success of dynamic_cast
    const KernelLookup& tmp = dynamic_cast<const KernelLookup&>(data);
    SVM* sc;
    if (score_)
      sc = new SVM(tmp,target,*score_,nof_inputs_);
    else
      sc = new SVM(tmp,target);

    //Copy those variables possible to modify from outside
    return sc;
  }

  void SVM::predict(const DataLookup2D& input, gslapi::matrix& prediction) const
  {
    assert(input.rows()==alpha_.size());
    prediction=gslapi::matrix(2,input.columns(),0);
    for (size_t i = 0; i<input.columns(); i++){
      for (size_t j = 0; i<input.rows(); i++)
        prediction(0,i) += target(j)*alpha_(i)*input(j,i);
      prediction(0,i) += bias_;
    }

    for (size_t i = 0; i<prediction.columns(); i++)
      prediction(1,i) = -prediction(0,i);
  }

  double SVM::predict(const DataLookup1D& x) const
  {
    double y=0;
    for (size_t i=0; i<alpha_.size(); i++)
      y += alpha_(i)*target_(i)*kernel_->element(x,i);

    return y+bias_;
  }

  double SVM::predict(const DataLookup1D& x, const DataLookup1D& w) const
  {
    double y=0;
    for (size_t i=0; i<alpha_.size(); i++)
      y += alpha_(i)*target_(i)*kernel_->element(x,w,i);

    return y+bias_;
  }

  bool SVM::train(void) 
  {
    // initializing variables for optimization
    assert(target_.size()==kernel_->rows());
    assert(target_.size()==alpha_.size());

    sample_.init(alpha_,tolerance_);
    gslapi::vector  E(target_.size(),0);
    for (size_t i=0; i<E.size(); i++) {
      for (size_t j=0; j<E.size(); j++) 
        E(i) += kernel_mod(i,j)*target(j)*alpha_(j);
      E(i)=E(i)-target(i);
    }
    assert(target_.size()==E.size());
    assert(target_.size()==sample_.n());

    unsigned long int epochs = 0;
    double alpha_new2;
    double alpha_new1;
    double u;
    double v;

    // Training loop
    while(choose(E)) {
      bounds(u,v);        
      double k = ( kernel_mod(sample_.value_first(), sample_.value_first()) + 
                   kernel_mod(sample_.value_second(), sample_.value_second()) - 
                   2*kernel_mod(sample_.value_first(), sample_.value_second()));
      
      double alpha_old1=alpha_(sample_.value_first());
      double alpha_old2=alpha_(sample_.value_second());

      alpha_new2 = ( alpha_(sample_.value_second()) + 
                     target(sample_.value_second())*
                     ( E(sample_.value_first())-E(sample_.value_second()) )/k );
      
      if (alpha_new2 > v)
        alpha_new2 = v;
      else if (alpha_new2<u)
        alpha_new2 = u;
      
      
      // Updating the alphas
      // if alpha is 'zero' make the sample a non-support vector
      if (alpha_new2 < tolerance_){
        sample_.nsv_second();
      }
      else{
        sample_.sv_second();
      }
      
      
      alpha_new1 = (alpha_(sample_.value_first()) + 
                    (target(sample_.value_first()) * 
                     target(sample_.value_second()) * 
                     (alpha_(sample_.value_second()) - alpha_new2) ));
            
      // if alpha is 'zero' make the sample a non-support vector
      if (alpha_new1 < tolerance_){
        sample_.nsv_first();
      }
      else
        sample_.sv_first();
      
      alpha_(sample_.value_first()) = alpha_new1;
      alpha_(sample_.value_second()) = alpha_new2;
      
      // update E vector
      // Peter, perhaps one should only update SVs, but what happens in choose?
      for (size_t i=0; i<E.size(); i++) {
        E(i)+=( kernel_mod(i,sample_.value_first())*
                target(sample_.value_first()) *
                (alpha_new1-alpha_old1) );
        E(i)+=( kernel_mod(i,sample_.value_second())*
                target(sample_.value_second()) *
                (alpha_new2-alpha_old2) );
      }
            
      epochs++;  
      if (epochs>max_epochs_){
        std::cerr << "WARNING: SVM: maximal number of epochs reached.\n";
        return false;
      }
    }
    
    trained_ = calculate_bias();
    return trained_;
  }


  bool SVM::choose(const theplu::gslapi::vector& E)
  {
    // First check for violation among SVs
    // E should be the same for all SVs
    // Choose that pair having largest violation/difference.
    sample_.update_second(0);
    sample_.update_first(0);
    if (sample_.nof_sv()>1){

      double max = E(sample_(0));
      double min = max;
      for (size_t i=1; i<sample_.nof_sv(); i++){ 
        assert(alpha_(sample_(i))>tolerance_);
        if (E(sample_(i)) > max){
          max = E(sample_(i));
          sample_.update_second(i);
        }
        else if (E(sample_(i))<min){
          min = E(sample_(i));
          sample_.update_first(i);
        }
      }
      assert(alpha_(sample_.value_first())>tolerance_);
      assert(alpha_(sample_.value_second())>tolerance_);

      
      if (E(sample_.value_second()) - E(sample_.value_first()) > 2*tolerance_){
        return true;
      }
      
      // If no violation check among non-support vectors
      sample_.shuffle();
      for (size_t i=sample_.nof_sv(); i<sample_.n();i++){
        if (target_.binary(sample_(i))){
          if(E(sample_(i)) < E(sample_.value_first()) - 2*tolerance_){
            sample_.update_second(i);
            return true;
          }
        }
        else{
          if(E(sample_(i)) > E(sample_.value_second()) + 2*tolerance_){
            sample_.update_first(i);
            return true;
          }
        }
      }
    }

    // if no support vectors - special case
    else{
      for (size_t i=0; i<sample_.n(); i++) {
        if (target_.binary(sample_(i))){
          for (size_t j=0; j<sample_.n(); j++) {
            if ( !target_.binary(sample_(j)) && 
                 E(sample_(i)) < E(sample_(j))+2*tolerance_ ){
              sample_.update_first(i);
              sample_.update_second(j);
              return true;
            }
          }
        }
      }
    }
    
    // If there is no violation then we should stop training
    return false;

  }
  
  
  void SVM::bounds( double& u, double& v) const
  {
    if (target(sample_.value_first())!=target(sample_.value_second())) {
      if (alpha_(sample_.value_second()) > alpha_(sample_.value_first())) {
        v = std::numeric_limits<double>::max();
        u = alpha_(sample_.value_second()) - alpha_(sample_.value_first());
      }
      else {
        v = (std::numeric_limits<double>::max() - 
             alpha_(sample_.value_first()) + 
             alpha_(sample_.value_second()));
        u = 0;
      }
    }
    else {       
      if (alpha_(sample_.value_second()) + alpha_(sample_.value_first()) > 
           std::numeric_limits<double>::max()) {
        u = (alpha_(sample_.value_second()) + alpha_(sample_.value_first()) - 
              std::numeric_limits<double>::max());
        v =  std::numeric_limits<double>::max();    
      }
      else {
        u = 0;
        v = alpha_(sample_.value_first()) + alpha_(sample_.value_second());
      }
    }
  }
  
  bool SVM::calculate_bias(void)
  {

    // calculating output without bias
    for (size_t i=0; i<output_.size(); i++) {
      output_(i)=0;
      for (size_t j=0; j<output_.size(); j++) 
        output_(i)+=alpha_(j)*target(j) * (*kernel_)(i,j);
    }

    if (!sample_.nof_sv()){
      std::cerr << "SVM::train() error: " 
                << "Cannot calculate bias because there is no support vector" 
                << std::endl;
      return false;
    }

    // For samples with alpha>0, we have: target*output=1-alpha/C
    bias_=0;
    for (size_t i=0; i<sample_.nof_sv(); i++) 
      bias_+= ( target(sample_(i)) * (1-alpha_(sample_(i))*C_inverse_) - 
                output_(sample_(i)) );
    bias_=bias_/sample_.nof_sv();
    for (size_t i=0; i<output_.size(); i++) 
      output_(i) += bias_;
      
    return true;
  }

  Index::Index(void)
    : nof_sv_(0), vec_(std::vector<size_t>(0))
  {
  }

  Index::Index(const size_t n)
    : nof_sv_(0), vec_(std::vector<size_t>(n))
  {
    for (size_t i=0; i<vec_.size(); i++)
      vec_[i]=i;
  }

  void Index::init(const gslapi::vector& alpha, const double tol)
  {
    nof_sv_=0;
    size_t nof_nsv=0;
    for (size_t i=0; i<alpha.size(); i++) 
      if (alpha(i)<tol){
        nof_nsv++;
        vec_[vec_.size()-nof_nsv]=i;
      }
      else{
        vec_[nof_sv_]=i;
        nof_sv_++;
      }
    assert(nof_sv_+nof_nsv==vec_.size());

  }

  void Index::sv_first(void)
  {
    // if already sv, do nothing
    if (index_first_<nof_sv())
      return;

    // swap elements
    if(index_second_==nof_sv_){
      index_second_=index_first_;
    }
    vec_[index_first_]=vec_[nof_sv_];
    vec_[nof_sv_]=value_first_;
    index_first_ = nof_sv_;

    nof_sv_++;

  }

  void Index::sv_second(void)
  {
    // if already sv, do nothing
    if (index_second_<nof_sv())
      return;

    // swap elements
    if(index_first_==nof_sv_){
      index_first_=index_second_;
    }

    vec_[index_second_]=vec_[nof_sv_];
    vec_[nof_sv_]=value_second_;
    index_second_=nof_sv_;

    nof_sv_++;
  }

  void Index::nsv_first(void)
  {
    // if already nsv, do nothing
    if ( !(index_first_<nof_sv()) )
      return;
    
    if(index_second_==nof_sv_-1)
      index_second_=index_first_;
    vec_[index_first_]=vec_[nof_sv_-1];
    vec_[nof_sv_-1]=value_first_;
    index_first_=nof_sv_-1;
    
    nof_sv_--;
  }

  void Index::nsv_second(void)
  {
    // if already nsv, do nothing
    if ( !(index_second_<nof_sv()) )
      return;

    if(index_first_==nof_sv_-1)
      index_first_=index_second_;
    vec_[index_second_]=vec_[nof_sv_-1];
    vec_[nof_sv_-1]=value_second_;
    index_second_ = nof_sv_-1;
    
    nof_sv_--;
  }


  void Index::shuffle(void)
  {
    random::DiscreteUniform a;
    random_shuffle(vec_.begin()+nof_sv_, vec_.end(), a); 
  }

  void Index::update_first(const size_t i)
  {
    assert(i<n());
    index_first_=i;
    value_first_=vec_[i];
  }

  void Index::update_second(const size_t i)
  {
    assert(i<n());
    index_second_=i;
    value_second_=vec_[i];
  }

}} // of namespace classifier and namespace theplu