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/LocallyLinearEmbedding.h> 00012 #include <shogun/lib/config.h> 00013 #ifdef HAVE_LAPACK 00014 #include <shogun/converter/EmbeddingConverter.h> 00015 #include <shogun/mathematics/arpack.h> 00016 #include <shogun/mathematics/lapack.h> 00017 #include <shogun/lib/FibonacciHeap.h> 00018 #include <shogun/base/DynArray.h> 00019 #include <shogun/mathematics/Math.h> 00020 #include <shogun/io/SGIO.h> 00021 #include <shogun/lib/Time.h> 00022 #include <shogun/distance/Distance.h> 00023 #include <shogun/lib/Signal.h> 00024 00025 #ifdef HAVE_PTHREAD 00026 #include <pthread.h> 00027 #endif 00028 00029 using namespace shogun; 00030 00031 #ifndef DOXYGEN_SHOULD_SKIP_THIS 00032 struct LINRECONSTRUCTION_THREAD_PARAM 00033 { 00035 int32_t idx_start; 00037 int32_t idx_stop; 00039 int32_t idx_step; 00041 int32_t m_k; 00043 int32_t dim; 00045 int32_t N; 00047 const int32_t* neighborhood_matrix; 00049 const float64_t* feature_matrix; 00051 float64_t* z_matrix; 00053 float64_t* covariance_matrix; 00055 float64_t* id_vector; 00057 float64_t* W_matrix; 00059 float64_t m_reconstruction_shift; 00060 }; 00061 00062 struct NEIGHBORHOOD_THREAD_PARAM 00063 { 00065 int32_t idx_start; 00067 int32_t idx_step; 00069 int32_t idx_stop; 00071 int32_t N; 00073 int32_t m_k; 00075 CFibonacciHeap* heap; 00077 const float64_t* distance_matrix; 00079 int32_t* neighborhood_matrix; 00080 }; 00081 #endif /* DOXYGEN_SHOULD_SKIP_THIS */ 00082 00083 CLocallyLinearEmbedding::CLocallyLinearEmbedding() : 00084 CEmbeddingConverter() 00085 { 00086 m_k = 10; 00087 m_max_k = 60; 00088 m_auto_k = false; 00089 m_nullspace_shift = -1e-9; 00090 m_reconstruction_shift = 1e-3; 00091 #ifdef HAVE_ARPACK 00092 m_use_arpack = true; 00093 #else 00094 m_use_arpack = false; 00095 #endif 00096 init(); 00097 } 00098 00099 void CLocallyLinearEmbedding::init() 00100 { 00101 m_parameters->add(&m_auto_k, "auto_k", 00102 "whether k should be determined automatically in range"); 00103 m_parameters->add(&m_k, "k", "number of neighbors"); 00104 m_parameters->add(&m_max_k, "max_k", 00105 "maximum number of neighbors used to compute optimal one"); 00106 m_parameters->add(&m_nullspace_shift, "nullspace_shift", 00107 "nullspace finding regularization shift"); 00108 m_parameters->add(&m_reconstruction_shift, "reconstruction_shift", 00109 "shift used to regularize reconstruction step"); 00110 m_parameters->add(&m_use_arpack, "use_arpack", 00111 "whether arpack is being used or not"); 00112 } 00113 00114 00115 CLocallyLinearEmbedding::~CLocallyLinearEmbedding() 00116 { 00117 } 00118 00119 void CLocallyLinearEmbedding::set_k(int32_t k) 00120 { 00121 ASSERT(k>0); 00122 m_k = k; 00123 } 00124 00125 int32_t CLocallyLinearEmbedding::get_k() const 00126 { 00127 return m_k; 00128 } 00129 00130 void CLocallyLinearEmbedding::set_max_k(int32_t max_k) 00131 { 00132 ASSERT(max_k>=m_k); 00133 m_max_k = max_k; 00134 } 00135 00136 int32_t CLocallyLinearEmbedding::get_max_k() const 00137 { 00138 return m_max_k; 00139 } 00140 00141 void CLocallyLinearEmbedding::set_auto_k(bool auto_k) 00142 { 00143 m_auto_k = auto_k; 00144 } 00145 00146 bool CLocallyLinearEmbedding::get_auto_k() const 00147 { 00148 return m_auto_k; 00149 } 00150 00151 void CLocallyLinearEmbedding::set_nullspace_shift(float64_t nullspace_shift) 00152 { 00153 m_nullspace_shift = nullspace_shift; 00154 } 00155 00156 float64_t CLocallyLinearEmbedding::get_nullspace_shift() const 00157 { 00158 return m_nullspace_shift; 00159 } 00160 00161 void CLocallyLinearEmbedding::set_reconstruction_shift(float64_t reconstruction_shift) 00162 { 00163 m_reconstruction_shift = reconstruction_shift; 00164 } 00165 00166 float64_t CLocallyLinearEmbedding::get_reconstruction_shift() const 00167 { 00168 return m_reconstruction_shift; 00169 } 00170 00171 void CLocallyLinearEmbedding::set_use_arpack(bool use_arpack) 00172 { 00173 m_use_arpack = use_arpack; 00174 } 00175 00176 bool CLocallyLinearEmbedding::get_use_arpack() const 00177 { 00178 return m_use_arpack; 00179 } 00180 00181 const char* CLocallyLinearEmbedding::get_name() const 00182 { 00183 return "LocallyLinearEmbedding"; 00184 } 00185 00186 CFeatures* CLocallyLinearEmbedding::apply(CFeatures* features) 00187 { 00188 ASSERT(features); 00189 // check features 00190 if (!(features->get_feature_class()==C_SIMPLE && 00191 features->get_feature_type()==F_DREAL)) 00192 { 00193 SG_ERROR("Given features are not of SimpleRealFeatures type.\n"); 00194 } 00195 // shorthand for simplefeatures 00196 CSimpleFeatures<float64_t>* simple_features = (CSimpleFeatures<float64_t>*) features; 00197 SG_REF(features); 00198 00199 // get and check number of vectors 00200 int32_t N = simple_features->get_num_vectors(); 00201 if (m_k>=N) 00202 SG_ERROR("Number of neighbors (%d) should be less than number of objects (%d).\n", 00203 m_k, N); 00204 00205 // compute distance matrix 00206 SG_DEBUG("Computing distance matrix\n"); 00207 ASSERT(m_distance); 00208 m_distance->init(simple_features,simple_features); 00209 SGMatrix<float64_t> distance_matrix = m_distance->get_distance_matrix(); 00210 m_distance->remove_lhs_and_rhs(); 00211 SG_DEBUG("Calculating neighborhood matrix\n"); 00212 SGMatrix<int32_t> neighborhood_matrix; 00213 00214 if (m_auto_k) 00215 { 00216 neighborhood_matrix = get_neighborhood_matrix(distance_matrix,m_max_k); 00217 m_k = estimate_k(simple_features,neighborhood_matrix); 00218 SG_DEBUG("Estimated k with value of %d\n",m_k); 00219 } 00220 else 00221 neighborhood_matrix = get_neighborhood_matrix(distance_matrix,m_k); 00222 00223 // init W (weight) matrix 00224 float64_t* W_matrix = distance_matrix.matrix; 00225 memset(W_matrix,0,sizeof(float64_t)*N*N); 00226 00227 // construct weight matrix 00228 SG_DEBUG("Constructing weight matrix\n"); 00229 SGMatrix<float64_t> weight_matrix = construct_weight_matrix(simple_features,W_matrix,neighborhood_matrix); 00230 neighborhood_matrix.destroy_matrix(); 00231 00232 // find null space of weight matrix 00233 SG_DEBUG("Finding nullspace\n"); 00234 SGMatrix<float64_t> new_feature_matrix = construct_embedding(weight_matrix,m_target_dim); 00235 weight_matrix.destroy_matrix(); 00236 00237 SG_UNREF(features); 00238 return (CFeatures*)(new CSimpleFeatures<float64_t>(new_feature_matrix)); 00239 } 00240 00241 int32_t CLocallyLinearEmbedding::estimate_k(CSimpleFeatures<float64_t>* simple_features, SGMatrix<int32_t> neighborhood_matrix) 00242 { 00243 int32_t right = m_max_k; 00244 int32_t left = m_k; 00245 int32_t left_third; 00246 int32_t right_third; 00247 ASSERT(right>=left); 00248 if (right==left) return left; 00249 int32_t dim; 00250 int32_t N; 00251 float64_t* feature_matrix= simple_features->get_feature_matrix(dim,N); 00252 float64_t* z_matrix = SG_MALLOC(float64_t,right*dim); 00253 float64_t* covariance_matrix = SG_MALLOC(float64_t,right*right); 00254 float64_t* resid_vector = SG_MALLOC(float64_t, right); 00255 float64_t* id_vector = SG_MALLOC(float64_t, right); 00256 while (right-left>2) 00257 { 00258 left_third = (left*2+right)/3; 00259 right_third = (right*2+left)/3; 00260 float64_t left_val = compute_reconstruction_error(left_third,dim,N,feature_matrix,z_matrix, 00261 covariance_matrix,resid_vector, 00262 id_vector,neighborhood_matrix); 00263 float64_t right_val = compute_reconstruction_error(right_third,dim,N,feature_matrix,z_matrix, 00264 covariance_matrix,resid_vector, 00265 id_vector,neighborhood_matrix); 00266 if (left_val<right_val) 00267 right = right_third; 00268 else 00269 left = left_third; 00270 } 00271 SG_FREE(z_matrix); 00272 SG_FREE(covariance_matrix); 00273 SG_FREE(resid_vector); 00274 SG_FREE(id_vector); 00275 return right; 00276 } 00277 00278 float64_t CLocallyLinearEmbedding::compute_reconstruction_error(int32_t k, int dim, int N, float64_t* feature_matrix, 00279 float64_t* z_matrix, float64_t* covariance_matrix, 00280 float64_t* resid_vector, float64_t* id_vector, 00281 SGMatrix<int32_t> neighborhood_matrix) 00282 { 00283 // todo parse params 00284 int32_t i,j; 00285 float64_t total_residual_norm = 0.0; 00286 for (i=CMath::random(0,20); i<N; i+=N/20) 00287 { 00288 for (j=0; j<k; j++) 00289 { 00290 cblas_dcopy(dim,feature_matrix+neighborhood_matrix[j*N+i]*dim,1,z_matrix+j*dim,1); 00291 cblas_daxpy(dim,-1.0,feature_matrix+i*dim,1,z_matrix+j*dim,1); 00292 } 00293 cblas_dgemm(CblasColMajor,CblasTrans,CblasNoTrans, 00294 k,k,dim, 00295 1.0,z_matrix,dim, 00296 z_matrix,dim, 00297 0.0,covariance_matrix,k); 00298 for (j=0; j<k; j++) 00299 { 00300 resid_vector[j] = 1.0; 00301 id_vector[j] = 1.0; 00302 } 00303 if (k>dim) 00304 { 00305 float64_t trace = 0.0; 00306 for (j=0; j<k; j++) 00307 trace += covariance_matrix[j*k+j]; 00308 for (j=0; j<m_k; j++) 00309 covariance_matrix[j*k+j] += m_reconstruction_shift*trace; 00310 } 00311 clapack_dposv(CblasColMajor,CblasLower,k,1,covariance_matrix,k,id_vector,k); 00312 float64_t norming=0.0; 00313 for (j=0; j<k; j++) 00314 norming += id_vector[j]; 00315 cblas_dscal(k,-1.0/norming,id_vector,1); 00316 cblas_dsymv(CblasColMajor,CblasLower,k,-1.0,covariance_matrix,k,id_vector,1,1.0,resid_vector,1); 00317 total_residual_norm += cblas_dnrm2(k,resid_vector,1); 00318 } 00319 return total_residual_norm/k; 00320 } 00321 00322 SGMatrix<float64_t> CLocallyLinearEmbedding::construct_weight_matrix(CSimpleFeatures<float64_t>* simple_features, 00323 float64_t* W_matrix, SGMatrix<int32_t> neighborhood_matrix) 00324 { 00325 int32_t N = simple_features->get_num_vectors(); 00326 int32_t dim = simple_features->get_num_features(); 00327 int32_t t; 00328 #ifdef HAVE_PTHREAD 00329 int32_t num_threads = parallel->get_num_threads(); 00330 ASSERT(num_threads>0); 00331 // allocate threads 00332 pthread_t* threads = SG_MALLOC(pthread_t, num_threads); 00333 LINRECONSTRUCTION_THREAD_PARAM* parameters = SG_MALLOC(LINRECONSTRUCTION_THREAD_PARAM, num_threads); 00334 pthread_attr_t attr; 00335 pthread_attr_init(&attr); 00336 pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE); 00337 #else 00338 int32_t num_threads = 1; 00339 #endif 00340 // init storages to be used 00341 float64_t* z_matrix = SG_MALLOC(float64_t, m_k*dim*num_threads); 00342 float64_t* covariance_matrix = SG_MALLOC(float64_t, m_k*m_k*num_threads); 00343 float64_t* id_vector = SG_MALLOC(float64_t, m_k*num_threads); 00344 00345 // get feature matrix 00346 SGMatrix<float64_t> feature_matrix = simple_features->get_feature_matrix(); 00347 00348 #ifdef HAVE_PTHREAD 00349 for (t=0; t<num_threads; t++) 00350 { 00351 parameters[t].idx_start = t; 00352 parameters[t].idx_step = num_threads; 00353 parameters[t].idx_stop = N; 00354 parameters[t].m_k = m_k; 00355 parameters[t].dim = dim; 00356 parameters[t].N = N; 00357 parameters[t].neighborhood_matrix = neighborhood_matrix.matrix; 00358 parameters[t].z_matrix = z_matrix+(m_k*dim)*t; 00359 parameters[t].feature_matrix = feature_matrix.matrix; 00360 parameters[t].covariance_matrix = covariance_matrix+(m_k*m_k)*t; 00361 parameters[t].id_vector = id_vector+m_k*t; 00362 parameters[t].W_matrix = W_matrix; 00363 parameters[t].m_reconstruction_shift = m_reconstruction_shift; 00364 pthread_create(&threads[t], &attr, run_linearreconstruction_thread, (void*)¶meters[t]); 00365 } 00366 for (t=0; t<num_threads; t++) 00367 pthread_join(threads[t], NULL); 00368 pthread_attr_destroy(&attr); 00369 SG_FREE(parameters); 00370 SG_FREE(threads); 00371 #else 00372 LINRECONSTRUCTION_THREAD_PARAM single_thread_param; 00373 single_thread_param.idx_start = 0; 00374 single_thread_param.idx_step = 1; 00375 single_thread_param.idx_stop = N; 00376 single_thread_param.m_k = m_k; 00377 single_thread_param.dim = dim; 00378 single_thread_param.N = N; 00379 single_thread_param.neighborhood_matrix = neighborhood_matrix.matrix; 00380 single_thread_param.z_matrix = z_matrix; 00381 single_thread_param.feature_matrix = feature_matrix.matrix; 00382 single_thread_param.covariance_matrix = covariance_matrix; 00383 single_thread_param.id_vector = id_vector; 00384 single_thread_param.W_matrix = W_matrix; 00385 single_thread_param.m_reconstruction_shift = m_reconstruction_shift; 00386 run_linearreconstruction_thread((void*)single_thread_param); 00387 #endif 00388 00389 // clean 00390 SG_FREE(id_vector); 00391 SG_FREE(z_matrix); 00392 SG_FREE(covariance_matrix); 00393 00394 return SGMatrix<float64_t>(W_matrix,N,N); 00395 } 00396 00397 SGMatrix<float64_t> CLocallyLinearEmbedding::construct_embedding(SGMatrix<float64_t> matrix,int dimension) 00398 { 00399 int i,j; 00400 ASSERT(matrix.num_cols==matrix.num_rows); 00401 int N = matrix.num_cols; 00402 // get eigenvectors with ARPACK or LAPACK 00403 int eigenproblem_status = 0; 00404 00405 float64_t* eigenvalues_vector = NULL; 00406 float64_t* eigenvectors = NULL; 00407 float64_t* nullspace_features = NULL; 00408 if (m_use_arpack) 00409 { 00410 #ifndef HAVE_ARPACK 00411 SG_ERROR("ARPACK is not supported in this configuration.\n"); 00412 #endif 00413 // using ARPACK (faster) 00414 eigenvalues_vector = SG_MALLOC(float64_t, dimension+1); 00415 #ifdef HAVE_ARPACK 00416 arpack_dsxupd(matrix.matrix,NULL,false,N,dimension+1,"LA",true,3,true,m_nullspace_shift,0.0, 00417 eigenvalues_vector,matrix.matrix,eigenproblem_status); 00418 #endif 00419 if (eigenproblem_status) 00420 SG_ERROR("ARPACK failed with code: %d", eigenproblem_status); 00421 nullspace_features = SG_MALLOC(float64_t, N*dimension); 00422 for (i=0; i<dimension; i++) 00423 { 00424 for (j=0; j<N; j++) 00425 nullspace_features[j*dimension+i] = matrix[j*(dimension+1)+i+1]; 00426 } 00427 SG_FREE(eigenvalues_vector); 00428 } 00429 else 00430 { 00431 // using LAPACK (slower) 00432 eigenvalues_vector = SG_MALLOC(float64_t, N); 00433 eigenvectors = SG_MALLOC(float64_t,(dimension+1)*N); 00434 wrap_dsyevr('V','U',N,matrix.matrix,N,2,dimension+2,eigenvalues_vector,eigenvectors,&eigenproblem_status); 00435 if (eigenproblem_status) 00436 SG_ERROR("LAPACK failed with code: %d", eigenproblem_status); 00437 nullspace_features = SG_MALLOC(float64_t, N*dimension); 00438 // LAPACKed eigenvectors 00439 for (i=0; i<dimension; i++) 00440 { 00441 for (j=0; j<N; j++) 00442 nullspace_features[j*dimension+i] = eigenvectors[i*N+j]; 00443 } 00444 SG_FREE(eigenvectors); 00445 SG_FREE(eigenvalues_vector); 00446 } 00447 return SGMatrix<float64_t>(nullspace_features,dimension,N); 00448 } 00449 00450 void* CLocallyLinearEmbedding::run_linearreconstruction_thread(void* p) 00451 { 00452 LINRECONSTRUCTION_THREAD_PARAM* parameters = (LINRECONSTRUCTION_THREAD_PARAM*)p; 00453 int32_t idx_start = parameters->idx_start; 00454 int32_t idx_step = parameters->idx_step; 00455 int32_t idx_stop = parameters->idx_stop; 00456 int32_t m_k = parameters->m_k; 00457 int32_t dim = parameters->dim; 00458 int32_t N = parameters->N; 00459 const int32_t* neighborhood_matrix = parameters->neighborhood_matrix; 00460 float64_t* z_matrix = parameters->z_matrix; 00461 const float64_t* feature_matrix = parameters->feature_matrix; 00462 float64_t* covariance_matrix = parameters->covariance_matrix; 00463 float64_t* id_vector = parameters->id_vector; 00464 float64_t* W_matrix = parameters->W_matrix; 00465 float64_t m_reconstruction_shift = parameters->m_reconstruction_shift; 00466 00467 int32_t i,j,k; 00468 float64_t norming,trace; 00469 00470 for (i=idx_start; i<idx_stop; i+=idx_step) 00471 { 00472 // compute local feature matrix containing neighbors of i-th vector 00473 // center features by subtracting i-th feature column 00474 for (j=0; j<m_k; j++) 00475 { 00476 cblas_dcopy(dim,feature_matrix+neighborhood_matrix[j*N+i]*dim,1,z_matrix+j*dim,1); 00477 cblas_daxpy(dim,-1.0,feature_matrix+i*dim,1,z_matrix+j*dim,1); 00478 } 00479 00480 // compute local covariance matrix 00481 cblas_dgemm(CblasColMajor,CblasTrans,CblasNoTrans, 00482 m_k,m_k,dim, 00483 1.0,z_matrix,dim, 00484 z_matrix,dim, 00485 0.0,covariance_matrix,m_k); 00486 00487 for (j=0; j<m_k; j++) 00488 id_vector[j] = 1.0; 00489 00490 // regularize in case of ill-posed system 00491 if (m_k>dim) 00492 { 00493 // compute tr(C) 00494 trace = 0.0; 00495 for (j=0; j<m_k; j++) 00496 trace += covariance_matrix[j*m_k+j]; 00497 00498 for (j=0; j<m_k; j++) 00499 covariance_matrix[j*m_k+j] += m_reconstruction_shift*trace; 00500 } 00501 00502 // solve system of linear equations: covariance_matrix * X = 1 00503 // covariance_matrix is a pos-def matrix 00504 clapack_dposv(CblasColMajor,CblasLower,m_k,1,covariance_matrix,m_k,id_vector,m_k); 00505 00506 // normalize weights 00507 norming=0.0; 00508 for (j=0; j<m_k; j++) 00509 norming += id_vector[j]; 00510 00511 cblas_dscal(m_k,1.0/norming,id_vector,1); 00512 00513 memset(covariance_matrix,0,sizeof(float64_t)*m_k*m_k); 00514 cblas_dger(CblasColMajor,m_k,m_k,1.0,id_vector,1,id_vector,1,covariance_matrix,m_k); 00515 00516 // put weights into W matrix 00517 W_matrix[N*i+i] += 1.0; 00518 for (j=0; j<m_k; j++) 00519 { 00520 W_matrix[N*i+neighborhood_matrix[j*N+i]] -= id_vector[j]; 00521 W_matrix[N*neighborhood_matrix[j*N+i]+i] -= id_vector[j]; 00522 } 00523 for (j=0; j<m_k; j++) 00524 { 00525 for (k=0; k<m_k; k++) 00526 W_matrix[N*neighborhood_matrix[j*N+i]+neighborhood_matrix[k*N+i]]+=covariance_matrix[j*m_k+k]; 00527 } 00528 } 00529 return NULL; 00530 } 00531 00532 SGMatrix<int32_t> CLocallyLinearEmbedding::get_neighborhood_matrix(SGMatrix<float64_t> distance_matrix, int32_t k) 00533 { 00534 int32_t t; 00535 int32_t N = distance_matrix.num_rows; 00536 // init matrix and heap to be used 00537 int32_t* neighborhood_matrix = SG_MALLOC(int32_t, N*k); 00538 #ifdef HAVE_PTHREAD 00539 int32_t num_threads = parallel->get_num_threads(); 00540 ASSERT(num_threads>0); 00541 NEIGHBORHOOD_THREAD_PARAM* parameters = SG_MALLOC(NEIGHBORHOOD_THREAD_PARAM, num_threads); 00542 pthread_t* threads = SG_MALLOC(pthread_t, num_threads); 00543 pthread_attr_t attr; 00544 pthread_attr_init(&attr); 00545 pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE); 00546 #else 00547 int32_t num_threads = 1; 00548 #endif 00549 CFibonacciHeap** heaps = SG_MALLOC(CFibonacciHeap*, num_threads); 00550 for (t=0; t<num_threads; t++) 00551 heaps[t] = new CFibonacciHeap(N); 00552 00553 #ifdef HAVE_PTHREAD 00554 for (t=0; t<num_threads; t++) 00555 { 00556 parameters[t].idx_start = t; 00557 parameters[t].idx_step = num_threads; 00558 parameters[t].idx_stop = N; 00559 parameters[t].m_k = k; 00560 parameters[t].N = N; 00561 parameters[t].heap = heaps[t]; 00562 parameters[t].neighborhood_matrix = neighborhood_matrix; 00563 parameters[t].distance_matrix = distance_matrix.matrix; 00564 pthread_create(&threads[t], &attr, run_neighborhood_thread, (void*)¶meters[t]); 00565 } 00566 for (t=0; t<num_threads; t++) 00567 pthread_join(threads[t], NULL); 00568 pthread_attr_destroy(&attr); 00569 SG_FREE(threads); 00570 SG_FREE(parameters); 00571 #else 00572 NEIGHBORHOOD_THREAD_PARAM single_thread_param; 00573 single_thread_param.idx_start = 0; 00574 single_thread_param.idx_step = 1; 00575 single_thread_param.idx_stop = N; 00576 single_thread_param.m_k = k; 00577 single_thread_param.N = N; 00578 single_thread_param.heap = heaps[0] 00579 single_thread_param.neighborhood_matrix = neighborhood_matrix; 00580 single_thread_param.distance_matrix = distance_matrix.matrix; 00581 run_neighborhood_thread((void*)&single_thread_param); 00582 #endif 00583 00584 for (t=0; t<num_threads; t++) 00585 delete heaps[t]; 00586 SG_FREE(heaps); 00587 00588 return SGMatrix<int32_t>(neighborhood_matrix,k,N); 00589 } 00590 00591 void* CLocallyLinearEmbedding::run_neighborhood_thread(void* p) 00592 { 00593 NEIGHBORHOOD_THREAD_PARAM* parameters = (NEIGHBORHOOD_THREAD_PARAM*)p; 00594 int32_t idx_start = parameters->idx_start; 00595 int32_t idx_step = parameters->idx_step; 00596 int32_t idx_stop = parameters->idx_stop; 00597 int32_t N = parameters->N; 00598 int32_t m_k = parameters->m_k; 00599 CFibonacciHeap* heap = parameters->heap; 00600 const float64_t* distance_matrix = parameters->distance_matrix; 00601 int32_t* neighborhood_matrix = parameters->neighborhood_matrix; 00602 00603 int32_t i,j; 00604 float64_t tmp; 00605 for (i=idx_start; i<idx_stop; i+=idx_step) 00606 { 00607 for (j=0; j<N; j++) 00608 { 00609 heap->insert(j,distance_matrix[i*N+j]); 00610 } 00611 00612 heap->extract_min(tmp); 00613 00614 for (j=0; j<m_k; j++) 00615 neighborhood_matrix[j*N+i] = heap->extract_min(tmp); 00616 00617 heap->clear(); 00618 } 00619 00620 return NULL; 00621 } 00622 #endif /* HAVE_LAPACK */