replaced uneccessary eigen3 typedefs with proper eigen3 type. fixed a…

…n off by 1 error in amari index. added jade example to example makefile. added links to shogun in ica notebook.

doc/ipython-notebooks/ica/bss_jade.ipynb

@@ -327,7 +327,7 @@
       "\n",
       "Cardoso, J. F., & Souloumiac, A. (1993). Blind beamforming for non-Gaussian signals. In IEE Proceedings F (Radar and Signal Processing) (Vol. 140, No. 6, pp. 362-370). IET Digital Library.\n",
       "\n",
-      "Shogun also has several other ICA algorithms including the Second Order Blind Identification (SOBI) algorithm, FFSep, JediSep, UWedgeSep and FastICA. All of the algorithms inherit from the ICAConverter base class and share some common methods for setting an intial guess for the mixing matrix, retrieving the final mixing matrix and getting/setting the number of iterations to run and the desired convergence tolerance. Some of the algorithms have additional getters for intermediate calculations, for example Jade has a method for returning the 4th order cumulant tensor while the \"Sep\" algorithms have a getter for the time lagged covariance matrices. Check out the source code on GitHub or the Shogun docs for more details!  "
+      "Shogun also has several other ICA algorithms including the Second Order Blind Identification (SOBI) algorithm, FFSep, JediSep, UWedgeSep and FastICA. All of the algorithms inherit from the ICAConverter base class and share some common methods for setting an intial guess for the mixing matrix, retrieving the final mixing matrix and getting/setting the number of iterations to run and the desired convergence tolerance. Some of the algorithms have additional getters for intermediate calculations, for example Jade has a method for returning the 4th order cumulant tensor while the \"Sep\" algorithms have a getter for the time lagged covariance matrices. Check out the source code on GitHub (https://github.com/shogun-toolbox/shogun) or the Shogun docs (http://www.shogun-toolbox.org/doc/en/latest/annotated.html) for more details!  "
      ]
     },
     {

examples/undocumented/libshogun/CMakeLists.txt

@@ -67,7 +67,8 @@ SET(EXAMPLES
 	converter_hessianlocallylinearembedding 
 	converter_kernellocallylinearembedding 
 	converter_multidimensionalscaling 
-	converter_isomap 
+	converter_isomap
+	converter_jade_bss  
 	converter_diffusionmaps 
 	converter_laplacianeigenmaps 
 	converter_neighborhoodpreservingembedding 

examples/undocumented/libshogun/Makefile

@@ -80,6 +80,7 @@ TARGETS = basic_minimal \
 		  converter_kernellocallylinearembedding \
 		  converter_multidimensionalscaling \
 		  converter_isomap \
+		  converter_jade_bss \
 		  converter_diffusionmaps \
 		  converter_laplacianeigenmaps \
 		  converter_neighborhoodpreservingembedding \

examples/undocumented/libshogun/converter_jade_bss.cpp

@@ -24,23 +24,19 @@
 #include <shogun/evaluation/ica/PermutationMatrix.h>
 #include <shogun/evaluation/ica/AmariIndex.h>
 
-using namespace Eigen;
-
-typedef Matrix< float64_t, Dynamic, Dynamic, ColMajor > EMatrix;
-typedef Matrix< float64_t, Dynamic, 1, ColMajor > EVector;
-
 using namespace shogun;
+using namespace Eigen;
 
 void test()
 {
 	// Generate sample data	
 	CMath::init_random(0); 
 	int n_samples = 2000;
-	EVector time(n_samples, true);
+	VectorXd time(n_samples, true);
 	time.setLinSpaced(n_samples,0,10);
 
 	// Source Signals
-	EMatrix S(2,n_samples);
+	MatrixXd S(2,n_samples);
 	for(int i = 0; i < n_samples; i++)
 	{
 		// Sin wave
@@ -53,14 +49,14 @@ void test()
 	}
 
 	// Standardize data
-	EVector avg = S.rowwise().sum() / n_samples;
-	EVector std = ((S.colwise() - avg).array().pow(2).rowwise().sum() / n_samples).array().sqrt();
+	VectorXd avg = S.rowwise().sum() / n_samples;
+	VectorXd std = ((S.colwise() - avg).array().pow(2).rowwise().sum() / n_samples).array().sqrt();
 	for(int i = 0; i < n_samples; i++)
 		S.col(i) = S.col(i).cwiseQuotient(std);
 
 	// Mixing Matrix
 	SGMatrix<float64_t> mixing_matrix(2,2);
-	Eigen::Map<EMatrix> A(mixing_matrix.matrix,2,2);
+	Map<MatrixXd> A(mixing_matrix.matrix,2,2);
 	A(0,0) = 1;    A(0,1) = 0.5;
 	A(1,0) = 0.5;  A(1,1) = 1;
 
@@ -69,7 +65,7 @@ void test()
 
 	// Mix signals
 	SGMatrix<float64_t> X(2,n_samples);
-	Eigen::Map<EMatrix> EX(X.matrix,2,n_samples);
+	Map<MatrixXd> EX(X.matrix,2,n_samples);
 	EX = A * S;
 	CDenseFeatures< float64_t >* mixed_signals = new CDenseFeatures< float64_t >(X);
 
@@ -81,13 +77,13 @@ void test()
 	SG_REF(signals);
 
 	// Close to a permutation matrix (with random scales)
-	Eigen::Map<EMatrix> EA(jade->get_mixing_matrix().matrix,2,2);
+	Map<MatrixXd> EA(jade->get_mixing_matrix().matrix,2,2);
 
 	std::cout << "Estimated Mixing Matrix:" << std::endl;
 	std::cout << EA << std::endl << std::endl;
 
 	SGMatrix<float64_t> P(2,2);
-	Eigen::Map<EMatrix> EP(P.matrix,2,2);
+	Eigen::Map<MatrixXd> EP(P.matrix,2,2);
 	EP = EA.inverse() * A;
 
 	bool isperm = is_permutation_matrix(P);

src/shogun/converter/ica/FFSep.cpp

@@ -17,14 +17,10 @@
 #include <shogun/mathematics/eigen3.h>
 #include <shogun/mathematics/ajd/FFDiag.h>
 
-using namespace Eigen;
-
-typedef Matrix< float64_t, Dynamic, 1, ColMajor > EVector;
-typedef Matrix< float64_t, Dynamic, Dynamic, ColMajor > EMatrix;
-
 using namespace shogun;
+using namespace Eigen;
 
-namespace { EMatrix cor(EMatrix x, int tau = 0, bool mean_flag = true); };
+namespace { MatrixXd cor(MatrixXd x, int tau = 0, bool mean_flag = true); };
 
 CFFSep::CFFSep() : CICAConverter()
 {		
@@ -71,7 +67,7 @@ CFeatures* CFFSep::apply(CFeatures* features)
 	int m = X.num_cols;
 	int N = m_tau.vlen;
 
-	Eigen::Map<EMatrix> EX(X.matrix,n,m);	
+	Map<MatrixXd> EX(X.matrix,n,m);	
 
 	// Compute Correlation Matrices
 	index_t * M_dims = SG_MALLOC(index_t, 3);
@@ -82,17 +78,17 @@ CFeatures* CFFSep::apply(CFeatures* features)
 
 	for (int t = 0; t < N; t++)
 	{
-		Eigen::Map<EMatrix> EM(m_covs.get_matrix(t),n,n);
+		Map<MatrixXd> EM(m_covs.get_matrix(t),n,n);
 		EM = cor(EX,m_tau[t]);
 	}
 
 	// Diagonalize
 	SGMatrix<float64_t> Q = CFFDiag::diagonalize(m_covs, m_mixing_matrix, tol, max_iter);
-	Eigen::Map<EMatrix> EQ(Q.matrix,n,n);
+	Map<MatrixXd> EQ(Q.matrix,n,n);
 
 	// Compute Mixing Matrix
 	m_mixing_matrix = SGMatrix<float64_t>(n,n);	
-	Eigen::Map<EMatrix> C(m_mixing_matrix.matrix,n,n);
+	Map<MatrixXd> C(m_mixing_matrix.matrix,n,n);
 	C = EQ.inverse();
 
 	// Normalize Estimated Mixing Matrix
@@ -108,25 +104,25 @@ CFeatures* CFFSep::apply(CFeatures* features)
 // Computing time delayed correlation matrix
 namespace 
 {
-	EMatrix cor(EMatrix x, int tau, bool mean_flag)
+	MatrixXd cor(MatrixXd x, int tau, bool mean_flag)
 	{	
 		int m = x.rows();
 		int n = x.cols();
 
 		// Center the data
 		if ( mean_flag )
 		{		
-			EVector mean = x.rowwise().sum();
+			VectorXd mean = x.rowwise().sum();
 			mean /= n;
 			x = x.colwise() - mean;
 		}
 
 		// Time-delayed Signal Matrix
-		EMatrix L = x.leftCols(n-tau);
-		EMatrix R = x.rightCols(n-tau);
+		MatrixXd L = x.leftCols(n-tau);
+		MatrixXd R = x.rightCols(n-tau);
 
 		// Compute Correlations
-		EMatrix K(m,m);
+		MatrixXd K(m,m);
 		K = (L * R.transpose()) / (n-tau);
 
 		// Symmetrize

src/shogun/converter/ica/FastICA.cpp

@@ -17,20 +17,16 @@
 #include <shogun/mathematics/Math.h>
 #include <shogun/mathematics/eigen3.h>
 
-using namespace Eigen;
-
-typedef Matrix< float64_t, Dynamic, 1, ColMajor > EVector;
-typedef Matrix< float64_t, Dynamic, Dynamic, ColMajor > EMatrix;
-
 using namespace shogun;
+using namespace Eigen;
 
 namespace {
 
-	EMatrix sym_decorrelation(EMatrix W)
+	MatrixXd sym_decorrelation(MatrixXd W)
 	{
-		EMatrix K = W * W.transpose();
+		MatrixXd K = W * W.transpose();
 
-		SelfAdjointEigenSolver<EMatrix> eig;
+		SelfAdjointEigenSolver<MatrixXd> eig;
 		eig.compute(K);
 
 		return ((eig.eigenvectors() * eig.eigenvalues().cwiseSqrt().asDiagonal().inverse()) * eig.eigenvectors().transpose()) * W;
@@ -86,21 +82,21 @@ CFeatures* CFastICA::apply(CFeatures* features)
 	int p = X.num_cols;
 	int m = n;
 
-	Eigen::Map<EMatrix> EX(X.matrix,n,p);
+	Map<MatrixXd> EX(X.matrix,n,p);
 
 	// Whiten
-	EMatrix K;
-	EMatrix WX;
+	MatrixXd K;
+	MatrixXd WX;
 	if (whiten)
 	{
-		EVector mean = (EX.rowwise().sum() / (float64_t)p);	
-		EMatrix SPX = EX.colwise() - mean;
+		VectorXd mean = (EX.rowwise().sum() / (float64_t)p);	
+		MatrixXd SPX = EX.colwise() - mean;
 
-		Eigen::JacobiSVD<EMatrix> svd;
+		Eigen::JacobiSVD<MatrixXd> svd;
 		svd.compute(SPX, Eigen::ComputeThinU);
 
-		EMatrix u = svd.matrixU();
-		EMatrix d = svd.singularValues();
+		MatrixXd u = svd.matrixU();
+		MatrixXd d = svd.singularValues();
 
 		// for matching numpy/scikit-learn
 		//u.rightCols(u.cols() - 1) *= -1;
@@ -132,20 +128,20 @@ CFeatures* CFastICA::apply(CFeatures* features)
 		}
 	}
 
-	Eigen::Map<EMatrix> W(m_mixing_matrix.matrix, m, m);
+	Map<MatrixXd> W(m_mixing_matrix.matrix, m, m);
 
 	W = sym_decorrelation(W);
 
 	int iter = 0;
 	float64_t lim = tol+1;
 	while (lim > tol && iter < max_iter)
 	{
-		EMatrix wtx = W * WX;
+		MatrixXd wtx = W * WX;
 
-		EMatrix gwtx = wtx.unaryExpr(std::ptr_fun(&gx));
-		EMatrix g_wtx = wtx.unaryExpr(std::ptr_fun(&g_x));
+		MatrixXd gwtx = wtx.unaryExpr(std::ptr_fun(&gx));
+		MatrixXd g_wtx = wtx.unaryExpr(std::ptr_fun(&g_x));
 
-		EMatrix W1 = (gwtx * WX.transpose()) / (float64_t)p - (g_wtx.rowwise().sum()/(float64_t)p).asDiagonal() * W;
+		MatrixXd W1 = (gwtx * WX.transpose()) / (float64_t)p - (g_wtx.rowwise().sum()/(float64_t)p).asDiagonal() * W;
 
 		W1 = sym_decorrelation(W1);
 

src/shogun/converter/ica/Jade.cpp

@@ -21,12 +21,8 @@
 #include <iostream>
 #endif
 
-using namespace Eigen;
-
-typedef Matrix< float64_t, Dynamic, 1, ColMajor > EVector;
-typedef Matrix< float64_t, Dynamic, Dynamic, ColMajor > EMatrix;
-
 using namespace shogun;
+using namespace Eigen;
 
 CJade::CJade() : CICAConverter()
 {	
@@ -59,21 +55,21 @@ CFeatures* CJade::apply(CFeatures* features)
 	int T = X.num_cols;
 	int m = n;
 
-	Eigen::Map<EMatrix> EX(X.matrix,n,T);
+	Map<MatrixXd> EX(X.matrix,n,T);
 
 	// Mean center X
-	EVector mean = (EX.rowwise().sum() / (float64_t)T);	
-	EMatrix SPX = EX.colwise() - mean;
+	VectorXd mean = (EX.rowwise().sum() / (float64_t)T);	
+	MatrixXd SPX = EX.colwise() - mean;
 
-	EMatrix cov = (SPX * SPX.transpose()) / (float64_t)T;
+	MatrixXd cov = (SPX * SPX.transpose()) / (float64_t)T;
 
 	#ifdef DEBUG_JADE
 	std::cout << "cov" << std::endl;
 	std::cout << cov << std::endl;
 	#endif
 
 	// Whitening & Projection onto signal subspace
-	SelfAdjointEigenSolver<EMatrix> eig;
+	SelfAdjointEigenSolver<MatrixXd> eig;
 	eig.compute(cov);
 
 	#ifdef DEBUG_JADE
@@ -85,8 +81,8 @@ CFeatures* CJade::apply(CFeatures* features)
 	#endif
 
 	// Scaling
-	EVector scales = eig.eigenvalues().cwiseSqrt();
-	EMatrix B = scales.cwiseInverse().asDiagonal() * eig.eigenvectors().transpose();
+	VectorXd scales = eig.eigenvalues().cwiseSqrt();
+	MatrixXd B = scales.cwiseInverse().asDiagonal() * eig.eigenvectors().transpose();
 
 	#ifdef DEBUG_JADE
 	std::cout << "whitener" << std::endl;	
@@ -100,12 +96,12 @@ CFeatures* CJade::apply(CFeatures* features)
 	int dimsymm = (m * ( m + 1)) / 2; // Dim. of the space of real symm matrices
 	int nbcm = dimsymm; //  number of cumulant matrices	
 	m_cumulant_matrix = SGMatrix<float64_t>(m,m*nbcm);	// Storage for cumulant matrices
-	Eigen::Map<EMatrix> CM(m_cumulant_matrix.matrix,m,m*nbcm); 
-	EMatrix R(m,m); R.setIdentity();
-	EMatrix Qij = EMatrix::Zero(m,m); // Temp for a cum. matrix
-	EVector Xim = EVector::Zero(m); // Temp
-	EVector Xjm = EVector::Zero(m); // Temp
-	EVector Xijm = EVector::Zero(m); // Temp
+	Map<MatrixXd> CM(m_cumulant_matrix.matrix,m,m*nbcm); 
+	MatrixXd R(m,m); R.setIdentity();
+	MatrixXd Qij = MatrixXd::Zero(m,m); // Temp for a cum. matrix
+	VectorXd Xim = VectorXd::Zero(m); // Temp
+	VectorXd Xjm = VectorXd::Zero(m); // Temp
+	VectorXd Xijm = VectorXd::Zero(m); // Temp
 	int Range = 0;
 
 	for (int im = 0; im < m; im++)
@@ -139,13 +135,13 @@ CFeatures* CJade::apply(CFeatures* features)
 
 	for (int i = 0; i < nbcm; i++)
 	{
-		Eigen::Map<EMatrix> EM(M.get_matrix(i),m,m);
+		Map<MatrixXd> EM(M.get_matrix(i),m,m);
 		EM = CM.block(0,i*m,m,m);
 	}
 
 	// Diagonalize
 	SGMatrix<float64_t> Q = CJADiagOrth::diagonalize(M, m_mixing_matrix, tol, max_iter);
-	Eigen::Map<EMatrix> EQ(Q.matrix,m,m);
+	Map<MatrixXd> EQ(Q.matrix,m,m);
 	EQ = -1 * EQ.inverse();
 
 	#ifdef DEBUG_JADE
@@ -155,11 +151,11 @@ CFeatures* CJade::apply(CFeatures* features)
 
 	// Separating matrix
 	SGMatrix<float64_t> sep_matrix = SGMatrix<float64_t>(m,m);
-	Eigen::Map<EMatrix> C(sep_matrix.matrix,m,m);
+	Map<MatrixXd> C(sep_matrix.matrix,m,m);
 	C = EQ.transpose() * B;
 
 	// Sort
-	EVector A = (B.inverse()*EQ).cwiseAbs2().colwise().sum();
+	VectorXd A = (B.inverse()*EQ).cwiseAbs2().colwise().sum();
 	bool swap = false;
 	do
 	{
@@ -180,7 +176,7 @@ CFeatures* CJade::apply(CFeatures* features)
 		C.row(j).swap( C.row(m-1-j) ); 
 
 	// Fix Signs
-	EVector signs = EVector::Zero(m);
+	VectorXd signs = VectorXd::Zero(m);
 	for (int i = 0; i < m; i++)
 	{
 		if( C(i,0) < 0 )
@@ -199,7 +195,7 @@ CFeatures* CJade::apply(CFeatures* features)
 	EX = C * EX;
 
 	m_mixing_matrix = SGMatrix<float64_t>(m,m);
-	Eigen::Map<EMatrix> Emixing_matrix(m_mixing_matrix.matrix,m,m);
+	Map<MatrixXd> Emixing_matrix(m_mixing_matrix.matrix,m,m);
 	Emixing_matrix = C.inverse();
 
 	return features;