00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include <shogun/converter/HessianLocallyLinearEmbedding.h>
00012 #ifdef HAVE_LAPACK
00013 #include <shogun/mathematics/lapack.h>
00014 #include <shogun/lib/common.h>
00015 #include <shogun/mathematics/Math.h>
00016 #include <shogun/io/SGIO.h>
00017 #include <shogun/distance/Distance.h>
00018 #include <shogun/lib/Signal.h>
00019
00020 #ifdef HAVE_PTHREAD
00021 #include <pthread.h>
00022 #endif
00023
00024 using namespace shogun;
00025
00026 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00027 struct HESSIANESTIMATION_THREAD_PARAM
00028 {
00030 int32_t idx_start;
00032 int32_t idx_step;
00034 int32_t idx_stop;
00036 int32_t m_k;
00038 int32_t dim;
00040 int32_t N;
00042 int32_t dp;
00044 int32_t target_dim;
00046 const int32_t* neighborhood_matrix;
00048 const float64_t* feature_matrix;
00050 float64_t* local_feature_matrix;
00052 float64_t* Yi_matrix;
00054 float64_t* mean_vector;
00056 float64_t* s_values_vector;
00058 float64_t* tau;
00060 int32_t tau_len;
00062 float64_t* w_sum_vector;
00064 float64_t* q_matrix;
00066 float64_t* W_matrix;
00067 #ifdef HAVE_PTHREAD
00068
00069 PTHREAD_LOCK_T* W_matrix_lock;
00070 #endif
00071 };
00072 #endif
00073
00074 CHessianLocallyLinearEmbedding::CHessianLocallyLinearEmbedding() :
00075 CLocallyLinearEmbedding()
00076 {
00077 }
00078
00079 CHessianLocallyLinearEmbedding::~CHessianLocallyLinearEmbedding()
00080 {
00081 }
00082
00083 const char* CHessianLocallyLinearEmbedding::get_name() const
00084 {
00085 return "HessianLocallyLinearEmbedding";
00086 };
00087
00088 SGMatrix<float64_t> CHessianLocallyLinearEmbedding::construct_weight_matrix(CSimpleFeatures<float64_t>* simple_features,float64_t* W_matrix,
00089 SGMatrix<int32_t> neighborhood_matrix)
00090 {
00091 int32_t N = simple_features->get_num_vectors();
00092 int32_t dim = simple_features->get_num_features();
00093 int32_t dp = m_target_dim*(m_target_dim+1)/2;
00094 if (m_k<(1+m_target_dim+dp))
00095 SG_ERROR("K parameter should have value greater than 1+target dimensionality+dp.\n");
00096 int32_t t;
00097 #ifdef HAVE_PTHREAD
00098 int32_t num_threads = parallel->get_num_threads();
00099 ASSERT(num_threads>0);
00100
00101 pthread_t* threads = SG_MALLOC(pthread_t, num_threads);
00102 HESSIANESTIMATION_THREAD_PARAM* parameters = SG_MALLOC(HESSIANESTIMATION_THREAD_PARAM, num_threads);
00103
00104 #else
00105 int32_t num_threads = 1;
00106 #endif
00107
00108
00109 float64_t* local_feature_matrix = SG_MALLOC(float64_t, m_k*dim*num_threads);
00110 float64_t* s_values_vector = SG_MALLOC(float64_t, dim*num_threads);
00111 int32_t tau_len = CMath::min((1+m_target_dim+dp), m_k);
00112 float64_t* tau = SG_MALLOC(float64_t, tau_len*num_threads);
00113 float64_t* mean_vector = SG_MALLOC(float64_t, dim*num_threads);
00114 float64_t* q_matrix = SG_MALLOC(float64_t, m_k*m_k*num_threads);
00115 float64_t* w_sum_vector = SG_MALLOC(float64_t, dp*num_threads);
00116 float64_t* Yi_matrix = SG_MALLOC(float64_t, m_k*(1+m_target_dim+dp)*num_threads);
00117
00118 SGMatrix<float64_t> feature_matrix = simple_features->get_feature_matrix();
00119
00120 #ifdef HAVE_PTHREAD
00121 PTHREAD_LOCK_T W_matrix_lock;
00122 pthread_attr_t attr;
00123 PTHREAD_LOCK_INIT(&W_matrix_lock);
00124 pthread_attr_init(&attr);
00125 pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
00126
00127 for (t=0; t<num_threads; t++)
00128 {
00129 parameters[t].idx_start = t;
00130 parameters[t].idx_step = num_threads;
00131 parameters[t].idx_stop = N;
00132 parameters[t].m_k = m_k;
00133 parameters[t].dim = dim;
00134 parameters[t].target_dim = m_target_dim;
00135 parameters[t].N = N;
00136 parameters[t].dp = dp;
00137 parameters[t].neighborhood_matrix = neighborhood_matrix.matrix;
00138 parameters[t].feature_matrix = feature_matrix.matrix;
00139 parameters[t].local_feature_matrix = local_feature_matrix + (m_k*dim)*t;
00140 parameters[t].Yi_matrix = Yi_matrix + (m_k*(1+m_target_dim+dp))*t;
00141 parameters[t].mean_vector = mean_vector + dim*t;
00142 parameters[t].s_values_vector = s_values_vector + dim*t;
00143 parameters[t].tau = tau+tau_len*t;
00144 parameters[t].tau_len = tau_len;
00145 parameters[t].w_sum_vector = w_sum_vector + dp*t;
00146 parameters[t].q_matrix = q_matrix + (m_k*m_k)*t;
00147 parameters[t].W_matrix = W_matrix;
00148 parameters[t].W_matrix_lock = &W_matrix_lock;
00149 pthread_create(&threads[t], &attr, run_hessianestimation_thread, (void*)¶meters[t]);
00150 }
00151 for (t=0; t<num_threads; t++)
00152 pthread_join(threads[t], NULL);
00153 PTHREAD_LOCK_DESTROY(&W_matrix_lock);
00154 SG_FREE(parameters);
00155 SG_FREE(threads);
00156 #else
00157 HESSIANESTIMATION_THREAD_PARAM single_thread_param;
00158 single_thread_param.idx_start = t;
00159 single_thread_param.idx_step = num_threads;
00160 single_thread_param.idx_stop = N;
00161 single_thread_param.m_k = m_k;
00162 single_thread_param.dim = dim;
00163 single_thread_param.target_dim = m_target_dim;
00164 single_thread_param.N = N;
00165 single_thread_param.dp = dp;
00166 single_thread_param.neighborhood_matrix = neighborhood_matrix.matrix;
00167 single_thread_param.feature_matrix = feature_matrix.matrix;
00168 single_thread_param.local_feature_matrix = local_feature_matrix;
00169 single_thread_param.Yi_matrix = Yi_matrix;
00170 single_thread_param.mean_vector = mean_vector;
00171 single_thread_param.s_values_vector = s_values_vector;
00172 single_thread_param.tau = tau;
00173 single_thread_param.tau_len = tau_len;
00174 single_thread_param.w_sum_vector = w_sum_vector;
00175 single_thread_param.q_matrix = q_matrix;
00176 single_thread_param.W_matrix = W_matrix;
00177 run_hessianestimation_thread((void*)&single_thread_param);
00178 #endif
00179
00180
00181 SG_FREE(Yi_matrix);
00182 SG_FREE(s_values_vector);
00183 SG_FREE(mean_vector);
00184 SG_FREE(tau);
00185 SG_FREE(w_sum_vector);
00186 SG_FREE(local_feature_matrix);
00187 SG_FREE(q_matrix);
00188
00189 return SGMatrix<float64_t>(W_matrix,N,N);
00190 }
00191
00192 void* CHessianLocallyLinearEmbedding::run_hessianestimation_thread(void* p)
00193 {
00194 HESSIANESTIMATION_THREAD_PARAM* parameters = (HESSIANESTIMATION_THREAD_PARAM*)p;
00195 int32_t idx_start = parameters->idx_start;
00196 int32_t idx_step = parameters->idx_step;
00197 int32_t idx_stop = parameters->idx_stop;
00198 int32_t m_k = parameters->m_k;
00199 int32_t dim = parameters->dim;
00200 int32_t N = parameters->N;
00201 int32_t dp = parameters->dp;
00202 int32_t target_dim = parameters->target_dim;
00203 const int32_t* neighborhood_matrix = parameters->neighborhood_matrix;
00204 const float64_t* feature_matrix = parameters->feature_matrix;
00205 float64_t* local_feature_matrix = parameters->local_feature_matrix;
00206 float64_t* Yi_matrix = parameters->Yi_matrix;
00207 float64_t* mean_vector = parameters->mean_vector;
00208 float64_t* s_values_vector = parameters->s_values_vector;
00209 float64_t* tau = parameters->tau;
00210 int32_t tau_len = parameters->tau_len;
00211 float64_t* w_sum_vector = parameters->w_sum_vector;
00212 float64_t* q_matrix = parameters->q_matrix;
00213 float64_t* W_matrix = parameters->W_matrix;
00214 #ifdef HAVE_PTHREAD
00215 PTHREAD_LOCK_T* W_matrix_lock = parameters->W_matrix_lock;
00216 #endif
00217
00218 int i,j,k,l;
00219
00220 for (i=idx_start; i<idx_stop; i+=idx_step)
00221 {
00222
00223 for (j=0; j<m_k; j++)
00224 Yi_matrix[j] = 1.0;
00225
00226
00227 memset(mean_vector,0,sizeof(float64_t)*dim);
00228
00229
00230 for (j=0; j<m_k; j++)
00231 {
00232 for (k=0; k<dim; k++)
00233 local_feature_matrix[j*dim+k] = feature_matrix[neighborhood_matrix[j*N+i]*dim+k];
00234
00235 cblas_daxpy(dim,1.0,local_feature_matrix+j*dim,1,mean_vector,1);
00236 }
00237
00238
00239 cblas_dscal(dim,1.0/m_k,mean_vector,1);
00240
00241
00242 for (j=0; j<m_k; j++)
00243 cblas_daxpy(dim,-1.0,mean_vector,1,local_feature_matrix+j*dim,1);
00244
00245 int32_t info = 0;
00246
00247 wrap_dgesvd('N','O', dim,m_k,local_feature_matrix,dim,
00248 s_values_vector,
00249 NULL,1, NULL,1, &info);
00250 ASSERT(info==0);
00251
00252
00253 for (j=0; j<target_dim; j++)
00254 {
00255 for (k=0; k<m_k; k++)
00256 Yi_matrix[(j+1)*m_k+k] = local_feature_matrix[k*dim+j];
00257 }
00258
00259 int32_t ct = 0;
00260
00261
00262 for (j=0; j<target_dim; j++)
00263 {
00264 for (k=0; k<target_dim-j; k++)
00265 {
00266 for (l=0; l<m_k; l++)
00267 {
00268 Yi_matrix[(ct+k+1+target_dim)*m_k+l] = Yi_matrix[(j+1)*m_k+l]*Yi_matrix[(j+k+1)*m_k+l];
00269 }
00270 }
00271 ct += ct + target_dim - j;
00272 }
00273
00274
00275 wrap_dgeqrf(m_k,(1+target_dim+dp),Yi_matrix,m_k,tau,&info);
00276 ASSERT(info==0);
00277 wrap_dorgqr(m_k,(1+target_dim+dp),tau_len,Yi_matrix,m_k,tau,&info);
00278 ASSERT(info==0);
00279
00280 float64_t* Pii = (Yi_matrix+m_k*(1+target_dim));
00281
00282 for (j=0; j<dp; j++)
00283 {
00284 w_sum_vector[j] = 0.0;
00285 for (k=0; k<m_k; k++)
00286 {
00287 w_sum_vector[j] += Pii[j*m_k+k];
00288 }
00289 if (w_sum_vector[j]<0.001)
00290 w_sum_vector[j] = 1.0;
00291 for (k=0; k<m_k; k++)
00292 Pii[j*m_k+k] /= w_sum_vector[j];
00293 }
00294
00295 cblas_dgemm(CblasColMajor,CblasNoTrans,CblasTrans,
00296 m_k,m_k,dp,
00297 1.0,Pii,m_k,
00298 Pii,m_k,
00299 0.0,q_matrix,m_k);
00300 #ifdef HAVE_PTHREAD
00301 PTHREAD_LOCK(W_matrix_lock);
00302 #endif
00303 for (j=0; j<m_k; j++)
00304 {
00305 for (k=0; k<m_k; k++)
00306 W_matrix[N*neighborhood_matrix[k*N+i]+neighborhood_matrix[j*N+i]] += q_matrix[j*m_k+k];
00307 }
00308 #ifdef HAVE_PTHREAD
00309 PTHREAD_UNLOCK(W_matrix_lock);
00310 #endif
00311 }
00312 return NULL;
00313 }
00314 #endif