269 lines
11 KiB
C++
269 lines
11 KiB
C++
// Copyright (C) 2012 Emanuele Cesena (emanuele.cesena@gmail.com), Davis E. King
|
|
// License: Boost Software License See LICENSE.txt for the full license.
|
|
#ifndef DLIB_SAMMoN_Hh_
|
|
#define DLIB_SAMMoN_Hh_
|
|
|
|
#include "sammon_abstract.h"
|
|
#include "../matrix.h"
|
|
#include "../algs.h"
|
|
#include "dpca.h"
|
|
#include <vector>
|
|
|
|
namespace dlib
|
|
{
|
|
|
|
class sammon_projection
|
|
{
|
|
|
|
public:
|
|
|
|
// ------------------------------------------------------------------------------------
|
|
|
|
template <typename matrix_type>
|
|
std::vector<matrix<double,0,1> > operator() (
|
|
const std::vector<matrix_type>& data,
|
|
const long num_dims
|
|
)
|
|
{
|
|
// make sure requires clause is not broken
|
|
DLIB_ASSERT(num_dims > 0,
|
|
"\t std::vector<matrix<double,0,1> > sammon_projection::operator()"
|
|
<< "\n\t Invalid inputs were given to this function."
|
|
<< "\n\t num_dims: " << num_dims
|
|
);
|
|
std::vector<matrix<double,0,1> > result; // projections
|
|
if (data.size() == 0)
|
|
{
|
|
return result;
|
|
}
|
|
|
|
#ifdef ENABLE_ASSERTS
|
|
DLIB_ASSERT(0 < num_dims && num_dims <= data[0].size(),
|
|
"\t std::vector<matrix<double,0,1> > sammon_projection::operator()"
|
|
<< "\n\t Invalid inputs were given to this function."
|
|
<< "\n\t data.size(): " << data.size()
|
|
<< "\n\t num_dims: " << num_dims
|
|
<< "\n\t data[0].size(): " << data[0].size()
|
|
);
|
|
for (unsigned long i = 0; i < data.size(); ++i)
|
|
{
|
|
DLIB_ASSERT(is_col_vector(data[i]) && data[i].size() == data[0].size(),
|
|
"\t std::vector<matrix<double,0,1> > sammon_projection::operator()"
|
|
<< "\n\t Invalid inputs were given to this function."
|
|
<< "\n\t data["<<i<<"].size(): " << data[i].size()
|
|
<< "\n\t data[0].size(): " << data[0].size()
|
|
<< "\n\t is_col_vector(data["<<i<<"]): " << is_col_vector(data[i])
|
|
);
|
|
}
|
|
#endif
|
|
|
|
double err; // error (discarded)
|
|
do_sammon_projection(data, num_dims, result, err);
|
|
return result;
|
|
}
|
|
|
|
// ------------------------------------------------------------------------------------
|
|
|
|
template <typename matrix_type>
|
|
void operator() (
|
|
const std::vector<matrix_type>& data,
|
|
const long num_dims,
|
|
std::vector<matrix<double,0,1> >& result,
|
|
double &err,
|
|
const unsigned long num_iters = 1000,
|
|
const double err_delta = 1.0e-9
|
|
)
|
|
{
|
|
// make sure requires clause is not broken
|
|
DLIB_ASSERT(num_dims > 0 && num_iters > 0 && err_delta > 0.0,
|
|
"\t std::vector<matrix<double,0,1> > sammon_projection::operator()"
|
|
<< "\n\t Invalid inputs were given to this function."
|
|
<< "\n\t data.size(): " << data.size()
|
|
<< "\n\t num_dims: " << num_dims
|
|
<< "\n\t num_iters: " << num_iters
|
|
<< "\n\t err_delta: " << err_delta
|
|
);
|
|
if (data.size() == 0)
|
|
{
|
|
result.clear();
|
|
err = 0;
|
|
return;
|
|
}
|
|
|
|
#ifdef ENABLE_ASSERTS
|
|
DLIB_ASSERT(0 < num_dims && num_dims <= data[0].size(),
|
|
"\t std::vector<matrix<double,0,1> > sammon_projection::operator()"
|
|
<< "\n\t Invalid inputs were given to this function."
|
|
<< "\n\t data.size(): " << data.size()
|
|
<< "\n\t num_dims: " << num_dims
|
|
<< "\n\t data[0].size(): " << data[0].size()
|
|
);
|
|
for (unsigned long i = 0; i < data.size(); ++i)
|
|
{
|
|
DLIB_ASSERT(is_col_vector(data[i]) && data[i].size() == data[0].size(),
|
|
"\t std::vector<matrix<double,0,1> > sammon_projection::operator()"
|
|
<< "\n\t Invalid inputs were given to this function."
|
|
<< "\n\t data["<<i<<"].size(): " << data[i].size()
|
|
<< "\n\t data[0].size(): " << data[0].size()
|
|
<< "\n\t is_col_vector(data["<<i<<"]): " << is_col_vector(data[i])
|
|
);
|
|
}
|
|
#endif
|
|
|
|
do_sammon_projection(data, num_dims, result, err, num_iters, err_delta);
|
|
}
|
|
|
|
// ----------------------------------------------------------------------------------------
|
|
// ----------------------------------------------------------------------------------------
|
|
|
|
private:
|
|
|
|
void compute_relative_distances(
|
|
matrix<double,0,1>& dist, // relative distances (output)
|
|
matrix<double,0,0>& data, // input data (matrix whose columns are the input vectors)
|
|
double eps_ratio = 1.0e-7 // to compute the minimum distance eps
|
|
)
|
|
/*!
|
|
requires
|
|
- dist.nc() == comb( data.nc(), 2 ), preallocated
|
|
- eps_ratio > 0
|
|
ensures
|
|
- dist[k] == lenght(data[i] - data[j]) for k = j(j-1)/2 + i
|
|
!*/
|
|
{
|
|
const long N = data.nc(); // num of points
|
|
double eps; // minimum distance, forced to avoid vectors collision
|
|
// computed at runtime as eps_ration * mean(vectors distances)
|
|
for (int k = 0, i = 1; i < N; ++i)
|
|
for (int j = 0; j < i; ++j)
|
|
dist(k++) = length(colm(data, i) - colm(data, j));
|
|
|
|
eps = eps_ratio * mean(dist);
|
|
dist = lowerbound(dist, eps);
|
|
}
|
|
|
|
// ----------------------------------------------------------------------------------------
|
|
|
|
template <typename matrix_type>
|
|
void do_sammon_projection(
|
|
const std::vector<matrix_type>& data, // input data
|
|
unsigned long num_dims, // dimension of the reduced space
|
|
std::vector<matrix<double,0,1> >& result, // projections (output)
|
|
double &err, // error (output)
|
|
unsigned long num_iters = 1000, // max num of iterations: stop condition
|
|
const double err_delta = 1.0e-9 // delta error: stop condition
|
|
)
|
|
/*!
|
|
requires
|
|
- matrix_type should be a kind of dlib::matrix<double,N,1>
|
|
- num_dims > 0
|
|
- num_iters > 0
|
|
- err_delta > 0
|
|
ensures
|
|
- result == a set of matrix<double,num_dims,1> objects that represent
|
|
the Sammon's projections of data vectors.
|
|
- err == the estimated error done in the projection, with the extra
|
|
property that err(at previous iteration) - err < err_delta
|
|
!*/
|
|
{
|
|
// other params
|
|
const double mf = 0.3; // magic factor
|
|
|
|
matrix<double> mdata; // input data as matrix
|
|
matrix<double> projs; // projected vectors, i.e. output data as matrix
|
|
|
|
// std::vector<matrix> -> matrix
|
|
mdata.set_size(data[0].size(), data.size());
|
|
for (unsigned int i = 0; i < data.size(); i++)
|
|
set_colm(mdata, i) = data[i];
|
|
|
|
const long N = mdata.nc(); // num of points
|
|
const long d = num_dims; // size of the reduced space
|
|
const long nd = N * (N - 1) / 2; // num of pairs of points = size of the distances vectors
|
|
|
|
matrix<double, 0, 1> dsij, inv_dsij; // d*_ij: pair-wise distances in the input space (and inverses)
|
|
dsij.set_size(nd, 1);
|
|
inv_dsij.set_size(nd, 1);
|
|
double ic; // 1.0 / sum of dsij
|
|
|
|
matrix<double, 0, 1> dij; // d_ij: pair-wise distances in the reduced space
|
|
dij.set_size(nd, 1);
|
|
|
|
matrix<double, 0, 0> dE, dE2, dtemp; // matrices representing error partial derivatives
|
|
dE.set_size(d, N);
|
|
dE2.set_size(d, N);
|
|
dtemp.set_size(d, N);
|
|
|
|
matrix<double, 0, 1> inv_dij, alpha; // utility vectors used to compute the partial derivatives
|
|
inv_dij.set_size(N, 1); // inv_dij is 1.0/dij, but we only need it column-wise
|
|
alpha.set_size(N, 1); // (slightly wasting a bit of computation)
|
|
// alpha = 1.0/dij - 1.0/dsij, again column-wise
|
|
|
|
|
|
// initialize projs with PCA
|
|
discriminant_pca<matrix<double> > dpca;
|
|
for (int i = 0; i < mdata.nc(); ++i)
|
|
{
|
|
dpca.add_to_total_variance(colm(mdata, i));
|
|
}
|
|
matrix<double> mat = dpca.dpca_matrix_of_size(num_dims);
|
|
projs = mat * mdata;
|
|
|
|
// compute dsij, inv_dsij and ic
|
|
compute_relative_distances(dsij, mdata);
|
|
inv_dsij = 1.0 / dsij;
|
|
ic = 1.0 / sum(dsij);
|
|
|
|
// compute dij and err
|
|
compute_relative_distances(dij, projs);
|
|
err = ic * sum(pointwise_multiply(squared(dij - dsij), inv_dsij));
|
|
|
|
// start iterating
|
|
while (num_iters--)
|
|
{
|
|
// compute dE, dE2 progressively column by column
|
|
for (int p = 0; p < N; ++p)
|
|
{
|
|
// compute
|
|
// - alpha_p, the column vector with 1/d_pj - 1/d*_pj
|
|
// - dtemp, the matrix with the p-th column repeated all along
|
|
//TODO: optimize constructions
|
|
for (int i = 0; i < N; ++i)
|
|
{
|
|
int pos = (i < p) ? p * (p - 1) / 2 + i : i * (i - 1) / 2 + p;
|
|
inv_dij(i) = (i == p) ? 0.0 : 1.0 / dij(pos);
|
|
alpha(i) = (i == p) ? 0.0 : inv_dij(i) - inv_dsij(pos);
|
|
set_colm(dtemp, i) = colm(projs, p);
|
|
}
|
|
|
|
dtemp -= projs;
|
|
set_colm(dE, p) = dtemp * alpha;
|
|
|
|
double sum_alpha = sum(alpha);
|
|
set_colm(dE2, p) = abs( sum_alpha + squared(dtemp) * cubed(inv_dij) );
|
|
}
|
|
|
|
|
|
// compute the update projections
|
|
projs += pointwise_multiply(dE, mf * reciprocal(dE2));
|
|
|
|
// compute new dij and error
|
|
compute_relative_distances(dij, projs);
|
|
double err_new = ic * sum( pointwise_multiply(squared(dij - dsij), inv_dsij) );
|
|
if (err - err_new < err_delta)
|
|
break;
|
|
err = err_new;
|
|
}
|
|
|
|
// matrix -> std::vector<matrix>
|
|
result.clear();
|
|
for (int i = 0; i < projs.nc(); ++i)
|
|
result.push_back(colm(projs, i));
|
|
}
|
|
|
|
};
|
|
|
|
} // namespace dlib
|
|
|
|
#endif // DLIB_SAMMoN_Hh_
|
|
|