30 using namespace shogun;
32 #ifndef DOXYGEN_SHOULD_SKIP_THIS
33 struct LINRECONSTRUCTION_THREAD_PARAM
44 const int32_t* neighborhood_matrix;
65 PTHREAD_LOCK_T* W_matrix_lock;
69 class LLE_COVERTREE_POINT
76 distance_matrix = dmatrix;
79 inline double distance(
const LLE_COVERTREE_POINT& p)
const
81 return distance_matrix[point_index*distance_matrix.num_rows+p.point_index];
84 inline bool operator==(
const LLE_COVERTREE_POINT& p)
const
86 return (p.point_index==point_index);
114 "whether k should be determined automatically in range",
MS_AVAILABLE);
117 "maximum number of neighbors used to compute optimal one",
MS_AVAILABLE);
119 "nullspace finding regularization shift");
195 return "LocallyLinearEmbedding";
205 SG_ERROR(
"Given features are not of SimpleRealFeatures type.\n");
214 SG_ERROR(
"Number of neighbors (%d) should be less than number of objects (%d).\n",
218 SG_DEBUG(
"Computing distance matrix\n");
226 SG_DEBUG(
"Calculating neighborhood matrix\n");
245 SG_DEBUG(
"Constructing weight matrix\n");
268 if (right==left)
return left;
278 left_third = (left*2+right)/3;
279 right_third = (right*2+left)/3;
281 covariance_matrix,resid_vector,
282 id_vector,neighborhood_matrix);
284 covariance_matrix,resid_vector,
285 id_vector,neighborhood_matrix);
286 if (left_val<right_val)
310 cblas_dcopy(dim,feature_matrix+neighborhood_matrix[i*
m_k+j]*dim,1,z_matrix+j*dim,1);
311 cblas_daxpy(dim,-1.0,feature_matrix+i*dim,1,z_matrix+j*dim,1);
313 cblas_dgemm(CblasColMajor,CblasTrans,CblasNoTrans,
317 0.0,covariance_matrix,k);
320 resid_vector[j] = 1.0;
327 trace += covariance_matrix[j*k+j];
328 for (j=0; j<
m_k; j++)
331 clapack_dposv(CblasColMajor,CblasLower,k,1,covariance_matrix,k,id_vector,k);
334 norming += id_vector[j];
335 cblas_dscal(k,-1.0/norming,id_vector,1);
336 cblas_dsymv(CblasColMajor,CblasLower,k,-1.0,covariance_matrix,k,id_vector,1,1.0,resid_vector,1);
337 total_residual_norm += cblas_dnrm2(k,resid_vector,1);
339 return total_residual_norm/k;
352 pthread_t* threads =
SG_MALLOC(pthread_t, num_threads);
353 LINRECONSTRUCTION_THREAD_PARAM* parameters =
SG_MALLOC(LINRECONSTRUCTION_THREAD_PARAM, num_threads);
355 pthread_attr_init(&attr);
356 pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
358 int32_t num_threads = 1;
369 PTHREAD_LOCK_T W_matrix_lock;
370 PTHREAD_LOCK_INIT(&W_matrix_lock);
371 for (t=0; t<num_threads; t++)
373 parameters[t].idx_start = t;
374 parameters[t].idx_step = num_threads;
375 parameters[t].idx_stop = N;
376 parameters[t].m_k =
m_k;
377 parameters[t].neighborhood_matrix = neighborhood_matrix.
matrix;
378 parameters[t].z_matrix = z_matrix+(
m_k*dim)*t;
379 parameters[t].feature_matrix = feature_matrix.
matrix;
380 parameters[t].covariance_matrix = covariance_matrix+(
m_k*
m_k)*t;
381 parameters[t].id_vector = id_vector+
m_k*t;
382 parameters[t].W_matrix = W_matrix;
384 parameters[t].dim = feature_matrix.
num_rows;
385 parameters[t].N = feature_matrix.
num_cols;
386 parameters[t].actual_k = neighborhood_matrix.
num_rows;
387 parameters[t].W_matrix_lock = &W_matrix_lock;
390 for (t=0; t<num_threads; t++)
391 pthread_join(threads[t], NULL);
392 pthread_attr_destroy(&attr);
396 LINRECONSTRUCTION_THREAD_PARAM single_thread_param;
397 single_thread_param.idx_start = 0;
398 single_thread_param.idx_step = 1;
399 single_thread_param.idx_stop = N;
400 single_thread_param.m_k =
m_k;
401 single_thread_param.neighborhood_matrix = neighborhood_matrix.
matrix;
402 single_thread_param.z_matrix = z_matrix;
403 single_thread_param.feature_matrix = feature_matrix.
matrix;
404 single_thread_param.covariance_matrix = covariance_matrix;
405 single_thread_param.id_vector = id_vector;
406 single_thread_param.W_matrix = W_matrix;
425 int eigenproblem_status = 0;
433 SG_ERROR(
"ARPACK is not supported in this configuration.\n");
438 arpack_dsxupd(matrix.
matrix,NULL,
false,N,dimension+1,
"LA",
true,3,
true,
false,
m_nullspace_shift,0.0,
439 eigenvalues_vector,matrix.
matrix,eigenproblem_status);
442 if (eigenproblem_status)
443 SG_ERROR(
"ARPACK failed with code: %d", eigenproblem_status);
445 for (i=0; i<dimension; i++)
448 nullspace_features[j*dimension+i] = matrix[j*(dimension+1)+i+1];
457 wrap_dsyevr(
'V',
'U',N,matrix.
matrix,N,2,dimension+2,eigenvalues_vector,eigenvectors,&eigenproblem_status);
458 if (eigenproblem_status)
459 SG_ERROR(
"LAPACK failed with code: %d", eigenproblem_status);
462 for (i=0; i<dimension; i++)
465 nullspace_features[j*dimension+i] = eigenvectors[i*N+j];
475 LINRECONSTRUCTION_THREAD_PARAM* parameters = (LINRECONSTRUCTION_THREAD_PARAM*)p;
476 int32_t idx_start = parameters->idx_start;
477 int32_t idx_step = parameters->idx_step;
478 int32_t idx_stop = parameters->idx_stop;
479 int32_t
m_k = parameters->m_k;
480 const int32_t* neighborhood_matrix = parameters->neighborhood_matrix;
481 float64_t* z_matrix = parameters->z_matrix;
482 const float64_t* feature_matrix = parameters->feature_matrix;
483 float64_t* covariance_matrix = parameters->covariance_matrix;
484 float64_t* id_vector = parameters->id_vector;
485 float64_t* W_matrix = parameters->W_matrix;
488 PTHREAD_LOCK_T* W_matrix_lock = parameters->W_matrix_lock;
492 int32_t dim = parameters->dim;
493 int32_t N = parameters->N;
494 int32_t actual_k = parameters->actual_k;
497 for (i=idx_start; i<idx_stop; i+=idx_step)
501 for (j=0; j<
m_k; j++)
503 cblas_dcopy(dim,feature_matrix+neighborhood_matrix[i*actual_k+j]*dim,1,z_matrix+j*dim,1);
504 cblas_daxpy(dim,-1.0,feature_matrix+i*dim,1,z_matrix+j*dim,1);
508 cblas_dgemm(CblasColMajor,CblasTrans,CblasNoTrans,
512 0.0,covariance_matrix,m_k);
514 for (j=0; j<
m_k; j++)
522 for (j=0; j<
m_k; j++)
523 trace += covariance_matrix[j*m_k+j];
525 for (j=0; j<
m_k; j++)
526 covariance_matrix[j*m_k+j] += m_reconstruction_shift*trace;
531 clapack_dposv(CblasColMajor,CblasLower,m_k,1,covariance_matrix,m_k,id_vector,m_k);
535 for (j=0; j<
m_k; j++)
536 norming += id_vector[j];
538 cblas_dscal(m_k,1.0/norming,id_vector,1);
540 memset(covariance_matrix,0,
sizeof(
float64_t)*m_k*m_k);
541 cblas_dger(CblasColMajor,m_k,m_k,1.0,id_vector,1,id_vector,1,covariance_matrix,m_k);
544 W_matrix[N*i+i] += 1.0;
546 PTHREAD_LOCK(W_matrix_lock);
548 for (j=0; j<
m_k; j++)
550 W_matrix[N*i+neighborhood_matrix[i*actual_k+j]] -= id_vector[j];
551 W_matrix[N*neighborhood_matrix[i*actual_k+j]+i] -= id_vector[j];
553 for (j=0; j<
m_k; j++)
555 for (k=0; k<
m_k; k++)
556 W_matrix[N*neighborhood_matrix[i*actual_k+j]+neighborhood_matrix[i*actual_k+k]]+=
557 covariance_matrix[j*m_k+k];
560 PTHREAD_UNLOCK(W_matrix_lock);
569 int32_t N = distance_matrix.
num_rows;
571 int32_t* neighborhood_matrix =
SG_MALLOC(int32_t, N*k);
578 coverTree->
insert(LLE_COVERTREE_POINT(i,distance_matrix));
582 std::vector<LLE_COVERTREE_POINT> neighbors =
584 for (std::size_t m=1; m<unsigned(k+1); m++)
585 neighborhood_matrix[i*k+m-1] = neighbors[m].point_index;