00001
00002
00003
00004
00005
00006
00007
00008
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
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
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
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
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
00211 memset(mean_vector,0,sizeof(float64_t)*dim);
00212
00213
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
00223 cblas_dscal(dim,1.0/m_k,mean_vector,1);
00224
00225
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
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
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
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