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