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

SHOGUN Machine Learning Toolbox - Documentation