```// Copyright (C) 2010  Davis E. King (davis@dlib.net)
#ifndef DLIB_SORT_BASIS_VECTORs_Hh_
#define DLIB_SORT_BASIS_VECTORs_Hh_

#include <vector>

#include "sort_basis_vectors_abstract.h"
#include "../matrix.h"
#include "../statistics.h"

namespace dlib
{

// ----------------------------------------------------------------------------------------

namespace bs_impl
{
template <typename EXP>
typename EXP::matrix_type invert (
const matrix_exp<EXP>& m
)
{
eigenvalue_decomposition<EXP> eig(make_symmetric(m));

typedef typename EXP::type scalar_type;
typedef typename EXP::mem_manager_type mm_type;

matrix<scalar_type,0,1,mm_type> vals = eig.get_real_eigenvalues();

const scalar_type max_eig = max(abs(vals));
const scalar_type thresh = max_eig*std::sqrt(std::numeric_limits<scalar_type>::epsilon());

// Since m might be singular or almost singular we need to do something about
// any very small eigenvalues.  So here we set the smallest eigenvalues to
// be equal to a large value to make the inversion stable.  We can't just set
// them to zero like in a normal pseudo-inverse since we want the resulting
// inverse matrix to be full rank.
for (long i = 0; i < vals.size(); ++i)
{
if (std::abs(vals(i)) < thresh)
vals(i) = max_eig;
}

// Build the inverse matrix.  This is basically a pseudo-inverse.
return make_symmetric(eig.get_pseudo_v()*diagm(reciprocal(vals))*trans(eig.get_pseudo_v()));
}

// ----------------------------------------------------------------------------------------

template <
typename kernel_type,
typename vect1_type,
typename vect2_type,
typename vect3_type
>
const std::vector<typename kernel_type::sample_type> sort_basis_vectors_impl (
const kernel_type& kern,
const vect1_type& samples,
const vect2_type& labels,
const vect3_type& basis,
double eps
)
{
DLIB_ASSERT(is_binary_classification_problem(samples, labels) &&
0 < eps && eps <= 1 &&
basis.size() > 0,
"\t void sort_basis_vectors()"
<< "\n\t Invalid arguments were given to this function."
<< "\n\t is_binary_classification_problem(samples, labels): " << is_binary_classification_problem(samples, labels)
<< "\n\t basis.size(): " << basis.size()
<< "\n\t eps:          " << eps
);

typedef typename kernel_type::scalar_type scalar_type;
typedef typename kernel_type::mem_manager_type mm_type;

typedef matrix<scalar_type,0,1,mm_type> col_matrix;
typedef matrix<scalar_type,0,0,mm_type> gen_matrix;

col_matrix c1_mean, c2_mean, temp, delta;

col_matrix weights;

running_covariance<gen_matrix> cov;

// compute the covariance matrix and the means of the two classes.
for (long i = 0; i < samples.size(); ++i)
{
temp = kernel_matrix(kern, basis, samples(i));
if (labels(i) > 0)
c1_mean += temp;
else
c2_mean += temp;
}

c1_mean /= sum(labels > 0);
c2_mean /= sum(labels < 0);

delta = c1_mean - c2_mean;

gen_matrix cov_inv = bs_impl::invert(cov.covariance());

matrix<long,0,1,mm_type> total_perm = trans(range(0, delta.size()-1));
matrix<long,0,1,mm_type> perm = total_perm;

std::vector<std::pair<scalar_type,long> > sorted_feats(delta.size());

long best_size = delta.size();
long misses = 0;
matrix<long,0,1,mm_type> best_total_perm = perm;

// Now we basically find fisher's linear discriminant over and over.  Each
// time sorting the features so that the most important ones pile up together.
weights = trans(chol(cov_inv))*delta;
while (true)
{

for (unsigned long i = 0; i < sorted_feats.size(); ++i)
sorted_feats[i] = make_pair(std::abs(weights(i)), i);

std::sort(sorted_feats.begin(), sorted_feats.end());

// make a permutation vector according to the sorting
for (long i = 0; i < perm.size(); ++i)
perm(i) = sorted_feats[i].second;

// Apply the permutation.  Doing this gives the same result as permuting all the
// features and then recomputing the delta and cov_inv from scratch.
cov_inv = subm(cov_inv,perm,perm);
delta = rowm(delta,perm);

// Record all the permutations we have done so we will know how the final
// weights match up with the original basis vectors when we are done.
total_perm = rowm(total_perm, perm);

// compute new Fisher weights for sorted features.
weights = trans(chol(cov_inv))*delta;

// Measure how many features it takes to account for eps% of the weights vector.
const scalar_type total_weight = length_squared(weights);
scalar_type weight_accum = 0;
long size = 0;
// figure out how to get eps% of the weights
for (long i = weights.size()-1; i >= 0; --i)
{
++size;
weight_accum += weights(i)*weights(i);
if (weight_accum/total_weight > eps)
break;
}

// loop until the best_size stops dropping
if (size < best_size)
{
misses = 0;
best_size = size;
best_total_perm = total_perm;
}
else
{
++misses;

// Give up once we have had 10 rounds where we didn't find a weights vector with
// a smaller concentration of good features.
if (misses >= 10)
break;
}

}

// make sure best_size isn't zero
if (best_size == 0)
best_size = 1;

std::vector<typename kernel_type::sample_type> sorted_basis;

// permute the basis so that it matches up with the contents of the best weights
sorted_basis.resize(best_size);
for (unsigned long i = 0; i < sorted_basis.size(); ++i)
{
// Note that we load sorted_basis backwards so that the most important
// basis elements come first.
sorted_basis[i] = basis(best_total_perm(basis.size()-i-1));
}

return sorted_basis;
}

}

// ----------------------------------------------------------------------------------------

template <
typename kernel_type,
typename vect1_type,
typename vect2_type,
typename vect3_type
>
const std::vector<typename kernel_type::sample_type> sort_basis_vectors (
const kernel_type& kern,
const vect1_type& samples,
const vect2_type& labels,
const vect3_type& basis,
double eps = 0.99
)
{
return bs_impl::sort_basis_vectors_impl(kern,
mat(samples),
mat(labels),
mat(basis),
eps);
}

// ----------------------------------------------------------------------------------------

}

#endif // DLIB_SORT_BASIS_VECTORs_Hh_

```