00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
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
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
00086
00087
00088
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
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
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
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
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
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
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
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
00335
00336 (*ptr_t) = t;
00337 (*ptr_aHa11) = aHa11;
00338 (*ptr_aHa22) = aHa22;
00339 (*ptr_History) = History;
00340
00341
00342 SG_FREE(Ha1);
00343 SG_FREE(Ha2);
00344
00345 return( exitflag );
00346 }
00347
00348
00349
00350
00351
00352
00353
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
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
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
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
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
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
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
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
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
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
00671
00672 (*ptr_t) = t;
00673 (*ptr_aHa11) = aHa11;
00674 (*ptr_aHa22) = aHa22;
00675 (*ptr_History) = History;
00676
00677
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 }