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

SHOGUN Machine Learning Toolbox - Documentation