Tapkee
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
diffusion_maps.hpp
Go to the documentation of this file.
1 /* This software is distributed under BSD 3-clause license (see LICENSE file).
2  *
3  * Copyright (c) 2012-2013 Sergey Lisitsyn
4  */
5 
6 #ifndef TAPKEE_DiffusionMapS_H_
7 #define TAPKEE_DiffusionMapS_H_
8 
9 /* Tapkee includes */
10 #include <tapkee/defines.hpp>
11 #include <tapkee/utils/time.hpp>
12 /* End of Tapke includes */
13 
14 namespace tapkee
15 {
16 namespace tapkee_internal
17 {
18 
35 template <class RandomAccessIterator, class DistanceCallback>
36 DenseSymmetricMatrix compute_diffusion_matrix(RandomAccessIterator begin, RandomAccessIterator end, DistanceCallback callback,
37  const IndexType timesteps, const ScalarType width)
38 {
39  timed_context context("Diffusion map matrix computation");
40 
41  const IndexType n_vectors = end-begin;
42  DenseSymmetricMatrix diffusion_matrix(n_vectors,n_vectors);
43  DenseVector p = DenseVector::Zero(n_vectors);
44 
46 
47  // compute gaussian kernel matrix
48 #pragma omp parallel shared(diffusion_matrix,begin,callback) default(none)
49  {
50  IndexType i_index_iter, j_index_iter;
51 #pragma omp for nowait
52  for (i_index_iter=0; i_index_iter<n_vectors; ++i_index_iter)
53  {
54  for (j_index_iter=i_index_iter; j_index_iter<n_vectors; ++j_index_iter)
55  {
56  ScalarType k = callback.distance(begin[i_index_iter],begin[j_index_iter]);
57  ScalarType gk = exp(-(k*k)/width);
58  diffusion_matrix(i_index_iter,j_index_iter) = gk;
59  diffusion_matrix(j_index_iter,i_index_iter) = gk;
60  }
61  }
62  }
63  // compute column sum vector
64  p = diffusion_matrix.colwise().sum();
65 
66  // compute full matrix as we need to compute sum later
67  for (IndexType i=0; i<n_vectors; i++)
68  for (IndexType j=0; j<n_vectors; j++)
69  diffusion_matrix(i,j) /= pow(p(i)*p(j),timesteps);
70 
71  // compute sqrt of column sum vector
72  p = diffusion_matrix.colwise().sum().cwiseSqrt();
73 
74  for (IndexType i=0; i<n_vectors; i++)
75  for (IndexType j=i; j<n_vectors; j++)
76  diffusion_matrix(i,j) /= p(i)*p(j);
77 
79 
80  return diffusion_matrix;
81 }
82 
83 } // End of namespace tapkee_internal
84 } // End of namespace tapkee
85 
86 #endif
double ScalarType
default scalar value (can be overrided with TAPKEE_CUSTOM_INTERNAL_NUMTYPE define) ...
Definition: types.hpp:15
DenseSymmetricMatrix compute_diffusion_matrix(RandomAccessIterator begin, RandomAccessIterator end, DistanceCallback callback, const IndexType timesteps, const ScalarType width)
Computes diffusion process matrix. Uses the following algorithm:
Eigen::Matrix< tapkee::ScalarType, Eigen::Dynamic, 1 > DenseVector
dense vector type (non-overridable)
Definition: types.hpp:21
int IndexType
indexing type (non-overridable) set to int for compatibility with OpenMP 2.0
Definition: types.hpp:19
#define UNRESTRICT_ALLOC
Definition: eigen3.hpp:27
#define RESTRICT_ALLOC
Definition: eigen3.hpp:26
tapkee::DenseMatrix DenseSymmetricMatrix
dense symmetric matrix (non-overridable, currently just dense matrix, can be improved later) ...
Definition: types.hpp:25