SHOGUN  4.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
SOBI.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2013 Kevin Hughes
8  */
9 
11 
13 
14 #ifdef HAVE_EIGEN3
15 
19 
20 using namespace shogun;
21 using namespace Eigen;
22 
23 namespace { MatrixXd cor(MatrixXd x, int tau = 0, bool mean_flag = true); };
24 
26 {
27  init();
28 }
29 
31 {
32  m_tau = SGVector<float64_t>(4);
33  m_tau[0]=0; m_tau[1]=1; m_tau[2]=2; m_tau[3]=3;
34 
35  m_covs = SGNDArray<float64_t>();
36 
37  SG_ADD(&m_tau, "tau", "tau vector", MS_AVAILABLE);
38 }
39 
41 {
42 }
43 
45 {
46  m_tau = tau;
47 }
48 
50 {
51  return m_tau;
52 }
53 
55 {
56  return m_covs;
57 }
58 
60 {
61  ASSERT(features);
62  SG_REF(features);
63 
64  SGMatrix<float64_t> X = ((CDenseFeatures<float64_t>*)features)->get_feature_matrix();
65 
66  int n = X.num_rows;
67  int m = X.num_cols;
68  int N = m_tau.vlen;
69 
70  Map<MatrixXd> EX(X.matrix,n,m);
71 
72  // Whitening or Sphering
73  MatrixXd M0 = cor(EX,int(m_tau[0]));
74  EigenSolver<MatrixXd> eig;
75  eig.compute(M0);
76  MatrixXd SPH = (eig.pseudoEigenvectors() * eig.pseudoEigenvalueMatrix().cwiseSqrt() * eig.pseudoEigenvectors ().transpose()).inverse();
77  MatrixXd spx = SPH*EX;
78 
79  // Compute Correlation Matrices
80  index_t * M_dims = SG_MALLOC(index_t, 3);
81  M_dims[0] = n;
82  M_dims[1] = n;
83  M_dims[2] = N;
84  m_covs = SGNDArray< float64_t >(M_dims, 3);
85 
86  for(int t = 0; t < N; t++)
87  {
88  Map<MatrixXd> EM(m_covs.get_matrix(t),n,n);
89  EM = cor(spx,m_tau[t]);
90  }
91 
92  // Diagonalize
94  Map<MatrixXd> EQ(Q.matrix,n,n);
95 
96  // Compute Mixing Matrix
99  C = SPH.inverse() * EQ.transpose();
100 
101  // Normalize Estimated Mixing Matrix
102  for(int t = 0; t < C.cols(); t++)
103  {
104  C.col(t) /= C.col(t).maxCoeff();
105  }
106 
107  // Unmix
108  EX = C.inverse() * EX;
109 
110  return features;
111 }
112 
113 // Computing time delayed correlation matrix
114 namespace
115 {
116  MatrixXd cor(MatrixXd x, int tau, bool mean_flag)
117  {
118  int m = x.rows();
119  int n = x.cols();
120 
121  // Center the data
122  if ( mean_flag )
123  {
124  VectorXd mean = x.rowwise().sum();
125  mean /= n;
126  x = x.colwise() - mean;
127  }
128 
129  // Time-delayed Signal Matrix
130  MatrixXd L = x.leftCols(n-tau);
131  MatrixXd R = x.rightCols(n-tau);
132 
133  // Compute Correlations
134  MatrixXd K(m,m);
135  K = (L * R.transpose()) / (n-tau);
136 
137  // Symmetrize
138  K = (K + K.transpose()) / 2.0;
139 
140  return K;
141  }
142 };
143 #endif // HAVE_EIGEN3
int32_t index_t
Definition: common.h:62
void set_tau(SGVector< float64_t > tau)
Definition: SOBI.cpp:44
static SGMatrix< float64_t > diagonalize(SGNDArray< float64_t > C, SGMatrix< float64_t > V0=SGMatrix< float64_t >(NULL, 0, 0, false), double eps=CMath::MACHINE_EPSILON, int itermax=200)
Definition: JADiagOrth.cpp:18
Definition: SGMatrix.h:20
index_t num_cols
Definition: SGMatrix.h:378
#define SG_REF(x)
Definition: SGObject.h:51
index_t num_rows
Definition: SGMatrix.h:376
T * get_matrix(index_t matIdx) const
Definition: SGNDArray.h:72
index_t vlen
Definition: SGVector.h:494
virtual CFeatures * apply(CFeatures *features)
Definition: SOBI.cpp:59
#define ASSERT(x)
Definition: SGIO.h:201
void init()
Definition: SOBI.cpp:30
SGVector< float64_t > get_tau() const
Definition: SOBI.cpp:49
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
SGNDArray< float64_t > get_covs() const
Definition: SOBI.cpp:54
class ICAConverter Base class for ICA algorithms
Definition: ICAConverter.h:27
virtual ~CSOBI()
Definition: SOBI.cpp:40
The class Features is the base class of all feature objects.
Definition: Features.h:68
SGMatrix< float64_t > m_mixing_matrix
Definition: ICAConverter.h:83
static void inverse(SGMatrix< float64_t > matrix)
inverses square matrix in-place
Definition: SGMatrix.cpp:885
#define SG_ADD(...)
Definition: SGObject.h:81

SHOGUN Machine Learning Toolbox - Documentation