SHOGUN  v1.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SVMOcas.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>
14 #include <shogun/lib/Time.h>
15 #include <shogun/base/Parameter.h>
16 #include <shogun/base/Parallel.h>
20 #include <shogun/features/Labels.h>
21 
22 using namespace shogun;
23 
26 {
27  init();
28 }
29 
30 CSVMOcas::CSVMOcas(E_SVM_TYPE type)
32 {
33  init();
34  method=type;
35 }
36 
38  float64_t C, CDotFeatures* traindat, CLabels* trainlab)
40 {
41  init();
42  C1=C;
43  C2=C;
44 
45  set_features(traindat);
46  set_labels(trainlab);
47 }
48 
49 
51 {
52 }
53 
55 {
56  SG_INFO("C=%f, epsilon=%f, bufsize=%d\n", get_C1(), get_epsilon(), bufsize);
57  SG_DEBUG("use_bias = %i\n", get_bias_enabled()) ;
58 
59  ASSERT(labels);
60  if (data)
61  {
62  if (!data->has_property(FP_DOT))
63  SG_ERROR("Specified features are not of type CDotFeatures\n");
64  set_features((CDotFeatures*) data);
65  }
68 
71  int32_t num_vec=features->get_num_vectors();
72 
73  if (num_vec!=lab.vlen || num_vec<=0)
74  SG_ERROR("num_vec=%d num_train_labels=%d\n", num_vec, lab.vlen);
75 
76  SG_FREE(w);
78  memset(w, 0, w_dim*sizeof(float64_t));
79 
80  SG_FREE(old_w);
82  memset(old_w, 0, w_dim*sizeof(float64_t));
83  bias=0;
84  old_bias=0;
85 
88  memset(cp_value, sizeof(float64_t*)*bufsize, 0);
89  cp_index=SG_MALLOC(uint32_t*, bufsize);
90  memset(cp_index, sizeof(float64_t*)*bufsize, 0);
91  cp_nz_dims=SG_MALLOC(uint32_t, bufsize);
93  memset(cp_bias, 0, sizeof(float64_t)*bufsize);
94 
95  float64_t TolAbs=0;
96  float64_t QPBound=0;
97  int32_t Method=0;
98  if (method == SVM_OCAS)
99  Method = 1;
100  ocas_return_value_T result = svm_ocas_solver( get_C1(), num_vec, get_epsilon(),
101  TolAbs, QPBound, get_max_train_time(), bufsize, Method,
108  this);
109 
110  SG_INFO("Ocas Converged after %d iterations\n"
111  "==================================\n"
112  "timing statistics:\n"
113  "output_time: %f s\n"
114  "sort_time: %f s\n"
115  "add_time: %f s\n"
116  "w_time: %f s\n"
117  "solver_time %f s\n"
118  "ocas_time %f s\n\n", result.nIter, result.output_time, result.sort_time,
119  result.add_time, result.w_time, result.qp_solver_time, result.ocas_time);
120 
122 
123  uint32_t num_cut_planes = result.nCutPlanes;
124 
125  SG_DEBUG("num_cut_planes=%d\n", num_cut_planes);
126  for (uint32_t i=0; i<num_cut_planes; i++)
127  {
128  SG_DEBUG("cp_value[%d]=%p\n", i, cp_value);
129  SG_FREE(cp_value[i]);
130  SG_DEBUG("cp_index[%d]=%p\n", i, cp_index);
131  SG_FREE(cp_index[i]);
132  }
133 
134  SG_FREE(cp_value);
135  cp_value=NULL;
136  SG_FREE(cp_index);
137  cp_index=NULL;
139  cp_nz_dims=NULL;
140  SG_FREE(cp_bias);
141  cp_bias=NULL;
142 
143  lab.free_vector();
144 
145  SG_FREE(old_w);
146  old_w=NULL;
147 
148  return true;
149 }
150 
151 /*----------------------------------------------------------------------------------
152  sq_norm_W = sparse_update_W( t ) does the following:
153 
154  W = oldW*(1-t) + t*W;
155  sq_norm_W = W'*W;
156 
157  ---------------------------------------------------------------------------------*/
159 {
160  float64_t sq_norm_W = 0;
161  CSVMOcas* o = (CSVMOcas*) ptr;
162  uint32_t nDim = (uint32_t) o->w_dim;
163  float64_t* W=o->w;
164  float64_t* oldW=o->old_w;
165 
166  for(uint32_t j=0; j <nDim; j++)
167  {
168  W[j] = oldW[j]*(1-t) + t*W[j];
169  sq_norm_W += W[j]*W[j];
170  }
171  o->bias=o->old_bias*(1-t) + t*o->bias;
172  sq_norm_W += CMath::sq(o->bias);
173 
174  return( sq_norm_W );
175 }
176 
177 /*----------------------------------------------------------------------------------
178  sparse_add_new_cut( new_col_H, new_cut, cut_length, nSel ) does the following:
179 
180  new_a = sum(data_X(:,find(new_cut ~=0 )),2);
181  new_col_H = [sparse_A(:,1:nSel)'*new_a ; new_a'*new_a];
182  sparse_A(:,nSel+1) = new_a;
183 
184  ---------------------------------------------------------------------------------*/
186  float64_t *new_col_H, uint32_t *new_cut, uint32_t cut_length,
187  uint32_t nSel, void* ptr)
188 {
189  CSVMOcas* o = (CSVMOcas*) ptr;
190  CDotFeatures* f = o->features;
191  uint32_t nDim=(uint32_t) o->w_dim;
192  float64_t* y = o->lab.vector;
193 
194  float64_t** c_val = o->cp_value;
195  uint32_t** c_idx = o->cp_index;
196  uint32_t* c_nzd = o->cp_nz_dims;
197  float64_t* c_bias = o->cp_bias;
198 
199  float64_t sq_norm_a;
200  uint32_t i, j, nz_dims;
201 
202  /* temporary vector */
203  float64_t* new_a = o->tmp_a_buf;
204  memset(new_a, 0, sizeof(float64_t)*nDim);
205 
206  for(i=0; i < cut_length; i++)
207  {
208  f->add_to_dense_vec(y[new_cut[i]], new_cut[i], new_a, nDim);
209 
210  if (o->use_bias)
211  c_bias[nSel]+=y[new_cut[i]];
212  }
213 
214  /* compute new_a'*new_a and count number of non-zerou dimensions */
215  nz_dims = 0;
216  sq_norm_a = CMath::sq(c_bias[nSel]);
217  for(j=0; j < nDim; j++ ) {
218  if(new_a[j] != 0) {
219  nz_dims++;
220  sq_norm_a += new_a[j]*new_a[j];
221  }
222  }
223 
224  /* sparsify new_a and insert it to the last column of sparse_A */
225  c_nzd[nSel] = nz_dims;
226  c_idx[nSel]=NULL;
227  c_val[nSel]=NULL;
228 
229  if(nz_dims > 0)
230  {
231  c_idx[nSel]=SG_MALLOC(uint32_t, nz_dims);
232  c_val[nSel]=SG_MALLOC(float64_t, nz_dims);
233 
234  uint32_t idx=0;
235  for(j=0; j < nDim; j++ )
236  {
237  if(new_a[j] != 0)
238  {
239  c_idx[nSel][idx] = j;
240  c_val[nSel][idx++] = new_a[j];
241  }
242  }
243  }
244 
245  new_col_H[nSel] = sq_norm_a;
246 
247  for(i=0; i < nSel; i++)
248  {
249  float64_t tmp = c_bias[nSel]*c_bias[i];
250  for(j=0; j < c_nzd[i]; j++)
251  tmp += new_a[c_idx[i][j]]*c_val[i][j];
252 
253  new_col_H[i] = tmp;
254  }
255  //CMath::display_vector(new_col_H, nSel+1, "new_col_H");
256  //CMath::display_vector((int32_t*) c_idx[nSel], (int32_t) nz_dims, "c_idx");
257  //CMath::display_vector((float64_t*) c_val[nSel], nz_dims, "c_val");
258  return 0;
259 }
260 
261 int CSVMOcas::sort(float64_t* vals, float64_t* data, uint32_t size)
262 {
263  CMath::qsort_index(vals, data, size);
264  return 0;
265 }
266 
267 /*----------------------------------------------------------------------
268  sparse_compute_output( output ) does the follwing:
269 
270  output = data_X'*W;
271  ----------------------------------------------------------------------*/
272 int CSVMOcas::compute_output(float64_t *output, void* ptr)
273 {
274  CSVMOcas* o = (CSVMOcas*) ptr;
275  CDotFeatures* f=o->features;
276  int32_t nData=f->get_num_vectors();
277 
278  float64_t* y = o->lab.vector;
279 
280  f->dense_dot_range(output, 0, nData, y, o->w, o->w_dim, 0.0);
281 
282  for (int32_t i=0; i<nData; i++)
283  output[i]+=y[i]*o->bias;
284  //CMath::display_vector(o->w, o->w_dim, "w");
285  //CMath::display_vector(output, nData, "out");
286  return 0;
287 }
288 
289 /*----------------------------------------------------------------------
290  sq_norm_W = compute_W( alpha, nSel ) does the following:
291 
292  oldW = W;
293  W = sparse_A(:,1:nSel)'*alpha;
294  sq_norm_W = W'*W;
295  dp_WoldW = W'*oldW';
296 
297  ----------------------------------------------------------------------*/
299  float64_t *sq_norm_W, float64_t *dp_WoldW, float64_t *alpha, uint32_t nSel,
300  void* ptr )
301 {
302  CSVMOcas* o = (CSVMOcas*) ptr;
303  uint32_t nDim= (uint32_t) o->w_dim;
304  CMath::swap(o->w, o->old_w);
305  float64_t* W=o->w;
306  float64_t* oldW=o->old_w;
307  memset(W, 0, sizeof(float64_t)*nDim);
309  float64_t bias=0;
310 
311  float64_t** c_val = o->cp_value;
312  uint32_t** c_idx = o->cp_index;
313  uint32_t* c_nzd = o->cp_nz_dims;
314  float64_t* c_bias = o->cp_bias;
315 
316  for(uint32_t i=0; i<nSel; i++)
317  {
318  uint32_t nz_dims = c_nzd[i];
319 
320  if(nz_dims > 0 && alpha[i] > 0)
321  {
322  for(uint32_t j=0; j < nz_dims; j++)
323  W[c_idx[i][j]] += alpha[i]*c_val[i][j];
324  }
325  bias += c_bias[i]*alpha[i];
326  }
327 
328  *sq_norm_W = CMath::dot(W,W, nDim) + CMath::sq(bias);
329  *dp_WoldW = CMath::dot(W,oldW, nDim) + bias*old_bias;
330  //SG_PRINT("nSel=%d sq_norm_W=%f dp_WoldW=%f\n", nSel, *sq_norm_W, *dp_WoldW);
331 
332  o->bias = bias;
333  o->old_bias = old_bias;
334 }
335 
336 void CSVMOcas::init()
337 {
338  use_bias=true;
339  bufsize=3000;
340  C1=1;
341  C2=1;
342 
343  epsilon=1e-3;
344  method=SVM_OCAS;
345  w=NULL;
346  old_w=NULL;
347  tmp_a_buf=NULL;
349  cp_value=NULL;
350  cp_index=NULL;
351  cp_nz_dims=NULL;
352  cp_bias=NULL;
353 
354  m_parameters->add(&C1, "C1", "Cost constant 1.");
355  m_parameters->add(&C2, "C2", "Cost constant 2.");
356  m_parameters->add(&use_bias, "use_bias",
357  "Indicates if bias is used.");
358  m_parameters->add(&epsilon, "epsilon", "Convergence precision.");
359  m_parameters->add(&bufsize, "bufsize", "Maximum number of cutting planes.");
360  m_parameters->add((machine_int_t*) &method, "method",
361  "SVMOcas solver type.");
362 }

SHOGUN Machine Learning Toolbox - Documentation