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

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

forked tests and fixed bugs

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