Tapkee
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
multidimensional_scaling.hpp
Go to the documentation of this file.
1 /* This software is distributed under BSD 3-clause license (see LICENSE file).
2  *
3  * Copyright (c) 2012-2013 Sergey Lisitsyn
4  */
5 
6 #ifndef TAPKEE_MDS_H_
7 #define TAPKEE_MDS_H_
8 
9 /* Tapkee includes */
10 #include <tapkee/defines.hpp>
11 #include <tapkee/utils/time.hpp>
12 /* End of Tapkee includes */
13 
14 namespace tapkee
15 {
16 namespace tapkee_internal
17 {
18 
19 template <class RandomAccessIterator>
20 Landmarks select_landmarks_random(RandomAccessIterator begin, RandomAccessIterator end, ScalarType ratio)
21 {
22  Landmarks landmarks;
23  landmarks.reserve(end-begin);
24  for (RandomAccessIterator iter=begin; iter!=end; ++iter)
25  landmarks.push_back(iter-begin);
26  tapkee::random_shuffle(landmarks.begin(),landmarks.end());
27  landmarks.erase(landmarks.begin() + static_cast<IndexType>(landmarks.size()*ratio),landmarks.end());
28  return landmarks;
29 }
30 
31 template <class RandomAccessIterator, class PairwiseCallback>
32 DenseSymmetricMatrix compute_distance_matrix(RandomAccessIterator begin, RandomAccessIterator /*end*/,
33  const Landmarks& landmarks, PairwiseCallback callback)
34 {
35  timed_context context("Multidimensional scaling distance matrix computation");
36 
37  const IndexType n_landmarks = landmarks.size();
38  DenseSymmetricMatrix distance_matrix(n_landmarks,n_landmarks);
39 
40 #pragma omp parallel shared(begin,landmarks,distance_matrix,callback) default(none)
41  {
42  IndexType i_index_iter,j_index_iter;
43 #pragma omp for nowait
44  for (i_index_iter=0; i_index_iter<n_landmarks; ++i_index_iter)
45  {
46  for (j_index_iter=i_index_iter; j_index_iter<n_landmarks; ++j_index_iter)
47  {
48  ScalarType d = callback.distance(begin[landmarks[i_index_iter]],begin[landmarks[j_index_iter]]);
49  d *= d;
50  distance_matrix(i_index_iter,j_index_iter) = d;
51  distance_matrix(j_index_iter,i_index_iter) = d;
52  }
53  }
54  }
55  return distance_matrix;
56 }
57 
58 template <class RandomAccessIterator, class PairwiseCallback>
59 DenseMatrix triangulate(RandomAccessIterator begin, RandomAccessIterator end, PairwiseCallback distance_callback,
60  const Landmarks& landmarks, const DenseVector& landmark_distances_squared,
61  EigendecompositionResult& landmarks_embedding, IndexType target_dimension)
62 {
63  timed_context context("Landmark triangulation");
64 
65  const IndexType n_vectors = end-begin;
66  const IndexType n_landmarks = landmarks.size();
67 
68  bool* to_process = new bool[n_vectors];
69  std::fill(to_process,to_process+n_vectors,true);
70 
71  DenseMatrix embedding(n_vectors,target_dimension);
72 
73  for (IndexType index_iter=0; index_iter<n_landmarks; ++index_iter)
74  {
75  to_process[landmarks[index_iter]] = false;
76  embedding.row(landmarks[index_iter]).noalias() = landmarks_embedding.first.row(index_iter);
77  }
78 
79  for (IndexType i=0; i<target_dimension; ++i)
80  landmarks_embedding.first.col(i).array() /= landmarks_embedding.second(i);
81 
82 #pragma omp parallel shared(begin,end,to_process,distance_callback,landmarks, \
83  landmarks_embedding,landmark_distances_squared,embedding) default(none)
84  {
85  DenseVector distances_to_landmarks(n_landmarks);
86  IndexType index_iter;
87 #pragma omp for nowait
88  for (index_iter=0; index_iter<n_vectors; ++index_iter)
89  {
90  if (!to_process[index_iter])
91  continue;
92 
93  for (IndexType i=0; i<n_landmarks; ++i)
94  {
95  ScalarType d = distance_callback.distance(begin[index_iter],begin[landmarks[i]]);
96  distances_to_landmarks(i) = d*d;
97  }
98  //distances_to_landmarks.array().square();
99 
100  distances_to_landmarks -= landmark_distances_squared;
101  embedding.row(index_iter).noalias() = -0.5*landmarks_embedding.first.transpose()*distances_to_landmarks;
102  }
103  }
104 
105  delete[] to_process;
106 
107  return embedding;
108 }
109 
110 template <class RandomAccessIterator, class PairwiseCallback>
111 DenseSymmetricMatrix compute_distance_matrix(RandomAccessIterator begin, RandomAccessIterator end,
112  PairwiseCallback callback)
113 {
114  timed_context context("Multidimensional scaling distance matrix computation");
115 
116  const IndexType n_vectors = end-begin;
117  DenseSymmetricMatrix distance_matrix(n_vectors,n_vectors);
118 
119 #pragma omp parallel shared(begin,distance_matrix,callback) default(none)
120  {
121  IndexType i_index_iter,j_index_iter;
122 #pragma omp for nowait
123  for (i_index_iter=0; i_index_iter<n_vectors; ++i_index_iter)
124  {
125  for (j_index_iter=i_index_iter; j_index_iter<n_vectors; ++j_index_iter)
126  {
127  ScalarType d = callback.distance(begin[i_index_iter],begin[j_index_iter]);
128  d *= d;
129  distance_matrix(i_index_iter,j_index_iter) = d;
130  distance_matrix(j_index_iter,i_index_iter) = d;
131  }
132  }
133  }
134  return distance_matrix;
135 }
136 
137 } // End of namespace tapkee_internal
138 } // End of namespace tapkee
139 
140 #endif
void random_shuffle(RAI first, RAI last)
Definition: random.hpp:58
Eigen::Matrix< tapkee::ScalarType, Eigen::Dynamic, Eigen::Dynamic > DenseMatrix
dense matrix type (non-overridable)
Definition: types.hpp:23
DenseSymmetricMatrix compute_distance_matrix(RandomAccessIterator begin, RandomAccessIterator, const Landmarks &landmarks, PairwiseCallback callback)
double ScalarType
default scalar value (can be overrided with TAPKEE_CUSTOM_INTERNAL_NUMTYPE define) ...
Definition: types.hpp:15
Eigen::Matrix< tapkee::ScalarType, Eigen::Dynamic, 1 > DenseVector
dense vector type (non-overridable)
Definition: types.hpp:21
DenseMatrix triangulate(RandomAccessIterator begin, RandomAccessIterator end, PairwiseCallback distance_callback, const Landmarks &landmarks, const DenseVector &landmark_distances_squared, EigendecompositionResult &landmarks_embedding, IndexType target_dimension)
int IndexType
indexing type (non-overridable) set to int for compatibility with OpenMP 2.0
Definition: types.hpp:19
TAPKEE_INTERNAL_PAIR< tapkee::DenseMatrix, tapkee::DenseVector > EigendecompositionResult
Definition: synonyms.hpp:41
Landmarks select_landmarks_random(RandomAccessIterator begin, RandomAccessIterator end, ScalarType ratio)
TAPKEE_INTERNAL_VECTOR< tapkee::IndexType > Landmarks
Definition: synonyms.hpp:42
tapkee::DenseMatrix DenseSymmetricMatrix
dense symmetric matrix (non-overridable, currently just dense matrix, can be improved later) ...
Definition: types.hpp:25