Merge pull request #1447 from pickle27/develop

added UWedgeSep Alg and python modular examples for all current ICA methods

examples/undocumented/python_modular/graphical/converter_ffsep_bss.py

@@ -0,0 +1,57 @@
+"""
+Blind Source Separation using the FFSep Algorithm with Shogun
+
+Based on the example from scikit-learn
+http://scikit-learn.org/
+
+Kevin Hughes 2013
+"""
+
+
+import numpy as np
+import pylab as pl
+
+from shogun.Features  import RealFeatures
+from shogun.Converter import FFSep
+
+# Generate sample data
+np.random.seed(0)
+n_samples = 2000
+time = np.linspace(0, 10, n_samples)
+
+# Source Signals
+s1 = np.sin(2 * time)  # sin wave
+s2 = np.sign(np.sin(3 * time))  # square wave
+S = np.c_[s1, s2]
+S += 0.2 * np.random.normal(size=S.shape)  # add noise
+
+# Standardize data
+S /= S.std(axis=0)  
+S = S.T
+
+# Mixing Matrix
+A = np.array([[1, 0.5], [0.5, 1]])
+
+# Mix Signals
+X = np.dot(A,S)
+mixed_signals = RealFeatures(X)
+
+# Separating
+ffsep = FFSep()
+signals = ffsep.apply(mixed_signals)
+S_ = signals.get_feature_matrix()
+A_ = ffsep.get_mixing_matrix();
+
+# Plot results
+pl.figure()
+pl.subplot(3, 1, 1)
+pl.plot(S.T)
+pl.title('True Sources')
+pl.subplot(3, 1, 2)
+pl.plot(X.T)
+pl.title('Mixed Sources')
+pl.subplot(3, 1, 3)
+pl.plot(S_.T)
+pl.title('Estimated Sources')
+pl.subplots_adjust(0.09, 0.04, 0.94, 0.94, 0.26, 0.36)
+pl.show()

examples/undocumented/python_modular/graphical/converter_jedi_bss.py

@@ -0,0 +1,59 @@
+"""
+Blind Source Separation using the JediSep Algorithm with Shogun
+
+Based on the example from scikit-learn
+http://scikit-learn.org/
+
+Kevin Hughes 2013
+"""
+
+
+import numpy as np
+import pylab as pl
+
+from shogun.Features  import RealFeatures
+from shogun.Converter import JediSep
+
+# Generate sample data
+np.random.seed(0)
+n_samples = 2000
+time = np.linspace(0, 10, n_samples)
+
+# Source Signals
+s1 = np.sin(2 * time)  # sin wave
+s2 = np.sign(np.sin(3 * time))  # square wave
+S = np.c_[s1, s2]
+S += 0.2 * np.random.normal(size=S.shape)  # add noise
+
+# Standardize data
+S /= S.std(axis=0)  
+S = S.T
+
+# Mixing Matrix
+A = np.array([[1, 0.5], [0.5, 1]])
+
+# Mix Signals
+X = np.dot(A,S)
+mixed_signals = RealFeatures(X)
+
+# Separating
+jedi = JediSep()
+signals = jedi.apply(mixed_signals)
+S_ = signals.get_feature_matrix()
+A_ = jedi.get_mixing_matrix();
+
+print A_
+
+# Plot results
+pl.figure()
+pl.subplot(3, 1, 1)
+pl.plot(S.T)
+pl.title('True Sources')
+pl.subplot(3, 1, 2)
+pl.plot(X.T)
+pl.title('Mixed Sources')
+pl.subplot(3, 1, 3)
+pl.plot(S_.T)
+pl.title('Estimated Sources')
+pl.subplots_adjust(0.09, 0.04, 0.94, 0.94, 0.26, 0.36)
+pl.show()

examples/undocumented/python_modular/graphical/converter_sobi_bss.py

