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