SHOGUN  6.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules
FFDiag.cpp
Go to the documentation of this file.
2 
3 
4 #include <shogun/base/init.h>
5 
8 
9 using namespace shogun;
10 using namespace Eigen;
11 
12 void getW(float64_t *C, int *ptN, int *ptK, float64_t *W);
13 
15  double eps, int itermax)
16 {
17  int n = C0.dims[0];
18  int K = C0.dims[2];
19 
20  index_t * C_dims = SG_MALLOC(index_t, 3);
21  C_dims[0] = C0.dims[0];
22  C_dims[1] = C0.dims[1];
23  C_dims[2] = C0.dims[2];
24  SGNDArray<float64_t> C(C_dims,3);
25  sg_memcpy(C.array, C0.array, C0.dims[0]*C0.dims[1]*C0.dims[2]*sizeof(float64_t));
26 
28  if (V0.num_rows == n && V0.num_cols == n)
29  V = V0.clone();
30  else
32 
33  MatrixXd Id(n,n); Id.setIdentity();
34  Map<MatrixXd> EV(V.matrix,n,n);
35 
36  float64_t inum = 0;
37  float64_t df = 1;
38  std::vector<float64_t> crit;
39  while (df > eps && inum < itermax)
40  {
41  MatrixXd W = MatrixXd::Zero(n,n);
42 
43  getW(C.get_matrix(0),
44  &n, &K,
45  W.data());
46 
47  W.transposeInPlace();
48  int e = CMath::ceil(log2(W.array().abs().rowwise().sum().maxCoeff()));
49  int s = std::max(0,e-1);
50  W /= pow(2,s);
51 
52  EV = (Id+W) * EV;
53  MatrixXd d = MatrixXd::Zero(EV.rows(),EV.cols());
54  d.diagonal() = VectorXd::Ones(EV.diagonalSize()).cwiseQuotient((EV * EV.transpose()).diagonal().cwiseSqrt());
55  EV = d * EV;
56 
57  for (int i = 0; i < K; i++)
58  {
59  Map<MatrixXd> Ci(C.get_matrix(i), n, n);
60  Map<MatrixXd> C0i(C0.get_matrix(i), n, n);
61  Ci = EV * C0i * EV.transpose();
62  }
63 
64  float64_t f = 0;
65  for (int i = 0; i < K; i++)
66  {
67  Map<MatrixXd> C0i(C0.get_matrix(i), n, n);
68  MatrixXd F = EV * C0i * EV.transpose();
69  f += (F.transpose() * F).diagonal().sum() - F.array().pow(2).matrix().diagonal().sum();
70  }
71 
72  crit.push_back(f);
73 
74  if (inum > 1)
75  df = CMath::abs(crit[inum-1]-crit[inum]);
76 
77  inum++;
78  }
79 
80  if (inum == itermax)
81  SG_SERROR("Convergence not reached\n")
82 
83  return V;
84 
85 }
86 
87 void getW(float64_t *C, int *ptN, int *ptK, float64_t *W)
88 {
89  int N=*ptN;
90  int K=*ptK;
91  int auxij,auxji,auxii,auxjj;
92  SGMatrix<float64_t> z(N, N);
93  SGMatrix<float64_t> y(N, N);
94 
95  z.zero();
96  y.zero();
97 
98  for (int i = 0; i < N; i++)
99  {
100  for (int j = 0; j < N; j++)
101  {
102  for (int k = 0; k < K; k++)
103  {
104  auxij = N*N*k+N*i+j;
105  auxji = N*N*k+N*j+i;
106  auxii = N*N*k+N*i+i;
107  auxjj = N*N*k+N*j+j;
108  z(i, j) += C[auxii]*C[auxjj];
109  y(i, j) += 0.5*C[auxjj]*(C[auxij]+C[auxji]);
110  }
111  }
112  }
113 
114  for (int i = 0; i < N-1; i++)
115  {
116  for (int j = i+1; j < N; j++)
117  {
118  auxij = N*i+j;
119  auxji = N*j+i;
120  W[auxij] = (z(j, i)*y(j, i) - z(i, i)*y(i, j))/(z(j, j)*z(i, i)-z(i, j)*z(i, j));
121  W[auxji] = (z(i, j)*y(i, j) - z(j, j)*y(j, i))/(z(j, j)*z(i, i)-z(i, j)*z(i, j));
122  }
123  }
124 
125  return;
126 }
SGMatrix< T > clone() const
Definition: SGMatrix.cpp:330
int32_t index_t
Definition: common.h:72
static float64_t ceil(float64_t d)
Definition: Math.h:411
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: FFDiag.cpp:14
Definition: SGMatrix.h:24
index_t num_cols
Definition: SGMatrix.h:465
index_t num_rows
Definition: SGMatrix.h:463
T * get_matrix(index_t matIdx) const
Definition: SGNDArray.h:72
double float64_t
Definition: common.h:60
void getW(float64_t *C, int *ptN, int *ptK, float64_t *W)
Definition: FFDiag.cpp:87
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
index_t * dims
Definition: SGNDArray.h:177
#define SG_SERROR(...)
Definition: SGIO.h:178
T max(const Container< T > &a)
static SGMatrix< T > create_identity_matrix(index_t size, T scale)
static T abs(T a)
Definition: Math.h:175

SHOGUN Machine Learning Toolbox - Documentation