LCOV - code coverage report
Current view: top level - routines - isomap.hpp (source / functions) Hit Total Coverage
Test: clean.info Lines: 79 79 100.0 %
Date: 2013-05-24 Functions: 5 21 23.8 %
Branches: 51 256 19.9 %

           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_Isomap_H_
       7                 :            : #define TAPKEE_Isomap_H_
       8                 :            : 
       9                 :            : /* Tapkee includes */
      10                 :            : #include <tapkee/defines.hpp>
      11                 :            : #include <tapkee/utils/fibonacci_heap.hpp>
      12                 :            : #include <tapkee/utils/reservable_priority_queue.hpp>
      13                 :            : #include <tapkee/utils/time.hpp>
      14                 :            : /* End of Tapkee includes */
      15                 :            : 
      16                 :            : #include <limits>
      17                 :            : 
      18                 :            : namespace tapkee
      19                 :            : {
      20                 :            : namespace tapkee_internal
      21                 :            : {
      22                 :            : 
      23                 :            : #ifdef TAPKEE_USE_PRIORITY_QUEUE
      24                 :            : typedef std::pair<IndexType,ScalarType> HeapElement;
      25                 :            : 
      26                 :            : struct HeapElementComparator
      27                 :            : {
      28                 :      24097 :         inline bool operator()(const HeapElement& l, const HeapElement& r) const 
      29                 :            :         {
      30                 :      24097 :                 return l.second > r.second;
      31                 :            :         }
      32                 :            : };
      33                 :            : #endif
      34                 :            : 
      35                 :            : //! Computes shortest distances (so-called geodesic distances)
      36                 :            : //! using Dijkstra algorithm.
      37                 :            : //!
      38                 :            : //! @param begin begin data iterator
      39                 :            : //! @param end end data iterator
      40                 :            : //! @param neighbors neighbors of each vector
      41                 :            : //! @param callback distance callback
      42                 :            : //!
      43                 :            : template <class RandomAccessIterator, class DistanceCallback>
      44                 :          1 : DenseSymmetricMatrix compute_shortest_distances_matrix(const RandomAccessIterator& begin, const RandomAccessIterator& end, 
      45                 :            :                 const Neighbors& neighbors, DistanceCallback callback)
      46                 :            : {
      47 [ +  - ][ +  - ]:          1 :         timed_context context("Distances shortest path relaxing");
         [ +  - ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      48                 :          1 :         const IndexType n_neighbors = neighbors[0].size();
      49                 :          1 :         const IndexType N = (end-begin);
      50                 :            : 
      51   [ +  -  #  #  :          1 :         DenseSymmetricMatrix shortest_distances(N,N);
             #  #  #  # ]
      52                 :            :         
      53                 :            : #pragma omp parallel shared(shortest_distances,neighbors,begin,callback) default(none)
      54                 :            :         {
      55 [ +  + ][ #  # ]:        154 :                 bool* f = new bool[N];
         [ #  # ][ #  # ]
      56 [ +  + ][ #  # ]:        168 :                 bool* s = new bool[N];
         [ #  # ][ #  # ]
      57                 :            :                 IndexType k;
      58                 :            : 
      59                 :            : #ifdef TAPKEE_USE_PRIORITY_QUEUE
      60                 :          4 :                 reservable_priority_queue<HeapElement,HeapElementComparator> heap(N);
      61                 :            : #else
      62                 :            :                 fibonacci_heap heap(N);
      63                 :            : #endif
      64                 :            : 
      65 [ +  + ][ #  # ]:          7 : #pragma omp for nowait
         [ #  # ][ #  # ]
      66                 :            :                 for (k=0; k<N; k++)
      67                 :            :                 {
      68                 :            :                         // fill s and f with false, fill shortest_D with infinity
      69 [ +  + ][ #  # ]:       2545 :                         for (IndexType j=0; j<N; j++)
         [ #  # ][ #  # ]
      70                 :            :                         {
      71                 :       2491 :                                 shortest_distances(k,j) = std::numeric_limits<DenseMatrix::Scalar>::max();
      72                 :       2488 :                                 s[j] = false;
      73                 :       2488 :                                 f[j] = false;
      74                 :            :                         }
      75                 :            :                         // set distance from k to k as zero
      76                 :         50 :                         shortest_distances(k,k) = 0.0;
      77                 :            : 
      78                 :            :                         // insert kth object to heap with zero distance and set f[k] true
      79                 :            : #ifdef TAPKEE_USE_PRIORITY_QUEUE
      80                 :         50 :                         HeapElement heap_element_of_self(k,0.0);
      81                 :         50 :                         heap.push(heap_element_of_self);
      82                 :            : #else
      83                 :            :                         heap.insert(k,0.0);
      84                 :            : #endif
      85                 :         50 :                         f[k] = true;
      86                 :            : 
      87                 :            :                         // while heap is not empty
      88 [ +  + ][ #  # ]:       3639 :                         while (!heap.empty())
         [ #  # ][ #  # ]
      89                 :            :                         {
      90                 :            :                                 // extract min and set (s)olution state as true and (f)rontier as false
      91                 :            : #ifdef TAPKEE_USE_PRIORITY_QUEUE
      92                 :       3556 :                                 int min_item = heap.top().first;
      93                 :       3550 :                                 ScalarType min_item_d = heap.top().second;
      94                 :       3559 :                                 heap.pop();
      95   [ +  +  #  #  :       3564 :                                 if (min_item_d > shortest_distances(k,min_item))
             #  #  #  # ]
      96                 :       1064 :                                         continue;
      97                 :            : #else
      98                 :            :                                 ScalarType tmp;
      99                 :            :                                 int min_item = heap.extract_min(tmp);
     100                 :            : #endif
     101                 :            : 
     102                 :       2493 :                                 s[min_item] = true;
     103                 :       2493 :                                 f[min_item] = false;
     104                 :            : 
     105                 :            :                                 // for-each edge (min_item->w)
     106 [ +  + ][ #  # ]:      26208 :                                 for (IndexType i=0; i<n_neighbors; i++)
         [ #  # ][ #  # ]
     107                 :            :                                 {
     108                 :            :                                         // get w idx
     109                 :      23683 :                                         int w = neighbors[min_item][i];
     110                 :            :                                         // if w is not in solution yet
     111   [ +  +  #  #  :      24622 :                                         if (s[w] == false)
             #  #  #  # ]
     112                 :            :                                         {
     113                 :            :                                                 // get distance from k to i through min_item
     114                 :      11087 :                                                 ScalarType dist = shortest_distances(k,min_item) + callback.distance(begin[min_item],begin[w]);
     115                 :            :                                                 // if distance can be relaxed
     116   [ +  +  #  #  :       9729 :                                                 if (dist < shortest_distances(k,w))
             #  #  #  # ]
     117                 :            :                                                 {
     118                 :            :                                                         // relax distance
     119                 :       3410 :                                                         shortest_distances(k,w) = dist;
     120                 :            : #ifdef TAPKEE_USE_PRIORITY_QUEUE
     121                 :       3466 :                                                         HeapElement relaxed_heap_element(w,dist);
     122                 :       3434 :                                                         heap.push(relaxed_heap_element);
     123                 :       3493 :                                                         f[w] = true;
     124                 :            : #else
     125                 :            :                                                         // if w is in (f)rontier
     126                 :            :                                                         if (f[w])
     127                 :            :                                                         {
     128                 :            :                                                                 // decrease distance in heap
     129                 :            :                                                                 heap.decrease_key(w, dist);
     130                 :            :                                                         }
     131                 :            :                                                         else
     132                 :            :                                                         {
     133                 :            :                                                                 // insert w to heap and set (f)rontier as true
     134                 :            :                                                                 heap.insert(w, dist);
     135                 :            :                                                                 f[w] = true;
     136                 :            :                                                         }
     137                 :            : #endif
     138                 :            :                                                 }
     139                 :            :                                         }
     140                 :            :                                 }
     141                 :            :                         }
     142 [ +  - ][ #  # ]:         50 :                         heap.clear();
         [ #  # ][ #  # ]
     143                 :            :                 }
     144                 :            : 
     145 [ -  + ][ #  # ]:          4 :                 delete[] s;
         [ #  # ][ #  # ]
     146 [ -  + ][ #  # ]:          8 :                 delete[] f;
         [ #  # ][ #  # ]
     147                 :            :         }
     148                 :          1 :         return shortest_distances;
     149                 :            : }
     150                 :            : 
     151                 :            : //! Computes shortest distances (so-called geodesic distances)
     152                 :            : //! using Dijkstra algorithm with landmarks.
     153                 :            : //!
     154                 :            : //! @param begin begin data iterator
     155                 :            : //! @param end end data iterator
     156                 :            : //! @param landmarks landmarks
     157                 :            : //! @param neighbors neighbors of each vector
     158                 :            : //! @param callback distance callback
     159                 :            : //!
     160                 :            : template <class RandomAccessIterator, class DistanceCallback>
     161                 :          1 : DenseMatrix compute_shortest_distances_matrix(const RandomAccessIterator& begin, const RandomAccessIterator& end, 
     162                 :            :                 const Landmarks& landmarks, const Neighbors& neighbors, DistanceCallback callback)
     163                 :            : {
     164 [ +  - ][ +  - ]:          1 :         timed_context context("Distances shortest path relaxing");
         [ +  - ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
     165                 :          1 :         const IndexType n_neighbors = neighbors[0].size();
     166                 :          1 :         const IndexType N = end-begin;
     167                 :          1 :         const IndexType N_landmarks = landmarks.size();
     168                 :            : 
     169   [ +  -  #  #  :          1 :         DenseMatrix shortest_distances(landmarks.size(),N);
             #  #  #  # ]
     170                 :            :         
     171                 :            : #pragma omp parallel shared(shortest_distances,begin,landmarks,neighbors,callback) default(none)
     172                 :            :         {
     173 [ +  + ][ #  # ]:        167 :                 bool* f = new bool[N];
         [ #  # ][ #  # ]
     174 [ +  + ][ #  # ]:        171 :                 bool* s = new bool[N];
         [ #  # ][ #  # ]
     175                 :            :                 IndexType k;
     176                 :            : 
     177                 :            : #ifdef TAPKEE_USE_PRIORITY_QUEUE
     178                 :          4 :                 reservable_priority_queue<HeapElement,HeapElementComparator> heap(N);
     179                 :            : #else
     180                 :            :                 fibonacci_heap heap(N);
     181                 :            : #endif
     182                 :            : 
     183 [ +  + ][ #  # ]:          5 : #pragma omp for nowait
         [ #  # ][ #  # ]
     184                 :            :                 for (k=0; k<N_landmarks; k++)
     185                 :            :                 {
     186                 :            :                         // fill s and f with false, fill shortest_D with infinity
     187 [ +  + ][ #  # ]:       1267 :                         for (IndexType j=0; j<N; j++)
         [ #  # ][ #  # ]
     188                 :            :                         {
     189                 :       1238 :                                 shortest_distances(k,j) = std::numeric_limits<DenseMatrix::Scalar>::max();
     190                 :       1237 :                                 s[j] = false;
     191                 :       1237 :                                 f[j] = false;
     192                 :            :                         }
     193                 :            :                         // set distance from k to k as zero
     194                 :         25 :                         shortest_distances(k,landmarks[k]) = 0.0;
     195                 :            : 
     196                 :            :                         // insert kth object to heap with zero distance and set f[k] true
     197                 :            : #ifdef TAPKEE_USE_PRIORITY_QUEUE
     198                 :         25 :                         HeapElement heap_element_of_self(landmarks[k],0.0);
     199                 :         25 :                         heap.push(heap_element_of_self);
     200                 :            : #else
     201                 :            :                         heap.insert(landmarks[k],0.0);
     202                 :            : #endif
     203                 :         25 :                         f[k] = true;
     204                 :            : 
     205                 :            :                         // while heap is not empty
     206 [ +  + ][ #  # ]:       1787 :                         while (!heap.empty())
         [ #  # ][ #  # ]
     207                 :            :                         {
     208                 :            :                                 // extract min and set (s)olution state as true and (f)rontier as false
     209                 :            : #ifdef TAPKEE_USE_PRIORITY_QUEUE
     210                 :       1854 :                                 int min_item = heap.top().first;
     211                 :       1852 :                                 ScalarType min_item_d = heap.top().second;
     212                 :       1852 :                                 heap.pop();
     213   [ +  +  #  #  :       1854 :                                 if (min_item_d > shortest_distances(k,min_item))
             #  #  #  # ]
     214                 :        607 :                                         continue;
     215                 :            : #else
     216                 :            :                                 ScalarType tmp;
     217                 :            :                                 int min_item = heap.extract_min(tmp);
     218                 :            : #endif
     219                 :            : 
     220                 :       1246 :                                 s[min_item] = true;
     221                 :       1246 :                                 f[min_item] = false;
     222                 :            : 
     223                 :            :                                 // for-each edge (min_item->w)
     224 [ +  + ][ #  # ]:      13000 :                                 for (IndexType i=0; i<n_neighbors; i++)
         [ #  # ][ #  # ]
     225                 :            :                                 {
     226                 :            :                                         // get w idx
     227                 :      11845 :                                         int w = neighbors[min_item][i];
     228                 :            :                                         // if w is not in solution yet
     229   [ +  +  #  #  :      12219 :                                         if (s[w] == false)
             #  #  #  # ]
     230                 :            :                                         {
     231                 :            :                                                 // get distance from k to i through min_item
     232                 :       5493 :                                                 ScalarType dist = shortest_distances(k,min_item) + callback.distance(begin[min_item],begin[w]);
     233                 :            :                                                 // if distance can be relaxed
     234   [ +  +  #  #  :       4746 :                                                 if (dist < shortest_distances(k,w))
             #  #  #  # ]
     235                 :            :                                                 {
     236                 :            :                                                         // relax distance
     237                 :       1750 :                                                         shortest_distances(k,w) = dist;
     238                 :            : #ifdef TAPKEE_USE_PRIORITY_QUEUE
     239                 :       1792 :                                                         HeapElement relaxed_heap_element(w,dist);
     240                 :       1764 :                                                         heap.push(relaxed_heap_element);
     241                 :       1822 :                                                         f[w] = true;
     242                 :            : #else
     243                 :            :                                                         // if w is in (f)rontier
     244                 :            :                                                         if (f[w])
     245                 :            :                                                         {
     246                 :            :                                                                 // decrease distance in heap
     247                 :            :                                                                 heap.decrease_key(w, dist);
     248                 :            :                                                         }
     249                 :            :                                                         else
     250                 :            :                                                         {
     251                 :            :                                                                 // insert w to heap and set (f)rontier as true
     252                 :            :                                                                 heap.insert(w, dist);
     253                 :            :                                                                 f[w] = true;
     254                 :            :                                                         }
     255                 :            : #endif
     256                 :            :                                                 }
     257                 :            :                                         }
     258                 :            :                                 }
     259                 :            :                         }
     260 [ +  + ][ #  # ]:         25 :                         heap.clear();
         [ #  # ][ #  # ]
     261                 :            :                 }
     262                 :            : 
     263 [ -  + ][ #  # ]:          4 :                 delete[] s;
         [ #  # ][ #  # ]
     264 [ -  + ][ #  # ]:          8 :                 delete[] f;
         [ #  # ][ #  # ]
     265                 :            :         }
     266                 :          1 :         return shortest_distances;
     267                 :            : }
     268                 :            : 
     269                 :            : }
     270                 :            : }
     271                 :            : 
     272                 :            : #endif

Generated by: LCOV version 1.9