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