GNPPLib.cpp

Go to the documentation of this file.
00001 /*-----------------------------------------------------------------------
00002  *
00003  * This program is free software; you can redistribute it and/or modify
00004  * it under the terms of the GNU General Public License as published by
00005  * the Free Software Foundation; either version 3 of the License, or
00006  * (at your option) any later version.
00007  *
00008  * Library of solvers for Generalized Nearest Point Problem (GNPP).
00009  *
00010  * Written (W) 1999-2008 Vojtech Franc, xfrancv@cmp.felk.cvut.cz
00011  * Copyright (C) 1999-2008 Center for Machine Perception, CTU FEL Prague 
00012  *
00013 -------------------------------------------------------------------- */
00014 
00015 #include <math.h>
00016 #include <limits.h>
00017 #include <shogun/lib/common.h>
00018 #include <shogun/io/SGIO.h>
00019 #include <shogun/mathematics/Math.h>
00020 
00021 #include <shogun/classifier/svm/GNPPLib.h>
00022 #include <shogun/kernel/Kernel.h>
00023 
00024 using namespace shogun;
00025 
00026 #define HISTORY_BUF 1000000
00027 
00028 #define MINUS_INF INT_MIN
00029 #define PLUS_INF  INT_MAX
00030 
00031 #define INDEX(ROW,COL,DIM) ((COL*DIM)+ROW)
00032 
00033 CGNPPLib::CGNPPLib()
00034 {
00035     SG_UNSTABLE("CGNPPLib::CGNPPLib()", "\n");
00036 
00037     kernel_columns = NULL;
00038     cache_index = NULL;
00039     first_kernel_inx = 0;
00040     Cache_Size = 0;
00041     m_num_data = 0;
00042     m_reg_const = 0.0;
00043     m_vector_y = NULL;
00044     m_kernel = NULL;
00045 }
00046 
00047 CGNPPLib::CGNPPLib(
00048     float64_t* vector_y, CKernel* kernel, int32_t num_data, float64_t reg_const)
00049 : CSGObject()
00050 {
00051   m_reg_const = reg_const;
00052   m_num_data = num_data;
00053   m_vector_y = vector_y;
00054   m_kernel = kernel;
00055 
00056   Cache_Size = ((int64_t) kernel->get_cache_size())*1024*1024/(sizeof(float64_t)*num_data);
00057   Cache_Size = CMath::min(Cache_Size, (int64_t) num_data);
00058 
00059   SG_INFO("using %d kernel cache lines\n", Cache_Size);
00060   ASSERT(Cache_Size>=2);
00061 
00062   /* allocates memory for kernel cache */
00063   kernel_columns = SG_MALLOC(float64_t*, Cache_Size);
00064   cache_index = SG_MALLOC(float64_t, Cache_Size);
00065 
00066   for(int32_t i = 0; i < Cache_Size; i++ ) 
00067   {
00068     kernel_columns[i] = SG_MALLOC(float64_t, num_data);
00069     cache_index[i] = -2;
00070   }
00071   first_kernel_inx = 0;
00072 
00073 }
00074 
00075 CGNPPLib::~CGNPPLib()
00076 {
00077   for(int32_t i = 0; i < Cache_Size; i++ ) 
00078       SG_FREE(kernel_columns[i]);
00079 
00080   SG_FREE(cache_index);
00081   SG_FREE(kernel_columns);
00082 }
00083 
00084 /* --------------------------------------------------------------
00085  QP solver based on Mitchell-Demyanov-Malozemov  algorithm.
00086 
00087  Usage: exitflag = gnpp_mdm( diag_H, vector_c, vector_y,
00088        dim, tmax, tolabs, tolrel, th, &alpha, &t, &aHa11, &aHa22, &History );
00089 -------------------------------------------------------------- */
00090 int8_t CGNPPLib::gnpp_mdm(float64_t *diag_H,
00091                        float64_t *vector_c,
00092                        float64_t *vector_y,
00093                        int32_t dim,
00094                        int32_t tmax,
00095                        float64_t tolabs,
00096                        float64_t tolrel,
00097                        float64_t th,
00098                        float64_t *alpha,
00099                        int32_t  *ptr_t,
00100                        float64_t *ptr_aHa11,
00101                        float64_t *ptr_aHa22,
00102                        float64_t **ptr_History,
00103                        int32_t verb)
00104 {
00105   float64_t LB;
00106   float64_t UB;
00107   float64_t aHa11, aHa12, aHa22, ac1, ac2;
00108   float64_t tmp;
00109   float64_t Huu, Huv, Hvv;
00110   float64_t min_beta1, max_beta1, min_beta2, max_beta2, beta;
00111   float64_t lambda;
00112   float64_t delta1, delta2;
00113   float64_t *History;
00114   float64_t *Ha1;
00115   float64_t *Ha2;
00116   float64_t *tmp_ptr;
00117   float64_t *col_u, *col_v;
00118   float64_t *col_v1, *col_v2;
00119   int64_t u1=0, u2=0;
00120   int64_t v1, v2;
00121   int64_t i;
00122   int64_t t;
00123   int64_t History_size;
00124   int8_t exitflag;
00125 
00126   /* ------------------------------------------------------------ */
00127   /* Initialization                                               */
00128   /* ------------------------------------------------------------ */
00129 
00130   Ha1 = SG_MALLOC(float64_t, dim);
00131   if( Ha1 == NULL ) SG_ERROR("Not enough memory.\n");
00132   Ha2 = SG_MALLOC(float64_t, dim);
00133   if( Ha2 == NULL ) SG_ERROR("Not enough memory.\n");
00134 
00135   History_size = (tmax < HISTORY_BUF ) ? tmax+1 : HISTORY_BUF;
00136   History = SG_MALLOC(float64_t, History_size*2);
00137   if( History == NULL ) SG_ERROR("Not enough memory.\n");
00138 
00139   /* inx1 = firts of find( y ==1 ), inx2 = firts of find( y ==2 ) */
00140   v1 = -1; v2 = -1; i = 0;
00141   while( (v1 == -1 || v2 == -1) && i < dim ) {
00142     if( v1 == -1 && vector_y[i] == 1 ) { v1 = i; }
00143     if( v2 == -1 && vector_y[i] == 2 ) { v2 = i; } 
00144     i++;
00145   }
00146 
00147   col_v1 = (float64_t*)get_col(v1,-1);
00148   col_v2 = (float64_t*)get_col(v2,v1);
00149   
00150   aHa12 = col_v1[v2];
00151   aHa11 = diag_H[v1];
00152   aHa22 = diag_H[v2];
00153   ac1 = vector_c[v1];
00154   ac2 = vector_c[v2];
00155 
00156   min_beta1 = PLUS_INF; min_beta2 = PLUS_INF;
00157   for( i = 0; i < dim; i++ ) 
00158   {
00159     alpha[i] = 0;
00160     Ha1[i] = col_v1[i];
00161     Ha2[i] = col_v2[i];
00162 
00163     beta = Ha1[i] + Ha2[i] + vector_c[i];
00164 
00165     if( vector_y[i] == 1 && min_beta1 > beta ) {
00166       u1 = i;
00167       min_beta1 = beta;
00168     }
00169 
00170     if( vector_y[i] == 2 && min_beta2 > beta ) {
00171       u2 = i;
00172       min_beta2 = beta;
00173     }
00174   }
00175 
00176   alpha[v1] = 1;
00177   alpha[v2] = 1;
00178 
00179   UB = 0.5*(aHa11 + 2*aHa12 + aHa22) + ac1 + ac2;
00180   LB = min_beta1 + min_beta2 - 0.5*(aHa11 + 2*aHa12 + aHa22);
00181 
00182   delta1 = Ha1[v1] + Ha2[v1] + vector_c[v1] - min_beta1;
00183   delta2 = Ha1[v2] + Ha2[v2] + vector_c[v2] - min_beta2;
00184 
00185   t = 0;
00186   History[INDEX(0,0,2)] = LB;
00187   History[INDEX(1,0,2)] = UB;
00188 
00189   if( verb ) {
00190     SG_PRINT("Init: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
00191       UB, LB, UB-LB,(UB-LB)/UB);
00192   }  
00193 
00194   /* Stopping conditions */
00195   if( UB-LB <= tolabs ) exitflag = 1;
00196   else if(UB-LB <= CMath::abs(UB)*tolrel ) exitflag = 2;
00197   else if(LB > th) exitflag = 3;
00198   else exitflag = -1;
00199 
00200   /* ------------------------------------------------------------ */
00201   /* Main optimization loop                                       */
00202   /* ------------------------------------------------------------ */
00203 
00204   while( exitflag == -1 ) 
00205   {
00206     t++;     
00207 
00208     if( delta1 > delta2 ) 
00209     {
00210       col_u = (float64_t*)get_col(u1,-1);
00211       col_v = (float64_t*)get_col(v1,u1);
00212 
00213       Huu = diag_H[u1];
00214       Hvv = diag_H[v1];
00215       Huv = col_u[v1];
00216 
00217       lambda = delta1/(alpha[v1]*(Huu - 2*Huv + Hvv ));
00218       lambda = CMath::min(1.0,lambda);
00219 
00220       tmp = lambda*alpha[v1];
00221 
00222       aHa11 = aHa11 + 2*tmp*(Ha1[u1]-Ha1[v1])+tmp*tmp*( Huu - 2*Huv + Hvv );
00223       aHa12 = aHa12 + tmp*(Ha2[u1]-Ha2[v1]);
00224       ac1 = ac1 + tmp*(vector_c[u1]-vector_c[v1]);
00225 
00226       alpha[u1] = alpha[u1] + tmp;
00227       alpha[v1] = alpha[v1] - tmp;
00228 
00229       min_beta1 = PLUS_INF; min_beta2 = PLUS_INF;
00230       max_beta1 = MINUS_INF; max_beta2 = MINUS_INF; 
00231       for( i = 0; i < dim; i ++ )
00232       {
00233          Ha1[i] = Ha1[i] + tmp*(col_u[i] - col_v[i]);
00234 
00235          beta = Ha1[i] + Ha2[i] + vector_c[i];
00236          if( vector_y[i] == 1 ) 
00237            {
00238              if( min_beta1 > beta ) { u1 = i; min_beta1 = beta; }
00239              if( max_beta1 < beta && alpha[i] > 0 ) { v1 = i; max_beta1 = beta; }
00240            }
00241          else
00242            {
00243              if( min_beta2 > beta ) { u2 = i; min_beta2 = beta; }
00244              if( max_beta2 < beta && alpha[i] > 0) { v2 = i; max_beta2 = beta; }
00245            }
00246       }
00247     }
00248     else
00249     {
00250       col_u = (float64_t*)get_col(u2,-1);
00251       col_v = (float64_t*)get_col(v2,u2);
00252 
00253       Huu = diag_H[u2];
00254       Hvv = diag_H[v2];
00255       Huv = col_u[v2];
00256   
00257       lambda = delta2/(alpha[v2]*( Huu - 2*Huv + Hvv ));
00258       lambda = CMath::min(1.0,lambda);
00259 
00260       tmp = lambda*alpha[v2];
00261       aHa22 = aHa22 + 2*tmp*( Ha2[u2]-Ha2[v2]) + tmp*tmp*( Huu - 2*Huv + Hvv);
00262       aHa12 = aHa12 + tmp*(Ha1[u2]-Ha1[v2]);
00263       ac2 = ac2 + tmp*( vector_c[u2]-vector_c[v2] );
00264 
00265       alpha[u2] = alpha[u2] + tmp;
00266       alpha[v2] = alpha[v2] - tmp;
00267 
00268       min_beta1 = PLUS_INF; min_beta2 = PLUS_INF;
00269       max_beta1 = MINUS_INF; max_beta2 = MINUS_INF; 
00270       for(i = 0; i < dim; i++ ) 
00271       {  
00272          Ha2[i] = Ha2[i] + tmp*( col_u[i] - col_v[i] );
00273 
00274          beta = Ha1[i] + Ha2[i] + vector_c[i];
00275 
00276          if( vector_y[i] == 1 ) 
00277          {
00278            if( min_beta1 > beta ) { u1 = i; min_beta1 = beta; }
00279            if( max_beta1 < beta && alpha[i] > 0 ) { v1 = i; max_beta1 = beta; }
00280          }
00281          else
00282          {
00283            if( min_beta2 > beta ) { u2 = i; min_beta2 = beta; }
00284            if( max_beta2 < beta && alpha[i] > 0) { v2 = i; max_beta2 = beta; }
00285          }
00286       }
00287     }
00288 
00289     UB = 0.5*(aHa11 + 2*aHa12 + aHa22) + ac1 + ac2;
00290     LB = min_beta1 + min_beta2 - 0.5*(aHa11 + 2*aHa12 + aHa22);
00291   
00292     delta1 = Ha1[v1] + Ha2[v1] + vector_c[v1] - min_beta1;
00293     delta2 = Ha1[v2] + Ha2[v2] + vector_c[v2] - min_beta2;
00294 
00295     /* Stopping conditions */
00296     if( UB-LB <= tolabs ) exitflag = 1; 
00297     else if( UB-LB <= CMath::abs(UB)*tolrel ) exitflag = 2;
00298     else if(LB > th) exitflag = 3;
00299     else if(t >= tmax) exitflag = 0; 
00300 
00301     if( verb && (t % verb) == 0) {
00302      SG_PRINT("%d: UB=%f,LB=%f,UB-LB=%f,(UB-LB)/|UB|=%f\n",
00303         t, UB, LB, UB-LB,(UB-LB)/UB); 
00304     }  
00305 
00306     /* Store selected values */
00307     if( t < History_size ) {
00308       History[INDEX(0,t,2)] = LB;
00309       History[INDEX(1,t,2)] = UB;
00310     }
00311     else {
00312       tmp_ptr = SG_MALLOC(float64_t, (History_size+HISTORY_BUF)*2);
00313       if( tmp_ptr == NULL ) SG_ERROR("Not enough memory.\n");
00314       for( i = 0; i < History_size; i++ ) {
00315         tmp_ptr[INDEX(0,i,2)] = History[INDEX(0,i,2)];
00316         tmp_ptr[INDEX(1,i,2)] = History[INDEX(1,i,2)];
00317       }
00318       tmp_ptr[INDEX(0,t,2)] = LB;
00319       tmp_ptr[INDEX(1,t,2)] = UB;
00320       
00321       History_size += HISTORY_BUF;
00322       SG_FREE(History);
00323       History = tmp_ptr;
00324     }
00325   }
00326 
00327   /* print info about last iteration*/
00328   if(verb && (t % verb) ) {
00329     SG_PRINT("Exit: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
00330       UB, LB, UB-LB,(UB-LB)/UB);
00331   }  
00332 
00333   /*------------------------------------------------------- */
00334   /* Set outputs                                            */
00335   /*------------------------------------------------------- */
00336   (*ptr_t) = t;
00337   (*ptr_aHa11) = aHa11;
00338   (*ptr_aHa22) = aHa22;
00339   (*ptr_History) = History;
00340 
00341   /* Free memory */
00342   SG_FREE(Ha1);
00343   SG_FREE(Ha2);
00344   
00345   return( exitflag ); 
00346 }
00347 
00348 
00349 /* --------------------------------------------------------------
00350  QP solver based on Improved MDM algorithm (u fixed v optimized)
00351 
00352  Usage: exitflag = gnpp_imdm( diag_H, vector_c, vector_y,
00353        dim, tmax, tolabs, tolrel, th, &alpha, &t, &aHa11, &aHa22, &History );
00354 -------------------------------------------------------------- */
00355 int8_t CGNPPLib::gnpp_imdm(float64_t *diag_H,
00356             float64_t *vector_c,
00357             float64_t *vector_y,
00358             int32_t dim,
00359             int32_t tmax,
00360             float64_t tolabs,
00361             float64_t tolrel,
00362             float64_t th,
00363             float64_t *alpha,
00364             int32_t  *ptr_t,
00365             float64_t *ptr_aHa11,
00366             float64_t *ptr_aHa22,
00367             float64_t **ptr_History,
00368             int32_t verb)
00369 {
00370   float64_t LB;
00371   float64_t UB;
00372   float64_t aHa11, aHa12, aHa22, ac1, ac2;
00373   float64_t tmp;
00374   float64_t Huu, Huv, Hvv;
00375   float64_t min_beta1, max_beta1, min_beta2, max_beta2, beta;
00376   float64_t lambda;
00377   float64_t delta1, delta2;
00378   float64_t improv, max_improv;
00379   float64_t *History;
00380   float64_t *Ha1;
00381   float64_t *Ha2;
00382   float64_t *tmp_ptr;
00383   float64_t *col_u, *col_v;
00384   float64_t *col_v1, *col_v2;
00385   int64_t u1=0, u2=0;
00386   int64_t v1, v2;
00387   int64_t i;
00388   int64_t t;
00389   int64_t History_size;
00390   int8_t exitflag;
00391   int8_t which_case;
00392 
00393   /* ------------------------------------------------------------ */
00394   /* Initialization                                               */
00395   /* ------------------------------------------------------------ */
00396 
00397   Ha1 = SG_MALLOC(float64_t, dim);
00398   if( Ha1 == NULL ) SG_ERROR("Not enough memory.\n");
00399   Ha2 = SG_MALLOC(float64_t, dim);
00400   if( Ha2 == NULL ) SG_ERROR("Not enough memory.\n");
00401 
00402   History_size = (tmax < HISTORY_BUF ) ? tmax+1 : HISTORY_BUF;
00403   History = SG_MALLOC(float64_t, History_size*2);
00404   if( History == NULL ) SG_ERROR("Not enough memory.\n");
00405 
00406   /* inx1 = firts of find( y ==1 ), inx2 = firts of find( y ==2 ) */
00407   v1 = -1; v2 = -1; i = 0;
00408   while( (v1 == -1 || v2 == -1) && i < dim ) {
00409     if( v1 == -1 && vector_y[i] == 1 ) { v1 = i; }
00410     if( v2 == -1 && vector_y[i] == 2 ) { v2 = i; } 
00411     i++;
00412   }
00413 
00414   col_v1 = (float64_t*)get_col(v1,-1);
00415   col_v2 = (float64_t*)get_col(v2,v1);
00416   
00417   aHa12 = col_v1[v2];
00418   aHa11 = diag_H[v1];
00419   aHa22 = diag_H[v2];
00420   ac1 = vector_c[v1];
00421   ac2 = vector_c[v2];
00422 
00423   min_beta1 = PLUS_INF; min_beta2 = PLUS_INF;
00424   for( i = 0; i < dim; i++ ) 
00425   {
00426     alpha[i] = 0;
00427     Ha1[i] = col_v1[i];
00428     Ha2[i] = col_v2[i];
00429 
00430     beta = Ha1[i] + Ha2[i] + vector_c[i];
00431 
00432     if( vector_y[i] == 1 && min_beta1 > beta ) {
00433       u1 = i;
00434       min_beta1 = beta;
00435     }
00436 
00437     if( vector_y[i] == 2 && min_beta2 > beta ) {
00438       u2 = i;
00439       min_beta2 = beta;
00440     }
00441   }
00442 
00443   alpha[v1] = 1;
00444   alpha[v2] = 1;
00445 
00446   UB = 0.5*(aHa11 + 2*aHa12 + aHa22) + ac1 + ac2;
00447   LB = min_beta1 + min_beta2 - 0.5*(aHa11 + 2*aHa12 + aHa22);
00448 
00449   delta1 = Ha1[v1] + Ha2[v1] + vector_c[v1] - min_beta1;
00450   delta2 = Ha1[v2] + Ha2[v2] + vector_c[v2] - min_beta2;
00451 
00452   t = 0;
00453   History[INDEX(0,0,2)] = LB;
00454   History[INDEX(1,0,2)] = UB;
00455 
00456   if( verb ) {
00457     SG_PRINT("Init: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
00458       UB, LB, UB-LB,(UB-LB)/UB);
00459   }  
00460 
00461   if( delta1 > delta2 ) 
00462   {
00463      which_case = 1;
00464      col_u = (float64_t*)get_col(u1,v1);
00465      col_v = col_v1;
00466   }
00467   else
00468   {
00469      which_case = 2;
00470      col_u = (float64_t*)get_col(u2,v2);
00471      col_v = col_v2;
00472   }
00473 
00474   /* Stopping conditions */
00475   if( UB-LB <= tolabs ) exitflag = 1;
00476   else if(UB-LB <= CMath::abs(UB)*tolrel ) exitflag = 2;
00477   else if(LB > th) exitflag = 3;
00478   else exitflag = -1;
00479 
00480   /* ------------------------------------------------------------ */
00481   /* Main optimization loop                                       */
00482   /* ------------------------------------------------------------ */
00483 
00484   while( exitflag == -1 ) 
00485   {
00486     t++;     
00487 
00488     if( which_case == 1 )
00489     {
00490       Huu = diag_H[u1];
00491       Hvv = diag_H[v1];
00492       Huv = col_u[v1];
00493 
00494       lambda = delta1/(alpha[v1]*(Huu - 2*Huv + Hvv ));
00495       lambda = CMath::min(1.0,lambda);
00496 
00497       tmp = lambda*alpha[v1];
00498 
00499       aHa11 = aHa11 + 2*tmp*(Ha1[u1]-Ha1[v1])+tmp*tmp*( Huu - 2*Huv + Hvv );
00500       aHa12 = aHa12 + tmp*(Ha2[u1]-Ha2[v1]);
00501       ac1 = ac1 + tmp*(vector_c[u1]-vector_c[v1]);
00502 
00503       alpha[u1] = alpha[u1] + tmp;
00504       alpha[v1] = alpha[v1] - tmp;
00505 
00506       min_beta1 = PLUS_INF; min_beta2 = PLUS_INF;
00507       max_beta1 = MINUS_INF; max_beta2 = MINUS_INF; 
00508       for( i = 0; i < dim; i ++ )
00509       {
00510          Ha1[i] = Ha1[i] + tmp*(col_u[i] - col_v[i]);
00511 
00512          beta = Ha1[i] + Ha2[i] + vector_c[i];
00513          if( vector_y[i] == 1 ) 
00514            {
00515              if( min_beta1 > beta ) { u1 = i; min_beta1 = beta; }
00516              if( max_beta1 < beta && alpha[i] > 0 ) { v1 = i; max_beta1 = beta; }
00517            }
00518          else
00519            {
00520              if( min_beta2 > beta ) { u2 = i; min_beta2 = beta; }
00521              if( max_beta2 < beta && alpha[i] > 0) { v2 = i; max_beta2 = beta; }
00522            }
00523       }
00524     }
00525     else
00526     {
00527       Huu = diag_H[u2];
00528       Hvv = diag_H[v2];
00529       Huv = col_u[v2];
00530   
00531       lambda = delta2/(alpha[v2]*( Huu - 2*Huv + Hvv ));
00532       lambda = CMath::min(1.0,lambda);
00533 
00534       tmp = lambda*alpha[v2];
00535       aHa22 = aHa22 + 2*tmp*( Ha2[u2]-Ha2[v2]) + tmp*tmp*( Huu - 2*Huv + Hvv);
00536       aHa12 = aHa12 + tmp*(Ha1[u2]-Ha1[v2]);
00537       ac2 = ac2 + tmp*( vector_c[u2]-vector_c[v2] );
00538 
00539       alpha[u2] = alpha[u2] + tmp;
00540       alpha[v2] = alpha[v2] - tmp;
00541 
00542       min_beta1 = PLUS_INF; min_beta2 = PLUS_INF;
00543       max_beta1 = MINUS_INF; max_beta2 = MINUS_INF; 
00544       for(i = 0; i < dim; i++ ) 
00545       {  
00546          Ha2[i] = Ha2[i] + tmp*( col_u[i] - col_v[i] );
00547 
00548          beta = Ha1[i] + Ha2[i] + vector_c[i];
00549 
00550          if( vector_y[i] == 1 ) 
00551          {
00552            if( min_beta1 > beta ) { u1 = i; min_beta1 = beta; }
00553            if( max_beta1 < beta && alpha[i] > 0 ) { v1 = i; max_beta1 = beta; }
00554          }
00555          else
00556          {
00557            if( min_beta2 > beta ) { u2 = i; min_beta2 = beta; }
00558            if( max_beta2 < beta && alpha[i] > 0) { v2 = i; max_beta2 = beta; }
00559          }
00560       }
00561     }
00562 
00563     UB = 0.5*(aHa11 + 2*aHa12 + aHa22) + ac1 + ac2;
00564     LB = min_beta1 + min_beta2 - 0.5*(aHa11 + 2*aHa12 + aHa22);
00565   
00566     delta1 = Ha1[v1] + Ha2[v1] + vector_c[v1] - min_beta1;
00567     delta2 = Ha1[v2] + Ha2[v2] + vector_c[v2] - min_beta2;
00568 
00569     if( delta1 > delta2 ) 
00570     {
00571        col_u = (float64_t*)get_col(u1,-1);
00572 
00573       /* search for optimal v while u is fixed */
00574       for( max_improv =  MINUS_INF, i = 0; i < dim; i++ ) {
00575 
00576         if( vector_y[i] == 1 && alpha[i] != 0 ) {
00577 
00578           beta = Ha1[i] + Ha2[i] + vector_c[i];
00579 
00580           if( beta >= min_beta1 ) {
00581 
00582             tmp = diag_H[u1] - 2*col_u[i] + diag_H[i];
00583             if( tmp != 0 ) {
00584               improv = (0.5*(beta-min_beta1)*(beta-min_beta1))/tmp;
00585 
00586               if( improv > max_improv ) {
00587                 max_improv = improv;
00588                 v1 = i;
00589               }
00590             }
00591           }
00592         }
00593       }
00594       col_v = (float64_t*)get_col(v1,u1);
00595       delta1 = Ha1[v1] + Ha2[v1] + vector_c[v1] - min_beta1;
00596       which_case = 1;
00597       
00598     }
00599     else
00600     {
00601        col_u = (float64_t*)get_col(u2,-1);
00602 
00603       /* search for optimal v while u is fixed */
00604       for( max_improv =  MINUS_INF, i = 0; i < dim; i++ ) {
00605 
00606         if( vector_y[i] == 2 && alpha[i] != 0 ) {
00607 
00608           beta = Ha1[i] + Ha2[i] + vector_c[i];
00609 
00610           if( beta >= min_beta2 ) {
00611 
00612             tmp = diag_H[u2] - 2*col_u[i] + diag_H[i];
00613             if( tmp != 0 ) {
00614               improv = (0.5*(beta-min_beta2)*(beta-min_beta2))/tmp;
00615 
00616               if( improv > max_improv ) {
00617                 max_improv = improv;
00618                 v2 = i;
00619               }
00620             }
00621           }
00622         }
00623       }
00624 
00625       col_v = (float64_t*)get_col(v2,u2);
00626       delta2 = Ha1[v2] + Ha2[v2] + vector_c[v2] - min_beta2;
00627       which_case = 2;
00628     }
00629     
00630 
00631     /* Stopping conditions */
00632     if( UB-LB <= tolabs ) exitflag = 1; 
00633     else if( UB-LB <= CMath::abs(UB)*tolrel ) exitflag = 2;
00634     else if(LB > th) exitflag = 3;
00635     else if(t >= tmax) exitflag = 0; 
00636 
00637     if( verb && (t % verb) == 0) {
00638      SG_PRINT("%d: UB=%f,LB=%f,UB-LB=%f,(UB-LB)/|UB|=%f\n",
00639         t, UB, LB, UB-LB,(UB-LB)/UB); 
00640     }  
00641 
00642     /* Store selected values */
00643     if( t < History_size ) {
00644       History[INDEX(0,t,2)] = LB;
00645       History[INDEX(1,t,2)] = UB;
00646     }
00647     else {
00648       tmp_ptr = SG_MALLOC(float64_t, (History_size+HISTORY_BUF)*2);
00649       if( tmp_ptr == NULL ) SG_ERROR("Not enough memory.\n");
00650       for( i = 0; i < History_size; i++ ) {
00651         tmp_ptr[INDEX(0,i,2)] = History[INDEX(0,i,2)];
00652         tmp_ptr[INDEX(1,i,2)] = History[INDEX(1,i,2)];
00653       }
00654       tmp_ptr[INDEX(0,t,2)] = LB;
00655       tmp_ptr[INDEX(1,t,2)] = UB;
00656       
00657       History_size += HISTORY_BUF;
00658       SG_FREE(History);
00659       History = tmp_ptr;
00660     }
00661   }
00662 
00663   /* print info about last iteration*/
00664   if(verb && (t % verb) ) {
00665     SG_PRINT("Exit: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
00666       UB, LB, UB-LB,(UB-LB)/UB);
00667   }  
00668 
00669   /*------------------------------------------------------- */
00670   /* Set outputs                                            */
00671   /*------------------------------------------------------- */
00672   (*ptr_t) = t;
00673   (*ptr_aHa11) = aHa11;
00674   (*ptr_aHa22) = aHa22;
00675   (*ptr_History) = History;
00676 
00677   /* Free memory */
00678   SG_FREE(Ha1);
00679   SG_FREE(Ha2);
00680   
00681   return( exitflag ); 
00682 }
00683 
00684 
00685 float64_t* CGNPPLib::get_col(int64_t a, int64_t b)
00686 {
00687   float64_t *col_ptr;
00688   float64_t y;
00689   int64_t i;
00690   int64_t inx;
00691 
00692   inx = -1;
00693   for( i=0; i < Cache_Size; i++ ) {
00694     if( cache_index[i] == a ) { inx = i; break; }
00695   }
00696     
00697   if( inx != -1 ) {
00698     col_ptr = kernel_columns[inx];
00699     return( col_ptr );
00700   }
00701    
00702   col_ptr = kernel_columns[first_kernel_inx];
00703   cache_index[first_kernel_inx] = a;
00704 
00705   first_kernel_inx++;
00706   if( first_kernel_inx >= Cache_Size ) first_kernel_inx = 0;
00707 
00708   y = m_vector_y[a];
00709   for( i=0; i < m_num_data; i++ ) {
00710     if( m_vector_y[i] == y )  
00711     {
00712       col_ptr[i] = 2*m_kernel->kernel(i,a); 
00713     }
00714     else 
00715     {
00716       col_ptr[i] = -2*m_kernel->kernel(i,a);
00717     }
00718   }
00719 
00720   col_ptr[a] = col_ptr[a] + m_reg_const;
00721 
00722   return( col_ptr );
00723 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation