SHOGUN  v1.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
LocalTangentSpaceAlignment.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
14 #include <shogun/lib/Time.h>
15 #include <shogun/lib/common.h>
17 #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 LTSA_THREAD_PARAM
29 {
31  int32_t idx_start;
33  int32_t idx_step;
35  int32_t idx_stop;
37  int32_t m_k;
39  int32_t target_dim;
41  int32_t dim;
43  int32_t N;
45  const int32_t* neighborhood_matrix;
47  float64_t* G_matrix;
49  float64_t* mean_vector;
51  float64_t* local_feature_matrix;
53  const float64_t* feature_matrix;
55  float64_t* s_values_vector;
57  float64_t* q_matrix;
59  float64_t* W_matrix;
60 #ifdef HAVE_PTHREAD
61 
62  PTHREAD_LOCK_T* W_matrix_lock;
63 #endif
64 };
65 #endif
66 
69 {
70 }
71 
73 {
74 }
75 
77 {
78  return "LocalTangentSpaceAlignment";
79 };
80 
82  SGMatrix<int32_t> neighborhood_matrix)
83 {
84  int32_t N = simple_features->get_num_vectors();
85  int32_t dim = simple_features->get_num_features();
86  int32_t t;
87 #ifdef HAVE_PTHREAD
88  int32_t num_threads = parallel->get_num_threads();
89  ASSERT(num_threads>0);
90  // allocate threads and params
91  pthread_t* threads = SG_MALLOC(pthread_t, num_threads);
92  LTSA_THREAD_PARAM* parameters = SG_MALLOC(LTSA_THREAD_PARAM, num_threads);
93 #else
94  int32_t num_threads = 1;
95 #endif
96 
97  // init matrices and norm factor to be used
98  float64_t* local_feature_matrix = SG_MALLOC(float64_t, m_k*dim*num_threads);
99  float64_t* mean_vector = SG_MALLOC(float64_t, dim*num_threads);
100  float64_t* q_matrix = SG_MALLOC(float64_t, m_k*m_k*num_threads);
101  float64_t* s_values_vector = SG_MALLOC(float64_t, dim*num_threads);
102  float64_t* G_matrix = SG_MALLOC(float64_t, m_k*(1+m_target_dim)*num_threads);
103 
104  // get feature matrix
105  SGMatrix<float64_t> feature_matrix = simple_features->get_feature_matrix();
106 
107 #ifdef HAVE_PTHREAD
108  PTHREAD_LOCK_T W_matrix_lock;
109  pthread_attr_t attr;
110  PTHREAD_LOCK_INIT(&W_matrix_lock);
111  pthread_attr_init(&attr);
112  pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
113 
114  for (t=0; t<num_threads; t++)
115  {
116  parameters[t].idx_start = t;
117  parameters[t].idx_step = num_threads;
118  parameters[t].idx_stop = N;
119  parameters[t].m_k = m_k;
120  parameters[t].target_dim = m_target_dim;
121  parameters[t].dim = dim;
122  parameters[t].N = N;
123  parameters[t].neighborhood_matrix = neighborhood_matrix.matrix;
124  parameters[t].G_matrix = G_matrix + (m_k*(1+m_target_dim))*t;
125  parameters[t].mean_vector = mean_vector + dim*t;
126  parameters[t].local_feature_matrix = local_feature_matrix + (m_k*dim)*t;
127  parameters[t].feature_matrix = feature_matrix.matrix;
128  parameters[t].s_values_vector = s_values_vector + dim*t;
129  parameters[t].q_matrix = q_matrix + (m_k*m_k)*t;
130  parameters[t].W_matrix = W_matrix;
131  parameters[t].W_matrix_lock = &W_matrix_lock;
132  pthread_create(&threads[t], &attr, run_ltsa_thread, (void*)&parameters[t]);
133  }
134  for (t=0; t<num_threads; t++)
135  pthread_join(threads[t], NULL);
136  PTHREAD_LOCK_DESTROY(&W_matrix_lock);
137  SG_FREE(parameters);
138  SG_FREE(threads);
139 #else
140  LTSA_THREAD_PARAM single_thread_param;
141  single_thread_param.idx_start = 0;
142  single_thread_param.idx_step = 1;
143  single_thread_param.idx_stop = N;
144  single_thread_param.m_k = m_k;
145  single_thread_param.target_dim = m_target_dim;
146  single_thread_param.dim = dim;
147  single_thread_param.N = N;
148  single_thread_param.neighborhood_matrix = neighborhood_matrix.matrix;
149  single_thread_param.G_matrix = G_matrix;
150  single_thread_param.mean_vector = mean_vector;
151  single_thread_param.local_feature_matrix = local_feature_matrix;
152  single_thread_param.feature_matrix = feature_matrix.matrix;
153  single_thread_param.s_values_vector = s_values_vector;
154  single_thread_param.q_matrix = q_matrix;
155  single_thread_param.W_matrix = W_matrix;
156  run_ltsa_thread((void*)&single_thread_param);
157 #endif
158 
159  // clean
160  SG_FREE(G_matrix);
161  SG_FREE(s_values_vector);
162  SG_FREE(mean_vector);
163  SG_FREE(local_feature_matrix);
164  SG_FREE(q_matrix);
165 
166  for (int32_t i=0; i<N; i++)
167  {
168  for (int32_t j=0; j<m_k; j++)
169  W_matrix[N*neighborhood_matrix[j*N+i]+neighborhood_matrix[j*N+i]] += 1.0;
170  }
171 
172  return SGMatrix<float64_t>(W_matrix,N,N);
173 }
174 
176 {
177  LTSA_THREAD_PARAM* parameters = (LTSA_THREAD_PARAM*)p;
178  int32_t idx_start = parameters->idx_start;
179  int32_t idx_step = parameters->idx_step;
180  int32_t idx_stop = parameters->idx_stop;
181  int32_t m_k = parameters->m_k;
182  int32_t target_dim = parameters->target_dim;
183  int32_t dim = parameters->dim;
184  int32_t N = parameters->N;
185  const int32_t* neighborhood_matrix = parameters->neighborhood_matrix;
186  float64_t* G_matrix = parameters->G_matrix;
187  float64_t* mean_vector = parameters->mean_vector;
188  float64_t* local_feature_matrix = parameters->local_feature_matrix;
189  const float64_t* feature_matrix = parameters->feature_matrix;
190  float64_t* s_values_vector = parameters->s_values_vector;
191  float64_t* q_matrix = parameters->q_matrix;
192  float64_t* W_matrix = parameters->W_matrix;
193 #ifdef HAVE_PTHREAD
194  PTHREAD_LOCK_T* W_matrix_lock = parameters->W_matrix_lock;
195 #endif
196 
197  int32_t i,j,k;
198 
199  for (j=0; j<m_k; j++)
200  G_matrix[j] = 1.0/CMath::sqrt((float64_t)m_k);
201 
202  for (i=idx_start; i<idx_stop; i+=idx_step)
203  {
204  // fill mean vector with zeros
205  memset(mean_vector,0,sizeof(float64_t)*dim);
206 
207  // compute local feature matrix containing neighbors of i-th vector
208  for (j=0; j<m_k; j++)
209  {
210  for (k=0; k<dim; k++)
211  local_feature_matrix[j*dim+k] = feature_matrix[neighborhood_matrix[j*N+i]*dim+k];
212 
213  cblas_daxpy(dim,1.0,local_feature_matrix+j*dim,1,mean_vector,1);
214  }
215 
216  // compute mean
217  cblas_dscal(dim,1.0/m_k,mean_vector,1);
218 
219  // center feature vectors by mean
220  for (j=0; j<m_k; j++)
221  cblas_daxpy(dim,-1.0,mean_vector,1,local_feature_matrix+j*dim,1);
222 
223  int32_t info = 0;
224  // find right eigenvectors of local_feature_matrix
225  wrap_dgesvd('N','O',dim,m_k,local_feature_matrix,dim,
226  s_values_vector,NULL,1, NULL,1,&info);
227  ASSERT(info==0);
228 
229  for (j=0; j<target_dim; j++)
230  {
231  for (k=0; k<m_k; k++)
232  G_matrix[(j+1)*m_k+k] = local_feature_matrix[k*dim+j];
233  }
234 
235  // compute GG'
236  cblas_dgemm(CblasColMajor,CblasNoTrans,CblasTrans,
237  m_k,m_k,1+target_dim,
238  1.0,G_matrix,m_k,
239  G_matrix,m_k,
240  0.0,q_matrix,m_k);
241 
242  // W[neighbors of i, neighbors of i] = I - GG'
243 #ifdef HAVE_PTHREAD
244  PTHREAD_LOCK(W_matrix_lock);
245 #endif
246  for (j=0; j<m_k; j++)
247  {
248  for (k=0; k<m_k; k++)
249  W_matrix[N*neighborhood_matrix[k*N+i]+neighborhood_matrix[j*N+i]] -= q_matrix[j*m_k+k];
250  }
251 #ifdef HAVE_PTHREAD
252  PTHREAD_UNLOCK(W_matrix_lock);
253 #endif
254  }
255  return NULL;
256 }
257 
258 #endif /* HAVE_LAPACK */

SHOGUN Machine Learning Toolbox - Documentation