// Copyright (C) 2009 Davis E. King (davis@dlib.net)
// License: Boost Software License See LICENSE.txt for the full license.
#include <dlib/matrix.h>
#include <sstream>
#include <string>
#include <cstdlib>
#include <ctime>
#include <vector>
#include "../stl_checked.h"
#include "../array.h"
#include "../rand.h"
#include <dlib/string.h>
#include "tester.h"
namespace
{
using namespace test;
using namespace dlib;
using namespace std;
logger dlog("test.matrix_eig");
dlib::rand rnd;
// ----------------------------------------------------------------------------------------
template <typename type>
const matrix<type> randm(long r, long c)
{
matrix<type> m(r,c);
for (long row = 0; row < m.nr(); ++row)
{
for (long col = 0; col < m.nc(); ++col)
{
m(row,col) = static_cast<type>(rnd.get_random_double());
}
}
return m;
}
template <typename type, long NR, long NC>
const matrix<type,NR,NC> randm()
{
matrix<type,NR,NC> m;
for (long row = 0; row < m.nr(); ++row)
{
for (long col = 0; col < m.nc(); ++col)
{
m(row,col) = static_cast<type>(rnd.get_random_double());
}
}
return m;
}
// ----------------------------------------------------------------------------------------
template <typename matrix_type, typename U>
void test_eigenvalue_impl ( const matrix_type& m, const eigenvalue_decomposition<U>& test )
{
typedef typename matrix_type::type type;
const type eps = 10*max(abs(m))*sqrt(std::numeric_limits<type>::epsilon());
dlog << LDEBUG << "test_eigenvalue(): " << m.nr() << " x " << m.nc() << " eps: " << eps;
print_spinner();
DLIB_TEST(test.dim() == m.nr());
// make sure all the various ways of asking for the eigenvalues are actually returning a
// consistent set of eigenvalues.
DLIB_TEST(equal(real(test.get_eigenvalues()), test.get_real_eigenvalues(), eps));
DLIB_TEST(equal(imag(test.get_eigenvalues()), test.get_imag_eigenvalues(), eps));
DLIB_TEST(equal(real(diag(test.get_d())), test.get_real_eigenvalues(), eps));
DLIB_TEST(equal(imag(diag(test.get_d())), test.get_imag_eigenvalues(), eps));
matrix<type> eig1 ( real_eigenvalues(m));
matrix<type> eig2 ( test.get_real_eigenvalues());
sort(&eig1(0), &eig1(0) + eig1.size());
sort(&eig2(0), &eig2(0) + eig2.size());
DLIB_TEST(max(abs(eig1 - eig2)) < eps);
const matrix<type> V = test.get_pseudo_v();
const matrix<type> D = test.get_pseudo_d();
const matrix<complex<type> > CV = test.get_v();
const matrix<complex<type> > CD = test.get_d();
const matrix<complex<type> > CM = complex_matrix(m, uniform_matrix<type>(m.nr(),m.nc(),0));
DLIB_TEST(V.nr() == test.dim());
DLIB_TEST(V.nc() == test.dim());
DLIB_TEST(D.nr() == test.dim());
DLIB_TEST(D.nc() == test.dim());
// CD is a diagonal matrix
DLIB_TEST(diagm(diag(CD)) == CD);
// verify that these things are actually eigenvalues and eigenvectors of m
DLIB_TEST_MSG(max(abs(m*V - V*D)) < eps, max(abs(m*V - V*D)) << " " << eps);
DLIB_TEST(max(norm(CM*CV - CV*CD)) < eps);
// if m is a symmetric matrix
if (max(abs(m-trans(m))) < 1e-5)
{
dlog << LTRACE << "m is symmetric";
// there aren't any imaginary eigenvalues
DLIB_TEST(max(abs(test.get_imag_eigenvalues())) < eps);
DLIB_TEST(diagm(diag(D)) == D);
// only check the determinant against the eigenvalues for small matrices
// because for huge ones the determinant might be so big it overflows a floating point number.
if (m.nr() < 50)
{
const type mdet = det(m);
DLIB_TEST_MSG(std::abs(prod(test.get_real_eigenvalues()) - mdet) < std::abs(mdet)*sqrt(std::numeric_limits<type>::epsilon()),
std::abs(prod(test.get_real_eigenvalues()) - mdet) <<" eps: " << std::abs(mdet)*sqrt(std::numeric_limits<type>::epsilon())
<< " mdet: "<< mdet << " prod(eig): " << prod(test.get_real_eigenvalues())
);
}
// V is orthogonal
DLIB_TEST(equal(V*trans(V), identity_matrix<type>(test.dim()), eps));
DLIB_TEST(equal(m , V*D*trans(V), eps));
}
else
{
dlog << LTRACE << "m is NOT symmetric";
DLIB_TEST_MSG(equal(m , V*D*inv(V), eps), max(abs(m - V*D*inv(V))));
}
}
// ----------------------------------------------------------------------------------------
template <typename matrix_type>
void test_eigenvalue ( const matrix_type& m )
{
typedef typename matrix_type::type type;
typedef typename matrix_type::mem_manager_type MM;
matrix<type,matrix_type::NR, matrix_type::NC, MM, row_major_layout> mr(m);
matrix<type,matrix_type::NR, matrix_type::NC, MM, column_major_layout> mc(m);
{
eigenvalue_decomposition<matrix_type> test(mr);
test_eigenvalue_impl(mr, test);
eigenvalue_decomposition<matrix_type> test_symm(make_symmetric(mr));
test_eigenvalue_impl(make_symmetric(mr), test_symm);
}
{
eigenvalue_decomposition<matrix_type> test(mc);
test_eigenvalue_impl(mc, test);
eigenvalue_decomposition<matrix_type> test_symm(make_symmetric(mc));
test_eigenvalue_impl(make_symmetric(mc), test_symm);
}
}
// ----------------------------------------------------------------------------------------
void matrix_test_double()
{
test_eigenvalue(10*randm<double>(1,1));
test_eigenvalue(10*randm<double>(2,2));
test_eigenvalue(10*randm<double>(3,3));
test_eigenvalue(10*randm<double>(4,4));
test_eigenvalue(10*randm<double>(15,15));
test_eigenvalue(10*randm<double>(150,150));
test_eigenvalue(10*randm<double,1,1>());
test_eigenvalue(10*randm<double,2,2>());
test_eigenvalue(10*randm<double,3,3>());
}
// ----------------------------------------------------------------------------------------
void matrix_test_float()
{
test_eigenvalue(10*randm<float>(1,1));
test_eigenvalue(10*randm<float>(2,2));
test_eigenvalue(10*randm<float>(3,3));
test_eigenvalue(10*randm<float>(4,4));
test_eigenvalue(10*randm<float>(15,15));
test_eigenvalue(10*randm<float>(50,50));
test_eigenvalue(10*randm<float,1,1>());
test_eigenvalue(10*randm<float,2,2>());
test_eigenvalue(10*randm<float,3,3>());
}
template <int dims>
void test_eigenvalue2()
{
for (int seed = 0; seed < 10; ++seed)
{
print_spinner();
matrix<double> H = gaussian_randm(dims,dims,seed);
H = H*trans(H);
eigenvalue_decomposition<matrix<double> > eig(H);
matrix<double> HH = eig.get_pseudo_v()*diagm(eig.get_real_eigenvalues())*trans(eig.get_pseudo_v());
DLIB_TEST_MSG(max(abs(H - HH))<1e-12, "dims: " << dims << " error: " << max(abs(H - HH)));
}
}
// ----------------------------------------------------------------------------------------
class matrix_tester : public tester
{
public:
matrix_tester (
) :
tester ("test_matrix_eig",
"Runs tests on the matrix eigen decomp component.")
{
//rnd.set_seed(cast_to_string(time(0)));
}
void perform_test (
)
{
dlog << LINFO << "seed string: " << rnd.get_seed();
dlog << LINFO << "begin testing with double";
matrix_test_double();
dlog << LINFO << "begin testing with float";
matrix_test_float();
test_eigenvalue2<10>();
test_eigenvalue2<11>();
test_eigenvalue2<3>();
test_eigenvalue2<2>();
test_eigenvalue2<1>();
}
} a;
}