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