source: trunk/c++_tools/classifier/SVM.cc @ 635

Last change on this file since 635 was 635, checked in by Peter, 16 years ago

closes #106 make_classifier is again taking DataLookup2D and Target making SupervisedClassifer? useful also in structure without SubsetGenerator?.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date ID
File size: 11.1 KB
Line 
1// $Id$
2
3#include <c++_tools/classifier/SVM.h>
4
5#include <c++_tools/classifier/DataLookup2D.h>
6#include <c++_tools/random/random.h>
7#include <c++_tools/statistics/Averager.h>
8#include <c++_tools/utility/matrix.h>
9#include <c++_tools/utility/vector.h>
10
11#include <algorithm>
12#include <cassert>
13#include <cctype>
14#include <cmath>
15#include <limits>
16#include <utility>
17#include <vector>
18
19
20namespace theplu {
21namespace classifier { 
22
23  SVM::SVM(const KernelLookup& kernel, const Target& target)
24    : SupervisedClassifier(target),
25      alpha_(target.size(),0),
26      bias_(0),
27      C_inverse_(0),
28      kernel_(&kernel),
29      margin_(0),
30      max_epochs_(100000),
31      output_(target.size(),0),
32      owner_(false),
33      sample_(target.size()),
34      trained_(false),
35      tolerance_(0.00000001)
36  {
37#ifndef NDEBUG
38    for (size_t i=0; i<alpha_.size(); i++) 
39      for (size_t j=0; j<alpha_.size(); j++) 
40        assert(kernel(i,j)==kernel(j,i));
41    for (size_t i=0; i<alpha_.size(); i++) 
42      for (size_t j=0; j<alpha_.size(); j++) 
43        assert((*kernel_)(i,j)==(*kernel_)(j,i));
44    for (size_t i = 0; i<kernel_->rows(); i++)
45      for (size_t j = 0; j<kernel_->columns(); j++)
46        if (std::isnan((*kernel_)(i,j)))
47          std::cerr << "SVM: Found nan in kernel: " << i << " " 
48                    << j << std::endl;
49#endif       
50
51  }
52
53  SVM::~SVM()
54  {
55    if (owner_)
56      delete kernel_;
57  }
58
59
60  void SVM::calculate_margin(void)
61  {
62    margin_ = 0;
63    for(size_t i = 0; i<alpha_.size(); ++i){
64      margin_ += alpha_(i)*target(i)*kernel_mod(i,i)*alpha_(i)*target(i);
65      for(size_t j = i+1; j<alpha_.size(); ++j)
66        margin_ += 2*alpha_(i)*target(i)*kernel_mod(i,j)*alpha_(j)*target(j);
67    }
68  }
69
70
71  SupervisedClassifier* SVM::make_classifier(const DataLookup2D& data,
72                                             const Target& target) const
73  {
74    const KernelLookup& kernel = 
75      dynamic_cast<const KernelLookup&>(data);
76
77    assert(data.rows()==data.columns());
78    assert(data.columns()==target.size());
79    SVM* sc;
80    sc = new SVM(kernel,target);
81
82
83    //Copy those variables possible to modify from outside
84    // Peter, in particular C
85    return sc;
86  }
87
88  void SVM::predict(const DataLookup2D& input, utility::matrix& prediction) const
89  {
90    // Peter, should check success of dynamic_cast
91    const KernelLookup input_kernel = dynamic_cast<const KernelLookup&>(input);
92
93    const KernelLookup* kernel_pointer;
94    kernel_pointer = &input_kernel;
95
96    assert(input.rows()==alpha_.size());
97    prediction = utility::matrix(2,input.columns(),0);
98    for (size_t i = 0; i<input.columns(); i++){
99      for (size_t j = 0; j<input.rows(); j++){
100        prediction(0,i) += target(j)*alpha_(j)*(*kernel_pointer)(j,i);
101        assert(target(j));
102      }
103      prediction(0,i) = margin_ * (prediction(0,i) + bias_);
104    }
105
106    for (size_t i = 0; i<prediction.columns(); i++)
107      prediction(1,i) = -prediction(0,i);
108   
109    assert(prediction(0,0));
110  }
111
112  double SVM::predict(const DataLookup1D& x) const
113  {
114    double y=0;
115    for (size_t i=0; i<alpha_.size(); i++)
116      y += alpha_(i)*target_(i)*kernel_->element(x,i);
117
118    return margin_*(y+bias_);
119  }
120
121  double SVM::predict(const DataLookupWeighted1D& x) const
122  {
123    double y=0;
124    for (size_t i=0; i<alpha_.size(); i++)
125      y += alpha_(i)*target_(i)*kernel_->element(x,i);
126
127    return margin_*(y+bias_);
128  }
129
130  bool SVM::train(void) 
131  {
132    // initializing variables for optimization
133    assert(target_.size()==kernel_->rows());
134    assert(target_.size()==alpha_.size());
135
136    sample_.init(alpha_,tolerance_);
137    utility::vector   E(target_.size(),0);
138    for (size_t i=0; i<E.size(); i++) {
139      for (size_t j=0; j<E.size(); j++) 
140        E(i) += kernel_mod(i,j)*target(j)*alpha_(j);
141      E(i)-=target(i);
142    }
143    assert(target_.size()==E.size());
144    assert(target_.size()==sample_.size());
145
146    unsigned long int epochs = 0;
147    double alpha_new2;
148    double alpha_new1;
149    double u;
150    double v;
151
152    // Training loop
153    while(choose(E)) {
154      bounds(u,v);       
155      double k = ( kernel_mod(sample_.value_first(), sample_.value_first()) + 
156                   kernel_mod(sample_.value_second(), sample_.value_second()) - 
157                   2*kernel_mod(sample_.value_first(), sample_.value_second()));
158     
159      double alpha_old1=alpha_(sample_.value_first());
160      double alpha_old2=alpha_(sample_.value_second());
161      alpha_new2 = ( alpha_(sample_.value_second()) + 
162                     target(sample_.value_second())*
163                     ( E(sample_.value_first())-E(sample_.value_second()) )/k );
164     
165      if (alpha_new2 > v)
166        alpha_new2 = v;
167      else if (alpha_new2<u)
168        alpha_new2 = u;
169     
170      // Updating the alphas
171      // if alpha is 'zero' make the sample a non-support vector
172      if (alpha_new2 < tolerance_){
173        sample_.nsv_second();
174      }
175      else{
176        sample_.sv_second();
177      }
178     
179     
180      alpha_new1 = (alpha_(sample_.value_first()) + 
181                    (target(sample_.value_first()) * 
182                     target(sample_.value_second()) * 
183                     (alpha_(sample_.value_second()) - alpha_new2) ));
184           
185      // if alpha is 'zero' make the sample a non-support vector
186      if (alpha_new1 < tolerance_){
187        sample_.nsv_first();
188      }
189      else
190        sample_.sv_first();
191     
192      alpha_(sample_.value_first()) = alpha_new1;
193      alpha_(sample_.value_second()) = alpha_new2;
194     
195      // update E vector
196      // Peter, perhaps one should only update SVs, but what happens in choose?
197      for (size_t i=0; i<E.size(); i++) {
198        E(i)+=( kernel_mod(i,sample_.value_first())*
199                target(sample_.value_first()) *
200                (alpha_new1-alpha_old1) );
201        E(i)+=( kernel_mod(i,sample_.value_second())*
202                target(sample_.value_second()) *
203                (alpha_new2-alpha_old2) );
204      }
205           
206      epochs++; 
207      if (epochs>max_epochs_){
208        std::cerr << "WARNING: SVM: maximal number of epochs reached.\n";
209        calculate_bias();
210        calculate_margin();
211        return false;
212      }
213    }
214    calculate_margin();
215    trained_ = calculate_bias();
216    return trained_;
217  }
218
219
220  bool SVM::choose(const theplu::utility::vector& E)
221  {
222    // First check for violation among SVs
223    // E should be the same for all SVs
224    // Choose that pair having largest violation/difference.
225    sample_.update_second(0);
226    sample_.update_first(0);
227    if (sample_.nof_sv()>1){
228
229      double max = E(sample_(0));
230      double min = max;
231      for (size_t i=1; i<sample_.nof_sv(); i++){ 
232        assert(alpha_(sample_(i))>tolerance_);
233        if (E(sample_(i)) > max){
234          max = E(sample_(i));
235          sample_.update_second(i);
236        }
237        else if (E(sample_(i))<min){
238          min = E(sample_(i));
239          sample_.update_first(i);
240        }
241      }
242      assert(alpha_(sample_.value_first())>tolerance_);
243      assert(alpha_(sample_.value_second())>tolerance_);
244
245      if (E(sample_.value_second()) - E(sample_.value_first()) > 2*tolerance_){
246        return true;
247      }
248     
249      // If no violation check among non-support vectors
250      sample_.shuffle();
251      for (size_t i=sample_.nof_sv(); i<sample_.size();i++){
252        if (target_.binary(sample_(i))){
253          if(E(sample_(i)) < E(sample_.value_first()) - 2*tolerance_){
254            sample_.update_second(i);
255            return true;
256          }
257        }
258        else{
259          if(E(sample_(i)) > E(sample_.value_second()) + 2*tolerance_){
260            sample_.update_first(i);
261            return true;
262          }
263        }
264      }
265    }
266
267    // if no support vectors - special case
268    else{
269      // to avoid getting stuck we shuffle
270      sample_.shuffle();
271      for (size_t i=0; i<sample_.size(); i++) {
272        if (target(sample_(i))==1){
273          for (size_t j=0; j<sample_.size(); j++) {
274            if ( target(sample_(j))==-1 && 
275                 E(sample_(i)) < E(sample_(j))+2*tolerance_ ){
276              sample_.update_first(i);
277              sample_.update_second(j);
278              return true;
279            }
280          }
281        }
282      }
283    }
284   
285    // If there is no violation then we should stop training
286    return false;
287
288  }
289 
290 
291  void SVM::bounds( double& u, double& v) const
292  {
293    if (target(sample_.value_first())!=target(sample_.value_second())) {
294      if (alpha_(sample_.value_second()) > alpha_(sample_.value_first())) {
295        v = std::numeric_limits<double>::max();
296        u = alpha_(sample_.value_second()) - alpha_(sample_.value_first());
297      }
298      else {
299        v = (std::numeric_limits<double>::max() - 
300             alpha_(sample_.value_first()) + 
301             alpha_(sample_.value_second()));
302        u = 0;
303      }
304    }
305    else {       
306      if (alpha_(sample_.value_second()) + alpha_(sample_.value_first()) > 
307           std::numeric_limits<double>::max()) {
308        u = (alpha_(sample_.value_second()) + alpha_(sample_.value_first()) - 
309              std::numeric_limits<double>::max());
310        v =  std::numeric_limits<double>::max();   
311      }
312      else {
313        u = 0;
314        v = alpha_(sample_.value_first()) + alpha_(sample_.value_second());
315      }
316    }
317  }
318 
319  bool SVM::calculate_bias(void)
320  {
321
322    // calculating output without bias
323    for (size_t i=0; i<output_.size(); i++) {
324      output_(i)=0;
325      for (size_t j=0; j<output_.size(); j++) 
326        output_(i)+=alpha_(j)*target(j) * (*kernel_)(i,j);
327    }
328
329    if (!sample_.nof_sv()){
330      std::cerr << "SVM::train() error: " 
331                << "Cannot calculate bias because there is no support vector" 
332                << std::endl;
333      return false;
334    }
335
336    // For samples with alpha>0, we have: target*output=1-alpha/C
337    bias_=0;
338    for (size_t i=0; i<sample_.nof_sv(); i++) 
339      bias_+= ( target(sample_(i)) * (1-alpha_(sample_(i))*C_inverse_) - 
340                output_(sample_(i)) );
341    bias_=bias_/sample_.nof_sv();
342    for (size_t i=0; i<output_.size(); i++) 
343      output_(i) += bias_;
344     
345    return true;
346  }
347
348  void SVM::set_C(const double C)
349  {
350    C_inverse_ = 1/C;
351  }
352
353  Index::Index(void)
354    : nof_sv_(0), vec_(std::vector<size_t>(0))
355  {
356  }
357
358  Index::Index(const size_t n)
359    : nof_sv_(0), vec_(std::vector<size_t>(n))
360  {
361    for (size_t i=0; i<vec_.size(); i++)
362      vec_[i]=i;
363  }
364
365  void Index::init(const utility::vector& alpha, const double tol)
366  {
367    nof_sv_=0;
368    size_t nof_nsv=0;
369    for (size_t i=0; i<alpha.size(); i++) 
370      if (alpha(i)<tol){
371        nof_nsv++;
372        vec_[vec_.size()-nof_nsv]=i;
373      }
374      else{
375        vec_[nof_sv_]=i;
376        nof_sv_++;
377      }
378    assert(nof_sv_+nof_nsv==vec_.size());
379
380  }
381
382  void Index::sv_first(void)
383  {
384    // if already sv, do nothing
385    if (index_first_<nof_sv())
386      return;
387
388    // swap elements
389    if(index_second_==nof_sv_){
390      index_second_=index_first_;
391    }
392    vec_[index_first_]=vec_[nof_sv_];
393    vec_[nof_sv_]=value_first_;
394    index_first_ = nof_sv_;
395
396    nof_sv_++;
397
398  }
399
400  void Index::sv_second(void)
401  {
402    // if already sv, do nothing
403    if (index_second_<nof_sv())
404      return;
405
406    // swap elements
407    if(index_first_==nof_sv_){
408      index_first_=index_second_;
409    }
410
411    vec_[index_second_]=vec_[nof_sv_];
412    vec_[nof_sv_]=value_second_;
413    index_second_=nof_sv_;
414
415    nof_sv_++;
416  }
417
418  void Index::nsv_first(void)
419  {
420    // if already nsv, do nothing
421    if ( !(index_first_<nof_sv()) )
422      return;
423   
424    if(index_second_==nof_sv_-1)
425      index_second_=index_first_;
426    vec_[index_first_]=vec_[nof_sv_-1];
427    vec_[nof_sv_-1]=value_first_;
428    index_first_=nof_sv_-1;
429   
430    nof_sv_--;
431  }
432
433  void Index::nsv_second(void)
434  {
435    // if already nsv, do nothing
436    if ( !(index_second_<nof_sv()) )
437      return;
438
439    if(index_first_==nof_sv_-1)
440      index_first_=index_second_;
441    vec_[index_second_]=vec_[nof_sv_-1];
442    vec_[nof_sv_-1]=value_second_;
443    index_second_ = nof_sv_-1;
444   
445    nof_sv_--;
446  }
447
448
449  void Index::shuffle(void)
450  {
451    random::DiscreteUniform a;
452    random_shuffle(vec_.begin()+nof_sv_, vec_.end(), a); 
453  }
454
455  void Index::update_first(const size_t i)
456  {
457    assert(i<size());
458    index_first_=i;
459    value_first_=vec_[i];
460  }
461
462  void Index::update_second(const size_t i)
463  {
464    assert(i<size());
465    index_second_=i;
466    value_second_=vec_[i];
467  }
468
469}} // of namespace classifier and namespace theplu
Note: See TracBrowser for help on using the repository browser.