SHOGUN
v3.2.0
|
00001 /* 00002 * This program is free software; you can redistribute it and/or modify 00003 * it under the terms of the GNU General Public License as published by 00004 * the Free Software Foundation; either version 3 of the License, or 00005 * (at your option) any later version. 00006 * 00007 * Written (W) 2013 Kevin Hughes 00008 */ 00009 00010 #include <shogun/converter/ica/Jade.h> 00011 00012 #include <shogun/features/DenseFeatures.h> 00013 00014 #ifdef HAVE_EIGEN3 00015 00016 #include <shogun/mathematics/Math.h> 00017 #include <shogun/mathematics/eigen3.h> 00018 #include <shogun/mathematics/ajd/JADiagOrth.h> 00019 00020 #ifdef DEBUG_JADE 00021 #include <iostream> 00022 #endif 00023 00024 using namespace shogun; 00025 using namespace Eigen; 00026 00027 CJade::CJade() : CICAConverter() 00028 { 00029 init(); 00030 } 00031 00032 void CJade::init() 00033 { 00034 m_cumulant_matrix = SGMatrix<float64_t>(); 00035 SG_ADD(&m_cumulant_matrix, "cumulant_matrix", "m_cumulant_matrix", MS_NOT_AVAILABLE); 00036 } 00037 00038 CJade::~CJade() 00039 { 00040 } 00041 00042 SGMatrix<float64_t> CJade::get_cumulant_matrix() const 00043 { 00044 return m_cumulant_matrix; 00045 } 00046 00047 CFeatures* CJade::apply(CFeatures* features) 00048 { 00049 ASSERT(features); 00050 SG_REF(features); 00051 00052 SGMatrix<float64_t> X = ((CDenseFeatures<float64_t>*)features)->get_feature_matrix(); 00053 00054 int n = X.num_rows; 00055 int T = X.num_cols; 00056 int m = n; 00057 00058 Map<MatrixXd> EX(X.matrix,n,T); 00059 00060 // Mean center X 00061 VectorXd mean = (EX.rowwise().sum() / (float64_t)T); 00062 MatrixXd SPX = EX.colwise() - mean; 00063 00064 MatrixXd cov = (SPX * SPX.transpose()) / (float64_t)T; 00065 00066 #ifdef DEBUG_JADE 00067 std::cout << "cov" << std::endl; 00068 std::cout << cov << std::endl; 00069 #endif 00070 00071 // Whitening & Projection onto signal subspace 00072 SelfAdjointEigenSolver<MatrixXd> eig; 00073 eig.compute(cov); 00074 00075 #ifdef DEBUG_JADE 00076 std::cout << "eigenvectors" << std::endl; 00077 std::cout << eig.eigenvectors() << std::endl; 00078 00079 std::cout << "eigenvalues" << std::endl; 00080 std::cout << eig.eigenvalues().asDiagonal() << std::endl; 00081 #endif 00082 00083 // Scaling 00084 VectorXd scales = eig.eigenvalues().cwiseSqrt(); 00085 MatrixXd B = scales.cwiseInverse().asDiagonal() * eig.eigenvectors().transpose(); 00086 00087 #ifdef DEBUG_JADE 00088 std::cout << "whitener" << std::endl; 00089 std::cout << B << std::endl; 00090 #endif 00091 00092 // Sphering 00093 SPX = B * SPX; 00094 00095 // Estimation of the cumulant matrices 00096 int dimsymm = (m * ( m + 1)) / 2; // Dim. of the space of real symm matrices 00097 int nbcm = dimsymm; // number of cumulant matrices 00098 m_cumulant_matrix = SGMatrix<float64_t>(m,m*nbcm); // Storage for cumulant matrices 00099 Map<MatrixXd> CM(m_cumulant_matrix.matrix,m,m*nbcm); 00100 MatrixXd R(m,m); R.setIdentity(); 00101 MatrixXd Qij = MatrixXd::Zero(m,m); // Temp for a cum. matrix 00102 VectorXd Xim = VectorXd::Zero(m); // Temp 00103 VectorXd Xjm = VectorXd::Zero(m); // Temp 00104 VectorXd Xijm = VectorXd::Zero(m); // Temp 00105 int Range = 0; 00106 00107 for (int im = 0; im < m; im++) 00108 { 00109 Xim = SPX.row(im); 00110 Xijm = Xim.cwiseProduct(Xim); 00111 Qij = SPX.cwiseProduct(Xijm.replicate(1,m).transpose()) * SPX.transpose() / (float)T - R - 2*R.col(im)*R.col(im).transpose(); 00112 CM.block(0,Range,m,m) = Qij; 00113 Range = Range + m; 00114 for (int jm = 0; jm < im; jm++) 00115 { 00116 Xjm = SPX.row(jm); 00117 Xijm = Xim.cwiseProduct(Xjm); 00118 Qij = SPX.cwiseProduct(Xijm.replicate(1,m).transpose()) * SPX.transpose() / (float)T - R.col(im)*R.col(jm).transpose() - R.col(jm)*R.col(im).transpose(); 00119 CM.block(0,Range,m,m) = sqrt(2)*Qij; 00120 Range = Range + m; 00121 } 00122 } 00123 00124 #ifdef DEBUG_JADE 00125 std::cout << "cumulatant matrices" << std::endl; 00126 std::cout << CM << std::endl; 00127 #endif 00128 00129 // Stack cumulant matrix into ND Array 00130 index_t * M_dims = SG_MALLOC(index_t, 3); 00131 M_dims[0] = m; 00132 M_dims[1] = m; 00133 M_dims[2] = nbcm; 00134 SGNDArray< float64_t > M(M_dims, 3); 00135 00136 for (int i = 0; i < nbcm; i++) 00137 { 00138 Map<MatrixXd> EM(M.get_matrix(i),m,m); 00139 EM = CM.block(0,i*m,m,m); 00140 } 00141 00142 // Diagonalize 00143 SGMatrix<float64_t> Q = CJADiagOrth::diagonalize(M, m_mixing_matrix, tol, max_iter); 00144 Map<MatrixXd> EQ(Q.matrix,m,m); 00145 EQ = -1 * EQ.inverse(); 00146 00147 #ifdef DEBUG_JADE 00148 std::cout << "diagonalizer" << std::endl; 00149 std::cout << EQ << std::endl; 00150 #endif 00151 00152 // Separating matrix 00153 SGMatrix<float64_t> sep_matrix = SGMatrix<float64_t>(m,m); 00154 Map<MatrixXd> C(sep_matrix.matrix,m,m); 00155 C = EQ.transpose() * B; 00156 00157 // Sort 00158 VectorXd A = (B.inverse()*EQ).cwiseAbs2().colwise().sum(); 00159 bool swap = false; 00160 do 00161 { 00162 swap = false; 00163 for (int j = 1; j < n; j++) 00164 { 00165 if ( A(j) < A(j-1) ) 00166 { 00167 std::swap(A(j),A(j-1)); 00168 C.col(j).swap(C.col(j-1)); 00169 swap = true; 00170 } 00171 } 00172 00173 } while(swap); 00174 00175 for (int j = 0; j < m/2; j++) 00176 C.row(j).swap( C.row(m-1-j) ); 00177 00178 // Fix Signs 00179 VectorXd signs = VectorXd::Zero(m); 00180 for (int i = 0; i < m; i++) 00181 { 00182 if( C(i,0) < 0 ) 00183 signs(i) = -1; 00184 else 00185 signs(i) = 1; 00186 } 00187 C = signs.asDiagonal() * C; 00188 00189 #ifdef DEBUG_JADE 00190 std::cout << "un mixing matrix" << std::endl; 00191 std::cout << C << std::endl; 00192 #endif 00193 00194 // Unmix 00195 EX = C * EX; 00196 00197 m_mixing_matrix = SGMatrix<float64_t>(m,m); 00198 Map<MatrixXd> Emixing_matrix(m_mixing_matrix.matrix,m,m); 00199 Emixing_matrix = C.inverse(); 00200 00201 return features; 00202 } 00203 00204 #endif // HAVE_EIGEN3