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