00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
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
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
00071
00072
00073
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
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
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
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
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
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
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
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
00320
00321 (*ptr_t) = t;
00322 (*ptr_aHa11) = aHa11;
00323 (*ptr_aHa22) = aHa22;
00324 (*ptr_History) = History;
00325
00326
00327 delete[] Ha1 ;
00328 delete[] Ha2;
00329
00330 return( exitflag );
00331 }
00332
00333
00334
00335
00336
00337
00338
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
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
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
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
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
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
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
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
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
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
00656
00657 (*ptr_t) = t;
00658 (*ptr_aHa11) = aHa11;
00659 (*ptr_aHa22) = aHa22;
00660 (*ptr_History) = History;
00661
00662
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