SHOGUN  v1.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Isomap.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2011 Sergey Lisitsyn
8  * Copyright (C) 2011 Berlin Institute of Technology and Max-Planck-Society
9  */
10 
12 #ifdef HAVE_LAPACK
13 #include <shogun/lib/common.h>
16 #include <shogun/io/SGIO.h>
17 #include <shogun/base/Parallel.h>
18 #include <shogun/lib/Signal.h>
19 
20 #ifdef HAVE_PTHREAD
21 #include <pthread.h>
22 #endif
23 
24 using namespace shogun;
25 
26 #ifndef DOXYGEN_SHOULD_SKIP_THIS
27 /* struct storing thread params
28  */
29 struct DIJKSTRA_THREAD_PARAM
30 {
32  CFibonacciHeap* heap;
35  const float64_t* edges_matrix;
38  const int32_t* edges_idx_matrix;
40  float64_t* shortest_D;
42  int32_t i_start;
44  int32_t i_stop;
46  int32_t i_step;
48  int32_t m_k;
50  bool* s;
52  bool* f;
53 };
54 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
55 
57 {
58  m_k = 3;
59 
60  init();
61 }
62 
64 {
65  m_parameters->add(&m_k, "k", "number of neighbors");
66 }
67 
69 {
70 }
71 
72 void CIsomap::set_k(int32_t k)
73 {
74  ASSERT(k>0);
75  m_k = k;
76 }
77 
78 int32_t CIsomap::get_k() const
79 {
80  return m_k;
81 }
82 
83 const char* CIsomap::get_name() const
84 {
85  return "Isomap";
86 }
87 
89 {
90  return isomap_distance(distance_matrix);
91 }
92 
94 {
95  int32_t N,t,i,j;
96  float64_t tmp;
97  N = D_matrix.num_cols;
98  if (D_matrix.num_cols!=D_matrix.num_rows)
99  {
100  D_matrix.destroy_matrix();
101  SG_ERROR("Given distance matrix is not square.\n");
102  }
103  if (m_k>=N)
104  {
105  D_matrix.destroy_matrix();
106  SG_ERROR("K parameter should be less than number of given vectors (k=%d, N=%d)\n", m_k, N);
107  }
108 
109  // cut by k-nearest neighbors
110  int32_t* edges_idx_matrix = SG_MALLOC(int32_t, N*m_k);
111  float64_t* edges_matrix = SG_MALLOC(float64_t, N*m_k);
112 
113  // query neighbors and edges to neighbors
114  CFibonacciHeap* heap = new CFibonacciHeap(N);
115  for (i=0; i<N; i++)
116  {
117  // insert distances to heap
118  for (j=0; j<N; j++)
119  heap->insert(j,D_matrix[i*N+j]);
120 
121  // extract nearest neighbor: the jth object itself
122  heap->extract_min(tmp);
123 
124  // extract m_k neighbors and distances
125  for (j=0; j<m_k; j++)
126  {
127  edges_idx_matrix[i*m_k+j] = heap->extract_min(tmp);
128  edges_matrix[i*m_k+j] = tmp;
129  }
130  // clear heap
131  heap->clear();
132  }
133  // cleanup
134  delete heap;
135 
136 #ifdef HAVE_PTHREAD
137 
138  // Parallel Dijkstra with Fibonacci Heap
139  int32_t num_threads = parallel->get_num_threads();
140  ASSERT(num_threads>0);
141  // allocate threads and thread parameters
142  pthread_t* threads = SG_MALLOC(pthread_t, num_threads);
143  DIJKSTRA_THREAD_PARAM* parameters = SG_MALLOC(DIJKSTRA_THREAD_PARAM, num_threads);
144  // allocate heaps
145  CFibonacciHeap** heaps = SG_MALLOC(CFibonacciHeap*, num_threads);
146  for (t=0; t<num_threads; t++)
147  heaps[t] = new CFibonacciHeap(N);
148 
149 #else
150  int32_t num_threads = 1;
151 #endif
152 
153  // allocate (s)olution
154  bool* s = SG_MALLOC(bool,N*num_threads);
155  // allocate (f)rontier
156  bool* f = SG_MALLOC(bool,N*num_threads);
157  // init matrix to store shortest distances
158  float64_t* shortest_D = D_matrix.matrix;
159 
160 #ifdef HAVE_PTHREAD
161 
162  pthread_attr_t attr;
163  pthread_attr_init(&attr);
164  pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
165  for (t=0; t<num_threads; t++)
166  {
167  parameters[t].i_start = t;
168  parameters[t].i_stop = N;
169  parameters[t].i_step = num_threads;
170  parameters[t].heap = heaps[t];
171  parameters[t].edges_matrix = edges_matrix;
172  parameters[t].edges_idx_matrix = edges_idx_matrix;
173  parameters[t].s = s+t*N;
174  parameters[t].f = f+t*N;
175  parameters[t].m_k = m_k;
176  parameters[t].shortest_D = shortest_D;
177  pthread_create(&threads[t], &attr, CIsomap::run_dijkstra_thread, (void*)&parameters[t]);
178  }
179  for (t=0; t<num_threads; t++)
180  pthread_join(threads[t], NULL);
181  pthread_attr_destroy(&attr);
182  for (t=0; t<num_threads; t++)
183  delete heaps[t];
184  SG_FREE(heaps);
185  SG_FREE(parameters);
186  SG_FREE(threads);
187 #else
188  D_THREAD_PARAM single_thread_param;
189  single_thread_param.i_start = 0;
190  single_thread_param.i_stop = N;
191  single_thread_param.i_step = 1;
192  single_thread_param.m_k = m_k;
193  single_thread_param.heap = new CFibonacciHeap(N);
194  single_thread_param.edges_matrix = edges_matrix;
195  single_thread_param.edges_idx_matrix = edges_idx_matrix;
196  single_thread_param.s = s;
197  single_thread_param.f = f;
198  single_thread_param.shortest_D = shortest_D;
199  run_dijkstra_thread((void*)&single_thread_param);
200  delete single_thread_param.heap;
201 #endif
202  // cleanup
203  SG_FREE(edges_matrix);
204  SG_FREE(edges_idx_matrix);
205  SG_FREE(s);
206  SG_FREE(f);
207 
208  return SGMatrix<float64_t>(shortest_D,N,N);
209 }
210 
212 {
213  // get parameters from p
214  DIJKSTRA_THREAD_PARAM* parameters = (DIJKSTRA_THREAD_PARAM*)p;
215  CFibonacciHeap* heap = parameters->heap;
216  int32_t i_start = parameters->i_start;
217  int32_t i_stop = parameters->i_stop;
218  int32_t i_step = parameters->i_step;
219  bool* s = parameters->s;
220  bool* f = parameters->f;
221  const float64_t* edges_matrix = parameters->edges_matrix;
222  const int32_t* edges_idx_matrix = parameters->edges_idx_matrix;
223  float64_t* shortest_D = parameters->shortest_D;
224  int32_t m_k = parameters->m_k;
225  int32_t k,j,i,min_item,w;
226  int32_t N = i_stop;
227 
228  // temporary vars
229  float64_t dist,tmp;
230 
231  // main loop
232  for (k=i_start; k<i_stop; k+=i_step)
233  {
234  // fill s and f with false, fill shortest_D with infinity
235  for (j=0; j<N; j++)
236  {
237  shortest_D[k*N+j] = CMath::ALMOST_INFTY;
238  s[j] = false;
239  f[j] = false;
240  }
241  // set distance from k to k as zero
242  shortest_D[k*N+k] = 0.0;
243 
244  // insert kth object to heap with zero distance and set f[k] true
245  heap->insert(k,0.0);
246  f[k] = true;
247 
248  // while heap is not empty
249  while (heap->get_num_nodes()>0)
250  {
251  // extract min and set (s)olution state as true and (f)rontier as false
252  min_item = heap->extract_min(tmp);
253  s[min_item] = true;
254  f[min_item] = false;
255 
256  // for-each edge (min_item->w)
257  for (i=0; i<m_k; i++)
258  {
259  // get w idx
260  w = edges_idx_matrix[min_item*m_k+i];
261  // if w is not in solution yet
262  if (s[w] == false)
263  {
264  // get distance from k to i through min_item
265  dist = shortest_D[k*N+min_item] + edges_matrix[min_item*m_k+i];
266  // if distance can be relaxed
267  if (dist < shortest_D[k*N+w])
268  {
269  // relax distance
270  shortest_D[k*N+w] = dist;
271  // if w is in (f)rontier
272  if (f[w])
273  {
274  // decrease distance in heap
275  heap->decrease_key(w, dist);
276  }
277  else
278  {
279  // insert w to heap and set (f)rontier as true
280  heap->insert(w, dist);
281  f[w] = true;
282  }
283  }
284  }
285  }
286  }
287  // clear heap to re-use
288  heap->clear();
289  }
290  return NULL;
291 }
292 
293 #endif /* HAVE_LAPACK */

SHOGUN Machine Learning Toolbox - Documentation