00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071 #include <stdio.h>
00072 #include <stdlib.h>
00073 #include <string.h>
00074 #include <math.h>
00075 #include <shogun/lib/external/gpdt.h>
00076 #include <shogun/io/SGIO.h>
00077
00078 namespace shogun
00079 {
00080 #define maxvpm 30000
00081 #define maxprojections 200
00082 #define in 8000
00083 #define alpha_max 1e10
00084 #define alpha_min 1e-10
00085
00086 extern uint32_t Randnext;
00087 #define ThRand (Randnext = Randnext * 1103515245L + 12345L)
00088 #define ThRandPos ((Randnext = Randnext * 1103515245L + 12345L) & 0x7fffffff)
00089
00090 int32_t InnerProjector(
00091 int32_t method, int32_t n, int32_t *iy, float64_t e, float64_t *qk, float64_t l,
00092 float64_t u, float64_t *x, float64_t &lambda);
00093
00094
00095
00096
00097
00098
00099
00100
00101 #define VPM_ADA
00102
00103
00104
00105
00106
00107
00108
00109
00110 int32_t gvpm(
00111 int32_t Projector, int32_t n, float32_t *vecA, float64_t *b, float64_t c,
00112 float64_t e, int32_t *iy, float64_t *x, float64_t tol, int32_t *ls,
00113 int32_t *proj)
00114 {
00115 int32_t i, j, iter, it, it2, luv, info;
00116 float64_t gd, max, normd, dAd, lam, lamnew, alpha, kktlam, ak, bk;
00117
00118 int32_t lscount = 0, projcount = 0;
00119 float64_t eps = 1.0e-16;
00120 float64_t DELTAsv, ProdDELTAsv;
00121 float64_t lam_ext;
00122
00123
00124 #ifdef VPM_ADA
00125 int32_t nc = 1, ia1 = -1;
00126 float64_t alpha1, alpha2;
00127 #endif
00128
00129
00130 #ifdef VARIABLES_ON_STACK
00131
00132 int32_t ipt[in], ipt2[in], uv[in];
00133 float64_t g[in], y[in], tempv[in], d[in], Ad[in], t[in];
00134
00135 #else
00136
00137 int32_t *ipt, *ipt2, *uv;
00138 float64_t *g, *y, *tempv, *d, *Ad, *t;
00139
00140
00141 ipt = SG_MALLOC(int32_t, n);
00142 ipt2 = SG_MALLOC(int32_t, n);
00143 uv = SG_MALLOC(int32_t, n);
00144 g = SG_MALLOC(float64_t, n);
00145 y = SG_MALLOC(float64_t, n);
00146 d = SG_MALLOC(float64_t, n);
00147 Ad = SG_MALLOC(float64_t, n);
00148 t = SG_MALLOC(float64_t, n);
00149 tempv = SG_MALLOC(float64_t, n);
00150
00151 #endif
00152
00153 DELTAsv = EPS_SV * c;
00154 if (tol <= 1.0e-5 || n <= 20)
00155 ProdDELTAsv = 0.0F;
00156 else
00157 ProdDELTAsv = EPS_SV * c;
00158
00159 for (i = 0; i < n; i++)
00160 tempv[i] = -x[i];
00161 lam_ext = 0.0;
00162 projcount += InnerProjector(Projector, n, iy, e, tempv, 0, c, x, lam_ext);
00163
00164
00165
00166 {
00167 float32_t *tempA;
00168
00169 it = 0;
00170 for (i = 0; i < n; i++)
00171 if (fabs(x[i]) > ProdDELTAsv*1e-2)
00172 ipt[it++] = i;
00173
00174 memset(t, 0, n*sizeof(float64_t));
00175 for (i = 0; i < it; i++)
00176 {
00177 tempA = vecA + ipt[i]*n;
00178 for (j = 0; j < n; j++)
00179 t[j] += (tempA[j] * x[ipt[i]]);
00180 }
00181 }
00182
00183 for (i = 0; i < n; i++)
00184 {
00185 g[i] = t[i] + b[i],
00186 y[i] = g[i] - x[i];
00187 }
00188
00189 projcount += InnerProjector(Projector, n, iy, e, y, 0, c, tempv, lam_ext);
00190
00191 max = alpha_min;
00192 for (i = 0; i < n; i++)
00193 {
00194 y[i] = tempv[i] - x[i];
00195 if (fabs(y[i]) > max)
00196 max = fabs(y[i]);
00197 }
00198
00199 if (max < c*tol*1e-3)
00200 {
00201 lscount = 0;
00202 iter = 0;
00203 goto Clean;
00204 }
00205
00206 alpha = 1.0 / max;
00207
00208 for (iter = 1; iter <= maxvpm; iter++)
00209 {
00210 for (i = 0; i < n; i++)
00211 tempv[i] = alpha*g[i] - x[i];
00212
00213 projcount += InnerProjector(Projector, n, iy, e, tempv, 0, c, y, lam_ext);
00214
00215 gd = 0.0;
00216 for (i = 0; i < n; i++)
00217 {
00218 d[i] = y[i] - x[i];
00219 gd += d[i] * g[i];
00220 }
00221
00222
00223
00224 {
00225 float32_t *tempA;
00226
00227 it = it2 = 0;
00228 for (i = 0; i < n; i++)
00229 if (fabs(d[i]) > (ProdDELTAsv*1.0e-2))
00230 ipt[it++] = i;
00231 for (i = 0; i < n; i++)
00232 if (fabs(y[i]) > ProdDELTAsv)
00233 ipt2[it2++] = i;
00234
00235 memset(Ad, 0, n*sizeof(float64_t));
00236 if (it < it2)
00237 {
00238 for (i = 0; i < it; i++)
00239 {
00240 tempA = vecA + ipt[i]*n;
00241 for (j = 0; j < n; j++)
00242 Ad[j] += (tempA[j] * d[ipt[i]]);
00243 }
00244 }
00245 else
00246 {
00247 for (i = 0; i < it2; i++)
00248 {
00249 tempA = vecA + ipt2[i]*n;
00250 for (j = 0; j < n; j++)
00251 Ad[j] += (tempA[j] * y[ipt2[i]]);
00252 }
00253 for (j = 0; j < n; j++)
00254 Ad[j] -= t[j];
00255 }
00256 }
00257
00258 normd = 0.0;
00259 for (i = 0; i < n; i++)
00260 normd += d[i] * d[i];
00261
00262 dAd = 0.0;
00263 for (i = 0; i < n; i++)
00264 dAd += d[i]*Ad[i];
00265
00266 if (dAd > eps*normd && gd < 0.0)
00267 {
00268 lam = lamnew = -gd/dAd;
00269 if (lam > 1.0 || lam < 0.0)
00270 lam = 1.0;
00271 else
00272 lscount++;
00273
00274 #ifdef VPM_ADA
00275
00276
00277
00278
00279 alpha1 = normd / dAd;
00280
00281
00282 alpha2 = 0.0;
00283 for (i = 0; i < n; i++)
00284 alpha2 += Ad[i] * Ad[i];
00285 alpha2 = dAd / alpha2;
00286
00287 if ( (nc > 2
00288 && (
00289 (ia1 == 1
00290 && (
00291 lamnew < 0.1 || (alpha1 > alpha && alpha2 < alpha)
00292 )
00293 )
00294 ||
00295 (ia1 == -1
00296 && (
00297 lamnew > 5.0 || (alpha1 > alpha && alpha2 < alpha)
00298 )
00299 )
00300 )
00301 )
00302 || nc > 9 )
00303 {
00304 ia1 = -ia1;
00305 nc = 0;
00306 }
00307
00308 if (ia1 == 1)
00309 alpha = alpha1;
00310 else
00311 alpha = alpha2;
00312
00313 if (alpha < alpha_min)
00314 alpha = alpha_min;
00315 else if (alpha > alpha_max)
00316 alpha = alpha_max;
00317
00318 nc++;
00319
00320 #else
00321
00322
00323
00324 if ((iter % 6) < 3)
00325 {
00326 alpha = 0.0;
00327 for (i = 0; i < n; i++)
00328 alpha += Ad[i] * Ad[i];
00329 alpha = dAd / alpha;
00330 }
00331 else
00332 {
00333 alpha = 0.0;
00334 for (i = 0; i < n; i++)
00335 alpha += d[i] * d[i];
00336 alpha = alpha / dAd;
00337 }
00338
00339 #endif
00340
00341 }
00342 else
00343 {
00344 lam = 1.0;
00345 alpha = alpha_max;
00346 }
00347
00348 for (i = 0; i < n; i++)
00349 {
00350 x[i] = x[i] + lam*d[i];
00351 t[i] = t[i] + lam*Ad[i];
00352 g[i] = b[i] + t[i];
00353 }
00354
00355
00356 bk = 0.0;
00357 ak = 0.0;
00358 for (i = 0; i < n; i++)
00359 {
00360 bk += x[i] * x[i];
00361 ak += d[i] * d[i];
00362 }
00363
00364 if (lam*sqrt(ak) < tol*10 * sqrt(bk))
00365 {
00366 it = 0;
00367 luv = 0;
00368 kktlam = 0.0;
00369 for (i = 0; i < n; i++)
00370 {
00371 if (x[i] > DELTAsv && x[i] < c-DELTAsv)
00372 {
00373 ipt[it++] = i;
00374 kktlam = kktlam - iy[i]*g[i];
00375 }
00376 else
00377 uv[luv++] = i;
00378 }
00379
00380 if (it == 0)
00381 {
00382 if (lam*sqrt(ak) < tol*0.5 * sqrt(bk))
00383 goto Clean;
00384 }
00385 else
00386 {
00387 kktlam = kktlam/it;
00388 info = 1;
00389 for (i = 0; i < it; i++)
00390 if (fabs(iy[ipt[i]]*g[ipt[i]]+kktlam) > tol)
00391 {
00392 info = 0;
00393 break;
00394 }
00395
00396 if (info == 1)
00397 for (i = 0; i < luv; i++)
00398 if (x[uv[i]] <= DELTAsv)
00399 {
00400 if (g[uv[i]] + kktlam*iy[uv[i]] < -tol)
00401 {
00402 info = 0;
00403 break;
00404 }
00405 }
00406 else
00407 {
00408 if (g[uv[i]] + kktlam*iy[uv[i]] > tol)
00409 {
00410 info = 0;
00411 break;
00412 }
00413 }
00414
00415 if (info == 1)
00416 goto Clean;
00417 }
00418 }
00419 }
00420
00421 SG_SWARNING("GVPM exits after maxvpm = %d iterations.\n", maxvpm);
00422
00423 Clean:
00424
00425
00426 #ifndef VARIABLES_ON_STACK
00427 SG_FREE(t);
00428 SG_FREE(uv);
00429 SG_FREE(ipt2);
00430 SG_FREE(ipt);
00431 SG_FREE(g);
00432 SG_FREE(y);
00433 SG_FREE(tempv);
00434 SG_FREE(d);
00435 SG_FREE(Ad);
00436 #endif
00437
00438 if (ls != NULL) *ls = lscount;
00439 if (proj != NULL) *proj = projcount;
00440 return(iter);
00441 }
00442
00443
00444
00445
00446
00447
00448 int32_t FletcherAlg2A(
00449 int32_t Projector, int32_t n, float32_t *vecA, float64_t *b, float64_t c,
00450 float64_t e, int32_t *iy, float64_t *x, float64_t tol, int32_t *ls,
00451 int32_t *proj)
00452 {
00453 int32_t i, j, iter, it, it2, luv, info, lscount = 0, projcount = 0;
00454 float64_t gd, max, ak, bk, akold, bkold, lamnew, alpha, kktlam, lam_ext;
00455 float64_t eps = 1.0e-16;
00456 float64_t DELTAsv, ProdDELTAsv;
00457
00458
00459 int32_t L, llast;
00460 float64_t fr, fbest, fv, fc, fv0;
00461
00462
00463 #ifdef VARIABLES_ON_STACK
00464
00465 int32_t ipt[in], ipt2[in], uv[in];
00466 float64_t g[in], y[in], tempv[in], d[in], Ad[in], t[in],
00467 xplus[in], tplus[in], sk[in], yk[in];
00468 #else
00469
00470 int32_t *ipt, *ipt2, *uv;
00471 float64_t *g, *y, *tempv, *d, *Ad, *t, *xplus, *tplus, *sk, *yk;
00472
00473
00474 ipt = SG_MALLOC(int32_t, n);
00475 ipt2 = SG_MALLOC(int32_t, n);
00476 uv = SG_MALLOC(int32_t, n);
00477 g = SG_MALLOC(float64_t, n);
00478 y = SG_MALLOC(float64_t, n);
00479 tempv = SG_MALLOC(float64_t, n);
00480 d = SG_MALLOC(float64_t, n);
00481 Ad = SG_MALLOC(float64_t, n);
00482 t = SG_MALLOC(float64_t, n);
00483 xplus = SG_MALLOC(float64_t, n);
00484 tplus = SG_MALLOC(float64_t, n);
00485 sk = SG_MALLOC(float64_t, n);
00486 yk = SG_MALLOC(float64_t, n);
00487
00488 #endif
00489
00490 DELTAsv = EPS_SV * c;
00491 if (tol <= 1.0e-5 || n <= 20)
00492 ProdDELTAsv = 0.0F;
00493 else
00494 ProdDELTAsv = EPS_SV * c;
00495
00496 for (i = 0; i < n; i++)
00497 tempv[i] = -x[i];
00498
00499 lam_ext = 0.0;
00500
00501 projcount += InnerProjector(Projector, n, iy, e, tempv, 0, c, x, lam_ext);
00502
00503
00504
00505 {
00506 float32_t *tempA;
00507
00508 it = 0;
00509 for (i = 0; i < n; i++)
00510 if (fabs(x[i]) > ProdDELTAsv)
00511 ipt[it++] = i;
00512
00513 memset(t, 0, n*sizeof(float64_t));
00514 for (i = 0; i < it; i++)
00515 {
00516 tempA = vecA + ipt[i] * n;
00517 for (j = 0; j < n; j++)
00518 t[j] += (tempA[j]*x[ipt[i]]);
00519 }
00520 }
00521
00522 for (i = 0; i < n; i++)
00523 {
00524 g[i] = t[i] + b[i],
00525 y[i] = g[i] - x[i];
00526 }
00527
00528 projcount += InnerProjector(Projector, n, iy, e, y, 0, c, tempv, lam_ext);
00529
00530 max = alpha_min;
00531 for (i = 0; i < n; i++)
00532 {
00533 y[i] = tempv[i] - x[i];
00534 if (fabs(y[i]) > max)
00535 max = fabs(y[i]);
00536 }
00537
00538 if (max < c*tol*1e-3)
00539 {
00540 lscount = 0;
00541 iter = 0;
00542 goto Clean;
00543 }
00544
00545 alpha = 1.0 / max;
00546
00547 fv0 = 0.0;
00548 for (i = 0; i < n; i++)
00549 fv0 += x[i] * (0.5*t[i] + b[i]);
00550
00551
00552 L = 2;
00553 fr = alpha_max;
00554 fbest = fv0;
00555 fc = fv0;
00556 llast = 0;
00557 akold = bkold = 0.0;
00558
00559 for (iter = 1; iter <= maxvpm; iter++)
00560 {
00561 for (i = 0; i < n; i++)
00562 tempv[i] = alpha*g[i] - x[i];
00563
00564 projcount += InnerProjector(Projector, n, iy, e, tempv, 0, c, y, lam_ext);
00565
00566 gd = 0.0;
00567 for (i = 0; i < n; i++)
00568 {
00569 d[i] = y[i] - x[i];
00570 gd += d[i] * g[i];
00571 }
00572
00573
00574 {
00575 float32_t *tempA;
00576
00577 it = it2 = 0;
00578 for (i = 0; i < n; i++)
00579 if (fabs(d[i]) > (ProdDELTAsv*1.0e-2))
00580 ipt[it++] = i;
00581 for (i = 0; i < n; i++)
00582 if (fabs(y[i]) > ProdDELTAsv)
00583 ipt2[it2++] = i;
00584
00585 memset(Ad, 0, n*sizeof(float64_t));
00586 if (it < it2)
00587 {
00588 for (i = 0; i < it; i++)
00589 {
00590 tempA = vecA + ipt[i]*n;
00591 for (j = 0; j < n; j++)
00592 Ad[j] += (tempA[j] * d[ipt[i]]);
00593 }
00594 }
00595 else
00596 {
00597 for (i = 0; i < it2; i++)
00598 {
00599 tempA = vecA + ipt2[i]*n;
00600 for (j = 0; j < n; j++)
00601 Ad[j] += (tempA[j] * y[ipt2[i]]);
00602 }
00603 for (j = 0; j < n; j++)
00604 Ad[j] -= t[j];
00605 }
00606 }
00607
00608 ak = 0.0;
00609 for (i = 0; i < n; i++)
00610 ak += d[i] * d[i];
00611
00612 bk = 0.0;
00613 for (i = 0; i < n; i++)
00614 bk += d[i]*Ad[i];
00615
00616 if (bk > eps*ak && gd < 0.0)
00617 lamnew = -gd/bk;
00618 else
00619 lamnew = 1.0;
00620
00621 fv = 0.0;
00622 for (i = 0; i < n; i++)
00623 {
00624 xplus[i] = x[i] + d[i];
00625 tplus[i] = t[i] + Ad[i];
00626 fv += xplus[i] * (0.5*tplus[i] + b[i]);
00627 }
00628
00629 if ((iter == 1 && fv >= fv0) || (iter > 1 && fv >= fr))
00630 {
00631 lscount++;
00632 fv = 0.0;
00633 for (i = 0; i < n; i++)
00634 {
00635 xplus[i] = x[i] + lamnew*d[i];
00636 tplus[i] = t[i] + lamnew*Ad[i];
00637 fv += xplus[i] * (0.5*tplus[i] + b[i]);
00638 }
00639 }
00640
00641 for (i = 0; i < n; i++)
00642 {
00643 sk[i] = xplus[i] - x[i];
00644 yk[i] = tplus[i] - t[i];
00645 x[i] = xplus[i];
00646 t[i] = tplus[i];
00647 g[i] = t[i] + b[i];
00648 }
00649
00650
00651
00652 if (fv < fbest)
00653 {
00654 fbest = fv;
00655 fc = fv;
00656 llast = 0;
00657 }
00658 else
00659 {
00660 fc = (fc > fv ? fc : fv);
00661 llast++;
00662 if (llast == L)
00663 {
00664 fr = fc;
00665 fc = fv;
00666 llast = 0;
00667 }
00668 }
00669
00670 ak = bk = 0.0;
00671 for (i = 0; i < n; i++)
00672 {
00673 ak += sk[i] * sk[i];
00674 bk += sk[i] * yk[i];
00675 }
00676
00677 if (bk < eps*ak)
00678 alpha = alpha_max;
00679 else
00680 {
00681 if (bkold < eps*akold)
00682 alpha = ak/bk;
00683 else
00684 alpha = (akold+ak)/(bkold+bk);
00685
00686 if (alpha > alpha_max)
00687 alpha = alpha_max;
00688 else if (alpha < alpha_min)
00689 alpha = alpha_min;
00690 }
00691
00692 akold = ak;
00693 bkold = bk;
00694
00695
00696
00697 bk = 0.0;
00698 for (i = 0; i < n; i++)
00699 bk += x[i] * x[i];
00700
00701 if (sqrt(ak) < tol*10 * sqrt(bk))
00702 {
00703 it = 0;
00704 luv = 0;
00705 kktlam = 0.0;
00706 for (i = 0; i < n; i++)
00707 {
00708 if ((x[i] > DELTAsv) && (x[i] < c-DELTAsv))
00709 {
00710 ipt[it++] = i;
00711 kktlam = kktlam - iy[i]*g[i];
00712 }
00713 else
00714 uv[luv++] = i;
00715 }
00716
00717 if (it == 0)
00718 {
00719 if (sqrt(ak) < tol*0.5 * sqrt(bk))
00720 goto Clean;
00721 }
00722 else
00723 {
00724
00725 kktlam = kktlam/it;
00726 info = 1;
00727 for (i = 0; i < it; i++)
00728 if ( fabs(iy[ipt[i]] * g[ipt[i]] + kktlam) > tol )
00729 {
00730 info = 0;
00731 break;
00732 }
00733
00734 if (info == 1)
00735 for (i = 0; i < luv; i++)
00736 if (x[uv[i]] <= DELTAsv)
00737 {
00738 if (g[uv[i]] + kktlam*iy[uv[i]] < -tol)
00739 {
00740 info = 0;
00741 break;
00742 }
00743 }
00744 else
00745 {
00746 if (g[uv[i]] + kktlam*iy[uv[i]] > tol)
00747 {
00748 info = 0;
00749 break;
00750 }
00751 }
00752
00753 if (info == 1)
00754 goto Clean;
00755 }
00756 }
00757 }
00758
00759 SG_SWARNING("Dai-Fletcher method exits after maxvpm = %d iterations.\n", maxvpm);
00760
00761 Clean:
00762
00763 #ifndef VARIABLES_ON_STACK
00764 SG_FREE(sk);
00765 SG_FREE(yk);
00766 SG_FREE(tplus);
00767 SG_FREE(xplus);
00768 SG_FREE(t);
00769 SG_FREE(uv);
00770 SG_FREE(ipt2);
00771 SG_FREE(ipt);
00772 SG_FREE(g);
00773 SG_FREE(y);
00774 SG_FREE(tempv);
00775 SG_FREE(d);
00776 SG_FREE(Ad);
00777 #endif
00778
00779 if (ls != NULL) *ls = lscount;
00780 if (proj != NULL) *proj = projcount;
00781 return(iter);
00782
00783 }
00784
00785
00786
00787
00788 int32_t gpm_solver(
00789 int32_t Solver, int32_t Projector, int32_t n, float32_t *A, float64_t *b,
00790 float64_t c, float64_t e, int32_t *iy, float64_t *x, float64_t tol,
00791 int32_t *ls, int32_t *proj)
00792 {
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820 if (Solver == SOLVER_FLETCHER)
00821 return FletcherAlg2A(Projector, n, A, b, c, e, iy, x, tol, ls, proj);
00822 else
00823 return gvpm(Projector, n, A, b, c, e, iy, x, tol, ls, proj);
00824 }
00825
00826
00827
00828
00829
00830
00831
00832 float64_t ProjectR(
00833 float64_t *x, int32_t n, float64_t lambda, int32_t *a, float64_t b,
00834 float64_t *c, float64_t l, float64_t u)
00835 {
00836 int32_t i;
00837 float64_t r = 0.0;
00838
00839 for (i = 0; i < n; i++)
00840 {
00841 x[i] = -c[i] + lambda*(float64_t)a[i];
00842 if (x[i] >= u) x[i] = u;
00843 else if (x[i] < l) x[i] = l;
00844 r += (float64_t)a[i]*x[i];
00845 }
00846
00847 return (r - b);
00848 }
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860 int32_t ProjectDai(
00861 int32_t n, int32_t *a, float64_t b, float64_t *c, float64_t l, float64_t u,
00862 float64_t *x, float64_t &lam_ext)
00863 {
00864 float64_t lambda, lambdal, lambdau, dlambda, lambda_new, tol_lam;
00865 float64_t r, rl, ru, s, tol_r;
00866 int32_t iter;
00867
00868 tol_lam = 1.0e-11;
00869 tol_r = 1.0e-10 * sqrt((u-l)*(float64_t)n);
00870 lambda = lam_ext;
00871 dlambda = 0.5;
00872 iter = 1;
00873 b = -b;
00874
00875
00876 r = ProjectR(x, n, lambda, a, b, c, l, u);
00877 if (fabs(r) < tol_r)
00878 return 0;
00879
00880 if (r < 0.0)
00881 {
00882 lambdal = lambda;
00883 rl = r;
00884 lambda = lambda + dlambda;
00885 r = ProjectR(x, n, lambda, a, b, c, l, u);
00886 while (r < 0.0)
00887 {
00888 lambdal = lambda;
00889 s = rl/r - 1.0;
00890 if (s < 0.1) s = 0.1;
00891 dlambda = dlambda + dlambda/s;
00892 lambda = lambda + dlambda;
00893 rl = r;
00894 r = ProjectR(x, n, lambda, a, b, c, l, u);
00895 }
00896 lambdau = lambda;
00897 ru = r;
00898 }
00899 else
00900 {
00901 lambdau = lambda;
00902 ru = r;
00903 lambda = lambda - dlambda;
00904 r = ProjectR(x, n, lambda, a, b, c, l, u);
00905 while (r > 0.0)
00906 {
00907 lambdau = lambda;
00908 s = ru/r - 1.0;
00909 if (s < 0.1) s = 0.1;
00910 dlambda = dlambda + dlambda/s;
00911 lambda = lambda - dlambda;
00912 ru = r;
00913 r = ProjectR(x, n, lambda, a, b, c, l, u);
00914 }
00915 lambdal = lambda;
00916 rl = r;
00917 }
00918
00919
00920
00921 s = 1.0 - rl/ru;
00922 dlambda = dlambda/s;
00923 lambda = lambdau - dlambda;
00924 r = ProjectR(x, n, lambda, a, b, c, l, u);
00925
00926 while ( fabs(r) > tol_r
00927 && dlambda > tol_lam * (1.0 + fabs(lambda))
00928 && iter < maxprojections )
00929 {
00930 iter++;
00931 if (r > 0.0)
00932 {
00933 if (s <= 2.0)
00934 {
00935 lambdau = lambda;
00936 ru = r;
00937 s = 1.0 - rl/ru;
00938 dlambda = (lambdau - lambdal) / s;
00939 lambda = lambdau - dlambda;
00940 }
00941 else
00942 {
00943 s = ru/r-1.0;
00944 if (s < 0.1) s = 0.1;
00945 dlambda = (lambdau - lambda) / s;
00946 lambda_new = 0.75*lambdal + 0.25*lambda;
00947 if (lambda_new < (lambda - dlambda))
00948 lambda_new = lambda - dlambda;
00949 lambdau = lambda;
00950 ru = r;
00951 lambda = lambda_new;
00952 s = (lambdau - lambdal) / (lambdau - lambda);
00953 }
00954 }
00955 else
00956 {
00957 if (s >= 2.0)
00958 {
00959 lambdal = lambda;
00960 rl = r;
00961 s = 1.0 - rl/ru;
00962 dlambda = (lambdau - lambdal) / s;
00963 lambda = lambdau - dlambda;
00964 }
00965 else
00966 {
00967 s = rl/r - 1.0;
00968 if (s < 0.1) s = 0.1;
00969 dlambda = (lambda-lambdal) / s;
00970 lambda_new = 0.75*lambdau + 0.25*lambda;
00971 if (lambda_new > (lambda + dlambda))
00972 lambda_new = lambda + dlambda;
00973 lambdal = lambda;
00974 rl = r;
00975 lambda = lambda_new;
00976 s = (lambdau - lambdal) / (lambdau-lambda);
00977 }
00978 }
00979 r = ProjectR(x, n, lambda, a, b, c, l, u);
00980 }
00981
00982 lam_ext = lambda;
00983 if (iter >= maxprojections)
00984 SG_SERROR("Projector exits after max iterations: %d\n", iter);
00985
00986 return (iter);
00987 }
00988
00989 #define SWAP(a,b) { register float64_t t=(a);(a)=(b);(b)=t; }
00990
00991
00992 float64_t quick_select(float64_t *arr, int32_t n)
00993 {
00994 int32_t low, high;
00995 int32_t median;
00996 int32_t middle, l, h;
00997
00998 low = 0;
00999 high = n-1;
01000 median = (low + high) / 2;
01001
01002 for (;;)
01003 {
01004 if (high <= low)
01005 return arr[median];
01006
01007 if (high == low + 1)
01008 {
01009 if (arr[low] > arr[high])
01010 SWAP(arr[low], arr[high]);
01011 return arr[median];
01012 }
01013
01014 middle = (low + high) / 2;
01015 if (arr[middle] > arr[high]) SWAP(arr[middle], arr[high]);
01016 if (arr[low] > arr[high]) SWAP(arr[low], arr[high]);
01017 if (arr[middle] > arr[low]) SWAP(arr[middle], arr[low]);
01018
01019 SWAP(arr[middle], arr[low+1]);
01020
01021 l = low + 1;
01022 h = high;
01023 for (;;)
01024 {
01025 do l++; while (arr[low] > arr[l]);
01026 do h--; while (arr[h] > arr[low]);
01027 if (h < l)
01028 break;
01029 SWAP(arr[l], arr[h]);
01030 }
01031
01032 SWAP(arr[low], arr[h]);
01033 if (h <= median)
01034 low = l;
01035 if (h >= median)
01036 high = h - 1;
01037 }
01038 }
01039
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049
01050
01051
01052 int32_t Pardalos(
01053 int32_t n, int32_t *iy, float64_t e, float64_t *qk, float64_t low,
01054 float64_t up, float64_t *x)
01055 {
01056 int32_t i, l, iter;
01057 int32_t luv, lxint;
01058 float64_t d, xmin, xmax, xmold, xmid, xx, ts, sw, s, s1, testsum;
01059
01060
01061 #ifdef VARIABLES_ON_STACK
01062 int32_t uv[in], uvt[in];
01063 float64_t xint[2*in+2], xint2[2*in+2], a[in], b[in], at[in], bt[in];
01064 float64_t newdia[in], newdt[in];
01065 #else
01066
01067 int32_t *uv, *uvt;
01068 float64_t *xint, *xint2, *a, *b, *at, *bt, *newdia, *newdt;
01069
01070
01071 uv = SG_MALLOC(int32_t, n);
01072 uvt = SG_MALLOC(int32_t, n);
01073 a = SG_MALLOC(float64_t, n);
01074 b = SG_MALLOC(float64_t, n);
01075 at = SG_MALLOC(float64_t, n);
01076 bt = SG_MALLOC(float64_t, n);
01077 newdia = SG_MALLOC(float64_t, n);
01078 newdt = SG_MALLOC(float64_t, n);
01079 xint = SG_MALLOC(float64_t, (2*n + 2));
01080 xint2 = SG_MALLOC(float64_t, (2*n + 2));
01081
01082 #endif
01083
01084 d = 0.0;
01085 for (i = 0; i < n; i++)
01086 d += iy[i] * qk[i];
01087 d = 0.5 * (d-e);
01088
01089 for (i = 0; i < n; i++)
01090 {
01091
01092
01093
01094 if (iy[i] > 0)
01095 {
01096 a[i] = ((qk[i] + low) * iy[i]) * 0.5;
01097 b[i] = ((up + qk[i]) * iy[i]) * 0.5;
01098 }
01099 else
01100 {
01101 b[i] = ((qk[i] + low) * iy[i]) * 0.5;
01102 a[i] = ((up + qk[i]) * iy[i]) * 0.5;
01103 }
01104 newdia[i] = (iy[i]*iy[i]);
01105 }
01106
01107 xmin = -1e33;
01108 xmax = 1e33;
01109
01110
01111 for (i = 0; i < n; i++)
01112 {
01113 uv[i] = i;
01114 xint[i] = a[i];
01115 xint[n+i] = b[i];
01116 }
01117
01118 xmid = xmin;
01119 xint[2*n] = xmin;
01120 xint[2*n+1] = xmax;
01121 ts = 0.0;
01122 sw = 0.0;
01123 luv = n;
01124 lxint = 2*n+2;
01125
01126 iter = 0;
01127 do {
01128 for (i = 0; i < luv; i++)
01129 {
01130 at[i] = a[uv[i]];
01131 bt[i] = b[uv[i]];
01132 newdt[i] = newdia[uv[i]];
01133 }
01134
01135 xmold = xmid;
01136 xmid = quick_select(xint, lxint);
01137 if (xmold == xmid)
01138 xmid = xint[(int32_t)(ThRandPos % lxint)];
01139
01140 s = ts;
01141 s1 = sw;
01142 for (i = 0; i < luv; i++)
01143 {
01144 if (xmid > bt[i])
01145 s += newdt[i]*bt[i];
01146 else if (xmid < at[i])
01147 s += newdt[i]*at[i];
01148 else
01149 s1 += newdt[i];
01150 }
01151
01152 testsum = s + s1*xmid;
01153 if (testsum <= (d+(1e-15)))
01154 xmin = xmid;
01155 if (testsum >= (d-(1e-15)))
01156 xmax = xmid;
01157
01158 l = 0;
01159 for (i = 0; i < lxint; i++)
01160 if((xint[i] >= xmin) && (xint[i] <= xmax))
01161 xint2[l++] = xint[i];
01162 lxint = l;
01163 memcpy(xint, xint2, lxint*sizeof(float64_t));
01164
01165 l = 0;
01166 for (i = 0; i < luv; i++)
01167 {
01168 if (xmin >= bt[i])
01169 ts += newdt[i]*bt[i];
01170 else if (xmax <= at[i])
01171 ts += newdt[i]*at[i];
01172 else if ((at[i] <= xmin) && (bt[i] >= xmax))
01173 sw += newdt[i];
01174 else
01175 uvt[l++] = uv[i];
01176 }
01177 luv = l;
01178 memcpy(uv, uvt, luv*sizeof(int32_t));
01179 iter++;
01180 } while(luv != 0 && iter < maxprojections);
01181
01182 if (sw == 0)
01183 xx = xmin;
01184 else
01185 xx = (d-ts) / sw;
01186
01187 for (i = 0; i < n; i++)
01188 {
01189 if (b[i] <= xmin)
01190 x[i] = b[i];
01191 else if (a[i] >= xmax)
01192 x[i] = a[i];
01193 else if ((a[i]<=xmin) && (xmax<=b[i]))
01194 x[i] = xx;
01195 else
01196 SG_SWARNING("Inner solver troubles...\n");
01197 }
01198
01199 for (i = 0; i < n; i++)
01200 x[i] = (2.0*x[i]*iy[i]-qk[i]);
01201
01202 #ifndef VARIABLES_ON_STACK
01203 SG_FREE(newdt);
01204 SG_FREE(newdia);
01205 SG_FREE(a);
01206 SG_FREE(b);
01207 SG_FREE(uvt);
01208 SG_FREE(uv);
01209 SG_FREE(bt);
01210 SG_FREE(at);
01211 SG_FREE(xint2);
01212 SG_FREE(xint);
01213 #endif
01214
01215 return(iter);
01216 }
01217
01218
01219
01220
01221 int32_t InnerProjector(
01222 int32_t method, int32_t n, int32_t *iy, float64_t e, float64_t *qk,
01223 float64_t l, float64_t u, float64_t *x, float64_t &lambda)
01224 {
01225 if (method == 0)
01226 return Pardalos(n, iy, e, qk, l, u, x);
01227 else
01228 return ProjectDai(n, iy, e, qk, l, u, x, lambda);
01229 }
01230 }
01231
01232
01233