SHOGUN
v2.0.0
|
00001 /* 00002 * This program is free software; you can redistribute it and/or modify 00003 * it under the terms of the GNU General Public License as published by 00004 * the Free Software Foundation; either version 3 of the License, or 00005 * (at your option) any later version. 00006 * 00007 * libocas.c: Implementation of the OCAS solver for training 00008 * linear SVM classifiers. 00009 * 00010 * Copyright (C) 2008 Vojtech Franc, xfrancv@cmp.felk.cvut.cz 00011 * Soeren Sonnenburg, soeren.sonnenburg@first.fraunhofer.de 00012 *-------------------------------------------------------------------- */ 00013 00014 #include <stdlib.h> 00015 #include <string.h> 00016 #include <math.h> 00017 #include <sys/time.h> 00018 #include <time.h> 00019 #include <stdio.h> 00020 #include <stdint.h> 00021 00022 #include <shogun/lib/external/libocas.h> 00023 #include <shogun/lib/external/libocas_common.h> 00024 #include <shogun/lib/external/libqp.h> 00025 00026 namespace shogun 00027 { 00028 #define MU 0.1 /* must be from (0,1> 1..means that OCAS becomes equivalent to CPA */ 00029 /* see paper Franc&Sonneburg JMLR 2009 */ 00030 00031 static const uint32_t QPSolverMaxIter = 10000000; 00032 00033 static float64_t *H; 00034 static uint32_t BufSize; 00035 00036 /*---------------------------------------------------------------------- 00037 Returns pointer at i-th column of Hessian matrix. 00038 ----------------------------------------------------------------------*/ 00039 static const float64_t *get_col( uint32_t i) 00040 { 00041 return( &H[ BufSize*i ] ); 00042 } 00043 00044 /*---------------------------------------------------------------------- 00045 Returns time of the day in seconds. 00046 ----------------------------------------------------------------------*/ 00047 static float64_t get_time() 00048 { 00049 struct timeval tv; 00050 if (gettimeofday(&tv, NULL)==0) 00051 return tv.tv_sec+((float64_t)(tv.tv_usec))/1e6; 00052 else 00053 return 0.0; 00054 } 00055 00056 00057 /*---------------------------------------------------------------------- 00058 Linear binary Ocas-SVM solver with additinal contraint enforceing 00059 a subset of weights (indices of the weights given by num_nnw/nnw_idx) 00060 to be non-negative. 00061 00062 ----------------------------------------------------------------------*/ 00063 ocas_return_value_T svm_ocas_solver_nnw( 00064 float64_t C, 00065 uint32_t nData, 00066 uint32_t num_nnw, 00067 uint32_t* nnw_idx, 00068 float64_t TolRel, 00069 float64_t TolAbs, 00070 float64_t QPBound, 00071 float64_t MaxTime, 00072 uint32_t _BufSize, 00073 uint8_t Method, 00074 int (*add_pw_constr)(uint32_t, uint32_t, void*), 00075 void (*clip_neg_W)(uint32_t, uint32_t*, void*), 00076 void (*compute_W)(float64_t*, float64_t*, float64_t*, uint32_t, void*), 00077 float64_t (*update_W)(float64_t, void*), 00078 int (*add_new_cut)(float64_t*, uint32_t*, uint32_t, uint32_t, void*), 00079 int (*compute_output)(float64_t*, void* ), 00080 int (*sort)(float64_t*, float64_t*, uint32_t), 00081 void (*ocas_print)(ocas_return_value_T), 00082 void* user_data) 00083 { 00084 ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; 00085 float64_t *b, *alpha, *diag_H; 00086 float64_t *output, *old_output; 00087 float64_t xi, sq_norm_W, QPSolverTolRel, dot_prod_WoldW, sq_norm_oldW; 00088 float64_t A0, B0, GradVal, t, t1, t2, *Ci, *Bi, *hpf, *hpb; 00089 float64_t start_time, ocas_start_time; 00090 uint32_t cut_length; 00091 uint32_t i, *new_cut; 00092 uint32_t *I; 00093 /* uint8_t S = 1; */ 00094 libqp_state_T qp_exitflag; 00095 00096 float64_t max_cp_norm; 00097 float64_t max_b; 00098 float64_t _C[2]; 00099 uint8_t S[2]; 00100 00101 ocas_start_time = get_time(); 00102 ocas.qp_solver_time = 0; 00103 ocas.output_time = 0; 00104 ocas.sort_time = 0; 00105 ocas.add_time = 0; 00106 ocas.w_time = 0; 00107 ocas.print_time = 0; 00108 00109 BufSize = _BufSize; 00110 00111 QPSolverTolRel = TolRel*0.5; 00112 00113 H=NULL; 00114 b=NULL; 00115 alpha=NULL; 00116 new_cut=NULL; 00117 I=NULL; 00118 diag_H=NULL; 00119 output=NULL; 00120 old_output=NULL; 00121 hpf=NULL; 00122 hpb = NULL; 00123 Ci=NULL; 00124 Bi=NULL; 00125 00126 /* Hessian matrix contains dot product of normal vectors of selected cutting planes */ 00127 H = (float64_t*)LIBOCAS_CALLOC(BufSize*BufSize,sizeof(float64_t)); 00128 if(H == NULL) 00129 { 00130 ocas.exitflag=-2; 00131 goto cleanup; 00132 } 00133 00134 /* bias of cutting planes */ 00135 b = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t)); 00136 if(b == NULL) 00137 { 00138 ocas.exitflag=-2; 00139 goto cleanup; 00140 } 00141 00142 alpha = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t)); 00143 if(alpha == NULL) 00144 { 00145 ocas.exitflag=-2; 00146 goto cleanup; 00147 } 00148 00149 /* indices of examples which define a new cut */ 00150 new_cut = (uint32_t*)LIBOCAS_CALLOC(nData,sizeof(uint32_t)); 00151 if(new_cut == NULL) 00152 { 00153 ocas.exitflag=-2; 00154 goto cleanup; 00155 } 00156 00157 I = (uint32_t*)LIBOCAS_CALLOC(BufSize,sizeof(uint32_t)); 00158 if(I == NULL) 00159 { 00160 ocas.exitflag=-2; 00161 goto cleanup; 00162 } 00163 00164 for(i=0; i< BufSize; i++) I[i] = 2; 00165 00166 diag_H = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t)); 00167 if(diag_H == NULL) 00168 { 00169 ocas.exitflag=-2; 00170 goto cleanup; 00171 } 00172 00173 output = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t)); 00174 if(output == NULL) 00175 { 00176 ocas.exitflag=-2; 00177 goto cleanup; 00178 } 00179 00180 old_output = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t)); 00181 if(old_output == NULL) 00182 { 00183 ocas.exitflag=-2; 00184 goto cleanup; 00185 } 00186 00187 /* array of hinge points used in line-serach */ 00188 hpf = (float64_t*) LIBOCAS_CALLOC(nData, sizeof(hpf[0])); 00189 if(hpf == NULL) 00190 { 00191 ocas.exitflag=-2; 00192 goto cleanup; 00193 } 00194 00195 hpb = (float64_t*) LIBOCAS_CALLOC(nData, sizeof(hpb[0])); 00196 if(hpb == NULL) 00197 { 00198 ocas.exitflag=-2; 00199 goto cleanup; 00200 } 00201 00202 /* vectors Ci, Bi are used in the line search procedure */ 00203 Ci = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t)); 00204 if(Ci == NULL) 00205 { 00206 ocas.exitflag=-2; 00207 goto cleanup; 00208 } 00209 00210 Bi = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t)); 00211 if(Bi == NULL) 00212 { 00213 ocas.exitflag=-2; 00214 goto cleanup; 00215 } 00216 00217 /* initial cutting planes implementing the non-negativity constraints on W*/ 00218 _C[0]=10000000.0; 00219 _C[1]=C; 00220 S[0]=1; 00221 S[1]=1; 00222 for(i=0; i < num_nnw; i++) 00223 { 00224 if(add_pw_constr(nnw_idx[i],i, user_data) != 0) 00225 { 00226 ocas.exitflag=-2; 00227 goto cleanup; 00228 } 00229 diag_H[i] = 1.0; 00230 H[LIBOCAS_INDEX(i,i,BufSize)] = 1.0; 00231 I[i] = 1; 00232 } 00233 00234 max_cp_norm = 1; 00235 max_b = 0; 00236 00237 /* */ 00238 ocas.nCutPlanes = num_nnw; 00239 ocas.exitflag = 0; 00240 ocas.nIter = 0; 00241 00242 /* Compute initial value of Q_P assuming that W is zero vector.*/ 00243 sq_norm_W = 0; 00244 xi = nData; 00245 ocas.Q_P = 0.5*sq_norm_W + C*xi; 00246 ocas.Q_D = 0; 00247 00248 /* Compute the initial cutting plane */ 00249 cut_length = nData; 00250 for(i=0; i < nData; i++) 00251 new_cut[i] = i; 00252 00253 ocas.trn_err = nData; 00254 ocas.ocas_time = get_time() - ocas_start_time; 00255 /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, Q_P-Q_D/abs(Q_P)=%f\n", 00256 ocas.nIter,cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P)); 00257 */ 00258 ocas_print(ocas); 00259 00260 /* main loop */ 00261 while( ocas.exitflag == 0 ) 00262 { 00263 ocas.nIter++; 00264 00265 /* append a new cut to the buffer and update H */ 00266 b[ocas.nCutPlanes] = -(float64_t)cut_length; 00267 00268 max_b = LIBOCAS_MAX(max_b,(float64_t)cut_length); 00269 00270 start_time = get_time(); 00271 00272 if(add_new_cut( &H[LIBOCAS_INDEX(0,ocas.nCutPlanes,BufSize)], new_cut, cut_length, ocas.nCutPlanes, user_data ) != 0) 00273 { 00274 ocas.exitflag=-2; 00275 goto cleanup; 00276 } 00277 00278 ocas.add_time += get_time() - start_time; 00279 00280 /* copy new added row: H(ocas.nCutPlanes,ocas.nCutPlanes,1:ocas.nCutPlanes-1) = H(1:ocas.nCutPlanes-1:ocas.nCutPlanes)' */ 00281 diag_H[ocas.nCutPlanes] = H[LIBOCAS_INDEX(ocas.nCutPlanes,ocas.nCutPlanes,BufSize)]; 00282 for(i=0; i < ocas.nCutPlanes; i++) { 00283 H[LIBOCAS_INDEX(ocas.nCutPlanes,i,BufSize)] = H[LIBOCAS_INDEX(i,ocas.nCutPlanes,BufSize)]; 00284 } 00285 00286 max_cp_norm = LIBOCAS_MAX(max_cp_norm, sqrt(diag_H[ocas.nCutPlanes])); 00287 00288 00289 ocas.nCutPlanes++; 00290 00291 /* call inner QP solver */ 00292 start_time = get_time(); 00293 00294 /* compute upper bound on sum of dual variables associated with the positivity constraints */ 00295 _C[0] = sqrt((float64_t)nData)*(sqrt(C)*sqrt(max_b) + C*max_cp_norm); 00296 00297 /* qp_exitflag = libqp_splx_solver(&get_col, diag_H, b, &C, I, &S, alpha, 00298 ocas.nCutPlanes, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBOCAS_PLUS_INF,0);*/ 00299 qp_exitflag = libqp_splx_solver(&get_col, diag_H, b, _C, I, S, alpha, 00300 ocas.nCutPlanes, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBOCAS_PLUS_INF,0); 00301 00302 ocas.qp_exitflag = qp_exitflag.exitflag; 00303 00304 ocas.qp_solver_time += get_time() - start_time; 00305 ocas.Q_D = -qp_exitflag.QP; 00306 00307 ocas.nNZAlpha = 0; 00308 for(i=0; i < ocas.nCutPlanes; i++) { 00309 if( alpha[i] != 0) ocas.nNZAlpha++; 00310 } 00311 00312 sq_norm_oldW = sq_norm_W; 00313 start_time = get_time(); 00314 compute_W( &sq_norm_W, &dot_prod_WoldW, alpha, ocas.nCutPlanes, user_data ); 00315 clip_neg_W(num_nnw, nnw_idx, user_data); 00316 ocas.w_time += get_time() - start_time; 00317 00318 /* select a new cut */ 00319 switch( Method ) 00320 { 00321 /* cutting plane algorithm implemented in SVMperf and BMRM */ 00322 case 0: 00323 00324 start_time = get_time(); 00325 if( compute_output( output, user_data ) != 0) 00326 { 00327 ocas.exitflag=-2; 00328 goto cleanup; 00329 } 00330 ocas.output_time += get_time()-start_time; 00331 00332 xi = 0; 00333 cut_length = 0; 00334 ocas.trn_err = 0; 00335 for(i=0; i < nData; i++) 00336 { 00337 if(output[i] <= 0) ocas.trn_err++; 00338 00339 if(output[i] <= 1) { 00340 xi += 1 - output[i]; 00341 new_cut[cut_length] = i; 00342 cut_length++; 00343 } 00344 } 00345 ocas.Q_P = 0.5*sq_norm_W + C*xi; 00346 00347 ocas.ocas_time = get_time() - ocas_start_time; 00348 00349 /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, 1-Q_D/Q_P=%f, nza=%4d, err=%.2f%%, qpf=%d\n", 00350 ocas.nIter,cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P), 00351 ocas.nNZAlpha, 100*(float64_t)ocas.trn_err/(float64_t)nData, ocas.qp_exitflag ); 00352 */ 00353 00354 start_time = get_time(); 00355 ocas_print(ocas); 00356 ocas.print_time += get_time() - start_time; 00357 00358 break; 00359 00360 00361 /* Ocas strategy */ 00362 case 1: 00363 00364 /* Linesearch */ 00365 A0 = sq_norm_W -2*dot_prod_WoldW + sq_norm_oldW; 00366 B0 = dot_prod_WoldW - sq_norm_oldW; 00367 00368 memcpy( old_output, output, sizeof(float64_t)*nData ); 00369 00370 start_time = get_time(); 00371 if( compute_output( output, user_data ) != 0) 00372 { 00373 ocas.exitflag=-2; 00374 goto cleanup; 00375 } 00376 ocas.output_time += get_time()-start_time; 00377 00378 uint32_t num_hp = 0; 00379 GradVal = B0; 00380 for(i=0; i< nData; i++) { 00381 00382 Ci[i] = C*(1-old_output[i]); 00383 Bi[i] = C*(old_output[i] - output[i]); 00384 00385 float64_t val; 00386 if(Bi[i] != 0) 00387 val = -Ci[i]/Bi[i]; 00388 else 00389 val = -LIBOCAS_PLUS_INF; 00390 00391 if (val>0) 00392 { 00393 /* hpi[num_hp] = i;*/ 00394 hpb[num_hp] = Bi[i]; 00395 hpf[num_hp] = val; 00396 num_hp++; 00397 } 00398 00399 if( (Bi[i] < 0 && val > 0) || (Bi[i] > 0 && val <= 0)) 00400 GradVal += Bi[i]; 00401 00402 } 00403 00404 t = 0; 00405 if( GradVal < 0 ) 00406 { 00407 start_time = get_time(); 00408 /* if( sort(hpf, hpi, num_hp) != 0)*/ 00409 if( sort(hpf, hpb, num_hp) != 0 ) 00410 { 00411 ocas.exitflag=-2; 00412 goto cleanup; 00413 } 00414 ocas.sort_time += get_time() - start_time; 00415 00416 float64_t t_new, GradVal_new; 00417 i = 0; 00418 while( GradVal < 0 && i < num_hp ) 00419 { 00420 t_new = hpf[i]; 00421 GradVal_new = GradVal + LIBOCAS_ABS(hpb[i]) + A0*(t_new-t); 00422 00423 if( GradVal_new >= 0 ) 00424 { 00425 t = t + GradVal*(t-t_new)/(GradVal_new - GradVal); 00426 } 00427 else 00428 { 00429 t = t_new; 00430 i++; 00431 } 00432 00433 GradVal = GradVal_new; 00434 } 00435 } 00436 00437 /* 00438 t = hpf[0] - 1; 00439 i = 0; 00440 GradVal = t*A0 + Bsum; 00441 while( GradVal < 0 && i < num_hp && hpf[i] < LIBOCAS_PLUS_INF ) { 00442 t = hpf[i]; 00443 Bsum = Bsum + LIBOCAS_ABS(Bi[hpi[i]]); 00444 GradVal = t*A0 + Bsum; 00445 i++; 00446 } 00447 */ 00448 t = LIBOCAS_MAX(t,0); /* just sanity check; t < 0 should not ocure */ 00449 00450 /* this guarantees that the new solution will not violate the positivity constraints on W */ 00451 t = LIBOCAS_MIN(t,1); 00452 00453 t1 = t; /* new (best so far) W */ 00454 t2 = t+MU*(1.0-t); /* new cutting plane */ 00455 /* t2 = t+(1.0-t)/10.0; */ 00456 00457 /* update W to be the best so far solution */ 00458 sq_norm_W = update_W( t1, user_data ); 00459 00460 /* select a new cut */ 00461 xi = 0; 00462 cut_length = 0; 00463 ocas.trn_err = 0; 00464 for(i=0; i < nData; i++ ) { 00465 00466 if( (old_output[i]*(1-t2) + t2*output[i]) <= 1 ) 00467 { 00468 new_cut[cut_length] = i; 00469 cut_length++; 00470 } 00471 00472 output[i] = old_output[i]*(1-t1) + t1*output[i]; 00473 00474 if( output[i] <= 1) xi += 1-output[i]; 00475 if( output[i] <= 0) ocas.trn_err++; 00476 00477 } 00478 00479 ocas.Q_P = 0.5*sq_norm_W + C*xi; 00480 00481 ocas.ocas_time = get_time() - ocas_start_time; 00482 00483 /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, 1-Q_D/Q_P=%f, nza=%4d, err=%.2f%%, qpf=%d\n", 00484 ocas.nIter, cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P), 00485 ocas.nNZAlpha, 100*(float64_t)ocas.trn_err/(float64_t)nData, ocas.qp_exitflag ); 00486 */ 00487 00488 start_time = get_time(); 00489 ocas_print(ocas); 00490 ocas.print_time += get_time() - start_time; 00491 00492 break; 00493 } 00494 00495 /* Stopping conditions */ 00496 if( ocas.Q_P - ocas.Q_D <= TolRel*LIBOCAS_ABS(ocas.Q_P)) ocas.exitflag = 1; 00497 if( ocas.Q_P - ocas.Q_D <= TolAbs) ocas.exitflag = 2; 00498 if( ocas.Q_P <= QPBound) ocas.exitflag = 3; 00499 if( MaxTime > 0 && ocas.ocas_time >= MaxTime) ocas.exitflag = 4; 00500 if(ocas.nCutPlanes >= BufSize) ocas.exitflag = -1; 00501 00502 } /* end of the main loop */ 00503 00504 cleanup: 00505 00506 LIBOCAS_FREE(H); 00507 LIBOCAS_FREE(b); 00508 LIBOCAS_FREE(alpha); 00509 LIBOCAS_FREE(new_cut); 00510 LIBOCAS_FREE(I); 00511 LIBOCAS_FREE(diag_H); 00512 LIBOCAS_FREE(output); 00513 LIBOCAS_FREE(old_output); 00514 LIBOCAS_FREE(hpf); 00515 /* LIBOCAS_FREE(hpi);*/ 00516 LIBOCAS_FREE(hpb); 00517 LIBOCAS_FREE(Ci); 00518 LIBOCAS_FREE(Bi); 00519 00520 ocas.ocas_time = get_time() - ocas_start_time; 00521 00522 return(ocas); 00523 } 00524 00525 00526 00527 /*---------------------------------------------------------------------- 00528 Linear binary Ocas-SVM solver. 00529 ----------------------------------------------------------------------*/ 00530 ocas_return_value_T svm_ocas_solver( 00531 float64_t C, 00532 uint32_t nData, 00533 float64_t TolRel, 00534 float64_t TolAbs, 00535 float64_t QPBound, 00536 float64_t MaxTime, 00537 uint32_t _BufSize, 00538 uint8_t Method, 00539 void (*compute_W)(float64_t*, float64_t*, float64_t*, uint32_t, void*), 00540 float64_t (*update_W)(float64_t, void*), 00541 int (*add_new_cut)(float64_t*, uint32_t*, uint32_t, uint32_t, void*), 00542 int (*compute_output)(float64_t*, void* ), 00543 int (*sort)(float64_t*, float64_t*, uint32_t), 00544 void (*ocas_print)(ocas_return_value_T), 00545 void* user_data) 00546 { 00547 ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; 00548 float64_t *b, *alpha, *diag_H; 00549 float64_t *output, *old_output; 00550 float64_t xi, sq_norm_W, QPSolverTolRel, dot_prod_WoldW, sq_norm_oldW; 00551 float64_t A0, B0, GradVal, t, t1, t2, *Ci, *Bi, *hpf, *hpb; 00552 float64_t start_time, ocas_start_time; 00553 uint32_t cut_length; 00554 uint32_t i, *new_cut; 00555 uint32_t *I; 00556 uint8_t S = 1; 00557 libqp_state_T qp_exitflag; 00558 00559 ocas_start_time = get_time(); 00560 ocas.qp_solver_time = 0; 00561 ocas.output_time = 0; 00562 ocas.sort_time = 0; 00563 ocas.add_time = 0; 00564 ocas.w_time = 0; 00565 ocas.print_time = 0; 00566 float64_t gap; 00567 00568 BufSize = _BufSize; 00569 00570 QPSolverTolRel = TolRel*0.5; 00571 00572 H=NULL; 00573 b=NULL; 00574 alpha=NULL; 00575 new_cut=NULL; 00576 I=NULL; 00577 diag_H=NULL; 00578 output=NULL; 00579 old_output=NULL; 00580 hpf=NULL; 00581 hpb = NULL; 00582 Ci=NULL; 00583 Bi=NULL; 00584 00585 /* Hessian matrix contains dot product of normal vectors of selected cutting planes */ 00586 H = (float64_t*)LIBOCAS_CALLOC(BufSize*BufSize,sizeof(float64_t)); 00587 if(H == NULL) 00588 { 00589 ocas.exitflag=-2; 00590 goto cleanup; 00591 } 00592 00593 /* bias of cutting planes */ 00594 b = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t)); 00595 if(b == NULL) 00596 { 00597 ocas.exitflag=-2; 00598 goto cleanup; 00599 } 00600 00601 alpha = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t)); 00602 if(alpha == NULL) 00603 { 00604 ocas.exitflag=-2; 00605 goto cleanup; 00606 } 00607 00608 /* indices of examples which define a new cut */ 00609 new_cut = (uint32_t*)LIBOCAS_CALLOC(nData,sizeof(uint32_t)); 00610 if(new_cut == NULL) 00611 { 00612 ocas.exitflag=-2; 00613 goto cleanup; 00614 } 00615 00616 I = (uint32_t*)LIBOCAS_CALLOC(BufSize,sizeof(uint32_t)); 00617 if(I == NULL) 00618 { 00619 ocas.exitflag=-2; 00620 goto cleanup; 00621 } 00622 00623 for(i=0; i< BufSize; i++) I[i] = 1; 00624 00625 diag_H = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t)); 00626 if(diag_H == NULL) 00627 { 00628 ocas.exitflag=-2; 00629 goto cleanup; 00630 } 00631 00632 output = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t)); 00633 if(output == NULL) 00634 { 00635 ocas.exitflag=-2; 00636 goto cleanup; 00637 } 00638 00639 old_output = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t)); 00640 if(old_output == NULL) 00641 { 00642 ocas.exitflag=-2; 00643 goto cleanup; 00644 } 00645 00646 /* array of hinge points used in line-serach */ 00647 hpf = (float64_t*) LIBOCAS_CALLOC(nData, sizeof(hpf[0])); 00648 if(hpf == NULL) 00649 { 00650 ocas.exitflag=-2; 00651 goto cleanup; 00652 } 00653 00654 hpb = (float64_t*) LIBOCAS_CALLOC(nData, sizeof(hpb[0])); 00655 if(hpb == NULL) 00656 { 00657 ocas.exitflag=-2; 00658 goto cleanup; 00659 } 00660 00661 /* vectors Ci, Bi are used in the line search procedure */ 00662 Ci = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t)); 00663 if(Ci == NULL) 00664 { 00665 ocas.exitflag=-2; 00666 goto cleanup; 00667 } 00668 00669 Bi = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t)); 00670 if(Bi == NULL) 00671 { 00672 ocas.exitflag=-2; 00673 goto cleanup; 00674 } 00675 00676 ocas.nCutPlanes = 0; 00677 ocas.exitflag = 0; 00678 ocas.nIter = 0; 00679 00680 /* Compute initial value of Q_P assuming that W is zero vector.*/ 00681 sq_norm_W = 0; 00682 xi = nData; 00683 ocas.Q_P = 0.5*sq_norm_W + C*xi; 00684 ocas.Q_D = 0; 00685 00686 /* Compute the initial cutting plane */ 00687 cut_length = nData; 00688 for(i=0; i < nData; i++) 00689 new_cut[i] = i; 00690 00691 gap=(ocas.Q_P-ocas.Q_D)/CMath::abs(ocas.Q_P); 00692 SG_SABS_PROGRESS(gap, -CMath::log10(gap), -CMath::log10(1), -CMath::log10(TolRel), 6); 00693 00694 ocas.trn_err = nData; 00695 ocas.ocas_time = get_time() - ocas_start_time; 00696 /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, Q_P-Q_D/abs(Q_P)=%f\n", 00697 ocas.nIter,cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P)); 00698 */ 00699 ocas_print(ocas); 00700 00701 /* main loop */ 00702 while( ocas.exitflag == 0 ) 00703 { 00704 ocas.nIter++; 00705 00706 /* append a new cut to the buffer and update H */ 00707 b[ocas.nCutPlanes] = -(float64_t)cut_length; 00708 00709 start_time = get_time(); 00710 00711 if(add_new_cut( &H[LIBOCAS_INDEX(0,ocas.nCutPlanes,BufSize)], new_cut, cut_length, ocas.nCutPlanes, user_data ) != 0) 00712 { 00713 ocas.exitflag=-2; 00714 goto cleanup; 00715 } 00716 00717 ocas.add_time += get_time() - start_time; 00718 00719 /* copy new added row: H(ocas.nCutPlanes,ocas.nCutPlanes,1:ocas.nCutPlanes-1) = H(1:ocas.nCutPlanes-1:ocas.nCutPlanes)' */ 00720 diag_H[ocas.nCutPlanes] = H[LIBOCAS_INDEX(ocas.nCutPlanes,ocas.nCutPlanes,BufSize)]; 00721 for(i=0; i < ocas.nCutPlanes; i++) { 00722 H[LIBOCAS_INDEX(ocas.nCutPlanes,i,BufSize)] = H[LIBOCAS_INDEX(i,ocas.nCutPlanes,BufSize)]; 00723 } 00724 00725 ocas.nCutPlanes++; 00726 00727 /* call inner QP solver */ 00728 start_time = get_time(); 00729 00730 qp_exitflag = libqp_splx_solver(&get_col, diag_H, b, &C, I, &S, alpha, 00731 ocas.nCutPlanes, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBOCAS_PLUS_INF,0); 00732 00733 ocas.qp_exitflag = qp_exitflag.exitflag; 00734 00735 ocas.qp_solver_time += get_time() - start_time; 00736 ocas.Q_D = -qp_exitflag.QP; 00737 00738 ocas.nNZAlpha = 0; 00739 for(i=0; i < ocas.nCutPlanes; i++) { 00740 if( alpha[i] != 0) ocas.nNZAlpha++; 00741 } 00742 00743 sq_norm_oldW = sq_norm_W; 00744 start_time = get_time(); 00745 compute_W( &sq_norm_W, &dot_prod_WoldW, alpha, ocas.nCutPlanes, user_data ); 00746 ocas.w_time += get_time() - start_time; 00747 00748 /* select a new cut */ 00749 switch( Method ) 00750 { 00751 /* cutting plane algorithm implemented in SVMperf and BMRM */ 00752 case 0: 00753 00754 start_time = get_time(); 00755 if( compute_output( output, user_data ) != 0) 00756 { 00757 ocas.exitflag=-2; 00758 goto cleanup; 00759 } 00760 ocas.output_time += get_time()-start_time; 00761 gap=(ocas.Q_P-ocas.Q_D)/CMath::abs(ocas.Q_P); 00762 SG_SABS_PROGRESS(gap, -CMath::log10(gap), -CMath::log10(1), -CMath::log10(TolRel), 6); 00763 00764 xi = 0; 00765 cut_length = 0; 00766 ocas.trn_err = 0; 00767 for(i=0; i < nData; i++) 00768 { 00769 if(output[i] <= 0) ocas.trn_err++; 00770 00771 if(output[i] <= 1) { 00772 xi += 1 - output[i]; 00773 new_cut[cut_length] = i; 00774 cut_length++; 00775 } 00776 } 00777 ocas.Q_P = 0.5*sq_norm_W + C*xi; 00778 00779 ocas.ocas_time = get_time() - ocas_start_time; 00780 00781 /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, 1-Q_D/Q_P=%f, nza=%4d, err=%.2f%%, qpf=%d\n", 00782 ocas.nIter,cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P), 00783 ocas.nNZAlpha, 100*(float64_t)ocas.trn_err/(float64_t)nData, ocas.qp_exitflag ); 00784 */ 00785 00786 start_time = get_time(); 00787 ocas_print(ocas); 00788 ocas.print_time += get_time() - start_time; 00789 00790 break; 00791 00792 00793 /* Ocas strategy */ 00794 case 1: 00795 00796 /* Linesearch */ 00797 A0 = sq_norm_W -2*dot_prod_WoldW + sq_norm_oldW; 00798 B0 = dot_prod_WoldW - sq_norm_oldW; 00799 00800 memcpy( old_output, output, sizeof(float64_t)*nData ); 00801 00802 start_time = get_time(); 00803 if( compute_output( output, user_data ) != 0) 00804 { 00805 ocas.exitflag=-2; 00806 goto cleanup; 00807 } 00808 ocas.output_time += get_time()-start_time; 00809 00810 uint32_t num_hp = 0; 00811 GradVal = B0; 00812 for(i=0; i< nData; i++) { 00813 00814 Ci[i] = C*(1-old_output[i]); 00815 Bi[i] = C*(old_output[i] - output[i]); 00816 00817 float64_t val; 00818 if(Bi[i] != 0) 00819 val = -Ci[i]/Bi[i]; 00820 else 00821 val = -LIBOCAS_PLUS_INF; 00822 00823 if (val>0) 00824 { 00825 /* hpi[num_hp] = i;*/ 00826 hpb[num_hp] = Bi[i]; 00827 hpf[num_hp] = val; 00828 num_hp++; 00829 } 00830 00831 if( (Bi[i] < 0 && val > 0) || (Bi[i] > 0 && val <= 0)) 00832 GradVal += Bi[i]; 00833 00834 } 00835 00836 t = 0; 00837 if( GradVal < 0 ) 00838 { 00839 start_time = get_time(); 00840 /* if( sort(hpf, hpi, num_hp) != 0)*/ 00841 if( sort(hpf, hpb, num_hp) != 0 ) 00842 { 00843 ocas.exitflag=-2; 00844 goto cleanup; 00845 } 00846 ocas.sort_time += get_time() - start_time; 00847 00848 float64_t t_new, GradVal_new; 00849 i = 0; 00850 while( GradVal < 0 && i < num_hp ) 00851 { 00852 t_new = hpf[i]; 00853 GradVal_new = GradVal + LIBOCAS_ABS(hpb[i]) + A0*(t_new-t); 00854 00855 if( GradVal_new >= 0 ) 00856 { 00857 t = t + GradVal*(t-t_new)/(GradVal_new - GradVal); 00858 } 00859 else 00860 { 00861 t = t_new; 00862 i++; 00863 } 00864 00865 GradVal = GradVal_new; 00866 } 00867 } 00868 00869 /* 00870 t = hpf[0] - 1; 00871 i = 0; 00872 GradVal = t*A0 + Bsum; 00873 while( GradVal < 0 && i < num_hp && hpf[i] < LIBOCAS_PLUS_INF ) { 00874 t = hpf[i]; 00875 Bsum = Bsum + LIBOCAS_ABS(Bi[hpi[i]]); 00876 GradVal = t*A0 + Bsum; 00877 i++; 00878 } 00879 */ 00880 t = LIBOCAS_MAX(t,0); /* just sanity check; t < 0 should not ocure */ 00881 00882 t1 = t; /* new (best so far) W */ 00883 t2 = t+MU*(1.0-t); /* new cutting plane */ 00884 /* t2 = t+(1.0-t)/10.0; */ 00885 00886 /* update W to be the best so far solution */ 00887 sq_norm_W = update_W( t1, user_data ); 00888 00889 /* select a new cut */ 00890 xi = 0; 00891 cut_length = 0; 00892 ocas.trn_err = 0; 00893 for(i=0; i < nData; i++ ) { 00894 00895 if( (old_output[i]*(1-t2) + t2*output[i]) <= 1 ) 00896 { 00897 new_cut[cut_length] = i; 00898 cut_length++; 00899 } 00900 00901 output[i] = old_output[i]*(1-t1) + t1*output[i]; 00902 00903 if( output[i] <= 1) xi += 1-output[i]; 00904 if( output[i] <= 0) ocas.trn_err++; 00905 00906 } 00907 00908 ocas.Q_P = 0.5*sq_norm_W + C*xi; 00909 00910 ocas.ocas_time = get_time() - ocas_start_time; 00911 00912 /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, 1-Q_D/Q_P=%f, nza=%4d, err=%.2f%%, qpf=%d\n", 00913 ocas.nIter, cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P), 00914 ocas.nNZAlpha, 100*(float64_t)ocas.trn_err/(float64_t)nData, ocas.qp_exitflag ); 00915 */ 00916 00917 start_time = get_time(); 00918 ocas_print(ocas); 00919 ocas.print_time += get_time() - start_time; 00920 00921 break; 00922 } 00923 00924 /* Stopping conditions */ 00925 if( ocas.Q_P - ocas.Q_D <= TolRel*LIBOCAS_ABS(ocas.Q_P)) ocas.exitflag = 1; 00926 if( ocas.Q_P - ocas.Q_D <= TolAbs) ocas.exitflag = 2; 00927 if( ocas.Q_P <= QPBound) ocas.exitflag = 3; 00928 if( MaxTime > 0 && ocas.ocas_time >= MaxTime) ocas.exitflag = 4; 00929 if(ocas.nCutPlanes >= BufSize) ocas.exitflag = -1; 00930 00931 } /* end of the main loop */ 00932 00933 cleanup: 00934 00935 LIBOCAS_FREE(H); 00936 LIBOCAS_FREE(b); 00937 LIBOCAS_FREE(alpha); 00938 LIBOCAS_FREE(new_cut); 00939 LIBOCAS_FREE(I); 00940 LIBOCAS_FREE(diag_H); 00941 LIBOCAS_FREE(output); 00942 LIBOCAS_FREE(old_output); 00943 LIBOCAS_FREE(hpf); 00944 /* LIBOCAS_FREE(hpi);*/ 00945 LIBOCAS_FREE(hpb); 00946 LIBOCAS_FREE(Ci); 00947 LIBOCAS_FREE(Bi); 00948 00949 ocas.ocas_time = get_time() - ocas_start_time; 00950 00951 return(ocas); 00952 } 00953 00954 00955 /*---------------------------------------------------------------------- 00956 Binary linear Ocas-SVM solver which allows using different C for each 00957 training example. 00958 ----------------------------------------------------------------------*/ 00959 ocas_return_value_T svm_ocas_solver_difC( 00960 float64_t *C, 00961 uint32_t nData, 00962 float64_t TolRel, 00963 float64_t TolAbs, 00964 float64_t QPBound, 00965 float64_t MaxTime, 00966 uint32_t _BufSize, 00967 uint8_t Method, 00968 void (*compute_W)(float64_t*, float64_t*, float64_t*, uint32_t, void*), 00969 float64_t (*update_W)(float64_t, void*), 00970 int (*add_new_cut)(float64_t*, uint32_t*, uint32_t, uint32_t, void*), 00971 int (*compute_output)(float64_t*, void* ), 00972 int (*sort)(float64_t*, float64_t*, uint32_t), 00973 void (*ocas_print)(ocas_return_value_T), 00974 void* user_data) 00975 { 00976 ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; 00977 float64_t *b, *alpha, *diag_H; 00978 float64_t *output, *old_output; 00979 float64_t xi, sq_norm_W, QPSolverTolRel, dot_prod_WoldW, sq_norm_oldW; 00980 float64_t A0, B0, GradVal, t, t1, t2, *Ci, *Bi, *hpf, *hpb; 00981 float64_t start_time, ocas_start_time; 00982 float64_t qp_b = 1.0; 00983 float64_t new_b; 00984 uint32_t cut_length; 00985 uint32_t i, *new_cut; 00986 uint32_t *I; 00987 uint8_t S = 1; 00988 libqp_state_T qp_exitflag; 00989 00990 ocas_start_time = get_time(); 00991 ocas.qp_solver_time = 0; 00992 ocas.output_time = 0; 00993 ocas.sort_time = 0; 00994 ocas.add_time = 0; 00995 ocas.w_time = 0; 00996 ocas.print_time = 0; 00997 00998 BufSize = _BufSize; 00999 01000 QPSolverTolRel = TolRel*0.5; 01001 01002 H=NULL; 01003 b=NULL; 01004 alpha=NULL; 01005 new_cut=NULL; 01006 I=NULL; 01007 diag_H=NULL; 01008 output=NULL; 01009 old_output=NULL; 01010 hpf=NULL; 01011 hpb = NULL; 01012 Ci=NULL; 01013 Bi=NULL; 01014 01015 /* Hessian matrix contains dot product of normal vectors of selected cutting planes */ 01016 H = (float64_t*)LIBOCAS_CALLOC(BufSize*BufSize,sizeof(float64_t)); 01017 if(H == NULL) 01018 { 01019 ocas.exitflag=-2; 01020 goto cleanup; 01021 } 01022 01023 /* bias of cutting planes */ 01024 b = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t)); 01025 if(b == NULL) 01026 { 01027 ocas.exitflag=-2; 01028 goto cleanup; 01029 } 01030 01031 alpha = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t)); 01032 if(alpha == NULL) 01033 { 01034 ocas.exitflag=-2; 01035 goto cleanup; 01036 } 01037 01038 /* indices of examples which define a new cut */ 01039 new_cut = (uint32_t*)LIBOCAS_CALLOC(nData,sizeof(uint32_t)); 01040 if(new_cut == NULL) 01041 { 01042 ocas.exitflag=-2; 01043 goto cleanup; 01044 } 01045 01046 I = (uint32_t*)LIBOCAS_CALLOC(BufSize,sizeof(uint32_t)); 01047 if(I == NULL) 01048 { 01049 ocas.exitflag=-2; 01050 goto cleanup; 01051 } 01052 01053 for(i=0; i< BufSize; i++) I[i] = 1; 01054 01055 diag_H = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t)); 01056 if(diag_H == NULL) 01057 { 01058 ocas.exitflag=-2; 01059 goto cleanup; 01060 } 01061 01062 output = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t)); 01063 if(output == NULL) 01064 { 01065 ocas.exitflag=-2; 01066 goto cleanup; 01067 } 01068 01069 old_output = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t)); 01070 if(old_output == NULL) 01071 { 01072 ocas.exitflag=-2; 01073 goto cleanup; 01074 } 01075 01076 /* array of hinge points used in line-serach */ 01077 hpf = (float64_t*) LIBOCAS_CALLOC(nData, sizeof(hpf[0])); 01078 if(hpf == NULL) 01079 { 01080 ocas.exitflag=-2; 01081 goto cleanup; 01082 } 01083 01084 hpb = (float64_t*) LIBOCAS_CALLOC(nData, sizeof(hpb[0])); 01085 if(hpb == NULL) 01086 { 01087 ocas.exitflag=-2; 01088 goto cleanup; 01089 } 01090 01091 /* vectors Ci, Bi are used in the line search procedure */ 01092 Ci = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t)); 01093 if(Ci == NULL) 01094 { 01095 ocas.exitflag=-2; 01096 goto cleanup; 01097 } 01098 01099 Bi = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t)); 01100 if(Bi == NULL) 01101 { 01102 ocas.exitflag=-2; 01103 goto cleanup; 01104 } 01105 01106 ocas.nCutPlanes = 0; 01107 ocas.exitflag = 0; 01108 ocas.nIter = 0; 01109 01110 /* Compute initial value of Q_P assuming that W is zero vector.*/ 01111 sq_norm_W = 0; 01112 xi = nData; 01113 /* ocas.Q_P = 0.5*sq_norm_W + C*xi;*/ 01114 ocas.Q_D = 0; 01115 01116 /* Compute the initial cutting plane */ 01117 cut_length = nData; 01118 new_b = 0; 01119 for(i=0; i < nData; i++) 01120 { 01121 new_cut[i] = i; 01122 new_b += C[i]; 01123 } 01124 01125 ocas.Q_P = 0.5*sq_norm_W + new_b; 01126 01127 01128 ocas.trn_err = nData; 01129 ocas.ocas_time = get_time() - ocas_start_time; 01130 /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, Q_P-Q_D/abs(Q_P)=%f\n", 01131 ocas.nIter,cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P)); 01132 */ 01133 ocas_print(ocas); 01134 01135 /* main loop */ 01136 while( ocas.exitflag == 0 ) 01137 { 01138 ocas.nIter++; 01139 01140 /* append a new cut to the buffer and update H */ 01141 /* b[ocas.nCutPlanes] = -(float64_t)cut_length*C;*/ 01142 b[ocas.nCutPlanes] = -new_b; 01143 01144 start_time = get_time(); 01145 01146 if(add_new_cut( &H[LIBOCAS_INDEX(0,ocas.nCutPlanes,BufSize)], new_cut, cut_length, ocas.nCutPlanes, user_data ) != 0) 01147 { 01148 ocas.exitflag=-2; 01149 goto cleanup; 01150 } 01151 01152 ocas.add_time += get_time() - start_time; 01153 01154 /* copy new added row: H(ocas.nCutPlanes,ocas.nCutPlanes,1:ocas.nCutPlanes-1) = H(1:ocas.nCutPlanes-1:ocas.nCutPlanes)' */ 01155 diag_H[ocas.nCutPlanes] = H[LIBOCAS_INDEX(ocas.nCutPlanes,ocas.nCutPlanes,BufSize)]; 01156 for(i=0; i < ocas.nCutPlanes; i++) { 01157 H[LIBOCAS_INDEX(ocas.nCutPlanes,i,BufSize)] = H[LIBOCAS_INDEX(i,ocas.nCutPlanes,BufSize)]; 01158 } 01159 01160 ocas.nCutPlanes++; 01161 01162 /* call inner QP solver */ 01163 start_time = get_time(); 01164 01165 /* qp_exitflag = libqp_splx_solver(&get_col, diag_H, b, &C, I, &S, alpha,*/ 01166 /* ocas.nCutPlanes, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBOCAS_PLUS_INF,0);*/ 01167 qp_exitflag = libqp_splx_solver(&get_col, diag_H, b, &qp_b, I, &S, alpha, 01168 ocas.nCutPlanes, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBOCAS_PLUS_INF,0); 01169 01170 ocas.qp_exitflag = qp_exitflag.exitflag; 01171 01172 ocas.qp_solver_time += get_time() - start_time; 01173 ocas.Q_D = -qp_exitflag.QP; 01174 01175 ocas.nNZAlpha = 0; 01176 for(i=0; i < ocas.nCutPlanes; i++) { 01177 if( alpha[i] != 0) ocas.nNZAlpha++; 01178 } 01179 01180 sq_norm_oldW = sq_norm_W; 01181 start_time = get_time(); 01182 compute_W( &sq_norm_W, &dot_prod_WoldW, alpha, ocas.nCutPlanes, user_data ); 01183 ocas.w_time += get_time() - start_time; 01184 01185 /* select a new cut */ 01186 switch( Method ) 01187 { 01188 /* cutting plane algorithm implemented in SVMperf and BMRM */ 01189 case 0: 01190 01191 start_time = get_time(); 01192 if( compute_output( output, user_data ) != 0) 01193 { 01194 ocas.exitflag=-2; 01195 goto cleanup; 01196 } 01197 ocas.output_time += get_time()-start_time; 01198 01199 xi = 0; 01200 cut_length = 0; 01201 ocas.trn_err = 0; 01202 new_b = 0; 01203 for(i=0; i < nData; i++) 01204 { 01205 if(output[i] <= 0) ocas.trn_err++; 01206 01207 /* if(output[i] <= 1) {*/ 01208 /* xi += 1 - output[i];*/ 01209 if(output[i] <= C[i]) { 01210 xi += C[i] - output[i]; 01211 new_cut[cut_length] = i; 01212 cut_length++; 01213 new_b += C[i]; 01214 } 01215 } 01216 /* ocas.Q_P = 0.5*sq_norm_W + C*xi;*/ 01217 ocas.Q_P = 0.5*sq_norm_W + xi; 01218 01219 ocas.ocas_time = get_time() - ocas_start_time; 01220 01221 /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, 1-Q_D/Q_P=%f, nza=%4d, err=%.2f%%, qpf=%d\n", 01222 ocas.nIter,cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P), 01223 ocas.nNZAlpha, 100*(float64_t)ocas.trn_err/(float64_t)nData, ocas.qp_exitflag ); 01224 */ 01225 01226 start_time = get_time(); 01227 ocas_print(ocas); 01228 ocas.print_time += get_time() - start_time; 01229 01230 break; 01231 01232 01233 /* Ocas strategy */ 01234 case 1: 01235 01236 /* Linesearch */ 01237 A0 = sq_norm_W -2*dot_prod_WoldW + sq_norm_oldW; 01238 B0 = dot_prod_WoldW - sq_norm_oldW; 01239 01240 memcpy( old_output, output, sizeof(float64_t)*nData ); 01241 01242 start_time = get_time(); 01243 if( compute_output( output, user_data ) != 0) 01244 { 01245 ocas.exitflag=-2; 01246 goto cleanup; 01247 } 01248 ocas.output_time += get_time()-start_time; 01249 01250 uint32_t num_hp = 0; 01251 GradVal = B0; 01252 for(i=0; i< nData; i++) { 01253 01254 /* Ci[i] = C*(1-old_output[i]);*/ 01255 /* Bi[i] = C*(old_output[i] - output[i]);*/ 01256 Ci[i] = (C[i]-old_output[i]); 01257 Bi[i] = old_output[i] - output[i]; 01258 01259 float64_t val; 01260 if(Bi[i] != 0) 01261 val = -Ci[i]/Bi[i]; 01262 else 01263 val = -LIBOCAS_PLUS_INF; 01264 01265 if (val>0) 01266 { 01267 /* hpi[num_hp] = i;*/ 01268 hpb[num_hp] = Bi[i]; 01269 hpf[num_hp] = val; 01270 num_hp++; 01271 } 01272 01273 if( (Bi[i] < 0 && val > 0) || (Bi[i] > 0 && val <= 0)) 01274 GradVal += Bi[i]; 01275 01276 } 01277 01278 t = 0; 01279 if( GradVal < 0 ) 01280 { 01281 start_time = get_time(); 01282 /* if( sort(hpf, hpi, num_hp) != 0)*/ 01283 if( sort(hpf, hpb, num_hp) != 0 ) 01284 { 01285 ocas.exitflag=-2; 01286 goto cleanup; 01287 } 01288 ocas.sort_time += get_time() - start_time; 01289 01290 float64_t t_new, GradVal_new; 01291 i = 0; 01292 while( GradVal < 0 && i < num_hp ) 01293 { 01294 t_new = hpf[i]; 01295 GradVal_new = GradVal + LIBOCAS_ABS(hpb[i]) + A0*(t_new-t); 01296 01297 if( GradVal_new >= 0 ) 01298 { 01299 t = t + GradVal*(t-t_new)/(GradVal_new - GradVal); 01300 } 01301 else 01302 { 01303 t = t_new; 01304 i++; 01305 } 01306 01307 GradVal = GradVal_new; 01308 } 01309 } 01310 01311 /* 01312 t = hpf[0] - 1; 01313 i = 0; 01314 GradVal = t*A0 + Bsum; 01315 while( GradVal < 0 && i < num_hp && hpf[i] < LIBOCAS_PLUS_INF ) { 01316 t = hpf[i]; 01317 Bsum = Bsum + LIBOCAS_ABS(Bi[hpi[i]]); 01318 GradVal = t*A0 + Bsum; 01319 i++; 01320 } 01321 */ 01322 t = LIBOCAS_MAX(t,0); /* just sanity check; t < 0 should not ocure */ 01323 01324 t1 = t; /* new (best so far) W */ 01325 t2 = t+(1.0-t)*MU; /* new cutting plane */ 01326 /* t2 = t+(1.0-t)/10.0; new cutting plane */ 01327 01328 /* update W to be the best so far solution */ 01329 sq_norm_W = update_W( t1, user_data ); 01330 01331 /* select a new cut */ 01332 xi = 0; 01333 cut_length = 0; 01334 ocas.trn_err = 0; 01335 new_b = 0; 01336 for(i=0; i < nData; i++ ) { 01337 01338 /* if( (old_output[i]*(1-t2) + t2*output[i]) <= 1 ) */ 01339 if( (old_output[i]*(1-t2) + t2*output[i]) <= C[i] ) 01340 { 01341 new_cut[cut_length] = i; 01342 cut_length++; 01343 new_b += C[i]; 01344 } 01345 01346 output[i] = old_output[i]*(1-t1) + t1*output[i]; 01347 01348 /* if( output[i] <= 1) xi += 1-output[i];*/ 01349 if( output[i] <= C[i]) xi += C[i]-output[i]; 01350 if( output[i] <= 0) ocas.trn_err++; 01351 01352 } 01353 01354 /* ocas.Q_P = 0.5*sq_norm_W + C*xi;*/ 01355 ocas.Q_P = 0.5*sq_norm_W + xi; 01356 01357 ocas.ocas_time = get_time() - ocas_start_time; 01358 01359 /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, 1-Q_D/Q_P=%f, nza=%4d, err=%.2f%%, qpf=%d\n", 01360 ocas.nIter, cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P), 01361 ocas.nNZAlpha, 100*(float64_t)ocas.trn_err/(float64_t)nData, ocas.qp_exitflag ); 01362 */ 01363 01364 start_time = get_time(); 01365 ocas_print(ocas); 01366 ocas.print_time += get_time() - start_time; 01367 01368 break; 01369 } 01370 01371 /* Stopping conditions */ 01372 if( ocas.Q_P - ocas.Q_D <= TolRel*LIBOCAS_ABS(ocas.Q_P)) ocas.exitflag = 1; 01373 if( ocas.Q_P - ocas.Q_D <= TolAbs) ocas.exitflag = 2; 01374 if( ocas.Q_P <= QPBound) ocas.exitflag = 3; 01375 if( MaxTime > 0 && ocas.ocas_time >= MaxTime) ocas.exitflag = 4; 01376 if(ocas.nCutPlanes >= BufSize) ocas.exitflag = -1; 01377 01378 } /* end of the main loop */ 01379 01380 cleanup: 01381 01382 LIBOCAS_FREE(H); 01383 LIBOCAS_FREE(b); 01384 LIBOCAS_FREE(alpha); 01385 LIBOCAS_FREE(new_cut); 01386 LIBOCAS_FREE(I); 01387 LIBOCAS_FREE(diag_H); 01388 LIBOCAS_FREE(output); 01389 LIBOCAS_FREE(old_output); 01390 LIBOCAS_FREE(hpf); 01391 /* LIBOCAS_FREE(hpi);*/ 01392 LIBOCAS_FREE(hpb); 01393 LIBOCAS_FREE(Ci); 01394 LIBOCAS_FREE(Bi); 01395 01396 ocas.ocas_time = get_time() - ocas_start_time; 01397 01398 return(ocas); 01399 } 01400 01401 01402 01403 /*---------------------------------------------------------------------- 01404 Multiclass SVM-Ocas solver 01405 ----------------------------------------------------------------------*/ 01406 01407 /* Helper function needed by the multi-class SVM linesearch. 01408 01409 - This function finds a simplified representation of a piece-wise linear function 01410 by splitting the domain into intervals and fining active terms for these intevals */ 01411 static void findactive(float64_t *Theta, float64_t *SortedA, uint32_t *nSortedA, float64_t *A, float64_t *B, int n, 01412 int (*sort)(float64_t*, float64_t*, uint32_t)) 01413 { 01414 float64_t tmp, theta; 01415 uint32_t i, j, idx, idx2 = 0, start; 01416 01417 sort(A,B,n); 01418 01419 tmp = B[0]; 01420 idx = 0; 01421 i = 0; 01422 while( i < (uint32_t)n-1 && A[i] == A[i+1]) 01423 { 01424 if( B[i+1] > B[idx] ) 01425 { 01426 idx = i+1; 01427 tmp = B[i+1]; 01428 } 01429 i++; 01430 } 01431 01432 (*nSortedA) = 1; 01433 SortedA[0] = A[idx]; 01434 01435 while(1) 01436 { 01437 start = idx + 1; 01438 while( start < (uint32_t)n && A[idx] == A[start]) 01439 start++; 01440 01441 theta = LIBOCAS_PLUS_INF; 01442 for(j=start; j < (uint32_t)n; j++) 01443 { 01444 tmp = (B[j] - B[idx])/(A[idx]-A[j]); 01445 if( tmp < theta) 01446 { 01447 theta = tmp; 01448 idx2 = j; 01449 } 01450 } 01451 01452 if( theta < LIBOCAS_PLUS_INF) 01453 { 01454 Theta[(*nSortedA) - 1] = theta; 01455 SortedA[(*nSortedA)] = A[idx2]; 01456 (*nSortedA)++; 01457 idx = idx2; 01458 } 01459 else 01460 return; 01461 } 01462 } 01463 01464 01465 /*---------------------------------------------------------------------- 01466 Multiclass linear OCAS-SVM solver. 01467 ----------------------------------------------------------------------*/ 01468 ocas_return_value_T msvm_ocas_solver( 01469 float64_t C, 01470 float64_t *data_y, 01471 uint32_t nY, 01472 uint32_t nData, 01473 float64_t TolRel, 01474 float64_t TolAbs, 01475 float64_t QPBound, 01476 float64_t MaxTime, 01477 uint32_t _BufSize, 01478 uint8_t Method, 01479 void (*compute_W)(float64_t*, float64_t*, float64_t*, uint32_t, void*), 01480 float64_t (*update_W)(float64_t, void*), 01481 int (*add_new_cut)(float64_t*, uint32_t*, uint32_t, void*), 01482 int (*compute_output)(float64_t*, void* ), 01483 int (*sort)(float64_t*, float64_t*, uint32_t), 01484 void (*ocas_print)(ocas_return_value_T), 01485 void* user_data) 01486 { 01487 ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; 01488 float64_t *b, *alpha, *diag_H; 01489 float64_t *output, *old_output; 01490 float64_t xi, sq_norm_W, QPSolverTolRel, QPSolverTolAbs, dot_prod_WoldW, sq_norm_oldW; 01491 float64_t A0, B0, t, t1, t2, R, tmp, element_b, x; 01492 float64_t *A, *B, *theta, *Theta, *sortedA, *Add; 01493 float64_t start_time, ocas_start_time, grad_sum, grad, min_x = 0, old_x, old_grad; 01494 uint32_t i, y, y2, ypred = 0, *new_cut, cnt1, cnt2, j, nSortedA, idx; 01495 uint32_t *I; 01496 uint8_t S = 1; 01497 libqp_state_T qp_exitflag; 01498 01499 ocas_start_time = get_time(); 01500 ocas.qp_solver_time = 0; 01501 ocas.output_time = 0; 01502 ocas.sort_time = 0; 01503 ocas.add_time = 0; 01504 ocas.w_time = 0; 01505 ocas.print_time = 0; 01506 01507 BufSize = _BufSize; 01508 01509 QPSolverTolRel = TolRel*0.5; 01510 QPSolverTolAbs = TolAbs*0.5; 01511 01512 H=NULL; 01513 b=NULL; 01514 alpha=NULL; 01515 new_cut=NULL; 01516 I=NULL; 01517 diag_H=NULL; 01518 output=NULL; 01519 old_output=NULL; 01520 A = NULL; 01521 B = NULL; 01522 theta = NULL; 01523 Theta = NULL; 01524 sortedA = NULL; 01525 Add = NULL; 01526 01527 /* Hessian matrix contains dot product of normal vectors of selected cutting planes */ 01528 H = (float64_t*)LIBOCAS_CALLOC(BufSize*BufSize,sizeof(float64_t)); 01529 if(H == NULL) 01530 { 01531 ocas.exitflag=-2; 01532 goto cleanup; 01533 } 01534 01535 /* bias of cutting planes */ 01536 b = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t)); 01537 if(b == NULL) 01538 { 01539 ocas.exitflag=-2; 01540 goto cleanup; 01541 } 01542 01543 alpha = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t)); 01544 if(alpha == NULL) 01545 { 01546 ocas.exitflag=-2; 01547 goto cleanup; 01548 } 01549 01550 /* indices of examples which define a new cut */ 01551 new_cut = (uint32_t*)LIBOCAS_CALLOC(nData,sizeof(uint32_t)); 01552 if(new_cut == NULL) 01553 { 01554 ocas.exitflag=-2; 01555 goto cleanup; 01556 } 01557 01558 I = (uint32_t*)LIBOCAS_CALLOC(BufSize,sizeof(uint32_t)); 01559 if(I == NULL) 01560 { 01561 ocas.exitflag=-2; 01562 goto cleanup; 01563 } 01564 01565 for(i=0; i< BufSize; i++) 01566 I[i] = 1; 01567 01568 diag_H = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t)); 01569 if(diag_H == NULL) 01570 { 01571 ocas.exitflag=-2; 01572 goto cleanup; 01573 } 01574 01575 output = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t)); 01576 if(output == NULL) 01577 { 01578 ocas.exitflag=-2; 01579 goto cleanup; 01580 } 01581 01582 old_output = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t)); 01583 if(old_output == NULL) 01584 { 01585 ocas.exitflag=-2; 01586 goto cleanup; 01587 } 01588 01589 /* auxciliary variables used in the linesearch */ 01590 A = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t)); 01591 if(A == NULL) 01592 { 01593 ocas.exitflag=-2; 01594 goto cleanup; 01595 } 01596 01597 B = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t)); 01598 if(B == NULL) 01599 { 01600 ocas.exitflag=-2; 01601 goto cleanup; 01602 } 01603 01604 theta = (float64_t*)LIBOCAS_CALLOC(nY,sizeof(float64_t)); 01605 if(theta == NULL) 01606 { 01607 ocas.exitflag=-2; 01608 goto cleanup; 01609 } 01610 01611 sortedA = (float64_t*)LIBOCAS_CALLOC(nY,sizeof(float64_t)); 01612 if(sortedA == NULL) 01613 { 01614 ocas.exitflag=-2; 01615 goto cleanup; 01616 } 01617 01618 Theta = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t)); 01619 if(Theta == NULL) 01620 { 01621 ocas.exitflag=-2; 01622 goto cleanup; 01623 } 01624 01625 Add = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t)); 01626 if(Add == NULL) 01627 { 01628 ocas.exitflag=-2; 01629 goto cleanup; 01630 } 01631 01632 /* Set initial values*/ 01633 ocas.nCutPlanes = 0; 01634 ocas.exitflag = 0; 01635 ocas.nIter = 0; 01636 ocas.Q_D = 0; 01637 ocas.trn_err = nData; 01638 R = (float64_t)nData; 01639 sq_norm_W = 0; 01640 element_b = (float64_t)nData; 01641 ocas.Q_P = 0.5*sq_norm_W + C*R; 01642 01643 /* initial cutting plane */ 01644 for(i=0; i < nData; i++) 01645 { 01646 y2 = (uint32_t)data_y[i]; 01647 01648 if(y2 > 0) 01649 new_cut[i] = 0; 01650 else 01651 new_cut[i] = 1; 01652 01653 } 01654 01655 ocas.ocas_time = get_time() - ocas_start_time; 01656 01657 start_time = get_time(); 01658 ocas_print(ocas); 01659 ocas.print_time += get_time() - start_time; 01660 01661 /* main loop of the OCAS */ 01662 while( ocas.exitflag == 0 ) 01663 { 01664 ocas.nIter++; 01665 01666 /* append a new cut to the buffer and update H */ 01667 b[ocas.nCutPlanes] = -(float64_t)element_b; 01668 01669 start_time = get_time(); 01670 01671 if(add_new_cut( &H[LIBOCAS_INDEX(0,ocas.nCutPlanes,BufSize)], new_cut, ocas.nCutPlanes, user_data ) != 0) 01672 { 01673 ocas.exitflag=-2; 01674 goto cleanup; 01675 } 01676 01677 ocas.add_time += get_time() - start_time; 01678 01679 /* copy newly appended row: H(ocas.nCutPlanes,ocas.nCutPlanes,1:ocas.nCutPlanes-1) = H(1:ocas.nCutPlanes-1:ocas.nCutPlanes)' */ 01680 diag_H[ocas.nCutPlanes] = H[LIBOCAS_INDEX(ocas.nCutPlanes,ocas.nCutPlanes,BufSize)]; 01681 for(i=0; i < ocas.nCutPlanes; i++) 01682 { 01683 H[LIBOCAS_INDEX(ocas.nCutPlanes,i,BufSize)] = H[LIBOCAS_INDEX(i,ocas.nCutPlanes,BufSize)]; 01684 } 01685 01686 ocas.nCutPlanes++; 01687 01688 /* call inner QP solver */ 01689 start_time = get_time(); 01690 01691 qp_exitflag = libqp_splx_solver(&get_col, diag_H, b, &C, I, &S, alpha, 01692 ocas.nCutPlanes, QPSolverMaxIter, QPSolverTolAbs, QPSolverTolRel, -LIBOCAS_PLUS_INF,0); 01693 01694 ocas.qp_exitflag = qp_exitflag.exitflag; 01695 01696 ocas.qp_solver_time += get_time() - start_time; 01697 ocas.Q_D = -qp_exitflag.QP; 01698 01699 ocas.nNZAlpha = 0; 01700 for(i=0; i < ocas.nCutPlanes; i++) 01701 if( alpha[i] != 0) ocas.nNZAlpha++; 01702 01703 sq_norm_oldW = sq_norm_W; 01704 start_time = get_time(); 01705 compute_W( &sq_norm_W, &dot_prod_WoldW, alpha, ocas.nCutPlanes, user_data ); 01706 ocas.w_time += get_time() - start_time; 01707 01708 /* select a new cut */ 01709 switch( Method ) 01710 { 01711 /* cutting plane algorithm implemented in SVMperf and BMRM */ 01712 case 0: 01713 01714 start_time = get_time(); 01715 if( compute_output( output, user_data ) != 0) 01716 { 01717 ocas.exitflag=-2; 01718 goto cleanup; 01719 } 01720 ocas.output_time += get_time()-start_time; 01721 01722 /* the following loop computes: */ 01723 element_b = 0.0; /* element_b = R(old_W) - g'*old_W */ 01724 R = 0; /* R(W) = sum_i max_y ( [[y != y_i]] + (w_y- w_y_i)'*x_i ) */ 01725 ocas.trn_err = 0; /* trn_err = sum_i [[y != y_i ]] */ 01726 /* new_cut[i] = argmax_i ( [[y != y_i]] + (w_y- w_y_i)'*x_i ) */ 01727 for(i=0; i < nData; i++) 01728 { 01729 y2 = (uint32_t)data_y[i]; 01730 01731 for(xi=-LIBOCAS_PLUS_INF, y=0; y < nY; y++) 01732 { 01733 if(y2 != y && xi < output[LIBOCAS_INDEX(y,i,nY)]) 01734 { 01735 xi = output[LIBOCAS_INDEX(y,i,nY)]; 01736 ypred = y; 01737 } 01738 } 01739 01740 if(xi >= output[LIBOCAS_INDEX(y2,i,nY)]) 01741 ocas.trn_err ++; 01742 01743 xi = LIBOCAS_MAX(0,xi+1-output[LIBOCAS_INDEX(y2,i,nY)]); 01744 R += xi; 01745 if(xi > 0) 01746 { 01747 element_b++; 01748 new_cut[i] = ypred; 01749 } 01750 else 01751 new_cut[i] = y2; 01752 } 01753 01754 ocas.Q_P = 0.5*sq_norm_W + C*R; 01755 01756 ocas.ocas_time = get_time() - ocas_start_time; 01757 01758 start_time = get_time(); 01759 ocas_print(ocas); 01760 ocas.print_time += get_time() - start_time; 01761 01762 break; 01763 01764 /* The OCAS solver */ 01765 case 1: 01766 memcpy( old_output, output, sizeof(float64_t)*nData*nY ); 01767 01768 start_time = get_time(); 01769 if( compute_output( output, user_data ) != 0) 01770 { 01771 ocas.exitflag=-2; 01772 goto cleanup; 01773 } 01774 ocas.output_time += get_time()-start_time; 01775 01776 A0 = sq_norm_W - 2*dot_prod_WoldW + sq_norm_oldW; 01777 B0 = dot_prod_WoldW - sq_norm_oldW; 01778 01779 for(i=0; i < nData; i++) 01780 { 01781 y2 = (uint32_t)data_y[i]; 01782 01783 for(y=0; y < nY; y++) 01784 { 01785 A[LIBOCAS_INDEX(y,i,nY)] = C*(output[LIBOCAS_INDEX(y,i,nY)] - old_output[LIBOCAS_INDEX(y,i,nY)] 01786 + old_output[LIBOCAS_INDEX(y2,i,nY)] - output[LIBOCAS_INDEX(y2,i,nY)]); 01787 B[LIBOCAS_INDEX(y,i,nY)] = C*(old_output[LIBOCAS_INDEX(y,i,nY)] - old_output[LIBOCAS_INDEX(y2,i,nY)] 01788 + (float64_t)(y != y2)); 01789 } 01790 } 01791 01792 /* linesearch */ 01793 /* new_x = msvm_linesearch_mex(A0,B0,AA*C,BB*C);*/ 01794 01795 grad_sum = B0; 01796 cnt1 = 0; 01797 cnt2 = 0; 01798 for(i=0; i < nData; i++) 01799 { 01800 findactive(theta,sortedA,&nSortedA,&A[i*nY],&B[i*nY],nY,sort); 01801 01802 idx = 0; 01803 while( idx < nSortedA-1 && theta[idx] < 0 ) 01804 idx++; 01805 01806 grad_sum += sortedA[idx]; 01807 01808 for(j=idx; j < nSortedA-1; j++) 01809 { 01810 Theta[cnt1] = theta[j]; 01811 cnt1++; 01812 } 01813 01814 for(j=idx+1; j < nSortedA; j++) 01815 { 01816 Add[cnt2] = -sortedA[j-1]+sortedA[j]; 01817 cnt2++; 01818 } 01819 } 01820 01821 start_time = get_time(); 01822 sort(Theta,Add,cnt1); 01823 ocas.sort_time += get_time() - start_time; 01824 01825 grad = grad_sum; 01826 if(grad >= 0) 01827 { 01828 min_x = 0; 01829 } 01830 else 01831 { 01832 old_x = 0; 01833 old_grad = grad; 01834 01835 for(i=0; i < cnt1; i++) 01836 { 01837 x = Theta[i]; 01838 01839 grad = x*A0 + grad_sum; 01840 01841 if(grad >=0) 01842 { 01843 01844 min_x = (grad*old_x - old_grad*x)/(grad - old_grad); 01845 01846 break; 01847 } 01848 else 01849 { 01850 grad_sum = grad_sum + Add[i]; 01851 01852 grad = x*A0 + grad_sum; 01853 if( grad >= 0) 01854 { 01855 min_x = x; 01856 break; 01857 } 01858 } 01859 01860 old_grad = grad; 01861 old_x = x; 01862 } 01863 } 01864 /* end of the linesearch which outputs min_x */ 01865 01866 t = min_x; 01867 t1 = t; /* new (best so far) W */ 01868 t2 = t+(1.0-t)*MU; /* new cutting plane */ 01869 /* t2 = t+(1.0-t)/10.0; */ 01870 01871 /* update W to be the best so far solution */ 01872 sq_norm_W = update_W( t1, user_data ); 01873 01874 /* the following code computes a new cutting plane: */ 01875 element_b = 0.0; /* element_b = R(old_W) - g'*old_W */ 01876 /* new_cut[i] = argmax_i ( [[y != y_i]] + (w_y- w_y_i)'*x_i ) */ 01877 for(i=0; i < nData; i++) 01878 { 01879 y2 = (uint32_t)data_y[i]; 01880 01881 for(xi=-LIBOCAS_PLUS_INF, y=0; y < nY; y++) 01882 { 01883 tmp = old_output[LIBOCAS_INDEX(y,i,nY)]*(1-t2) + t2*output[LIBOCAS_INDEX(y,i,nY)]; 01884 if(y2 != y && xi < tmp) 01885 { 01886 xi = tmp; 01887 ypred = y; 01888 } 01889 } 01890 01891 tmp = old_output[LIBOCAS_INDEX(y2,i,nY)]*(1-t2) + t2*output[LIBOCAS_INDEX(y2,i,nY)]; 01892 xi = LIBOCAS_MAX(0,xi+1-tmp); 01893 if(xi > 0) 01894 { 01895 element_b++; 01896 new_cut[i] = ypred; 01897 } 01898 else 01899 new_cut[i] = y2; 01900 } 01901 01902 /* compute Risk, class. error and update outputs to correspond to the new W */ 01903 ocas.trn_err = 0; /* trn_err = sum_i [[y != y_i ]] */ 01904 R = 0; 01905 for(i=0; i < nData; i++) 01906 { 01907 y2 = (uint32_t)data_y[i]; 01908 01909 for(tmp=-LIBOCAS_PLUS_INF, y=0; y < nY; y++) 01910 { 01911 output[LIBOCAS_INDEX(y,i,nY)] = old_output[LIBOCAS_INDEX(y,i,nY)]*(1-t1) + t1*output[LIBOCAS_INDEX(y,i,nY)]; 01912 01913 if(y2 != y && tmp < output[LIBOCAS_INDEX(y,i,nY)]) 01914 { 01915 ypred = y; 01916 tmp = output[LIBOCAS_INDEX(y,i,nY)]; 01917 } 01918 } 01919 01920 R += LIBOCAS_MAX(0,1+tmp - output[LIBOCAS_INDEX(y2,i,nY)]); 01921 if( tmp >= output[LIBOCAS_INDEX(y2,i,nY)]) 01922 ocas.trn_err ++; 01923 } 01924 01925 ocas.Q_P = 0.5*sq_norm_W + C*R; 01926 01927 01928 /* get time and print status */ 01929 ocas.ocas_time = get_time() - ocas_start_time; 01930 01931 start_time = get_time(); 01932 ocas_print(ocas); 01933 ocas.print_time += get_time() - start_time; 01934 01935 break; 01936 01937 } 01938 01939 /* Stopping conditions */ 01940 if( ocas.Q_P - ocas.Q_D <= TolRel*LIBOCAS_ABS(ocas.Q_P)) ocas.exitflag = 1; 01941 if( ocas.Q_P - ocas.Q_D <= TolAbs) ocas.exitflag = 2; 01942 if( ocas.Q_P <= QPBound) ocas.exitflag = 3; 01943 if( MaxTime > 0 && ocas.ocas_time >= MaxTime) ocas.exitflag = 4; 01944 if(ocas.nCutPlanes >= BufSize) ocas.exitflag = -1; 01945 01946 } /* end of the main loop */ 01947 01948 cleanup: 01949 01950 LIBOCAS_FREE(H); 01951 LIBOCAS_FREE(b); 01952 LIBOCAS_FREE(alpha); 01953 LIBOCAS_FREE(new_cut); 01954 LIBOCAS_FREE(I); 01955 LIBOCAS_FREE(diag_H); 01956 LIBOCAS_FREE(output); 01957 LIBOCAS_FREE(old_output); 01958 LIBOCAS_FREE(A); 01959 LIBOCAS_FREE(B); 01960 LIBOCAS_FREE(theta); 01961 LIBOCAS_FREE(Theta); 01962 LIBOCAS_FREE(sortedA); 01963 LIBOCAS_FREE(Add); 01964 01965 ocas.ocas_time = get_time() - ocas_start_time; 01966 01967 return(ocas); 01968 } 01969 } 01970 01971