SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Jade.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/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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation