1 // Copyright (C) 2012 Emanuele Cesena (emanuele.cesena@gmail.com), Davis E. King 2 // License: Boost Software License See LICENSE.txt for the full license. 3 #ifndef DLIB_SAMMoN_Hh_ 4 #define DLIB_SAMMoN_Hh_ 5 6 #include "sammon_abstract.h" 7 #include "../matrix.h" 8 #include "../algs.h" 9 #include "dpca.h" 10 #include <vector> 11 12 namespace dlib 13 { 14 15 class sammon_projection 16 { 17 18 public: 19 20 // ------------------------------------------------------------------------------------ 21 22 template <typename matrix_type> operator()23 std::vector<matrix<double,0,1> > operator() ( 24 const std::vector<matrix_type>& data, 25 const long num_dims 26 ) 27 { 28 // make sure requires clause is not broken 29 DLIB_ASSERT(num_dims > 0, 30 "\t std::vector<matrix<double,0,1> > sammon_projection::operator()" 31 << "\n\t Invalid inputs were given to this function." 32 << "\n\t num_dims: " << num_dims 33 ); 34 std::vector<matrix<double,0,1> > result; // projections 35 if (data.size() == 0) 36 { 37 return result; 38 } 39 40 #ifdef ENABLE_ASSERTS 41 DLIB_ASSERT(0 < num_dims && num_dims <= data[0].size(), 42 "\t std::vector<matrix<double,0,1> > sammon_projection::operator()" 43 << "\n\t Invalid inputs were given to this function." 44 << "\n\t data.size(): " << data.size() 45 << "\n\t num_dims: " << num_dims 46 << "\n\t data[0].size(): " << data[0].size() 47 ); 48 for (unsigned long i = 0; i < data.size(); ++i) 49 { 50 DLIB_ASSERT(is_col_vector(data[i]) && data[i].size() == data[0].size(), 51 "\t std::vector<matrix<double,0,1> > sammon_projection::operator()" 52 << "\n\t Invalid inputs were given to this function." 53 << "\n\t data["<<i<<"].size(): " << data[i].size() 54 << "\n\t data[0].size(): " << data[0].size() 55 << "\n\t is_col_vector(data["<<i<<"]): " << is_col_vector(data[i]) 56 ); 57 } 58 #endif 59 60 double err; // error (discarded) 61 do_sammon_projection(data, num_dims, result, err); 62 return result; 63 } 64 65 // ------------------------------------------------------------------------------------ 66 67 template <typename matrix_type> operator()68 void operator() ( 69 const std::vector<matrix_type>& data, 70 const long num_dims, 71 std::vector<matrix<double,0,1> >& result, 72 double &err, 73 const unsigned long num_iters = 1000, 74 const double err_delta = 1.0e-9 75 ) 76 { 77 // make sure requires clause is not broken 78 DLIB_ASSERT(num_dims > 0 && num_iters > 0 && err_delta > 0.0, 79 "\t std::vector<matrix<double,0,1> > sammon_projection::operator()" 80 << "\n\t Invalid inputs were given to this function." 81 << "\n\t data.size(): " << data.size() 82 << "\n\t num_dims: " << num_dims 83 << "\n\t num_iters: " << num_iters 84 << "\n\t err_delta: " << err_delta 85 ); 86 if (data.size() == 0) 87 { 88 result.clear(); 89 err = 0; 90 return; 91 } 92 93 #ifdef ENABLE_ASSERTS 94 DLIB_ASSERT(0 < num_dims && num_dims <= data[0].size(), 95 "\t std::vector<matrix<double,0,1> > sammon_projection::operator()" 96 << "\n\t Invalid inputs were given to this function." 97 << "\n\t data.size(): " << data.size() 98 << "\n\t num_dims: " << num_dims 99 << "\n\t data[0].size(): " << data[0].size() 100 ); 101 for (unsigned long i = 0; i < data.size(); ++i) 102 { 103 DLIB_ASSERT(is_col_vector(data[i]) && data[i].size() == data[0].size(), 104 "\t std::vector<matrix<double,0,1> > sammon_projection::operator()" 105 << "\n\t Invalid inputs were given to this function." 106 << "\n\t data["<<i<<"].size(): " << data[i].size() 107 << "\n\t data[0].size(): " << data[0].size() 108 << "\n\t is_col_vector(data["<<i<<"]): " << is_col_vector(data[i]) 109 ); 110 } 111 #endif 112 113 do_sammon_projection(data, num_dims, result, err, num_iters, err_delta); 114 } 115 116 // ---------------------------------------------------------------------------------------- 117 // ---------------------------------------------------------------------------------------- 118 119 private: 120 121 void compute_relative_distances( 122 matrix<double,0,1>& dist, // relative distances (output) 123 matrix<double,0,0>& data, // input data (matrix whose columns are the input vectors) 124 double eps_ratio = 1.0e-7 // to compute the minimum distance eps 125 ) 126 /*! 127 requires 128 - dist.nc() == comb( data.nc(), 2 ), preallocated 129 - eps_ratio > 0 130 ensures 131 - dist[k] == lenght(data[i] - data[j]) for k = j(j-1)/2 + i 132 !*/ 133 { 134 const long N = data.nc(); // num of points 135 double eps; // minimum distance, forced to avoid vectors collision 136 // computed at runtime as eps_ration * mean(vectors distances) 137 for (int k = 0, i = 1; i < N; ++i) 138 for (int j = 0; j < i; ++j) 139 dist(k++) = length(colm(data, i) - colm(data, j)); 140 141 eps = eps_ratio * mean(dist); 142 dist = lowerbound(dist, eps); 143 } 144 145 // ---------------------------------------------------------------------------------------- 146 147 template <typename matrix_type> 148 void do_sammon_projection( 149 const std::vector<matrix_type>& data, // input data 150 unsigned long num_dims, // dimension of the reduced space 151 std::vector<matrix<double,0,1> >& result, // projections (output) 152 double &err, // error (output) 153 unsigned long num_iters = 1000, // max num of iterations: stop condition 154 const double err_delta = 1.0e-9 // delta error: stop condition 155 ) 156 /*! 157 requires 158 - matrix_type should be a kind of dlib::matrix<double,N,1> 159 - num_dims > 0 160 - num_iters > 0 161 - err_delta > 0 162 ensures 163 - result == a set of matrix<double,num_dims,1> objects that represent 164 the Sammon's projections of data vectors. 165 - err == the estimated error done in the projection, with the extra 166 property that err(at previous iteration) - err < err_delta 167 !*/ 168 { 169 // other params 170 const double mf = 0.3; // magic factor 171 172 matrix<double> mdata; // input data as matrix 173 matrix<double> projs; // projected vectors, i.e. output data as matrix 174 175 // std::vector<matrix> -> matrix 176 mdata.set_size(data[0].size(), data.size()); 177 for (unsigned int i = 0; i < data.size(); i++) 178 set_colm(mdata, i) = data[i]; 179 180 const long N = mdata.nc(); // num of points 181 const long d = num_dims; // size of the reduced space 182 const long nd = N * (N - 1) / 2; // num of pairs of points = size of the distances vectors 183 184 matrix<double, 0, 1> dsij, inv_dsij; // d*_ij: pair-wise distances in the input space (and inverses) 185 dsij.set_size(nd, 1); 186 inv_dsij.set_size(nd, 1); 187 double ic; // 1.0 / sum of dsij 188 189 matrix<double, 0, 1> dij; // d_ij: pair-wise distances in the reduced space 190 dij.set_size(nd, 1); 191 192 matrix<double, 0, 0> dE, dE2, dtemp; // matrices representing error partial derivatives 193 dE.set_size(d, N); 194 dE2.set_size(d, N); 195 dtemp.set_size(d, N); 196 197 matrix<double, 0, 1> inv_dij, alpha; // utility vectors used to compute the partial derivatives 198 inv_dij.set_size(N, 1); // inv_dij is 1.0/dij, but we only need it column-wise 199 alpha.set_size(N, 1); // (slightly wasting a bit of computation) 200 // alpha = 1.0/dij - 1.0/dsij, again column-wise 201 202 203 // initialize projs with PCA 204 discriminant_pca<matrix<double> > dpca; 205 for (int i = 0; i < mdata.nc(); ++i) 206 { 207 dpca.add_to_total_variance(colm(mdata, i)); 208 } 209 matrix<double> mat = dpca.dpca_matrix_of_size(num_dims); 210 projs = mat * mdata; 211 212 // compute dsij, inv_dsij and ic 213 compute_relative_distances(dsij, mdata); 214 inv_dsij = 1.0 / dsij; 215 ic = 1.0 / sum(dsij); 216 217 // compute dij and err 218 compute_relative_distances(dij, projs); 219 err = ic * sum(pointwise_multiply(squared(dij - dsij), inv_dsij)); 220 221 // start iterating 222 while (num_iters--) 223 { 224 // compute dE, dE2 progressively column by column 225 for (int p = 0; p < N; ++p) 226 { 227 // compute 228 // - alpha_p, the column vector with 1/d_pj - 1/d*_pj 229 // - dtemp, the matrix with the p-th column repeated all along 230 //TODO: optimize constructions 231 for (int i = 0; i < N; ++i) 232 { 233 int pos = (i < p) ? p * (p - 1) / 2 + i : i * (i - 1) / 2 + p; 234 inv_dij(i) = (i == p) ? 0.0 : 1.0 / dij(pos); 235 alpha(i) = (i == p) ? 0.0 : inv_dij(i) - inv_dsij(pos); 236 set_colm(dtemp, i) = colm(projs, p); 237 } 238 239 dtemp -= projs; 240 set_colm(dE, p) = dtemp * alpha; 241 242 double sum_alpha = sum(alpha); 243 set_colm(dE2, p) = abs( sum_alpha + squared(dtemp) * cubed(inv_dij) ); 244 } 245 246 247 // compute the update projections 248 projs += pointwise_multiply(dE, mf * reciprocal(dE2)); 249 250 // compute new dij and error 251 compute_relative_distances(dij, projs); 252 double err_new = ic * sum( pointwise_multiply(squared(dij - dsij), inv_dsij) ); 253 if (err - err_new < err_delta) 254 break; 255 err = err_new; 256 } 257 258 // matrix -> std::vector<matrix> 259 result.clear(); 260 for (int i = 0; i < projs.nc(); ++i) 261 result.push_back(colm(projs, i)); 262 } 263 264 }; 265 266 } // namespace dlib 267 268 #endif // DLIB_SAMMoN_Hh_ 269 270