00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #include <shogun/features/Labels.h>
00013 #include <shogun/mathematics/Math.h>
00014 #include <shogun/lib/DynamicArray.h>
00015 #include <shogun/lib/Time.h>
00016 #include <shogun/base/Parallel.h>
00017 #include <shogun/machine/Machine.h>
00018 #include <shogun/classifier/svm/libocas.h>
00019 #include <shogun/classifier/svm/WDSVMOcas.h>
00020 #include <shogun/features/StringFeatures.h>
00021 #include <shogun/features/Alphabet.h>
00022 #include <shogun/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 : CMachine(), 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 : CMachine(), 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 : CMachine(), 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::apply()
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, apply(i));
00114
00115 return output;
00116 }
00117
00118 return NULL;
00119 }
00120
00121 CLabels* CWDSVMOcas::apply(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 apply();
00134 }
00135
00136 int32_t CWDSVMOcas::set_wd_weights()
00137 {
00138 ASSERT(degree>0 && degree<=8);
00139 SG_FREE(wd_weights);
00140 wd_weights=SG_MALLOC(float32_t, degree);
00141 SG_FREE(w_offsets);
00142 w_offsets=SG_MALLOC(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_machine(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 SGVector<float64_t> labvec=labels->get_labels();
00177 lab=labvec.vector;
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, labvec.vlen);
00188 ASSERT(num_vec==labvec.vlen);
00189 ASSERT(num_vec>0);
00190
00191 SG_FREE(w);
00192 w=SG_MALLOC(float32_t, w_dim);
00193 memset(w, 0, w_dim*sizeof(float32_t));
00194
00195 SG_FREE(old_w);
00196 old_w=SG_MALLOC(float32_t, w_dim);
00197 memset(old_w, 0, w_dim*sizeof(float32_t));
00198 bias=0;
00199 old_bias=0;
00200
00201 cuts=SG_MALLOC(float32_t*, bufsize);
00202 memset(cuts, 0, sizeof(*cuts)*bufsize);
00203 cp_bias=SG_MALLOC(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 SG_FREE(cuts[i]);
00244 SG_FREE(cuts);
00245
00246 lab=NULL;
00247 labvec.free_vector();
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=SG_MALLOC(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 SG_FREE(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=SG_MALLOC(wdocas_thread_params_add, o->parallel->get_num_threads());
00358 pthread_t* threads=SG_MALLOC(pthread_t, o->parallel->get_num_threads());
00359 float32_t* new_a=SG_MALLOC(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 step= string_length/o->parallel->get_num_threads();
00365
00366 if (step<1)
00367 {
00368 nthreads=string_length-1;
00369 step=1;
00370 }
00371
00372 for (t=0; t<nthreads; t++)
00373 {
00374 params_add[t].wdocas=o;
00375
00376 params_add[t].new_a=new_a;
00377 params_add[t].new_cut=new_cut;
00378 params_add[t].start = step*t;
00379 params_add[t].end = step*(t+1);
00380 params_add[t].cut_length = cut_length;
00381
00382 if (pthread_create(&threads[t], NULL, &CWDSVMOcas::add_new_cut_helper, (void*)¶ms_add[t]) != 0)
00383 {
00384 nthreads=t;
00385 SG_SWARNING("thread creation failed\n");
00386 break;
00387 }
00388 }
00389
00390 params_add[t].wdocas=o;
00391
00392 params_add[t].new_a=new_a;
00393 params_add[t].new_cut=new_cut;
00394 params_add[t].start = step*t;
00395 params_add[t].end = string_length;
00396 params_add[t].cut_length = cut_length;
00397 add_new_cut_helper(¶ms_add[t]);
00398
00399
00400 for (t=0; t<nthreads; t++)
00401 {
00402 if (pthread_join(threads[t], NULL) != 0)
00403 SG_SWARNING( "pthread_join failed\n");
00404
00405
00406
00407
00408
00409 }
00410
00411 for(i=0; i < cut_length; i++)
00412 {
00413 if (o->use_bias)
00414 c_bias[nSel]+=o->lab[new_cut[i]];
00415 }
00416
00417
00418 for(i=0; i < nSel; i++)
00419 new_col_H[i] = CMath::dot(new_a, cuts[i], nDim) + c_bias[nSel]*c_bias[i];
00420 new_col_H[nSel] = CMath::dot(new_a, new_a, nDim) + CMath::sq(c_bias[nSel]);
00421
00422 cuts[nSel]=new_a;
00423
00424
00425
00426 SG_FREE(threads);
00427 SG_FREE(params_add);
00428
00429 return 0;
00430 }
00431
00432 int CWDSVMOcas::sort( float64_t* vals, float64_t* data, uint32_t size)
00433 {
00434 CMath::qsort_index(vals, data, size);
00435 return 0;
00436 }
00437
00438
00439
00440
00441
00442
00443 void* CWDSVMOcas::compute_output_helper(void* ptr)
00444 {
00445 wdocas_thread_params_output* p = (wdocas_thread_params_output*) ptr;
00446 CWDSVMOcas* o = p->wdocas;
00447 int32_t start = p->start;
00448 int32_t end = p->end;
00449 float32_t* out = p->out;
00450 float64_t* output = p->output;
00451 int32_t* val = p->val;
00452
00453 CStringFeatures<uint8_t>* f=o->get_features();
00454
00455 int32_t degree = o->degree;
00456 int32_t string_length = o->string_length;
00457 int32_t alphabet_size = o->alphabet_size;
00458 int32_t* w_offsets = o->w_offsets;
00459 float32_t* wd_weights = o->wd_weights;
00460 float32_t* w= o->w;
00461
00462 float64_t* y = o->lab;
00463 float64_t normalization_const = o->normalization_const;
00464
00465
00466 for (int32_t j=0; j<string_length; j++)
00467 {
00468 int32_t offs=o->w_dim_single_char*j;
00469 for (int32_t i=start ; i<end; i++)
00470 val[i]=0;
00471
00472 int32_t lim=CMath::min(degree, string_length-j);
00473 int32_t len;
00474
00475 for (int32_t k=0; k<lim; k++)
00476 {
00477 bool free_vec;
00478 uint8_t* vec=f->get_feature_vector(j+k, len, free_vec);
00479 float32_t wd = wd_weights[k];
00480
00481 for (int32_t i=start; i<end; i++)
00482 {
00483 val[i]=val[i]*alphabet_size + vec[i];
00484 out[i]+=wd*w[offs+val[i]];
00485 }
00486
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 offs+=w_offsets[k];
00523 f->free_feature_vector(vec, j+k, free_vec);
00524 }
00525 }
00526
00527 for (int32_t i=start; i<end; i++)
00528 output[i]=y[i]*o->bias + out[i]*y[i]/normalization_const;
00529
00530
00531
00532 return NULL;
00533 }
00534
00535 int CWDSVMOcas::compute_output( float64_t *output, void* ptr )
00536 {
00537 CWDSVMOcas* o = (CWDSVMOcas*) ptr;
00538 int32_t nData=o->num_vec;
00539 wdocas_thread_params_output* params_output=SG_MALLOC(wdocas_thread_params_output, o->parallel->get_num_threads());
00540 pthread_t* threads = SG_MALLOC(pthread_t, o->parallel->get_num_threads());
00541
00542 float32_t* out=SG_MALLOC(float32_t, nData);
00543 int32_t* val=SG_MALLOC(int32_t, nData);
00544 memset(out, 0, sizeof(float32_t)*nData);
00545
00546 int32_t t;
00547 int32_t nthreads=o->parallel->get_num_threads()-1;
00548 int32_t step= nData/o->parallel->get_num_threads();
00549
00550 if (step<1)
00551 {
00552 nthreads=nData-1;
00553 step=1;
00554 }
00555
00556 for (t=0; t<nthreads; t++)
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 = step*(t+1);
00564
00565
00566 if (pthread_create(&threads[t], NULL, &CWDSVMOcas::compute_output_helper, (void*)¶ms_output[t]) != 0)
00567 {
00568 nthreads=t;
00569 SG_SWARNING("thread creation failed\n");
00570 break;
00571 }
00572 }
00573
00574 params_output[t].wdocas=o;
00575 params_output[t].output=output;
00576 params_output[t].out=out;
00577 params_output[t].val=val;
00578 params_output[t].start = step*t;
00579 params_output[t].end = nData;
00580 compute_output_helper(¶ms_output[t]);
00581
00582
00583 for (t=0; t<nthreads; t++)
00584 {
00585 if (pthread_join(threads[t], NULL) != 0)
00586 SG_SWARNING( "pthread_join failed\n");
00587 }
00588 SG_FREE(threads);
00589 SG_FREE(params_output);
00590 SG_FREE(val);
00591 SG_FREE(out);
00592 return 0;
00593 }
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603 void CWDSVMOcas::compute_W(
00604 float64_t *sq_norm_W, float64_t *dp_WoldW, float64_t *alpha, uint32_t nSel,
00605 void* ptr)
00606 {
00607 CWDSVMOcas* o = (CWDSVMOcas*) ptr;
00608 uint32_t nDim= (uint32_t) o->w_dim;
00609 CMath::swap(o->w, o->old_w);
00610 float32_t* W=o->w;
00611 float32_t* oldW=o->old_w;
00612 float32_t** cuts=o->cuts;
00613 memset(W, 0, sizeof(float32_t)*nDim);
00614 float64_t* c_bias = o->cp_bias;
00615 float64_t old_bias=o->bias;
00616 float64_t bias=0;
00617
00618 for (uint32_t i=0; i<nSel; i++)
00619 {
00620 if (alpha[i] > 0)
00621 CMath::vec1_plus_scalar_times_vec2(W, (float32_t) alpha[i], cuts[i], nDim);
00622
00623 bias += c_bias[i]*alpha[i];
00624 }
00625
00626 *sq_norm_W = CMath::dot(W,W, nDim) +CMath::sq(bias);
00627 *dp_WoldW = CMath::dot(W,oldW, nDim) + bias*old_bias;;
00628
00629
00630 o->bias = bias;
00631 o->old_bias = old_bias;
00632 }