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/SOBI.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 using namespace shogun; 00021 using namespace Eigen; 00022 00023 namespace { MatrixXd cor(MatrixXd x, int tau = 0, bool mean_flag = true); }; 00024 00025 CSOBI::CSOBI() : CICAConverter() 00026 { 00027 init(); 00028 } 00029 00030 void CSOBI::init() 00031 { 00032 m_tau = SGVector<float64_t>(4); 00033 m_tau[0]=0; m_tau[1]=1; m_tau[2]=2; m_tau[3]=3; 00034 00035 m_covs = SGNDArray<float64_t>(); 00036 00037 SG_ADD(&m_tau, "tau", "tau vector", MS_AVAILABLE); 00038 } 00039 00040 CSOBI::~CSOBI() 00041 { 00042 } 00043 00044 void CSOBI::set_tau(SGVector<float64_t> tau) 00045 { 00046 m_tau = tau; 00047 } 00048 00049 SGVector<float64_t> CSOBI::get_tau() const 00050 { 00051 return m_tau; 00052 } 00053 00054 SGNDArray<float64_t> CSOBI::get_covs() const 00055 { 00056 return m_covs; 00057 } 00058 00059 CFeatures* CSOBI::apply(CFeatures* features) 00060 { 00061 ASSERT(features); 00062 SG_REF(features); 00063 00064 SGMatrix<float64_t> X = ((CDenseFeatures<float64_t>*)features)->get_feature_matrix(); 00065 00066 int n = X.num_rows; 00067 int m = X.num_cols; 00068 int N = m_tau.vlen; 00069 00070 Map<MatrixXd> EX(X.matrix,n,m); 00071 00072 // Whitening or Sphering 00073 MatrixXd M0 = cor(EX,int(m_tau[0])); 00074 EigenSolver<MatrixXd> eig; 00075 eig.compute(M0); 00076 MatrixXd SPH = (eig.pseudoEigenvectors() * eig.pseudoEigenvalueMatrix().cwiseSqrt() * eig.pseudoEigenvectors ().transpose()).inverse(); 00077 MatrixXd spx = SPH*EX; 00078 00079 // Compute Correlation Matrices 00080 index_t * M_dims = SG_MALLOC(index_t, 3); 00081 M_dims[0] = n; 00082 M_dims[1] = n; 00083 M_dims[2] = N; 00084 m_covs = SGNDArray< float64_t >(M_dims, 3); 00085 00086 for(int t = 0; t < N; t++) 00087 { 00088 Map<MatrixXd> EM(m_covs.get_matrix(t),n,n); 00089 EM = cor(spx,m_tau[t]); 00090 } 00091 00092 // Diagonalize 00093 SGMatrix<float64_t> Q = CJADiagOrth::diagonalize(m_covs); 00094 Map<MatrixXd> EQ(Q.matrix,n,n); 00095 00096 // Compute Mixing Matrix 00097 m_mixing_matrix = SGMatrix<float64_t>(n,n); 00098 Map<MatrixXd> C(m_mixing_matrix.matrix,n,n); 00099 C = SPH.inverse() * EQ.transpose(); 00100 00101 // Normalize Estimated Mixing Matrix 00102 for(int t = 0; t < C.cols(); t++) 00103 { 00104 C.col(t) /= C.col(t).maxCoeff(); 00105 } 00106 00107 // Unmix 00108 EX = C.inverse() * EX; 00109 00110 return features; 00111 } 00112 00113 // Computing time delayed correlation matrix 00114 namespace 00115 { 00116 MatrixXd cor(MatrixXd x, int tau, bool mean_flag) 00117 { 00118 int m = x.rows(); 00119 int n = x.cols(); 00120 00121 // Center the data 00122 if ( mean_flag ) 00123 { 00124 VectorXd mean = x.rowwise().sum(); 00125 mean /= n; 00126 x = x.colwise() - mean; 00127 } 00128 00129 // Time-delayed Signal Matrix 00130 MatrixXd L = x.leftCols(n-tau); 00131 MatrixXd R = x.rightCols(n-tau); 00132 00133 // Compute Correlations 00134 MatrixXd K(m,m); 00135 K = (L * R.transpose()) / (n-tau); 00136 00137 // Symmetrize 00138 K = (K + K.transpose()) / 2.0; 00139 00140 return K; 00141 } 00142 }; 00143 #endif // HAVE_EIGEN3