00001
00002
00003
00004
00005
00006
00007
00008
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/lib/SGMatrix.h>
00018 #include <shogun/mathematics/Math.h>
00019 #include <shogun/io/SGIO.h>
00020 #include <shogun/distance/Distance.h>
00021 #include <shogun/lib/Signal.h>
00022 #include <shogun/base/Parallel.h>
00023
00024 #ifdef HAVE_PTHREAD
00025 #include <pthread.h>
00026 #endif
00027
00028 using namespace shogun;
00029
00030 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00031 struct KLTSA_THREAD_PARAM
00032 {
00034 int32_t idx_start;
00036 int32_t idx_step;
00038 int32_t idx_stop;
00040 int32_t m_k;
00042 int32_t target_dim;
00044 int32_t N;
00046 const int32_t* neighborhood_matrix;
00048 const float64_t* kernel_matrix;
00050 float64_t* local_gram_matrix;
00052 float64_t* ev_vector;
00054 float64_t* G_matrix;
00056 float64_t* W_matrix;
00057 #ifdef HAVE_PTHREAD
00058
00059 PTHREAD_LOCK_T* W_matrix_lock;
00060 #endif
00061 };
00062 #endif
00063
00064 CKernelLocalTangentSpaceAlignment::CKernelLocalTangentSpaceAlignment() :
00065 CKernelLocallyLinearEmbedding()
00066 {
00067 }
00068
00069 CKernelLocalTangentSpaceAlignment::CKernelLocalTangentSpaceAlignment(CKernel* kernel) :
00070 CKernelLocallyLinearEmbedding(kernel)
00071 {
00072 }
00073
00074 CKernelLocalTangentSpaceAlignment::~CKernelLocalTangentSpaceAlignment()
00075 {
00076 }
00077
00078 const char* CKernelLocalTangentSpaceAlignment::get_name() const
00079 {
00080 return "KernelLocalTangentSpaceAlignment";
00081 };
00082
00083 SGMatrix<float64_t> CKernelLocalTangentSpaceAlignment::construct_weight_matrix(SGMatrix<float64_t> kernel_matrix,
00084 SGMatrix<int32_t> neighborhood_matrix)
00085 {
00086 int32_t N = kernel_matrix.num_cols;
00087 int32_t t;
00088 #ifdef HAVE_PTHREAD
00089 int32_t num_threads = parallel->get_num_threads();
00090 ASSERT(num_threads>0);
00091
00092 pthread_t* threads = SG_MALLOC(pthread_t, num_threads);
00093 KLTSA_THREAD_PARAM* parameters = SG_MALLOC(KLTSA_THREAD_PARAM, num_threads);
00094 #else
00095 int32_t num_threads = 1;
00096 #endif
00097
00098
00099 float64_t* local_gram_matrix = SG_MALLOC(float64_t, m_k*m_k*num_threads);
00100 float64_t* G_matrix = SG_MALLOC(float64_t, m_k*(1+m_target_dim)*num_threads);
00101 float64_t* W_matrix = SG_CALLOC(float64_t, N*N);
00102 float64_t* ev_vector = SG_MALLOC(float64_t, m_k*num_threads);
00103
00104 #ifdef HAVE_PTHREAD
00105 PTHREAD_LOCK_T W_matrix_lock;
00106 pthread_attr_t attr;
00107 PTHREAD_LOCK_INIT(&W_matrix_lock);
00108 pthread_attr_init(&attr);
00109 pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
00110
00111 for (t=0; t<num_threads; t++)
00112 {
00113 KLTSA_THREAD_PARAM params = {t,num_threads,N,m_k,m_target_dim,N,neighborhood_matrix.matrix,
00114 kernel_matrix.matrix,local_gram_matrix+(m_k*m_k)*t,ev_vector+m_k*t,
00115 G_matrix+(m_k*(1+m_target_dim))*t,W_matrix,&W_matrix_lock};
00116 parameters[t] = params;
00117 pthread_create(&threads[t], &attr, run_kltsa_thread, (void*)¶meters[t]);
00118 }
00119 for (t=0; t<num_threads; t++)
00120 pthread_join(threads[t], NULL);
00121 PTHREAD_LOCK_DESTROY(&W_matrix_lock);
00122 SG_FREE(parameters);
00123 SG_FREE(threads);
00124 #else
00125 KLTSA_THREAD_PARAM single_thread_param = {0,1,N,m_k,m_target_dim,N,neighborhood_matrix.matrix,
00126 kernel_matrix.matrix,local_gram_matrix,ev_vector,
00127 G_matrix,W_matrix};
00128 run_kltsa_thread((void*)&single_thread_param);
00129 #endif
00130
00131
00132 SG_FREE(local_gram_matrix);
00133 SG_FREE(ev_vector);
00134 SG_FREE(G_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[i*m_k+j]+neighborhood_matrix[i*m_k+j]] += 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[i*m_k+j]*N+neighborhood_matrix[i*m_k+k]];
00175 }
00176
00177 SGMatrix<float64_t>::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[i*m_k+k]+neighborhood_matrix[i*m_k+j]] -= 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