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/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
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
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
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