25 #define HISTORY_BUF 1000000
27 #define MINUS_INF INT_MIN
28 #define PLUS_INF INT_MAX
30 #define INDEX(ROW,COL,DIM) ((COL*DIM)+ROW)
109 float64_t min_beta1, max_beta1, min_beta2, max_beta2, beta;
122 int64_t History_size;
130 if( Ha1 == NULL )
SG_ERROR(
"Not enough memory.\n")
132 if( Ha2 == NULL )
SG_ERROR(
"Not enough memory.\n")
135 History = SG_MALLOC(
float64_t, History_size*2);
136 if( History == NULL )
SG_ERROR(
"Not enough memory.\n")
139 v1 = -1; v2 = -1; i = 0;
140 while( (v1 == -1 || v2 == -1) && i < dim ) {
141 if( v1 == -1 && vector_y[i] == 1 ) { v1 = i; }
142 if( v2 == -1 && vector_y[i] == 2 ) { v2 = i; }
156 for( i = 0; i < dim; i++ )
162 beta = Ha1[i] + Ha2[i] + vector_c[i];
164 if( vector_y[i] == 1 && min_beta1 > beta ) {
169 if( vector_y[i] == 2 && min_beta2 > beta ) {
178 UB = 0.5*(aHa11 + 2*aHa12 + aHa22) + ac1 + ac2;
179 LB = min_beta1 + min_beta2 - 0.5*(aHa11 + 2*aHa12 + aHa22);
181 delta1 = Ha1[v1] + Ha2[v1] + vector_c[v1] - min_beta1;
182 delta2 = Ha1[v2] + Ha2[v2] + vector_c[v2] - min_beta2;
185 History[
INDEX(0,0,2)] = LB;
186 History[
INDEX(1,0,2)] = UB;
189 SG_PRINT(
"Init: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
190 UB, LB, UB-LB,(UB-LB)/UB);
194 if( UB-LB <= tolabs ) exitflag = 1;
195 else if(UB-LB <=
CMath::abs(UB)*tolrel ) exitflag = 2;
196 else if(LB > th) exitflag = 3;
203 while( exitflag == -1 )
207 if( delta1 > delta2 )
216 lambda = delta1/(alpha[v1]*(Huu - 2*Huv + Hvv ));
219 tmp = lambda*alpha[v1];
221 aHa11 = aHa11 + 2*tmp*(Ha1[u1]-Ha1[v1])+tmp*tmp*( Huu - 2*Huv + Hvv );
222 aHa12 = aHa12 + tmp*(Ha2[u1]-Ha2[v1]);
223 ac1 = ac1 + tmp*(vector_c[u1]-vector_c[v1]);
225 alpha[u1] = alpha[u1] + tmp;
226 alpha[v1] = alpha[v1] - tmp;
230 for( i = 0; i < dim; i ++ )
232 Ha1[i] = Ha1[i] + tmp*(col_u[i] - col_v[i]);
234 beta = Ha1[i] + Ha2[i] + vector_c[i];
235 if( vector_y[i] == 1 )
237 if( min_beta1 > beta ) { u1 = i; min_beta1 = beta; }
238 if( max_beta1 < beta && alpha[i] > 0 ) { v1 = i; max_beta1 = beta; }
242 if( min_beta2 > beta ) { u2 = i; min_beta2 = beta; }
243 if( max_beta2 < beta && alpha[i] > 0) { v2 = i; max_beta2 = beta; }
256 lambda = delta2/(alpha[v2]*( Huu - 2*Huv + Hvv ));
259 tmp = lambda*alpha[v2];
260 aHa22 = aHa22 + 2*tmp*( Ha2[u2]-Ha2[v2]) + tmp*tmp*( Huu - 2*Huv + Hvv);
261 aHa12 = aHa12 + tmp*(Ha1[u2]-Ha1[v2]);
262 ac2 = ac2 + tmp*( vector_c[u2]-vector_c[v2] );
264 alpha[u2] = alpha[u2] + tmp;
265 alpha[v2] = alpha[v2] - tmp;
269 for(i = 0; i < dim; i++ )
271 Ha2[i] = Ha2[i] + tmp*( col_u[i] - col_v[i] );
273 beta = Ha1[i] + Ha2[i] + vector_c[i];
275 if( vector_y[i] == 1 )
277 if( min_beta1 > beta ) { u1 = i; min_beta1 = beta; }
278 if( max_beta1 < beta && alpha[i] > 0 ) { v1 = i; max_beta1 = beta; }
282 if( min_beta2 > beta ) { u2 = i; min_beta2 = beta; }
283 if( max_beta2 < beta && alpha[i] > 0) { v2 = i; max_beta2 = beta; }
288 UB = 0.5*(aHa11 + 2*aHa12 + aHa22) + ac1 + ac2;
289 LB = min_beta1 + min_beta2 - 0.5*(aHa11 + 2*aHa12 + aHa22);
291 delta1 = Ha1[v1] + Ha2[v1] + vector_c[v1] - min_beta1;
292 delta2 = Ha1[v2] + Ha2[v2] + vector_c[v2] - min_beta2;
295 if( UB-LB <= tolabs ) exitflag = 1;
296 else if( UB-LB <=
CMath::abs(UB)*tolrel ) exitflag = 2;
297 else if(LB > th) exitflag = 3;
298 else if(t >= tmax) exitflag = 0;
300 if( verb && (t % verb) == 0) {
301 SG_PRINT(
"%d: UB=%f,LB=%f,UB-LB=%f,(UB-LB)/|UB|=%f\n",
302 t, UB, LB, UB-LB,(UB-LB)/UB);
306 if( t < History_size ) {
307 History[
INDEX(0,t,2)] = LB;
308 History[
INDEX(1,t,2)] = UB;
312 if( tmp_ptr == NULL )
SG_ERROR(
"Not enough memory.\n")
313 for( i = 0; i < History_size; i++ ) {
314 tmp_ptr[
INDEX(0,i,2)] = History[
INDEX(0,i,2)];
315 tmp_ptr[
INDEX(1,i,2)] = History[
INDEX(1,i,2)];
317 tmp_ptr[
INDEX(0,t,2)] = LB;
318 tmp_ptr[
INDEX(1,t,2)] = UB;
327 if(verb && (t % verb) ) {
328 SG_PRINT(
"Exit: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
329 UB, LB, UB-LB,(UB-LB)/UB);
336 (*ptr_aHa11) = aHa11;
337 (*ptr_aHa22) = aHa22;
338 (*ptr_History) = History;
374 float64_t min_beta1, max_beta1, min_beta2, max_beta2, beta;
388 int64_t History_size;
397 if( Ha1 == NULL )
SG_ERROR(
"Not enough memory.\n")
399 if( Ha2 == NULL )
SG_ERROR(
"Not enough memory.\n")
402 History = SG_MALLOC(
float64_t, History_size*2);
403 if( History == NULL )
SG_ERROR(
"Not enough memory.\n")
406 v1 = -1; v2 = -1; i = 0;
407 while( (v1 == -1 || v2 == -1) && i < dim ) {
408 if( v1 == -1 && vector_y[i] == 1 ) { v1 = i; }
409 if( v2 == -1 && vector_y[i] == 2 ) { v2 = i; }
423 for( i = 0; i < dim; i++ )
429 beta = Ha1[i] + Ha2[i] + vector_c[i];
431 if( vector_y[i] == 1 && min_beta1 > beta ) {
436 if( vector_y[i] == 2 && min_beta2 > beta ) {
445 UB = 0.5*(aHa11 + 2*aHa12 + aHa22) + ac1 + ac2;
446 LB = min_beta1 + min_beta2 - 0.5*(aHa11 + 2*aHa12 + aHa22);
448 delta1 = Ha1[v1] + Ha2[v1] + vector_c[v1] - min_beta1;
449 delta2 = Ha1[v2] + Ha2[v2] + vector_c[v2] - min_beta2;
452 History[
INDEX(0,0,2)] = LB;
453 History[
INDEX(1,0,2)] = UB;
456 SG_PRINT(
"Init: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
457 UB, LB, UB-LB,(UB-LB)/UB);
460 if( delta1 > delta2 )
474 if( UB-LB <= tolabs ) exitflag = 1;
475 else if(UB-LB <=
CMath::abs(UB)*tolrel ) exitflag = 2;
476 else if(LB > th) exitflag = 3;
483 while( exitflag == -1 )
487 if( which_case == 1 )
493 lambda = delta1/(alpha[v1]*(Huu - 2*Huv + Hvv ));
496 tmp = lambda*alpha[v1];
498 aHa11 = aHa11 + 2*tmp*(Ha1[u1]-Ha1[v1])+tmp*tmp*( Huu - 2*Huv + Hvv );
499 aHa12 = aHa12 + tmp*(Ha2[u1]-Ha2[v1]);
500 ac1 = ac1 + tmp*(vector_c[u1]-vector_c[v1]);
502 alpha[u1] = alpha[u1] + tmp;
503 alpha[v1] = alpha[v1] - tmp;
507 for( i = 0; i < dim; i ++ )
509 Ha1[i] = Ha1[i] + tmp*(col_u[i] - col_v[i]);
511 beta = Ha1[i] + Ha2[i] + vector_c[i];
512 if( vector_y[i] == 1 )
514 if( min_beta1 > beta ) { u1 = i; min_beta1 = beta; }
515 if( max_beta1 < beta && alpha[i] > 0 ) { v1 = i; max_beta1 = beta; }
519 if( min_beta2 > beta ) { u2 = i; min_beta2 = beta; }
520 if( max_beta2 < beta && alpha[i] > 0) { v2 = i; max_beta2 = beta; }
530 lambda = delta2/(alpha[v2]*( Huu - 2*Huv + Hvv ));
533 tmp = lambda*alpha[v2];
534 aHa22 = aHa22 + 2*tmp*( Ha2[u2]-Ha2[v2]) + tmp*tmp*( Huu - 2*Huv + Hvv);
535 aHa12 = aHa12 + tmp*(Ha1[u2]-Ha1[v2]);
536 ac2 = ac2 + tmp*( vector_c[u2]-vector_c[v2] );
538 alpha[u2] = alpha[u2] + tmp;
539 alpha[v2] = alpha[v2] - tmp;
543 for(i = 0; i < dim; i++ )
545 Ha2[i] = Ha2[i] + tmp*( col_u[i] - col_v[i] );
547 beta = Ha1[i] + Ha2[i] + vector_c[i];
549 if( vector_y[i] == 1 )
551 if( min_beta1 > beta ) { u1 = i; min_beta1 = beta; }
552 if( max_beta1 < beta && alpha[i] > 0 ) { v1 = i; max_beta1 = beta; }
556 if( min_beta2 > beta ) { u2 = i; min_beta2 = beta; }
557 if( max_beta2 < beta && alpha[i] > 0) { v2 = i; max_beta2 = beta; }
562 UB = 0.5*(aHa11 + 2*aHa12 + aHa22) + ac1 + ac2;
563 LB = min_beta1 + min_beta2 - 0.5*(aHa11 + 2*aHa12 + aHa22);
565 delta1 = Ha1[v1] + Ha2[v1] + vector_c[v1] - min_beta1;
566 delta2 = Ha1[v2] + Ha2[v2] + vector_c[v2] - min_beta2;
568 if( delta1 > delta2 )
573 for( max_improv =
MINUS_INF, i = 0; i < dim; i++ ) {
575 if( vector_y[i] == 1 && alpha[i] != 0 ) {
577 beta = Ha1[i] + Ha2[i] + vector_c[i];
579 if( beta >= min_beta1 ) {
581 tmp = diag_H[u1] - 2*col_u[i] + diag_H[i];
583 improv = (0.5*(beta-min_beta1)*(beta-min_beta1))/tmp;
585 if( improv > max_improv ) {
594 delta1 = Ha1[v1] + Ha2[v1] + vector_c[v1] - min_beta1;
603 for( max_improv =
MINUS_INF, i = 0; i < dim; i++ ) {
605 if( vector_y[i] == 2 && alpha[i] != 0 ) {
607 beta = Ha1[i] + Ha2[i] + vector_c[i];
609 if( beta >= min_beta2 ) {
611 tmp = diag_H[u2] - 2*col_u[i] + diag_H[i];
613 improv = (0.5*(beta-min_beta2)*(beta-min_beta2))/tmp;
615 if( improv > max_improv ) {
625 delta2 = Ha1[v2] + Ha2[v2] + vector_c[v2] - min_beta2;
631 if( UB-LB <= tolabs ) exitflag = 1;
632 else if( UB-LB <=
CMath::abs(UB)*tolrel ) exitflag = 2;
633 else if(LB > th) exitflag = 3;
634 else if(t >= tmax) exitflag = 0;
636 if( verb && (t % verb) == 0) {
637 SG_PRINT(
"%d: UB=%f,LB=%f,UB-LB=%f,(UB-LB)/|UB|=%f\n",
638 t, UB, LB, UB-LB,(UB-LB)/UB);
642 if( t < History_size ) {
643 History[
INDEX(0,t,2)] = LB;
644 History[
INDEX(1,t,2)] = UB;
648 if( tmp_ptr == NULL )
SG_ERROR(
"Not enough memory.\n")
649 for( i = 0; i < History_size; i++ ) {
650 tmp_ptr[
INDEX(0,i,2)] = History[
INDEX(0,i,2)];
651 tmp_ptr[
INDEX(1,i,2)] = History[
INDEX(1,i,2)];
653 tmp_ptr[
INDEX(0,t,2)] = LB;
654 tmp_ptr[
INDEX(1,t,2)] = UB;
663 if(verb && (t % verb) ) {
664 SG_PRINT(
"Exit: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
665 UB, LB, UB-LB,(UB-LB)/UB);
672 (*ptr_aHa11) = aHa11;
673 (*ptr_aHa22) = aHa22;
674 (*ptr_History) = History;
#define INDEX(ROW, COL, DIM)
float64_t ** kernel_columns
float64_t kernel(int32_t idx_a, int32_t idx_b)
int8_t gnpp_mdm(float64_t *diag_H, float64_t *vector_c, float64_t *vector_y, int32_t dim, int32_t tmax, float64_t tolabs, float64_t tolrel, float64_t th, float64_t *alpha, int32_t *ptr_t, float64_t *ptr_aHa11, float64_t *ptr_aHa22, float64_t **ptr_History, int32_t verb)
Class SGObject is the base class of all shogun objects.
float64_t * get_col(int64_t a, int64_t b)
all of classes and functions are contained in the shogun namespace
int8_t gnpp_imdm(float64_t *diag_H, float64_t *vector_c, float64_t *vector_y, int32_t dim, int32_t tmax, float64_t tolabs, float64_t tolrel, float64_t th, float64_t *alpha, int32_t *ptr_t, float64_t *ptr_aHa11, float64_t *ptr_aHa22, float64_t **ptr_History, int32_t verb)
#define SG_UNSTABLE(func,...)