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

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

added test for target
redesign crossSplitter
added two class function in Target

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