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