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/UWedgeSep.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/UWedge.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 CUWedgeSep::CUWedgeSep() : CICAConverter() 00026 { 00027 init(); 00028 } 00029 00030 void CUWedgeSep::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 CUWedgeSep::~CUWedgeSep() 00041 { 00042 } 00043 00044 void CUWedgeSep::set_tau(SGVector<float64_t> tau) 00045 { 00046 m_tau = tau; 00047 } 00048 00049 SGVector<float64_t> CUWedgeSep::get_tau() const 00050 { 00051 return m_tau; 00052 } 00053 00054 SGNDArray<float64_t> CUWedgeSep::get_covs() const 00055 { 00056 return m_covs; 00057 } 00058 00059 CFeatures* CUWedgeSep::apply(CFeatures* features) 00060 { 00061 REQUIRE(features, "features is null"); 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 // Compute Correlation Matrices 00073 index_t * M_dims = SG_MALLOC(index_t, 3); 00074 M_dims[0] = n; 00075 M_dims[1] = n; 00076 M_dims[2] = N; 00077 m_covs = SGNDArray< float64_t >(M_dims, 3); 00078 00079 for (int t = 0; t < N; t++) 00080 { 00081 Map<MatrixXd> EM(m_covs.get_matrix(t),n,n); 00082 EM = cor(EX,m_tau[t]); 00083 } 00084 00085 // Diagonalize 00086 SGMatrix<float64_t> Q = CUWedge::diagonalize(m_covs, m_mixing_matrix, tol, max_iter); 00087 Map<MatrixXd> EQ(Q.matrix,n,n); 00088 00089 // Compute Mixing Matrix 00090 m_mixing_matrix = SGMatrix<float64_t>(n,n); 00091 Map<MatrixXd> C(m_mixing_matrix.matrix,n,n); 00092 C = EQ.inverse(); 00093 00094 // Normalize Estimated Mixing Matrix 00095 for (int t = 0; t < C.cols(); t++) 00096 C.col(t) /= C.col(t).maxCoeff(); 00097 00098 // Unmix 00099 EX = C.inverse() * EX; 00100 00101 return features; 00102 } 00103 00104 // Computing time delayed correlation matrix 00105 namespace 00106 { 00107 MatrixXd cor(MatrixXd x, int tau, bool mean_flag) 00108 { 00109 int m = x.rows(); 00110 int n = x.cols(); 00111 00112 // Center the data 00113 if ( mean_flag ) 00114 { 00115 VectorXd mean = x.rowwise().sum(); 00116 mean /= n; 00117 x = x.colwise() - mean; 00118 } 00119 00120 // Time-delayed Signal Matrix 00121 MatrixXd L = x.leftCols(n-tau); 00122 MatrixXd R = x.rightCols(n-tau); 00123 00124 // Compute Correlations 00125 MatrixXd K(m,m); 00126 K = (L * R.transpose()) / (n-tau); 00127 00128 // Symmetrize 00129 K = (K + K.transpose()) / 2.0; 00130 00131 return K; 00132 } 00133 }; 00134 00135 #endif // HAVE_EIGEN3