LCOV - code coverage report
Current view: top level - routines - diffusion_maps.hpp (source / functions) Hit Total Coverage
Test: clean.info Lines: 21 22 95.5 %
Date: 2013-05-24 Functions: 1 9 11.1 %
Branches: 33 217 15.2 %

           Branch data     Line data    Source code
       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                 :            : 
      19                 :            : //! Computes diffusion process matrix. Uses the following algorithm:
      20                 :            : //!
      21                 :            : //! <ol>
      22                 :            : //! <li> Compute matrix \f$ K \f$ such as \f$ K_{i,j} = \exp\left(-\frac{d(x_i,x_j)^2}{w}\right) \f$.
      23                 :            : //! <li> Compute sum vector \f$ p = \sum_i K_{i,j}\f$.
      24                 :            : //! <li> Modify \f$ K \f$ with \f$ K_{i,j} = K_{i,j} / (p_i  p_j)^t \f$.
      25                 :            : //! <li> Compute sum vector \f$ p = \sum_i K_{i,j}\f$ again.
      26                 :            : //! <li> Normalize \f$ K \f$ with \f$ K_{i,j} = K_{i,j} / (p_i p_j) \f$.
      27                 :            : //! </ol>
      28                 :            : //!
      29                 :            : //! @param begin begin data iterator
      30                 :            : //! @param end end data iterator
      31                 :            : //! @param callback distance callback
      32                 :            : //! @param timesteps number of timesteps \f$ t \f$ of diffusion process
      33                 :            : //! @param width width \f$ w \f$ of the gaussian kernel
      34                 :            : //!
      35                 :            : template <class RandomAccessIterator, class DistanceCallback>
      36                 :          1 : DenseSymmetricMatrix compute_diffusion_matrix(RandomAccessIterator begin, RandomAccessIterator end, DistanceCallback callback, 
      37                 :            :                                               const IndexType timesteps, const ScalarType width)
      38                 :            : {
      39 [ +  - ][ +  - ]:          1 :         timed_context context("Diffusion map matrix computation");
         [ +  - ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      40                 :            : 
      41                 :          1 :         const IndexType n_vectors = end-begin;
      42   [ +  -  #  #  :          1 :         DenseSymmetricMatrix diffusion_matrix(n_vectors,n_vectors);
             #  #  #  # ]
      43 [ +  - ][ +  - ]:          1 :         DenseVector p = DenseVector::Zero(n_vectors);
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      44                 :            : 
      45                 :          1 :         RESTRICT_ALLOC;
      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 [ #  # ][ #  # ]:          0 : #pragma omp for nowait
         [ #  # ][ #  # ]
                    [ + ]
      52                 :            :                 for (i_index_iter=0; i_index_iter<n_vectors; ++i_index_iter)
      53                 :            :                 {
      54 [ +  + ][ +  + ]:       1226 :                         for (j_index_iter=i_index_iter; j_index_iter<n_vectors; ++j_index_iter)
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      55                 :            :                         {
      56                 :       1208 :                                 ScalarType k = callback.distance(begin[i_index_iter],begin[j_index_iter]);
      57                 :       1116 :                                 ScalarType gk = exp(-(k*k)/width);
      58                 :       1116 :                                 diffusion_matrix(i_index_iter,j_index_iter) = gk;
      59                 :       1182 :                                 diffusion_matrix(j_index_iter,i_index_iter) = gk;
      60                 :            :                         }
      61                 :            :                 }
      62                 :            :         }
      63                 :            :         // compute column sum vector
      64 [ +  - ][ +  - ]:          1 :         p = diffusion_matrix.colwise().sum();
           [ +  -  #  # ]
                 [ #  # ]
           [ #  #  #  # ]
                 [ #  # ]
           [ #  #  #  # ]
         [ #  # ][ #  # ]
      65                 :            : 
      66                 :            :         // compute full matrix as we need to compute sum later
      67 [ +  + ][ #  # ]:         51 :         for (IndexType i=0; i<n_vectors; i++)
         [ #  # ][ #  # ]
      68 [ +  + ][ #  # ]:       2550 :                 for (IndexType j=0; j<n_vectors; j++)
         [ #  # ][ #  # ]
      69 [ +  - ][ +  - ]:       2500 :                         diffusion_matrix(i,j) /= pow(p(i)*p(j),timesteps);
         [ +  - ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      70                 :            : 
      71                 :            :         // compute sqrt of column sum vector
      72 [ +  - ][ +  - ]:          1 :         p = diffusion_matrix.colwise().sum().cwiseSqrt();
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      73                 :            :         
      74 [ +  + ][ #  # ]:         51 :         for (IndexType i=0; i<n_vectors; i++)
         [ #  # ][ #  # ]
      75 [ +  + ][ #  # ]:       1325 :                 for (IndexType j=i; j<n_vectors; j++)
         [ #  # ][ #  # ]
      76 [ +  - ][ +  - ]:       1275 :                         diffusion_matrix(i,j) /= p(i)*p(j);
         [ +  - ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      77                 :            : 
      78                 :          1 :         UNRESTRICT_ALLOC;
      79                 :            : 
      80 [ +  - ][ #  # ]:          1 :         return diffusion_matrix;
         [ #  # ][ #  # ]
      81                 :            : }
      82                 :            : 
      83                 :            : } // End of namespace tapkee_internal
      84                 :            : } // End of namespace tapkee
      85                 :            : 
      86                 :            : #endif

Generated by: LCOV version 1.9