@@ -0,0 +1,57 @@
+"""
+Blind Source Separation using the SOBI Algorithm with Shogun
+
+Based on the example from scikit-learn
+http://scikit-learn.org/
+
+Kevin Hughes 2013
+"""
+
+
+import numpy as np
+import pylab as pl
+
+from shogun.Features  import RealFeatures
+from shogun.Converter import SOBI
+
+# Generate sample data
+np.random.seed(0)
+n_samples = 2000
+time = np.linspace(0, 10, n_samples)
+
+# Source Signals
+s1 = np.sin(2 * time)  # sin wave
+s2 = np.sign(np.sin(3 * time))  # square wave
+S = np.c_[s1, s2]
+S += 0.2 * np.random.normal(size=S.shape)  # add noise
+
+# Standardize data
+S /= S.std(axis=0)  
+S = S.T
+
+# Mixing Matrix
+A = np.array([[1, 0.5], [0.5, 1]])
+
+# Mix Signals
+X = np.dot(A,S)
+mixed_signals = RealFeatures(X)
+
+# Separating
+sobi = SOBI()
+signals = sobi.apply(mixed_signals)
+S_ = signals.get_feature_matrix()
+A_ = sobi.get_mixing_matrix();
+
+# Plot results
+pl.figure()
+pl.subplot(3, 1, 1)
+pl.plot(S.T)
+pl.title('True Sources')
+pl.subplot(3, 1, 2)
+pl.plot(X.T)
+pl.title('Mixed Sources')
+pl.subplot(3, 1, 3)
+pl.plot(S_.T)
+pl.title('Estimated Sources')
+pl.subplots_adjust(0.09, 0.04, 0.94, 0.94, 0.26, 0.36)
+pl.show()

examples/undocumented/python_modular/graphical/converter_uwedge_bss.py

@@ -0,0 +1,56 @@
+"""
+Blind Source Separation using the UWedgeSep Algorithm with Shogun
+
+Kevin Hughes 2013
+"""
+
+
+import numpy as np
+import pylab as pl
+
+from shogun.Features  import RealFeatures
+from shogun.Converter import UWedgeSep
+
+# Generate sample data
+np.random.seed(0)
+n_samples = 2000
+time = np.linspace(0, 1, n_samples)
+
+# Source Signals
+s1 = np.sin(2*3.14*55*time)  # sin wave
+s2 = np.cos(2*3.14*100*time) # cos wave
+S = np.c_[s1, s2]
+S += 0.1 * np.random.normal(size=S.shape)  # add noise
+
+# Standardize data
+S /= S.std(axis=0)  
+S = S.T
+
+# Mixing Matrix
+A = np.array([[1, 0.85], [0.55, 1]])
+#print A
+
+# Mix Signals
+X = np.dot(A,S)
+mixed_signals = RealFeatures(X)
+
+# Separating
+uwedge = UWedgeSep()
+signals = uwedge.apply(mixed_signals)
+S_ = signals.get_feature_matrix()
+A_ = uwedge.get_mixing_matrix();
+#print A_
+
+# Plot results
+pl.figure()
+pl.subplot(3, 1, 1)
+pl.plot(S.T)
+pl.title('True Sources')
+pl.subplot(3, 1, 2)
+pl.plot(X.T)
+pl.title('Mixed Sources')
+pl.subplot(3, 1, 3)
+pl.plot(S_.T)
+pl.title('Estimated Sources')
+pl.subplots_adjust(0.09, 0.04, 0.94, 0.94, 0.26, 0.36)
+pl.show()

src/interfaces/modular/Converter.i

@@ -29,6 +29,7 @@
 %rename(SOBI) CSOBI;
 %rename(FFSep) CFFSep;
 %rename(JediSep) CJediSep;
+%rename(UWedgeSep) CUWedgeSep;
 
 %newobject shogun::CEmbeddingConverter::apply;
 %newobject shogun::*::embed_kernel;
@@ -56,3 +57,4 @@
 %include <shogun/converter/ica/SOBI.h>
 %include <shogun/converter/ica/FFSep.h>
 %include <shogun/converter/ica/JediSep.h>
+%include <shogun/converter/ica/UWedgeSep.h>

src/interfaces/modular/Converter_includes.i

@@ -21,4 +21,5 @@
 #include <shogun/converter/ica/SOBI.h>
 #include <shogun/converter/ica/FFSep.h>
 #include <shogun/converter/ica/JediSep.h>
