SHOGUN  6.0.0
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
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
20 #include <iostream>
21 #endif
22
23 using namespace shogun;
24 using namespace Eigen;
25
27 {
28  init();
29 }
30
32 {
33  m_cumulant_matrix = SGMatrix<float64_t>();
35 }
36
38 {
39 }
40
42 {
43  return m_cumulant_matrix;
44 }
45
47 {
48  ASSERT(features);
49  SG_REF(features);
50
51  SGMatrix<float64_t> X = ((CDenseFeatures<float64_t>*)features)->get_feature_matrix();
52
53  int n = X.num_rows;
54  int T = X.num_cols;
55  int m = n;
56
57  Map<MatrixXd> EX(X.matrix,n,T);
58
59  // Mean center X
60  VectorXd mean = (EX.rowwise().sum() / (float64_t)T);
61  MatrixXd SPX = EX.colwise() - mean;
62
63  MatrixXd cov = (SPX * SPX.transpose()) / (float64_t)T;
64
66  std::cout << "cov" << std::endl;
67  std::cout << cov << std::endl;
68  #endif
69
70  // Whitening & Projection onto signal subspace
72  eig.compute(cov);
73
75  std::cout << "eigenvectors" << std::endl;
76  std::cout << eig.eigenvectors() << std::endl;
77
78  std::cout << "eigenvalues" << std::endl;
79  std::cout << eig.eigenvalues().asDiagonal() << std::endl;
80  #endif
81
82  // Scaling
83  VectorXd scales = eig.eigenvalues().cwiseSqrt();
84  MatrixXd B = scales.cwiseInverse().asDiagonal() * eig.eigenvectors().transpose();
85
87  std::cout << "whitener" << std::endl;
88  std::cout << B << std::endl;
89  #endif
90
91  // Sphering
92  SPX = B * SPX;
93
94  // Estimation of the cumulant matrices
95  int dimsymm = (m * ( m + 1)) / 2; // Dim. of the space of real symm matrices
96  int nbcm = dimsymm; // number of cumulant matrices
97  m_cumulant_matrix = SGMatrix<float64_t>(m,m*nbcm); // Storage for cumulant matrices
98  Map<MatrixXd> CM(m_cumulant_matrix.matrix,m,m*nbcm);
99  MatrixXd R(m,m); R.setIdentity();
100  MatrixXd Qij = MatrixXd::Zero(m,m); // Temp for a cum. matrix
101  VectorXd Xim = VectorXd::Zero(m); // Temp
102  VectorXd Xjm = VectorXd::Zero(m); // Temp
103  VectorXd Xijm = VectorXd::Zero(m); // Temp
104  int Range = 0;
105
106  for (int im = 0; im < m; im++)
107  {
108  Xim = SPX.row(im);
109  Xijm = Xim.cwiseProduct(Xim);
110  Qij = SPX.cwiseProduct(Xijm.replicate(1,m).transpose()) * SPX.transpose() / (float)T - R - 2*R.col(im)*R.col(im).transpose();
111  CM.block(0,Range,m,m) = Qij;
112  Range = Range + m;
113  for (int jm = 0; jm < im; jm++)
114  {
115  Xjm = SPX.row(jm);
116  Xijm = Xim.cwiseProduct(Xjm);
117  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();
118  CM.block(0,Range,m,m) = sqrt(2)*Qij;
119  Range = Range + m;
120  }
121  }
122
124  std::cout << "cumulatant matrices" << std::endl;
125  std::cout << CM << std::endl;
126  #endif
127
128  // Stack cumulant matrix into ND Array
129  index_t * M_dims = SG_MALLOC(index_t, 3);
130  M_dims[0] = m;
131  M_dims[1] = m;
132  M_dims[2] = nbcm;
133  SGNDArray< float64_t > M(M_dims, 3);
134
135  for (int i = 0; i < nbcm; i++)
136  {
137  Map<MatrixXd> EM(M.get_matrix(i),m,m);
138  EM = CM.block(0,i*m,m,m);
139  }
140
141  // Diagonalize
143  Map<MatrixXd> EQ(Q.matrix,m,m);
144  EQ = -1 * EQ.inverse();
145
147  std::cout << "diagonalizer" << std::endl;
148  std::cout << EQ << std::endl;
149  #endif
150
151  // Separating matrix
152  SGMatrix<float64_t> sep_matrix = SGMatrix<float64_t>(m,m);
153  Map<MatrixXd> C(sep_matrix.matrix,m,m);
154  C = EQ.transpose() * B;
155
156  // Sort
157  VectorXd A = (B.inverse()*EQ).cwiseAbs2().colwise().sum();
158  bool swap = false;
159  do
160  {
161  swap = false;
162  for (int j = 1; j < n; j++)
163  {
164  if ( A(j) < A(j-1) )
165  {
166  std::swap(A(j),A(j-1));
167  C.col(j).swap(C.col(j-1));
168  swap = true;
169  }
170  }
171
172  } while(swap);
173
174  for (int j = 0; j < m/2; j++)
175  C.row(j).swap( C.row(m-1-j) );
176
177  // Fix Signs
178  VectorXd signs = VectorXd::Zero(m);
179  for (int i = 0; i < m; i++)
180  {
181  if( C(i,0) < 0 )
182  signs(i) = -1;
183  else
184  signs(i) = 1;
185  }
186  C = signs.asDiagonal() * C;
187
189  std::cout << "un mixing matrix" << std::endl;
190  std::cout << C << std::endl;
191  #endif
192
193  // Unmix
194  EX = C * EX;
195
197  Map<MatrixXd> Emixing_matrix(m_mixing_matrix.matrix,m,m);
198  Emixing_matrix = C.inverse();
199
200  return features;
201 }
202
int32_t index_t
Definition: common.h:72
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: SGMatrix.h:24
index_t num_cols
Definition: SGMatrix.h:465
std::enable_if<!std::is_same< T, complex128_t >::value, float64_t >::type mean(const Container< T > &a)
void init()
#define SG_REF(x)
Definition: SGObject.h:52
index_t num_rows
Definition: SGMatrix.h:463
T * get_matrix(index_t matIdx) const
Definition: SGNDArray.h:72
#define ASSERT(x)
Definition: SGIO.h:200
double float64_t
Definition: common.h:60
virtual CFeatures * apply(CFeatures *features)
SGMatrix< float64_t > get_cumulant_matrix() const
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
class ICAConverter Base class for ICA algorithms
Definition: ICAConverter.h:26
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:972