SHOGUN
v3.2.0
|
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*)¶ms[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