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

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

re-introduced prediction in SVM

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