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