00001
00002
00003
00004
00005
00006
00007
00008
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
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
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
00196 CSimpleFeatures<float64_t>* simple_features = (CSimpleFeatures<float64_t>*) features;
00197 SG_REF(features);
00198
00199
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
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
00224 float64_t* W_matrix = distance_matrix.matrix;
00225 memset(W_matrix,0,sizeof(float64_t)*N*N);
00226
00227
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
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
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
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
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
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
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
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
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
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
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
00473
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
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
00491 if (m_k>dim)
00492 {
00493
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
00503
00504 clapack_dposv(CblasColMajor,CblasLower,m_k,1,covariance_matrix,m_k,id_vector,m_k);
00505
00506
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
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
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