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-2008 Soeren Sonnenburg
00009  * Copyright (C) 2007-2008 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/SparseLinearClassifier.h"
00017 #include "classifier/svm/SVMOcas.h"
00018 #include "classifier/svm/libocas.h"
00019 #include "features/SparseFeatures.h"
00020 #include "features/Labels.h"
00021 
00022 CSVMOcas::CSVMOcas(E_SVM_TYPE type)
00023 : CSparseLinearClassifier(), use_bias(false), bufsize(3000), C1(1), C2(1),
00024     epsilon(1e-3), method(type)
00025 {
00026     w=NULL;
00027     old_w=NULL;
00028 }
00029 
00030 CSVMOcas::CSVMOcas(
00031     float64_t C, CSparseFeatures<float64_t>* traindat, CLabels* trainlab)
00032 : CSparseLinearClassifier(), use_bias(false), bufsize(3000), C1(C), C2(C),
00033     epsilon(1e-3)
00034 {
00035     w=NULL;
00036     old_w=NULL;
00037     method=SVM_OCAS;
00038     set_features(traindat);
00039     set_labels(trainlab);
00040 }
00041 
00042 
00043 CSVMOcas::~CSVMOcas()
00044 {
00045 }
00046 
00047 bool CSVMOcas::train()
00048 {
00049     SG_INFO("C=%f, epsilon=%f, bufsize=%d\n", get_C1(), get_epsilon(), bufsize);
00050 
00051     ASSERT(labels);
00052     ASSERT(get_features());
00053     ASSERT(labels->is_two_class_labeling());
00054 
00055     int32_t num_train_labels=0;
00056     lab=labels->get_labels(num_train_labels);
00057     w_dim=features->get_num_features();
00058     int32_t num_vec=features->get_num_vectors();
00059 
00060     ASSERT(num_vec==num_train_labels);
00061     ASSERT(num_vec>0);
00062 
00063     delete[] w;
00064     w=new float64_t[w_dim];
00065     memset(w, 0, w_dim*sizeof(float64_t));
00066 
00067     delete[] old_w;
00068     old_w=new float64_t[w_dim];
00069     memset(old_w, 0, w_dim*sizeof(float64_t));
00070     bias=0;
00071 
00072     tmp_a_buf=new float64_t[w_dim];
00073     cp_value=new float64_t*[bufsize];
00074     cp_index=new uint32_t*[bufsize];
00075     cp_nz_dims=new uint32_t[bufsize];
00076 
00077     float64_t TolAbs=0;
00078     float64_t QPBound=0;
00079     int32_t Method=0;
00080     if (method == SVM_OCAS)
00081         Method = 1;
00082     ocas_return_value_T result = svm_ocas_solver( get_C1(), num_vec, get_epsilon(),
00083             TolAbs, QPBound, bufsize, Method, 
00084             &CSVMOcas::compute_W,
00085             &CSVMOcas::update_W, 
00086             &CSVMOcas::add_new_cut, 
00087             &CSVMOcas::compute_output,
00088             &CSVMOcas::sort,
00089             this);
00090 
00091     SG_INFO("Ocas Converged after %d iterations\n"
00092             "==================================\n"
00093             "timing statistics:\n"
00094             "output_time: %f s\n"
00095             "sort_time: %f s\n"
00096             "add_time: %f s\n"
00097             "w_time: %f s\n"
00098             "solver_time %f s\n"
00099             "ocas_time %f s\n\n", result.nIter, result.output_time, result.sort_time,
00100             result.add_time, result.w_time, result.solver_time, result.ocas_time);
00101 
00102     delete[] tmp_a_buf;
00103 
00104     uint32_t num_cut_planes = result.nCutPlanes;
00105 
00106     for (uint32_t i=0; i<num_cut_planes; i++)
00107     {
00108         delete[] cp_value[i];
00109         delete[] cp_index[i];
00110     }
00111 
00112     delete[] cp_value;
00113     cp_value=NULL;
00114     delete[] cp_index;
00115     cp_index=NULL;
00116     delete[] cp_nz_dims;
00117     cp_nz_dims=NULL;
00118 
00119     delete[] lab;
00120     lab=NULL;
00121 
00122     delete[] old_w;
00123     old_w=NULL;
00124 
00125     return true;
00126 }
00127 
00128 /*----------------------------------------------------------------------------------
00129   sq_norm_W = sparse_update_W( t ) does the following:
00130 
00131   W = oldW*(1-t) + t*W;
00132   sq_norm_W = W'*W;
00133 
00134   ---------------------------------------------------------------------------------*/
00135 float64_t CSVMOcas::update_W( float64_t t, void* ptr )
00136 {
00137   float64_t sq_norm_W = 0;         
00138   CSVMOcas* o = (CSVMOcas*) ptr;
00139   uint32_t nDim = (uint32_t) o->w_dim;
00140   float64_t* W=o->w;
00141   float64_t* oldW=o->old_w;
00142 
00143   for(uint32_t j=0; j <nDim; j++)
00144   {
00145       W[j] = oldW[j]*(1-t) + t*W[j];
00146       sq_norm_W += W[j]*W[j];
00147   }          
00148 
00149   return( sq_norm_W );
00150 }
00151 
00152 /*----------------------------------------------------------------------------------
00153   sparse_add_new_cut( new_col_H, new_cut, cut_length, nSel ) does the following:
00154 
00155     new_a = sum(data_X(:,find(new_cut ~=0 )),2);
00156     new_col_H = [sparse_A(:,1:nSel)'*new_a ; new_a'*new_a];
00157     sparse_A(:,nSel+1) = new_a;
00158 
00159   ---------------------------------------------------------------------------------*/
00160 void CSVMOcas::add_new_cut(
00161     float64_t *new_col_H, uint32_t *new_cut, uint32_t cut_length,
00162     uint32_t nSel, void* ptr)
00163 {
00164     CSVMOcas* o = (CSVMOcas*) ptr;
00165     CSparseFeatures<float64_t>* f = o->get_features();
00166     uint32_t nDim=(uint32_t) o->w_dim;
00167     float64_t* y = o->lab;
00168 
00169     float64_t** c_val = o->cp_value;
00170     uint32_t** c_idx = o->cp_index;
00171     uint32_t* c_nzd = o->cp_nz_dims;
00172 
00173     float64_t sq_norm_a;
00174     uint32_t i, j, nz_dims;
00175 
00176     /* temporary vector */
00177     float64_t* new_a = o->tmp_a_buf;
00178     memset(new_a, 0, sizeof(float64_t)*nDim);
00179 
00180     for(i=0; i < cut_length; i++) 
00181         f->add_to_dense_vec(y[new_cut[i]], new_cut[i], new_a, nDim);
00182 
00183     /* compute new_a'*new_a and count number of non-zerou dimensions */
00184     nz_dims = 0; 
00185     sq_norm_a = 0;
00186     for(j=0; j < nDim; j++ ) {
00187         if(new_a[j] != 0) {
00188             nz_dims++;
00189             sq_norm_a += new_a[j]*new_a[j];
00190         }
00191     }
00192 
00193     /* sparsify new_a and insert it to the last column of sparse_A */
00194     c_nzd[nSel] = nz_dims;
00195     if(nz_dims > 0)
00196     {
00197         c_idx[nSel]=new uint32_t[nz_dims];
00198         c_val[nSel]=new float64_t[nz_dims];
00199 
00200         uint32_t idx=0;
00201         for(j=0; j < nDim; j++ )
00202         {
00203             if(new_a[j] != 0)
00204             {
00205                 c_idx[nSel][idx] = j;
00206                 c_val[nSel][idx++] = new_a[j];
00207             }
00208         }
00209     }
00210 
00211     new_col_H[nSel] = sq_norm_a;
00212 
00213     for(i=0; i < nSel; i++)
00214     {
00215         float64_t tmp = 0;
00216         for(j=0; j < c_nzd[i]; j++)
00217             tmp += new_a[c_idx[i][j]]*c_val[i][j];
00218 
00219         new_col_H[i] = tmp;
00220     }
00221     //CMath::display_vector(new_col_H, nSel+1, "new_col_H");
00222     //CMath::display_vector((int32_t*) c_idx[nSel], (int32_t) nz_dims, "c_idx");
00223     //CMath::display_vector((float64_t*) c_val[nSel], nz_dims, "c_val");
00224 }
00225 
00226 void CSVMOcas::sort(float64_t* vals, uint32_t* idx, uint32_t size)
00227 {
00228     CMath::qsort_index(vals, idx, size);
00229 }
00230 
00231 /*----------------------------------------------------------------------
00232   sparse_compute_output( output ) does the follwing:
00233 
00234   output = data_X'*W;
00235   ----------------------------------------------------------------------*/
00236 void CSVMOcas::compute_output(float64_t *output, void* ptr)
00237 {
00238     CSVMOcas* o = (CSVMOcas*) ptr;
00239     CSparseFeatures<float64_t>* f=o->get_features();
00240     int32_t nData=f->get_num_vectors();
00241 
00242     float64_t* y = o->lab;
00243     f->dense_dot_range(output, 0, nData, y, o->w, o->w_dim, 0.0);
00244     //CMath::display_vector(o->w, o->w_dim, "w");
00245     //CMath::display_vector(output, nData, "out");
00246 }
00247 
00248 /*----------------------------------------------------------------------
00249   sq_norm_W = compute_W( alpha, nSel ) does the following:
00250 
00251   oldW = W;
00252   W = sparse_A(:,1:nSel)'*alpha;
00253   sq_norm_W = W'*W;
00254   dp_WoldW = W'*oldW';
00255 
00256   ----------------------------------------------------------------------*/
00257 void CSVMOcas::compute_W(
00258     float64_t *sq_norm_W, float64_t *dp_WoldW, float64_t *alpha, uint32_t nSel,
00259     void* ptr )
00260 {
00261     CSVMOcas* o = (CSVMOcas*) ptr;
00262     uint32_t nDim= (uint32_t) o->w_dim;
00263     CMath::swap(o->w, o->old_w);
00264     float64_t* W=o->w;
00265     float64_t* oldW=o->old_w;
00266     memset(W, 0, sizeof(float64_t)*nDim);
00267 
00268     float64_t** c_val = o->cp_value;
00269     uint32_t** c_idx = o->cp_index;
00270     uint32_t* c_nzd = o->cp_nz_dims;
00271 
00272     for(uint32_t i=0; i<nSel; i++)
00273     {
00274         uint32_t nz_dims = c_nzd[i];
00275 
00276         if(nz_dims > 0 && alpha[i] > 0)
00277         {
00278             for(uint32_t j=0; j < nz_dims; j++)
00279                 W[c_idx[i][j]] += alpha[i]*c_val[i][j];
00280         }
00281     }
00282 
00283     *sq_norm_W = CMath::dot(W,W, nDim);
00284     *dp_WoldW = CMath::dot(W,oldW, nDim);;
00285     //SG_PRINT("nSel=%d sq_norm_W=%f dp_WoldW=%f\n", nSel, *sq_norm_W, *dp_WoldW);
00286 }

SHOGUN Machine Learning Toolbox - Documentation