SHOGUN  4.2.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 
18 
19 using namespace shogun;
20 using namespace Eigen;
21 
22 namespace { MatrixXd cor(MatrixXd x, int tau = 0, bool mean_flag = true); };
23 
25 {
26  init();
27 }
28 
30 {
31  m_tau = SGVector<float64_t>(4);
32  m_tau[0]=0; m_tau[1]=1; m_tau[2]=2; m_tau[3]=3;
33 
34  m_covs = SGNDArray<float64_t>();
35 
36  SG_ADD(&m_tau, "tau", "tau vector", MS_AVAILABLE);
37 }
38 
40 {
41 }
42 
44 {
45  m_tau = tau;
46 }
47 
49 {
50  return m_tau;
51 }
52 
54 {
55  return m_covs;
56 }
57 
59 {
60  ASSERT(features);
61  SG_REF(features);
62 
63  SGMatrix<float64_t> X = ((CDenseFeatures<float64_t>*)features)->get_feature_matrix();
64 
65  int n = X.num_rows;
66  int m = X.num_cols;
67  int N = m_tau.vlen;
68 
69  Map<MatrixXd> EX(X.matrix,n,m);
70 
71  // Whitening or Sphering
72  MatrixXd M0 = cor(EX,int(m_tau[0]));
73  EigenSolver<MatrixXd> eig;
74  eig.compute(M0);
75  MatrixXd SPH = (eig.pseudoEigenvectors() * eig.pseudoEigenvalueMatrix().cwiseSqrt() * eig.pseudoEigenvectors ().transpose()).inverse();
76  MatrixXd spx = SPH*EX;
77 
78  // Compute Correlation Matrices
79  index_t * M_dims = SG_MALLOC(index_t, 3);
80  M_dims[0] = n;
81  M_dims[1] = n;
82  M_dims[2] = N;
83  m_covs = SGNDArray< float64_t >(M_dims, 3);
84 
85  for(int t = 0; t < N; t++)
86  {
87  Map<MatrixXd> EM(m_covs.get_matrix(t),n,n);
88  EM = cor(spx,m_tau[t]);
89  }
90 
91  // Diagonalize
93  Map<MatrixXd> EQ(Q.matrix,n,n);
94 
95  // Compute Mixing Matrix
98  C = SPH.inverse() * EQ.transpose();
99 
100  // Normalize Estimated Mixing Matrix
101  for(int t = 0; t < C.cols(); t++)
102  {
103  C.col(t) /= C.col(t).maxCoeff();
104  }
105 
106  // Unmix
107  EX = C.inverse() * EX;
108 
109  return features;
110 }
111 
112 // Computing time delayed correlation matrix
113 namespace
114 {
115  MatrixXd cor(MatrixXd x, int tau, bool mean_flag)
116  {
117  int m = x.rows();
118  int n = x.cols();
119 
120  // Center the data
121  if ( mean_flag )
122  {
123  VectorXd mean = x.rowwise().sum();
124  mean /= n;
125  x = x.colwise() - mean;
126  }
127 
128  // Time-delayed Signal Matrix
129  MatrixXd L = x.leftCols(n-tau);
130  MatrixXd R = x.rightCols(n-tau);
131 
132  // Compute Correlations
133  MatrixXd K(m,m);
134  K = (L * R.transpose()) / (n-tau);
135 
136  // Symmetrize
137  K = (K + K.transpose()) / 2.0;
138 
139  return K;
140  }
141 };
int32_t index_t
Definition: common.h:62
void set_tau(SGVector< float64_t > tau)
Definition: SOBI.cpp:43
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:17
Definition: SGMatrix.h:20
index_t num_cols
Definition: SGMatrix.h:376
#define SG_REF(x)
Definition: SGObject.h:54
index_t num_rows
Definition: SGMatrix.h:374
T * get_matrix(index_t matIdx) const
Definition: SGNDArray.h:72
index_t vlen
Definition: SGVector.h:492
virtual CFeatures * apply(CFeatures *features)
Definition: SOBI.cpp:58
#define ASSERT(x)
Definition: SGIO.h:201
void init()
Definition: SOBI.cpp:29
SGVector< float64_t > get_tau() const
Definition: SOBI.cpp:48
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:53
class ICAConverter Base class for ICA algorithms
Definition: ICAConverter.h:26
virtual ~CSOBI()
Definition: SOBI.cpp:39
The class Features is the base class of all feature objects.
Definition: Features.h:68
SGMatrix< float64_t > m_mixing_matrix
Definition: ICAConverter.h:82
static void inverse(SGMatrix< float64_t > matrix)
inverses square matrix in-place
Definition: SGMatrix.cpp:881
#define SG_ADD(...)
Definition: SGObject.h:84

SHOGUN Machine Learning Toolbox - Documentation