Tapkee
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
isomap.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_Isomap_H_
7 #define TAPKEE_Isomap_H_
8 
9 /* Tapkee includes */
10 #include <tapkee/defines.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  inline bool operator()(const HeapElement& l, const HeapElement& r) const
29  {
30  return l.second > r.second;
31  }
32 };
33 #endif
34 
43 template <class RandomAccessIterator, class DistanceCallback>
44 DenseSymmetricMatrix compute_shortest_distances_matrix(const RandomAccessIterator& begin, const RandomAccessIterator& end,
45  const Neighbors& neighbors, DistanceCallback callback)
46 {
47  timed_context context("Distances shortest path relaxing");
48  const IndexType n_neighbors = neighbors[0].size();
49  const IndexType N = (end-begin);
50 
51  DenseSymmetricMatrix shortest_distances(N,N);
52 
53 #pragma omp parallel shared(shortest_distances,neighbors,begin,callback) default(none)
54  {
55  bool* f = new bool[N];
56  bool* s = new bool[N];
57  IndexType k;
58 
59 #ifdef TAPKEE_USE_PRIORITY_QUEUE
61 #else
62  fibonacci_heap heap(N);
63 #endif
64 
65 #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  for (IndexType j=0; j<N; j++)
70  {
71  shortest_distances(k,j) = std::numeric_limits<DenseMatrix::Scalar>::max();
72  s[j] = false;
73  f[j] = false;
74  }
75  // set distance from k to k as zero
76  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  HeapElement heap_element_of_self(k,0.0);
81  heap.push(heap_element_of_self);
82 #else
83  heap.insert(k,0.0);
84 #endif
85  f[k] = true;
86 
87  // while heap is not empty
88  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  int min_item = heap.top().first;
93  ScalarType min_item_d = heap.top().second;
94  heap.pop();
95  if (min_item_d > shortest_distances(k,min_item))
96  continue;
97 #else
98  ScalarType tmp;
99  int min_item = heap.extract_min(tmp);
100 #endif
101 
102  s[min_item] = true;
103  f[min_item] = false;
104 
105  // for-each edge (min_item->w)
106  for (IndexType i=0; i<n_neighbors; i++)
107  {
108  // get w idx
109  int w = neighbors[min_item][i];
110  // if w is not in solution yet
111  if (s[w] == false)
112  {
113  // get distance from k to i through min_item
114  ScalarType dist = shortest_distances(k,min_item) + callback.distance(begin[min_item],begin[w]);
115  // if distance can be relaxed
116  if (dist < shortest_distances(k,w))
117  {
118  // relax distance
119  shortest_distances(k,w) = dist;
120 #ifdef TAPKEE_USE_PRIORITY_QUEUE
121  HeapElement relaxed_heap_element(w,dist);
122  heap.push(relaxed_heap_element);
123  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  heap.clear();
143  }
144 
145  delete[] s;
146  delete[] f;
147  }
148  return shortest_distances;
149 }
150 
160 template <class RandomAccessIterator, class DistanceCallback>
161 DenseMatrix compute_shortest_distances_matrix(const RandomAccessIterator& begin, const RandomAccessIterator& end,
162  const Landmarks& landmarks, const Neighbors& neighbors, DistanceCallback callback)
163 {
164  timed_context context("Distances shortest path relaxing");
165  const IndexType n_neighbors = neighbors[0].size();
166  const IndexType N = end-begin;
167  const IndexType N_landmarks = landmarks.size();
168 
169  DenseMatrix shortest_distances(landmarks.size(),N);
170 
171 #pragma omp parallel shared(shortest_distances,begin,landmarks,neighbors,callback) default(none)
172  {
173  bool* f = new bool[N];
174  bool* s = new bool[N];
175  IndexType k;
176 
177 #ifdef TAPKEE_USE_PRIORITY_QUEUE
179 #else
180  fibonacci_heap heap(N);
181 #endif
182 
183 #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  for (IndexType j=0; j<N; j++)
188  {
189  shortest_distances(k,j) = std::numeric_limits<DenseMatrix::Scalar>::max();
190  s[j] = false;
191  f[j] = false;
192  }
193  // set distance from k to k as zero
194  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  HeapElement heap_element_of_self(landmarks[k],0.0);
199  heap.push(heap_element_of_self);
200 #else
201  heap.insert(landmarks[k],0.0);
202 #endif
203  f[k] = true;
204 
205  // while heap is not empty
206  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  int min_item = heap.top().first;
211  ScalarType min_item_d = heap.top().second;
212  heap.pop();
213  if (min_item_d > shortest_distances(k,min_item))
214  continue;
215 #else
216  ScalarType tmp;
217  int min_item = heap.extract_min(tmp);
218 #endif
219 
220  s[min_item] = true;
221  f[min_item] = false;
222 
223  // for-each edge (min_item->w)
224  for (IndexType i=0; i<n_neighbors; i++)
225  {
226  // get w idx
227  int w = neighbors[min_item][i];
228  // if w is not in solution yet
229  if (s[w] == false)
230  {
231  // get distance from k to i through min_item
232  ScalarType dist = shortest_distances(k,min_item) + callback.distance(begin[min_item],begin[w]);
233  // if distance can be relaxed
234  if (dist < shortest_distances(k,w))
235  {
236  // relax distance
237  shortest_distances(k,w) = dist;
238 #ifdef TAPKEE_USE_PRIORITY_QUEUE
239  HeapElement relaxed_heap_element(w,dist);
240  heap.push(relaxed_heap_element);
241  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  heap.clear();
261  }
262 
263  delete[] s;
264  delete[] f;
265  }
266  return shortest_distances;
267 }
268 
269 }
270 }
271 
272 #endif
Eigen::Matrix< tapkee::ScalarType, Eigen::Dynamic, Eigen::Dynamic > DenseMatrix
dense matrix type (non-overridable)
Definition: types.hpp:23
double ScalarType
default scalar value (can be overrided with TAPKEE_CUSTOM_INTERNAL_NUMTYPE define) ...
Definition: types.hpp:15
DenseSymmetricMatrix compute_shortest_distances_matrix(const RandomAccessIterator &begin, const RandomAccessIterator &end, const Neighbors &neighbors, DistanceCallback callback)
Computes shortest distances (so-called geodesic distances) using Dijkstra algorithm.
Definition: isomap.hpp:44
int IndexType
indexing type (non-overridable) set to int for compatibility with OpenMP 2.0
Definition: types.hpp:19
void decrease_key(int index, ScalarType &key)
void insert(int index, ScalarType key)
TAPKEE_INTERNAL_VECTOR< tapkee::tapkee_internal::LocalNeighbors > Neighbors
Definition: synonyms.hpp:40
the class fibonacci_heap, a fibonacci heap. Generally used by Isomap for Dijkstra heap algorithm ...
TAPKEE_INTERNAL_VECTOR< tapkee::IndexType > Landmarks
Definition: synonyms.hpp:42
tapkee::DenseMatrix DenseSymmetricMatrix
dense symmetric matrix (non-overridable, currently just dense matrix, can be improved later) ...
Definition: types.hpp:25