Tapkee
multidimensional_scaling.hpp
Go to the documentation of this file.
00001 /* This software is distributed under BSD 3-clause license (see LICENSE file).
00002  *
00003  * Copyright (c) 2012-2013 Sergey Lisitsyn
00004  */
00005 
00006 #ifndef TAPKEE_MDS_H_
00007 #define TAPKEE_MDS_H_
00008 
00009 /* Tapkee includes */
00010 #include <tapkee/defines.hpp>
00011 #include <tapkee/utils/time.hpp>
00012 /* End of Tapkee includes */
00013 
00014 namespace tapkee
00015 {
00016 namespace tapkee_internal
00017 {
00018 
00019 template <class RandomAccessIterator>
00020 Landmarks select_landmarks_random(RandomAccessIterator begin, RandomAccessIterator end, ScalarType ratio)
00021 {
00022     Landmarks landmarks;
00023     landmarks.reserve(end-begin);
00024     for (RandomAccessIterator iter=begin; iter!=end; ++iter)
00025         landmarks.push_back(iter-begin);
00026     tapkee::random_shuffle(landmarks.begin(),landmarks.end());
00027     landmarks.erase(landmarks.begin() + static_cast<IndexType>(landmarks.size()*ratio),landmarks.end());
00028     return landmarks;
00029 }
00030 
00031 template <class RandomAccessIterator, class PairwiseCallback>
00032 DenseSymmetricMatrix compute_distance_matrix(RandomAccessIterator begin, RandomAccessIterator /*end*/,
00033                                              const Landmarks& landmarks, PairwiseCallback callback)
00034 {
00035     timed_context context("Multidimensional scaling distance matrix computation");
00036 
00037     const IndexType n_landmarks = landmarks.size();
00038     DenseSymmetricMatrix distance_matrix(n_landmarks,n_landmarks);
00039 
00040 #pragma omp parallel shared(begin,landmarks,distance_matrix,callback) default(none)
00041     {
00042         IndexType i_index_iter,j_index_iter;
00043 #pragma omp for nowait
00044         for (i_index_iter=0; i_index_iter<n_landmarks; ++i_index_iter)
00045         {
00046             for (j_index_iter=i_index_iter; j_index_iter<n_landmarks; ++j_index_iter)
00047             {
00048                 ScalarType d = callback.distance(begin[landmarks[i_index_iter]],begin[landmarks[j_index_iter]]);
00049                 d *= d;
00050                 distance_matrix(i_index_iter,j_index_iter) = d;
00051                 distance_matrix(j_index_iter,i_index_iter) = d;
00052             }
00053         }
00054     }
00055     return distance_matrix;
00056 }
00057 
00058 template <class RandomAccessIterator, class PairwiseCallback>
00059 DenseMatrix triangulate(RandomAccessIterator begin, RandomAccessIterator end, PairwiseCallback distance_callback,
00060                         const Landmarks& landmarks, const DenseVector& landmark_distances_squared, 
00061                         EigendecompositionResult& landmarks_embedding, IndexType target_dimension)
00062 {
00063     timed_context context("Landmark triangulation");
00064     
00065     const IndexType n_vectors = end-begin;
00066     const IndexType n_landmarks = landmarks.size();
00067 
00068     bool* to_process = new bool[n_vectors];
00069     std::fill(to_process,to_process+n_vectors,true);
00070     
00071     DenseMatrix embedding(n_vectors,target_dimension);
00072 
00073     for (IndexType index_iter=0; index_iter<n_landmarks; ++index_iter)
00074     {
00075         to_process[landmarks[index_iter]] = false;
00076         embedding.row(landmarks[index_iter]).noalias() = landmarks_embedding.first.row(index_iter);
00077     }
00078 
00079     for (IndexType i=0; i<target_dimension; ++i)
00080         landmarks_embedding.first.col(i).array() /= landmarks_embedding.second(i);
00081 
00082 #pragma omp parallel shared(begin,end,to_process,distance_callback,landmarks, \
00083         landmarks_embedding,landmark_distances_squared,embedding) default(none)
00084     {
00085         DenseVector distances_to_landmarks(n_landmarks);
00086         IndexType index_iter;
00087 #pragma omp for nowait
00088         for (index_iter=0; index_iter<n_vectors; ++index_iter)
00089         {
00090             if (!to_process[index_iter])
00091                 continue;
00092 
00093             for (IndexType i=0; i<n_landmarks; ++i)
00094             {
00095                 ScalarType d = distance_callback.distance(begin[index_iter],begin[landmarks[i]]);
00096                 distances_to_landmarks(i) = d*d;
00097             }
00098             //distances_to_landmarks.array().square();
00099 
00100             distances_to_landmarks -= landmark_distances_squared;
00101             embedding.row(index_iter).noalias() = -0.5*landmarks_embedding.first.transpose()*distances_to_landmarks;
00102         }
00103     }
00104 
00105     delete[] to_process;
00106 
00107     return embedding;
00108 }
00109 
00110 template <class RandomAccessIterator, class PairwiseCallback>
00111 DenseSymmetricMatrix compute_distance_matrix(RandomAccessIterator begin, RandomAccessIterator end, 
00112                                              PairwiseCallback callback)
00113 {
00114     timed_context context("Multidimensional scaling distance matrix computation");
00115 
00116     const IndexType n_vectors = end-begin;
00117     DenseSymmetricMatrix distance_matrix(n_vectors,n_vectors);
00118 
00119 #pragma omp parallel shared(begin,distance_matrix,callback) default(none)
00120     {
00121         IndexType i_index_iter,j_index_iter;
00122 #pragma omp for nowait
00123         for (i_index_iter=0; i_index_iter<n_vectors; ++i_index_iter)
00124         {
00125             for (j_index_iter=i_index_iter; j_index_iter<n_vectors; ++j_index_iter)
00126             {
00127                 ScalarType d = callback.distance(begin[i_index_iter],begin[j_index_iter]);
00128                 d *= d;
00129                 distance_matrix(i_index_iter,j_index_iter) = d;
00130                 distance_matrix(j_index_iter,i_index_iter) = d;
00131             }
00132         }
00133     }
00134     return distance_matrix;
00135 }
00136 
00137 } // End of namespace tapkee_internal
00138 } // End of namespace tapkee
00139 
00140 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines