26 using namespace shogun;
28 #ifndef DOXYGEN_SHOULD_SKIP_THIS
29 struct TRIANGULATION_THREAD_PARAM
54 const int32_t* lmk_idxs;
56 const bool* to_process;
89 SG_ERROR(
"Number of landmarks should be greater than 3 to make triangulation possible while %d given.",
111 return "MultidimensionalScaling";
133 return distance_matrix;
153 int32_t N = distance_matrix.
num_cols;
165 distance_matrix[i*N+j] = dsq;
166 distance_matrix[j*N+i] = dsq;
172 distance_matrix[i*N+i] *= -0.5;
173 for (j=i+1; j<N; j++)
175 distance_matrix[i*N+j] *= -0.5;
176 distance_matrix[j*N+i] *= -0.5;
184 int eigenproblem_status = 0;
189 arpack_dsxupd(distance_matrix.
matrix,NULL,
false,N,m_target_dim,
"LM",
false,1,
false,0.0,0.0,
190 eigenvalues_vector,replace_feature_matrix,eigenproblem_status);
192 ASSERT(eigenproblem_status == 0);
197 for (i=0; i<m_target_dim/2; i++)
199 tmp = replace_feature_matrix[j*m_target_dim+i];
200 replace_feature_matrix[j*m_target_dim+i] =
201 replace_feature_matrix[j*m_target_dim+(m_target_dim-i-1)];
202 replace_feature_matrix[j*m_target_dim+(m_target_dim-i-1)] = tmp;
206 for (i=0; i<m_target_dim/2; i++)
208 tmp = eigenvalues_vector[i];
209 eigenvalues_vector[i] = eigenvalues_vector[m_target_dim-i-1];
210 eigenvalues_vector[m_target_dim-i-1] = tmp;
217 replace_feature_matrix[j*m_target_dim+i] *=
229 wrap_dsyevr(
'V',
'U',N,distance_matrix.
matrix,N,N-m_target_dim+1,N,eigenvalues_vector,eigenvectors,&eigenproblem_status);
231 ASSERT(eigenproblem_status==0);
249 replace_feature_matrix[j*m_target_dim+i] =
261 SG_WARNING(
"Embedding is not consistent (got neg eigenvalues): features %d-%d are wrong",
276 int32_t total_N = distance_matrix.
num_cols;
279 SG_ERROR(
"Number of landmarks (%d) should be greater than 3 for proper triangulation.\n",
284 SG_ERROR(
"Number of landmarks (%d) should be less than total number of vectors (%d).\n",
292 for (i=0; i<lmk_N; i++)
294 for (j=0; j<lmk_N; j++)
295 lmk_dist_matrix[i*lmk_N+j] =
296 distance_matrix[lmk_idxs.
vector[i]*total_N+lmk_idxs.
vector[j]];
303 for (i=0; i<lmk_N; i++)
305 for (j=0; j<lmk_N; j++)
306 mean_sq_dist_vector[i] +=
CMath::sq(lmk_dist_matrix[i*lmk_N+j]);
308 mean_sq_dist_vector[i] /= lmk_N;
318 for (i=0; i<lmk_N; i++)
321 new_feature_matrix[lmk_idxs.
vector[i]*m_target_dim+j] =
322 lmk_feature_matrix[i*m_target_dim+j];
327 for (i=0; i<lmk_N; i++)
335 bool* to_process =
SG_MALLOC(
bool, total_N);
336 for (j=0; j<total_N; j++)
337 to_process[j] =
true;
338 for (j=0; j<lmk_N; j++)
339 to_process[lmk_idxs.
vector[j]] =
false;
346 pthread_t* threads =
SG_MALLOC(pthread_t, num_threads);
347 TRIANGULATION_THREAD_PARAM* parameters =
SG_MALLOC(TRIANGULATION_THREAD_PARAM, num_threads);
349 pthread_attr_init(&attr);
350 pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
353 for (t=0; t<num_threads; t++)
355 parameters[t].idx_start = t;
356 parameters[t].idx_stop = total_N;
357 parameters[t].idx_step = num_threads;
358 parameters[t].lmk_N = lmk_N;
359 parameters[t].total_N = total_N;
361 parameters[t].current_dist_to_lmks = current_dist_to_lmks+t*lmk_N;
362 parameters[t].distance_matrix = distance_matrix.
matrix;
363 parameters[t].mean_sq_dist_vector = mean_sq_dist_vector;
364 parameters[t].lmk_idxs = lmk_idxs.
vector;
365 parameters[t].lmk_feature_matrix = lmk_feature_matrix.
matrix;
366 parameters[t].new_feature_matrix = new_feature_matrix;
367 parameters[t].to_process = to_process;
371 for (t=0; t<num_threads; t++)
372 pthread_join(threads[t], NULL);
373 pthread_attr_destroy(&attr);
379 TRIANGULATION_THREAD_PARAM single_thread_param;
380 single_thread_param.idx_start = 0;
381 single_thread_param.idx_stop = total_N;
382 single_thread_param.idx_step = 1;
383 single_thread_param.lmk_N = lmk_N;
384 single_thread_param.total_N = total_N;
386 single_thread_param.current_dist_to_lmks = current_dist_to_lmks;
387 single_thread_param.distance_matrix = distance_matrix.
matrix;
388 single_thread_param.mean_sq_dist_vector = mean_sq_dist_vector;
389 single_thread_param.lmk_idxs = lmk_idxs.
vector;
390 single_thread_param.lmk_feature_matrix = lmk_feature_matrix.
matrix;
391 single_thread_param.new_feature_matrix = new_feature_matrix;
392 single_thread_param.to_process = to_process;
408 TRIANGULATION_THREAD_PARAM* parameters = (TRIANGULATION_THREAD_PARAM*)p;
409 int32_t idx_start = parameters->idx_start;
410 int32_t idx_step = parameters->idx_step;
411 int32_t idx_stop = parameters->idx_stop;
412 const int32_t* lmk_idxs = parameters->lmk_idxs;
413 const float64_t* distance_matrix = parameters->distance_matrix;
414 const float64_t* mean_sq_dist_vector = parameters->mean_sq_dist_vector;
415 float64_t* current_dist_to_lmks = parameters->current_dist_to_lmks;
417 int32_t lmk_N = parameters->lmk_N;
418 int32_t total_N = parameters->total_N;
419 const bool* to_process = parameters->to_process;
420 float64_t* lmk_feature_matrix = parameters->lmk_feature_matrix;
421 float64_t* new_feature_matrix = parameters->new_feature_matrix;
424 for (i=idx_start; i<idx_stop; i+=idx_step)
431 for (k=0; k<lmk_N; k++)
433 current_dist_to_lmks[k] =
434 CMath::sq(distance_matrix[i*total_N+lmk_idxs[k]]) -
435 mean_sq_dist_vector[k];
438 cblas_dgemv(CblasColMajor,CblasNoTrans,
440 -0.5,lmk_feature_matrix,m_target_dim,
441 current_dist_to_lmks,1,
442 0.0,(new_feature_matrix+i*m_target_dim),1);
450 int32_t* idxs =
SG_MALLOC(int32_t, total_count);
452 int32_t* permuted_idxs =
SG_MALLOC(int32_t, count);
455 for (i=0; i<total_count; i++)
457 for (i=0; i<count; i++)
458 permuted_idxs[i] = idxs[i];
459 for (i=count; i<total_count; i++)
463 permuted_idxs[rnd] = idxs[i];