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

SHOGUN Machine Learning Toolbox - Documentation