SHOGUN  5.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules
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 
13 #ifdef USE_GPL_SHOGUN
14 
15 #include <shogun/labels/Labels.h>
17 #include <shogun/lib/Time.h>
18 #include <shogun/base/Parameter.h>
19 #include <shogun/base/Parallel.h>
22 #include <shogun/labels/Labels.h>
24 
25 using namespace shogun;
26 
27 CSVMOcas::CSVMOcas()
29 {
30  init();
31 }
32 
33 CSVMOcas::CSVMOcas(E_SVM_TYPE type)
35 {
36  init();
37  method=type;
38 }
39 
40 CSVMOcas::CSVMOcas(
41  float64_t C, CDotFeatures* traindat, CLabels* trainlab)
43 {
44  init();
45  C1=C;
46  C2=C;
47 
48  set_features(traindat);
49  set_labels(trainlab);
50 }
51 
52 
53 CSVMOcas::~CSVMOcas()
54 {
55 }
56 
57 bool CSVMOcas::train_machine(CFeatures* data)
58 {
59  SG_INFO("C=%f, epsilon=%f, bufsize=%d\n", get_C1(), get_epsilon(), bufsize)
60  SG_DEBUG("use_bias = %i\n", get_bias_enabled())
61 
62  ASSERT(m_labels)
63  ASSERT(m_labels->get_label_type() == LT_BINARY)
64  if (data)
65  {
66  if (!data->has_property(FP_DOT))
67  SG_ERROR("Specified features are not of type CDotFeatures\n")
68  set_features((CDotFeatures*) data);
69  }
70  ASSERT(features)
71 
72  int32_t num_vec=features->get_num_vectors();
73  lab = SGVector<float64_t>(num_vec);
74  for (int32_t i=0; i<num_vec; i++)
75  lab[i] = ((CBinaryLabels*)m_labels)->get_label(i);
76 
77  w=SGVector<float64_t>(features->get_dim_feature_space());
78  w.zero();
79 
80  if (num_vec!=lab.vlen || num_vec<=0)
81  SG_ERROR("num_vec=%d num_train_labels=%d\n", num_vec, lab.vlen)
82 
83  SG_FREE(old_w);
84  old_w=SG_CALLOC(float64_t, w.vlen);
85  bias=0;
86  old_bias=0;
87 
88  tmp_a_buf=SG_CALLOC(float64_t, w.vlen);
89  cp_value=SG_CALLOC(float64_t*, bufsize);
90  cp_index=SG_CALLOC(uint32_t*, bufsize);
91  cp_nz_dims=SG_CALLOC(uint32_t, bufsize);
92  cp_bias=SG_CALLOC(float64_t, bufsize);
93 
94  float64_t TolAbs=0;
95  float64_t QPBound=0;
96  int32_t Method=0;
97  if (method == SVM_OCAS)
98  Method = 1;
99  ocas_return_value_T result = svm_ocas_solver( get_C1(), num_vec, get_epsilon(),
100  TolAbs, QPBound, get_max_train_time(), bufsize, Method,
101  &CSVMOcas::compute_W,
102  &CSVMOcas::update_W,
103  &CSVMOcas::add_new_cut,
104  &CSVMOcas::compute_output,
105  &CSVMOcas::sort,
106  &CSVMOcas::print,
107  this);
108 
109  SG_INFO("Ocas Converged after %d iterations\n"
110  "==================================\n"
111  "timing statistics:\n"
112  "output_time: %f s\n"
113  "sort_time: %f s\n"
114  "add_time: %f s\n"
115  "w_time: %f s\n"
116  "solver_time %f s\n"
117  "ocas_time %f s\n\n", result.nIter, result.output_time, result.sort_time,
118  result.add_time, result.w_time, result.qp_solver_time, result.ocas_time);
119 
120  SG_FREE(tmp_a_buf);
121 
122  primal_objective = result.Q_P;
123 
124  uint32_t num_cut_planes = result.nCutPlanes;
125 
126  SG_DEBUG("num_cut_planes=%d\n", num_cut_planes)
127  for (uint32_t i=0; i<num_cut_planes; i++)
128  {
129  SG_DEBUG("cp_value[%d]=%p\n", i, cp_value)
130  SG_FREE(cp_value[i]);
131  SG_DEBUG("cp_index[%d]=%p\n", i, cp_index)
132  SG_FREE(cp_index[i]);
133  }
134 
135  SG_FREE(cp_value);
136  cp_value=NULL;
137  SG_FREE(cp_index);
138  cp_index=NULL;
139  SG_FREE(cp_nz_dims);
140  cp_nz_dims=NULL;
141  SG_FREE(cp_bias);
142  cp_bias=NULL;
143 
144  SG_FREE(old_w);
145  old_w=NULL;
146 
147  return true;
148 }
149 
150 /*----------------------------------------------------------------------------------
151  sq_norm_W = sparse_update_W( t ) does the following:
152 
153  W = oldW*(1-t) + t*W;
154  sq_norm_W = W'*W;
155 
156  ---------------------------------------------------------------------------------*/
157 float64_t CSVMOcas::update_W( float64_t t, void* ptr )
158 {
159  float64_t sq_norm_W = 0;
160  CSVMOcas* o = (CSVMOcas*) ptr;
161  uint32_t nDim = (uint32_t) o->w.vlen;
162  float64_t* W=o->w.vector;
163  float64_t* oldW=o->old_w;
164 
165  for(uint32_t j=0; j <nDim; j++)
166  {
167  W[j] = oldW[j]*(1-t) + t*W[j];
168  sq_norm_W += W[j]*W[j];
169  }
170  o->bias=o->old_bias*(1-t) + t*o->bias;
171  sq_norm_W += CMath::sq(o->bias);
172 
173  return( sq_norm_W );
174 }
175 
176 /*----------------------------------------------------------------------------------
177  sparse_add_new_cut( new_col_H, new_cut, cut_length, nSel ) does the following:
178 
179  new_a = sum(data_X(:,find(new_cut ~=0 )),2);
180  new_col_H = [sparse_A(:,1:nSel)'*new_a ; new_a'*new_a];
181  sparse_A(:,nSel+1) = new_a;
182 
183  ---------------------------------------------------------------------------------*/
184 int CSVMOcas::add_new_cut(
185  float64_t *new_col_H, uint32_t *new_cut, uint32_t cut_length,
186  uint32_t nSel, void* ptr)
187 {
188  CSVMOcas* o = (CSVMOcas*) ptr;
189  CDotFeatures* f = o->features;
190  uint32_t nDim=(uint32_t) o->w.vlen;
191  float64_t* y = o->lab.vector;
192 
193  float64_t** c_val = o->cp_value;
194  uint32_t** c_idx = o->cp_index;
195  uint32_t* c_nzd = o->cp_nz_dims;
196  float64_t* c_bias = o->cp_bias;
197 
198  float64_t sq_norm_a;
199  uint32_t i, j, nz_dims;
200 
201  /* temporary vector */
202  float64_t* new_a = o->tmp_a_buf;
203  memset(new_a, 0, sizeof(float64_t)*nDim);
204 
205  for(i=0; i < cut_length; i++)
206  {
207  f->add_to_dense_vec(y[new_cut[i]], new_cut[i], new_a, nDim);
208 
209  if (o->use_bias)
210  c_bias[nSel]+=y[new_cut[i]];
211  }
212 
213  /* compute new_a'*new_a and count number of non-zerou dimensions */
214  nz_dims = 0;
215  sq_norm_a = CMath::sq(c_bias[nSel]);
216  for(j=0; j < nDim; j++ ) {
217  if(new_a[j] != 0) {
218  nz_dims++;
219  sq_norm_a += new_a[j]*new_a[j];
220  }
221  }
222 
223  /* sparsify new_a and insert it to the last column of sparse_A */
224  c_nzd[nSel] = nz_dims;
225  c_idx[nSel]=NULL;
226  c_val[nSel]=NULL;
227 
228  if(nz_dims > 0)
229  {
230  c_idx[nSel]=SG_MALLOC(uint32_t, nz_dims);
231  c_val[nSel]=SG_MALLOC(float64_t, nz_dims);
232 
233  uint32_t idx=0;
234  for(j=0; j < nDim; j++ )
235  {
236  if(new_a[j] != 0)
237  {
238  c_idx[nSel][idx] = j;
239  c_val[nSel][idx++] = new_a[j];
240  }
241  }
242  }
243 
244  new_col_H[nSel] = sq_norm_a;
245 
246  for(i=0; i < nSel; i++)
247  {
248  float64_t tmp = c_bias[nSel]*c_bias[i];
249  for(j=0; j < c_nzd[i]; j++)
250  tmp += new_a[c_idx[i][j]]*c_val[i][j];
251 
252  new_col_H[i] = tmp;
253  }
254  //CMath::display_vector(new_col_H, nSel+1, "new_col_H");
255  //CMath::display_vector((int32_t*) c_idx[nSel], (int32_t) nz_dims, "c_idx");
256  //CMath::display_vector((float64_t*) c_val[nSel], nz_dims, "c_val");
257  return 0;
258 }
259 
260 int CSVMOcas::sort(float64_t* vals, float64_t* data, uint32_t size)
261 {
262  CMath::qsort_index(vals, data, size);
263  return 0;
264 }
265 
266 /*----------------------------------------------------------------------
267  sparse_compute_output( output ) does the follwing:
268 
269  output = data_X'*W;
270  ----------------------------------------------------------------------*/
271 int CSVMOcas::compute_output(float64_t *output, void* ptr)
272 {
273  CSVMOcas* o = (CSVMOcas*) ptr;
274  CDotFeatures* f=o->features;
275  int32_t nData=f->get_num_vectors();
276 
277  float64_t* y = o->lab.vector;
278 
279  f->dense_dot_range(output, 0, nData, y, o->w.vector, o->w.vlen, 0.0);
280 
281  for (int32_t i=0; i<nData; i++)
282  output[i]+=y[i]*o->bias;
283  //CMath::display_vector(o->w, o->w.vlen, "w");
284  //CMath::display_vector(output, nData, "out");
285  return 0;
286 }
287 
288 /*----------------------------------------------------------------------
289  sq_norm_W = compute_W( alpha, nSel ) does the following:
290 
291  oldW = W;
292  W = sparse_A(:,1:nSel)'*alpha;
293  sq_norm_W = W'*W;
294  dp_WoldW = W'*oldW';
295 
296  ----------------------------------------------------------------------*/
297 void CSVMOcas::compute_W(
298  float64_t *sq_norm_W, float64_t *dp_WoldW, float64_t *alpha, uint32_t nSel,
299  void* ptr )
300 {
301  CSVMOcas* o = (CSVMOcas*) ptr;
302  uint32_t nDim= (uint32_t) o->w.vlen;
303  CMath::swap(o->w.vector, o->old_w);
304  float64_t* W=o->w.vector;
305  float64_t* oldW=o->old_w;
306  memset(W, 0, sizeof(float64_t)*nDim);
307  float64_t old_bias=o->bias;
308  float64_t bias=0;
309 
310  float64_t** c_val = o->cp_value;
311  uint32_t** c_idx = o->cp_index;
312  uint32_t* c_nzd = o->cp_nz_dims;
313  float64_t* c_bias = o->cp_bias;
314 
315  for(uint32_t i=0; i<nSel; i++)
316  {
317  uint32_t nz_dims = c_nzd[i];
318 
319  if(nz_dims > 0 && alpha[i] > 0)
320  {
321  for(uint32_t j=0; j < nz_dims; j++)
322  W[c_idx[i][j]] += alpha[i]*c_val[i][j];
323  }
324  bias += c_bias[i]*alpha[i];
325  }
326 
327  *sq_norm_W = CMath::dot(W,W, nDim) + CMath::sq(bias);
328  *dp_WoldW = CMath::dot(W,oldW, nDim) + bias*old_bias;
329  //SG_PRINT("nSel=%d sq_norm_W=%f dp_WoldW=%f\n", nSel, *sq_norm_W, *dp_WoldW)
330 
331  o->bias = bias;
332  o->old_bias = old_bias;
333 }
334 
335 void CSVMOcas::init()
336 {
337  use_bias=true;
338  bufsize=3000;
339  C1=1;
340  C2=1;
341 
342  epsilon=1e-3;
343  method=SVM_OCAS;
344  old_w=NULL;
345  tmp_a_buf=NULL;
346  cp_value=NULL;
347  cp_index=NULL;
348  cp_nz_dims=NULL;
349  cp_bias=NULL;
350 
351  primal_objective = 0.0;
352 
353  m_parameters->add(&C1, "C1", "Cost constant 1.");
354  m_parameters->add(&C2, "C2", "Cost constant 2.");
355  m_parameters->add(&use_bias, "use_bias",
356  "Indicates if bias is used.");
357  m_parameters->add(&epsilon, "epsilon", "Convergence precision.");
358  m_parameters->add(&bufsize, "bufsize", "Maximum number of cutting planes.");
359  m_parameters->add((machine_int_t*) &method, "method",
360  "SVMOcas solver type.");
361 }
362 
363 float64_t CSVMOcas::compute_primal_objective() const
364 {
365  return primal_objective;
366 }
367 
368 #endif //USE_GPL_SHOGUN
virtual void dense_dot_range(float64_t *output, int32_t start, int32_t stop, float64_t *alphas, float64_t *vec, int32_t dim, float64_t b)
Definition: DotFeatures.cpp:67
#define SG_INFO(...)
Definition: SGIO.h:118
binary labels +1/-1
Definition: LabelTypes.h:18
The class Labels models labels, i.e. class assignments of objects.
Definition: Labels.h:43
float64_t C1
Definition: SVMLin.h:109
static void qsort_index(T1 *output, T2 *index, uint32_t size)
Definition: Math.h:2202
static T sq(T x)
Definition: Math.h:450
virtual int32_t get_num_vectors() const =0
#define SG_ERROR(...)
Definition: SGIO.h:129
virtual void add_to_dense_vec(float64_t alpha, int32_t vec_idx1, float64_t *vec2, int32_t vec2_len, bool abs_val=false)=0
Features that support dot products among other operations.
Definition: DotFeatures.h:44
#define ASSERT(x)
Definition: SGIO.h:201
void print(CJLCoverTreePoint &p)
shogun vector
double float64_t
Definition: common.h:50
virtual void set_features(CDotFeatures *feat)
Class LinearMachine is a generic interface for all kinds of linear machines like classifiers.
Definition: LinearMachine.h:63
static float64_t dot(const bool *v1, const bool *v2, int32_t n)
Compute dot product between v1 and v2 (blas optimized)
Definition: Math.h:627
#define SG_DEBUG(...)
Definition: SGIO.h:107
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
int machine_int_t
Definition: common.h:59
The class Features is the base class of all feature objects.
Definition: Features.h:68
Binary Labels for binary classification.
Definition: BinaryLabels.h:37
static void swap(T &a, T &b)
Definition: Math.h:438
virtual void set_labels(CLabels *lab)
Definition: Machine.cpp:65
float64_t C2
Definition: SVMLin.h:111

SHOGUN Machine Learning Toolbox - Documentation