SHOGUN
v3.2.0
|
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 }