00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #include "features/Labels.h"
00013 #include "lib/Mathematics.h"
00014 #include "lib/DynamicArray.h"
00015 #include "lib/Time.h"
00016 #include "base/Parallel.h"
00017 #include "classifier/Classifier.h"
00018 #include "classifier/svm/libocas.h"
00019 #include "classifier/svm/WDSVMOcas.h"
00020 #include "features/StringFeatures.h"
00021 #include "features/Alphabet.h"
00022 #include "features/Labels.h"
00023
00024 using namespace shogun;
00025
00026 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00027 struct wdocas_thread_params_output
00028 {
00029 float32_t* out;
00030 int32_t* val;
00031 float64_t* output;
00032 CWDSVMOcas* wdocas;
00033 int32_t start;
00034 int32_t end;
00035 };
00036
00037 struct wdocas_thread_params_add
00038 {
00039 CWDSVMOcas* wdocas;
00040 float32_t* new_a;
00041 uint32_t* new_cut;
00042 int32_t start;
00043 int32_t end;
00044 uint32_t cut_length;
00045 };
00046 #endif // DOXYGEN_SHOULD_SKIP_THIS
00047
00048 CWDSVMOcas::CWDSVMOcas(E_SVM_TYPE type)
00049 : CClassifier(), use_bias(false), bufsize(3000), C1(1), C2(1),
00050 epsilon(1e-3), method(type)
00051 {
00052 w=NULL;
00053 old_w=NULL;
00054 features=NULL;
00055 degree=6;
00056 from_degree=40;
00057 wd_weights=NULL;
00058 w_offsets=NULL;
00059 normalization_const=1.0;
00060 }
00061
00062 CWDSVMOcas::CWDSVMOcas(
00063 float64_t C, int32_t d, int32_t from_d, CStringFeatures<uint8_t>* traindat,
00064 CLabels* trainlab)
00065 : CClassifier(), use_bias(false), bufsize(3000), C1(C), C2(C), epsilon(1e-3),
00066 degree(d), from_degree(from_d)
00067 {
00068 w=NULL;
00069 old_w=NULL;
00070 method=SVM_OCAS;
00071 features=traindat;
00072 set_labels(trainlab);
00073 wd_weights=NULL;
00074 w_offsets=NULL;
00075 normalization_const=1.0;
00076 }
00077
00078
00079 CWDSVMOcas::~CWDSVMOcas()
00080 {
00081 }
00082
00083 CLabels* CWDSVMOcas::classify()
00084 {
00085 set_wd_weights();
00086 set_normalization_const();
00087
00088 if (features)
00089 {
00090 int32_t num=features->get_num_vectors();
00091 ASSERT(num>0);
00092
00093 CLabels* output=new CLabels(num);
00094 SG_REF(output);
00095
00096 for (int32_t i=0; i<num; i++)
00097 output->set_label(i, classify_example(i));
00098
00099 return output;
00100 }
00101
00102 return NULL;
00103 }
00104
00105 CLabels* CWDSVMOcas::classify(CFeatures* data)
00106 {
00107 if (!data)
00108 SG_ERROR("No features specified\n");
00109
00110 if (data->get_feature_class() != C_STRING ||
00111 data->get_feature_type() != F_BYTE)
00112 {
00113 SG_ERROR("Features not of class string type byte\n");
00114 }
00115
00116 set_features((CStringFeatures<uint8_t>*) data);
00117 return classify();
00118 }
00119
00120 int32_t CWDSVMOcas::set_wd_weights()
00121 {
00122 ASSERT(degree>0 && degree<=8);
00123 delete[] wd_weights;
00124 wd_weights=new float32_t[degree];
00125 delete[] w_offsets;
00126 w_offsets=new int32_t[degree];
00127 int32_t w_dim_single_c=0;
00128
00129 for (int32_t i=0; i<degree; i++)
00130 {
00131 w_offsets[i]=CMath::pow(alphabet_size, i+1);
00132 wd_weights[i]=sqrt(2.0*(from_degree-i)/(from_degree*(from_degree+1)));
00133 w_dim_single_c+=w_offsets[i];
00134 }
00135 return w_dim_single_c;
00136 }
00137
00138 bool CWDSVMOcas::train(CFeatures* data)
00139 {
00140 SG_INFO("C=%f, epsilon=%f, bufsize=%d\n", get_C1(), get_epsilon(), bufsize);
00141
00142 ASSERT(labels);
00143 if (data)
00144 {
00145 if (data->get_feature_class() != C_STRING ||
00146 data->get_feature_type() != F_BYTE)
00147 {
00148 SG_ERROR("Features not of class string type byte\n");
00149 }
00150 set_features((CStringFeatures<uint8_t>*) data);
00151 }
00152
00153 ASSERT(get_features());
00154 ASSERT(labels->is_two_class_labeling());
00155 CAlphabet* alphabet=get_features()->get_alphabet();
00156 ASSERT(alphabet && alphabet->get_alphabet()==RAWDNA);
00157
00158 alphabet_size=alphabet->get_num_symbols();
00159 string_length=features->get_num_vectors();
00160 int32_t num_train_labels=0;
00161 lab=labels->get_labels(num_train_labels);
00162
00163 w_dim_single_char=set_wd_weights();
00164
00165 SG_DEBUG("w_dim_single_char=%d\n", w_dim_single_char);
00166 w_dim=string_length*w_dim_single_char;
00167 SG_DEBUG("cutting plane has %d dims\n", w_dim);
00168 num_vec=get_features()->get_max_vector_length();
00169
00170 set_normalization_const();
00171 SG_INFO("num_vec: %d num_lab: %d\n", num_vec, num_train_labels);
00172 ASSERT(num_vec==num_train_labels);
00173 ASSERT(num_vec>0);
00174
00175 delete[] w;
00176 w=new float32_t[w_dim];
00177 memset(w, 0, w_dim*sizeof(float32_t));
00178
00179 delete[] old_w;
00180 old_w=new float32_t[w_dim];
00181 memset(old_w, 0, w_dim*sizeof(float32_t));
00182 bias=0;
00183 old_bias=0;
00184
00185 cuts=new float32_t*[bufsize];
00186 memset(cuts, 0, sizeof(*cuts)*bufsize);
00187 cp_bias=new float64_t[bufsize];
00188 memset(cp_bias, 0, sizeof(float64_t)*bufsize);
00189
00191
00192
00193
00194
00195
00196
00197
00198
00200 float64_t TolAbs=0;
00201 float64_t QPBound=0;
00202 uint8_t Method=0;
00203 if (method == SVM_OCAS)
00204 Method = 1;
00205 ocas_return_value_T result = svm_ocas_solver( get_C1(), num_vec, get_epsilon(),
00206 TolAbs, QPBound, bufsize, Method,
00207 &CWDSVMOcas::compute_W,
00208 &CWDSVMOcas::update_W,
00209 &CWDSVMOcas::add_new_cut,
00210 &CWDSVMOcas::compute_output,
00211 &CWDSVMOcas::sort,
00212 this);
00213
00214 SG_INFO("Ocas Converged after %d iterations\n"
00215 "==================================\n"
00216 "timing statistics:\n"
00217 "output_time: %f s\n"
00218 "sort_time: %f s\n"
00219 "add_time: %f s\n"
00220 "w_time: %f s\n"
00221 "solver_time %f s\n"
00222 "ocas_time %f s\n\n", result.nIter, result.output_time, result.sort_time,
00223 result.add_time, result.w_time, result.solver_time, result.ocas_time);
00224
00225 for (int32_t i=bufsize-1; i>=0; i--)
00226 delete[] cuts[i];
00227 delete[] cuts;
00228
00229 delete[] lab;
00230 lab=NULL;
00231
00232 SG_UNREF(alphabet);
00233
00234 return true;
00235 }
00236
00237
00238
00239
00240
00241
00242
00243
00244 float64_t CWDSVMOcas::update_W( float64_t t, void* ptr )
00245 {
00246 float64_t sq_norm_W = 0;
00247 CWDSVMOcas* o = (CWDSVMOcas*) ptr;
00248 uint32_t nDim = (uint32_t) o->w_dim;
00249 float32_t* W=o->w;
00250 float32_t* oldW=o->old_w;
00251 float64_t bias=o->bias;
00252 float64_t old_bias=bias;
00253
00254 for(uint32_t j=0; j <nDim; j++)
00255 {
00256 W[j] = oldW[j]*(1-t) + t*W[j];
00257 sq_norm_W += W[j]*W[j];
00258 }
00259
00260 bias=old_bias*(1-t) + t*bias;
00261 sq_norm_W += CMath::sq(bias);
00262
00263 o->bias=bias;
00264 o->old_bias=old_bias;
00265
00266 return( sq_norm_W );
00267 }
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277 void* CWDSVMOcas::add_new_cut_helper( void* ptr)
00278 {
00279 wdocas_thread_params_add* p = (wdocas_thread_params_add*) ptr;
00280 CWDSVMOcas* o = p->wdocas;
00281 int32_t start = p->start;
00282 int32_t end = p->end;
00283 int32_t string_length = o->string_length;
00284
00285 uint32_t cut_length=p->cut_length;
00286 uint32_t* new_cut=p->new_cut;
00287 int32_t* w_offsets = o->w_offsets;
00288 float64_t* y = o->lab;
00289 int32_t alphabet_size = o->alphabet_size;
00290 float32_t* wd_weights = o->wd_weights;
00291 int32_t degree = o->degree;
00292 CStringFeatures<uint8_t>* f = o->features;
00293 float64_t normalization_const = o->normalization_const;
00294
00295
00296 float32_t* new_a = p->new_a;
00297
00298
00299
00300 int32_t* val=new int32_t[cut_length];
00301 for (int32_t j=start; j<end; j++)
00302 {
00303 int32_t offs=o->w_dim_single_char*j;
00304 memset(val,0,sizeof(int32_t)*cut_length);
00305 int32_t lim=CMath::min(degree, string_length-j);
00306 int32_t len;
00307
00308 for (int32_t k=0; k<lim; k++)
00309 {
00310 bool free_vec;
00311 uint8_t* vec = f->get_feature_vector(j+k, len, free_vec);
00312 float32_t wd = wd_weights[k]/normalization_const;
00313
00314 for(uint32_t i=0; i < cut_length; i++)
00315 {
00316 val[i]=val[i]*alphabet_size + vec[new_cut[i]];
00317 new_a[offs+val[i]]+=wd * y[new_cut[i]];
00318 }
00319 offs+=w_offsets[k];
00320 f->free_feature_vector(vec, j+k, free_vec);
00321 }
00322 }
00323
00324
00325 delete[] val;
00326 return NULL;
00327 }
00328
00329 void CWDSVMOcas::add_new_cut(
00330 float64_t *new_col_H, uint32_t *new_cut, uint32_t cut_length,
00331 uint32_t nSel, void* ptr)
00332 {
00333 CWDSVMOcas* o = (CWDSVMOcas*) ptr;
00334 int32_t string_length = o->string_length;
00335 uint32_t nDim=(uint32_t) o->w_dim;
00336 float32_t** cuts=o->cuts;
00337 float64_t* c_bias = o->cp_bias;
00338
00339 uint32_t i;
00340 wdocas_thread_params_add* params_add=new wdocas_thread_params_add[o->parallel->get_num_threads()];
00341 pthread_t* threads=new pthread_t[o->parallel->get_num_threads()];
00342 float32_t* new_a=new float32_t[nDim];
00343 memset(new_a, 0, sizeof(float32_t)*nDim);
00344
00345 int32_t t;
00346 int32_t nthreads=o->parallel->get_num_threads()-1;
00347 int32_t end=0;
00348 int32_t step= string_length/o->parallel->get_num_threads();
00349
00350 if (step<1)
00351 {
00352 nthreads=string_length-1;
00353 step=1;
00354 }
00355
00356 for (t=0; t<nthreads; t++)
00357 {
00358 params_add[t].wdocas=o;
00359
00360 params_add[t].new_a=new_a;
00361 params_add[t].new_cut=new_cut;
00362 params_add[t].start = step*t;
00363 params_add[t].end = step*(t+1);
00364 params_add[t].cut_length = cut_length;
00365
00366 if (pthread_create(&threads[t], NULL, &CWDSVMOcas::add_new_cut_helper, (void*)¶ms_add[t]) != 0)
00367 {
00368 nthreads=t;
00369 SG_SWARNING("thread creation failed\n");
00370 break;
00371 }
00372 end=params_add[t].end;
00373 }
00374
00375 params_add[t].wdocas=o;
00376
00377 params_add[t].new_a=new_a;
00378 params_add[t].new_cut=new_cut;
00379 params_add[t].start = step*t;
00380 params_add[t].end = string_length;
00381 params_add[t].cut_length = cut_length;
00382 add_new_cut_helper(¶ms_add[t]);
00383
00384
00385 for (t=0; t<nthreads; t++)
00386 {
00387 if (pthread_join(threads[t], NULL) != 0)
00388 SG_SWARNING( "pthread_join failed\n");
00389
00390
00391
00392
00393
00394 }
00395
00396 for(i=0; i < cut_length; i++)
00397 {
00398 if (o->use_bias)
00399 c_bias[nSel]+=o->lab[new_cut[i]];
00400 }
00401
00402
00403 for(i=0; i < nSel; i++)
00404 new_col_H[i] = CMath::dot(new_a, cuts[i], nDim) + c_bias[nSel]*c_bias[i];
00405 new_col_H[nSel] = CMath::dot(new_a, new_a, nDim) + CMath::sq(c_bias[nSel]);
00406
00407 cuts[nSel]=new_a;
00408
00409
00410
00411 delete[] threads;
00412 delete[] params_add;
00413 }
00414
00415 void CWDSVMOcas::sort( float64_t* vals, uint32_t* idx, uint32_t size)
00416 {
00417 CMath::qsort_index(vals, idx, size);
00418 }
00419
00420
00421
00422
00423
00424
00425 void* CWDSVMOcas::compute_output_helper(void* ptr)
00426 {
00427 wdocas_thread_params_output* p = (wdocas_thread_params_output*) ptr;
00428 CWDSVMOcas* o = p->wdocas;
00429 int32_t start = p->start;
00430 int32_t end = p->end;
00431 float32_t* out = p->out;
00432 float64_t* output = p->output;
00433 int32_t* val = p->val;
00434
00435 CStringFeatures<uint8_t>* f=o->get_features();
00436
00437 int32_t degree = o->degree;
00438 int32_t string_length = o->string_length;
00439 int32_t alphabet_size = o->alphabet_size;
00440 int32_t* w_offsets = o->w_offsets;
00441 float32_t* wd_weights = o->wd_weights;
00442 float32_t* w= o->w;
00443
00444 float64_t* y = o->lab;
00445 float64_t normalization_const = o->normalization_const;
00446
00447
00448 for (int32_t j=0; j<string_length; j++)
00449 {
00450 int32_t offs=o->w_dim_single_char*j;
00451 for (int32_t i=start ; i<end; i++)
00452 val[i]=0;
00453
00454 int32_t lim=CMath::min(degree, string_length-j);
00455 int32_t len;
00456
00457 for (int32_t k=0; k<lim; k++)
00458 {
00459 bool free_vec;
00460 uint8_t* vec=f->get_feature_vector(j+k, len, free_vec);
00461 float32_t wd = wd_weights[k];
00462
00463 for (int32_t i=start; i<end; i++)
00464 {
00465 val[i]=val[i]*alphabet_size + vec[i];
00466 out[i]+=wd*w[offs+val[i]];
00467 }
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504 offs+=w_offsets[k];
00505 f->free_feature_vector(vec, j+k, free_vec);
00506 }
00507 }
00508
00509 for (int32_t i=start; i<end; i++)
00510 output[i]=y[i]*o->bias + out[i]*y[i]/normalization_const;
00511
00512
00513
00514 return NULL;
00515 }
00516
00517 void CWDSVMOcas::compute_output( float64_t *output, void* ptr )
00518 {
00519 CWDSVMOcas* o = (CWDSVMOcas*) ptr;
00520 int32_t nData=o->num_vec;
00521 wdocas_thread_params_output* params_output=new wdocas_thread_params_output[o->parallel->get_num_threads()];
00522 pthread_t* threads = new pthread_t[o->parallel->get_num_threads()];
00523
00524 float32_t* out=new float32_t[nData];
00525 int32_t* val=new int32_t[nData];
00526 memset(out, 0, sizeof(float32_t)*nData);
00527
00528 int32_t t;
00529 int32_t nthreads=o->parallel->get_num_threads()-1;
00530 int32_t end=0;
00531 int32_t step= nData/o->parallel->get_num_threads();
00532
00533 if (step<1)
00534 {
00535 nthreads=nData-1;
00536 step=1;
00537 }
00538
00539 for (t=0; t<nthreads; t++)
00540 {
00541 params_output[t].wdocas=o;
00542 params_output[t].output=output;
00543 params_output[t].out=out;
00544 params_output[t].val=val;
00545 params_output[t].start = step*t;
00546 params_output[t].end = step*(t+1);
00547
00548
00549 if (pthread_create(&threads[t], NULL, &CWDSVMOcas::compute_output_helper, (void*)¶ms_output[t]) != 0)
00550 {
00551 nthreads=t;
00552 SG_SWARNING("thread creation failed\n");
00553 break;
00554 }
00555 end=params_output[t].end;
00556 }
00557
00558 params_output[t].wdocas=o;
00559 params_output[t].output=output;
00560 params_output[t].out=out;
00561 params_output[t].val=val;
00562 params_output[t].start = step*t;
00563 params_output[t].end = nData;
00564 compute_output_helper(¶ms_output[t]);
00565
00566
00567 for (t=0; t<nthreads; t++)
00568 {
00569 if (pthread_join(threads[t], NULL) != 0)
00570 SG_SWARNING( "pthread_join failed\n");
00571 }
00572 delete[] threads;
00573 delete[] params_output;
00574 delete[] val;
00575 delete[] out;
00576 }
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586 void CWDSVMOcas::compute_W(
00587 float64_t *sq_norm_W, float64_t *dp_WoldW, float64_t *alpha, uint32_t nSel,
00588 void* ptr)
00589 {
00590 CWDSVMOcas* o = (CWDSVMOcas*) ptr;
00591 uint32_t nDim= (uint32_t) o->w_dim;
00592 CMath::swap(o->w, o->old_w);
00593 float32_t* W=o->w;
00594 float32_t* oldW=o->old_w;
00595 float32_t** cuts=o->cuts;
00596 memset(W, 0, sizeof(float32_t)*nDim);
00597 float64_t* c_bias = o->cp_bias;
00598 float64_t old_bias=o->bias;
00599 float64_t bias=0;
00600
00601 for (uint32_t i=0; i<nSel; i++)
00602 {
00603 if (alpha[i] > 0)
00604 CMath::vec1_plus_scalar_times_vec2(W, (float32_t) alpha[i], cuts[i], nDim);
00605
00606 bias += c_bias[i]*alpha[i];
00607 }
00608
00609 *sq_norm_W = CMath::dot(W,W, nDim) +CMath::sq(bias);
00610 *dp_WoldW = CMath::dot(W,oldW, nDim) + bias*old_bias;;
00611
00612
00613 o->bias = bias;
00614 o->old_bias = old_bias;
00615 }