SHOGUN
v2.0.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 const int32_t* neighborhood_matrix; 00043 float64_t* G_matrix; 00045 float64_t* mean_vector; 00047 float64_t* local_feature_matrix; 00049 const float64_t* feature_matrix; 00051 float64_t* s_values_vector; 00053 float64_t* q_matrix; 00055 float64_t* W_matrix; 00057 int32_t N; 00059 int32_t dim; 00061 int32_t actual_k; 00062 #ifdef HAVE_PTHREAD 00063 00064 PTHREAD_LOCK_T* W_matrix_lock; 00065 #endif 00066 }; 00067 #endif 00068 00069 CLocalTangentSpaceAlignment::CLocalTangentSpaceAlignment() : 00070 CLocallyLinearEmbedding() 00071 { 00072 } 00073 00074 CLocalTangentSpaceAlignment::~CLocalTangentSpaceAlignment() 00075 { 00076 } 00077 00078 const char* CLocalTangentSpaceAlignment::get_name() const 00079 { 00080 return "LocalTangentSpaceAlignment"; 00081 }; 00082 00083 SGMatrix<float64_t> CLocalTangentSpaceAlignment::construct_weight_matrix(CDenseFeatures<float64_t>* simple_features, float64_t* W_matrix, 00084 SGMatrix<int32_t> neighborhood_matrix) 00085 { 00086 int32_t N = simple_features->get_num_vectors(); 00087 int32_t dim = simple_features->get_num_features(); 00088 int32_t t; 00089 #ifdef HAVE_PTHREAD 00090 int32_t num_threads = parallel->get_num_threads(); 00091 ASSERT(num_threads>0); 00092 // allocate threads and params 00093 pthread_t* threads = SG_MALLOC(pthread_t, num_threads); 00094 LTSA_THREAD_PARAM* parameters = SG_MALLOC(LTSA_THREAD_PARAM, num_threads); 00095 #else 00096 int32_t num_threads = 1; 00097 #endif 00098 00099 // init matrices and norm factor to be used 00100 float64_t* local_feature_matrix = SG_MALLOC(float64_t, m_k*dim*num_threads); 00101 float64_t* mean_vector = SG_MALLOC(float64_t, dim*num_threads); 00102 float64_t* q_matrix = SG_MALLOC(float64_t, m_k*m_k*num_threads); 00103 float64_t* s_values_vector = SG_MALLOC(float64_t, dim*num_threads); 00104 float64_t* G_matrix = SG_MALLOC(float64_t, m_k*(1+m_target_dim)*num_threads); 00105 00106 // get feature matrix 00107 SGMatrix<float64_t> feature_matrix = simple_features->get_feature_matrix(); 00108 00109 #ifdef HAVE_PTHREAD 00110 PTHREAD_LOCK_T W_matrix_lock; 00111 pthread_attr_t attr; 00112 PTHREAD_LOCK_INIT(&W_matrix_lock); 00113 pthread_attr_init(&attr); 00114 pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE); 00115 00116 for (t=0; t<num_threads; t++) 00117 { 00118 parameters[t].idx_start = t; 00119 parameters[t].idx_step = num_threads; 00120 parameters[t].idx_stop = N; 00121 parameters[t].m_k = m_k; 00122 parameters[t].target_dim = m_target_dim; 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 parameters[t].N = N; 00133 parameters[t].dim = dim; 00134 parameters[t].actual_k = neighborhood_matrix.num_rows; 00135 pthread_create(&threads[t], &attr, run_ltsa_thread, (void*)¶meters[t]); 00136 } 00137 for (t=0; t<num_threads; t++) 00138 pthread_join(threads[t], NULL); 00139 PTHREAD_LOCK_DESTROY(&W_matrix_lock); 00140 SG_FREE(parameters); 00141 SG_FREE(threads); 00142 #else 00143 LTSA_THREAD_PARAM single_thread_param; 00144 single_thread_param.idx_start = 0; 00145 single_thread_param.idx_step = 1; 00146 single_thread_param.idx_stop = N; 00147 single_thread_param.m_k = m_k; 00148 single_thread_param.target_dim = m_target_dim; 00149 single_thread_param.neighborhood_matrix = neighborhood_matrix.matrix; 00150 single_thread_param.G_matrix = G_matrix; 00151 single_thread_param.mean_vector = mean_vector; 00152 single_thread_param.local_feature_matrix = local_feature_matrix; 00153 single_thread_param.feature_matrix = feature_matrix; 00154 single_thread_param.s_values_vector = s_values_vector; 00155 single_thread_param.q_matrix = q_matrix; 00156 single_thread_param.W_matrix = W_matrix; 00157 single_thread_param.N = N; 00158 single_thread_param.dim = dim; 00159 single_thread_param.actual_k = neighborhood_matrix.num_rows; 00160 run_ltsa_thread((void*)&single_thread_param); 00161 #endif 00162 00163 // clean 00164 SG_FREE(G_matrix); 00165 SG_FREE(s_values_vector); 00166 SG_FREE(mean_vector); 00167 SG_FREE(local_feature_matrix); 00168 SG_FREE(q_matrix); 00169 00170 int32_t actual_k = neighborhood_matrix.num_rows; 00171 for (int32_t i=0; i<N; i++) 00172 { 00173 for (int32_t j=0; j<m_k; j++) 00174 W_matrix[N*neighborhood_matrix[i*actual_k+j]+neighborhood_matrix[i*actual_k+j]] += 1.0; 00175 } 00176 00177 return SGMatrix<float64_t>(W_matrix,N,N); 00178 } 00179 00180 void* CLocalTangentSpaceAlignment::run_ltsa_thread(void* p) 00181 { 00182 LTSA_THREAD_PARAM* parameters = (LTSA_THREAD_PARAM*)p; 00183 int32_t idx_start = parameters->idx_start; 00184 int32_t idx_step = parameters->idx_step; 00185 int32_t idx_stop = parameters->idx_stop; 00186 int32_t m_k = parameters->m_k; 00187 int32_t target_dim = parameters->target_dim; 00188 const int32_t* neighborhood_matrix = parameters->neighborhood_matrix; 00189 float64_t* G_matrix = parameters->G_matrix; 00190 float64_t* mean_vector = parameters->mean_vector; 00191 float64_t* local_feature_matrix = parameters->local_feature_matrix; 00192 const float64_t* feature_matrix = parameters->feature_matrix; 00193 float64_t* s_values_vector = parameters->s_values_vector; 00194 float64_t* q_matrix = parameters->q_matrix; 00195 float64_t* W_matrix = parameters->W_matrix; 00196 #ifdef HAVE_PTHREAD 00197 PTHREAD_LOCK_T* W_matrix_lock = parameters->W_matrix_lock; 00198 #endif 00199 00200 int32_t i,j,k; 00201 int32_t N = parameters->N; 00202 int32_t dim = parameters->dim; 00203 int32_t actual_k = parameters->actual_k; 00204 00205 for (j=0; j<m_k; j++) 00206 G_matrix[j] = 1.0/CMath::sqrt((float64_t)m_k); 00207 00208 for (i=idx_start; i<idx_stop; i+=idx_step) 00209 { 00210 // fill mean vector with zeros 00211 memset(mean_vector,0,sizeof(float64_t)*dim); 00212 00213 // compute local feature matrix containing neighbors of i-th vector 00214 for (j=0; j<m_k; j++) 00215 { 00216 for (k=0; k<dim; k++) 00217 local_feature_matrix[j*dim+k] = feature_matrix[neighborhood_matrix[i*actual_k+j]*dim+k]; 00218 00219 cblas_daxpy(dim,1.0,local_feature_matrix+j*dim,1,mean_vector,1); 00220 } 00221 00222 // compute mean 00223 cblas_dscal(dim,1.0/m_k,mean_vector,1); 00224 00225 // center feature vectors by mean 00226 for (j=0; j<m_k; j++) 00227 cblas_daxpy(dim,-1.0,mean_vector,1,local_feature_matrix+j*dim,1); 00228 00229 int32_t info = 0; 00230 // find right eigenvectors of local_feature_matrix 00231 wrap_dgesvd('N','O',dim,m_k,local_feature_matrix,dim, 00232 s_values_vector,NULL,1, NULL,1,&info); 00233 ASSERT(info==0); 00234 00235 for (j=0; j<target_dim; j++) 00236 { 00237 for (k=0; k<m_k; k++) 00238 G_matrix[(j+1)*m_k+k] = local_feature_matrix[k*dim+j]; 00239 } 00240 00241 // compute GG' 00242 cblas_dgemm(CblasColMajor,CblasNoTrans,CblasTrans, 00243 m_k,m_k,1+target_dim, 00244 1.0,G_matrix,m_k, 00245 G_matrix,m_k, 00246 0.0,q_matrix,m_k); 00247 00248 // W[neighbors of i, neighbors of i] = I - GG' 00249 #ifdef HAVE_PTHREAD 00250 PTHREAD_LOCK(W_matrix_lock); 00251 #endif 00252 for (j=0; j<m_k; j++) 00253 { 00254 for (k=0; k<m_k; k++) 00255 W_matrix[N*neighborhood_matrix[i*actual_k+k]+neighborhood_matrix[i*actual_k+j]] -= q_matrix[j*m_k+k]; 00256 } 00257 #ifdef HAVE_PTHREAD 00258 PTHREAD_UNLOCK(W_matrix_lock); 00259 #endif 00260 } 00261 return NULL; 00262 } 00263 00264 #endif /* HAVE_LAPACK */