SHOGUN  v1.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
KernelLocallyLinearEmbedding.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
15 #include <shogun/lib/common.h>
16 #include <shogun/base/DynArray.h>
18 #include <shogun/io/SGIO.h>
19 #include <shogun/lib/Signal.h>
20 
21 #ifdef HAVE_PTHREAD
22 #include <pthread.h>
23 #endif
24 
25 using namespace shogun;
26 
27 #ifndef DOXYGEN_SHOULD_SKIP_THIS
28 struct LK_RECONSTRUCTION_THREAD_PARAM
29 {
31  int32_t idx_start;
33  int32_t idx_stop;
35  int32_t idx_step;
37  int32_t m_k;
39  int32_t N;
41  const int32_t* neighborhood_matrix;
43  float64_t* local_gram_matrix;
45  const float64_t* kernel_matrix;
47  float64_t* id_vector;
49  float64_t reconstruction_shift;
51  float64_t* W_matrix;
52 };
53 
54 struct K_NEIGHBORHOOD_THREAD_PARAM
55 {
57  int32_t idx_start;
59  int32_t idx_step;
61  int32_t idx_stop;
63  int32_t N;
65  int32_t m_k;
67  CFibonacciHeap* heap;
69  const float64_t* kernel_matrix;
71  int32_t* neighborhood_matrix;
72 };
73 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
74 
77 {
78 }
79 
82 {
83  set_kernel(kernel);
84 }
85 
87 {
88  return "KernelLocallyLinearEmbedding";
89 };
90 
92 {
93 }
94 
96 {
97  ASSERT(features);
98  SG_REF(features);
99 
100  // get dimensionality and number of vectors of data
101  int32_t N = features->get_num_vectors();
102  if (m_k>=N)
103  SG_ERROR("Number of neighbors (%d) should be less than number of objects (%d).\n",
104  m_k, N);
105 
106  // compute kernel matrix
107  ASSERT(m_kernel);
108  m_kernel->init(features,features);
110  m_kernel->cleanup();
111  SG_UNREF(features);
112  return (CFeatures*)embedding;
113 }
114 
116 {
117  SGMatrix<float64_t> kernel_matrix = kernel->get_kernel_matrix();
118  SGMatrix<int32_t> neighborhood_matrix = get_neighborhood_matrix(kernel_matrix,m_k);
119 
120  SGMatrix<float64_t> M_matrix = construct_weight_matrix(kernel_matrix,neighborhood_matrix);
121  neighborhood_matrix.destroy_matrix();
122 
124  M_matrix.destroy_matrix();
125 
126  return new CSimpleFeatures<float64_t>(nullspace);
127 }
128 
130  SGMatrix<int32_t> neighborhood_matrix)
131 {
132  int32_t N = kernel_matrix.num_cols;
133  // loop variables
134  int32_t t;
135 #ifdef HAVE_PTHREAD
136  int32_t num_threads = parallel->get_num_threads();
137  ASSERT(num_threads>0);
138  // allocate threads
139  pthread_t* threads = SG_MALLOC(pthread_t, num_threads);
140  LK_RECONSTRUCTION_THREAD_PARAM* parameters = SG_MALLOC(LK_RECONSTRUCTION_THREAD_PARAM, num_threads);
141  pthread_attr_t attr;
142  pthread_attr_init(&attr);
143  pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
144 #else
145  int32_t num_threads = 1;
146 #endif
147  float64_t* W_matrix = SG_CALLOC(float64_t, N*N);
148  // init matrices and norm factor to be used
149  float64_t* local_gram_matrix = SG_MALLOC(float64_t, m_k*m_k*num_threads);
150  float64_t* id_vector = SG_MALLOC(float64_t, m_k*num_threads);
151 
152 #ifdef HAVE_PTHREAD
153  for (t=0; t<num_threads; t++)
154  {
155  parameters[t].idx_start = t;
156  parameters[t].idx_step = num_threads;
157  parameters[t].idx_stop = N;
158  parameters[t].m_k = m_k;
159  parameters[t].N = N;
160  parameters[t].neighborhood_matrix = neighborhood_matrix.matrix;
161  parameters[t].kernel_matrix = kernel_matrix.matrix;
162  parameters[t].local_gram_matrix = local_gram_matrix+(m_k*m_k)*t;
163  parameters[t].id_vector = id_vector+m_k*t;
164  parameters[t].W_matrix = W_matrix;
165  parameters[t].reconstruction_shift = m_reconstruction_shift;
166  pthread_create(&threads[t], &attr, run_linearreconstruction_thread, (void*)&parameters[t]);
167  }
168  for (t=0; t<num_threads; t++)
169  pthread_join(threads[t], NULL);
170  pthread_attr_destroy(&attr);
171  SG_FREE(parameters);
172  SG_FREE(threads);
173 #else
174  LK_RECONSTRUCTION_THREAD_PARAM single_thread_param;
175  single_thread_param.idx_start = 0;
176  single_thread_param.idx_step = 1;
177  single_thread_param.idx_stop = N;
178  single_thread_param.m_k = m_k;
179  single_thread_param.N = N;
180  single_thread_param.neighborhood_matrix = neighborhood_matrix.matrix;
181  single_thread_param.local_gram_matrix = local_gram_matrix;
182  single_thread_param.kernel_matrix = kernel_matrix.matrix;
183  single_thread_param.id_vector = id_vector;
184  single_thread_param.W_matrix = W_matrix;
185  run_linearreconstruction_thread((void*)single_thread_param);
186 #endif
187 
188  // clean
189  SG_FREE(id_vector);
190  SG_FREE(local_gram_matrix);
191 
192  return SGMatrix<float64_t>(W_matrix,N,N);
193 }
194 
196 {
197  LK_RECONSTRUCTION_THREAD_PARAM* parameters = (LK_RECONSTRUCTION_THREAD_PARAM*)p;
198  int32_t idx_start = parameters->idx_start;
199  int32_t idx_step = parameters->idx_step;
200  int32_t idx_stop = parameters->idx_stop;
201  int32_t m_k = parameters->m_k;
202  int32_t N = parameters->N;
203  const int32_t* neighborhood_matrix = parameters->neighborhood_matrix;
204  float64_t* local_gram_matrix = parameters->local_gram_matrix;
205  const float64_t* kernel_matrix = parameters->kernel_matrix;
206  float64_t* id_vector = parameters->id_vector;
207  float64_t* W_matrix = parameters->W_matrix;
208  float64_t reconstruction_shift = parameters->reconstruction_shift;
209 
210  int32_t i,j,k;
211  float64_t norming,trace;
212 
213  for (i=idx_start; i<idx_stop; i+=idx_step)
214  {
215  for (j=0; j<m_k; j++)
216  {
217  for (k=0; k<m_k; k++)
218  local_gram_matrix[j*m_k+k] =
219  kernel_matrix[i*N+i] -
220  kernel_matrix[i*N+neighborhood_matrix[j*N+i]] -
221  kernel_matrix[i*N+neighborhood_matrix[k*N+i]] +
222  kernel_matrix[neighborhood_matrix[j*N+i]*N+neighborhood_matrix[k*N+i]];
223  }
224 
225  for (j=0; j<m_k; j++)
226  id_vector[j] = 1.0;
227 
228  // compute tr(C)
229  trace = 0.0;
230  for (j=0; j<m_k; j++)
231  trace += local_gram_matrix[j*m_k+j];
232 
233  // regularize gram matrix
234  for (j=0; j<m_k; j++)
235  local_gram_matrix[j*m_k+j] += reconstruction_shift*trace;
236 
237  clapack_dposv(CblasColMajor,CblasLower,m_k,1,local_gram_matrix,m_k,id_vector,m_k);
238 
239  // normalize weights
240  norming=0.0;
241  for (j=0; j<m_k; j++)
242  norming += id_vector[j];
243 
244  cblas_dscal(m_k,1.0/norming,id_vector,1);
245 
246  memset(local_gram_matrix,0,sizeof(float64_t)*m_k*m_k);
247  cblas_dger(CblasColMajor,m_k,m_k,1.0,id_vector,1,id_vector,1,local_gram_matrix,m_k);
248 
249  // put weights into W matrix
250  W_matrix[N*i+i] += 1.0;
251  for (j=0; j<m_k; j++)
252  {
253  W_matrix[N*i+neighborhood_matrix[j*N+i]] -= id_vector[j];
254  W_matrix[N*neighborhood_matrix[j*N+i]+i] -= id_vector[j];
255  }
256  for (j=0; j<m_k; j++)
257  {
258  for (k=0; k<m_k; k++)
259  W_matrix[N*neighborhood_matrix[j*N+i]+neighborhood_matrix[k*N+i]]+=local_gram_matrix[j*m_k+k];
260  }
261  }
262  return NULL;
263 }
264 
266 {
267  int32_t t;
268  int32_t N = kernel_matrix.num_cols;
269  // init matrix and heap to be used
270  int32_t* neighborhood_matrix = SG_MALLOC(int32_t, N*k);
271 #ifdef HAVE_PTHREAD
272  int32_t num_threads = parallel->get_num_threads();
273  ASSERT(num_threads>0);
274  K_NEIGHBORHOOD_THREAD_PARAM* parameters = SG_MALLOC(K_NEIGHBORHOOD_THREAD_PARAM, num_threads);
275  pthread_t* threads = SG_MALLOC(pthread_t, num_threads);
276  pthread_attr_t attr;
277  pthread_attr_init(&attr);
278  pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
279 #else
280  int32_t num_threads = 1;
281 #endif
282  CFibonacciHeap** heaps = SG_MALLOC(CFibonacciHeap*, num_threads);
283  for (t=0; t<num_threads; t++)
284  heaps[t] = new CFibonacciHeap(N);
285 
286 #ifdef HAVE_PTHREAD
287  for (t=0; t<num_threads; t++)
288  {
289  parameters[t].idx_start = t;
290  parameters[t].idx_step = num_threads;
291  parameters[t].idx_stop = N;
292  parameters[t].m_k = k;
293  parameters[t].N = N;
294  parameters[t].heap = heaps[t];
295  parameters[t].neighborhood_matrix = neighborhood_matrix;
296  parameters[t].kernel_matrix = kernel_matrix.matrix;
297  pthread_create(&threads[t], &attr, run_neighborhood_thread, (void*)&parameters[t]);
298  }
299  for (t=0; t<num_threads; t++)
300  pthread_join(threads[t], NULL);
301  pthread_attr_destroy(&attr);
302  SG_FREE(threads);
303  SG_FREE(parameters);
304 #else
305  K_NEIGHBORHOOD_THREAD_PARAM single_thread_param;
306  single_thread_param.idx_start = 0;
307  single_thread_param.idx_step = 1;
308  single_thread_param.idx_stop = N;
309  single_thread_param.m_k = k;
310  single_thread_param.N = N;
311  single_thread_param.heap = heaps[0]
312  single_thread_param.neighborhood_matrix = neighborhood_matrix;
313  single_thread_param.kernel_matrix = kernel_matrix.matrix;
314  run_neighborhood_thread((void*)&single_thread_param);
315 #endif
316 
317  for (t=0; t<num_threads; t++)
318  delete heaps[t];
319  SG_FREE(heaps);
320 
321  return SGMatrix<int32_t>(neighborhood_matrix,k,N);
322 }
323 
325 {
326  K_NEIGHBORHOOD_THREAD_PARAM* parameters = (K_NEIGHBORHOOD_THREAD_PARAM*)p;
327  int32_t idx_start = parameters->idx_start;
328  int32_t idx_step = parameters->idx_step;
329  int32_t idx_stop = parameters->idx_stop;
330  int32_t N = parameters->N;
331  int32_t m_k = parameters->m_k;
332  CFibonacciHeap* heap = parameters->heap;
333  const float64_t* kernel_matrix = parameters->kernel_matrix;
334  int32_t* neighborhood_matrix = parameters->neighborhood_matrix;
335 
336  int32_t i,j;
337  float64_t tmp;
338  for (i=idx_start; i<idx_stop; i+=idx_step)
339  {
340  for (j=0; j<N; j++)
341  {
342  heap->insert(j,kernel_matrix[i*N+i]-2*kernel_matrix[i*N+j]+kernel_matrix[j*N+j]);
343  }
344 
345  heap->extract_min(tmp);
346 
347  for (j=0; j<m_k; j++)
348  neighborhood_matrix[j*N+i] = heap->extract_min(tmp);
349 
350  heap->clear();
351  }
352 
353  return NULL;
354 }
355 #endif /* HAVE_LAPACK */

SHOGUN Machine Learning Toolbox - Documentation