Tapkee
isomap.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_Isomap_H_
00007 #define TAPKEE_Isomap_H_
00008 
00009 /* Tapkee includes */
00010 #include <tapkee/defines.hpp>
00011 #include <tapkee/utils/fibonacci_heap.hpp>
00012 #include <tapkee/utils/reservable_priority_queue.hpp>
00013 #include <tapkee/utils/time.hpp>
00014 /* End of Tapkee includes */
00015 
00016 #include <limits>
00017 
00018 namespace tapkee
00019 {
00020 namespace tapkee_internal
00021 {
00022 
00023 #ifdef TAPKEE_USE_PRIORITY_QUEUE
00024 typedef std::pair<IndexType,ScalarType> HeapElement;
00025 
00026 struct HeapElementComparator
00027 {
00028     inline bool operator()(const HeapElement& l, const HeapElement& r) const 
00029     {
00030         return l.second > r.second;
00031     }
00032 };
00033 #endif
00034 
00043 template <class RandomAccessIterator, class DistanceCallback>
00044 DenseSymmetricMatrix compute_shortest_distances_matrix(const RandomAccessIterator& begin, const RandomAccessIterator& end, 
00045         const Neighbors& neighbors, DistanceCallback callback)
00046 {
00047     timed_context context("Distances shortest path relaxing");
00048     const IndexType n_neighbors = neighbors[0].size();
00049     const IndexType N = (end-begin);
00050 
00051     DenseSymmetricMatrix shortest_distances(N,N);
00052     
00053 #pragma omp parallel shared(shortest_distances,neighbors,begin,callback) default(none)
00054     {
00055         bool* f = new bool[N];
00056         bool* s = new bool[N];
00057         IndexType k;
00058 
00059 #ifdef TAPKEE_USE_PRIORITY_QUEUE
00060         reservable_priority_queue<HeapElement,HeapElementComparator> heap(N);
00061 #else
00062         fibonacci_heap heap(N);
00063 #endif
00064 
00065 #pragma omp for nowait
00066         for (k=0; k<N; k++)
00067         {
00068             // fill s and f with false, fill shortest_D with infinity
00069             for (IndexType j=0; j<N; j++)
00070             {
00071                 shortest_distances(k,j) = std::numeric_limits<DenseMatrix::Scalar>::max();
00072                 s[j] = false;
00073                 f[j] = false;
00074             }
00075             // set distance from k to k as zero
00076             shortest_distances(k,k) = 0.0;
00077 
00078             // insert kth object to heap with zero distance and set f[k] true
00079 #ifdef TAPKEE_USE_PRIORITY_QUEUE
00080             HeapElement heap_element_of_self(k,0.0);
00081             heap.push(heap_element_of_self);
00082 #else
00083             heap.insert(k,0.0);
00084 #endif
00085             f[k] = true;
00086 
00087             // while heap is not empty
00088             while (!heap.empty())
00089             {
00090                 // extract min and set (s)olution state as true and (f)rontier as false
00091 #ifdef TAPKEE_USE_PRIORITY_QUEUE
00092                 int min_item = heap.top().first;
00093                 ScalarType min_item_d = heap.top().second;
00094                 heap.pop();
00095                 if (min_item_d > shortest_distances(k,min_item))
00096                     continue;
00097 #else
00098                 ScalarType tmp;
00099                 int min_item = heap.extract_min(tmp);
00100 #endif
00101 
00102                 s[min_item] = true;
00103                 f[min_item] = false;
00104 
00105                 // for-each edge (min_item->w)
00106                 for (IndexType i=0; i<n_neighbors; i++)
00107                 {
00108                     // get w idx
00109                     int w = neighbors[min_item][i];
00110                     // if w is not in solution yet
00111                     if (s[w] == false)
00112                     {
00113                         // get distance from k to i through min_item
00114                         ScalarType dist = shortest_distances(k,min_item) + callback.distance(begin[min_item],begin[w]);
00115                         // if distance can be relaxed
00116                         if (dist < shortest_distances(k,w))
00117                         {
00118                             // relax distance
00119                             shortest_distances(k,w) = dist;
00120 #ifdef TAPKEE_USE_PRIORITY_QUEUE
00121                             HeapElement relaxed_heap_element(w,dist);
00122                             heap.push(relaxed_heap_element);
00123                             f[w] = true;
00124 #else
00125                             // if w is in (f)rontier
00126                             if (f[w])
00127                             {
00128                                 // decrease distance in heap
00129                                 heap.decrease_key(w, dist);
00130                             }
00131                             else
00132                             {
00133                                 // insert w to heap and set (f)rontier as true
00134                                 heap.insert(w, dist);
00135                                 f[w] = true;
00136                             }
00137 #endif
00138                         }
00139                     }
00140                 }
00141             }
00142             heap.clear();
00143         }
00144 
00145         delete[] s;
00146         delete[] f;
00147     }
00148     return shortest_distances;
00149 }
00150 
00160 template <class RandomAccessIterator, class DistanceCallback>
00161 DenseMatrix compute_shortest_distances_matrix(const RandomAccessIterator& begin, const RandomAccessIterator& end, 
00162         const Landmarks& landmarks, const Neighbors& neighbors, DistanceCallback callback)
00163 {
00164     timed_context context("Distances shortest path relaxing");
00165     const IndexType n_neighbors = neighbors[0].size();
00166     const IndexType N = end-begin;
00167     const IndexType N_landmarks = landmarks.size();
00168 
00169     DenseMatrix shortest_distances(landmarks.size(),N);
00170     
00171 #pragma omp parallel shared(shortest_distances,begin,landmarks,neighbors,callback) default(none)
00172     {
00173         bool* f = new bool[N];
00174         bool* s = new bool[N];
00175         IndexType k;
00176 
00177 #ifdef TAPKEE_USE_PRIORITY_QUEUE
00178         reservable_priority_queue<HeapElement,HeapElementComparator> heap(N);
00179 #else
00180         fibonacci_heap heap(N);
00181 #endif
00182 
00183 #pragma omp for nowait
00184         for (k=0; k<N_landmarks; k++)
00185         {
00186             // fill s and f with false, fill shortest_D with infinity
00187             for (IndexType j=0; j<N; j++)
00188             {
00189                 shortest_distances(k,j) = std::numeric_limits<DenseMatrix::Scalar>::max();
00190                 s[j] = false;
00191                 f[j] = false;
00192             }
00193             // set distance from k to k as zero
00194             shortest_distances(k,landmarks[k]) = 0.0;
00195 
00196             // insert kth object to heap with zero distance and set f[k] true
00197 #ifdef TAPKEE_USE_PRIORITY_QUEUE
00198             HeapElement heap_element_of_self(landmarks[k],0.0);
00199             heap.push(heap_element_of_self);
00200 #else
00201             heap.insert(landmarks[k],0.0);
00202 #endif
00203             f[k] = true;
00204 
00205             // while heap is not empty
00206             while (!heap.empty())
00207             {
00208                 // extract min and set (s)olution state as true and (f)rontier as false
00209 #ifdef TAPKEE_USE_PRIORITY_QUEUE
00210                 int min_item = heap.top().first;
00211                 ScalarType min_item_d = heap.top().second;
00212                 heap.pop();
00213                 if (min_item_d > shortest_distances(k,min_item))
00214                     continue;
00215 #else
00216                 ScalarType tmp;
00217                 int min_item = heap.extract_min(tmp);
00218 #endif
00219 
00220                 s[min_item] = true;
00221                 f[min_item] = false;
00222 
00223                 // for-each edge (min_item->w)
00224                 for (IndexType i=0; i<n_neighbors; i++)
00225                 {
00226                     // get w idx
00227                     int w = neighbors[min_item][i];
00228                     // if w is not in solution yet
00229                     if (s[w] == false)
00230                     {
00231                         // get distance from k to i through min_item
00232                         ScalarType dist = shortest_distances(k,min_item) + callback.distance(begin[min_item],begin[w]);
00233                         // if distance can be relaxed
00234                         if (dist < shortest_distances(k,w))
00235                         {
00236                             // relax distance
00237                             shortest_distances(k,w) = dist;
00238 #ifdef TAPKEE_USE_PRIORITY_QUEUE
00239                             HeapElement relaxed_heap_element(w,dist);
00240                             heap.push(relaxed_heap_element);
00241                             f[w] = true;
00242 #else
00243                             // if w is in (f)rontier
00244                             if (f[w])
00245                             {
00246                                 // decrease distance in heap
00247                                 heap.decrease_key(w, dist);
00248                             }
00249                             else
00250                             {
00251                                 // insert w to heap and set (f)rontier as true
00252                                 heap.insert(w, dist);
00253                                 f[w] = true;
00254                             }
00255 #endif
00256                         }
00257                     }
00258                 }
00259             }
00260             heap.clear();
00261         }
00262 
00263         delete[] s;
00264         delete[] f;
00265     }
00266     return shortest_distances;
00267 }
00268 
00269 }
00270 }
00271 
00272 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines