00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "classifier/svm/gmnpmkl.h"
00012
00013 lpwrapper::lpwrapper()
00014 {
00015 lpwrappertype=-1;
00016 }
00017
00018 lpwrapper::~lpwrapper()
00019 {
00020
00021 }
00022
00023 void lpwrapper::setup(const int32_t numkernels)
00024 {
00025 throw ShogunException("void lpwrapper::setup(...): not implemented in derived class");
00026 }
00027
00028 void lpwrapper::addconstraint(const ::std::vector<float64_t> & normw2,
00029 const float64_t sumofpositivealphas)
00030 {
00031 throw ShogunException("void lpwrapper::addconstraint(...): not implemented in derived class");
00032
00033 }
00034
00035 void lpwrapper::computeweights(std::vector<float64_t> & weights2)
00036 {
00037 throw ShogunException("void lpwrapper::computeweights(...): not implemented in derived class");
00038
00039 }
00040
00041
00042
00043 glpkwrapper::glpkwrapper()
00044 {
00045 lpwrappertype=0;
00046
00047 #if defined(USE_GLPK)
00048 linearproblem=NULL;
00049 #endif
00050 }
00051 glpkwrapper::~glpkwrapper()
00052 {
00053 #if defined(USE_GLPK)
00054 if (linearproblem!=NULL)
00055 {
00056 glp_delete_prob(linearproblem);
00057 linearproblem=NULL;
00058 }
00059 printf("deleting glpk linprob struct\n");
00060
00061 #endif
00062 }
00063
00064 glpkwrapper glpkwrapper::operator=(glpkwrapper & gl)
00065 {
00066 throw ShogunException(" glpkwrapper glpkwrapper::operator=(...): must not be called, glpk structure is currently not copiable");
00067
00068 }
00069 glpkwrapper::glpkwrapper(glpkwrapper & gl)
00070 {
00071 throw ShogunException(" glpkwrapper::glpkwrapper(glpkwrapper & gl): must not be called, glpk structure is currently not copiable");
00072
00073 }
00074
00075 glpkwrapper4CGMNPMKL::glpkwrapper4CGMNPMKL()
00076 {
00077 numkernels=0;
00078 }
00079
00080 glpkwrapper4CGMNPMKL::~glpkwrapper4CGMNPMKL()
00081 {
00082
00083 }
00084
00085 void glpkwrapper4CGMNPMKL::setup(const int32_t numkernels2)
00086 {
00087 #if defined(USE_GLPK)
00088 numkernels=numkernels2;
00089 if(numkernels<=1)
00090 {
00091
00092
00093 char bla[1000];
00094 sprintf(bla,"void glpkwrapper::setup(const int32_tnumkernels): input numkernels out of bounds: %d\n",numkernels);
00095 throw ShogunException(bla);
00096 }
00097
00098 if(NULL==linearproblem)
00099 {
00100 linearproblem=glp_create_prob();
00101 }
00102
00103 glp_set_obj_dir(linearproblem, GLP_MAX);
00104
00105 glp_add_cols(linearproblem,1+numkernels);
00106
00107
00108 glp_set_col_name(linearproblem,1,"theta");
00109 glp_set_col_bnds(linearproblem,1,GLP_FR,0.0,0.0);
00110 glp_set_obj_coef(linearproblem,1,1.0);
00111
00112
00113 int32_t offset=2;
00114 for(int32_t i=0; i<numkernels;++i)
00115 {
00116
00117
00118 char bla[100];
00119 sprintf(bla,"beta_%d",i);
00120 glp_set_col_name(linearproblem,offset+i,bla);
00121 glp_set_col_bnds(linearproblem,offset+i,GLP_DB,0.0,1.0);
00122 glp_set_obj_coef(linearproblem,offset+i,0.0);
00123
00124 }
00125
00126
00127
00128
00129 glp_add_rows(linearproblem,1);
00130 glp_set_row_name(linearproblem,1,"betas_sumupto_one");
00131
00132 int32_t*betainds(NULL);
00133 betainds=new int[1+numkernels];
00134 for(int32_t i=0; i<numkernels;++i)
00135 {
00136 betainds[1+i]=2+i;
00137 }
00138
00139 float64_t *betacoeffs(NULL);
00140 betacoeffs=new float64_t[1+numkernels];
00141
00142 for(int32_t i=0; i<numkernels;++i)
00143 {
00144 betacoeffs[1+i]=1;
00145 }
00146
00147 glp_set_mat_row(linearproblem,1,numkernels, betainds,betacoeffs);
00148 glp_set_row_bnds(linearproblem,1,GLP_FX,1.0,1.0);
00149
00150 delete[] betainds;
00151 betainds=NULL;
00152
00153 delete[] betacoeffs;
00154 betacoeffs=NULL;
00155 #else
00156 SG_ERROR("glpk.h from GNU glpk not included at compile time necessary in gmnpmkl.cpp/.h");
00157 #endif
00158
00159 }
00160
00161 void glpkwrapper4CGMNPMKL::addconstraint(const ::std::vector<float64_t> & normw2,
00162 const float64_t sumofpositivealphas)
00163 {
00164 #if defined(USE_GLPK)
00165
00166 ASSERT((int)normw2.size()==numkernels);
00167 ASSERT(sumofpositivealphas>=0);
00168
00169 glp_add_rows(linearproblem,1);
00170
00171 int32_t curconstraint=glp_get_num_rows(linearproblem);
00172
00173
00174
00175 char bla[100];
00176 sprintf(bla,"constraintno_%d",curconstraint);
00177 glp_set_row_name(linearproblem,curconstraint,bla);
00178
00179 int32_t *betainds(NULL);
00180 betainds=new int[1+1+numkernels];
00181
00182 betainds[1]=1;
00183 for(int32_t i=0; i<numkernels;++i)
00184 {
00185 betainds[2+i]=2+i;
00186 }
00187
00188 float64_t *betacoeffs(NULL);
00189 betacoeffs=new float64_t[1+1+numkernels];
00190
00191 betacoeffs[1]=-1;
00192
00193
00194 for(int32_t i=0; i<numkernels;++i)
00195 {
00196 betacoeffs[2+i]=0.5*normw2[i];
00197 }
00198 glp_set_mat_row(linearproblem,curconstraint,1+numkernels, betainds,betacoeffs);
00199 glp_set_row_bnds(linearproblem,curconstraint,GLP_LO,sumofpositivealphas,sumofpositivealphas);
00200
00201 delete[] betainds;
00202 betainds=NULL;
00203
00204 delete[] betacoeffs;
00205 betacoeffs=NULL;
00206
00207
00208 #else
00209 SG_ERROR("glpk.h from GNU glpk not included at compile time necessary in gmnpmkl.cpp/.h");
00210 #endif
00211 }
00212
00213 void glpkwrapper4CGMNPMKL::computeweights(std::vector<float64_t> & weights2)
00214 {
00215 #if defined(USE_GLPK)
00216 weights2.resize(numkernels);
00217
00218 glp_simplex(linearproblem,NULL);
00219
00220 float64_t sum=0;
00221 for(int32_t i=0; i< numkernels;++i)
00222 {
00223 weights2[i]=glp_get_col_prim(linearproblem, i+2);
00224 weights2[i]= ::std::max(0.0, ::std::min(1.0,weights2[i]));
00225 sum+= weights2[i];
00226
00227 }
00228
00229 if(sum>0)
00230 {
00231 for(int32_t i=0; i< numkernels;++i)
00232 {
00233 weights2[i]/=sum;
00234 }
00235 }
00236 else
00237 {
00238
00239
00240
00241 char bla[1000];
00242 sprintf(bla,"void glpkwrapper::computeweights(std::vector<float64_t> & weights2): sum of weights nonpositive %f\n",sum);
00243 throw ShogunException(bla);
00244 }
00245
00246
00247
00248 #else
00249 SG_ERROR("glpk.h from GNU glpk not included at compile time necessary in gmnpmkl.cpp/.h");
00250 #endif
00251 }
00252
00253
00254 CGMNPMKL::CGMNPMKL()
00255 : CMultiClassSVM(ONE_VS_REST)
00256 {
00257 svm=NULL;
00258 lpw=NULL;
00259
00260 lpwrappertype=0;
00261 thresh=0.01;
00262 maxiters=999;
00263
00264 numdat=0;
00265 numcl=0;
00266 numker=0;
00267 }
00268
00269 CGMNPMKL::CGMNPMKL(float64_t C, CKernel* k, CLabels* lab)
00270 : CMultiClassSVM(ONE_VS_REST, C, k, lab)
00271 {
00272 svm=NULL;
00273 lpw=NULL;
00274
00275 lpwrappertype=0;
00276 thresh=0.01;
00277 maxiters=999;
00278
00279 }
00280
00281
00282 CGMNPMKL::~CGMNPMKL()
00283 {
00284 SG_UNREF(svm);
00285 svm=NULL;
00286 delete lpw;
00287 lpw=NULL;
00288 }
00289
00290 void CGMNPMKL::lpsetup(const int32_t numkernels)
00291 {
00292 numker=numkernels;
00293 ASSERT(numker>0);
00294 switch(lpwrappertype)
00295 {
00296 case 0:
00297 if(lpw!=NULL)
00298 {
00299 delete lpw;
00300 }
00301 lpw=new glpkwrapper4CGMNPMKL;
00302 lpw->setup(numkernels);
00303 break;
00304
00305 default:
00306 {
00307
00308
00309
00310
00311 char bla[1000];
00312 sprintf(bla,"CGMNPMKL::setup(const int32_tnumkernels): unknown value for lpwrappertype %d\n ",lpwrappertype);
00313 throw ShogunException(bla);
00314 }
00315 break;
00316 }
00317 }
00318
00319 void CGMNPMKL::initsvm()
00320 {
00321 ASSERT(labels);
00322
00323 SG_UNREF(svm);
00324 svm=new CGMNPSVM;
00325 SG_REF(svm);
00326
00327 svm->set_C(get_C1(),get_C2());
00328 svm->set_epsilon(epsilon);
00329
00330
00331
00332 int32_t numlabels;
00333 float64_t * lb=labels->get_labels ( numlabels);
00334 ASSERT(numlabels>0);
00335 numdat=numlabels;
00336 numcl=labels->get_num_classes();
00337
00338 CLabels* newlab=new CLabels(lb, labels->get_num_labels() );
00339 delete[] lb;
00340 lb=NULL;
00341
00342 svm->set_labels(newlab);
00343
00344 newlab=NULL;
00345 }
00346
00347 void CGMNPMKL::init()
00348 {
00349
00350 if(NULL==kernel)
00351 {
00352
00353 throw ShogunException("CGMNPMKL::init(): the set kernel is NULL ");
00354 }
00355 else
00356 {
00357 if(kernel->get_kernel_type()!=K_COMBINED)
00358 {
00359
00360
00361
00362
00363 char bla[1000];
00364 sprintf(bla,"CGMNPMKL::init(): given kernel is not of type K_COMBINED %d \n",kernel->get_kernel_type());
00365 throw ShogunException(bla);
00366 }
00367
00368 numker=dynamic_cast<CCombinedKernel *>(kernel)->get_num_subkernels();
00369
00370 lpsetup(numker);
00371
00372
00373 }
00374 }
00375
00376
00377 bool CGMNPMKL::evaluatefinishcriterion(const int32_t numberofsilpiterations)
00378 {
00379 if((maxiters>0)&&(numberofsilpiterations>=maxiters))
00380 {
00381 return(true);
00382 }
00383
00384 if(weightshistory.size()>1)
00385 {
00386 std::vector<float64_t> wold,wnew;
00387
00388 wold=weightshistory[ weightshistory.size()-2 ];
00389 wnew=weightshistory.back();
00390 float64_t delta=0;
00391
00392 ASSERT(wold.size()==wnew.size());
00393
00394 for(size_t i=0;i< wnew.size();++i)
00395 {
00396 delta+=(wold[i]-wnew[i])*(wold[i]-wnew[i]);
00397 }
00398 delta=sqrt(delta);
00399
00400 SG_SPRINT( "CGMNPMKL::evaluatefinishcriterion(): L2 norm of changes= %f, required for termination by member variables thresh= %f \n",delta, thresh);
00401
00402 if( (delta < thresh)&&(numberofsilpiterations>=1) )
00403 {
00404 return(true);
00405 }
00406 }
00407
00408 return(false);
00409 }
00410
00411 void CGMNPMKL::addingweightsstep( const std::vector<float64_t> & curweights)
00412 {
00413
00414
00415 if(weightshistory.size()>2)
00416 {
00417 weightshistory.erase(weightshistory.begin());
00418 }
00419
00420 float64_t* weights(NULL);
00421 weights=new float64_t[curweights.size()];
00422 std::copy(curweights.begin(),curweights.end(),weights);
00423
00424 kernel->set_subkernel_weights( weights, curweights.size());
00425 delete[] weights;
00426 weights=NULL;
00427
00428 initsvm();
00429
00430
00431 if (numdat!=kernel->get_num_vec_lhs())
00432 {
00433 SG_ERROR("numdat (%ld) does not match kernel->get_num_vec_lhs() (%ld)\n",
00434 numdat, kernel->get_num_vec_lhs());
00435 }
00436
00437 if (numdat!=kernel->get_num_vec_rhs())
00438 {
00439 SG_ERROR("numdat (%ld) does not match kernel->get_num_vec_rhs() (%ld)\n",
00440 numdat, kernel->get_num_vec_rhs());
00441 }
00442
00443 svm->set_kernel(kernel);
00444 svm->train();
00445
00446 float64_t sumofsignfreealphas=getsumofsignfreealphas();
00447 int32_t numkernels=dynamic_cast<CCombinedKernel *>(kernel)->get_num_subkernels();
00448
00449
00450 std::vector<float64_t> normw2(numkernels);
00451 for(int32_t ind=0; ind < numkernels; ++ind )
00452 {
00453 normw2[ind]=getsquarenormofprimalcoefficients( ind );
00454 }
00455
00456
00457 lpw->addconstraint(normw2,sumofsignfreealphas);
00458 }
00459
00460 float64_t CGMNPMKL::getsumofsignfreealphas()
00461 {
00462
00463
00464 std::vector<int> trainlabels2(labels->get_num_labels());
00465 int32_t tmpint;
00466 int32_t * lab=labels->get_int_labels ( tmpint);
00467 std::copy(lab,lab+labels->get_num_labels(), trainlabels2.begin());
00468 delete[] lab;
00469 lab=NULL;
00470
00471
00472 ASSERT(trainlabels2.size()>0);
00473 float64_t sum=0;
00474
00475 for(int32_t nc=0; nc< labels->get_num_classes();++nc)
00476 {
00477 CSVM * sm=svm->get_svm(nc);
00478
00479 float64_t bia=sm->get_bias();
00480 sum+= bia*bia;
00481
00482 SG_UNREF(sm);
00483 }
00484
00485 ::std::vector< ::std::vector<float64_t> > basealphas;
00486 svm->getbasealphas( basealphas);
00487
00488 for(size_t lb=0; lb< trainlabels2.size();++lb)
00489 {
00490 for(int32_t nc=0; nc< labels->get_num_classes();++nc)
00491 {
00492
00493 CSVM * sm=svm->get_svm(nc);
00494
00495
00496 if((int)nc!=trainlabels2[lb])
00497 {
00498 CSVM * sm2=svm->get_svm(trainlabels2[lb]);
00499
00500 float64_t bia1=sm2->get_bias();
00501 float64_t bia2=sm->get_bias();
00502 SG_UNREF(sm2);
00503
00504 sum+= -basealphas[nc][lb]*(bia1-bia2-1);
00505 }
00506
00507 SG_UNREF(sm);
00508
00509
00510 }
00511 }
00512
00513 return(sum);
00514 }
00515
00516 float64_t CGMNPMKL::getsquarenormofprimalcoefficients(
00517 const int32_t ind)
00518 {
00519
00520
00521 CKernel * ker=dynamic_cast<CCombinedKernel *>(kernel)->get_kernel(ind);
00522
00523 float64_t tmp=0;
00524
00525 for(int32_t classindex=0; classindex< labels->get_num_classes();++classindex)
00526 {
00527 CSVM * sm=svm->get_svm(classindex);
00528
00529 for (int32_t i=0; i < sm->get_num_support_vectors(); ++i)
00530 {
00531 float64_t alphai=sm->get_alpha(i);
00532 int32_t svindi= sm->get_support_vector(i);
00533
00534 for (int32_t k=0; k < sm->get_num_support_vectors(); ++k)
00535 {
00536 float64_t alphak=sm->get_alpha(k);
00537 int32_t svindk=sm->get_support_vector(k);
00538
00539 tmp+=alphai*ker->kernel(svindi,svindk)
00540 *alphak;
00541
00542 }
00543 }
00544 SG_UNREF(sm);
00545 }
00546 SG_UNREF(ker);
00547 ker=NULL;
00548
00549 return(tmp);
00550 }
00551
00552
00553 bool CGMNPMKL::train()
00554 {
00555 init();
00556 weightshistory.clear();
00557
00558 int32_t numkernels=dynamic_cast<CCombinedKernel *>(kernel)->get_num_subkernels();
00559
00560 ::std::vector<float64_t> curweights(numkernels,1.0/numkernels);
00561 weightshistory.push_back(curweights);
00562
00563 SG_PRINT("initial weights in silp \n");
00564 for(size_t i=0; i< curweights.size();++i)
00565 {
00566 SG_PRINT("%f ",curweights[i]);
00567 }
00568 SG_PRINT("\n");
00569
00570
00571 addingweightsstep(curweights);
00572
00573
00574 int32_t numberofsilpiterations=0;
00575 bool final=false;
00576 while(false==final)
00577 {
00578
00579 curweights.clear();
00580 lpw->computeweights(curweights);
00581 weightshistory.push_back(curweights);
00582
00583 SG_SPRINT("SILP iteration %d weights in silp \n",numberofsilpiterations);
00584 for(size_t i=0; i< curweights.size();++i)
00585 {
00586 SG_SPRINT("%f ",curweights[i]);
00587 }
00588 SG_SPRINT("\n");
00589
00590 final=evaluatefinishcriterion(numberofsilpiterations);
00591 ++numberofsilpiterations;
00592
00593 addingweightsstep(curweights);
00594
00595 }
00596
00597
00598
00599 ASSERT(numcl>=1);
00600 create_multiclass_svm(numcl);
00601
00602 for (int32_t i=0; i<numcl; i++)
00603 {
00604 CSVM* osvm=svm->get_svm(i);
00605 CSVM* nsvm=new CSVM(osvm->get_num_support_vectors());
00606
00607 for (int32_t k=0; k<osvm->get_num_support_vectors() ; k++)
00608 {
00609 nsvm->set_alpha(k, osvm->get_alpha(k) );
00610 nsvm->set_support_vector(k,osvm->get_support_vector(k) );
00611 }
00612 nsvm->set_bias(osvm->get_bias() );
00613 set_svm(i, nsvm);
00614
00615 SG_UNREF(osvm);
00616 osvm=NULL;
00617 }
00618
00619 SG_UNREF(svm);
00620 svm=NULL;
00621 return(true);
00622 }
00623
00624
00625
00626
00627 float64_t* CGMNPMKL::getsubkernelweights(int32_t & numweights)
00628 {
00629 if((numker<=0 )||weightshistory.empty() )
00630 {
00631 numweights=0;
00632 return(NULL);
00633 }
00634
00635 std::vector<float64_t> subkerw=weightshistory.back();
00636 ASSERT(numker=subkerw.size());
00637 numweights=numker;
00638
00639 float64_t* res=new float64_t[numker];
00640 std::copy(subkerw.begin(), subkerw.end(),res);
00641 return(res);
00642 }