+#include <shogun/converter/ica/UWedgeSep.h>
 %}

src/shogun/converter/ica/UWedgeSep.cpp

@@ -0,0 +1,140 @@
+/*
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Written (W) 2013 Kevin Hughes
+ */
+
+#include <shogun/converter/ica/UWedgeSep.h>
+
+#include <shogun/features/DenseFeatures.h>
+
+#ifdef HAVE_EIGEN3
+
+#include <shogun/mathematics/Math.h>
+#include <shogun/mathematics/eigen3.h>
+#include <shogun/mathematics/ajd/UWedge.h>
+
+using namespace Eigen;
+
+typedef Matrix< float64_t, Dynamic, 1, ColMajor > EVector;
+typedef Matrix< float64_t, Dynamic, Dynamic, ColMajor > EMatrix;
+
+using namespace shogun;
+
+namespace { EMatrix cor(EMatrix x, int tau = 0, bool mean_flag = true); };
+
+CUWedgeSep::CUWedgeSep() : CConverter()
+{
+	m_tau = SGVector<float64_t>(4); 
+	m_tau[0]=0; m_tau[1]=1; m_tau[2]=2; m_tau[3]=3;
+
+	init();
+}
+
+void CUWedgeSep::init()
+{
+	m_mixing_matrix = SGMatrix<float64_t>();
+
+	SG_ADD(&m_tau, "tau", "tau vector", MS_AVAILABLE);
+	SG_ADD(&m_mixing_matrix, "mixing_matrix", "m_mixing_matrix", MS_NOT_AVAILABLE);
+}
+
+CUWedgeSep::~CUWedgeSep()
+{
+}
+
+void CUWedgeSep::set_tau(SGVector<float64_t> tau)
+{
+	m_tau = tau;
+}
+
+SGVector<float64_t> CUWedgeSep::get_tau() const
+{
+	return m_tau;
+}
+
+SGMatrix<float64_t> CUWedgeSep::get_mixing_matrix() const
+{
+	return m_mixing_matrix;
+}
+
+CFeatures* CUWedgeSep::apply(CFeatures* features)
+{
+	REQUIRE(features, "features is null");	
+	SG_REF(features);
+
+	SGMatrix<float64_t> X = ((CDenseFeatures<float64_t>*)features)->get_feature_matrix();
+
+	int n = X.num_rows;
+	int m = X.num_cols;
+	int N = m_tau.vlen;
+
+	Eigen::Map<EMatrix> EX(X.matrix,n,m);	
+
+	// Compute Correlation Matrices
+	index_t * M_dims = SG_MALLOC(index_t, 3);
+	M_dims[0] = n;
+	M_dims[1] = n;
+	M_dims[2] = N;
+	SGNDArray< float64_t > M(M_dims, 3);
+
+	for (int t = 0; t < N; t++)
+	{
+		Eigen::Map<EMatrix> EM(M.get_matrix(t),n,n);
+		EM = cor(EX,m_tau[t]);
+	}
+
+	// Diagonalize
+	SGMatrix<float64_t> Q = CUWedge::diagonalize(M);
+	Eigen::Map<EMatrix> 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);
+	C = EQ.inverse();
+
+	// Normalize Estimated Mixing Matrix
+	for (int t = 0; t < C.cols(); t++)	
+		C.col(t) /= C.col(t).maxCoeff();
+
+	// Unmix
+	EX = C.inverse() * EX;
+
+	return features;
+}
+
+// Computing time delayed correlation matrix
+namespace 
+{
+	EMatrix cor(EMatrix 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();
+			mean /= n;
+			x = x.colwise() - mean;
+		}
+
+		// Time-delayed Signal Matrix
+		EMatrix L = x.leftCols(n-tau);
+		EMatrix R = x.rightCols(n-tau);
+
+		// Compute Correlations
+		EMatrix K(m,m);
+		K = (L * R.transpose()) / (n-tau);
+
+		// Symmetrize
+		K = (K + K.transpose()) / 2.0;
+
+		return K;
+	}
+};
+
+#endif // HAVE_EIGEN3