00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063 #include "classifier/svm/gmnplib.h"
00064 #include "lib/Mathematics.h"
00065
00066 #include <string.h>
00067 #include <limits.h>
00068
00069 #define HISTORY_BUF 1000000
00070
00071 #define MINUS_INF INT_MIN
00072 #define PLUS_INF INT_MAX
00073
00074 #define INDEX(ROW,COL,DIM) ((COL*DIM)+ROW)
00075 #define KDELTA(A,B) (A==B)
00076 #define KDELTA4(A1,A2,A3,A4) ((A1==A2)||(A1==A3)||(A1==A4)||(A2==A3)||(A2==A4)||(A3==A4))
00077
00078 CGMNPLib::CGMNPLib(DREAL* vector_y, CKernel* kernel, INT num_data, INT num_virt_data, INT num_classes, DREAL reg_const)
00079 : CSGObject()
00080 {
00081 m_num_classes=num_classes;
00082 m_num_virt_data=num_virt_data;
00083 m_reg_const = reg_const;
00084 m_num_data = num_data;
00085 m_vector_y = vector_y;
00086 m_kernel = kernel;
00087
00088 Cache_Size = ((LONG) kernel->get_cache_size())*1024*1024/(sizeof(DREAL)*num_data);
00089 Cache_Size = CMath::min(Cache_Size, (LONG) num_data);
00090
00091 SG_INFO("using %d kernel cache lines\n", Cache_Size);
00092 ASSERT(Cache_Size>=2);
00093
00094
00095 kernel_columns = new DREAL*[Cache_Size];
00096 cache_index = new DREAL[Cache_Size];
00097
00098 for(INT i = 0; i < Cache_Size; i++ )
00099 {
00100 kernel_columns[i] = new DREAL[num_data];
00101 cache_index[i] = -2;
00102 }
00103 first_kernel_inx = 0;
00104
00105
00106
00107 for(INT i = 0; i < 3; i++ )
00108 {
00109 virt_columns[i] = new DREAL[num_virt_data];
00110 }
00111 first_virt_inx = 0;
00112
00113 diag_H = new DREAL[num_virt_data];
00114
00115 for(INT i = 0; i < num_virt_data; i++ )
00116 diag_H[i] = kernel_fce(i,i);
00117 }
00118
00119 CGMNPLib::~CGMNPLib()
00120 {
00121 for(INT i = 0; i < Cache_Size; i++ )
00122 delete[] kernel_columns[i];
00123
00124 for(INT i = 0; i < 3; i++ )
00125 delete[] virt_columns[i];
00126
00127 delete[] cache_index;
00128 delete[] kernel_columns;
00129
00130 delete[] diag_H;
00131 }
00132
00133
00134
00135
00136
00137 DREAL* CGMNPLib::get_kernel_col( INT a )
00138 {
00139 double *col_ptr;
00140 INT i;
00141 INT inx;
00142
00143 inx = -1;
00144 for( i=0; i < Cache_Size; i++ ) {
00145 if( cache_index[i] == a ) { inx = i; break; }
00146 }
00147
00148 if( inx != -1 ) {
00149 col_ptr = kernel_columns[inx];
00150 return( col_ptr );
00151 }
00152
00153 col_ptr = kernel_columns[first_kernel_inx];
00154 cache_index[first_kernel_inx] = a;
00155
00156 first_kernel_inx++;
00157 if( first_kernel_inx >= Cache_Size ) first_kernel_inx = 0;
00158
00159 for( i=0; i < m_num_data; i++ ) {
00160 col_ptr[i] = m_kernel->kernel(i,a);
00161 }
00162
00163 return( col_ptr );
00164 }
00165
00166
00167
00168
00169
00170 void CGMNPLib::get_indices2( INT *index, INT *c, INT i )
00171 {
00172 *index = i / (m_num_classes-1);
00173
00174 *c= (i % (m_num_classes-1))+1;
00175 if( *c>= m_vector_y[ *index ]) (*c)++;
00176
00177 return;
00178 }
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188 DREAL* CGMNPLib::get_col( INT a, INT b )
00189 {
00190 INT i;
00191 double *col_ptr;
00192 double *ker_ptr;
00193 double value;
00194 INT i1,c1,i2,c2;
00195
00196 col_ptr = virt_columns[first_virt_inx++];
00197 if( first_virt_inx >= 3 ) first_virt_inx = 0;
00198
00199 get_indices2( &i1, &c1, a );
00200 ker_ptr = (double*) get_kernel_col( i1 );
00201
00202 for( i=0; i < m_num_virt_data; i++ ) {
00203 get_indices2( &i2, &c2, i );
00204
00205 if( KDELTA4(m_vector_y[i1],m_vector_y[i2],c1,c2) ) {
00206 value = (+KDELTA(m_vector_y[i1],m_vector_y[i2])
00207 -KDELTA(m_vector_y[i1],c2)
00208 -KDELTA(m_vector_y[i2],c1)
00209 +KDELTA(c1,c2)
00210 )*(ker_ptr[i2]+1);
00211 }
00212 else
00213 {
00214 value = 0;
00215 }
00216
00217 if(a==i) value += m_reg_const;
00218
00219 col_ptr[i] = value;
00220 }
00221
00222 return( col_ptr );
00223 }
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236 int CGMNPLib::gmnp_imdm(double *vector_c,
00237 INT dim,
00238 INT tmax,
00239 double tolabs,
00240 double tolrel,
00241 double th,
00242 double *alpha,
00243 INT *ptr_t,
00244 double **ptr_History,
00245 INT verb)
00246 {
00247 double LB;
00248 double UB;
00249 double aHa, ac;
00250 double tmp, tmp1;
00251 double Huu, Huv, Hvv;
00252 double min_beta, beta;
00253 double max_improv, improv;
00254 double lambda;
00255 double *History;
00256 double *Ha;
00257 double *tmp_ptr;
00258 double *col_u, *col_v;
00259 INT u=0;
00260 INT v=0;
00261 INT new_u=0;
00262 INT i;
00263 INT t;
00264 INT History_size;
00265 int exitflag;
00266
00267
00268
00269
00270
00271 Ha = new double[dim];
00272 if( Ha == NULL ) SG_ERROR("Not enough memory.");
00273
00274 History_size = (tmax < HISTORY_BUF ) ? tmax+1 : HISTORY_BUF;
00275 History = new double[History_size*2];
00276 if( History == NULL ) SG_ERROR("Not enough memory.");
00277
00278
00279 for( tmp1 = PLUS_INF, i = 0; i < dim; i++ ) {
00280 tmp = 0.5*diag_H[i] + vector_c[i];
00281 if( tmp1 > tmp) {
00282 tmp1 = tmp;
00283 v = i;
00284 }
00285 }
00286
00287 col_v = (double*)get_col(v,-1);
00288
00289 for( min_beta = PLUS_INF, i = 0; i < dim; i++ )
00290 {
00291 alpha[i] = 0;
00292 Ha[i] = col_v[i];
00293
00294 beta = Ha[i] + vector_c[i];
00295 if( beta < min_beta ) {
00296 min_beta = beta;
00297 u = i;
00298 }
00299 }
00300
00301 alpha[v] = 1;
00302 aHa = diag_H[v];
00303 ac = vector_c[v];
00304
00305 UB = 0.5*aHa + ac;
00306 LB = min_beta - 0.5*aHa;
00307
00308 t = 0;
00309 History[INDEX(0,0,2)] = LB;
00310 History[INDEX(1,0,2)] = UB;
00311
00312 if( verb ) {
00313 SG_PRINT("Init: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
00314 UB, LB, UB-LB,(UB-LB)/UB);
00315 }
00316
00317
00318 if( UB-LB <= tolabs ) exitflag = 1;
00319 else if(UB-LB <= CMath::abs(UB)*tolrel ) exitflag = 2;
00320 else if(LB > th ) exitflag = 3;
00321 else exitflag = -1;
00322
00323
00324
00325
00326
00327 col_u = (double*)get_col(u,-1);
00328 while( exitflag == -1 )
00329 {
00330 t++;
00331
00332 col_v = (double*)get_col(v,u);
00333
00334
00335 Huu = diag_H[u];
00336 Hvv = diag_H[v];
00337 Huv = col_u[v];
00338
00339 lambda = (Ha[v]-Ha[u]+vector_c[v]-vector_c[u])/(alpha[v]*(Huu-2*Huv+Hvv));
00340 if( lambda < 0 ) lambda = 0; else if (lambda > 1) lambda = 1;
00341
00342 aHa = aHa + 2*alpha[v]*lambda*(Ha[u]-Ha[v])+
00343 lambda*lambda*alpha[v]*alpha[v]*(Huu-2*Huv+Hvv);
00344
00345 ac = ac + lambda*alpha[v]*(vector_c[u]-vector_c[v]);
00346
00347 tmp = alpha[v];
00348 alpha[u]=alpha[u]+lambda*alpha[v];
00349 alpha[v]=alpha[v]-lambda*alpha[v];
00350
00351 UB = 0.5*aHa + ac;
00352
00353
00354 for( min_beta = PLUS_INF, i = 0; i < dim; i++ )
00355 {
00356 Ha[i] = Ha[i] + lambda*tmp*(col_u[i] - col_v[i]);
00357
00358 beta = Ha[i]+ vector_c[i];
00359
00360 if( beta < min_beta )
00361 {
00362 new_u = i;
00363 min_beta = beta;
00364 }
00365 }
00366
00367 LB = min_beta - 0.5*aHa;
00368 u = new_u;
00369 col_u = (double*)get_col(u,-1);
00370
00371
00372 for( max_improv = MINUS_INF, i = 0; i < dim; i++ ) {
00373
00374 if( alpha[i] != 0 ) {
00375 beta = Ha[i] + vector_c[i];
00376
00377 if( beta >= min_beta ) {
00378
00379 tmp = diag_H[u] - 2*col_u[i] + diag_H[i];
00380 if( tmp != 0 ) {
00381 improv = (0.5*(beta-min_beta)*(beta-min_beta))/tmp;
00382
00383 if( improv > max_improv ) {
00384 max_improv = improv;
00385 v = i;
00386 }
00387 }
00388 }
00389 }
00390 }
00391
00392
00393 if( UB-LB <= tolabs ) exitflag = 1;
00394 else if( UB-LB <= CMath::abs(UB)*tolrel ) exitflag = 2;
00395 else if(LB > th ) exitflag = 3;
00396 else if(t >= tmax) exitflag = 0;
00397
00398
00399 if(verb && (t % verb) == 0 ) {
00400 SG_PRINT("%d: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
00401 t, UB, LB, UB-LB,(UB-LB)/UB);
00402 }
00403
00404
00405 if( t < History_size ) {
00406 History[INDEX(0,t,2)] = LB;
00407 History[INDEX(1,t,2)] = UB;
00408 }
00409 else {
00410 tmp_ptr = new double[(History_size+HISTORY_BUF)*2];
00411 if( tmp_ptr == NULL ) SG_ERROR("Not enough memory.");
00412 for( i = 0; i < History_size; i++ ) {
00413 tmp_ptr[INDEX(0,i,2)] = History[INDEX(0,i,2)];
00414 tmp_ptr[INDEX(1,i,2)] = History[INDEX(1,i,2)];
00415 }
00416 tmp_ptr[INDEX(0,t,2)] = LB;
00417 tmp_ptr[INDEX(1,t,2)] = UB;
00418
00419 History_size += HISTORY_BUF;
00420 delete[] History;
00421 History = tmp_ptr;
00422 }
00423 }
00424
00425
00426 if(verb && (t % verb) ) {
00427 SG_PRINT("exit: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
00428 UB, LB, UB-LB,(UB-LB)/UB);
00429 }
00430
00431
00432
00433
00434
00435 (*ptr_t) = t;
00436 (*ptr_History) = History;
00437
00438
00439 delete[] Ha;
00440
00441 return( exitflag );
00442 }
00443
00444
00445
00446
00447
00448 double CGMNPLib::kernel_fce( INT a, INT b )
00449 {
00450 double value;
00451 INT i1,c1,i2,c2;
00452
00453 get_indices2( &i1, &c1, a );
00454 get_indices2( &i2, &c2, b );
00455
00456 if( KDELTA4(m_vector_y[i1],m_vector_y[i2],c1,c2) ) {
00457 value = (+KDELTA(m_vector_y[i1],m_vector_y[i2])
00458 -KDELTA(m_vector_y[i1],c2)
00459 -KDELTA(m_vector_y[i2],c1)
00460 +KDELTA(c1,c2)
00461 )*(m_kernel->kernel( i1, i2 )+1);
00462 }
00463 else
00464 {
00465 value = 0;
00466 }
00467
00468 if(a==b) value += m_reg_const;
00469
00470 return( value );
00471 }