SVMOcas.cpp

Go to the documentation of this file.
00001 /*
00002  * This program is free software; you can redistribute it and/or modify
00003  * it under the terms of the GNU General Public License as published by
00004  * the Free Software Foundation; either version 3 of the License, or
00005  * (at your option) any later version.
00006  *
00007  * Written (W) 2007-2008 Vojtech Franc
00008  * Written (W) 2007-2009 Soeren Sonnenburg
00009  * Copyright (C) 2007-2009 Fraunhofer Institute FIRST and Max-Planck-Society
00010  */
00011 
00012 #include "features/Labels.h"
00013 #include "lib/Mathematics.h"
00014 #include "lib/Time.h"
00015 #include "base/Parallel.h"
00016 #include "classifier/LinearClassifier.h"
00017 #include "classifier/svm/SVMOcas.h"
00018 #include "classifier/svm/libocas.h"
00019 #include "features/DotFeatures.h"
00020 #include "features/Labels.h"
00021 
00022 using namespace shogun;
00023 
00024 CSVMOcas::CSVMOcas(E_SVM_TYPE type)
00025 : CLinearClassifier(), use_bias(true), bufsize(3000), C1(1), C2(1),
00026     epsilon(1e-3), method(type)
00027 {
00028     w=NULL;
00029     old_w=NULL;
00030 }
00031 
00032 CSVMOcas::CSVMOcas(
00033     float64_t C, CDotFeatures* traindat, CLabels* trainlab)
00034 : CLinearClassifier(), use_bias(true), bufsize(3000), C1(C), C2(C),
00035     epsilon(1e-3)
00036 {
00037     w=NULL;
00038     old_w=NULL;
00039     method=SVM_OCAS;
00040     set_features(traindat);
00041     set_labels(trainlab);
00042 }
00043 
00044 
00045 CSVMOcas::~CSVMOcas()
00046 {
00047 }
00048 
00049 bool CSVMOcas::train(CFeatures* data)
00050 {
00051     SG_INFO("C=%f, epsilon=%f, bufsize=%d\n", get_C1(), get_epsilon(), bufsize);
00052     SG_DEBUG("use_bias = %i\n", get_bias_enabled()) ;
00053 
00054     ASSERT(labels);
00055     if (data)
00056     {
00057         if (!data->has_property(FP_DOT))
00058             SG_ERROR("Specified features are not of type CDotFeatures\n");
00059         set_features((CDotFeatures*) data);
00060     }
00061     ASSERT(features);
00062     ASSERT(labels->is_two_class_labeling());
00063 
00064     int32_t num_train_labels=0;
00065     lab=labels->get_labels(num_train_labels);
00066     w_dim=features->get_dim_feature_space();
00067     int32_t num_vec=features->get_num_vectors();
00068 
00069     ASSERT(num_vec==num_train_labels);
00070     ASSERT(num_vec>0);
00071 
00072     delete[] w;
00073     w=new float64_t[w_dim];
00074     memset(w, 0, w_dim*sizeof(float64_t));
00075 
00076     delete[] old_w;
00077     old_w=new float64_t[w_dim];
00078     memset(old_w, 0, w_dim*sizeof(float64_t));
00079     bias=0;
00080     old_bias=0;
00081 
00082     tmp_a_buf=new float64_t[w_dim];
00083     cp_value=new float64_t*[bufsize];
00084     cp_index=new uint32_t*[bufsize];
00085     cp_nz_dims=new uint32_t[bufsize];
00086     cp_bias=new float64_t[bufsize];
00087     memset(cp_bias, 0, sizeof(float64_t)*bufsize);
00088 
00089     float64_t TolAbs=0;
00090     float64_t QPBound=0;
00091     int32_t Method=0;
00092     if (method == SVM_OCAS)
00093         Method = 1;
00094     ocas_return_value_T result = svm_ocas_solver( get_C1(), num_vec, get_epsilon(),
00095             TolAbs, QPBound, bufsize, Method, 
00096             &CSVMOcas::compute_W,
00097             &CSVMOcas::update_W, 
00098             &CSVMOcas::add_new_cut, 
00099             &CSVMOcas::compute_output,
00100             &CSVMOcas::sort,
00101             this);
00102 
00103     SG_INFO("Ocas Converged after %d iterations\n"
00104             "==================================\n"
00105             "timing statistics:\n"
00106             "output_time: %f s\n"
00107             "sort_time: %f s\n"
00108             "add_time: %f s\n"
00109             "w_time: %f s\n"
00110             "solver_time %f s\n"
00111             "ocas_time %f s\n\n", result.nIter, result.output_time, result.sort_time,
00112             result.add_time, result.w_time, result.solver_time, result.ocas_time);
00113 
00114     delete[] tmp_a_buf;
00115 
00116     uint32_t num_cut_planes = result.nCutPlanes;
00117 
00118     SG_DEBUG("num_cut_planes=%d\n", num_cut_planes);
00119     for (uint32_t i=0; i<num_cut_planes; i++)
00120     {
00121         SG_DEBUG("cp_value[%d]=%p\n", i, cp_value);
00122         delete[] cp_value[i];
00123         SG_DEBUG("cp_index[%d]=%p\n", i, cp_index);
00124         delete[] cp_index[i];
00125     }
00126 
00127     delete[] cp_value;
00128     cp_value=NULL;
00129     delete[] cp_index;
00130     cp_index=NULL;
00131     delete[] cp_nz_dims;
00132     cp_nz_dims=NULL;
00133     delete[] cp_bias;
00134     cp_bias=NULL;
00135 
00136     delete[] lab;
00137     lab=NULL;
00138 
00139     delete[] old_w;
00140     old_w=NULL;
00141 
00142     return true;
00143 }
00144 
00145 /*----------------------------------------------------------------------------------
00146   sq_norm_W = sparse_update_W( t ) does the following:
00147 
00148   W = oldW*(1-t) + t*W;
00149   sq_norm_W = W'*W;
00150 
00151   ---------------------------------------------------------------------------------*/
00152 float64_t CSVMOcas::update_W( float64_t t, void* ptr )
00153 {
00154   float64_t sq_norm_W = 0;         
00155   CSVMOcas* o = (CSVMOcas*) ptr;
00156   uint32_t nDim = (uint32_t) o->w_dim;
00157   float64_t* W=o->w;
00158   float64_t* oldW=o->old_w;
00159 
00160   for(uint32_t j=0; j <nDim; j++)
00161   {
00162       W[j] = oldW[j]*(1-t) + t*W[j];
00163       sq_norm_W += W[j]*W[j];
00164   }          
00165   o->bias=o->old_bias*(1-t) + t*o->bias;
00166   sq_norm_W += CMath::sq(o->bias);
00167 
00168   return( sq_norm_W );
00169 }
00170 
00171 /*----------------------------------------------------------------------------------
00172   sparse_add_new_cut( new_col_H, new_cut, cut_length, nSel ) does the following:
00173 
00174     new_a = sum(data_X(:,find(new_cut ~=0 )),2);
00175     new_col_H = [sparse_A(:,1:nSel)'*new_a ; new_a'*new_a];
00176     sparse_A(:,nSel+1) = new_a;
00177 
00178   ---------------------------------------------------------------------------------*/
00179 void CSVMOcas::add_new_cut(
00180     float64_t *new_col_H, uint32_t *new_cut, uint32_t cut_length,
00181     uint32_t nSel, void* ptr)
00182 {
00183     CSVMOcas* o = (CSVMOcas*) ptr;
00184     CDotFeatures* f = o->features;
00185     uint32_t nDim=(uint32_t) o->w_dim;
00186     float64_t* y = o->lab;
00187 
00188     float64_t** c_val = o->cp_value;
00189     uint32_t** c_idx = o->cp_index;
00190     uint32_t* c_nzd = o->cp_nz_dims;
00191     float64_t* c_bias = o->cp_bias;
00192 
00193     float64_t sq_norm_a;
00194     uint32_t i, j, nz_dims;
00195 
00196     /* temporary vector */
00197     float64_t* new_a = o->tmp_a_buf;
00198     memset(new_a, 0, sizeof(float64_t)*nDim);
00199 
00200     for(i=0; i < cut_length; i++) 
00201     {
00202         f->add_to_dense_vec(y[new_cut[i]], new_cut[i], new_a, nDim);
00203 
00204         if (o->use_bias)
00205             c_bias[nSel]+=y[new_cut[i]];
00206     }
00207 
00208     /* compute new_a'*new_a and count number of non-zerou dimensions */
00209     nz_dims = 0; 
00210     sq_norm_a = CMath::sq(c_bias[nSel]);
00211     for(j=0; j < nDim; j++ ) {
00212         if(new_a[j] != 0) {
00213             nz_dims++;
00214             sq_norm_a += new_a[j]*new_a[j];
00215         }
00216     }
00217 
00218     /* sparsify new_a and insert it to the last column of sparse_A */
00219     c_nzd[nSel] = nz_dims;
00220     c_idx[nSel]=NULL;
00221     c_val[nSel]=NULL;
00222 
00223     if(nz_dims > 0)
00224     {
00225         c_idx[nSel]=new uint32_t[nz_dims];
00226         c_val[nSel]=new float64_t[nz_dims];
00227 
00228         uint32_t idx=0;
00229         for(j=0; j < nDim; j++ )
00230         {
00231             if(new_a[j] != 0)
00232             {
00233                 c_idx[nSel][idx] = j;
00234                 c_val[nSel][idx++] = new_a[j];
00235             }
00236         }
00237     }
00238 
00239     new_col_H[nSel] = sq_norm_a;
00240 
00241     for(i=0; i < nSel; i++)
00242     {
00243         float64_t tmp = c_bias[nSel]*c_bias[i];
00244         for(j=0; j < c_nzd[i]; j++)
00245             tmp += new_a[c_idx[i][j]]*c_val[i][j];
00246 
00247         new_col_H[i] = tmp;
00248     }
00249     //CMath::display_vector(new_col_H, nSel+1, "new_col_H");
00250     //CMath::display_vector((int32_t*) c_idx[nSel], (int32_t) nz_dims, "c_idx");
00251     //CMath::display_vector((float64_t*) c_val[nSel], nz_dims, "c_val");
00252 }
00253 
00254 void CSVMOcas::sort(float64_t* vals, uint32_t* idx, uint32_t size)
00255 {
00256     CMath::qsort_index(vals, idx, size);
00257 }
00258 
00259 /*----------------------------------------------------------------------
00260   sparse_compute_output( output ) does the follwing:
00261 
00262   output = data_X'*W;
00263   ----------------------------------------------------------------------*/
00264 void CSVMOcas::compute_output(float64_t *output, void* ptr)
00265 {
00266     CSVMOcas* o = (CSVMOcas*) ptr;
00267     CDotFeatures* f=o->features;
00268     int32_t nData=f->get_num_vectors();
00269 
00270     float64_t* y = o->lab;
00271 
00272     f->dense_dot_range(output, 0, nData, y, o->w, o->w_dim, 0.0);
00273 
00274     for (int32_t i=0; i<nData; i++)
00275         output[i]+=y[i]*o->bias;
00276     //CMath::display_vector(o->w, o->w_dim, "w");
00277     //CMath::display_vector(output, nData, "out");
00278 }
00279 
00280 /*----------------------------------------------------------------------
00281   sq_norm_W = compute_W( alpha, nSel ) does the following:
00282 
00283   oldW = W;
00284   W = sparse_A(:,1:nSel)'*alpha;
00285   sq_norm_W = W'*W;
00286   dp_WoldW = W'*oldW';
00287 
00288   ----------------------------------------------------------------------*/
00289 void CSVMOcas::compute_W(
00290     float64_t *sq_norm_W, float64_t *dp_WoldW, float64_t *alpha, uint32_t nSel,
00291     void* ptr )
00292 {
00293     CSVMOcas* o = (CSVMOcas*) ptr;
00294     uint32_t nDim= (uint32_t) o->w_dim;
00295     CMath::swap(o->w, o->old_w);
00296     float64_t* W=o->w;
00297     float64_t* oldW=o->old_w;
00298     memset(W, 0, sizeof(float64_t)*nDim);
00299     float64_t old_bias=o->bias;
00300     float64_t bias=0;
00301 
00302     float64_t** c_val = o->cp_value;
00303     uint32_t** c_idx = o->cp_index;
00304     uint32_t* c_nzd = o->cp_nz_dims;
00305     float64_t* c_bias = o->cp_bias;
00306 
00307     for(uint32_t i=0; i<nSel; i++)
00308     {
00309         uint32_t nz_dims = c_nzd[i];
00310 
00311         if(nz_dims > 0 && alpha[i] > 0)
00312         {
00313             for(uint32_t j=0; j < nz_dims; j++)
00314                 W[c_idx[i][j]] += alpha[i]*c_val[i][j];
00315         }
00316         bias += c_bias[i]*alpha[i];
00317     }
00318 
00319     *sq_norm_W = CMath::dot(W,W, nDim) + CMath::sq(bias);
00320     *dp_WoldW = CMath::dot(W,oldW, nDim) + bias*old_bias;
00321     //SG_PRINT("nSel=%d sq_norm_W=%f dp_WoldW=%f\n", nSel, *sq_norm_W, *dp_WoldW);
00322     
00323     o->bias = bias;
00324     o->old_bias = old_bias;
00325 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation