KernelLocalTangentSpaceAlignment.cpp

Go to the documentation of this file.
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/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     // allocate threads and params
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     // init matrices and norm factor to be used
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*)&parameters[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     // clean
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 /* HAVE_LAPACK */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation