SHOGUN  v1.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
WDSVMOcas.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2007-2008 Vojtech Franc
8  * Written (W) 2007-2009 Soeren Sonnenburg
9  * Copyright (C) 2007-2009 Fraunhofer Institute FIRST and Max-Planck-Society
10  */
11 
12 #include <shogun/features/Labels.h>
15 #include <shogun/lib/Time.h>
16 #include <shogun/base/Parallel.h>
17 #include <shogun/machine/Machine.h>
22 #include <shogun/features/Labels.h>
23 
24 using namespace shogun;
25 
26 #ifndef DOXYGEN_SHOULD_SKIP_THIS
27 struct wdocas_thread_params_output
28 {
29  float32_t* out;
30  int32_t* val;
31  float64_t* output;
32  CWDSVMOcas* wdocas;
33  int32_t start;
34  int32_t end;
35 };
36 
37 struct wdocas_thread_params_add
38 {
39  CWDSVMOcas* wdocas;
40  float32_t* new_a;
41  uint32_t* new_cut;
42  int32_t start;
43  int32_t end;
44  uint32_t cut_length;
45 };
46 #endif // DOXYGEN_SHOULD_SKIP_THIS
47 
49 : CMachine(), use_bias(false), bufsize(3000), C1(1), C2(1),
50  epsilon(1e-3), method(SVM_OCAS)
51 {
52  SG_UNSTABLE("CWDSVMOcas::CWDSVMOcas()", "\n");
53 
54  w=NULL;
55  old_w=NULL;
56  features=NULL;
57  degree=6;
58  from_degree=40;
59  wd_weights=NULL;
60  w_offsets=NULL;
62 }
63 
64 CWDSVMOcas::CWDSVMOcas(E_SVM_TYPE type)
65 : CMachine(), use_bias(false), bufsize(3000), C1(1), C2(1),
66  epsilon(1e-3), method(type)
67 {
68  w=NULL;
69  old_w=NULL;
70  features=NULL;
71  degree=6;
72  from_degree=40;
73  wd_weights=NULL;
74  w_offsets=NULL;
76 }
77 
79  float64_t C, int32_t d, int32_t from_d, CStringFeatures<uint8_t>* traindat,
80  CLabels* trainlab)
81 : CMachine(), use_bias(false), bufsize(3000), C1(C), C2(C), epsilon(1e-3),
82  degree(d), from_degree(from_d)
83 {
84  w=NULL;
85  old_w=NULL;
86  method=SVM_OCAS;
87  features=traindat;
88  set_labels(trainlab);
89  wd_weights=NULL;
90  w_offsets=NULL;
92 }
93 
94 
96 {
97 }
98 
100 {
101  set_wd_weights();
103 
104  if (features)
105  {
106  int32_t num=features->get_num_vectors();
107  ASSERT(num>0);
108 
109  CLabels* output=new CLabels(num);
110  SG_REF(output);
111 
112  for (int32_t i=0; i<num; i++)
113  output->set_label(i, apply(i));
114 
115  return output;
116  }
117 
118  return NULL;
119 }
120 
122 {
123  if (!data)
124  SG_ERROR("No features specified\n");
125 
126  if (data->get_feature_class() != C_STRING ||
127  data->get_feature_type() != F_BYTE)
128  {
129  SG_ERROR("Features not of class string type byte\n");
130  }
131 
133  return apply();
134 }
135 
137 {
138  ASSERT(degree>0 && degree<=8);
142  w_offsets=SG_MALLOC(int32_t, degree);
143  int32_t w_dim_single_c=0;
144 
145  for (int32_t i=0; i<degree; i++)
146  {
148  wd_weights[i]=sqrt(2.0*(from_degree-i)/(from_degree*(from_degree+1)));
149  w_dim_single_c+=w_offsets[i];
150  }
151  return w_dim_single_c;
152 }
153 
155 {
156  SG_INFO("C=%f, epsilon=%f, bufsize=%d\n", get_C1(), get_epsilon(), bufsize);
157 
158  ASSERT(labels);
159  if (data)
160  {
161  if (data->get_feature_class() != C_STRING ||
162  data->get_feature_type() != F_BYTE)
163  {
164  SG_ERROR("Features not of class string type byte\n");
165  }
167  }
168 
169  ASSERT(get_features());
171  CAlphabet* alphabet=get_features()->get_alphabet();
172  ASSERT(alphabet && alphabet->get_alphabet()==RAWDNA);
173 
174  alphabet_size=alphabet->get_num_symbols();
177  lab=labvec.vector;
178 
180  //CMath::display_vector(wd_weights, degree, "wd_weights");
181  SG_DEBUG("w_dim_single_char=%d\n", w_dim_single_char);
183  SG_DEBUG("cutting plane has %d dims\n", w_dim);
185 
187  SG_INFO("num_vec: %d num_lab: %d\n", num_vec, labvec.vlen);
188  ASSERT(num_vec==labvec.vlen);
189  ASSERT(num_vec>0);
190 
191  SG_FREE(w);
193  memset(w, 0, w_dim*sizeof(float32_t));
194 
195  SG_FREE(old_w);
197  memset(old_w, 0, w_dim*sizeof(float32_t));
198  bias=0;
199  old_bias=0;
200 
202  memset(cuts, 0, sizeof(*cuts)*bufsize);
203  cp_bias=SG_MALLOC(float64_t, bufsize);
204  memset(cp_bias, 0, sizeof(float64_t)*bufsize);
205 
207  /*float64_t* tmp = SG_MALLOC(float64_t, num_vec);
208  float64_t start=CTime::get_curtime();
209  CMath::random_vector(w, w_dim, (float32_t) 0, (float32_t) 1000);
210  compute_output(tmp, this);
211  start=CTime::get_curtime()-start;
212  SG_PRINT("timing:%f\n", start);
213  SG_FREE(tmp);
214  exit(1);*/
216  float64_t TolAbs=0;
217  float64_t QPBound=0;
218  uint8_t Method=0;
219  if (method == SVM_OCAS)
220  Method = 1;
221  ocas_return_value_T result = svm_ocas_solver( get_C1(), num_vec, get_epsilon(),
222  TolAbs, QPBound, get_max_train_time(), bufsize, Method,
229  this);
230 
231  SG_INFO("Ocas Converged after %d iterations\n"
232  "==================================\n"
233  "timing statistics:\n"
234  "output_time: %f s\n"
235  "sort_time: %f s\n"
236  "add_time: %f s\n"
237  "w_time: %f s\n"
238  "solver_time %f s\n"
239  "ocas_time %f s\n\n", result.nIter, result.output_time, result.sort_time,
240  result.add_time, result.w_time, result.qp_solver_time, result.ocas_time);
241 
242  for (int32_t i=bufsize-1; i>=0; i--)
243  SG_FREE(cuts[i]);
244  SG_FREE(cuts);
245 
246  lab=NULL;
247  labvec.free_vector();
248 
249  SG_UNREF(alphabet);
250 
251  return true;
252 }
253 
254 /*----------------------------------------------------------------------------------
255  sq_norm_W = sparse_update_W( t ) does the following:
256 
257  W = oldW*(1-t) + t*W;
258  sq_norm_W = W'*W;
259 
260  ---------------------------------------------------------------------------------*/
262 {
263  float64_t sq_norm_W = 0;
264  CWDSVMOcas* o = (CWDSVMOcas*) ptr;
265  uint32_t nDim = (uint32_t) o->w_dim;
266  float32_t* W=o->w;
267  float32_t* oldW=o->old_w;
268  float64_t bias=o->bias;
270 
271  for(uint32_t j=0; j <nDim; j++)
272  {
273  W[j] = oldW[j]*(1-t) + t*W[j];
274  sq_norm_W += W[j]*W[j];
275  }
276 
277  bias=old_bias*(1-t) + t*bias;
278  sq_norm_W += CMath::sq(bias);
279 
280  o->bias=bias;
281  o->old_bias=old_bias;
282 
283  return( sq_norm_W );
284 }
285 
286 /*----------------------------------------------------------------------------------
287  sparse_add_new_cut( new_col_H, new_cut, cut_length, nSel ) does the following:
288 
289  new_a = sum(data_X(:,find(new_cut ~=0 )),2);
290  new_col_H = [sparse_A(:,1:nSel)'*new_a ; new_a'*new_a];
291  sparse_A(:,nSel+1) = new_a;
292 
293  ---------------------------------------------------------------------------------*/
295 {
296  wdocas_thread_params_add* p = (wdocas_thread_params_add*) ptr;
297  CWDSVMOcas* o = p->wdocas;
298  int32_t start = p->start;
299  int32_t end = p->end;
300  int32_t string_length = o->string_length;
301  //uint32_t nDim=(uint32_t) o->w_dim;
302  uint32_t cut_length=p->cut_length;
303  uint32_t* new_cut=p->new_cut;
304  int32_t* w_offsets = o->w_offsets;
305  float64_t* y = o->lab;
306  int32_t alphabet_size = o->alphabet_size;
308  int32_t degree = o->degree;
311 
312  // temporary vector
313  float32_t* new_a = p->new_a;
314  //float32_t* new_a = SG_MALLOC(float32_t, nDim);
315  //memset(new_a, 0, sizeof(float32_t)*nDim);
316 
317  int32_t* val=SG_MALLOC(int32_t, cut_length);
318  for (int32_t j=start; j<end; j++)
319  {
320  int32_t offs=o->w_dim_single_char*j;
321  memset(val,0,sizeof(int32_t)*cut_length);
322  int32_t lim=CMath::min(degree, string_length-j);
323  int32_t len;
324 
325  for (int32_t k=0; k<lim; k++)
326  {
327  bool free_vec;
328  uint8_t* vec = f->get_feature_vector(j+k, len, free_vec);
329  float32_t wd = wd_weights[k]/normalization_const;
330 
331  for(uint32_t i=0; i < cut_length; i++)
332  {
333  val[i]=val[i]*alphabet_size + vec[new_cut[i]];
334  new_a[offs+val[i]]+=wd * y[new_cut[i]];
335  }
336  offs+=w_offsets[k];
337  f->free_feature_vector(vec, j+k, free_vec);
338  }
339  }
340 
341  //p->new_a=new_a;
342  SG_FREE(val);
343  return NULL;
344 }
345 
347  float64_t *new_col_H, uint32_t *new_cut, uint32_t cut_length,
348  uint32_t nSel, void* ptr)
349 {
350  CWDSVMOcas* o = (CWDSVMOcas*) ptr;
351  int32_t string_length = o->string_length;
352  uint32_t nDim=(uint32_t) o->w_dim;
353  float32_t** cuts=o->cuts;
354  float64_t* c_bias = o->cp_bias;
355 
356  uint32_t i;
357  wdocas_thread_params_add* params_add=SG_MALLOC(wdocas_thread_params_add, o->parallel->get_num_threads());
358  pthread_t* threads=SG_MALLOC(pthread_t, o->parallel->get_num_threads());
359  float32_t* new_a=SG_MALLOC(float32_t, nDim);
360  memset(new_a, 0, sizeof(float32_t)*nDim);
361 
362  int32_t t;
363  int32_t nthreads=o->parallel->get_num_threads()-1;
364  int32_t step= string_length/o->parallel->get_num_threads();
365 
366  if (step<1)
367  {
368  nthreads=string_length-1;
369  step=1;
370  }
371 
372  for (t=0; t<nthreads; t++)
373  {
374  params_add[t].wdocas=o;
375  //params_add[t].new_a=NULL;
376  params_add[t].new_a=new_a;
377  params_add[t].new_cut=new_cut;
378  params_add[t].start = step*t;
379  params_add[t].end = step*(t+1);
380  params_add[t].cut_length = cut_length;
381 
382  if (pthread_create(&threads[t], NULL, &CWDSVMOcas::add_new_cut_helper, (void*)&params_add[t]) != 0)
383  {
384  nthreads=t;
385  SG_SWARNING("thread creation failed\n");
386  break;
387  }
388  }
389 
390  params_add[t].wdocas=o;
391  //params_add[t].new_a=NULL;
392  params_add[t].new_a=new_a;
393  params_add[t].new_cut=new_cut;
394  params_add[t].start = step*t;
395  params_add[t].end = string_length;
396  params_add[t].cut_length = cut_length;
397  add_new_cut_helper(&params_add[t]);
398  //float32_t* new_a=params_add[t].new_a;
399 
400  for (t=0; t<nthreads; t++)
401  {
402  if (pthread_join(threads[t], NULL) != 0)
403  SG_SWARNING( "pthread_join failed\n");
404 
405  //float32_t* a=params_add[t].new_a;
406  //for (i=0; i<nDim; i++)
407  // new_a[i]+=a[i];
408  //SG_FREE(a);
409  }
410 
411  for(i=0; i < cut_length; i++)
412  {
413  if (o->use_bias)
414  c_bias[nSel]+=o->lab[new_cut[i]];
415  }
416 
417  // insert new_a into the last column of sparse_A
418  for(i=0; i < nSel; i++)
419  new_col_H[i] = CMath::dot(new_a, cuts[i], nDim) + c_bias[nSel]*c_bias[i];
420  new_col_H[nSel] = CMath::dot(new_a, new_a, nDim) + CMath::sq(c_bias[nSel]);
421 
422  cuts[nSel]=new_a;
423  //CMath::display_vector(new_col_H, nSel+1, "new_col_H");
424  //CMath::display_vector(cuts[nSel], nDim, "cut[nSel]");
425  //
426  SG_FREE(threads);
427  SG_FREE(params_add);
428 
429  return 0;
430 }
431 
432 int CWDSVMOcas::sort( float64_t* vals, float64_t* data, uint32_t size)
433 {
434  CMath::qsort_index(vals, data, size);
435  return 0;
436 }
437 
438 /*----------------------------------------------------------------------
439  sparse_compute_output( output ) does the follwing:
440 
441  output = data_X'*W;
442  ----------------------------------------------------------------------*/
444 {
445  wdocas_thread_params_output* p = (wdocas_thread_params_output*) ptr;
446  CWDSVMOcas* o = p->wdocas;
447  int32_t start = p->start;
448  int32_t end = p->end;
449  float32_t* out = p->out;
450  float64_t* output = p->output;
451  int32_t* val = p->val;
452 
454 
455  int32_t degree = o->degree;
456  int32_t string_length = o->string_length;
457  int32_t alphabet_size = o->alphabet_size;
458  int32_t* w_offsets = o->w_offsets;
460  float32_t* w= o->w;
461 
462  float64_t* y = o->lab;
464 
465 
466  for (int32_t j=0; j<string_length; j++)
467  {
468  int32_t offs=o->w_dim_single_char*j;
469  for (int32_t i=start ; i<end; i++)
470  val[i]=0;
471 
472  int32_t lim=CMath::min(degree, string_length-j);
473  int32_t len;
474 
475  for (int32_t k=0; k<lim; k++)
476  {
477  bool free_vec;
478  uint8_t* vec=f->get_feature_vector(j+k, len, free_vec);
479  float32_t wd = wd_weights[k];
480 
481  for (int32_t i=start; i<end; i++) // quite fast 1.9s
482  {
483  val[i]=val[i]*alphabet_size + vec[i];
484  out[i]+=wd*w[offs+val[i]];
485  }
486 
487  /*for (int32_t i=0; i<nData/4; i++) // slowest 2s
488  {
489  uint32_t x=((uint32_t*) vec)[i];
490  int32_t ii=4*i;
491  val[ii]=val[ii]*alphabet_size + (x&255);
492  val[ii+1]=val[ii+1]*alphabet_size + ((x>>8)&255);
493  val[ii+2]=val[ii+2]*alphabet_size + ((x>>16)&255);
494  val[ii+3]=val[ii+3]*alphabet_size + (x>>24);
495  out[ii]+=wd*w[offs+val[ii]];
496  out[ii+1]+=wd*w[offs+val[ii+1]];
497  out[ii+2]+=wd*w[offs+val[ii+2]];
498  out[ii+3]+=wd*w[offs+val[ii+3]];
499  }*/
500 
501  /*for (int32_t i=0; i<nData>>3; i++) // fastest on 64bit: 1.5s
502  {
503  uint64_t x=((uint64_t*) vec)[i];
504  int32_t ii=i<<3;
505  val[ii]=val[ii]*alphabet_size + (x&255);
506  val[ii+1]=val[ii+1]*alphabet_size + ((x>>8)&255);
507  val[ii+2]=val[ii+2]*alphabet_size + ((x>>16)&255);
508  val[ii+3]=val[ii+3]*alphabet_size + ((x>>24)&255);
509  val[ii+4]=val[ii+4]*alphabet_size + ((x>>32)&255);
510  val[ii+5]=val[ii+5]*alphabet_size + ((x>>40)&255);
511  val[ii+6]=val[ii+6]*alphabet_size + ((x>>48)&255);
512  val[ii+7]=val[ii+7]*alphabet_size + (x>>56);
513  out[ii]+=wd*w[offs+val[ii]];
514  out[ii+1]+=wd*w[offs+val[ii+1]];
515  out[ii+2]+=wd*w[offs+val[ii+2]];
516  out[ii+3]+=wd*w[offs+val[ii+3]];
517  out[ii+4]+=wd*w[offs+val[ii+4]];
518  out[ii+5]+=wd*w[offs+val[ii+5]];
519  out[ii+6]+=wd*w[offs+val[ii+6]];
520  out[ii+7]+=wd*w[offs+val[ii+7]];
521  }*/
522  offs+=w_offsets[k];
523  f->free_feature_vector(vec, j+k, free_vec);
524  }
525  }
526 
527  for (int32_t i=start; i<end; i++)
528  output[i]=y[i]*o->bias + out[i]*y[i]/normalization_const;
529 
530  //CMath::display_vector(o->w, o->w_dim, "w");
531  //CMath::display_vector(output, nData, "out");
532  return NULL;
533 }
534 
535 int CWDSVMOcas::compute_output( float64_t *output, void* ptr )
536 {
537  CWDSVMOcas* o = (CWDSVMOcas*) ptr;
538  int32_t nData=o->num_vec;
539  wdocas_thread_params_output* params_output=SG_MALLOC(wdocas_thread_params_output, o->parallel->get_num_threads());
540  pthread_t* threads = SG_MALLOC(pthread_t, o->parallel->get_num_threads());
541 
542  float32_t* out=SG_MALLOC(float32_t, nData);
543  int32_t* val=SG_MALLOC(int32_t, nData);
544  memset(out, 0, sizeof(float32_t)*nData);
545 
546  int32_t t;
547  int32_t nthreads=o->parallel->get_num_threads()-1;
548  int32_t step= nData/o->parallel->get_num_threads();
549 
550  if (step<1)
551  {
552  nthreads=nData-1;
553  step=1;
554  }
555 
556  for (t=0; t<nthreads; t++)
557  {
558  params_output[t].wdocas=o;
559  params_output[t].output=output;
560  params_output[t].out=out;
561  params_output[t].val=val;
562  params_output[t].start = step*t;
563  params_output[t].end = step*(t+1);
564 
565  //SG_SPRINT("t=%d start=%d end=%d output=%p\n", t, params_output[t].start, params_output[t].end, params_output[t].output);
566  if (pthread_create(&threads[t], NULL, &CWDSVMOcas::compute_output_helper, (void*)&params_output[t]) != 0)
567  {
568  nthreads=t;
569  SG_SWARNING("thread creation failed\n");
570  break;
571  }
572  }
573 
574  params_output[t].wdocas=o;
575  params_output[t].output=output;
576  params_output[t].out=out;
577  params_output[t].val=val;
578  params_output[t].start = step*t;
579  params_output[t].end = nData;
580  compute_output_helper(&params_output[t]);
581  //SG_SPRINT("t=%d start=%d end=%d output=%p\n", t, params_output[t].start, params_output[t].end, params_output[t].output);
582 
583  for (t=0; t<nthreads; t++)
584  {
585  if (pthread_join(threads[t], NULL) != 0)
586  SG_SWARNING( "pthread_join failed\n");
587  }
588  SG_FREE(threads);
589  SG_FREE(params_output);
590  SG_FREE(val);
591  SG_FREE(out);
592  return 0;
593 }
594 /*----------------------------------------------------------------------
595  sq_norm_W = compute_W( alpha, nSel ) does the following:
596 
597  oldW = W;
598  W = sparse_A(:,1:nSel)'*alpha;
599  sq_norm_W = W'*W;
600  dp_WoldW = W'*oldW';
601 
602  ----------------------------------------------------------------------*/
604  float64_t *sq_norm_W, float64_t *dp_WoldW, float64_t *alpha, uint32_t nSel,
605  void* ptr)
606 {
607  CWDSVMOcas* o = (CWDSVMOcas*) ptr;
608  uint32_t nDim= (uint32_t) o->w_dim;
609  CMath::swap(o->w, o->old_w);
610  float32_t* W=o->w;
611  float32_t* oldW=o->old_w;
612  float32_t** cuts=o->cuts;
613  memset(W, 0, sizeof(float32_t)*nDim);
614  float64_t* c_bias = o->cp_bias;
616  float64_t bias=0;
617 
618  for (uint32_t i=0; i<nSel; i++)
619  {
620  if (alpha[i] > 0)
621  CMath::vec1_plus_scalar_times_vec2(W, (float32_t) alpha[i], cuts[i], nDim);
622 
623  bias += c_bias[i]*alpha[i];
624  }
625 
626  *sq_norm_W = CMath::dot(W,W, nDim) +CMath::sq(bias);
627  *dp_WoldW = CMath::dot(W,oldW, nDim) + bias*old_bias;;
628  //SG_PRINT("nSel=%d sq_norm_W=%f dp_WoldW=%f\n", nSel, *sq_norm_W, *dp_WoldW);
629 
630  o->bias = bias;
631  o->old_bias = old_bias;
632 }

SHOGUN Machine Learning Toolbox - Documentation