gmnpmkl.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) 2009 Alexander Binder
00008  * Copyright (C) 2009 Fraunhofer Institute FIRST and Max-Planck-Society
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 //TODO: check for correctness of
00085 void glpkwrapper4CGMNPMKL::setup(const int32_t numkernels2)
00086 {
00087 #if defined(USE_GLPK)
00088     numkernels=numkernels2;
00089     if(numkernels<=1)
00090     {
00091         //std::ostringstream helper;
00092         //helper << "void glpkwrapper::setup(const int32_tnumkernels): input numkernels out of bounds: "<<numkernels <<std::endl;
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); // one for theta (objectivelike), all others for betas=weights
00106 
00107     //set up theta
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     //set up betas
00113     int32_t offset=2;
00114     for(int32_t i=0; i<numkernels;++i)
00115     {
00116         //std::ostringstream helper;
00117         //helper << "beta_"<<i;
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     // objective is maximize theta over { beta and theta } subject to constraints
00127 
00128     //set sumupconstraint32_t/sum_l \beta_l=1
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; // coefficient for theta stays zero, therefore start at 2 not at 1 !
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     //std::ostringstream helper;
00174     //helper << "constraintno_"<<curconstraint;
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; // coefficient for theta stays zero, therefore start at 2 not at 1 !
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     //addedconstraint=true;
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); // standard parameters
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         //std::ostringstream helper;
00239         //helper << "void glpkwrapper::computeweights(std::vector<float64_t> & weights2): sum of weights nonpositive "<<sum <<std::endl;
00240         //throw ShogunException(helper.str().c_str());
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             //std::ostringstream helper;
00308             //helper << "CGMNPMKL::setup(const int32_tnumkernels): unknown value for lpwrappertype "<<lpwrappertype <<std::endl;
00309             //throw ShogunException(helper.str().c_str());
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     //CGNMPSVM->set_labels(get_labels());
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             //std::ostringstream helper;
00360             //helper << "CGMNPMKL::init(): given kernel is not of type K_COMBINED "<<k->get_kernel_type() <<std::endl;
00361             //throw ShogunException(helper.str().c_str());
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     //weightshistory.push_back(curweights);
00414     //error prone?
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     //number of labels equal to number of features?
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     //add a constraint
00457     lpw->addconstraint(normw2,sumofsignfreealphas);
00458 }
00459 
00460 float64_t CGMNPMKL::getsumofsignfreealphas()
00461 {
00462     //returns \sum_y b_y^2-\sum_i \sum_{ y \neq y_i} \alpha_{iy}(b_{y_i}-b_y-1)
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     // alphas are already correctly transformed!
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);// svmstuff[classindex].alphas[i];
00532             int32_t svindi= sm->get_support_vector(i); //svmstuff[classindex].svind[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     } // while(false==final)
00596 
00597 
00598     //set alphas, bias, support vecs
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 }

SHOGUN Machine Learning Toolbox - Documentation