00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #include <stdlib.h>
00018 #include <stdio.h>
00019 #include <time.h>
00020 #include <math.h>
00021 #include <shogun/lib/slep/q1/epph.h>
00022
00023 #define delta 1e-8
00024
00025 #define innerIter 1000
00026 #define outerIter 1000
00027
00028 void eplb(double * x, double *root, int * steps, double * v,int n, double z, double lambda0)
00029 {
00030
00031 int i, j, flag=0;
00032 int rho_1, rho_2, rho, rho_T, rho_S;
00033 int V_i_b, V_i_e, V_i;
00034 double lambda_1, lambda_2, lambda_T, lambda_S, lambda;
00035 double s_1, s_2, s, s_T, s_S, v_max, temp;
00036 double f_lambda_1, f_lambda_2, f_lambda, f_lambda_T, f_lambda_S;
00037 int iter_step=0;
00038
00039
00040
00041
00042
00043 if (z< 0){
00044 printf("\n z should be nonnegative!");
00045 return;
00046 }
00047
00048 V_i=0;
00049 if (v[0] !=0){
00050 rho_1=1;
00051 s_1=x[V_i]=v_max=fabs(v[0]);
00052 V_i++;
00053 }
00054 else{
00055 rho_1=0;
00056 s_1=v_max=0;
00057 }
00058
00059 for (i=1;i<n; i++){
00060 if (v[i]!=0){
00061 x[V_i]=fabs(v[i]); s_1+= x[V_i]; rho_1++;
00062
00063 if (x[V_i] > v_max)
00064 v_max=x[V_i];
00065 V_i++;
00066 }
00067 }
00068
00069
00070 if (s_1 <= z){
00071 flag=1; lambda=0;
00072 for(i=0;i<n;i++){
00073 x[i]=v[i];
00074 }
00075 *root=lambda;
00076 *steps=iter_step;
00077 return;
00078 }
00079
00080 lambda_1=0; lambda_2=v_max;
00081 f_lambda_1=s_1 -z;
00082
00083 rho_2=0; s_2=0; f_lambda_2=-z;
00084 V_i_b=0; V_i_e=V_i-1;
00085
00086 lambda=lambda0;
00087 if ( (lambda<lambda_2) && (lambda> lambda_1) ){
00088
00089
00090
00091
00092
00093 i=V_i_b; j=V_i_e; rho=0; s=0;
00094 while (i <= j){
00095 while( (i <= V_i_e) && (x[i] <= lambda) ){
00096 i++;
00097 }
00098 while( (j>=V_i_b) && (x[j] > lambda) ){
00099 s+=x[j];
00100 j--;
00101 }
00102 if (i<j){
00103 s+=x[i];
00104
00105 temp=x[i]; x[i]=x[j]; x[j]=temp;
00106 i++; j--;
00107 }
00108 }
00109
00110 rho=V_i_e-j; rho+=rho_2; s+=s_2;
00111 f_lambda=s-rho*lambda-z;
00112
00113 if ( fabs(f_lambda)< delta ){
00114 flag=1;
00115 }
00116
00117 if (f_lambda <0){
00118 lambda_2=lambda; s_2=s; rho_2=rho; f_lambda_2=f_lambda;
00119
00120 V_i_e=j; V_i=V_i_e-V_i_b+1;
00121 }
00122 else{
00123 lambda_1=lambda; rho_1=rho; s_1=s; f_lambda_1=f_lambda;
00124
00125 V_i_b=i; V_i=V_i_e-V_i_b+1;
00126 }
00127
00128 if (V_i==0){
00129
00130
00131
00132
00133 lambda=(s - z)/ rho;
00134 flag=1;
00135 }
00136
00137
00138
00139
00140
00141 }
00142
00143 while (!flag){
00144 iter_step++;
00145
00146
00147 lambda_T=lambda_1 + f_lambda_1 /rho_1;
00148 if(rho_2 !=0){
00149 if (lambda_2 + f_lambda_2 /rho_2 > lambda_T)
00150 lambda_T=lambda_2 + f_lambda_2 /rho_2;
00151 }
00152
00153
00154 lambda_S=lambda_2 - f_lambda_2 *(lambda_2-lambda_1)/(f_lambda_2-f_lambda_1);
00155
00156 if (fabs(lambda_T-lambda_S) <= delta){
00157 lambda=lambda_T; flag=1;
00158 break;
00159 }
00160
00161
00162 lambda=(lambda_T+lambda_S)/2;
00163
00164 s_T=s_S=s=0;
00165 rho_T=rho_S=rho=0;
00166 i=V_i_b; j=V_i_e;
00167 while (i <= j){
00168 while( (i <= V_i_e) && (x[i] <= lambda) ){
00169 if (x[i]> lambda_T){
00170 s_T+=x[i]; rho_T++;
00171 }
00172 i++;
00173 }
00174 while( (j>=V_i_b) && (x[j] > lambda) ){
00175 if (x[j] > lambda_S){
00176 s_S+=x[j]; rho_S++;
00177 }
00178 else{
00179 s+=x[j]; rho++;
00180 }
00181 j--;
00182 }
00183 if (i<j){
00184 if (x[i] > lambda_S){
00185 s_S+=x[i]; rho_S++;
00186 }
00187 else{
00188 s+=x[i]; rho++;
00189 }
00190
00191 if (x[j]> lambda_T){
00192 s_T+=x[j]; rho_T++;
00193 }
00194
00195 temp=x[i]; x[i]=x[j]; x[j]=temp;
00196 i++; j--;
00197 }
00198 }
00199
00200 s_S+=s_2; rho_S+=rho_2;
00201 s+=s_S; rho+=rho_S;
00202 s_T+=s; rho_T+=rho;
00203 f_lambda_S=s_S-rho_S*lambda_S-z;
00204 f_lambda=s-rho*lambda-z;
00205 f_lambda_T=s_T-rho_T*lambda_T-z;
00206
00207
00208
00209 if ( fabs(f_lambda)< delta ){
00210
00211 flag=1;
00212 break;
00213 }
00214 if ( fabs(f_lambda_S)< delta ){
00215
00216 lambda=lambda_S; flag=1;
00217 break;
00218 }
00219 if ( fabs(f_lambda_T)< delta ){
00220
00221 lambda=lambda_T; flag=1;
00222 break;
00223 }
00224
00225
00226
00227
00228
00229
00230
00231 if (f_lambda <0){
00232 lambda_2=lambda; s_2=s; rho_2=rho;
00233 f_lambda_2=f_lambda;
00234
00235 lambda_1=lambda_T; s_1=s_T; rho_1=rho_T;
00236 f_lambda_1=f_lambda_T;
00237
00238 V_i_e=j; i=V_i_b;
00239 while (i <= j){
00240 while( (i <= V_i_e) && (x[i] <= lambda_T) ){
00241 i++;
00242 }
00243 while( (j>=V_i_b) && (x[j] > lambda_T) ){
00244 j--;
00245 }
00246 if (i<j){
00247 x[j]=x[i];
00248 i++; j--;
00249 }
00250 }
00251 V_i_b=i; V_i=V_i_e-V_i_b+1;
00252 }
00253 else{
00254 lambda_1=lambda; s_1=s; rho_1=rho;
00255 f_lambda_1=f_lambda;
00256
00257 lambda_2=lambda_S; s_2=s_S; rho_2=rho_S;
00258 f_lambda_2=f_lambda_S;
00259
00260 V_i_b=i; j=V_i_e;
00261 while (i <= j){
00262 while( (i <= V_i_e) && (x[i] <= lambda_S) ){
00263 i++;
00264 }
00265 while( (j>=V_i_b) && (x[j] > lambda_S) ){
00266 j--;
00267 }
00268 if (i<j){
00269 x[i]=x[j];
00270 i++; j--;
00271 }
00272 }
00273 V_i_e=j; V_i=V_i_e-V_i_b+1;
00274 }
00275
00276 if (V_i==0){
00277 lambda=(s - z)/ rho; flag=1;
00278
00279 break;
00280 }
00281 }
00282
00283
00284 for(i=0;i<n;i++){
00285 if (v[i] > lambda)
00286 x[i]=v[i]-lambda;
00287 else
00288 if (v[i]< -lambda)
00289 x[i]=v[i]+lambda;
00290 else
00291 x[i]=0;
00292 }
00293 *root=lambda;
00294 *steps=iter_step;
00295 }
00296
00297 void epp1(double *x, double *v, int n, double rho)
00298 {
00299 int i;
00300
00301
00302
00303
00304
00305 for(i=0;i<n;i++){
00306 if (fabs(v[i])<=rho)
00307 x[i]=0;
00308 else
00309 if (v[i]< -rho)
00310 x[i]=v[i]+rho;
00311 else
00312 x[i]=v[i]-rho;
00313 }
00314 }
00315
00316 void epp2(double *x, double *v, int n, double rho)
00317 {
00318 int i;
00319 double v2=0, ratio;
00320
00321
00322
00323
00324
00325 for(i=0; i< n; i++){
00326 v2+=v[i]*v[i];
00327 }
00328 v2=sqrt(v2);
00329
00330 if (rho >= v2)
00331 for(i=0;i<n;i++)
00332 x[i]=0;
00333 else{
00334 ratio= (v2-rho) /v2;
00335 for(i=0;i<n;i++)
00336 x[i]=v[i]*ratio;
00337 }
00338 }
00339
00340 void eppInf(double *x, double * c, int * iter_step, double *v, int n, double rho, double c0)
00341 {
00342 int i, steps;
00343
00344
00345
00346
00347
00348 eplb(x, c, &steps, v, n, rho, c0);
00349
00350 for(i=0; i< n; i++){
00351 x[i]=v[i]-x[i];
00352 }
00353 iter_step[0]=steps;
00354 iter_step[1]=0;
00355 }
00356
00357 void zerofind(double *root, int * iterStep, double v, double p, double c, double x0)
00358 {
00359 double x, f, fprime, p1=p-1, pp;
00360 int step=0;
00361
00362
00363 if (v==0){
00364 *root=0; *iterStep=0; return;
00365 }
00366
00367 if (c==0){
00368 *root=v; * iterStep=0; return;
00369 }
00370
00371
00372 if ( (x0 <v) && (x0>0) )
00373 x=x0;
00374 else
00375 x=v;
00376
00377
00378 pp=pow(x, p1);
00379 f= x + c* pp -v;
00380
00381
00382
00383
00384
00385 while (1){
00386 step++;
00387
00388 fprime=1 + c* p1 * pp / x;
00389
00390
00391
00392
00393 x = x- f/fprime;
00394
00395
00396
00397
00398
00399 if (p>2){
00400 if (x>v){
00401 x=v;
00402 }
00403 }
00404 else{
00405 if ( (x<0) || (x>v)){
00406 x=1e-30;
00407
00408 f= x+c* pow(x,p1)-v;
00409
00410 if (f>0){
00411
00412
00413
00414
00415
00416 *root=x;
00417 * iterStep=step;
00418
00419 break;
00420 }
00421 }
00422 }
00423
00424
00425
00426
00427 pp=pow(x, p1);
00428 f= x + c* pp -v;
00429
00430
00431
00432
00433 if ( fabs(f) <= delta){
00434 *root=x;
00435 * iterStep=step;
00436 break;
00437 }
00438
00439 if (step>=innerIter){
00440 printf("\n The number of steps exceed %d, in finding the root for f(x)= x + c x^{p-1} - v, 0< x< v.", innerIter);
00441 printf("\n If you meet with this problem, please contact Jun Liu (j.liu@asu.edu). Thanks!");
00442 return;
00443 }
00444
00445 }
00446
00447
00448
00449
00450 }
00451
00452 double norm(double * v, double p, int n)
00453 {
00454 int i;
00455 double t=0;
00456
00457
00458
00459
00460
00461
00462
00463 for(i=0;i<n;i++)
00464 t+=pow(v[i], p);
00465
00466 return( pow(t, 1/p) );
00467 };
00468
00469 void eppO(double *x, double * cc, int * iter_step, double *v, int n, double rho, double p)
00470 {
00471 int i, *flag, bisStep, newtonStep=0, totoalStep=0;
00472 double vq=0, epsilon, vmax=0, vmin=1e10;
00473 double q=1/(1-1/p), c, c1, c2, root, f, xp;
00474
00475 double x_diff=0;
00476 double temp;
00477 int p_n=1;
00478
00479 flag=(int *)malloc(sizeof(int)*n);
00480
00481
00482
00483
00484
00485
00486
00487
00488 for(i=0; i< n; i++){
00489
00490 x[i]=0;
00491
00492 if (v[i]==0)
00493 flag[i]=0;
00494 else
00495 {
00496 if (v[i]>0)
00497 flag[i]=0;
00498 else
00499 {
00500 flag[i]=1;
00501 v[i]=-v[i];
00502
00503
00504 }
00505
00506 vq+=pow(v[i], q);
00507
00508
00509 if (v[i]>vmax)
00510 vmax=v[i];
00511
00512 if (v[i]<vmin)
00513 vmin=v[i];
00514 }
00515 }
00516 vq=pow(vq, 1/q);
00517
00518
00519
00520
00521 if (rho >= vq){
00522 *cc=0;
00523 iter_step[0]=iter_step[1]=0;
00524
00525
00526 for(i=0;i<n;i++){
00527 if (flag[i])
00528 v[i]=-v[i];
00529 }
00530
00531 free(flag);
00532 return;
00533 }
00534
00535
00536
00537
00538
00539 epsilon=(vq -rho)/ vq;
00540 if (p>2){
00541
00542 if ( log((1-epsilon) * vmin) - (p-1) * log( epsilon* vmin ) >= 709 )
00543 {
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553 for(i=0;i<n;i++){
00554 if (flag[i])
00555 v[i]=-v[i];
00556 }
00557
00558 eppInf(x, cc, iter_step, v, n, rho, 0);
00559
00560 free(flag);
00561 return;
00562 }
00563
00564 c1= (1-epsilon) * vmax / pow(epsilon* vmax, p-1);
00565 c2= (1-epsilon) * vmin / pow(epsilon* vmin, p-1);
00566 }
00567 else{
00568
00569 c2= (1-epsilon) * vmax / pow(epsilon* vmax, p-1);
00570 c1= (1-epsilon) * vmin / pow(epsilon* vmin, p-1);
00571 }
00572
00573
00574
00575
00576
00577
00578 if (fabs(c1-c2) <= delta){
00579 c=c1;
00580 }
00581 else
00582 c=(c1+c2)/2;
00583
00584
00585 bisStep =0;
00586
00587 while(1){
00588 bisStep++;
00589
00590
00591 x_diff=0;
00592 for(i=0;i<n;i++){
00593 zerofind(&root, &newtonStep, v[i], p, c, x[i]);
00594
00595 temp=fabs(root-x[i]);
00596 if (x_diff< temp )
00597 x_diff=temp;
00598
00599 x[i]=root;
00600 totoalStep+=newtonStep;
00601 }
00602
00603 xp=norm(x, p, n);
00604
00605 f=rho * pow(xp, 1-p) - c;
00606
00607 if ( fabs(f)<=delta || fabs(c1-c2)<=delta )
00608 break;
00609 else{
00610 if (f>0){
00611 if ( (x_diff <=delta) && (p_n==0) )
00612 break;
00613
00614 c1=c; p_n=1;
00615 }
00616 else{
00617
00618 if ( (x_diff <=delta) && (p_n==1) )
00619 break;
00620
00621 c2=c; p_n=0;
00622 }
00623 }
00624 c=(c1+c2)/2;
00625
00626 if (bisStep>=outerIter){
00627
00628
00629 if ( fabs(c1-c2) <=delta * c2 )
00630 break;
00631 else{
00632 printf("\n The number of bisection steps exceed %d.", outerIter);
00633 printf("\n c1=%e, c2=%e, x_diff=%e, f=%e",c1,c2,x_diff,f);
00634 printf("\n If you meet with this problem, please contact Jun Liu (j.liu@asu.edu). Thanks!");
00635
00636 return;
00637 }
00638 }
00639
00640
00641
00642
00643 }
00644
00645
00646
00647
00648
00649 for(i=0;i<n;i++){
00650 if (flag[i]){
00651 x[i]=-x[i];
00652 v[i]=-v[i];
00653 }
00654 }
00655 free(flag);
00656
00657 *cc=c;
00658
00659 iter_step[0]=bisStep;
00660 iter_step[1]=totoalStep;
00661 }
00662
00663 void epp(double *x, double * c, int * iter_step, double * v, int n, double rho, double p, double c0){
00664 if (rho <0){
00665 printf("\n rho should be non-negative!");
00666 exit(1);
00667 }
00668
00669 if (p==1){
00670 epp1(x, v, n, rho);
00671 *c=0;
00672 iter_step[0]=iter_step[1]=0;
00673 }
00674 else
00675 if (p==2){
00676 epp2(x, v, n, rho);
00677 *c=0;
00678 iter_step[0]=iter_step[1]=0;
00679 }
00680 else
00681 if (p>=1e6)
00682 eppInf(x, c, iter_step, v, n, rho, c0);
00683 else
00684 eppO(x, c, iter_step, v, n, rho, p);
00685 }
00686