SHOGUN
v1.1.0
|
00001 /* 00002 * This program is free software; you can redistribute it and/or modify 00003 * it under the terms of the GNU General Public License as published by 00004 * the Free Software Foundation; either version 3 of the License, or 00005 * (at your option) any later version. 00006 * 00007 * Written (W) 2011 Sergey Lisitsyn 00008 * Copyright (C) 2011 Berlin Institute of Technology and Max-Planck-Society 00009 */ 00010 00011 #include <shogun/converter/LocalTangentSpaceAlignment.h> 00012 #ifdef HAVE_LAPACK 00013 #include <shogun/mathematics/lapack.h> 00014 #include <shogun/lib/Time.h> 00015 #include <shogun/lib/common.h> 00016 #include <shogun/mathematics/Math.h> 00017 #include <shogun/io/SGIO.h> 00018 #include <shogun/distance/Distance.h> 00019 #include <shogun/lib/Signal.h> 00020 00021 #ifdef HAVE_PTHREAD 00022 #include <pthread.h> 00023 #endif 00024 00025 using namespace shogun; 00026 00027 #ifndef DOXYGEN_SHOULD_SKIP_THIS 00028 struct LTSA_THREAD_PARAM 00029 { 00031 int32_t idx_start; 00033 int32_t idx_step; 00035 int32_t idx_stop; 00037 int32_t m_k; 00039 int32_t target_dim; 00041 int32_t dim; 00043 int32_t N; 00045 const int32_t* neighborhood_matrix; 00047 float64_t* G_matrix; 00049 float64_t* mean_vector; 00051 float64_t* local_feature_matrix; 00053 const float64_t* feature_matrix; 00055 float64_t* s_values_vector; 00057 float64_t* q_matrix; 00059 float64_t* W_matrix; 00060 #ifdef HAVE_PTHREAD 00061 00062 PTHREAD_LOCK_T* W_matrix_lock; 00063 #endif 00064 }; 00065 #endif 00066 00067 CLocalTangentSpaceAlignment::CLocalTangentSpaceAlignment() : 00068 CLocallyLinearEmbedding() 00069 { 00070 } 00071 00072 CLocalTangentSpaceAlignment::~CLocalTangentSpaceAlignment() 00073 { 00074 } 00075 00076 const char* CLocalTangentSpaceAlignment::get_name() const 00077 { 00078 return "LocalTangentSpaceAlignment"; 00079 }; 00080 00081 SGMatrix<float64_t> CLocalTangentSpaceAlignment::construct_weight_matrix(CSimpleFeatures<float64_t>* simple_features, float64_t* W_matrix, 00082 SGMatrix<int32_t> neighborhood_matrix) 00083 { 00084 int32_t N = simple_features->get_num_vectors(); 00085 int32_t dim = simple_features->get_num_features(); 00086 int32_t t; 00087 #ifdef HAVE_PTHREAD 00088 int32_t num_threads = parallel->get_num_threads(); 00089 ASSERT(num_threads>0); 00090 // allocate threads and params 00091 pthread_t* threads = SG_MALLOC(pthread_t, num_threads); 00092 LTSA_THREAD_PARAM* parameters = SG_MALLOC(LTSA_THREAD_PARAM, num_threads); 00093 #else 00094 int32_t num_threads = 1; 00095 #endif 00096 00097 // init matrices and norm factor to be used 00098 float64_t* local_feature_matrix = SG_MALLOC(float64_t, m_k*dim*num_threads); 00099 float64_t* mean_vector = SG_MALLOC(float64_t, dim*num_threads); 00100 float64_t* q_matrix = SG_MALLOC(float64_t, m_k*m_k*num_threads); 00101 float64_t* s_values_vector = SG_MALLOC(float64_t, dim*num_threads); 00102 float64_t* G_matrix = SG_MALLOC(float64_t, m_k*(1+m_target_dim)*num_threads); 00103 00104 // get feature matrix 00105 SGMatrix<float64_t> feature_matrix = simple_features->get_feature_matrix(); 00106 00107 #ifdef HAVE_PTHREAD 00108 PTHREAD_LOCK_T W_matrix_lock; 00109 pthread_attr_t attr; 00110 PTHREAD_LOCK_INIT(&W_matrix_lock); 00111 pthread_attr_init(&attr); 00112 pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE); 00113 00114 for (t=0; t<num_threads; t++) 00115 { 00116 parameters[t].idx_start = t; 00117 parameters[t].idx_step = num_threads; 00118 parameters[t].idx_stop = N; 00119 parameters[t].m_k = m_k; 00120 parameters[t].target_dim = m_target_dim; 00121 parameters[t].dim = dim; 00122 parameters[t].N = N; 00123 parameters[t].neighborhood_matrix = neighborhood_matrix.matrix; 00124 parameters[t].G_matrix = G_matrix + (m_k*(1+m_target_dim))*t; 00125 parameters[t].mean_vector = mean_vector + dim*t; 00126 parameters[t].local_feature_matrix = local_feature_matrix + (m_k*dim)*t; 00127 parameters[t].feature_matrix = feature_matrix.matrix; 00128 parameters[t].s_values_vector = s_values_vector + dim*t; 00129 parameters[t].q_matrix = q_matrix + (m_k*m_k)*t; 00130 parameters[t].W_matrix = W_matrix; 00131 parameters[t].W_matrix_lock = &W_matrix_lock; 00132 pthread_create(&threads[t], &attr, run_ltsa_thread, (void*)¶meters[t]); 00133 } 00134 for (t=0; t<num_threads; t++) 00135 pthread_join(threads[t], NULL); 00136 PTHREAD_LOCK_DESTROY(&W_matrix_lock); 00137 SG_FREE(parameters); 00138 SG_FREE(threads); 00139 #else 00140 LTSA_THREAD_PARAM single_thread_param; 00141 single_thread_param.idx_start = 0; 00142 single_thread_param.idx_step = 1; 00143 single_thread_param.idx_stop = N; 00144 single_thread_param.m_k = m_k; 00145 single_thread_param.target_dim = m_target_dim; 00146 single_thread_param.dim = dim; 00147 single_thread_param.N = N; 00148 single_thread_param.neighborhood_matrix = neighborhood_matrix.matrix; 00149 single_thread_param.G_matrix = G_matrix; 00150 single_thread_param.mean_vector = mean_vector; 00151 single_thread_param.local_feature_matrix = local_feature_matrix; 00152 single_thread_param.feature_matrix = feature_matrix.matrix; 00153 single_thread_param.s_values_vector = s_values_vector; 00154 single_thread_param.q_matrix = q_matrix; 00155 single_thread_param.W_matrix = W_matrix; 00156 run_ltsa_thread((void*)&single_thread_param); 00157 #endif 00158 00159 // clean 00160 SG_FREE(G_matrix); 00161 SG_FREE(s_values_vector); 00162 SG_FREE(mean_vector); 00163 SG_FREE(local_feature_matrix); 00164 SG_FREE(q_matrix); 00165 00166 for (int32_t i=0; i<N; i++) 00167 { 00168 for (int32_t j=0; j<m_k; j++) 00169 W_matrix[N*neighborhood_matrix[j*N+i]+neighborhood_matrix[j*N+i]] += 1.0; 00170 } 00171 00172 return SGMatrix<float64_t>(W_matrix,N,N); 00173 } 00174 00175 void* CLocalTangentSpaceAlignment::run_ltsa_thread(void* p) 00176 { 00177 LTSA_THREAD_PARAM* parameters = (LTSA_THREAD_PARAM*)p; 00178 int32_t idx_start = parameters->idx_start; 00179 int32_t idx_step = parameters->idx_step; 00180 int32_t idx_stop = parameters->idx_stop; 00181 int32_t m_k = parameters->m_k; 00182 int32_t target_dim = parameters->target_dim; 00183 int32_t dim = parameters->dim; 00184 int32_t N = parameters->N; 00185 const int32_t* neighborhood_matrix = parameters->neighborhood_matrix; 00186 float64_t* G_matrix = parameters->G_matrix; 00187 float64_t* mean_vector = parameters->mean_vector; 00188 float64_t* local_feature_matrix = parameters->local_feature_matrix; 00189 const float64_t* feature_matrix = parameters->feature_matrix; 00190 float64_t* s_values_vector = parameters->s_values_vector; 00191 float64_t* q_matrix = parameters->q_matrix; 00192 float64_t* W_matrix = parameters->W_matrix; 00193 #ifdef HAVE_PTHREAD 00194 PTHREAD_LOCK_T* W_matrix_lock = parameters->W_matrix_lock; 00195 #endif 00196 00197 int32_t i,j,k; 00198 00199 for (j=0; j<m_k; j++) 00200 G_matrix[j] = 1.0/CMath::sqrt((float64_t)m_k); 00201 00202 for (i=idx_start; i<idx_stop; i+=idx_step) 00203 { 00204 // fill mean vector with zeros 00205 memset(mean_vector,0,sizeof(float64_t)*dim); 00206 00207 // compute local feature matrix containing neighbors of i-th vector 00208 for (j=0; j<m_k; j++) 00209 { 00210 for (k=0; k<dim; k++) 00211 local_feature_matrix[j*dim+k] = feature_matrix[neighborhood_matrix[j*N+i]*dim+k]; 00212 00213 cblas_daxpy(dim,1.0,local_feature_matrix+j*dim,1,mean_vector,1); 00214 } 00215 00216 // compute mean 00217 cblas_dscal(dim,1.0/m_k,mean_vector,1); 00218 00219 // center feature vectors by mean 00220 for (j=0; j<m_k; j++) 00221 cblas_daxpy(dim,-1.0,mean_vector,1,local_feature_matrix+j*dim,1); 00222 00223 int32_t info = 0; 00224 // find right eigenvectors of local_feature_matrix 00225 wrap_dgesvd('N','O',dim,m_k,local_feature_matrix,dim, 00226 s_values_vector,NULL,1, NULL,1,&info); 00227 ASSERT(info==0); 00228 00229 for (j=0; j<target_dim; j++) 00230 { 00231 for (k=0; k<m_k; k++) 00232 G_matrix[(j+1)*m_k+k] = local_feature_matrix[k*dim+j]; 00233 } 00234 00235 // compute GG' 00236 cblas_dgemm(CblasColMajor,CblasNoTrans,CblasTrans, 00237 m_k,m_k,1+target_dim, 00238 1.0,G_matrix,m_k, 00239 G_matrix,m_k, 00240 0.0,q_matrix,m_k); 00241 00242 // W[neighbors of i, neighbors of i] = I - GG' 00243 #ifdef HAVE_PTHREAD 00244 PTHREAD_LOCK(W_matrix_lock); 00245 #endif 00246 for (j=0; j<m_k; j++) 00247 { 00248 for (k=0; k<m_k; k++) 00249 W_matrix[N*neighborhood_matrix[k*N+i]+neighborhood_matrix[j*N+i]] -= q_matrix[j*m_k+k]; 00250 } 00251 #ifdef HAVE_PTHREAD 00252 PTHREAD_UNLOCK(W_matrix_lock); 00253 #endif 00254 } 00255 return NULL; 00256 } 00257 00258 #endif /* HAVE_LAPACK */