Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #ifndef EPSGLASSO_SLEP
00018 #define EPSGLASSO_SLEP
00019
00020 #include <stdlib.h>
00021 #include <stdio.h>
00022 #include <time.h>
00023 #include <math.h>
00024 #include <shogun/lib/slep/q1/epph.h>
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 void epsgLasso(double *X, double * normx, double * V, int k, int n, double lambda1, double lambda2, int pflag){
00038 int i, j, *iter_step, nn=n*k, m;
00039 double *v, *x;
00040 double normValue,c0=0, c;
00041
00042 v=(double *)malloc(sizeof(double)*n);
00043 x=(double *)malloc(sizeof(double)*n);
00044 iter_step=(int *)malloc(sizeof(int)*2);
00045
00046
00047
00048
00049 normx[0]=normx[1]=0;
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066 for (i=0;i<nn;i++){
00067 if (V[i]< -lambda1)
00068 X[i]=V[i] + lambda1;
00069 else
00070 if (V[i]> lambda1)
00071 X[i]=V[i] - lambda1;
00072 else
00073 X[i]=0;
00074 }
00075
00076
00077
00078
00079
00080 if (pflag==2){
00081 for(i=0; i<k; i++){
00082
00083
00084
00085
00086 normValue=0;
00087
00088 m=n%5;
00089 for(j=0;j<m;j++){
00090 v[j]=X[i + j*k];
00091 }
00092 for(j=m;j<n;j+=5){
00093 v[j ]=X[i + j*k];
00094 v[j+1]=X[i + (j+1)*k ];
00095 v[j+2]=X[i + (j+2)*k];
00096 v[j+3]=X[i + (j+3)*k];
00097 v[j+4]=X[i + (j+4)*k];
00098 }
00099
00100 m=n%5;
00101 for(j=0;j<m;j++){
00102 normValue+=v[j]*v[j];
00103 }
00104 for(j=m;j<n;j+=5){
00105 normValue+=v[j]*v[j]+
00106 v[j+1]*v[j+1]+
00107 v[j+2]*v[j+2]+
00108 v[j+3]*v[j+3]+
00109 v[j+4]*v[j+4];
00110 }
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120 normValue=sqrt(normValue);
00121
00122 if (normValue<= lambda2){
00123 for(j=0; j<n; j++)
00124 X[i + j*k]=0;
00125
00126
00127 }
00128 else{
00129
00130 normx[1]+=normValue-lambda2;
00131
00132
00133 normValue=(normValue-lambda2)/normValue;
00134
00135 m=n%5;
00136 for(j=0;j<m;j++){
00137 X[i + j*k]*=normValue;
00138 normx[0]+=fabs(X[i + j*k]);
00139 }
00140 for(j=m; j<n;j+=5){
00141 X[i + j*k]*=normValue;
00142 X[i + (j+1)*k]*=normValue;
00143 X[i + (j+2)*k]*=normValue;
00144 X[i + (j+3)*k]*=normValue;
00145 X[i + (j+4)*k]*=normValue;
00146
00147 normx[0]+=fabs(X[i + j*k])+
00148 fabs(X[i + (j+1)*k])+
00149 fabs(X[i + (j+2)*k])+
00150 fabs(X[i + (j+3)*k])+
00151 fabs(X[i + (j+4)*k]);
00152 }
00153
00154
00155
00156
00157
00158 }
00159 }
00160 }
00161 else{
00162 for(i=0; i<k; i++){
00163
00164
00165
00166
00167 normValue=0;
00168 for(j=0; j<n; j++){
00169 v[j]=X[i + j*k];
00170
00171 normValue+=fabs(v[j]);
00172 }
00173
00174 if (normValue<= lambda2){
00175 for(j=0; j<n; j++)
00176 X[i + j*k]=0;
00177 }
00178 else{
00179 eplb(x, &c, iter_step, v, n, lambda2, c0);
00180
00181 for(j=0; j<n; j++){
00182 if (X[i + j*k] > c)
00183 X[i + j*k]=c;
00184 else
00185 if (X[i + j*k]<-c)
00186 X[i + j*k]=-c;
00187 }
00188 }
00189 }
00190 }
00191
00192
00193 free(v);
00194 free(x);
00195 free(iter_step);
00196 }
00197 #endif
00198