Changeset 323
- Timestamp:
- May 26, 2005, 10:54:30 PM (18 years ago)
- Location:
- trunk
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/lib/svm/SVM.cc
r307 r323 10 10 #include <c++_tools/utility/random_singleton.h> 11 11 12 #include <iostream> 13 #include <algorithm> 12 14 #include <cmath> 13 #include <limits .h>15 #include <limits> 14 16 #include <utility> 15 17 #include <vector> … … 19 21 namespace svm { 20 22 23 Index::Index(void) 24 : nof_sv_(0), vec_(std::vector<size_t>(0)) 25 { 26 } 27 28 Index::Index(const size_t n) 29 : nof_sv_(0), vec_(std::vector<size_t>(n)) 30 { 31 for (size_t i=0; i<vec_.size(); i++) 32 vec_[i]=i; 33 } 34 35 void Index::init(const gslapi::vector& alpha, const double tol) 36 { 37 nof_sv_=0; 38 size_t nof_nsv=0; 39 for (size_t i=0; i<alpha.size(); i++) 40 if (alpha(i)<tol){ 41 nof_nsv++; 42 vec_[vec_.size()-nof_nsv]=i; 43 } 44 else{ 45 vec_[nof_sv_]=i; 46 nof_sv_++; 47 } 48 assert(nof_sv_+nof_nsv==vec_.size()); 49 50 } 51 52 void Index::sv_first(void) 53 { 54 // if already sv, do nothing 55 if (index_first_<nof_sv()) 56 return; 57 58 // swap elements 59 if(index_second_==nof_sv_){ 60 index_second_=index_first_; 61 } 62 vec_[index_first_]=vec_[nof_sv_]; 63 vec_[nof_sv_]=value_first_; 64 index_first_ = nof_sv_; 65 66 nof_sv_++; 67 68 } 69 70 void Index::sv_second(void) 71 { 72 // if already sv, do nothing 73 if (index_second_<nof_sv()) 74 return; 75 76 // swap elements 77 if(index_first_==nof_sv_){ 78 index_first_=index_second_; 79 } 80 81 vec_[index_second_]=vec_[nof_sv_]; 82 vec_[nof_sv_]=value_second_; 83 index_second_=nof_sv_; 84 85 nof_sv_++; 86 } 87 88 void Index::nsv_first(void) 89 { 90 // if already nsv, do nothing 91 if ( !(index_first_<nof_sv()) ) 92 return; 93 94 if(index_second_==nof_sv_-1) 95 index_second_=index_first_; 96 vec_[index_first_]=vec_[nof_sv_-1]; 97 vec_[nof_sv_-1]=value_first_; 98 index_first_=nof_sv_-1; 99 100 nof_sv_--; 101 } 102 103 void Index::nsv_second(void) 104 { 105 // if already nsv, do nothing 106 if ( !(index_second_<nof_sv()) ) 107 return; 108 109 if(index_first_==nof_sv_-1) 110 index_first_=index_second_; 111 vec_[index_second_]=vec_[nof_sv_-1]; 112 vec_[nof_sv_-1]=value_second_; 113 index_second_ = nof_sv_-1; 114 115 nof_sv_--; 116 } 117 118 119 void Index::shuffle(void) 120 { 121 utility::my_uniform_rng a; 122 random_shuffle(vec_.begin()+nof_sv_, vec_.end(), a); 123 } 124 125 void Index::update_first(const size_t i) 126 { 127 assert(i<n()); 128 index_first_=i; 129 value_first_=vec_[i]; 130 } 131 132 void Index::update_second(const size_t i) 133 { 134 assert(i<n()); 135 index_second_=i; 136 value_second_=vec_[i]; 137 } 138 21 139 SVM::SVM(const Kernel_MEV& kernel, 22 const gslapi::vector& target ,23 const std::vector<size_t>& train_set)140 const gslapi::vector& target) 141 24 142 : alpha_(target.size(),0), 25 143 bias_(0), 26 c_(0),144 C_inverse_(0), 27 145 kernel_(kernel), 28 146 max_epochs_(10000000), 147 output_(target.size(),0), 148 sample_(target.size()), 29 149 target_(target), 30 150 trained_(false), 31 train_set_(train_set), 32 tolerance_(0.000001) 151 tolerance_(0.00000001) 33 152 34 153 { 35 154 if (max_epochs_>ULONG_MAX) 36 155 max_epochs_=ULONG_MAX; 37 if (!train_set_.size()) 38 for (size_t i=0; i<target_.size(); i++) 39 train_set_.push_back(i); 40 } 41 42 43 44 std::pair<size_t,size_t> 45 SVM::choose(const theplu::gslapi::vector& E, 46 const theplu::gslapi::vector& target_train, 47 const theplu::gslapi::vector& alpha_train, 48 bool& stop_condition) 49 { 50 std::pair<size_t, size_t> index(0,0); 51 52 // First check for violation among support vectors 53 double max = -10000000; 54 double min = 10000000; 55 // Peter, avoid looping over all alphas 56 for (size_t i=0; i<E.size(); i++){ 57 if (alpha_train(i)){ 58 if (E(i)>max){ 59 max = E(i); 60 index.second = i; 61 } 62 if (E(i)<min){ 63 min = E(i); 64 index.first = i; 65 } 66 } 67 } 68 //index = E.minmax_index(support_vectors); 69 if (E(index.second) - E(index.first) > 2*tolerance_) 70 return index; 71 72 // If no violation check among non-support vectors 73 std::vector<size_t> tmp; 74 tmp.resize(E.size()); 75 for (size_t i=0; i<E.size(); i++) 76 tmp[i]=i; 77 utility::my_uniform_rng a; 78 random_shuffle(tmp.begin(), tmp.end(), a); 79 for (size_t i=0; i<tmp.size(); i++){ 80 if (target_train(tmp[i])==1){ 81 if ( E(tmp[i]) < E(index.first) - 2*tolerance_){ 82 index.second = tmp[i]; 83 return index; 84 } 85 } 86 else if (target_train(tmp[i])==-1){ 87 if ( E(tmp[i]) > E(index.second) + 2*tolerance_){ 88 index.first = tmp[i]; 89 return index; 90 } 91 } 92 } 93 // If there is no violation then we should stop 94 stop_condition = true; 95 return index; 96 } 97 98 99 100 bool SVM::train() 101 { 102 using namespace std; 156 } 157 158 159 bool SVM::train(void) 160 { 103 161 // initializing variables for optimization 104 gslapi::matrix kernel_train(train_set_.size(),train_set_.size()); 105 gslapi::vector alpha_train(train_set_.size()); 106 gslapi::vector target_train(train_set_.size()); 107 108 // copy submatrix to train on. Peter, might be unnecessary. 109 for (unsigned int i=0; i<train_set_.size(); i++) { 110 target_train(i) = target_(train_set_[i]); 111 for (unsigned int j=0; j<train_set_.size(); j++) 112 kernel_train(i,j) = kernel_(train_set_[i],train_set_[j]); 113 } 114 // Modify the diagonal of the kernel matrix (soft margin) 115 if (c_>0){ 116 double d = 1/c_; 117 for (unsigned int i=0; i<kernel_train.columns(); i++) 118 kernel_train(i,i) += d; 119 } 120 121 gslapi::vector E = ( kernel_train*target_train.mul(alpha_train) 122 - target_train ) ; 123 124 double upper_bound = pow(10., 32); //Peter, should be static... 162 assert(target_.size()==kernel_.size()); 163 assert(target_.size()==alpha_.size()); 164 165 sample_.init(alpha_,tolerance_); 166 gslapi::vector E(target_.size()); 167 for (size_t i=0; i<E.size(); i++) { 168 E(i)=0; 169 for (size_t j=0; j<E.size(); j++) 170 E(i) += kernel_mod(i,j)*target_(j)*alpha_(j); 171 E(i)=E(i)-target_(i); 172 } 173 assert(target_.size()==E.size()); 174 assert(target_.size()==sample_.n()); 175 125 176 unsigned long int epochs = 0; 126 177 double alpha_new2; … … 128 179 double u; 129 180 double v; 130 bool stop_condition = false; 181 131 182 // Training loop 132 while(!stop_condition) { 133 pair<size_t, size_t> index = choose(E, target_train, 134 alpha_train,stop_condition); 135 // Defining upper and lower bounds for alpha 136 if (target_train[index.first]!=target_train[index.second]) { 137 if (alpha_train[index.second] > alpha_train[index.first]) { 138 v = upper_bound; 139 u = alpha_train[index.second] - alpha_train[index.first]; 140 } 141 else { 142 v = (upper_bound - alpha_train[index.first] + 143 alpha_train[index.second]); 144 u = 0; 145 } 146 } 147 else { 148 if (alpha_train[index.second] + alpha_train[index.first] > 149 upper_bound) { 150 u = (alpha_train[index.second] + alpha_train[index.first] - 151 upper_bound); 152 v = upper_bound; 153 } 154 else { 155 u = 0; 156 v = alpha_train[index.first] + alpha_train[index.second]; 157 } 158 } 159 160 double k = ( kernel_train(index.first, index.first) + 161 kernel_train(index.second, index.second) - 162 2*kernel_train(index.first, index.second)); 163 164 alpha_new2 = (alpha_train[index.second] + target_train[index.second]* 165 (E[index.first]-E[index.second])/k ); 166 183 while(choose(E)) { 184 bounds(u,v); 185 double k = ( kernel_mod(sample_.value_first(), sample_.value_first()) + 186 kernel_mod(sample_.value_second(), sample_.value_second()) - 187 2*kernel_mod(sample_.value_first(), sample_.value_second())); 188 189 double alpha_old1=alpha_(sample_.value_first()); 190 double alpha_old2=alpha_(sample_.value_second()); 191 192 alpha_new2 = ( alpha_(sample_.value_second()) + 193 target_(sample_.value_second())* 194 ( E(sample_.value_first())-E(sample_.value_second()) )/k ); 195 167 196 if (alpha_new2 > v) 168 197 alpha_new2 = v; 169 198 else if (alpha_new2<u) 170 199 alpha_new2 = u; 171 172 // Updating the alphas and some more 173 if (alpha_new2 < tolerance_) 174 alpha_new2=0; 175 176 alpha_new1 = (alpha_train[index.first] + 177 (target_train[index.first] * target_train[index.second] * 178 (alpha_train[index.second] - alpha_new2)) ); 200 201 202 // Updating the alphas 203 // if alpha is 'zero' make the sample a non-support vector 204 if (alpha_new2 < tolerance_){ 205 sample_.nsv_second(); 206 } 207 else{ 208 sample_.sv_second(); 209 } 210 211 212 alpha_new1 = (alpha_(sample_.value_first()) + 213 (target_(sample_.value_first()) * 214 target_(sample_.value_second()) * 215 (alpha_(sample_.value_second()) - alpha_new2) )); 216 217 // if alpha is 'zero' make the sample a non-support vector 179 218 if (alpha_new1 < tolerance_){ 180 alpha_new1=0; 181 } 182 183 alpha_train[index.first] = alpha_new1; 184 alpha_train[index.second] = alpha_new2; 185 186 E = kernel_train*target_train.mul(alpha_train) - target_train; 187 gslapi::vector ones(alpha_train.size(),1); 188 219 sample_.nsv_first(); 220 } 221 else 222 sample_.sv_first(); 223 224 alpha_(sample_.value_first()) = alpha_new1; 225 alpha_(sample_.value_second()) = alpha_new2; 226 227 // update E vector 228 // Peter, perhaps one should only update SVs, but what happens in choose? 229 for (size_t i=0; i<E.size(); i++) { 230 E(i)+=( kernel_mod(i,sample_.value_first())* 231 target_(sample_.value_first()) * 232 (alpha_new1-alpha_old1) ); 233 E(i)+=( kernel_mod(i,sample_.value_second())* 234 target_(sample_.value_second()) * 235 (alpha_new2-alpha_old2) ); 236 } 237 189 238 epochs++; 190 if (epochs>max_epochs_) {239 if (epochs>max_epochs_) 191 240 return false; 192 } 193 } 194 195 // Calculating the bias from support vectors. E = output - bias - target 196 statistics::Averager bias; 197 for (size_t i=0; i<E.size(); i++) 198 if (alpha_train(i)) 199 bias.add(-E(i)); 200 bias_ = bias.mean(); 201 202 trained_= true; 203 alpha_.set_all(0); 204 for (unsigned int i=0; i<train_set_.size(); i++) 205 alpha_(train_set_[i]) = alpha_train(i); 206 207 // training performed succesfully 208 return true; 209 210 } 211 241 } 242 243 trained_ = calculate_bias(); 244 return trained_; 245 } 246 247 248 bool SVM::choose(const theplu::gslapi::vector& E) 249 { 250 //std::cout << "e choose\n" ; 251 // First check for violation among SVs 252 // E should be the same for all SVs 253 // Choose that pair having largest violation/difference. 254 sample_.update_second(0); 255 sample_.update_first(0); 256 if (sample_.nof_sv()>1){ 257 //std::cout << "there is SVs\n"; 258 double max = E(sample_(0)); 259 double min = max; 260 for (size_t i=1; i<sample_.nof_sv(); i++){ 261 assert(alpha_(sample_(i))>tolerance_); 262 if (E(sample_(i)) > max){ 263 max = E(sample_(i)); 264 sample_.update_second(i); 265 } 266 else if (E(sample_(i))<min){ 267 min = E(sample_(i)); 268 sample_.update_first(i); 269 } 270 } 271 assert(alpha_(sample_.value_first())>tolerance_); 272 assert(alpha_(sample_.value_second())>tolerance_); 273 274 275 if (E(sample_.value_second()) - E(sample_.value_first()) > 2*tolerance_){ 276 return true; 277 } 278 279 // If no violation check among non-support vectors 280 //std::cout << "no violation SVs\n"; 281 282 sample_.shuffle(); 283 //std::cout << "randomized\n"; 284 285 for (size_t i=sample_.nof_sv(); i<sample_.n();i++){ 286 //std::cout << "nr: " << i << std::endl; 287 if (target_(sample_(i))==1){ 288 if(E(sample_(i)) < E(sample_.value_first()) - 2*tolerance_){ 289 sample_.update_second(i); 290 return true; 291 } 292 } 293 else{ 294 if(E(sample_(i)) > E(sample_.value_second()) + 2*tolerance_){ 295 sample_.update_first(i); 296 return true; 297 } 298 } 299 } 300 } 301 302 // if no support vectors - special case 303 else{ 304 for (size_t i=0; i<sample_.n(); i++) { 305 if (target_(sample_(i))==1){ 306 for (size_t j=0; j<sample_.n(); j++) { 307 if ( target_(sample_(j))==-1 && 308 E(sample_(i)) < E(sample_(j))+2*tolerance_ ){ 309 sample_.update_first(i); 310 sample_.update_second(j); 311 return true; 312 } 313 } 314 } 315 } 316 } 317 318 //std::cout << "Done!" << std::endl; 319 // If there is no violation then we should stop training 320 return false; 321 322 } 212 323 213 214 324 325 void SVM::bounds( double& u, double& v) const 326 { 327 if (target_(sample_.value_first())!=target_(sample_.value_second())) { 328 if (alpha_(sample_.value_second()) > alpha_(sample_.value_first())) { 329 v = std::numeric_limits<double>::max(); 330 u = alpha_(sample_.value_second()) - alpha_(sample_.value_first()); 331 } 332 else { 333 v = (std::numeric_limits<double>::max() - 334 alpha_(sample_.value_first()) + 335 alpha_(sample_.value_second())); 336 u = 0; 337 } 338 } 339 else { 340 if (alpha_(sample_.value_second()) + alpha_(sample_.value_first()) > 341 std::numeric_limits<double>::max()) { 342 u = (alpha_(sample_.value_second()) + alpha_(sample_.value_first()) - 343 std::numeric_limits<double>::max()); 344 v = std::numeric_limits<double>::max(); 345 } 346 else { 347 u = 0; 348 v = alpha_(sample_.value_first()) + alpha_(sample_.value_second()); 349 } 350 } 351 } 352 353 bool SVM::calculate_bias(void) 354 { 355 356 // calculating output without bias 357 for (size_t i=0; i<output_.size(); i++) { 358 output_(i)=0; 359 for (size_t j=0; j<output_.size(); j++) 360 output_(i)+=alpha_(j)*target_(j)*kernel_(i,j); 361 } 362 363 if (!sample_.nof_sv()){ 364 std::cerr << "SVM::train() error: " 365 << "Cannot calculate bias because there is no support vector" 366 << std::endl; 367 return false; 368 } 369 370 // For samples with alpha>0, we have: target*output=1-alpha/C 371 bias_=0; 372 for (size_t i=0; i<sample_.nof_sv(); i++) 373 bias_+= ( target_(sample_(i)) * (1-alpha_(sample_(i))*C_inverse_) - 374 output_(sample_(i)) ); 375 bias_=bias_/sample_.nof_sv(); 376 for (size_t i=0; i<output_.size(); i++) 377 output_(i) += bias_; 378 379 return true; 380 } 215 381 216 382 }} // of namespace svm and namespace theplu -
trunk/lib/svm/SVM.h
r318 r323 13 13 namespace theplu { 14 14 namespace svm { 15 /// 15 16 // @internal Class keeping track of which samples are support vectors and 17 // not. The first nof_sv elements in the vector are indices of the 18 // support vectors 19 // 20 class Index 21 { 22 23 public: 24 //Default Contructor 25 Index(); 26 27 // 28 Index(const size_t); 29 30 //not implemented 31 Index(const Index&); 32 33 // @return index_first 34 inline size_t index_first(void) const 35 { assert(index_first_<n()); return index_first_; } 36 37 // @return index_second 38 inline size_t index_second(void) const { assert(index_second_<n());return index_second_; } 39 40 // synch the object against alpha 41 void init(const gslapi::vector& alpha, const double); 42 43 // @return nof samples 44 inline size_t n(void) const { return vec_.size(); } 45 46 // @return nof support vectors 47 inline size_t nof_sv(void) const { return nof_sv_; } 48 49 // making first to an nsv. If already sv, nothing happens. 50 void nsv_first(void); 51 52 // making first to an nsv. If already sv, nothing happens. 53 void nsv_second(void); 54 55 // randomizes the nsv part of vector and sets index_first to 56 // nof_sv_ (the first nsv) 57 void shuffle(void); 58 59 // making first to a sv. If already sv, nothing happens. 60 void sv_first(void); 61 62 // making first to a sv. If already sv, nothing happens. 63 void sv_second(void); 64 65 // 66 void update_first(const size_t); 67 68 // 69 void update_second(const size_t); 70 71 // @return value_first 72 inline size_t value_first(void) const { assert(value_first_<n()); return value_first_; } 73 74 // @return const ref value_second 75 inline size_t value_second(void) const { assert(value_first_<n()); return value_second_; } 76 77 inline size_t operator()(size_t i) const { 78 assert(i<n()); assert(vec_[i]<n()); return vec_[i]; } 79 80 private: 81 size_t index_first_; 82 size_t index_second_; 83 size_t nof_sv_; 84 std::vector<size_t> vec_; 85 size_t value_first_; // vec_[index_first_] exists for fast access 86 size_t value_second_; // vec_[index_second_] exists for fast access 87 88 }; 89 90 /// 16 91 /// Class for SVM using Keerthi's second modification of Platt's SMO. 17 /// @todo speed it up.18 92 /// 19 93 class SVM … … 21 95 22 96 public: 23 /// 24 /// Constructor taking the kernel matrix and the target vector as input 97 /// 98 /// Default constructor (not implemented) 99 /// 100 SVM(); 101 102 /// 103 /// Constructor taking the kernel and the target vector as input. 25 104 /// 26 105 SVM(const Kernel_MEV&, 27 const gslapi::vector&, 28 const std::vector<size_t>& = std::vector<size_t>()); 106 const gslapi::vector&); 107 108 /// 109 /// @todo Copy constructor. 110 /// 111 SVM(const SVM&); 29 112 30 113 /// 31 /// Function returns \f$\alpha\f$ 32 /// 33 inline gslapi::vector get_alpha(void) const { return alpha_; } 34 35 /// 36 /// Function returns the C-parameter 37 /// 38 inline double get_c(void) const { return c_; } 39 40 /// 114 /// The C-parameter is the balance term (see train()). A very 115 /// large C means the training will be focused on getting samples 116 /// correctly classified, with risk for overfitting and poor 117 /// generalisation. A too small C will result in a training where 118 /// misclassifications are not penalized. C is weighted with 119 /// respect to the size, so \f$ n_+C_+ = n_-C_- \f$, meaning a 120 /// misclassificaion of the smaller group is penalized 121 /// harder. This balance is equivalent to the one occuring for 122 /// regression with regularisation, or ANN-training with a 123 /// weight-decay term. Default is C set to infinity. 124 /// 125 /// @returns mean of vector \f$ C_i \f$ 126 /// 127 inline double C(void) const { return 1/C_inverse_; } 128 129 /// 130 /// @return \f$\alpha\f$ 131 /// 132 inline const gslapi::vector& alpha(void) const { return alpha_; } 133 134 /// 135 /// Default is max_epochs set to 10,000,000. 136 /// 41 137 /// @return number of maximal epochs 42 138 /// 43 139 inline long int max_epochs(void) const {return max_epochs_;} 44 140 45 /// 46 /// Changing number of maximal epochs 47 /// 48 inline void max_epochs(const unsigned long int d) {max_epochs_=d;} 49 50 /// 141 /// 142 /// The output is calculated as \f$ o_i = \sum \alpha_j t_j K_{ij} 143 /// + bias \f$, where \f$ t \f$ is the target. 144 /// 51 145 /// @return output 52 /// @todo53 theplu::gslapi::vector output(void) const;54 55 /// 56 /// Changing the C-parameter 57 58 inline void set_c(const double c) {c_ = c;}146 /// 147 inline const theplu::gslapi::vector& output(void) const { return output_; } 148 149 /// 150 /// @todo Function sets \f$ \alpha=0 \f$ and makes SVM untrained. 151 /// 152 void reset(void); 59 153 60 154 /// 61 155 /// Training the SVM following Platt's SMO, with Keerti's 62 /// modifacation. However the complete kernel is stored in 63 /// memory. The reason for this is speed. When number of samples N 64 /// is large this is not possible since the memory cost for the 65 /// kernel scales N^2. In that case one should follow the SMO and 66 /// calculate the kernel elements sequentially. Minimizing \f$ 67 /// \frac{1}{2}\sum 156 /// modifacation. Minimizing \f$ \frac{1}{2}\sum 68 157 /// y_iy_j\alpha_i\alpha_j(K_{ij}+\frac{1}{C_i}\delta_{ij}) \f$, 69 /// which corresponds to minimizing \f$ \sum w_i^2+\sum 70 /// C_i\xi_i^2 \f$ 71 /// 72 73 bool train(void); 74 158 /// which corresponds to minimizing \f$ \sum w_i^2+\sum C_i\xi_i^2 159 /// \f$. train_set defines which samples to use for the 160 /// training. Default is to use all samples. 161 /// 162 bool train(); 163 164 75 165 76 166 private: 77 167 gslapi::vector alpha_; 78 168 double bias_; 79 double c_;80 const Kernel_MEV kernel_;169 double C_inverse_; 170 const Kernel_MEV& kernel_; 81 171 unsigned long int max_epochs_; 82 const gslapi::vector& target_; 172 gslapi::vector output_; 173 Index sample_; 174 const gslapi::vector& target_; 83 175 bool trained_; 84 std::vector<size_t> train_set_;85 176 double tolerance_; 86 177 178 179 /// 180 /// Calculates bounds for alpha2 181 /// 182 void bounds(double&, double&) const; 183 184 /// 185 /// @todo Function calculating bias 186 /// 187 /// @return true if successful 188 /// 189 bool calculate_bias(void); 190 87 191 /// 88 192 /// Private function choosing which two elements that should be 89 193 /// updated. First checking for the biggest violation (output - target = 90 194 /// 0) among support vectors (alpha!=0). If no violation was found check 91 /// for sequentially among the other samples. If no violation there as 92 /// well, stop_condition is fullfilled. 93 /// 94 std::pair<size_t, size_t> choose(const theplu::gslapi::vector&, 95 const theplu::gslapi::vector&, 96 const theplu::gslapi::vector&, 97 bool&); 98 99 195 /// sequentially among the other samples. If no violation there as 196 /// well training is completed 197 /// 198 /// @return true if a pair of samples that violate the conditions 199 /// can be found 200 /// 201 bool choose(const theplu::gslapi::vector&); 202 203 /// 204 /// @return kernel modified in with diagonal term (soft margin) 205 /// 206 inline double kernel_mod(const size_t i, const size_t j) const 207 { return i!=j ? kernel_(i,j) : kernel_(i,j) + C_inverse_; } 100 208 101 209 }; 102 210 103 104 105 211 /// 212 /// @todo The output operator for the SVM class. 213 /// 214 std::ostream& operator<< (std::ostream& s, const SVM&); 215 216 106 217 }} // of namespace svm and namespace theplu 107 218 -
trunk/test/svm_test.cc
r307 r323 14 14 using namespace theplu; 15 15 16 int main() 17 16 int main(const int argc,const char* argv[]) 18 17 { 19 18 20 // ifstream is("data/nm_kernel.txt"); 21 // theplu::gslapi::matrix kernel_matrix(is); 22 // is.close(); 23 24 // is.open("data/nm_target_bin.txt"); 25 // theplu::gslapi::vector target(is); 26 // is.close(); 19 bool print = (argc>1 && argv[1]==std::string("-p")); 20 bool ok = true; 27 21 28 // is.open("data/nm_alpha_linear_matlab.txt"); 29 // theplu::gslapi::vector alpha_matlab(is); 30 // is.close(); 22 gslapi::matrix data2(2,3); 23 data2(0,0)=0; 24 data2(1,0)=0; 25 data2(0,1)=0; 26 data2(1,1)=1; 27 data2(0,2)=1; 28 data2(1,2)=0; 29 gslapi::vector target2(3); 30 target2(0)=-1; 31 target2(1)=1; 32 target2(2)=1; 33 svm::KernelFunction* kf2 = new svm::PolynomialKernelFunction(); 34 svm::Kernel_MEV kernel2(data2,*kf2); 35 assert(kernel2.size()==3); 36 svm::SVM svm2(kernel2, target2); 37 svm2.train(); 31 38 32 // theplu::cpptools::SVM svm(kernel_matrix, target); 33 // svm.train(); 39 if (svm2.alpha()*target2){ 40 std::cerr << "condition not fullfilled" << std::endl; 41 return -1; 42 } 34 43 35 // theplu::gslapi::vector alpha = svm.get_alpha(); 44 if (svm2.alpha()(1)!=2 || svm2.alpha()(2)!=2){ 45 std::cerr << "wrong alpha" << std::endl; 46 std::cerr << "alpha: " << svm2.alpha() << std::endl; 47 std::cerr << "expected: 4 2 2" << std::endl; 48 49 return -1; 50 } 51 52 53 54 std::ifstream is("data/nm_data_centralized.txt"); 55 gslapi::matrix transposed_data(is); 56 is.close(); 57 // Because how the kernel is treated is changed, data must be transposed. 58 gslapi::matrix data=transposed_data; 59 60 svm::KernelFunction* kf = new svm::PolynomialKernelFunction(); 61 svm::Kernel_SEV kernel(data,*kf); 62 63 64 is.open("data/nm_target_bin.txt"); 65 theplu::gslapi::vector target(is); 66 is.close(); 67 68 is.open("data/nm_alpha_linear_matlab.txt"); 69 theplu::gslapi::vector alpha_matlab(is); 70 is.close(); 71 72 theplu::svm::SVM svm(kernel, target); 73 if (!svm.train()){ 74 ok=false; 75 if (print) 76 std::cerr << "Training failured" << std::endl; 77 } 78 79 theplu::gslapi::vector alpha = svm.alpha(); 36 80 37 // // Comparing alpha to alpha_matlab 38 // theplu::gslapi::vector diff_alpha = alpha - alpha_matlab; 39 // if (diff_alpha*diff_alpha> pow(10.0,-10)){ 40 // cerr << "Difference to matlab alphas too large/n"; 41 // return -1; 42 // } 43 // // Comparing output to target 44 // theplu::gslapi::vector output = svm.get_output(); 45 // double slack = 0; 46 // for (unsigned int i=0; i<target.size(); i++){ 47 // if (output[i]*target[i] < 1){ 48 // slack += 1 - output[i]*target[i]; 49 // } 50 // } 51 // if (slack > pow(10.0,-10)){ 52 // cerr << "Difference to matlab alphas too large/n"; 53 // return -1; 54 // } 81 // Comparing alpha to alpha_matlab 82 theplu::gslapi::vector diff_alpha = alpha - alpha_matlab; 83 if (diff_alpha*diff_alpha> 1e-10 ){ 84 if (print) 85 std::cerr << "Difference to matlab alphas too large\n"; 86 ok=false; 87 } 88 89 // Comparing output to target 90 theplu::gslapi::vector output = svm.output(); 91 double slack = 0; 92 for (unsigned int i=0; i<target.size(); i++){ 93 if (output[i]*target[i] < 1){ 94 slack += 1 - output[i]*target[i]; 95 } 96 } 97 if (slack > 1e-7){ 98 if (print){ 99 std::cerr << "Slack too large. Is the bias correct?\n"; 100 std::cerr << "slack: " << slack << std::endl; 101 } 102 ok = false; 103 } 55 104 56 // testing on non solvable problem57 gslapi::matrix data(3,1);58 data(0,0)=-1;59 data(1,0)=10;60 data(2,0)=0;61 62 gslapi::vector target(3,1.0);63 target(2)=-1;64 svm::KernelFunction* kf =65 new svm::PolynomialKernelFunction();66 svm::Kernel_SEV kernel(data,*kf);67 //theplu::gslapi::matrix m=kernel.get();68 //std::cout << "kernel:\n" << m << "\n";69 svm::SVM svm(kernel,target);70 // svm.set_c(1);71 //svm.train();72 105 73 return 0; 106 if(ok) 107 return 0; 108 return -1; 74 109 75 110 }
Note: See TracChangeset
for help on using the changeset viewer.