29 using namespace shogun;
31 #ifndef DOXYGEN_SHOULD_SKIP_THIS
32 struct LINRECONSTRUCTION_THREAD_PARAM
47 const int32_t* neighborhood_matrix;
62 struct NEIGHBORHOOD_THREAD_PARAM
79 int32_t* neighborhood_matrix;
102 "whether k should be determined automatically in range");
105 "maximum number of neighbors used to compute optimal one");
107 "nullspace finding regularization shift");
109 "shift used to regularize reconstruction step");
111 "whether arpack is being used or not");
183 return "LocallyLinearEmbedding";
193 SG_ERROR(
"Given features are not of SimpleRealFeatures type.\n");
202 SG_ERROR(
"Number of neighbors (%d) should be less than number of objects (%d).\n",
206 SG_DEBUG(
"Computing distance matrix\n");
211 SG_DEBUG(
"Calculating neighborhood matrix\n");
225 memset(W_matrix,0,
sizeof(
float64_t)*N*N);
228 SG_DEBUG(
"Constructing weight matrix\n");
248 if (right==left)
return left;
258 left_third = (left*2+right)/3;
259 right_third = (right*2+left)/3;
261 covariance_matrix,resid_vector,
262 id_vector,neighborhood_matrix);
264 covariance_matrix,resid_vector,
265 id_vector,neighborhood_matrix);
266 if (left_val<right_val)
290 cblas_dcopy(dim,feature_matrix+neighborhood_matrix[j*N+i]*dim,1,z_matrix+j*dim,1);
291 cblas_daxpy(dim,-1.0,feature_matrix+i*dim,1,z_matrix+j*dim,1);
293 cblas_dgemm(CblasColMajor,CblasTrans,CblasNoTrans,
297 0.0,covariance_matrix,k);
300 resid_vector[j] = 1.0;
307 trace += covariance_matrix[j*k+j];
308 for (j=0; j<
m_k; j++)
311 clapack_dposv(CblasColMajor,CblasLower,k,1,covariance_matrix,k,id_vector,k);
314 norming += id_vector[j];
315 cblas_dscal(k,-1.0/norming,id_vector,1);
316 cblas_dsymv(CblasColMajor,CblasLower,k,-1.0,covariance_matrix,k,id_vector,1,1.0,resid_vector,1);
317 total_residual_norm += cblas_dnrm2(k,resid_vector,1);
319 return total_residual_norm/k;
332 pthread_t* threads =
SG_MALLOC(pthread_t, num_threads);
333 LINRECONSTRUCTION_THREAD_PARAM* parameters =
SG_MALLOC(LINRECONSTRUCTION_THREAD_PARAM, num_threads);
335 pthread_attr_init(&attr);
336 pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
338 int32_t num_threads = 1;
349 for (t=0; t<num_threads; t++)
351 parameters[t].idx_start = t;
352 parameters[t].idx_step = num_threads;
353 parameters[t].idx_stop = N;
354 parameters[t].m_k =
m_k;
355 parameters[t].dim = dim;
357 parameters[t].neighborhood_matrix = neighborhood_matrix.
matrix;
358 parameters[t].z_matrix = z_matrix+(
m_k*dim)*t;
359 parameters[t].feature_matrix = feature_matrix.
matrix;
360 parameters[t].covariance_matrix = covariance_matrix+(
m_k*
m_k)*t;
361 parameters[t].id_vector = id_vector+
m_k*t;
362 parameters[t].W_matrix = W_matrix;
366 for (t=0; t<num_threads; t++)
367 pthread_join(threads[t], NULL);
368 pthread_attr_destroy(&attr);
372 LINRECONSTRUCTION_THREAD_PARAM single_thread_param;
373 single_thread_param.idx_start = 0;
374 single_thread_param.idx_step = 1;
375 single_thread_param.idx_stop = N;
376 single_thread_param.m_k =
m_k;
377 single_thread_param.dim = dim;
378 single_thread_param.N = N;
379 single_thread_param.neighborhood_matrix = neighborhood_matrix.
matrix;
380 single_thread_param.z_matrix = z_matrix;
381 single_thread_param.feature_matrix = feature_matrix.
matrix;
382 single_thread_param.covariance_matrix = covariance_matrix;
383 single_thread_param.id_vector = id_vector;
384 single_thread_param.W_matrix = W_matrix;
403 int eigenproblem_status = 0;
411 SG_ERROR(
"ARPACK is not supported in this configuration.\n");
416 arpack_dsxupd(matrix.
matrix,NULL,
false,N,dimension+1,
"LA",
true,3,
true,
m_nullspace_shift,0.0,
417 eigenvalues_vector,matrix.
matrix,eigenproblem_status);
419 if (eigenproblem_status)
420 SG_ERROR(
"ARPACK failed with code: %d", eigenproblem_status);
422 for (i=0; i<dimension; i++)
425 nullspace_features[j*dimension+i] = matrix[j*(dimension+1)+i+1];
434 wrap_dsyevr(
'V',
'U',N,matrix.
matrix,N,2,dimension+2,eigenvalues_vector,eigenvectors,&eigenproblem_status);
435 if (eigenproblem_status)
436 SG_ERROR(
"LAPACK failed with code: %d", eigenproblem_status);
439 for (i=0; i<dimension; i++)
442 nullspace_features[j*dimension+i] = eigenvectors[i*N+j];
452 LINRECONSTRUCTION_THREAD_PARAM* parameters = (LINRECONSTRUCTION_THREAD_PARAM*)p;
453 int32_t idx_start = parameters->idx_start;
454 int32_t idx_step = parameters->idx_step;
455 int32_t idx_stop = parameters->idx_stop;
456 int32_t
m_k = parameters->m_k;
457 int32_t dim = parameters->dim;
458 int32_t N = parameters->N;
459 const int32_t* neighborhood_matrix = parameters->neighborhood_matrix;
460 float64_t* z_matrix = parameters->z_matrix;
461 const float64_t* feature_matrix = parameters->feature_matrix;
462 float64_t* covariance_matrix = parameters->covariance_matrix;
463 float64_t* id_vector = parameters->id_vector;
464 float64_t* W_matrix = parameters->W_matrix;
470 for (i=idx_start; i<idx_stop; i+=idx_step)
474 for (j=0; j<
m_k; j++)
476 cblas_dcopy(dim,feature_matrix+neighborhood_matrix[j*N+i]*dim,1,z_matrix+j*dim,1);
477 cblas_daxpy(dim,-1.0,feature_matrix+i*dim,1,z_matrix+j*dim,1);
481 cblas_dgemm(CblasColMajor,CblasTrans,CblasNoTrans,
485 0.0,covariance_matrix,m_k);
487 for (j=0; j<
m_k; j++)
495 for (j=0; j<
m_k; j++)
496 trace += covariance_matrix[j*m_k+j];
498 for (j=0; j<
m_k; j++)
499 covariance_matrix[j*m_k+j] += m_reconstruction_shift*trace;
504 clapack_dposv(CblasColMajor,CblasLower,m_k,1,covariance_matrix,m_k,id_vector,m_k);
508 for (j=0; j<
m_k; j++)
509 norming += id_vector[j];
511 cblas_dscal(m_k,1.0/norming,id_vector,1);
513 memset(covariance_matrix,0,
sizeof(
float64_t)*m_k*m_k);
514 cblas_dger(CblasColMajor,m_k,m_k,1.0,id_vector,1,id_vector,1,covariance_matrix,m_k);
517 W_matrix[N*i+i] += 1.0;
518 for (j=0; j<
m_k; j++)
520 W_matrix[N*i+neighborhood_matrix[j*N+i]] -= id_vector[j];
521 W_matrix[N*neighborhood_matrix[j*N+i]+i] -= id_vector[j];
523 for (j=0; j<
m_k; j++)
525 for (k=0; k<
m_k; k++)
526 W_matrix[N*neighborhood_matrix[j*N+i]+neighborhood_matrix[k*N+i]]+=covariance_matrix[j*m_k+k];
535 int32_t N = distance_matrix.
num_rows;
537 int32_t* neighborhood_matrix =
SG_MALLOC(int32_t, N*k);
541 NEIGHBORHOOD_THREAD_PARAM* parameters =
SG_MALLOC(NEIGHBORHOOD_THREAD_PARAM, num_threads);
542 pthread_t* threads =
SG_MALLOC(pthread_t, num_threads);
544 pthread_attr_init(&attr);
545 pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
547 int32_t num_threads = 1;
549 CFibonacciHeap** heaps =
SG_MALLOC(CFibonacciHeap*, num_threads);
550 for (t=0; t<num_threads; t++)
551 heaps[t] =
new CFibonacciHeap(N);
554 for (t=0; t<num_threads; t++)
556 parameters[t].idx_start = t;
557 parameters[t].idx_step = num_threads;
558 parameters[t].idx_stop = N;
559 parameters[t].m_k = k;
561 parameters[t].heap = heaps[t];
562 parameters[t].neighborhood_matrix = neighborhood_matrix;
563 parameters[t].distance_matrix = distance_matrix.
matrix;
566 for (t=0; t<num_threads; t++)
567 pthread_join(threads[t], NULL);
568 pthread_attr_destroy(&attr);
572 NEIGHBORHOOD_THREAD_PARAM single_thread_param;
573 single_thread_param.idx_start = 0;
574 single_thread_param.idx_step = 1;
575 single_thread_param.idx_stop = N;
576 single_thread_param.m_k = k;
577 single_thread_param.N = N;
578 single_thread_param.heap = heaps[0]
579 single_thread_param.neighborhood_matrix = neighborhood_matrix;
580 single_thread_param.distance_matrix = distance_matrix.
matrix;
584 for (t=0; t<num_threads; t++)
593 NEIGHBORHOOD_THREAD_PARAM* parameters = (NEIGHBORHOOD_THREAD_PARAM*)p;
594 int32_t idx_start = parameters->idx_start;
595 int32_t idx_step = parameters->idx_step;
596 int32_t idx_stop = parameters->idx_stop;
597 int32_t N = parameters->N;
598 int32_t
m_k = parameters->m_k;
599 CFibonacciHeap* heap = parameters->heap;
600 const float64_t* distance_matrix = parameters->distance_matrix;
601 int32_t* neighborhood_matrix = parameters->neighborhood_matrix;
605 for (i=idx_start; i<idx_stop; i+=idx_step)
609 heap->insert(j,distance_matrix[i*N+j]);
612 heap->extract_min(tmp);
614 for (j=0; j<
m_k; j++)
615 neighborhood_matrix[j*N+i] = heap->extract_min(tmp);