00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
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/classifier/svm/libocas.h>
00023 #include <shogun/classifier/svm/libocas_common.h>
00024 #include <shogun/classifier/svm/libqp.h>
00025
00026 namespace shogun
00027 {
00028 #define LAMBDA 0.1
00029
00030 static const uint32_t QPSolverMaxIter = 10000000;
00031
00032 static float64_t *H;
00033 static uint32_t BufSize;
00034
00035
00036
00037
00038 static const float64_t *get_col( uint32_t i)
00039 {
00040 return( &H[ BufSize*i ] );
00041 }
00042
00043
00044
00045
00046 static float64_t get_time()
00047 {
00048 struct timeval tv;
00049 if (gettimeofday(&tv, NULL)==0)
00050 return tv.tv_sec+((float64_t)(tv.tv_usec))/1e6;
00051 else
00052 return 0.0;
00053 }
00054
00055
00056
00057
00058 ocas_return_value_T svm_ocas_solver(
00059 float64_t C,
00060 uint32_t nData,
00061 float64_t TolRel,
00062 float64_t TolAbs,
00063 float64_t QPBound,
00064 float64_t MaxTime,
00065 uint32_t _BufSize,
00066 uint8_t Method,
00067 void (*compute_W)(float64_t*, float64_t*, float64_t*, uint32_t, void*),
00068 float64_t (*update_W)(float64_t, void*),
00069 int (*add_new_cut)(float64_t*, uint32_t*, uint32_t, uint32_t, void*),
00070 int (*compute_output)(float64_t*, void* ),
00071 int (*sort)(float64_t*, float64_t*, uint32_t),
00072 void (*ocas_print)(ocas_return_value_T),
00073 void* user_data)
00074 {
00075 ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
00076 float64_t *b, *alpha, *diag_H;
00077 float64_t *output, *old_output;
00078 float64_t xi, sq_norm_W, QPSolverTolRel, dot_prod_WoldW, sq_norm_oldW;
00079 float64_t A0, B0, GradVal, t, t1, t2, *Ci, *Bi, *hpf, *hpb;
00080 float64_t start_time, ocas_start_time;
00081 uint32_t cut_length;
00082 uint32_t i, *new_cut;
00083 uint32_t *I;
00084 uint8_t S = 1;
00085 libqp_state_T qp_exitflag;
00086
00087 ocas_start_time = get_time();
00088 ocas.qp_solver_time = 0;
00089 ocas.output_time = 0;
00090 ocas.sort_time = 0;
00091 ocas.add_time = 0;
00092 ocas.w_time = 0;
00093 ocas.print_time = 0;
00094 float64_t gap;
00095
00096 BufSize = _BufSize;
00097
00098 QPSolverTolRel = TolRel*0.5;
00099
00100 H=NULL;
00101 b=NULL;
00102 alpha=NULL;
00103 new_cut=NULL;
00104 I=NULL;
00105 diag_H=NULL;
00106 output=NULL;
00107 old_output=NULL;
00108 hpf=NULL;
00109 hpb = NULL;
00110 Ci=NULL;
00111 Bi=NULL;
00112
00113
00114 H = (float64_t*)LIBOCAS_CALLOC(BufSize*BufSize,sizeof(float64_t));
00115 if(H == NULL)
00116 {
00117 ocas.exitflag=-2;
00118 goto cleanup;
00119 }
00120
00121
00122 b = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
00123 if(b == NULL)
00124 {
00125 ocas.exitflag=-2;
00126 goto cleanup;
00127 }
00128
00129 alpha = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
00130 if(alpha == NULL)
00131 {
00132 ocas.exitflag=-2;
00133 goto cleanup;
00134 }
00135
00136
00137 new_cut = (uint32_t*)LIBOCAS_CALLOC(nData,sizeof(uint32_t));
00138 if(new_cut == NULL)
00139 {
00140 ocas.exitflag=-2;
00141 goto cleanup;
00142 }
00143
00144 I = (uint32_t*)LIBOCAS_CALLOC(BufSize,sizeof(uint32_t));
00145 if(I == NULL)
00146 {
00147 ocas.exitflag=-2;
00148 goto cleanup;
00149 }
00150
00151 for(i=0; i< BufSize; i++) I[i] = 1;
00152
00153 diag_H = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
00154 if(diag_H == NULL)
00155 {
00156 ocas.exitflag=-2;
00157 goto cleanup;
00158 }
00159
00160 output = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
00161 if(output == NULL)
00162 {
00163 ocas.exitflag=-2;
00164 goto cleanup;
00165 }
00166
00167 old_output = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
00168 if(old_output == NULL)
00169 {
00170 ocas.exitflag=-2;
00171 goto cleanup;
00172 }
00173
00174
00175 hpf = (float64_t*) LIBOCAS_CALLOC(nData, sizeof(hpf[0]));
00176 if(hpf == NULL)
00177 {
00178 ocas.exitflag=-2;
00179 goto cleanup;
00180 }
00181
00182 hpb = (float64_t*) LIBOCAS_CALLOC(nData, sizeof(hpb[0]));
00183 if(hpb == NULL)
00184 {
00185 ocas.exitflag=-2;
00186 goto cleanup;
00187 }
00188
00189
00190 Ci = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
00191 if(Ci == NULL)
00192 {
00193 ocas.exitflag=-2;
00194 goto cleanup;
00195 }
00196
00197 Bi = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
00198 if(Bi == NULL)
00199 {
00200 ocas.exitflag=-2;
00201 goto cleanup;
00202 }
00203
00204 ocas.nCutPlanes = 0;
00205 ocas.exitflag = 0;
00206 ocas.nIter = 0;
00207
00208
00209 sq_norm_W = 0;
00210 xi = nData;
00211 ocas.Q_P = 0.5*sq_norm_W + C*xi;
00212 ocas.Q_D = 0;
00213
00214
00215 cut_length = nData;
00216 for(i=0; i < nData; i++)
00217 new_cut[i] = i;
00218
00219 gap=(ocas.Q_P-ocas.Q_D)/CMath::abs(ocas.Q_P);
00220 SG_SABS_PROGRESS(gap, -CMath::log10(gap), -CMath::log10(1), -CMath::log10(TolRel), 6);
00221
00222 ocas.trn_err = nData;
00223 ocas.ocas_time = get_time() - ocas_start_time;
00224
00225
00226
00227 ocas_print(ocas);
00228
00229
00230 while( ocas.exitflag == 0 )
00231 {
00232 ocas.nIter++;
00233
00234
00235 b[ocas.nCutPlanes] = -(float64_t)cut_length;
00236
00237 start_time = get_time();
00238
00239 if(add_new_cut( &H[LIBOCAS_INDEX(0,ocas.nCutPlanes,BufSize)], new_cut, cut_length, ocas.nCutPlanes, user_data ) != 0)
00240 {
00241 ocas.exitflag=-2;
00242 goto cleanup;
00243 }
00244
00245 ocas.add_time += get_time() - start_time;
00246
00247
00248 diag_H[ocas.nCutPlanes] = H[LIBOCAS_INDEX(ocas.nCutPlanes,ocas.nCutPlanes,BufSize)];
00249 for(i=0; i < ocas.nCutPlanes; i++) {
00250 H[LIBOCAS_INDEX(ocas.nCutPlanes,i,BufSize)] = H[LIBOCAS_INDEX(i,ocas.nCutPlanes,BufSize)];
00251 }
00252
00253 ocas.nCutPlanes++;
00254
00255
00256 start_time = get_time();
00257
00258 qp_exitflag = libqp_splx_solver(&get_col, diag_H, b, &C, I, &S, alpha,
00259 ocas.nCutPlanes, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBOCAS_PLUS_INF,0);
00260
00261 ocas.qp_exitflag = qp_exitflag.exitflag;
00262
00263 ocas.qp_solver_time += get_time() - start_time;
00264 ocas.Q_D = -qp_exitflag.QP;
00265
00266 ocas.nNZAlpha = 0;
00267 for(i=0; i < ocas.nCutPlanes; i++) {
00268 if( alpha[i] != 0) ocas.nNZAlpha++;
00269 }
00270
00271 sq_norm_oldW = sq_norm_W;
00272 start_time = get_time();
00273 compute_W( &sq_norm_W, &dot_prod_WoldW, alpha, ocas.nCutPlanes, user_data );
00274 ocas.w_time += get_time() - start_time;
00275
00276
00277 switch( Method )
00278 {
00279
00280 case 0:
00281
00282 start_time = get_time();
00283 if( compute_output( output, user_data ) != 0)
00284 {
00285 ocas.exitflag=-2;
00286 goto cleanup;
00287 }
00288 ocas.output_time += get_time()-start_time;
00289
00290 xi = 0;
00291 cut_length = 0;
00292 ocas.trn_err = 0;
00293 for(i=0; i < nData; i++)
00294 {
00295 if(output[i] <= 0) ocas.trn_err++;
00296
00297 if(output[i] <= 1) {
00298 xi += 1 - output[i];
00299 new_cut[cut_length] = i;
00300 cut_length++;
00301 }
00302 }
00303 ocas.Q_P = 0.5*sq_norm_W + C*xi;
00304
00305 gap=(ocas.Q_P-ocas.Q_D)/CMath::abs(ocas.Q_P);
00306 SG_SABS_PROGRESS(gap, -CMath::log10(gap), -CMath::log10(1), -CMath::log10(TolRel), 6);
00307
00308
00309
00310
00311
00312 start_time = get_time();
00313 ocas_print(ocas);
00314 ocas.print_time += get_time() - start_time;
00315
00316 break;
00317
00318
00319
00320 case 1:
00321
00322
00323 A0 = sq_norm_W -2*dot_prod_WoldW + sq_norm_oldW;
00324 B0 = dot_prod_WoldW - sq_norm_oldW;
00325
00326 memcpy( old_output, output, sizeof(float64_t)*nData );
00327
00328 start_time = get_time();
00329 if( compute_output( output, user_data ) != 0)
00330 {
00331 ocas.exitflag=-2;
00332 goto cleanup;
00333 }
00334 ocas.output_time += get_time()-start_time;
00335
00336 uint32_t num_hp = 0;
00337 GradVal = B0;
00338 for(i=0; i< nData; i++) {
00339
00340 Ci[i] = C*(1-old_output[i]);
00341 Bi[i] = C*(old_output[i] - output[i]);
00342
00343 float64_t val;
00344 if(Bi[i] != 0)
00345 val = -Ci[i]/Bi[i];
00346 else
00347 val = -LIBOCAS_PLUS_INF;
00348
00349 if (val>0)
00350 {
00351
00352 hpb[num_hp] = Bi[i];
00353 hpf[num_hp] = val;
00354 num_hp++;
00355 }
00356
00357 if( (Bi[i] < 0 && val > 0) || (Bi[i] > 0 && val <= 0))
00358 GradVal += Bi[i];
00359
00360 }
00361
00362 t = 0;
00363 if( GradVal < 0 )
00364 {
00365 start_time = get_time();
00366
00367 if( sort(hpf, hpb, num_hp) != 0 )
00368 {
00369 ocas.exitflag=-2;
00370 goto cleanup;
00371 }
00372 ocas.sort_time += get_time() - start_time;
00373
00374 float64_t t_new, GradVal_new;
00375 i = 0;
00376 while( GradVal < 0 && i < num_hp )
00377 {
00378 t_new = hpf[i];
00379 GradVal_new = GradVal + LIBOCAS_ABS(hpb[i]) + A0*(t_new-t);
00380
00381 if( GradVal_new >= 0 )
00382 {
00383 t = t + GradVal*(t-t_new)/(GradVal_new - GradVal);
00384 }
00385 else
00386 {
00387 t = t_new;
00388 i++;
00389 }
00390
00391 GradVal = GradVal_new;
00392 }
00393 }
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406 t = LIBOCAS_MAX(t,0);
00407
00408 t1 = t;
00409 t2 = t+LAMBDA*(1.0-t);
00410
00411
00412
00413 sq_norm_W = update_W( t1, user_data );
00414
00415
00416 xi = 0;
00417 cut_length = 0;
00418 ocas.trn_err = 0;
00419 for(i=0; i < nData; i++ ) {
00420
00421 if( (old_output[i]*(1-t2) + t2*output[i]) <= 1 )
00422 {
00423 new_cut[cut_length] = i;
00424 cut_length++;
00425 }
00426
00427 output[i] = old_output[i]*(1-t1) + t1*output[i];
00428
00429 if( output[i] <= 1) xi += 1-output[i];
00430 if( output[i] <= 0) ocas.trn_err++;
00431
00432 }
00433
00434 ocas.Q_P = 0.5*sq_norm_W + C*xi;
00435
00436 ocas.ocas_time = get_time() - ocas_start_time;
00437
00438
00439
00440
00441
00442
00443 start_time = get_time();
00444 ocas_print(ocas);
00445 ocas.print_time += get_time() - start_time;
00446
00447 break;
00448 }
00449
00450
00451 if( ocas.Q_P - ocas.Q_D <= TolRel*LIBOCAS_ABS(ocas.Q_P)) ocas.exitflag = 1;
00452 if( ocas.Q_P - ocas.Q_D <= TolAbs) ocas.exitflag = 2;
00453 if( ocas.Q_P <= QPBound) ocas.exitflag = 3;
00454 if( ocas.ocas_time >= MaxTime) ocas.exitflag = 4;
00455 if(ocas.nCutPlanes >= BufSize) ocas.exitflag = -1;
00456
00457 }
00458
00459 cleanup:
00460
00461 LIBOCAS_FREE(H);
00462 LIBOCAS_FREE(b);
00463 LIBOCAS_FREE(alpha);
00464 LIBOCAS_FREE(new_cut);
00465 LIBOCAS_FREE(I);
00466 LIBOCAS_FREE(diag_H);
00467 LIBOCAS_FREE(output);
00468 LIBOCAS_FREE(old_output);
00469 LIBOCAS_FREE(hpf);
00470
00471 LIBOCAS_FREE(hpb);
00472 LIBOCAS_FREE(Ci);
00473 LIBOCAS_FREE(Bi);
00474
00475 ocas.ocas_time = get_time() - ocas_start_time;
00476
00477 return(ocas);
00478 }
00479
00480
00481
00482
00483
00484
00485 ocas_return_value_T svm_ocas_solver_difC(
00486 float64_t *C,
00487 uint32_t nData,
00488 float64_t TolRel,
00489 float64_t TolAbs,
00490 float64_t QPBound,
00491 float64_t MaxTime,
00492 uint32_t _BufSize,
00493 uint8_t Method,
00494 void (*compute_W)(float64_t*, float64_t*, float64_t*, uint32_t, void*),
00495 float64_t (*update_W)(float64_t, void*),
00496 int (*add_new_cut)(float64_t*, uint32_t*, uint32_t, uint32_t, void*),
00497 int (*compute_output)(float64_t*, void* ),
00498 int (*sort)(float64_t*, float64_t*, uint32_t),
00499 void (*ocas_print)(ocas_return_value_T),
00500 void* user_data)
00501 {
00502 ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
00503 float64_t *b, *alpha, *diag_H;
00504 float64_t *output, *old_output;
00505 float64_t xi, sq_norm_W, QPSolverTolRel, dot_prod_WoldW, sq_norm_oldW;
00506 float64_t A0, B0, GradVal, t, t1, t2, *Ci, *Bi, *hpf, *hpb;
00507 float64_t start_time, ocas_start_time;
00508 float64_t qp_b = 1.0;
00509 float64_t new_b;
00510 uint32_t cut_length;
00511 uint32_t i, *new_cut;
00512 uint32_t *I;
00513 uint8_t S = 1;
00514 libqp_state_T qp_exitflag;
00515
00516 ocas_start_time = get_time();
00517 ocas.qp_solver_time = 0;
00518 ocas.output_time = 0;
00519 ocas.sort_time = 0;
00520 ocas.add_time = 0;
00521 ocas.w_time = 0;
00522 ocas.print_time = 0;
00523
00524 BufSize = _BufSize;
00525
00526 QPSolverTolRel = TolRel*0.5;
00527
00528 H=NULL;
00529 b=NULL;
00530 alpha=NULL;
00531 new_cut=NULL;
00532 I=NULL;
00533 diag_H=NULL;
00534 output=NULL;
00535 old_output=NULL;
00536 hpf=NULL;
00537 hpb = NULL;
00538 Ci=NULL;
00539 Bi=NULL;
00540
00541
00542 H = (float64_t*)LIBOCAS_CALLOC(BufSize*BufSize,sizeof(float64_t));
00543 if(H == NULL)
00544 {
00545 ocas.exitflag=-2;
00546 goto cleanup;
00547 }
00548
00549
00550 b = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
00551 if(b == NULL)
00552 {
00553 ocas.exitflag=-2;
00554 goto cleanup;
00555 }
00556
00557 alpha = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
00558 if(alpha == NULL)
00559 {
00560 ocas.exitflag=-2;
00561 goto cleanup;
00562 }
00563
00564
00565 new_cut = (uint32_t*)LIBOCAS_CALLOC(nData,sizeof(uint32_t));
00566 if(new_cut == NULL)
00567 {
00568 ocas.exitflag=-2;
00569 goto cleanup;
00570 }
00571
00572 I = (uint32_t*)LIBOCAS_CALLOC(BufSize,sizeof(uint32_t));
00573 if(I == NULL)
00574 {
00575 ocas.exitflag=-2;
00576 goto cleanup;
00577 }
00578
00579 for(i=0; i< BufSize; i++) I[i] = 1;
00580
00581 diag_H = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
00582 if(diag_H == NULL)
00583 {
00584 ocas.exitflag=-2;
00585 goto cleanup;
00586 }
00587
00588 output = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
00589 if(output == NULL)
00590 {
00591 ocas.exitflag=-2;
00592 goto cleanup;
00593 }
00594
00595 old_output = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
00596 if(old_output == NULL)
00597 {
00598 ocas.exitflag=-2;
00599 goto cleanup;
00600 }
00601
00602
00603 hpf = (float64_t*) LIBOCAS_CALLOC(nData, sizeof(hpf[0]));
00604 if(hpf == NULL)
00605 {
00606 ocas.exitflag=-2;
00607 goto cleanup;
00608 }
00609
00610 hpb = (float64_t*) LIBOCAS_CALLOC(nData, sizeof(hpb[0]));
00611 if(hpb == NULL)
00612 {
00613 ocas.exitflag=-2;
00614 goto cleanup;
00615 }
00616
00617
00618 Ci = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
00619 if(Ci == NULL)
00620 {
00621 ocas.exitflag=-2;
00622 goto cleanup;
00623 }
00624
00625 Bi = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
00626 if(Bi == NULL)
00627 {
00628 ocas.exitflag=-2;
00629 goto cleanup;
00630 }
00631
00632 ocas.nCutPlanes = 0;
00633 ocas.exitflag = 0;
00634 ocas.nIter = 0;
00635
00636
00637 sq_norm_W = 0;
00638 xi = nData;
00639
00640 ocas.Q_D = 0;
00641
00642
00643 cut_length = nData;
00644 new_b = 0;
00645 for(i=0; i < nData; i++)
00646 {
00647 new_cut[i] = i;
00648 new_b += C[i];
00649 }
00650
00651 ocas.Q_P = 0.5*sq_norm_W + new_b;
00652
00653
00654 ocas.trn_err = nData;
00655 ocas.ocas_time = get_time() - ocas_start_time;
00656
00657
00658
00659 ocas_print(ocas);
00660
00661
00662 while( ocas.exitflag == 0 )
00663 {
00664 ocas.nIter++;
00665
00666
00667
00668 b[ocas.nCutPlanes] = -new_b;
00669
00670 start_time = get_time();
00671
00672 if(add_new_cut( &H[LIBOCAS_INDEX(0,ocas.nCutPlanes,BufSize)], new_cut, cut_length, ocas.nCutPlanes, user_data ) != 0)
00673 {
00674 ocas.exitflag=-2;
00675 goto cleanup;
00676 }
00677
00678 ocas.add_time += get_time() - start_time;
00679
00680
00681 diag_H[ocas.nCutPlanes] = H[LIBOCAS_INDEX(ocas.nCutPlanes,ocas.nCutPlanes,BufSize)];
00682 for(i=0; i < ocas.nCutPlanes; i++) {
00683 H[LIBOCAS_INDEX(ocas.nCutPlanes,i,BufSize)] = H[LIBOCAS_INDEX(i,ocas.nCutPlanes,BufSize)];
00684 }
00685
00686 ocas.nCutPlanes++;
00687
00688
00689 start_time = get_time();
00690
00691
00692
00693 qp_exitflag = libqp_splx_solver(&get_col, diag_H, b, &qp_b, I, &S, alpha,
00694 ocas.nCutPlanes, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBOCAS_PLUS_INF,0);
00695
00696 ocas.qp_exitflag = qp_exitflag.exitflag;
00697
00698 ocas.qp_solver_time += get_time() - start_time;
00699 ocas.Q_D = -qp_exitflag.QP;
00700
00701 ocas.nNZAlpha = 0;
00702 for(i=0; i < ocas.nCutPlanes; i++) {
00703 if( alpha[i] != 0) ocas.nNZAlpha++;
00704 }
00705
00706 sq_norm_oldW = sq_norm_W;
00707 start_time = get_time();
00708 compute_W( &sq_norm_W, &dot_prod_WoldW, alpha, ocas.nCutPlanes, user_data );
00709 ocas.w_time += get_time() - start_time;
00710
00711
00712 switch( Method )
00713 {
00714
00715 case 0:
00716
00717 start_time = get_time();
00718 if( compute_output( output, user_data ) != 0)
00719 {
00720 ocas.exitflag=-2;
00721 goto cleanup;
00722 }
00723 ocas.output_time += get_time()-start_time;
00724
00725 xi = 0;
00726 cut_length = 0;
00727 ocas.trn_err = 0;
00728 new_b = 0;
00729 for(i=0; i < nData; i++)
00730 {
00731 if(output[i] <= 0) ocas.trn_err++;
00732
00733
00734
00735 if(output[i] <= C[i]) {
00736 xi += C[i] - output[i];
00737 new_cut[cut_length] = i;
00738 cut_length++;
00739 new_b += C[i];
00740 }
00741 }
00742
00743 ocas.Q_P = 0.5*sq_norm_W + xi;
00744
00745 ocas.ocas_time = get_time() - ocas_start_time;
00746
00747
00748
00749
00750
00751
00752 start_time = get_time();
00753 ocas_print(ocas);
00754 ocas.print_time += get_time() - start_time;
00755
00756 break;
00757
00758
00759
00760 case 1:
00761
00762
00763 A0 = sq_norm_W -2*dot_prod_WoldW + sq_norm_oldW;
00764 B0 = dot_prod_WoldW - sq_norm_oldW;
00765
00766 memcpy( old_output, output, sizeof(float64_t)*nData );
00767
00768 start_time = get_time();
00769 if( compute_output( output, user_data ) != 0)
00770 {
00771 ocas.exitflag=-2;
00772 goto cleanup;
00773 }
00774 ocas.output_time += get_time()-start_time;
00775
00776 uint32_t num_hp = 0;
00777 GradVal = B0;
00778 for(i=0; i< nData; i++) {
00779
00780
00781
00782 Ci[i] = (C[i]-old_output[i]);
00783 Bi[i] = old_output[i] - output[i];
00784
00785 float64_t val;
00786 if(Bi[i] != 0)
00787 val = -Ci[i]/Bi[i];
00788 else
00789 val = -LIBOCAS_PLUS_INF;
00790
00791 if (val>0)
00792 {
00793
00794 hpb[num_hp] = Bi[i];
00795 hpf[num_hp] = val;
00796 num_hp++;
00797 }
00798
00799 if( (Bi[i] < 0 && val > 0) || (Bi[i] > 0 && val <= 0))
00800 GradVal += Bi[i];
00801
00802 }
00803
00804 t = 0;
00805 if( GradVal < 0 )
00806 {
00807 start_time = get_time();
00808
00809 if( sort(hpf, hpb, num_hp) != 0 )
00810 {
00811 ocas.exitflag=-2;
00812 goto cleanup;
00813 }
00814 ocas.sort_time += get_time() - start_time;
00815
00816 float64_t t_new, GradVal_new;
00817 i = 0;
00818 while( GradVal < 0 && i < num_hp )
00819 {
00820 t_new = hpf[i];
00821 GradVal_new = GradVal + LIBOCAS_ABS(hpb[i]) + A0*(t_new-t);
00822
00823 if( GradVal_new >= 0 )
00824 {
00825 t = t + GradVal*(t-t_new)/(GradVal_new - GradVal);
00826 }
00827 else
00828 {
00829 t = t_new;
00830 i++;
00831 }
00832
00833 GradVal = GradVal_new;
00834 }
00835 }
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848 t = LIBOCAS_MAX(t,0);
00849
00850 t1 = t;
00851 t2 = t+(1.0-t)*LAMBDA;
00852
00853
00854
00855 sq_norm_W = update_W( t1, user_data );
00856
00857
00858 xi = 0;
00859 cut_length = 0;
00860 ocas.trn_err = 0;
00861 new_b = 0;
00862 for(i=0; i < nData; i++ ) {
00863
00864
00865 if( (old_output[i]*(1-t2) + t2*output[i]) <= C[i] )
00866 {
00867 new_cut[cut_length] = i;
00868 cut_length++;
00869 new_b += C[i];
00870 }
00871
00872 output[i] = old_output[i]*(1-t1) + t1*output[i];
00873
00874
00875 if( output[i] <= C[i]) xi += C[i]-output[i];
00876 if( output[i] <= 0) ocas.trn_err++;
00877
00878 }
00879
00880
00881 ocas.Q_P = 0.5*sq_norm_W + xi;
00882
00883 ocas.ocas_time = get_time() - ocas_start_time;
00884
00885
00886
00887
00888
00889
00890 start_time = get_time();
00891 ocas_print(ocas);
00892 ocas.print_time += get_time() - start_time;
00893
00894 break;
00895 }
00896
00897
00898 if( ocas.Q_P - ocas.Q_D <= TolRel*LIBOCAS_ABS(ocas.Q_P)) ocas.exitflag = 1;
00899 if( ocas.Q_P - ocas.Q_D <= TolAbs) ocas.exitflag = 2;
00900 if( ocas.Q_P <= QPBound) ocas.exitflag = 3;
00901 if( ocas.ocas_time >= MaxTime) ocas.exitflag = 4;
00902 if(ocas.nCutPlanes >= BufSize) ocas.exitflag = -1;
00903
00904 }
00905
00906 cleanup:
00907
00908 LIBOCAS_FREE(H);
00909 LIBOCAS_FREE(b);
00910 LIBOCAS_FREE(alpha);
00911 LIBOCAS_FREE(new_cut);
00912 LIBOCAS_FREE(I);
00913 LIBOCAS_FREE(diag_H);
00914 LIBOCAS_FREE(output);
00915 LIBOCAS_FREE(old_output);
00916 LIBOCAS_FREE(hpf);
00917
00918 LIBOCAS_FREE(hpb);
00919 LIBOCAS_FREE(Ci);
00920 LIBOCAS_FREE(Bi);
00921
00922 ocas.ocas_time = get_time() - ocas_start_time;
00923
00924 return(ocas);
00925 }
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937 static void findactive(float64_t *Theta, float64_t *SortedA, uint32_t *nSortedA, float64_t *A, float64_t *B, int n,
00938 int (*sort)(float64_t*, float64_t*, uint32_t))
00939 {
00940 float64_t tmp, theta;
00941 int32_t i, j, idx, idx2 = 0, start;
00942
00943 sort(A,B,n);
00944
00945 tmp = B[0];
00946 idx = 0;
00947 i = 0;
00948 while( i < n-1 && A[i] == A[i+1])
00949 {
00950 if( B[i+1] > B[idx] )
00951 {
00952 idx = i+1;
00953 tmp = B[i+1];
00954 }
00955 i++;
00956 }
00957
00958 (*nSortedA) = 1;
00959 SortedA[0] = A[idx];
00960
00961 while(1)
00962 {
00963 start = idx + 1;
00964 while( start < n && A[idx] == A[start])
00965 start++;
00966
00967 theta = LIBOCAS_PLUS_INF;
00968 for(j=start; j < n; j++)
00969 {
00970 tmp = (B[j] - B[idx])/(A[idx]-A[j]);
00971 if( tmp < theta)
00972 {
00973 theta = tmp;
00974 idx2 = j;
00975 }
00976 }
00977
00978 if( theta < LIBOCAS_PLUS_INF)
00979 {
00980 Theta[(*nSortedA) - 1] = theta;
00981 SortedA[(*nSortedA)] = A[idx2];
00982 (*nSortedA)++;
00983 idx = idx2;
00984 }
00985 else
00986 return;
00987 }
00988 }
00989
00990
00991
00992
00993
00994 ocas_return_value_T msvm_ocas_solver(
00995 float64_t C,
00996 float64_t *data_y,
00997 uint32_t nY,
00998 uint32_t nData,
00999 float64_t TolRel,
01000 float64_t TolAbs,
01001 float64_t QPBound,
01002 float64_t MaxTime,
01003 uint32_t _BufSize,
01004 uint8_t Method,
01005 void (*compute_W)(float64_t*, float64_t*, float64_t*, uint32_t, void*),
01006 float64_t (*update_W)(float64_t, void*),
01007 int (*add_new_cut)(float64_t*, uint32_t*, uint32_t, void*),
01008 int (*compute_output)(float64_t*, void* ),
01009 int (*sort)(float64_t*, float64_t*, uint32_t),
01010 void (*ocas_print)(ocas_return_value_T),
01011 void* user_data)
01012 {
01013 ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
01014 float64_t *b, *alpha, *diag_H;
01015 float64_t *output, *old_output;
01016 float64_t xi, sq_norm_W, QPSolverTolRel, QPSolverTolAbs, dot_prod_WoldW, sq_norm_oldW;
01017 float64_t A0, B0, t, t1, t2, R, tmp, element_b, x;
01018 float64_t *A, *B, *theta, *Theta, *sortedA, *Add;
01019 float64_t start_time, ocas_start_time, grad_sum, grad, min_x = 0, old_x, old_grad;
01020 uint32_t i, y, y2, ypred = 0, *new_cut, cnt1, cnt2, j, nSortedA, idx;
01021 uint32_t *I;
01022 uint8_t S = 1;
01023 libqp_state_T qp_exitflag;
01024
01025 ocas_start_time = get_time();
01026 ocas.qp_solver_time = 0;
01027 ocas.output_time = 0;
01028 ocas.sort_time = 0;
01029 ocas.add_time = 0;
01030 ocas.w_time = 0;
01031 ocas.print_time = 0;
01032
01033 BufSize = _BufSize;
01034
01035 QPSolverTolRel = TolRel*0.5;
01036 QPSolverTolAbs = TolAbs*0.5;
01037
01038 H=NULL;
01039 b=NULL;
01040 alpha=NULL;
01041 new_cut=NULL;
01042 I=NULL;
01043 diag_H=NULL;
01044 output=NULL;
01045 old_output=NULL;
01046 A = NULL;
01047 B = NULL;
01048 theta = NULL;
01049 Theta = NULL;
01050 sortedA = NULL;
01051 Add = NULL;
01052
01053
01054 H = (float64_t*)LIBOCAS_CALLOC(BufSize*BufSize,sizeof(float64_t));
01055 if(H == NULL)
01056 {
01057 ocas.exitflag=-2;
01058 goto cleanup;
01059 }
01060
01061
01062 b = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
01063 if(b == NULL)
01064 {
01065 ocas.exitflag=-2;
01066 goto cleanup;
01067 }
01068
01069 alpha = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
01070 if(alpha == NULL)
01071 {
01072 ocas.exitflag=-2;
01073 goto cleanup;
01074 }
01075
01076
01077 new_cut = (uint32_t*)LIBOCAS_CALLOC(nData,sizeof(uint32_t));
01078 if(new_cut == NULL)
01079 {
01080 ocas.exitflag=-2;
01081 goto cleanup;
01082 }
01083
01084 I = (uint32_t*)LIBOCAS_CALLOC(BufSize,sizeof(uint32_t));
01085 if(I == NULL)
01086 {
01087 ocas.exitflag=-2;
01088 goto cleanup;
01089 }
01090
01091 for(i=0; i< BufSize; i++)
01092 I[i] = 1;
01093
01094 diag_H = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
01095 if(diag_H == NULL)
01096 {
01097 ocas.exitflag=-2;
01098 goto cleanup;
01099 }
01100
01101 output = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t));
01102 if(output == NULL)
01103 {
01104 ocas.exitflag=-2;
01105 goto cleanup;
01106 }
01107
01108 old_output = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t));
01109 if(old_output == NULL)
01110 {
01111 ocas.exitflag=-2;
01112 goto cleanup;
01113 }
01114
01115
01116 A = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t));
01117 if(A == NULL)
01118 {
01119 ocas.exitflag=-2;
01120 goto cleanup;
01121 }
01122
01123 B = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t));
01124 if(B == NULL)
01125 {
01126 ocas.exitflag=-2;
01127 goto cleanup;
01128 }
01129
01130 theta = (float64_t*)LIBOCAS_CALLOC(nY,sizeof(float64_t));
01131 if(theta == NULL)
01132 {
01133 ocas.exitflag=-2;
01134 goto cleanup;
01135 }
01136
01137 sortedA = (float64_t*)LIBOCAS_CALLOC(nY,sizeof(float64_t));
01138 if(sortedA == NULL)
01139 {
01140 ocas.exitflag=-2;
01141 goto cleanup;
01142 }
01143
01144 Theta = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t));
01145 if(Theta == NULL)
01146 {
01147 ocas.exitflag=-2;
01148 goto cleanup;
01149 }
01150
01151 Add = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t));
01152 if(Add == NULL)
01153 {
01154 ocas.exitflag=-2;
01155 goto cleanup;
01156 }
01157
01158
01159 ocas.nCutPlanes = 0;
01160 ocas.exitflag = 0;
01161 ocas.nIter = 0;
01162 ocas.Q_D = 0;
01163 ocas.trn_err = nData;
01164 R = (float64_t)nData;
01165 sq_norm_W = 0;
01166 element_b = (float64_t)nData;
01167 ocas.Q_P = 0.5*sq_norm_W + C*R;
01168
01169
01170 for(i=0; i < nData; i++)
01171 {
01172 y2 = (uint32_t)data_y[i]-1;
01173
01174 if(y2 > 0)
01175 new_cut[i] = 0;
01176 else
01177 new_cut[i] = 1;
01178
01179 }
01180
01181 ocas.ocas_time = get_time() - ocas_start_time;
01182
01183 start_time = get_time();
01184 ocas_print(ocas);
01185 ocas.print_time += get_time() - start_time;
01186
01187
01188 while( ocas.exitflag == 0 )
01189 {
01190 ocas.nIter++;
01191
01192
01193 b[ocas.nCutPlanes] = -(float64_t)element_b;
01194
01195 start_time = get_time();
01196
01197 if(add_new_cut( &H[LIBOCAS_INDEX(0,ocas.nCutPlanes,BufSize)], new_cut, ocas.nCutPlanes, user_data ) != 0)
01198 {
01199 ocas.exitflag=-2;
01200 goto cleanup;
01201 }
01202
01203 ocas.add_time += get_time() - start_time;
01204
01205
01206 diag_H[ocas.nCutPlanes] = H[LIBOCAS_INDEX(ocas.nCutPlanes,ocas.nCutPlanes,BufSize)];
01207 for(i=0; i < ocas.nCutPlanes; i++)
01208 {
01209 H[LIBOCAS_INDEX(ocas.nCutPlanes,i,BufSize)] = H[LIBOCAS_INDEX(i,ocas.nCutPlanes,BufSize)];
01210 }
01211
01212 ocas.nCutPlanes++;
01213
01214
01215 start_time = get_time();
01216
01217 qp_exitflag = libqp_splx_solver(&get_col, diag_H, b, &C, I, &S, alpha,
01218 ocas.nCutPlanes, QPSolverMaxIter, QPSolverTolAbs, QPSolverTolRel, -LIBOCAS_PLUS_INF,0);
01219
01220 ocas.qp_exitflag = qp_exitflag.exitflag;
01221
01222 ocas.qp_solver_time += get_time() - start_time;
01223 ocas.Q_D = -qp_exitflag.QP;
01224
01225 ocas.nNZAlpha = 0;
01226 for(i=0; i < ocas.nCutPlanes; i++)
01227 if( alpha[i] != 0) ocas.nNZAlpha++;
01228
01229 sq_norm_oldW = sq_norm_W;
01230 start_time = get_time();
01231 compute_W( &sq_norm_W, &dot_prod_WoldW, alpha, ocas.nCutPlanes, user_data );
01232 ocas.w_time += get_time() - start_time;
01233
01234
01235 switch( Method )
01236 {
01237
01238 case 0:
01239
01240 start_time = get_time();
01241 if( compute_output( output, user_data ) != 0)
01242 {
01243 ocas.exitflag=-2;
01244 goto cleanup;
01245 }
01246 ocas.output_time += get_time()-start_time;
01247
01248
01249 element_b = 0.0;
01250 R = 0;
01251 ocas.trn_err = 0;
01252
01253 for(i=0; i < nData; i++)
01254 {
01255 y2 = (uint32_t)data_y[i]-1;
01256
01257 for(xi=-LIBOCAS_PLUS_INF, y=0; y < nY; y++)
01258 {
01259 if(y2 != y && xi < output[LIBOCAS_INDEX(y,i,nY)])
01260 {
01261 xi = output[LIBOCAS_INDEX(y,i,nY)];
01262 ypred = y;
01263 }
01264 }
01265
01266 if(xi >= output[LIBOCAS_INDEX(y2,i,nY)])
01267 ocas.trn_err ++;
01268
01269 xi = LIBOCAS_MAX(0,xi+1-output[LIBOCAS_INDEX(y2,i,nY)]);
01270 R += xi;
01271 if(xi > 0)
01272 {
01273 element_b++;
01274 new_cut[i] = ypred;
01275 }
01276 else
01277 new_cut[i] = y2;
01278 }
01279
01280 ocas.Q_P = 0.5*sq_norm_W + C*R;
01281
01282 ocas.ocas_time = get_time() - ocas_start_time;
01283
01284 start_time = get_time();
01285 ocas_print(ocas);
01286 ocas.print_time += get_time() - start_time;
01287
01288 break;
01289
01290
01291 case 1:
01292 memcpy( old_output, output, sizeof(float64_t)*nData*nY );
01293
01294 start_time = get_time();
01295 if( compute_output( output, user_data ) != 0)
01296 {
01297 ocas.exitflag=-2;
01298 goto cleanup;
01299 }
01300 ocas.output_time += get_time()-start_time;
01301
01302 A0 = sq_norm_W - 2*dot_prod_WoldW + sq_norm_oldW;
01303 B0 = dot_prod_WoldW - sq_norm_oldW;
01304
01305 for(i=0; i < nData; i++)
01306 {
01307 y2 = (uint32_t)data_y[i]-1;
01308
01309 for(y=0; y < nY; y++)
01310 {
01311 A[LIBOCAS_INDEX(y,i,nY)] = C*(output[LIBOCAS_INDEX(y,i,nY)] - old_output[LIBOCAS_INDEX(y,i,nY)]
01312 + old_output[LIBOCAS_INDEX(y2,i,nY)] - output[LIBOCAS_INDEX(y2,i,nY)]);
01313 B[LIBOCAS_INDEX(y,i,nY)] = C*(old_output[LIBOCAS_INDEX(y,i,nY)] - old_output[LIBOCAS_INDEX(y2,i,nY)]
01314 + (float64_t)(y != y2));
01315 }
01316 }
01317
01318
01319
01320
01321 grad_sum = B0;
01322 cnt1 = 0;
01323 cnt2 = 0;
01324 for(i=0; i < nData; i++)
01325 {
01326 findactive(theta,sortedA,&nSortedA,&A[i*nY],&B[i*nY],nY,sort);
01327
01328 idx = 0;
01329 while( idx < nSortedA-1 && theta[idx] < 0 )
01330 idx++;
01331
01332 grad_sum += sortedA[idx];
01333
01334 for(j=idx; j < nSortedA-1; j++)
01335 {
01336 Theta[cnt1] = theta[j];
01337 cnt1++;
01338 }
01339
01340 for(j=idx+1; j < nSortedA; j++)
01341 {
01342 Add[cnt2] = -sortedA[j-1]+sortedA[j];
01343 cnt2++;
01344 }
01345 }
01346
01347 start_time = get_time();
01348 sort(Theta,Add,cnt1);
01349 ocas.sort_time += get_time() - start_time;
01350
01351 grad = grad_sum;
01352 if(grad >= 0)
01353 {
01354 min_x = 0;
01355 }
01356 else
01357 {
01358 old_x = 0;
01359 old_grad = grad;
01360
01361 for(i=0; i < cnt1; i++)
01362 {
01363 x = Theta[i];
01364
01365 grad = x*A0 + grad_sum;
01366
01367 if(grad >=0)
01368 {
01369
01370 min_x = (grad*old_x - old_grad*x)/(grad - old_grad);
01371
01372 break;
01373 }
01374 else
01375 {
01376 grad_sum = grad_sum + Add[i];
01377
01378 grad = x*A0 + grad_sum;
01379 if( grad >= 0)
01380 {
01381 min_x = x;
01382 break;
01383 }
01384 }
01385
01386 old_grad = grad;
01387 old_x = x;
01388 }
01389 }
01390
01391
01392 t = min_x;
01393 t1 = t;
01394 t2 = t+(1.0-t)*LAMBDA;
01395
01396
01397
01398 sq_norm_W = update_W( t1, user_data );
01399
01400
01401 element_b = 0.0;
01402
01403 for(i=0; i < nData; i++)
01404 {
01405 y2 = (uint32_t)data_y[i]-1;
01406
01407 for(xi=-LIBOCAS_PLUS_INF, y=0; y < nY; y++)
01408 {
01409 tmp = old_output[LIBOCAS_INDEX(y,i,nY)]*(1-t2) + t2*output[LIBOCAS_INDEX(y,i,nY)];
01410 if(y2 != y && xi < tmp)
01411 {
01412 xi = tmp;
01413 ypred = y;
01414 }
01415 }
01416
01417 tmp = old_output[LIBOCAS_INDEX(y2,i,nY)]*(1-t2) + t2*output[LIBOCAS_INDEX(y2,i,nY)];
01418 xi = LIBOCAS_MAX(0,xi+1-tmp);
01419 if(xi > 0)
01420 {
01421 element_b++;
01422 new_cut[i] = ypred;
01423 }
01424 else
01425 new_cut[i] = y2;
01426 }
01427
01428
01429 ocas.trn_err = 0;
01430 R = 0;
01431 for(i=0; i < nData; i++)
01432 {
01433 y2 = (uint32_t)data_y[i]-1;
01434
01435 for(tmp=-LIBOCAS_PLUS_INF, y=0; y < nY; y++)
01436 {
01437 output[LIBOCAS_INDEX(y,i,nY)] = old_output[LIBOCAS_INDEX(y,i,nY)]*(1-t1) + t1*output[LIBOCAS_INDEX(y,i,nY)];
01438
01439 if(y2 != y && tmp < output[LIBOCAS_INDEX(y,i,nY)])
01440 {
01441 ypred = y;
01442 tmp = output[LIBOCAS_INDEX(y,i,nY)];
01443 }
01444 }
01445
01446 R += LIBOCAS_MAX(0,1+tmp - output[LIBOCAS_INDEX(y2,i,nY)]);
01447 if( tmp >= output[LIBOCAS_INDEX(y2,i,nY)])
01448 ocas.trn_err ++;
01449 }
01450
01451 ocas.Q_P = 0.5*sq_norm_W + C*R;
01452
01453
01454
01455 ocas.ocas_time = get_time() - ocas_start_time;
01456
01457 start_time = get_time();
01458 ocas_print(ocas);
01459 ocas.print_time += get_time() - start_time;
01460
01461 break;
01462
01463 }
01464
01465
01466 if( ocas.Q_P - ocas.Q_D <= TolRel*LIBOCAS_ABS(ocas.Q_P)) ocas.exitflag = 1;
01467 if( ocas.Q_P - ocas.Q_D <= TolAbs) ocas.exitflag = 2;
01468 if( ocas.Q_P <= QPBound) ocas.exitflag = 3;
01469 if( ocas.ocas_time >= MaxTime) ocas.exitflag = 4;
01470 if(ocas.nCutPlanes >= BufSize) ocas.exitflag = -1;
01471
01472 }
01473
01474 cleanup:
01475
01476 LIBOCAS_FREE(H);
01477 LIBOCAS_FREE(b);
01478 LIBOCAS_FREE(alpha);
01479 LIBOCAS_FREE(new_cut);
01480 LIBOCAS_FREE(I);
01481 LIBOCAS_FREE(diag_H);
01482 LIBOCAS_FREE(output);
01483 LIBOCAS_FREE(old_output);
01484 LIBOCAS_FREE(A);
01485 LIBOCAS_FREE(B);
01486 LIBOCAS_FREE(theta);
01487 LIBOCAS_FREE(Theta);
01488 LIBOCAS_FREE(sortedA);
01489 LIBOCAS_FREE(Add);
01490
01491 ocas.ocas_time = get_time() - ocas_start_time;
01492
01493 return(ocas);
01494 }
01495 }