SHOGUN  5.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  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  float64_t z[N][N];
93  float64_t y[N][N];
94 
95  for (int i = 0; i < N; i++)
96  {
97  for (int j = 0; j < N; j++)
98  {
99  z[i][j] = 0;
100  y[i][j] = 0;
101  }
102  }
103 
104  for (int i = 0; i < N; i++)
105  {
106  for (int j = 0; j < N; j++)
107  {
108  for (int k = 0; k < K; k++)
109  {
110  auxij = N*N*k+N*i+j;
111  auxji = N*N*k+N*j+i;
112  auxii = N*N*k+N*i+i;
113  auxjj = N*N*k+N*j+j;
114  z[i][j] += C[auxii]*C[auxjj];
115  y[i][j] += 0.5*C[auxjj]*(C[auxij]+C[auxji]);
116  }
117  }
118  }
119 
120  for (int i = 0; i < N-1; i++)
121  {
122  for (int j = i+1; j < N; j++)
123  {
124  auxij = N*i+j;
125  auxji = N*j+i;
126  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]);
127  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]);
128  }
129  }
130 
131  return;
132 }
int32_t index_t
Definition: common.h:62
static float64_t ceil(float64_t d)
Definition: Math.h:416
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:20
SGMatrix< T > clone()
Definition: SGMatrix.cpp:256
index_t num_cols
Definition: SGMatrix.h:376
index_t num_rows
Definition: SGMatrix.h:374
T * get_matrix(index_t matIdx) const
Definition: SGNDArray.h:72
double float64_t
Definition: common.h:50
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:179
Matrix::Scalar max(Matrix m)
Definition: Redux.h:68
static SGMatrix< T > create_identity_matrix(index_t size, T scale)
static T abs(T a)
Definition: Math.h:179

SHOGUN Machine Learning Toolbox - Documentation