SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SVRLight.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) 1999-2009 Soeren Sonnenburg
00008  * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society
00009  */
00010 
00011 #include <shogun/lib/config.h>
00012 
00013 #ifdef USE_SVMLIGHT
00014 
00015 #include <shogun/io/SGIO.h>
00016 #include <shogun/mathematics/lapack.h>
00017 #include <shogun/lib/Signal.h>
00018 #include <shogun/mathematics/Math.h>
00019 #include <shogun/regression/svr/SVRLight.h>
00020 #include <shogun/machine/KernelMachine.h>
00021 #include <shogun/kernel/CombinedKernel.h>
00022 #include <shogun/labels/RegressionLabels.h>
00023 
00024 #include <unistd.h>
00025 
00026 #ifdef USE_CPLEX
00027 extern "C" {
00028 #include <ilcplex/cplex.h>
00029 }
00030 #endif
00031 
00032 #include <shogun/base/Parallel.h>
00033 
00034 #ifdef HAVE_PTHREAD
00035 #include <pthread.h>
00036 #endif
00037 
00038 using namespace shogun;
00039 
00040 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00041 struct S_THREAD_PARAM_SVRLIGHT
00042 {
00043     float64_t* lin;
00044     int32_t start, end;
00045     int32_t* active2dnum;
00046     int32_t* docs;
00047     CKernel* kernel;
00048     int32_t num_vectors;
00049 };
00050 #endif // DOXYGEN_SHOULD_SKIP_THIS
00051 
00052 CSVRLight::CSVRLight(float64_t C, float64_t eps, CKernel* k, CLabels* lab)
00053 : CSVMLight(C, k, lab)
00054 {
00055     set_tube_epsilon(eps);
00056 }
00057 
00058 CSVRLight::CSVRLight()
00059 : CSVMLight()
00060 {
00061 }
00062 
00064 CSVRLight::~CSVRLight()
00065 {
00066 }
00067 
00068 EMachineType CSVRLight::get_classifier_type()
00069 {
00070     return CT_SVRLIGHT;
00071 }
00072 
00073 bool CSVRLight::train_machine(CFeatures* data)
00074 {
00075     //certain setup params
00076     verbosity=1;
00077     init_margin=0.15;
00078     init_iter=500;
00079     precision_violations=0;
00080     opt_precision=DEF_PRECISION;
00081 
00082     strcpy (learn_parm->predfile, "");
00083     learn_parm->biased_hyperplane=1;
00084     learn_parm->sharedslack=0;
00085     learn_parm->remove_inconsistent=0;
00086     learn_parm->skip_final_opt_check=1;
00087     learn_parm->svm_maxqpsize=get_qpsize();
00088     learn_parm->svm_newvarsinqp=learn_parm->svm_maxqpsize-1;
00089     learn_parm->maxiter=100000;
00090     learn_parm->svm_iter_to_shrink=100;
00091     learn_parm->svm_c=get_C1();
00092     learn_parm->transduction_posratio=0.33;
00093     learn_parm->svm_costratio=get_C2()/get_C1();
00094     learn_parm->svm_costratio_unlab=1.0;
00095     learn_parm->svm_unlabbound=1E-5;
00096     learn_parm->epsilon_crit=epsilon; // GU: better decrease it ... ??
00097     learn_parm->epsilon_a=1E-15;
00098     learn_parm->compute_loo=0;
00099     learn_parm->rho=1.0;
00100     learn_parm->xa_depth=0;
00101 
00102     if (!kernel)
00103     {
00104         SG_ERROR("SVR_light can not proceed without kernel!\n")
00105         return false ;
00106     }
00107 
00108     if (!m_labels)
00109     {
00110         SG_ERROR("SVR_light can not proceed without labels!\n")
00111         return false;
00112     }
00113 
00114     if (data)
00115     {
00116         if (m_labels->get_num_labels() != data->get_num_vectors())
00117             SG_ERROR("Number of training vectors does not match number of labels\n")
00118         kernel->init(data, data);
00119     }
00120 
00121     if (kernel->has_property(KP_LINADD) && get_linadd_enabled())
00122         kernel->clear_normal();
00123 
00124     // output some info
00125     SG_DEBUG("qpsize = %i\n", learn_parm->svm_maxqpsize)
00126     SG_DEBUG("epsilon = %1.1e\n", learn_parm->epsilon_crit)
00127     SG_DEBUG("kernel->has_property(KP_LINADD) = %i\n", kernel->has_property(KP_LINADD))
00128     SG_DEBUG("kernel->has_property(KP_KERNCOMBINATION) = %i\n", kernel->has_property(KP_KERNCOMBINATION))
00129     SG_DEBUG("get_linadd_enabled() = %i\n", get_linadd_enabled())
00130     SG_DEBUG("kernel->get_num_subkernels() = %i\n", kernel->get_num_subkernels())
00131 
00132     use_kernel_cache = !((kernel->get_kernel_type() == K_CUSTOM) ||
00133                          (get_linadd_enabled() && kernel->has_property(KP_LINADD)));
00134 
00135     SG_DEBUG("use_kernel_cache = %i\n", use_kernel_cache)
00136 
00137     // train the svm
00138     svr_learn();
00139 
00140     // brain damaged svm light work around
00141     create_new_model(model->sv_num-1);
00142     set_bias(-model->b);
00143     for (int32_t i=0; i<model->sv_num-1; i++)
00144     {
00145         set_alpha(i, model->alpha[i+1]);
00146         set_support_vector(i, model->supvec[i+1]);
00147     }
00148 
00149     if (kernel->has_property(KP_LINADD) && get_linadd_enabled())
00150         kernel->clear_normal() ;
00151 
00152     return true ;
00153 }
00154 
00155 void CSVRLight::svr_learn()
00156 {
00157     int32_t *inconsistent, i, j;
00158     int32_t upsupvecnum;
00159     float64_t maxdiff, *lin, *c, *a;
00160     int32_t iterations;
00161     float64_t *xi_fullset; /* buffer for storing xi on full sample in loo */
00162     float64_t *a_fullset;  /* buffer for storing alpha on full sample in loo */
00163     TIMING timing_profile;
00164     SHRINK_STATE shrink_state;
00165     int32_t* label;
00166     int32_t* docs;
00167 
00168     ASSERT(m_labels)
00169     int32_t totdoc=m_labels->get_num_labels();
00170     num_vectors=totdoc;
00171 
00172     // set up regression problem in standard form
00173     docs=SG_MALLOC(int32_t, 2*totdoc);
00174     label=SG_MALLOC(int32_t, 2*totdoc);
00175     c = SG_MALLOC(float64_t, 2*totdoc);
00176 
00177   for(i=0;i<totdoc;i++) {
00178       docs[i]=i;
00179       j=2*totdoc-1-i;
00180       label[i]=+1;
00181       c[i]=((CRegressionLabels*) m_labels)->get_label(i);
00182       docs[j]=j;
00183       label[j]=-1;
00184       c[j]=((CRegressionLabels*) m_labels)->get_label(i);
00185   }
00186   totdoc*=2;
00187 
00188   //prepare kernel cache for regression (i.e. cachelines are twice of current size)
00189   kernel->resize_kernel_cache( kernel->get_cache_size(), true);
00190 
00191   if (kernel->get_kernel_type() == K_COMBINED)
00192   {
00193       CCombinedKernel* k = (CCombinedKernel*) kernel;
00194 
00195       for (index_t k_idx=0; k_idx<k->get_num_kernels(); k_idx++)
00196       {
00197           CKernel* kn = k->get_kernel(k_idx);
00198           kn->resize_kernel_cache( kernel->get_cache_size(), true);
00199           SG_UNREF(kn);
00200       }
00201   }
00202 
00203   timing_profile.time_kernel=0;
00204   timing_profile.time_opti=0;
00205   timing_profile.time_shrink=0;
00206   timing_profile.time_update=0;
00207   timing_profile.time_model=0;
00208   timing_profile.time_check=0;
00209   timing_profile.time_select=0;
00210 
00211     SG_FREE(W);
00212     W=NULL;
00213 
00214     if (kernel->has_property(KP_KERNCOMBINATION) && callback)
00215     {
00216         W = SG_MALLOC(float64_t, totdoc*kernel->get_num_subkernels());
00217         for (i=0; i<totdoc*kernel->get_num_subkernels(); i++)
00218             W[i]=0;
00219     }
00220 
00221     /* make sure -n value is reasonable */
00222     if((learn_parm->svm_newvarsinqp < 2)
00223             || (learn_parm->svm_newvarsinqp > learn_parm->svm_maxqpsize)) {
00224         learn_parm->svm_newvarsinqp=learn_parm->svm_maxqpsize;
00225     }
00226 
00227     init_shrink_state(&shrink_state,totdoc,(int32_t)MAXSHRINK);
00228 
00229     inconsistent = SG_MALLOC(int32_t, totdoc);
00230     a = SG_MALLOC(float64_t, totdoc);
00231     a_fullset = SG_MALLOC(float64_t, totdoc);
00232     xi_fullset = SG_MALLOC(float64_t, totdoc);
00233     lin = SG_MALLOC(float64_t, totdoc);
00234     learn_parm->svm_cost = SG_MALLOC(float64_t, totdoc);
00235     if (m_linear_term.vlen>0)
00236         learn_parm->eps=get_linear_term_array();
00237     else
00238     {
00239         learn_parm->eps=SG_MALLOC(float64_t, totdoc);      /* equivalent regression epsilon for classification */
00240         SGVector<float64_t>::fill_vector(learn_parm->eps, totdoc, tube_epsilon);
00241     }
00242 
00243     SG_FREE(model->supvec);
00244     SG_FREE(model->alpha);
00245     SG_FREE(model->index);
00246     model->supvec = SG_MALLOC(int32_t, totdoc+2);
00247     model->alpha = SG_MALLOC(float64_t, totdoc+2);
00248     model->index = SG_MALLOC(int32_t, totdoc+2);
00249 
00250     model->at_upper_bound=0;
00251     model->b=0;
00252     model->supvec[0]=0;  /* element 0 reserved and empty for now */
00253     model->alpha[0]=0;
00254     model->totdoc=totdoc;
00255 
00256     model->kernel=kernel;
00257 
00258     model->sv_num=1;
00259     model->loo_error=-1;
00260     model->loo_recall=-1;
00261     model->loo_precision=-1;
00262     model->xa_error=-1;
00263     model->xa_recall=-1;
00264     model->xa_precision=-1;
00265 
00266   for(i=0;i<totdoc;i++) {    /* various inits */
00267     inconsistent[i]=0;
00268     a[i]=0;
00269     lin[i]=0;
00270 
00271         if(label[i] > 0) {
00272             learn_parm->svm_cost[i]=learn_parm->svm_c*learn_parm->svm_costratio*
00273                 fabs((float64_t)label[i]);
00274         }
00275         else if(label[i] < 0) {
00276             learn_parm->svm_cost[i]=learn_parm->svm_c*fabs((float64_t)label[i]);
00277         }
00278         else
00279             ASSERT(false)
00280     }
00281 
00282     if(verbosity==1) {
00283         SG_DEBUG("Optimizing...\n")
00284     }
00285 
00286     /* train the svm */
00287         SG_DEBUG("num_train: %d\n", totdoc)
00288   iterations=optimize_to_convergence(docs,label,totdoc,
00289                      &shrink_state,inconsistent,a,lin,
00290                      c,&timing_profile,
00291                      &maxdiff,(int32_t)-1,
00292                      (int32_t)1);
00293 
00294 
00295     if(verbosity>=1) {
00296         SG_DONE()
00297         SG_INFO("(%ld iterations)\n",iterations)
00298         SG_INFO("Optimization finished (maxdiff=%.8f).\n",maxdiff)
00299         SG_INFO("obj = %.16f, rho = %.16f\n",get_objective(),model->b)
00300 
00301         upsupvecnum=0;
00302 
00303         SG_DEBUG("num sv: %d\n", model->sv_num)
00304         for(i=1;i<model->sv_num;i++)
00305         {
00306             if(fabs(model->alpha[i]) >=
00307                     (learn_parm->svm_cost[model->supvec[i]]-
00308                      learn_parm->epsilon_a))
00309                 upsupvecnum++;
00310         }
00311         SG_INFO("Number of SV: %d (including %d at upper bound)\n",
00312                 model->sv_num-1,upsupvecnum);
00313     }
00314 
00315   /* this makes sure the model we return does not contain pointers to the
00316      temporary documents */
00317   for(i=1;i<model->sv_num;i++) {
00318     j=model->supvec[i];
00319     if(j >= (totdoc/2)) {
00320       j=totdoc-j-1;
00321     }
00322     model->supvec[i]=j;
00323   }
00324 
00325   shrink_state_cleanup(&shrink_state);
00326     SG_FREE(label);
00327     SG_FREE(inconsistent);
00328     SG_FREE(c);
00329     SG_FREE(a);
00330     SG_FREE(a_fullset);
00331     SG_FREE(xi_fullset);
00332     SG_FREE(lin);
00333     SG_FREE(learn_parm->svm_cost);
00334     SG_FREE(docs);
00335 }
00336 
00337 float64_t CSVRLight::compute_objective_function(
00338     float64_t *a, float64_t *lin, float64_t *c, float64_t* eps, int32_t *label,
00339     int32_t totdoc)
00340 {
00341   /* calculate value of objective function */
00342   float64_t criterion=0;
00343 
00344   for(int32_t i=0;i<totdoc;i++)
00345       criterion+=(eps[i]-(float64_t)label[i]*c[i])*a[i]+0.5*a[i]*label[i]*lin[i];
00346 
00347   /* float64_t check=0;
00348   for(int32_t i=0;i<totdoc;i++)
00349   {
00350       check+=a[i]*eps-a[i]*label[i]*c[i];
00351       for(int32_t j=0;j<totdoc;j++)
00352           check+= 0.5*a[i]*label[i]*a[j]*label[j]*compute_kernel(i,j);
00353 
00354   }
00355 
00356   SG_INFO("REGRESSION OBJECTIVE %f vs. CHECK %f (diff %f)\n", criterion, check, criterion-check) */
00357 
00358   return(criterion);
00359 }
00360 
00361 void* CSVRLight::update_linear_component_linadd_helper(void *params_)
00362 {
00363     S_THREAD_PARAM_SVRLIGHT * params = (S_THREAD_PARAM_SVRLIGHT*) params_ ;
00364 
00365     int32_t jj=0, j=0 ;
00366 
00367     for(jj=params->start;(jj<params->end) && (j=params->active2dnum[jj])>=0;jj++)
00368         params->lin[j]+=params->kernel->compute_optimized(CSVRLight::regression_fix_index2(params->docs[j], params->num_vectors));
00369 
00370     return NULL ;
00371 }
00372 
00373 int32_t CSVRLight::regression_fix_index(int32_t i)
00374 {
00375     if (i>=num_vectors)
00376         i=2*num_vectors-1-i;
00377 
00378     return i;
00379 }
00380 
00381 int32_t CSVRLight::regression_fix_index2(
00382         int32_t i, int32_t num_vectors)
00383 {
00384     if (i>=num_vectors)
00385         i=2*num_vectors-1-i;
00386 
00387     return i;
00388 }
00389 
00390 float64_t CSVRLight::compute_kernel(int32_t i, int32_t j)
00391 {
00392     i=regression_fix_index(i);
00393     j=regression_fix_index(j);
00394     return kernel->kernel(i, j);
00395 }
00396 
00397 void CSVRLight::update_linear_component(
00398     int32_t* docs, int32_t* label, int32_t *active2dnum, float64_t *a,
00399     float64_t *a_old, int32_t *working2dnum, int32_t totdoc, float64_t *lin,
00400     float64_t *aicache, float64_t* c)
00401      /* keep track of the linear component */
00402      /* lin of the gradient etc. by updating */
00403      /* based on the change of the variables */
00404      /* in the current working set */
00405 {
00406     register int32_t i=0,ii=0,j=0,jj=0;
00407 
00408     if (kernel->has_property(KP_LINADD) && get_linadd_enabled())
00409     {
00410         if (callback)
00411         {
00412             update_linear_component_mkl_linadd(docs, label, active2dnum, a, a_old, working2dnum,
00413                                                totdoc,  lin, aicache, c) ;
00414         }
00415         else
00416         {
00417             kernel->clear_normal();
00418 
00419             int32_t num_working=0;
00420             for(ii=0;(i=working2dnum[ii])>=0;ii++) {
00421                 if(a[i] != a_old[i]) {
00422                     kernel->add_to_normal(regression_fix_index(docs[i]), (a[i]-a_old[i])*(float64_t)label[i]);
00423                     num_working++;
00424                 }
00425             }
00426 
00427             if (num_working>0)
00428             {
00429                 if (parallel->get_num_threads() < 2)
00430                 {
00431                     for(jj=0;(j=active2dnum[jj])>=0;jj++) {
00432                         lin[j]+=kernel->compute_optimized(regression_fix_index(docs[j]));
00433                     }
00434                 }
00435 #ifdef HAVE_PTHREAD
00436                 else
00437                 {
00438                     int32_t num_elem = 0 ;
00439                     for(jj=0;(j=active2dnum[jj])>=0;jj++) num_elem++ ;
00440 
00441                     pthread_t* threads = SG_MALLOC(pthread_t, parallel->get_num_threads()-1);
00442                     S_THREAD_PARAM_SVRLIGHT* params = SG_MALLOC(S_THREAD_PARAM_SVRLIGHT, parallel->get_num_threads()-1);
00443                     int32_t start = 0 ;
00444                     int32_t step = num_elem/parallel->get_num_threads() ;
00445                     int32_t end = step ;
00446 
00447                     for (int32_t t=0; t<parallel->get_num_threads()-1; t++)
00448                     {
00449                         params[t].kernel = kernel ;
00450                         params[t].lin = lin ;
00451                         params[t].docs = docs ;
00452                         params[t].active2dnum=active2dnum ;
00453                         params[t].start = start ;
00454                         params[t].end = end ;
00455                         params[t].num_vectors=num_vectors ;
00456 
00457                         start=end ;
00458                         end+=step ;
00459                         pthread_create(&threads[t], NULL, update_linear_component_linadd_helper, (void*)&params[t]) ;
00460                     }
00461 
00462                     for(jj=params[parallel->get_num_threads()-2].end;(j=active2dnum[jj])>=0;jj++) {
00463                         lin[j]+=kernel->compute_optimized(regression_fix_index(docs[j]));
00464                     }
00465                     void* ret;
00466                     for (int32_t t=0; t<parallel->get_num_threads()-1; t++)
00467                         pthread_join(threads[t], &ret) ;
00468 
00469                     SG_FREE(params);
00470                     SG_FREE(threads);
00471                 }
00472 #endif
00473             }
00474         }
00475     }
00476     else
00477     {
00478         if (callback)
00479         {
00480             update_linear_component_mkl(docs, label, active2dnum,
00481                     a, a_old, working2dnum, totdoc, lin, aicache, c) ;
00482         }
00483         else {
00484             for(jj=0;(i=working2dnum[jj])>=0;jj++) {
00485                 if(a[i] != a_old[i]) {
00486                     kernel->get_kernel_row(i,active2dnum,aicache);
00487                     for(ii=0;(j=active2dnum[ii])>=0;ii++)
00488                         lin[j]+=(a[i]-a_old[i])*aicache[j]*(float64_t)label[i];
00489                 }
00490             }
00491         }
00492     }
00493 }
00494 
00495 void CSVRLight::update_linear_component_mkl(
00496     int32_t* docs, int32_t* label, int32_t *active2dnum, float64_t *a,
00497     float64_t *a_old, int32_t *working2dnum, int32_t totdoc, float64_t *lin,
00498     float64_t *aicache, float64_t* c)
00499 {
00500     int32_t num         = totdoc;
00501     int32_t num_weights = -1;
00502     int32_t num_kernels = kernel->get_num_subkernels() ;
00503     const float64_t* old_beta  = kernel->get_subkernel_weights(num_weights);
00504 
00505     ASSERT(num_weights==num_kernels)
00506 
00507     if ((kernel->get_kernel_type()==K_COMBINED) &&
00508              (!((CCombinedKernel*)kernel)->get_append_subkernel_weights()))// for combined kernel
00509     {
00510         CCombinedKernel* k = (CCombinedKernel*) kernel;
00511 
00512         int32_t n = 0, i, j ;
00513 
00514         for (index_t k_idx=0; k_idx<k->get_num_kernels(); k_idx++)
00515         {
00516             CKernel* kn = k->get_kernel(k_idx);
00517             for(i=0;i<num;i++)
00518             {
00519                 if(a[i] != a_old[i])
00520                 {
00521                     kn->get_kernel_row(i,NULL,aicache, true);
00522                     for(j=0;j<num;j++)
00523                         W[j*num_kernels+n]+=(a[i]-a_old[i])*aicache[regression_fix_index(j)]*(float64_t)label[i];
00524                 }
00525             }
00526             SG_UNREF(kn);
00527             n++ ;
00528         }
00529     }
00530     else // hope the kernel is fast ...
00531     {
00532         float64_t* w_backup = SG_MALLOC(float64_t, num_kernels);
00533         float64_t* w1 = SG_MALLOC(float64_t, num_kernels);
00534 
00535         // backup and set to zero
00536         for (int32_t i=0; i<num_kernels; i++)
00537         {
00538             w_backup[i] = old_beta[i] ;
00539             w1[i]=0.0 ;
00540         }
00541         for (int32_t n=0; n<num_kernels; n++)
00542         {
00543             w1[n]=1.0 ;
00544             kernel->set_subkernel_weights(SGVector<float64_t>(w1, num_weights)) ;
00545 
00546             for(int32_t i=0;i<num;i++)
00547             {
00548                 if(a[i] != a_old[i])
00549                 {
00550                     for(int32_t j=0;j<num;j++)
00551                         W[j*num_kernels+n]+=(a[i]-a_old[i])*compute_kernel(i,j)*(float64_t)label[i];
00552                 }
00553             }
00554             w1[n]=0.0 ;
00555         }
00556 
00557         // restore old weights
00558         kernel->set_subkernel_weights(SGVector<float64_t>(w_backup,num_weights));
00559 
00560         SG_FREE(w_backup);
00561         SG_FREE(w1);
00562     }
00563 
00564     call_mkl_callback(a, label, lin, c, totdoc);
00565 }
00566 
00567 
00568 void CSVRLight::update_linear_component_mkl_linadd(
00569     int32_t* docs, int32_t* label, int32_t *active2dnum, float64_t *a,
00570     float64_t *a_old, int32_t *working2dnum, int32_t totdoc, float64_t *lin,
00571     float64_t *aicache, float64_t* c)
00572 {
00573     // kernel with LP_LINADD property is assumed to have
00574     // compute_by_subkernel functions
00575     int32_t num_weights = -1;
00576     int32_t num_kernels = kernel->get_num_subkernels() ;
00577     const float64_t* old_beta   = kernel->get_subkernel_weights(num_weights);
00578 
00579     ASSERT(num_weights==num_kernels)
00580 
00581     float64_t* w_backup=SG_MALLOC(float64_t, num_kernels);
00582     float64_t* w1=SG_MALLOC(float64_t, num_kernels);
00583 
00584     // backup and set to one
00585     for (int32_t i=0; i<num_kernels; i++)
00586     {
00587         w_backup[i] = old_beta[i] ;
00588         w1[i]=1.0 ;
00589     }
00590     // set the kernel weights
00591     kernel->set_subkernel_weights(SGVector<float64_t>(w1, num_weights));
00592 
00593     // create normal update (with changed alphas only)
00594     kernel->clear_normal();
00595     for(int32_t ii=0, i=0;(i=working2dnum[ii])>=0;ii++) {
00596         if(a[i] != a_old[i]) {
00597             kernel->add_to_normal(regression_fix_index(docs[i]), (a[i]-a_old[i])*(float64_t)label[i]);
00598         }
00599     }
00600 
00601     // determine contributions of different kernels
00602     for (int32_t i=0; i<num_vectors; i++)
00603         kernel->compute_by_subkernel(i,&W[i*num_kernels]) ;
00604 
00605     // restore old weights
00606     kernel->set_subkernel_weights(SGVector<float64_t>(w_backup,num_weights));
00607 
00608     call_mkl_callback(a, label, lin, c, totdoc);
00609 }
00610 
00611 void CSVRLight::call_mkl_callback(float64_t* a, int32_t* label, float64_t* lin, float64_t* c, int32_t totdoc)
00612 {
00613     int32_t num = totdoc;
00614     int32_t num_kernels = kernel->get_num_subkernels() ;
00615     float64_t sumalpha = 0;
00616     float64_t* sumw=SG_MALLOC(float64_t, num_kernels);
00617 
00618     for (int32_t i=0; i<num; i++)
00619         sumalpha-=a[i]*(learn_parm->eps[i]-label[i]*c[i]);
00620 
00621 #ifdef HAVE_LAPACK
00622     int nk = (int) num_kernels; // calling external lib
00623     double* alphay  = SG_MALLOC(double, num);
00624     for (int32_t i=0; i<num; i++)
00625         alphay[i]=a[i]*label[i];
00626 
00627     for (int32_t i=0; i<num_kernels; i++)
00628         sumw[i]=0;
00629 
00630     cblas_dgemv(CblasColMajor, CblasNoTrans, nk, (int) num, 0.5, (double*) W,
00631         nk, (double*) alphay, 1, 1.0, (double*) sumw, 1);
00632 
00633     SG_FREE(alphay);
00634 #else
00635     for (int32_t d=0; d<num_kernels; d++)
00636     {
00637         sumw[d]=0;
00638         for(int32_t i=0; i<num; i++)
00639             sumw[d] += 0.5*a[i]*label[i]*W[i*num_kernels+d];
00640     }
00641 #endif
00642 
00643     if (callback)
00644         mkl_converged=callback(mkl, sumw, sumalpha);
00645 
00646     const float64_t* new_beta   = kernel->get_subkernel_weights(num_kernels);
00647 
00648     // update lin
00649 #ifdef HAVE_LAPACK
00650     cblas_dgemv(CblasColMajor, CblasTrans, nk, (int) num, 1.0, (double*) W,
00651         nk, (double*) new_beta, 1, 0.0, (double*) lin, 1);
00652 #else
00653     for(int32_t i=0; i<num; i++)
00654         lin[i]=0 ;
00655     for (int32_t d=0; d<num_kernels; d++)
00656         if (new_beta[d]!=0)
00657             for(int32_t i=0; i<num; i++)
00658                 lin[i] += new_beta[d]*W[i*num_kernels+d] ;
00659 #endif
00660 
00661 
00662     SG_FREE(sumw);
00663 }
00664 
00665 
00666 void CSVRLight::reactivate_inactive_examples(
00667     int32_t* label, float64_t *a, SHRINK_STATE *shrink_state, float64_t *lin,
00668     float64_t *c, int32_t totdoc, int32_t iteration, int32_t *inconsistent,
00669     int32_t* docs, float64_t *aicache, float64_t *maxdiff)
00670      /* Make all variables active again which had been removed by
00671         shrinking. */
00672      /* Computes lin for those variables from scratch. */
00673 {
00674   register int32_t i=0,j,ii=0,jj,t,*changed2dnum,*inactive2dnum;
00675   int32_t *changed,*inactive;
00676   register float64_t *a_old,dist;
00677   float64_t ex_c,target;
00678 
00679   if (kernel->has_property(KP_LINADD) && get_linadd_enabled()) { /* special linear case */
00680       a_old=shrink_state->last_a;
00681 
00682       kernel->clear_normal();
00683       int32_t num_modified=0;
00684       for(i=0;i<totdoc;i++) {
00685           if(a[i] != a_old[i]) {
00686               kernel->add_to_normal(regression_fix_index(docs[i]), ((a[i]-a_old[i])*(float64_t)label[i]));
00687               a_old[i]=a[i];
00688               num_modified++;
00689           }
00690       }
00691 
00692       if (num_modified>0)
00693       {
00694           for(i=0;i<totdoc;i++) {
00695               if(!shrink_state->active[i]) {
00696                   lin[i]=shrink_state->last_lin[i]+kernel->compute_optimized(regression_fix_index(docs[i]));
00697               }
00698               shrink_state->last_lin[i]=lin[i];
00699           }
00700       }
00701   }
00702   else
00703   {
00704       changed=SG_MALLOC(int32_t, totdoc);
00705       changed2dnum=SG_MALLOC(int32_t, totdoc+11);
00706       inactive=SG_MALLOC(int32_t, totdoc);
00707       inactive2dnum=SG_MALLOC(int32_t, totdoc+11);
00708       for(t=shrink_state->deactnum-1;(t>=0) && shrink_state->a_history[t];t--) {
00709           if(verbosity>=2) {
00710               SG_INFO("%ld..",t)
00711           }
00712           a_old=shrink_state->a_history[t];
00713           for(i=0;i<totdoc;i++) {
00714               inactive[i]=((!shrink_state->active[i])
00715                       && (shrink_state->inactive_since[i] == t));
00716               changed[i]= (a[i] != a_old[i]);
00717           }
00718           compute_index(inactive,totdoc,inactive2dnum);
00719           compute_index(changed,totdoc,changed2dnum);
00720 
00721           for(ii=0;(i=changed2dnum[ii])>=0;ii++) {
00722               CKernelMachine::kernel->get_kernel_row(i,inactive2dnum,aicache);
00723               for(jj=0;(j=inactive2dnum[jj])>=0;jj++)
00724                   lin[j]+=(a[i]-a_old[i])*aicache[j]*(float64_t)label[i];
00725           }
00726       }
00727       SG_FREE(changed);
00728       SG_FREE(changed2dnum);
00729       SG_FREE(inactive);
00730       SG_FREE(inactive2dnum);
00731   }
00732 
00733   (*maxdiff)=0;
00734   for(i=0;i<totdoc;i++) {
00735     shrink_state->inactive_since[i]=shrink_state->deactnum-1;
00736     if(!inconsistent[i]) {
00737       dist=(lin[i]-model->b)*(float64_t)label[i];
00738       target=-(learn_parm->eps[i]-(float64_t)label[i]*c[i]);
00739       ex_c=learn_parm->svm_cost[i]-learn_parm->epsilon_a;
00740       if((a[i]>learn_parm->epsilon_a) && (dist > target)) {
00741     if((dist-target)>(*maxdiff))  /* largest violation */
00742       (*maxdiff)=dist-target;
00743       }
00744       else if((a[i]<ex_c) && (dist < target)) {
00745     if((target-dist)>(*maxdiff))  /* largest violation */
00746       (*maxdiff)=target-dist;
00747       }
00748       if((a[i]>(0+learn_parm->epsilon_a))
00749      && (a[i]<ex_c)) {
00750     shrink_state->active[i]=1;                         /* not at bound */
00751       }
00752       else if((a[i]<=(0+learn_parm->epsilon_a)) && (dist < (target+learn_parm->epsilon_shrink))) {
00753     shrink_state->active[i]=1;
00754       }
00755       else if((a[i]>=ex_c)
00756           && (dist > (target-learn_parm->epsilon_shrink))) {
00757     shrink_state->active[i]=1;
00758       }
00759       else if(learn_parm->sharedslack) { /* make all active when sharedslack */
00760     shrink_state->active[i]=1;
00761       }
00762     }
00763   }
00764   if (use_kernel_cache) { /* update history for non-linear */
00765       for(i=0;i<totdoc;i++) {
00766           (shrink_state->a_history[shrink_state->deactnum-1])[i]=a[i];
00767       }
00768       for(t=shrink_state->deactnum-2;(t>=0) && shrink_state->a_history[t];t--) {
00769           SG_FREE(shrink_state->a_history[t]);
00770           shrink_state->a_history[t]=0;
00771       }
00772   }
00773 }
00774 #endif //USE_SVMLIGHT
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation