SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SOBI.cpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation