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/KernelLocalTangentSpaceAlignment.h> 00012 #ifdef HAVE_LAPACK 00013 #include <shogun/mathematics/lapack.h> 00014 #include <shogun/mathematics/arpack.h> 00015 #include <shogun/lib/Time.h> 00016 #include <shogun/lib/common.h> 00017 #include <shogun/mathematics/Math.h> 00018 #include <shogun/io/SGIO.h> 00019 #include <shogun/distance/Distance.h> 00020 #include <shogun/lib/Signal.h> 00021 #include <shogun/base/Parallel.h> 00022 00023 #ifdef HAVE_PTHREAD 00024 #include <pthread.h> 00025 #endif 00026 00027 using namespace shogun; 00028 00029 #ifndef DOXYGEN_SHOULD_SKIP_THIS 00030 struct KLTSA_THREAD_PARAM 00031 { 00033 int32_t idx_start; 00035 int32_t idx_step; 00037 int32_t idx_stop; 00039 int32_t m_k; 00041 int32_t target_dim; 00043 int32_t N; 00045 const int32_t* neighborhood_matrix; 00047 const float64_t* kernel_matrix; 00049 float64_t* local_gram_matrix; 00051 float64_t* ev_vector; 00053 float64_t* G_matrix; 00055 float64_t* W_matrix; 00056 #ifdef HAVE_PTHREAD 00057 00058 PTHREAD_LOCK_T* W_matrix_lock; 00059 #endif 00060 }; 00061 #endif 00062 00063 CKernelLocalTangentSpaceAlignment::CKernelLocalTangentSpaceAlignment() : 00064 CKernelLocallyLinearEmbedding() 00065 { 00066 } 00067 00068 CKernelLocalTangentSpaceAlignment::CKernelLocalTangentSpaceAlignment(CKernel* kernel) : 00069 CKernelLocallyLinearEmbedding(kernel) 00070 { 00071 } 00072 00073 CKernelLocalTangentSpaceAlignment::~CKernelLocalTangentSpaceAlignment() 00074 { 00075 } 00076 00077 const char* CKernelLocalTangentSpaceAlignment::get_name() const 00078 { 00079 return "KernelLocalTangentSpaceAlignment"; 00080 }; 00081 00082 SGMatrix<float64_t> CKernelLocalTangentSpaceAlignment::construct_weight_matrix(SGMatrix<float64_t> kernel_matrix, 00083 SGMatrix<int32_t> neighborhood_matrix) 00084 { 00085 int32_t N = kernel_matrix.num_cols; 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 KLTSA_THREAD_PARAM* parameters = SG_MALLOC(KLTSA_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_gram_matrix = SG_MALLOC(float64_t, m_k*m_k*num_threads); 00099 float64_t* G_matrix = SG_MALLOC(float64_t, m_k*(1+m_target_dim)*num_threads); 00100 float64_t* W_matrix = SG_CALLOC(float64_t, N*N); 00101 float64_t* ev_vector = SG_MALLOC(float64_t, m_k*num_threads); 00102 00103 #ifdef HAVE_PTHREAD 00104 PTHREAD_LOCK_T W_matrix_lock; 00105 pthread_attr_t attr; 00106 PTHREAD_LOCK_INIT(&W_matrix_lock); 00107 pthread_attr_init(&attr); 00108 pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE); 00109 00110 for (t=0; t<num_threads; t++) 00111 { 00112 KLTSA_THREAD_PARAM params = {t,num_threads,N,m_k,m_target_dim,N,neighborhood_matrix.matrix, 00113 kernel_matrix.matrix,local_gram_matrix+(m_k*m_k)*t,ev_vector+m_k*t, 00114 G_matrix+(m_k*(1+m_target_dim))*t,W_matrix,&W_matrix_lock}; 00115 parameters[t] = params; 00116 pthread_create(&threads[t], &attr, run_kltsa_thread, (void*)¶meters[t]); 00117 } 00118 for (t=0; t<num_threads; t++) 00119 pthread_join(threads[t], NULL); 00120 PTHREAD_LOCK_DESTROY(&W_matrix_lock); 00121 SG_FREE(parameters); 00122 SG_FREE(threads); 00123 #else 00124 KLTSA_THREAD_PARAM single_thread_param = {0,1,N,m_k,m_target_dim,neighborhood_matrix.matrix, 00125 kernel_matrix.matrix,local_gram_matrix,ev_vector, 00126 G_matrix,W_matrix}; 00127 run_kltsa_thread((void*)&single_thread_param); 00128 #endif 00129 00130 // clean 00131 SG_FREE(local_gram_matrix); 00132 SG_FREE(ev_vector); 00133 SG_FREE(G_matrix); 00134 kernel_matrix.destroy_matrix(); 00135 00136 for (int32_t i=0; i<N; i++) 00137 { 00138 for (int32_t j=0; j<m_k; j++) 00139 W_matrix[N*neighborhood_matrix[j*N+i]+neighborhood_matrix[j*N+i]] += 1.0; 00140 } 00141 00142 return SGMatrix<float64_t>(W_matrix,N,N); 00143 } 00144 00145 void* CKernelLocalTangentSpaceAlignment::run_kltsa_thread(void* p) 00146 { 00147 KLTSA_THREAD_PARAM* parameters = (KLTSA_THREAD_PARAM*)p; 00148 int32_t idx_start = parameters->idx_start; 00149 int32_t idx_step = parameters->idx_step; 00150 int32_t idx_stop = parameters->idx_stop; 00151 int32_t m_k = parameters->m_k; 00152 int32_t target_dim = parameters->target_dim; 00153 int32_t N = parameters->N; 00154 const int32_t* neighborhood_matrix = parameters->neighborhood_matrix; 00155 const float64_t* kernel_matrix = parameters->kernel_matrix; 00156 float64_t* ev_vector = parameters->ev_vector; 00157 float64_t* G_matrix = parameters->G_matrix; 00158 float64_t* local_gram_matrix = parameters->local_gram_matrix; 00159 float64_t* W_matrix = parameters->W_matrix; 00160 #ifdef HAVE_PTHREAD 00161 PTHREAD_LOCK_T* W_matrix_lock = parameters->W_matrix_lock; 00162 #endif 00163 00164 int32_t i,j,k; 00165 00166 for (j=0; j<m_k; j++) 00167 G_matrix[j] = 1.0/CMath::sqrt((float64_t)m_k); 00168 00169 for (i=idx_start; i<idx_stop; i+=idx_step) 00170 { 00171 for (j=0; j<m_k; j++) 00172 { 00173 for (k=0; k<m_k; k++) 00174 local_gram_matrix[j*m_k+k] = kernel_matrix[neighborhood_matrix[j*N+i]*N+neighborhood_matrix[k*N+i]]; 00175 } 00176 00177 CMath::center_matrix(local_gram_matrix,m_k,m_k); 00178 00179 int32_t info = 0; 00180 wrap_dsyevr('V','U',m_k,local_gram_matrix,m_k,m_k-target_dim+1,m_k,ev_vector,G_matrix+m_k,&info); 00181 ASSERT(info==0); 00182 00183 cblas_dgemm(CblasColMajor,CblasNoTrans,CblasTrans, 00184 m_k,m_k,1+target_dim, 00185 1.0,G_matrix,m_k, 00186 G_matrix,m_k, 00187 0.0,local_gram_matrix,m_k); 00188 00189 #ifdef HAVE_PTHREAD 00190 PTHREAD_LOCK(W_matrix_lock); 00191 #endif 00192 for (j=0; j<m_k; j++) 00193 { 00194 for (k=0; k<m_k; k++) 00195 W_matrix[N*neighborhood_matrix[k*N+i]+neighborhood_matrix[j*N+i]] -= local_gram_matrix[j*m_k+k]; 00196 } 00197 #ifdef HAVE_PTHREAD 00198 PTHREAD_UNLOCK(W_matrix_lock); 00199 #endif 00200 } 00201 return NULL; 00202 } 00203 00204 #endif /* HAVE_LAPACK */