SHOGUN
v3.2.0
|
00001 /* This program is free software: you can redistribute it and/or modify 00002 * it under the terms of the GNU General Public License as published by 00003 * the Free Software Foundation, either version 3 of the License, or 00004 * (at your option) any later version. 00005 * 00006 * This program is distributed in the hope that it will be useful, 00007 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00008 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00009 * GNU General Public License for more details. 00010 * 00011 * You should have received a copy of the GNU General Public License 00012 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00013 * 00014 * Copyright (C) 2009 - 2012 Jun Liu and Jieping Ye 00015 */ 00016 00017 #ifndef EPSP_SLEP 00018 #define EPSP_SLEP 00019 00020 #include <stdlib.h> 00021 #include <stdio.h> 00022 #include <time.h> 00023 #include <math.h> 00024 00025 #define delta 1e-12 00026 00027 /* 00028 Euclidean Projection onto the simplex (epsp) 00029 00030 min 1/2 ||x- y||_2^2 00031 s.t. ||x||_1 = z, x >=0 00032 00033 which is converted to the following zero finding problem 00034 00035 f(lambda)= sum( max( x-lambda,0) )-z=0 00036 00037 Usage: 00038 [x, lambda, iter_step]=epsp(y, n, z, lambda0); 00039 00040 */ 00041 00042 void epsp(double * x, double *root, int * stepsp, double * v, 00043 int n, double z, double lambda0) 00044 { 00045 00046 int i, j, flag=0; 00047 int rho_1, rho_2, rho, rho_T, rho_S; 00048 int V_i_b, V_i_e, V_i; 00049 double lambda_1, lambda_2, lambda_T, lambda_S, lambda; 00050 double s_1, s_2, s, s_T, s_S, v_max, temp; 00051 double f_lambda_1, f_lambda_2, f_lambda, f_lambda_T, f_lambda_S; 00052 int iter_step=0; 00053 00054 00055 00056 if (z< 0){ 00057 printf("\n z should be nonnegative!"); 00058 exit(-1); 00059 } 00060 00061 00062 /* 00063 * find the maximal value in v 00064 */ 00065 v_max=v[0]; 00066 for (i=1;i<n; i++){ 00067 if (v[i] > v_max) 00068 v_max=v[i]; 00069 } 00070 00071 00072 lambda_1=v_max - z; lambda_2=v_max; 00073 /* 00074 * copy v to x 00075 * compute f_lambda_1, rho_1, s_1 00076 */ 00077 V_i=0;s_1=0; rho_1=0; 00078 for (i=0;i<n;i++){ 00079 if (v[i] > lambda_1){ 00080 x[V_i]=v[i]; 00081 00082 s_1+=x[V_i]; rho_1++; 00083 00084 V_i++; 00085 } 00086 } 00087 f_lambda_1=s_1-rho_1* lambda_1 -z; 00088 00089 rho_2=0; s_2=0; f_lambda_2=-z; 00090 V_i_b=0; V_i_e=V_i-1; 00091 00092 lambda=lambda0; 00093 if ( (lambda<lambda_2) && (lambda> lambda_1) ){ 00094 /*------------------------------------------------------------------- 00095 Initialization with the root 00096 *------------------------------------------------------------------- 00097 */ 00098 00099 i=V_i_b; j=V_i_e; rho=0; s=0; 00100 while (i <= j){ 00101 while( (i <= V_i_e) && (x[i] <= lambda) ){ 00102 i++; 00103 } 00104 while( (j>=V_i_b) && (x[j] > lambda) ){ 00105 s+=x[j]; 00106 j--; 00107 } 00108 if (i<j){ 00109 s+=x[i]; 00110 00111 temp=x[i]; x[i]=x[j]; x[j]=temp; 00112 i++; j--; 00113 } 00114 } 00115 00116 rho=V_i_e-j; rho+=rho_2; s+=s_2; 00117 f_lambda=s-rho*lambda-z; 00118 00119 if ( fabs(f_lambda)< delta ){ 00120 flag=1; 00121 } 00122 00123 if (f_lambda <0){ 00124 lambda_2=lambda; s_2=s; rho_2=rho; f_lambda_2=f_lambda; 00125 00126 V_i_e=j; V_i=V_i_e-V_i_b+1; 00127 } 00128 else{ 00129 lambda_1=lambda; rho_1=rho; s_1=s; f_lambda_1=f_lambda; 00130 00131 V_i_b=i; V_i=V_i_e-V_i_b+1; 00132 } 00133 00134 if (V_i==0){ 00135 /*printf("\n rho=%d, rho_1=%d, rho_2=%d",rho, rho_1, rho_2);*/ 00136 00137 /*printf("\n V_i=%d",V_i);*/ 00138 00139 lambda=(s - z)/ rho; 00140 flag=1; 00141 } 00142 /*------------------------------------------------------------------- 00143 End of initialization 00144 *-------------------------------------------------------------------- 00145 */ 00146 00147 }/* end of if(!flag) */ 00148 00149 while (!flag){ 00150 iter_step++; 00151 00152 /* compute lambda_T */ 00153 lambda_T=lambda_1 + f_lambda_1 /rho_1; 00154 if(rho_2 !=0){ 00155 if (lambda_2 + f_lambda_2 /rho_2 > lambda_T) 00156 lambda_T=lambda_2 + f_lambda_2 /rho_2; 00157 } 00158 00159 /* compute lambda_S */ 00160 lambda_S=lambda_2 - f_lambda_2 *(lambda_2-lambda_1)/(f_lambda_2-f_lambda_1); 00161 00162 if (fabs(lambda_T-lambda_S) <= delta){ 00163 lambda=lambda_T; flag=1; 00164 break; 00165 } 00166 00167 /* set lambda as the middle point of lambda_T and lambda_S */ 00168 lambda=(lambda_T+lambda_S)/2; 00169 00170 s_T=s_S=s=0; 00171 rho_T=rho_S=rho=0; 00172 i=V_i_b; j=V_i_e; 00173 while (i <= j){ 00174 while( (i <= V_i_e) && (x[i] <= lambda) ){ 00175 if (x[i]> lambda_T){ 00176 s_T+=x[i]; rho_T++; 00177 } 00178 i++; 00179 } 00180 while( (j>=V_i_b) && (x[j] > lambda) ){ 00181 if (x[j] > lambda_S){ 00182 s_S+=x[j]; rho_S++; 00183 } 00184 else{ 00185 s+=x[j]; rho++; 00186 } 00187 j--; 00188 } 00189 if (i<j){ 00190 if (x[i] > lambda_S){ 00191 s_S+=x[i]; rho_S++; 00192 } 00193 else{ 00194 s+=x[i]; rho++; 00195 } 00196 00197 if (x[j]> lambda_T){ 00198 s_T+=x[j]; rho_T++; 00199 } 00200 00201 temp=x[i]; x[i]=x[j]; x[j]=temp; 00202 i++; j--; 00203 } 00204 } 00205 00206 s_S+=s_2; rho_S+=rho_2; 00207 s+=s_S; rho+=rho_S; 00208 s_T+=s; rho_T+=rho; 00209 f_lambda_S=s_S-rho_S*lambda_S-z; 00210 f_lambda=s-rho*lambda-z; 00211 f_lambda_T=s_T-rho_T*lambda_T-z; 00212 00213 /*printf("\n %d & %d & %5.6f & %5.6f & %5.6f & %5.6f & %5.6f \\\\ \n \\hline ", iter_step, V_i, lambda_1, lambda_T, lambda, lambda_S, lambda_2);*/ 00214 00215 if ( fabs(f_lambda)< delta ){ 00216 /*printf("\n lambda");*/ 00217 flag=1; 00218 break; 00219 } 00220 if ( fabs(f_lambda_S)< delta ){ 00221 /* printf("\n lambda_S");*/ 00222 lambda=lambda_S; flag=1; 00223 break; 00224 } 00225 if ( fabs(f_lambda_T)< delta ){ 00226 /* printf("\n lambda_T");*/ 00227 lambda=lambda_T; flag=1; 00228 break; 00229 } 00230 00231 /* 00232 printf("\n\n f_lambda_1=%5.6f, f_lambda_2=%5.6f, f_lambda=%5.6f",f_lambda_1,f_lambda_2, f_lambda); 00233 printf("\n lambda_1=%5.6f, lambda_2=%5.6f, lambda=%5.6f",lambda_1, lambda_2, lambda); 00234 printf("\n rho_1=%d, rho_2=%d, rho=%d ",rho_1, rho_2, rho); 00235 */ 00236 00237 if (f_lambda <0){ 00238 lambda_2=lambda; s_2=s; rho_2=rho; 00239 f_lambda_2=f_lambda; 00240 00241 lambda_1=lambda_T; s_1=s_T; rho_1=rho_T; 00242 f_lambda_1=f_lambda_T; 00243 00244 V_i_e=j; i=V_i_b; 00245 while (i <= j){ 00246 while( (i <= V_i_e) && (x[i] <= lambda_T) ){ 00247 i++; 00248 } 00249 while( (j>=V_i_b) && (x[j] > lambda_T) ){ 00250 j--; 00251 } 00252 if (i<j){ 00253 x[j]=x[i]; 00254 i++; j--; 00255 } 00256 } 00257 V_i_b=i; V_i=V_i_e-V_i_b+1; 00258 } 00259 else{ 00260 lambda_1=lambda; s_1=s; rho_1=rho; 00261 f_lambda_1=f_lambda; 00262 00263 lambda_2=lambda_S; s_2=s_S; rho_2=rho_S; 00264 f_lambda_2=f_lambda_S; 00265 00266 V_i_b=i; j=V_i_e; 00267 while (i <= j){ 00268 while( (i <= V_i_e) && (x[i] <= lambda_S) ){ 00269 i++; 00270 } 00271 while( (j>=V_i_b) && (x[j] > lambda_S) ){ 00272 j--; 00273 } 00274 if (i<j){ 00275 x[i]=x[j]; 00276 i++; j--; 00277 } 00278 } 00279 V_i_e=j; V_i=V_i_e-V_i_b+1; 00280 } 00281 00282 if (V_i==0){ 00283 lambda=(s - z)/ rho; flag=1; 00284 /*printf("\n V_i=0, lambda=%5.6f",lambda);*/ 00285 break; 00286 } 00287 }/* end of while */ 00288 00289 00290 for(i=0;i<n;i++){ 00291 if (v[i] > lambda) 00292 x[i]=v[i]-lambda; 00293 else 00294 x[i]=0; 00295 } 00296 *root=lambda; 00297 *stepsp=iter_step; 00298 } 00299 #endif /* ----- #ifndef EPSP_SLEP ----- */ 00300