SHOGUN
v3.2.0
|
00001 #ifdef HAVE_EIGEN3 00002 00003 #include <shogun/mathematics/ajd/JediDiag.h> 00004 00005 #include <shogun/base/init.h> 00006 00007 #include <shogun/mathematics/Math.h> 00008 #include <shogun/mathematics/eigen3.h> 00009 00010 using namespace shogun; 00011 using namespace Eigen; 00012 00013 void sweepjedi(float64_t *C, int *pMatSize, int *pMatNumber, 00014 float64_t *s_max, float64_t *sh_max, float64_t *A); 00015 00016 void iterJDI(float64_t *C, int *pMatSize, int *pMatNumber, int *ptn,int *ptm, 00017 float64_t *s_max, float64_t *sh_max, float64_t *A); 00018 00019 SGMatrix<float64_t> CJediDiag::diagonalize(SGNDArray<float64_t> C, SGMatrix<float64_t> V0, 00020 double eps, int itermax) 00021 { 00022 int d = C.dims[0]; 00023 int L = C.dims[2]; 00024 00025 SGMatrix<float64_t> V; 00026 if (V0.num_rows == d && V0.num_cols == d) 00027 V = V0.clone(); 00028 else 00029 V = SGMatrix<float64_t>::create_identity_matrix(d,1); 00030 00031 int iter = 0; 00032 float64_t sh_max = 1; 00033 float64_t s_max = 1; 00034 while (iter < itermax && ( (sh_max>eps) || (s_max>eps) )) 00035 { 00036 sh_max = 0; 00037 s_max = 0; 00038 sweepjedi(C.array, 00039 &d, &L, 00040 &s_max, &sh_max, 00041 V.matrix); 00042 iter++; 00043 } 00044 00045 if (iter == itermax) 00046 SG_SERROR("Convergence not reached\n") 00047 00048 Map<MatrixXd> EV(V.matrix,d,d); 00049 EV = EV.inverse(); 00050 00051 return V; 00052 } 00053 00054 void sweepjedi(float64_t *C, int *pMatSize, int *pMatNumber, 00055 float64_t *s_max, float64_t *sh_max, float64_t *A) 00056 { 00057 for (int n=2;n<=*pMatSize;n++) 00058 for (int m=1;m<n;m++) 00059 iterJDI(C, pMatSize, pMatNumber, &n, &m, s_max, sh_max, A); 00060 00061 int MS=*pMatSize; 00062 int MN=*pMatNumber; 00063 float64_t col_norm[MS]; 00064 00065 for (int i=0;i<MS;i++) 00066 { 00067 col_norm[i] = 0; 00068 00069 for (int j=0;j<MS;j++) 00070 col_norm[i] += A[i*MS+j]*A[i*MS+j]; 00071 00072 col_norm[i] = sqrt(col_norm[i]); 00073 } 00074 00075 float64_t daux=1; 00076 float64_t d[MS]; 00077 00078 for (int i=0;i<MS;i++) 00079 daux*=col_norm[i]; 00080 00081 daux=pow((float64_t)daux,(float64_t) 1/MS); 00082 00083 for (int i=0;i<MS;i++) 00084 d[i]=daux/col_norm[i]; 00085 00086 for (int i=0;i<MS;i++) 00087 for (int j=0;j<MS;j++) 00088 A[j*MS+i] *= d[j]; 00089 00090 for (int k=0;k<MN;k++) 00091 { 00092 for (int i=0;i<MS;i++) 00093 { 00094 for (int j=0;j<MS;j++) 00095 C[k*MS*MS+i*MS+j] *= (1/d[i])*(1/d[j]); 00096 } 00097 } 00098 } 00099 00100 void iterJDI(float64_t *C, int *pMatSize, int *pMatNumber, int *ptn,int *ptm, 00101 float64_t *s_max, float64_t *sh_max, float64_t *A) 00102 { 00103 int n=*ptn-1; 00104 int m=*ptm-1; 00105 int MN=*pMatNumber; 00106 int MS=*pMatSize; 00107 00108 float64_t tmm[MN]; 00109 float64_t tnn[MN]; 00110 float64_t tmn[MN]; 00111 for (int i = 0, d3 = 0; i < MN; i++, d3+=MS*MS) 00112 { 00113 tmm[i] = C[d3+m*MS+m]; 00114 tnn[i] = C[d3+n*MS+n]; 00115 tmn[i] = C[d3+n*MS+m]; 00116 } 00117 00118 // here we evd 00119 float64_t G[MN][3]; 00120 float64_t evectors[9], evalues[3]; 00121 for (int i = 0; i < MN; i++) 00122 { 00123 G[i][0] = 0.5*(tmm[i]+tnn[i]); 00124 G[i][1] = 0.5*(tmm[i]-tnn[i]); 00125 G[i][2] = tmn[i]; 00126 } 00127 float64_t GG[9]; 00128 for (int i = 0; i < 3; i++) 00129 { 00130 for (int j = 0; j <= i; j++) 00131 { 00132 GG[3*j+i] = 0; 00133 00134 for (int k = 0; k < MN; k++) 00135 GG[3*j+i] += G[k][i]*G[k][j]; 00136 00137 GG[3*i+j] = GG[3*j+i]; 00138 } 00139 } 00140 00141 for (int i = 0; i < 3; i++) 00142 GG[3*i] *= -1; 00143 00144 Map<MatrixXd> EGG(GG,3,3); 00145 00146 EigenSolver<MatrixXd> eig; 00147 eig.compute(EGG); 00148 00149 MatrixXd eigenvectors = eig.pseudoEigenvectors(); 00150 VectorXd eigenvalues = eig.pseudoEigenvalueMatrix().diagonal(); 00151 00152 eigenvectors.col(0) = eigenvectors.col(0).normalized(); 00153 eigenvectors.col(1) = eigenvectors.col(1).normalized(); 00154 eigenvectors.col(2) = eigenvectors.col(2).normalized(); 00155 00156 memcpy(evectors, eigenvectors.data(), 9*sizeof(float64_t)); 00157 memcpy(evalues, eigenvalues.data(), 3*sizeof(float64_t)); 00158 00159 float64_t tmp_evec[3],tmp_eval; 00160 if (fabs(evalues[1])<fabs(evalues[2])) 00161 { 00162 tmp_eval = evalues[1]; 00163 evalues[1] = evalues[2]; 00164 evalues[2] = tmp_eval; 00165 memcpy(tmp_evec,&evectors[3],3*sizeof(float64_t)); 00166 memcpy(&evectors[3],&evectors[6],3*sizeof(float64_t)); 00167 memcpy(&evectors[6],tmp_evec,3*sizeof(float64_t)); 00168 } 00169 if (fabs(evalues[0])<fabs(evalues[1])) 00170 { 00171 tmp_eval = evalues[0]; 00172 evalues[0] = evalues[1]; 00173 evalues[1] = tmp_eval; 00174 memcpy(tmp_evec,evectors,3*sizeof(float64_t)); 00175 memcpy(evectors,&evectors[3],3*sizeof(float64_t)); 00176 memcpy(&evectors[3],tmp_evec,3*sizeof(float64_t)); 00177 } 00178 if (fabs(evalues[1])<fabs(evalues[2])) 00179 { 00180 tmp_eval = evalues[1]; 00181 evalues[1] = evalues[2]; 00182 evalues[2] = tmp_eval; 00183 memcpy(tmp_evec,&evectors[3],3*sizeof(float64_t)); 00184 memcpy(&evectors[3],&evectors[6],3*sizeof(float64_t)); 00185 memcpy(&evectors[6],tmp_evec,3*sizeof(float64_t)); 00186 } 00187 00188 float64_t aux[9]; 00189 float64_t diagaux[3]; 00190 float64_t v[3]; 00191 for (int i = 0; i < 9; i++) 00192 aux[i] = evectors[i]; 00193 00194 for (int i = 0; i < 3; i++) 00195 aux[3*i] *= -1; 00196 00197 for (int i = 0; i < 3; i++) 00198 { 00199 diagaux[i] = 0; 00200 00201 for (int j = 0; j < 3; j++) 00202 diagaux[i] += evectors[3*i+j] * aux[3*i+j]; 00203 00204 diagaux[i] = 1/sqrt(fabs(diagaux[i])); 00205 } 00206 00207 int ind_min=2; 00208 00209 for (int i = 0; i < 3; i++) 00210 v[i] = evectors[3*ind_min+i]*diagaux[ind_min]; 00211 00212 if (v[2]<0) 00213 for (int i = 0; i < 3; i++) 00214 v[i]*=-1; 00215 00216 float64_t ch = sqrt( (1+sqrt(1+v[0]*v[0]))/2 ); 00217 float64_t sh = v[0]/(2*ch); 00218 float64_t c = sqrt( ( 1 + v[2]/sqrt(1+v[0]*v[0]) )/2 ); 00219 float64_t s = -v[1]/( 2*c*sqrt( (1+v[0]*v[0]) ) ); 00220 *s_max = std::max(*s_max,fabs(s)); 00221 *sh_max = std::max(*sh_max,fabs(sh)); 00222 00223 float64_t rm11=c*ch - s*sh; 00224 float64_t rm12=c*sh - s*ch; 00225 float64_t rm21=c*sh + s*ch; 00226 float64_t rm22=c*ch + s*sh; 00227 00228 float64_t h_slice1,h_slice2; 00229 float64_t buf[MS][MN]; 00230 00231 for (int i = 0; i < MS ;i++) 00232 { 00233 for (int j = 0, d3 = 0; j < MN; j++, d3+=MS*MS) 00234 { 00235 h_slice1 = C[d3+i*MS+m]; 00236 h_slice2 = C[d3+i*MS+n]; 00237 buf[i][j] = rm11*h_slice1 + rm21*h_slice2; 00238 C[d3+i*MS+n] = rm12*h_slice1 + rm22*h_slice2; 00239 C[d3+i*MS+m] = buf[i][j]; 00240 } 00241 } 00242 00243 for (int i = 0; i < MS; i++) 00244 { 00245 for (int j = 0, d3 = 0; j < MN; j++, d3+=MS*MS) 00246 C[d3+m*MS+i] = buf[i][j]; 00247 } 00248 for (int i = 0; i < MS; i++) 00249 { 00250 for (int j = 0, d3 = 0; j < MN; j++, d3+=MS*MS) 00251 C[d3+n*MS+i] = C[d3+i*MS+n]; 00252 } 00253 for (int i = 0; i < MN; i++) 00254 { 00255 C[i*MS*MS+m*MS+m] = (rm11*rm11)*tmm[i]+(2*rm11*rm21)*tmn[i]+(rm21*rm21)*tnn[i]; 00256 C[i*MS*MS+n*MS+m] = (rm11*rm12)*tmm[i]+(rm11*rm22+rm12*rm21)*tmn[i]+(rm21*rm22)*tnn[i]; 00257 C[i*MS*MS+m*MS+n] = C[i*MS*MS+n*MS+m]; 00258 C[i*MS*MS+n*MS+n] = (rm12*rm12)*tmm[i]+(2*rm12*rm22)*tmn[i]+(rm22*rm22)*tnn[i]; 00259 } 00260 00261 float64_t col1; 00262 float64_t col2; 00263 00264 for (int i = 0; i < MS; i++) 00265 { 00266 col1 = A[m*MS+i]; 00267 col2 = A[n*MS+i]; 00268 A[m*MS+i] = rm22*col1 - rm12*col2; 00269 A[n*MS+i] = -rm21*col1 + rm11*col2; 00270 } 00271 } 00272 #endif //HAVE_EIGEN3