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
00064
00065
00066
00067
00068
00069
00070 #include <math.h>
00071 #include <stdlib.h>
00072 #include <stdio.h>
00073 #include <string.h>
00074 #include <stdint.h>
00075 #include <limits.h>
00076
00077 #include "classifier/svm/libocas_common.h"
00078 #include "classifier/svm/qpssvmlib.h"
00079
00080
00081
00082
00083
00084
00085
00086 int8_t qpssvm_solver(const void* (*get_col)(uint32_t),
00087 float64_t *diag_H,
00088 float64_t *f,
00089 float64_t b,
00090 uint16_t *I,
00091 float64_t *x,
00092 uint32_t n,
00093 uint32_t tmax,
00094 float64_t tolabs,
00095 float64_t tolrel,
00096 float64_t *QP,
00097 float64_t *QD,
00098 uint32_t verb)
00099 {
00100 float64_t *x_nequ;
00101 float64_t *d;
00102 float64_t *col_u, *col_v;
00103 float64_t LB;
00104 float64_t UB;
00105 float64_t tmp;
00106 float64_t improv;
00107 float64_t tmp_num;
00108 float64_t tmp_den=0;
00109 float64_t tau=0;
00110 float64_t delta;
00111 float64_t yu;
00112 uint32_t *inx;
00113 uint32_t *nk;
00114 uint32_t m;
00115 uint32_t t;
00116 uint32_t u=0;
00117 uint32_t v=0;
00118 uint32_t k;
00119 uint32_t i, j;
00120 int8_t exitflag;
00121
00122
00123
00124
00125
00126
00127 x_nequ=NULL;
00128 inx=NULL;
00129 nk=NULL;
00130 d=NULL;
00131
00132
00133 for( i=0, m=0; i < n; i++ ) m = CMath::max(m,(uint32_t) I[i]);
00134
00135
00136 x_nequ = (float64_t*) OCAS_CALLOC(m, sizeof(float64_t));
00137 if( x_nequ == NULL )
00138 {
00139 exitflag=-2;
00140 goto cleanup;
00141 }
00142
00143
00144 inx = (uint32_t*) OCAS_CALLOC(m*n, sizeof(uint32_t));
00145 if( inx == NULL )
00146 {
00147 exitflag=-2;
00148 goto cleanup;
00149 }
00150
00151 nk = (uint32_t*) OCAS_CALLOC(m, sizeof(uint32_t));
00152 if( nk == NULL )
00153 {
00154 exitflag=-2;
00155 goto cleanup;
00156 }
00157
00158 for( i=0; i < m; i++ ) x_nequ[i] = b;
00159 for( i=0; i < n; i++ ) {
00160 k = I[i]-1;
00161 x_nequ[k] -= x[i];
00162 inx[INDEX2(nk[k],k,n)] = i;
00163 nk[k]++;
00164 }
00165
00166
00167 d = (float64_t*) OCAS_CALLOC(n, sizeof(float64_t));
00168 if( d == NULL )
00169 {
00170 exitflag=-2;
00171 goto cleanup;
00172 }
00173
00174
00175 for( i=0; i < n; i++ ) {
00176 if( x[i] > 0 ) {
00177 col_u = (float64_t*)get_col(i);
00178 for( j=0; j < n; j++ ) {
00179 d[j] += col_u[j]*x[i];
00180 }
00181 }
00182 }
00183 for( i=0; i < n; i++ ) d[i] += f[i];
00184
00185
00186
00187 for( i=0, UB = 0, LB=0; i < n; i++) {
00188 UB += x[i]*(f[i]+d[i]);
00189 LB += x[i]*(f[i]-d[i]);
00190 }
00191 UB = 0.5*UB;
00192 LB = 0.5*LB;
00193
00194
00195
00196
00197
00198
00199
00200
00201 for( i=0; i < m; i++ ) {
00202 for( j=0, tmp = OCAS_PLUS_INF; j < nk[i]; j++ ) {
00203 tmp = CMath::min(tmp, d[inx[INDEX2(j,i,n)]]);
00204 }
00205 if( tmp < 0) LB += b*tmp;
00206 }
00207
00208 exitflag = 0;
00209 t = 0;
00210
00211
00212 while( (exitflag == 0) && (t < tmax))
00213 {
00214 t++;
00215
00216 exitflag = 1;
00217 for( k=0; k < m; k++ )
00218 {
00219
00220
00221
00222
00223
00224 for( j=0, tmp = OCAS_PLUS_INF, delta = 0; j < nk[k]; j++ ) {
00225 i = inx[INDEX2(j,k,n)];
00226 delta += x[i]*d[i];
00227 if( tmp > d[i] ) {
00228 tmp = d[i];
00229 u = i;
00230 }
00231 }
00232
00233
00234 if( d[u] < 0) yu = b; else yu = 0;
00235
00236
00237 delta -= yu*d[u];
00238
00239 if( delta > tolabs/m && delta > tolrel*CMath::abs(UB)/m)
00240 {
00241 exitflag = 0;
00242
00243 if( yu > 0 )
00244 {
00245 col_u = (float64_t*)get_col(u);
00246
00247 improv = -OCAS_PLUS_INF;
00248 for( j=0; j < nk[k]; j++ ) {
00249 i = inx[INDEX2(j,k,n)];
00250
00251
00252
00253 if(x[i] > 0) {
00254
00255 tmp_num = x[i]*(d[i] - d[u]);
00256 tmp_den = x[i]*x[i]*(diag_H[u] - 2*col_u[i] + diag_H[i]);
00257 if( tmp_den > 0 ) {
00258 if( tmp_num < tmp_den ) {
00259 tmp = tmp_num*tmp_num / tmp_den;
00260 } else {
00261 tmp = tmp_num - 0.5 * tmp_den;
00262 }
00263 }
00264 if( tmp > improv ) {
00265 improv = tmp;
00266 tau = CMath::min(1.0,tmp_num/tmp_den);
00267 v = i;
00268 }
00269 }
00270 }
00271
00272 tmp_num = -x_nequ[k]*d[u];
00273 if( tmp_num > 0 ) {
00274 tmp_den = x_nequ[k]*x_nequ[k]*diag_H[u];
00275 if( tmp_den > 0 ) {
00276 if( tmp_num < tmp_den ) {
00277 tmp = tmp_num*tmp_num / tmp_den;
00278 } else {
00279 tmp = tmp_num - 0.5 * tmp_den;
00280 }
00281 }
00282 } else {
00283 tmp = -OCAS_PLUS_INF;
00284 }
00285
00286 if( tmp > improv ) {
00287 tau = CMath::min(1.0,tmp_num/tmp_den);
00288 for( i = 0; i < n; i++ ) {
00289 d[i] += x_nequ[k]*tau*col_u[i];
00290 }
00291 x[u] += tau*x_nequ[k];
00292 x_nequ[k] -= tau*x_nequ[k];
00293
00294 } else {
00295
00296
00297 col_v = (float64_t*)get_col(v);
00298 for( i = 0; i < n; i++ ) {
00299 d[i] += x[v]*tau*(col_u[i]-col_v[i]);
00300 }
00301
00302 x[u] += tau*x[v];
00303 x[v] -= tau*x[v];
00304 }
00305 }
00306 else
00307 {
00308 improv = -OCAS_PLUS_INF;
00309 for( j=0; j < nk[k]; j++ ) {
00310 i = inx[INDEX2(j,k,n)];
00311
00312
00313
00314 if( x[i] > 0 && d[i] > 0) {
00315
00316 tmp_num = x[i]*d[i];
00317 tmp_den = x[i]*x[i]*diag_H[i];
00318 if( tmp_den > 0 ) {
00319 if( tmp_num < tmp_den ) {
00320 tmp = tmp_num*tmp_num / tmp_den;
00321 } else {
00322 tmp = tmp_num - 0.5 * tmp_den;
00323 }
00324 }
00325 if( tmp > improv ) {
00326 improv = tmp;
00327 tau = CMath::min(1.0,tmp_num/tmp_den);
00328 v = i;
00329 }
00330 }
00331 }
00332
00333
00334 col_v = (float64_t*)get_col(v);
00335 for( i = 0; i < n; i++ ) {
00336 d[i] -= x[v]*tau*col_v[i];
00337 }
00338
00339 x_nequ[k] += tau*x[v];
00340 x[v] -= tau*x[v];
00341 }
00342
00343 UB = UB - improv;
00344 }
00345
00346 }
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356 for( i=0, UB = 0, LB=0; i < n; i++) {
00357 UB += x[i]*(f[i]+d[i]);
00358 LB += x[i]*(f[i]-d[i]);
00359 }
00360 UB = 0.5*UB;
00361 LB = 0.5*LB;
00362
00363 for( k=0; k < m; k++ ) {
00364 for( j=0,tmp = OCAS_PLUS_INF; j < nk[k]; j++ ) {
00365 i = inx[INDEX2(j,k,n)];
00366
00367 tmp = CMath::min(tmp, d[i]);
00368 }
00369 if( tmp < 0) LB += b*tmp;
00370 }
00371
00372 if( verb > 0 && (exitflag > 0 || (t % verb)==0 ))
00373 {
00374 float64_t gap=(UB!=0) ? (UB-LB)/CMath::abs(UB) : 0;
00375 SG_SABS_PROGRESS(gap, -CMath::log10(gap), -CMath::log10(1), -CMath::log10(tolrel), 6);
00376 }
00377
00378 }
00379
00380
00381 if( UB-LB < tolabs ) exitflag = 1;
00382 else if(UB-LB < CMath::abs(UB)*tolrel ) exitflag = 2;
00383 else exitflag = 0;
00384
00385
00386
00387
00388 *QP = UB;
00389 *QD = LB;
00390
00391
00392
00393
00394 cleanup:
00395 OCAS_FREE( d );
00396 OCAS_FREE( inx );
00397 OCAS_FREE( nk );
00398 OCAS_FREE( x_nequ );
00399
00400 return( exitflag );
00401
00402 }
00403