LCOV - code coverage report
Current view: top level - routines - laplacian_eigenmaps.hpp (source / functions) Hit Total Coverage
Test: clean.info Lines: 35 35 100.0 %
Date: 2013-05-24 Functions: 3 9 33.3 %
Branches: 86 364 23.6 %

           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_LaplacianEigenmaps_H_
       7                 :            : #define TAPKEE_LaplacianEigenmaps_H_
       8                 :            : 
       9                 :            : /* Tapkee includes */
      10                 :            : #include <tapkee/defines.hpp>
      11                 :            : #include <tapkee/utils/time.hpp>
      12                 :            : /* End of Tapkee includes */
      13                 :            : 
      14                 :            : namespace tapkee
      15                 :            : {
      16                 :            : namespace tapkee_internal
      17                 :            : {
      18                 :            : 
      19                 :            : //! Computes laplacian of neighborhood graph.
      20                 :            : //!
      21                 :            : //! Follows the algorithm described below:
      22                 :            : //! <ul>
      23                 :            : //! <li> For each vector compute gaussian exp of distances to its neighbor vectors and 
      24                 :            : //!      put it to sparse matrix \f$ L_{i,N_i(j)} = \exp\left( - \frac{d(x_i,x_{N_i(j)})^2}{w} \right) \f$.
      25                 :            : //! <li> Symmetrize matrix \f$ L \f$ with \f$ L_{i,j} = \max (L_{i,j}, L_{j,i}) \f$ to
      26                 :            : //!      make neighborhood relationship symmetric.
      27                 :            : //! <li> Compute sum vector \f$ D = \sum_{i} L_{i,j} \f$.
      28                 :            : //! <li> Modify \f$ L = D - L \f$.
      29                 :            : //! <li> Output matrix sparse matrix \f$ L \f$ and diagonal matrix built of vector \f$ D \f$.
      30                 :            : //! </ul>
      31                 :            : //!
      32                 :            : //! @param begin begin data iterator
      33                 :            : //! @param end end data iterator
      34                 :            : //! @param neighbors neighbors of each vector
      35                 :            : //! @param callback distance callback
      36                 :            : //! @param width width \f$ w \f$ of the gaussian kernel
      37                 :            : //!
      38                 :            : template<class RandomAccessIterator, class DistanceCallback>
      39                 :          3 : Laplacian compute_laplacian(RandomAccessIterator begin, 
      40                 :            :                         RandomAccessIterator end,const Neighbors& neighbors, 
      41                 :            :                         DistanceCallback callback, ScalarType width)
      42                 :            : {
      43                 :          3 :         SparseTriplets sparse_triplets;
      44                 :            : 
      45 [ +  - ][ +  - ]:          3 :         timed_context context("Laplacian computation");
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      46                 :          3 :         const IndexType k = neighbors[0].size();
      47   [ +  -  +  -  :          3 :         sparse_triplets.reserve((k+1)*(end-begin));
             #  #  #  # ]
      48                 :            : 
      49 [ +  - ][ +  - ]:          3 :         DenseVector D = DenseVector::Zero(end-begin);
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      50 [ +  + ][ +  + ]:        123 :         for (RandomAccessIterator iter=begin; iter!=end; ++iter)
         [ #  # ][ #  # ]
      51                 :            :         {
      52                 :        120 :                 const LocalNeighbors& current_neighbors = neighbors[iter-begin];
      53                 :            : 
      54 [ +  + ][ +  + ]:       1220 :                 for (IndexType i=0; i<k; ++i)
         [ #  # ][ #  # ]
      55                 :            :                 {
      56 [ #  # ][ #  # ]:       1100 :                         ScalarType distance = callback.distance(*iter,begin[current_neighbors[i]]);
                 [ +  - ]
      57                 :       1100 :                         ScalarType heat = exp(-distance*distance/width);
      58   [ #  #  +  - ]:       1100 :                         D(iter-begin) += heat;
         [ #  # ][ #  # ]
                 [ +  - ]
      59 [ +  - ][ +  - ]:       1100 :                         D(current_neighbors[i]) += heat;
         [ #  # ][ #  # ]
      60 [ +  - ][ +  - ]:       1100 :                         sparse_triplets.push_back(SparseTriplet(current_neighbors[i],(iter-begin),-heat));
         [ #  # ][ #  # ]
      61 [ +  - ][ +  - ]:       1100 :                         sparse_triplets.push_back(SparseTriplet((iter-begin),current_neighbors[i],-heat));
         [ #  # ][ #  # ]
      62                 :            :                 }
      63                 :            :         }
      64 [ +  + ][ +  + ]:        123 :         for (IndexType i=0; i<static_cast<IndexType>(end-begin); ++i)
         [ #  # ][ #  # ]
      65 [ +  - ][ +  - ]:        120 :                 sparse_triplets.push_back(SparseTriplet(i,i,D(i)));
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      66                 :            : 
      67                 :            : #ifdef EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
      68                 :            :         Eigen::DynamicSparseMatrix<ScalarType> dynamic_weight_matrix(end-begin,end-begin);
      69                 :            :         dynamic_weight_matrix.reserve(sparse_triplets.size());
      70                 :            :         for (SparseTriplets::const_iterator it=sparse_triplets.begin(); it!=sparse_triplets.end(); ++it)
      71                 :            :                 dynamic_weight_matrix.coeffRef(it->col(),it->row()) += it->value();
      72                 :            :         SparseWeightMatrix weight_matrix(dynamic_weight_matrix);
      73                 :            : #else
      74 [ +  - ][ +  - ]:          3 :         SparseWeightMatrix weight_matrix(end-begin,end-begin);
         [ #  # ][ #  # ]
      75 [ +  - ][ +  - ]:          3 :         weight_matrix.setFromTriplets(sparse_triplets.begin(),sparse_triplets.end());
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      76                 :            : #endif
      77                 :            : 
      78 [ +  - ][ +  - ]:          3 :         return Laplacian(weight_matrix,DenseDiagonalMatrix(D));
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      79                 :            : }
      80                 :            : 
      81                 :            : template<class RandomAccessIterator, class FeatureVectorCallback>
      82                 :          1 : DenseSymmetricMatrixPair construct_locality_preserving_eigenproblem(SparseWeightMatrix& L,
      83                 :            :                 DenseDiagonalMatrix& D, RandomAccessIterator begin, RandomAccessIterator end, FeatureVectorCallback feature_vector_callback,
      84                 :            :                 IndexType dimension)
      85                 :            : {
      86 [ +  - ][ +  - ]:          1 :         timed_context context("Constructing LPP eigenproblem");
         [ +  - ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
      87                 :            : 
      88         [ +  - ]:          1 :         DenseSymmetricMatrix lhs = DenseSymmetricMatrix::Zero(dimension,dimension);
           [ +  -  #  # ]
           [ #  #  #  # ]
                 [ #  # ]
      89 [ +  - ][ +  - ]:          1 :         DenseSymmetricMatrix rhs = DenseSymmetricMatrix::Zero(dimension,dimension);
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      90                 :            : 
      91 [ +  - ][ #  # ]:          1 :         DenseVector rank_update_vector_i(dimension);
                 [ #  # ]
      92 [ +  - ][ #  # ]:          1 :         DenseVector rank_update_vector_j(dimension);
                 [ #  # ]
      93 [ +  + ][ #  # ]:         51 :         for (RandomAccessIterator iter=begin; iter!=end; ++iter)
                 [ #  # ]
      94                 :            :         {
      95 [ +  - ][ #  # ]:         50 :                 feature_vector_callback.vector(*iter,rank_update_vector_i);
                 [ #  # ]
      96 [ +  - ][ +  - ]:         50 :                 rhs.selfadjointView<Eigen::Upper>().rankUpdate(rank_update_vector_i,D.diagonal()(iter-begin));
         [ +  - ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
      97                 :            :         }
      98                 :            : 
      99 [ +  + ][ #  # ]:         51 :         for (int i=0; i<L.outerSize(); ++i)
                 [ #  # ]
     100                 :            :         {
     101 [ +  - ][ +  + ]:        702 :                 for (SparseWeightMatrix::InnerIterator it(L,i); it; ++it)
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
     102                 :            :                 {
     103 [ +  - ][ #  # ]:        652 :                         feature_vector_callback.vector(begin[it.row()],rank_update_vector_i);
                 [ #  # ]
     104 [ +  - ][ #  # ]:        652 :                         feature_vector_callback.vector(begin[it.col()],rank_update_vector_j);
                 [ #  # ]
     105 [ +  - ][ +  - ]:        652 :                         lhs.selfadjointView<Eigen::Upper>().rankUpdate(rank_update_vector_i, rank_update_vector_j, it.value());
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
     106                 :            :                 }
     107                 :            :         }
     108                 :            : 
     109 [ +  - ][ +  - ]:          1 :         return DenseSymmetricMatrixPair(lhs,rhs);
         [ +  - ][ +  - ]
         [ +  - ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
     110                 :            : }
     111                 :            : 
     112                 :            : }
     113                 :            : }
     114                 :            : 
     115                 :            : #endif

Generated by: LCOV version 1.9