49 if (gettimeofday(&tv, NULL)==0)
50 return tv.tv_sec+((
float64_t)(tv.tv_usec))/1e6;
69 int (*add_new_cut)(
float64_t*, uint32_t*, uint32_t, uint32_t,
void*),
70 int (*compute_output)(
float64_t*,
void* ),
72 void (*ocas_print)(ocas_return_value_T),
75 ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
78 float64_t xi, sq_norm_W, QPSolverTolRel, dot_prod_WoldW, sq_norm_oldW;
79 float64_t A0, B0, GradVal, t, t1, t2, *Ci, *Bi, *hpf, *hpb;
85 libqp_state_T qp_exitflag;
88 ocas.qp_solver_time = 0;
98 QPSolverTolRel = TolRel*0.5;
137 new_cut = (uint32_t*)LIBOCAS_CALLOC(nData,
sizeof(uint32_t));
144 I = (uint32_t*)LIBOCAS_CALLOC(BufSize,
sizeof(uint32_t));
151 for(i=0; i<
BufSize; i++) I[i] = 1;
168 if(old_output == NULL)
175 hpf = (
float64_t*) LIBOCAS_CALLOC(nData,
sizeof(hpf[0]));
182 hpb = (
float64_t*) LIBOCAS_CALLOC(nData,
sizeof(hpb[0]));
211 ocas.Q_P = 0.5*sq_norm_W + C*xi;
216 for(i=0; i < nData; i++)
222 ocas.trn_err = nData;
223 ocas.ocas_time =
get_time() - ocas_start_time;
230 while( ocas.exitflag == 0 )
235 b[ocas.nCutPlanes] = -(
float64_t)cut_length;
239 if(add_new_cut( &
H[LIBOCAS_INDEX(0,ocas.nCutPlanes,BufSize)], new_cut, cut_length, ocas.nCutPlanes, user_data ) != 0)
245 ocas.add_time +=
get_time() - start_time;
248 diag_H[ocas.nCutPlanes] =
H[LIBOCAS_INDEX(ocas.nCutPlanes,ocas.nCutPlanes,BufSize)];
249 for(i=0; i < ocas.nCutPlanes; i++) {
250 H[LIBOCAS_INDEX(ocas.nCutPlanes,i,BufSize)] =
H[LIBOCAS_INDEX(i,ocas.nCutPlanes,BufSize)];
259 ocas.nCutPlanes,
QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBOCAS_PLUS_INF,0);
261 ocas.qp_exitflag = qp_exitflag.exitflag;
263 ocas.qp_solver_time +=
get_time() - start_time;
264 ocas.Q_D = -qp_exitflag.QP;
267 for(i=0; i < ocas.nCutPlanes; i++) {
268 if( alpha[i] != 0) ocas.nNZAlpha++;
271 sq_norm_oldW = sq_norm_W;
273 compute_W( &sq_norm_W, &dot_prod_WoldW, alpha, ocas.nCutPlanes, user_data );
274 ocas.w_time +=
get_time() - start_time;
283 if( compute_output( output, user_data ) != 0)
288 ocas.output_time +=
get_time()-start_time;
293 for(i=0; i < nData; i++)
295 if(output[i] <= 0) ocas.trn_err++;
299 new_cut[cut_length] = i;
303 ocas.Q_P = 0.5*sq_norm_W + C*xi;
314 ocas.print_time +=
get_time() - start_time;
323 A0 = sq_norm_W -2*dot_prod_WoldW + sq_norm_oldW;
324 B0 = dot_prod_WoldW - sq_norm_oldW;
326 memcpy( old_output, output,
sizeof(
float64_t)*nData );
329 if( compute_output( output, user_data ) != 0)
334 ocas.output_time +=
get_time()-start_time;
338 for(i=0; i< nData; i++) {
340 Ci[i] = C*(1-old_output[i]);
341 Bi[i] = C*(old_output[i] - output[i]);
347 val = -LIBOCAS_PLUS_INF;
357 if( (Bi[i] < 0 && val > 0) || (Bi[i] > 0 && val <= 0))
367 if( sort(hpf, hpb, num_hp) != 0 )
372 ocas.sort_time +=
get_time() - start_time;
376 while( GradVal < 0 && i < num_hp )
379 GradVal_new = GradVal + LIBOCAS_ABS(hpb[i]) + A0*(t_new-t);
381 if( GradVal_new >= 0 )
383 t = t + GradVal*(t-t_new)/(GradVal_new - GradVal);
391 GradVal = GradVal_new;
406 t = LIBOCAS_MAX(t,0);
413 sq_norm_W = update_W( t1, user_data );
419 for(i=0; i < nData; i++ ) {
421 if( (old_output[i]*(1-t2) + t2*output[i]) <= 1 )
423 new_cut[cut_length] = i;
427 output[i] = old_output[i]*(1-t1) + t1*output[i];
429 if( output[i] <= 1) xi += 1-output[i];
430 if( output[i] <= 0) ocas.trn_err++;
434 ocas.Q_P = 0.5*sq_norm_W + C*xi;
436 ocas.ocas_time =
get_time() - ocas_start_time;
445 ocas.print_time +=
get_time() - start_time;
451 if( ocas.Q_P - ocas.Q_D <= TolRel*LIBOCAS_ABS(ocas.Q_P)) ocas.exitflag = 1;
452 if( ocas.Q_P - ocas.Q_D <= TolAbs) ocas.exitflag = 2;
453 if( ocas.Q_P <= QPBound) ocas.exitflag = 3;
454 if( ocas.ocas_time >= MaxTime) ocas.exitflag = 4;
455 if(ocas.nCutPlanes >= BufSize) ocas.exitflag = -1;
464 LIBOCAS_FREE(new_cut);
466 LIBOCAS_FREE(diag_H);
467 LIBOCAS_FREE(output);
468 LIBOCAS_FREE(old_output);
475 ocas.ocas_time =
get_time() - ocas_start_time;
496 int (*add_new_cut)(
float64_t*, uint32_t*, uint32_t, uint32_t,
void*),
497 int (*compute_output)(
float64_t*,
void* ),
499 void (*ocas_print)(ocas_return_value_T),
502 ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
505 float64_t xi, sq_norm_W, QPSolverTolRel, dot_prod_WoldW, sq_norm_oldW;
506 float64_t A0, B0, GradVal, t, t1, t2, *Ci, *Bi, *hpf, *hpb;
511 uint32_t i, *new_cut;
514 libqp_state_T qp_exitflag;
517 ocas.qp_solver_time = 0;
518 ocas.output_time = 0;
526 QPSolverTolRel = TolRel*0.5;
565 new_cut = (uint32_t*)LIBOCAS_CALLOC(nData,
sizeof(uint32_t));
572 I = (uint32_t*)LIBOCAS_CALLOC(BufSize,
sizeof(uint32_t));
579 for(i=0; i<
BufSize; i++) I[i] = 1;
596 if(old_output == NULL)
603 hpf = (
float64_t*) LIBOCAS_CALLOC(nData,
sizeof(hpf[0]));
610 hpb = (
float64_t*) LIBOCAS_CALLOC(nData,
sizeof(hpb[0]));
645 for(i=0; i < nData; i++)
651 ocas.Q_P = 0.5*sq_norm_W + new_b;
654 ocas.trn_err = nData;
655 ocas.ocas_time =
get_time() - ocas_start_time;
662 while( ocas.exitflag == 0 )
668 b[ocas.nCutPlanes] = -new_b;
672 if(add_new_cut( &
H[LIBOCAS_INDEX(0,ocas.nCutPlanes,BufSize)], new_cut, cut_length, ocas.nCutPlanes, user_data ) != 0)
678 ocas.add_time +=
get_time() - start_time;
681 diag_H[ocas.nCutPlanes] =
H[LIBOCAS_INDEX(ocas.nCutPlanes,ocas.nCutPlanes,BufSize)];
682 for(i=0; i < ocas.nCutPlanes; i++) {
683 H[LIBOCAS_INDEX(ocas.nCutPlanes,i,BufSize)] =
H[LIBOCAS_INDEX(i,ocas.nCutPlanes,BufSize)];
694 ocas.nCutPlanes,
QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBOCAS_PLUS_INF,0);
696 ocas.qp_exitflag = qp_exitflag.exitflag;
698 ocas.qp_solver_time +=
get_time() - start_time;
699 ocas.Q_D = -qp_exitflag.QP;
702 for(i=0; i < ocas.nCutPlanes; i++) {
703 if( alpha[i] != 0) ocas.nNZAlpha++;
706 sq_norm_oldW = sq_norm_W;
708 compute_W( &sq_norm_W, &dot_prod_WoldW, alpha, ocas.nCutPlanes, user_data );
709 ocas.w_time +=
get_time() - start_time;
718 if( compute_output( output, user_data ) != 0)
723 ocas.output_time +=
get_time()-start_time;
729 for(i=0; i < nData; i++)
731 if(output[i] <= 0) ocas.trn_err++;
735 if(output[i] <= C[i]) {
736 xi += C[i] - output[i];
737 new_cut[cut_length] = i;
743 ocas.Q_P = 0.5*sq_norm_W + xi;
745 ocas.ocas_time =
get_time() - ocas_start_time;
754 ocas.print_time +=
get_time() - start_time;
763 A0 = sq_norm_W -2*dot_prod_WoldW + sq_norm_oldW;
764 B0 = dot_prod_WoldW - sq_norm_oldW;
766 memcpy( old_output, output,
sizeof(
float64_t)*nData );
769 if( compute_output( output, user_data ) != 0)
774 ocas.output_time +=
get_time()-start_time;
778 for(i=0; i< nData; i++) {
782 Ci[i] = (C[i]-old_output[i]);
783 Bi[i] = old_output[i] - output[i];
789 val = -LIBOCAS_PLUS_INF;
799 if( (Bi[i] < 0 && val > 0) || (Bi[i] > 0 && val <= 0))
809 if( sort(hpf, hpb, num_hp) != 0 )
814 ocas.sort_time +=
get_time() - start_time;
818 while( GradVal < 0 && i < num_hp )
821 GradVal_new = GradVal + LIBOCAS_ABS(hpb[i]) + A0*(t_new-t);
823 if( GradVal_new >= 0 )
825 t = t + GradVal*(t-t_new)/(GradVal_new - GradVal);
833 GradVal = GradVal_new;
848 t = LIBOCAS_MAX(t,0);
855 sq_norm_W = update_W( t1, user_data );
862 for(i=0; i < nData; i++ ) {
865 if( (old_output[i]*(1-t2) + t2*output[i]) <= C[i] )
867 new_cut[cut_length] = i;
872 output[i] = old_output[i]*(1-t1) + t1*output[i];
875 if( output[i] <= C[i]) xi += C[i]-output[i];
876 if( output[i] <= 0) ocas.trn_err++;
881 ocas.Q_P = 0.5*sq_norm_W + xi;
883 ocas.ocas_time =
get_time() - ocas_start_time;
892 ocas.print_time +=
get_time() - start_time;
898 if( ocas.Q_P - ocas.Q_D <= TolRel*LIBOCAS_ABS(ocas.Q_P)) ocas.exitflag = 1;
899 if( ocas.Q_P - ocas.Q_D <= TolAbs) ocas.exitflag = 2;
900 if( ocas.Q_P <= QPBound) ocas.exitflag = 3;
901 if( ocas.ocas_time >= MaxTime) ocas.exitflag = 4;
902 if(ocas.nCutPlanes >= BufSize) ocas.exitflag = -1;
911 LIBOCAS_FREE(new_cut);
913 LIBOCAS_FREE(diag_H);
914 LIBOCAS_FREE(output);
915 LIBOCAS_FREE(old_output);
922 ocas.ocas_time =
get_time() - ocas_start_time;
941 int32_t i, j, idx, idx2 = 0, start;
948 while( i < n-1 && A[i] == A[i+1])
950 if( B[i+1] > B[idx] )
964 while( start < n && A[idx] == A[start])
967 theta = LIBOCAS_PLUS_INF;
968 for(j=start; j < n; j++)
970 tmp = (B[j] - B[idx])/(A[idx]-A[j]);
978 if( theta < LIBOCAS_PLUS_INF)
980 Theta[(*nSortedA) - 1] = theta;
981 SortedA[(*nSortedA)] = A[idx2];
1007 int (*add_new_cut)(
float64_t*, uint32_t*, uint32_t,
void*),
1008 int (*compute_output)(
float64_t*,
void* ),
1010 void (*ocas_print)(ocas_return_value_T),
1013 ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1016 float64_t xi, sq_norm_W, QPSolverTolRel, QPSolverTolAbs, dot_prod_WoldW, sq_norm_oldW;
1017 float64_t A0, B0, t, t1, t2, R, tmp, element_b, x;
1018 float64_t *A, *B, *theta, *Theta, *sortedA, *Add;
1019 float64_t start_time, ocas_start_time, grad_sum, grad, min_x = 0, old_x, old_grad;
1020 uint32_t i, y, y2, ypred = 0, *new_cut, cnt1, cnt2, j, nSortedA, idx;
1023 libqp_state_T qp_exitflag;
1026 ocas.qp_solver_time = 0;
1027 ocas.output_time = 0;
1031 ocas.print_time = 0;
1035 QPSolverTolRel = TolRel*0.5;
1036 QPSolverTolAbs = TolAbs*0.5;
1077 new_cut = (uint32_t*)LIBOCAS_CALLOC(nData,
sizeof(uint32_t));
1084 I = (uint32_t*)LIBOCAS_CALLOC(BufSize,
sizeof(uint32_t));
1109 if(old_output == NULL)
1159 ocas.nCutPlanes = 0;
1163 ocas.trn_err = nData;
1167 ocas.Q_P = 0.5*sq_norm_W + C*R;
1170 for(i=0; i < nData; i++)
1172 y2 = (uint32_t)data_y[i]-1;
1181 ocas.ocas_time =
get_time() - ocas_start_time;
1185 ocas.print_time +=
get_time() - start_time;
1188 while( ocas.exitflag == 0 )
1193 b[ocas.nCutPlanes] = -(
float64_t)element_b;
1197 if(add_new_cut( &
H[LIBOCAS_INDEX(0,ocas.nCutPlanes,BufSize)], new_cut, ocas.nCutPlanes, user_data ) != 0)
1203 ocas.add_time +=
get_time() - start_time;
1206 diag_H[ocas.nCutPlanes] =
H[LIBOCAS_INDEX(ocas.nCutPlanes,ocas.nCutPlanes,BufSize)];
1207 for(i=0; i < ocas.nCutPlanes; i++)
1209 H[LIBOCAS_INDEX(ocas.nCutPlanes,i,BufSize)] =
H[LIBOCAS_INDEX(i,ocas.nCutPlanes,BufSize)];
1218 ocas.nCutPlanes,
QPSolverMaxIter, QPSolverTolAbs, QPSolverTolRel, -LIBOCAS_PLUS_INF,0);
1220 ocas.qp_exitflag = qp_exitflag.exitflag;
1222 ocas.qp_solver_time +=
get_time() - start_time;
1223 ocas.Q_D = -qp_exitflag.QP;
1226 for(i=0; i < ocas.nCutPlanes; i++)
1227 if( alpha[i] != 0) ocas.nNZAlpha++;
1229 sq_norm_oldW = sq_norm_W;
1231 compute_W( &sq_norm_W, &dot_prod_WoldW, alpha, ocas.nCutPlanes, user_data );
1232 ocas.w_time +=
get_time() - start_time;
1241 if( compute_output( output, user_data ) != 0)
1246 ocas.output_time +=
get_time()-start_time;
1253 for(i=0; i < nData; i++)
1255 y2 = (uint32_t)data_y[i]-1;
1257 for(xi=-LIBOCAS_PLUS_INF, y=0; y < nY; y++)
1259 if(y2 != y && xi < output[LIBOCAS_INDEX(y,i,nY)])
1261 xi = output[LIBOCAS_INDEX(y,i,nY)];
1266 if(xi >= output[LIBOCAS_INDEX(y2,i,nY)])
1269 xi = LIBOCAS_MAX(0,xi+1-output[LIBOCAS_INDEX(y2,i,nY)]);
1280 ocas.Q_P = 0.5*sq_norm_W + C*R;
1282 ocas.ocas_time =
get_time() - ocas_start_time;
1286 ocas.print_time +=
get_time() - start_time;
1292 memcpy( old_output, output,
sizeof(
float64_t)*nData*nY );
1295 if( compute_output( output, user_data ) != 0)
1300 ocas.output_time +=
get_time()-start_time;
1302 A0 = sq_norm_W - 2*dot_prod_WoldW + sq_norm_oldW;
1303 B0 = dot_prod_WoldW - sq_norm_oldW;
1305 for(i=0; i < nData; i++)
1307 y2 = (uint32_t)data_y[i]-1;
1309 for(y=0; y < nY; y++)
1311 A[LIBOCAS_INDEX(y,i,nY)] = C*(output[LIBOCAS_INDEX(y,i,nY)] - old_output[LIBOCAS_INDEX(y,i,nY)]
1312 + old_output[LIBOCAS_INDEX(y2,i,nY)] - output[LIBOCAS_INDEX(y2,i,nY)]);
1313 B[LIBOCAS_INDEX(y,i,nY)] = C*(old_output[LIBOCAS_INDEX(y,i,nY)] - old_output[LIBOCAS_INDEX(y2,i,nY)]
1324 for(i=0; i < nData; i++)
1326 findactive(theta,sortedA,&nSortedA,&A[i*nY],&B[i*nY],nY,sort);
1329 while( idx < nSortedA-1 && theta[idx] < 0 )
1332 grad_sum += sortedA[idx];
1334 for(j=idx; j < nSortedA-1; j++)
1336 Theta[cnt1] = theta[j];
1340 for(j=idx+1; j < nSortedA; j++)
1342 Add[cnt2] = -sortedA[j-1]+sortedA[j];
1348 sort(Theta,Add,cnt1);
1349 ocas.sort_time +=
get_time() - start_time;
1361 for(i=0; i < cnt1; i++)
1365 grad = x*A0 + grad_sum;
1370 min_x = (grad*old_x - old_grad*x)/(grad - old_grad);
1376 grad_sum = grad_sum + Add[i];
1378 grad = x*A0 + grad_sum;
1398 sq_norm_W = update_W( t1, user_data );
1403 for(i=0; i < nData; i++)
1405 y2 = (uint32_t)data_y[i]-1;
1407 for(xi=-LIBOCAS_PLUS_INF, y=0; y < nY; y++)
1409 tmp = old_output[LIBOCAS_INDEX(y,i,nY)]*(1-t2) + t2*output[LIBOCAS_INDEX(y,i,nY)];
1410 if(y2 != y && xi < tmp)
1417 tmp = old_output[LIBOCAS_INDEX(y2,i,nY)]*(1-t2) + t2*output[LIBOCAS_INDEX(y2,i,nY)];
1418 xi = LIBOCAS_MAX(0,xi+1-tmp);
1431 for(i=0; i < nData; i++)
1433 y2 = (uint32_t)data_y[i]-1;
1435 for(tmp=-LIBOCAS_PLUS_INF, y=0; y < nY; y++)
1437 output[LIBOCAS_INDEX(y,i,nY)] = old_output[LIBOCAS_INDEX(y,i,nY)]*(1-t1) + t1*output[LIBOCAS_INDEX(y,i,nY)];
1439 if(y2 != y && tmp < output[LIBOCAS_INDEX(y,i,nY)])
1442 tmp = output[LIBOCAS_INDEX(y,i,nY)];
1446 R += LIBOCAS_MAX(0,1+tmp - output[LIBOCAS_INDEX(y2,i,nY)]);
1447 if( tmp >= output[LIBOCAS_INDEX(y2,i,nY)])
1451 ocas.Q_P = 0.5*sq_norm_W + C*R;
1455 ocas.ocas_time =
get_time() - ocas_start_time;
1459 ocas.print_time +=
get_time() - start_time;
1466 if( ocas.Q_P - ocas.Q_D <= TolRel*LIBOCAS_ABS(ocas.Q_P)) ocas.exitflag = 1;
1467 if( ocas.Q_P - ocas.Q_D <= TolAbs) ocas.exitflag = 2;
1468 if( ocas.Q_P <= QPBound) ocas.exitflag = 3;
1469 if( ocas.ocas_time >= MaxTime) ocas.exitflag = 4;
1470 if(ocas.nCutPlanes >= BufSize) ocas.exitflag = -1;
1478 LIBOCAS_FREE(alpha);
1479 LIBOCAS_FREE(new_cut);
1481 LIBOCAS_FREE(diag_H);
1482 LIBOCAS_FREE(output);
1483 LIBOCAS_FREE(old_output);
1486 LIBOCAS_FREE(theta);
1487 LIBOCAS_FREE(Theta);
1488 LIBOCAS_FREE(sortedA);
1491 ocas.ocas_time =
get_time() - ocas_start_time;