Tapkee
diffusion_maps.hpp
Go to the documentation of this file.
00001 /* This software is distributed under BSD 3-clause license (see LICENSE file).
00002  *
00003  * Copyright (c) 2012-2013 Sergey Lisitsyn
00004  */
00005 
00006 #ifndef TAPKEE_DiffusionMapS_H_
00007 #define TAPKEE_DiffusionMapS_H_
00008 
00009 /* Tapkee includes */
00010 #include <tapkee/defines.hpp>
00011 #include <tapkee/utils/time.hpp>
00012 /* End of Tapke includes */
00013 
00014 namespace tapkee
00015 {
00016 namespace tapkee_internal
00017 {
00018 
00035 template <class RandomAccessIterator, class DistanceCallback>
00036 DenseSymmetricMatrix compute_diffusion_matrix(RandomAccessIterator begin, RandomAccessIterator end, DistanceCallback callback, 
00037                                               const IndexType timesteps, const ScalarType width)
00038 {
00039     timed_context context("Diffusion map matrix computation");
00040 
00041     const IndexType n_vectors = end-begin;
00042     DenseSymmetricMatrix diffusion_matrix(n_vectors,n_vectors);
00043     DenseVector p = DenseVector::Zero(n_vectors);
00044 
00045     RESTRICT_ALLOC;
00046 
00047     // compute gaussian kernel matrix
00048 #pragma omp parallel shared(diffusion_matrix,begin,callback) default(none)
00049     {
00050         IndexType i_index_iter, j_index_iter;
00051 #pragma omp for nowait
00052         for (i_index_iter=0; i_index_iter<n_vectors; ++i_index_iter)
00053         {
00054             for (j_index_iter=i_index_iter; j_index_iter<n_vectors; ++j_index_iter)
00055             {
00056                 ScalarType k = callback.distance(begin[i_index_iter],begin[j_index_iter]);
00057                 ScalarType gk = exp(-(k*k)/width);
00058                 diffusion_matrix(i_index_iter,j_index_iter) = gk;
00059                 diffusion_matrix(j_index_iter,i_index_iter) = gk;
00060             }
00061         }
00062     }
00063     // compute column sum vector
00064     p = diffusion_matrix.colwise().sum();
00065 
00066     // compute full matrix as we need to compute sum later
00067     for (IndexType i=0; i<n_vectors; i++)
00068         for (IndexType j=0; j<n_vectors; j++)
00069             diffusion_matrix(i,j) /= pow(p(i)*p(j),timesteps);
00070 
00071     // compute sqrt of column sum vector
00072     p = diffusion_matrix.colwise().sum().cwiseSqrt();
00073     
00074     for (IndexType i=0; i<n_vectors; i++)
00075         for (IndexType j=i; j<n_vectors; j++)
00076             diffusion_matrix(i,j) /= p(i)*p(j);
00077 
00078     UNRESTRICT_ALLOC;
00079 
00080     return diffusion_matrix;
00081 }
00082 
00083 } // End of namespace tapkee_internal
00084 } // End of namespace tapkee
00085 
00086 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines