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