root/library/tests/square_mat_stress.cpp @ 468

#include "../bdm/math/square_mat.h"
#include "../bdm/math/chmat.h"
#include "base/user_info.h"
#include "square_mat_point.h"
#include "UnitTest++.h"
#include "TestReporterStdout.h"
#include <iostream>
#include <iomanip>
#include <stdlib.h>
#include <string.h>

using std::cout;
using std::cerr;
using std::endl;

using bdm::UIFile;
using bdm::UI;

const char *agenda_file_name = "agenda.cfg";
double epsilon = 0.00001;
bool fast = false;

namespace bdm {
UIREGISTER(square_mat_point);
}

namespace UnitTest
{

// can't include mat_checks.h because CheckClose is different in this file
extern bool AreClose(const itpp::vec &expected, const itpp::vec &actual,
                     double tolerance);

extern bool AreClose(const itpp::mat &expected, const itpp::mat &actual,
                     double tolerance);

void CheckClose(TestResults &results, const itpp::mat &expected,
                const itpp::mat &actual, double tolerance,
                TestDetails const& details) {
    if (!AreClose(expected, actual, tolerance)) { 
        MemoryOutStream stream;
        stream << "failed at " << expected.rows()
               << " x " << expected.cols();

        results.OnTestFailure(details, stream.GetText());
    }
}

}

typedef void (*FTestMatrix)(int, square_mat_point *);

template<typename TMatrix>
void test_matrix(int index, square_mat_point *point) {
    Real_Timer tt;
        
    cout << "agenda[" << index << "]:" << endl;
    mat A = point->get_matrix();
    int sz = A.rows();
    CHECK_EQUAL(A.cols(), sz);

    tt.tic();
    TMatrix sqmat(A);
    double elapsed = tt.toc();
    cout << "ctor(" << sz << " x " << sz << "): " << elapsed << " s" << endl;

    tt.tic();
    mat res = sqmat.to_mat();
    elapsed = tt.toc();

    if (!fast) {
        CHECK_CLOSE(A, res, epsilon);
    }

    cout << "to_mat: " << elapsed << " s" << endl;

    vec v = point->get_vector();
    double w = point->get_scalar();
    TMatrix sqmat2 = sqmat;
        
    tt.tic();
    sqmat2.opupdt(v, w);
    elapsed = tt.toc();

    if (!fast) {
        mat expA = A + w * outer_product(v, v);
        CHECK_CLOSE(expA, sqmat2.to_mat(), epsilon);
    }

    cout << "opupdt: " << elapsed << " s" << endl;

    TMatrix invmat(sz);

    tt.tic();
    sqmat.inv(invmat);
    elapsed = tt.toc();

    mat invA;
    if (!fast) {
        invA = inv(A);
        CHECK_CLOSE(invA, invmat.to_mat(), epsilon);
    }

    cout << "inv: " << elapsed << " s" << endl;

    tt.tic();
    double ld = sqmat.logdet();
    elapsed = tt.toc();

    if (!fast) {
        double d = det(A);
        CHECK_CLOSE(log(d), ld, epsilon);
    }

    cout << "logdet: " << elapsed << " s" << endl;

    tt.tic();
    double q = sqmat.qform(ones(sz));
    elapsed = tt.toc();

    if (!fast) {
        CHECK_CLOSE(sumsum(A), q, epsilon);
    }

    cout << "qform(1): " << elapsed << " s" << endl;

    tt.tic();
    q = sqmat.qform(v);
    elapsed = tt.toc();

    if (!fast) {
        double r = (A * v) * v;
        CHECK_CLOSE(r, q, epsilon);
    }

    cout << "qform(v): " << elapsed << " s" << endl;

    tt.tic();
    q = sqmat.invqform(v);
    elapsed = tt.toc();

    if (!fast) {
        double r = (invA * v) * v;
        CHECK_CLOSE(r, q, epsilon);
    }

    cout << "invqform: " << elapsed << " s" << endl;

    TMatrix twice = sqmat;

    tt.tic();
    twice += sqmat;
    elapsed = tt.toc();

    if (!fast) {
        res = 2 * A;
        CHECK_CLOSE(res, twice.to_mat(), epsilon);
    }

    cout << "+=: " << elapsed << " s" << endl;

    sqmat2 = sqmat;

    tt.tic();
    sqmat2.mult_sym(A);
    elapsed = tt.toc();

    if (!fast) {
        res = (A * A) * A.T();
        CHECK_CLOSE(res, sqmat2.to_mat(), epsilon);
    }

    cout << "^2: " << elapsed << " s" << endl;
}

void test_agenda(FTestMatrix test) {
    UIFile fag(agenda_file_name);
    Array<square_mat_point *> mag;
    UI::get(mag, fag, "agenda");
    int sz = mag.size();
    CHECK(sz > 0);
    for (int i = 0; i < sz; ++i) {
        test(i, mag(i));
    }

    for (int i = 0; i < sz; ++i) {
        square_mat_point *p = mag(i);
        mag(i) = 0;
        delete p;
    }
}

SUITE(ldmat) {
    TEST(agenda) {
        test_agenda(test_matrix<ldmat>);
    }
}

SUITE(fsqmat) {
    TEST(agenda) {
        test_agenda(test_matrix<fsqmat>);
    }
}

SUITE(chmat) {
    TEST(agenda) {
        test_agenda(test_matrix<chmat>);
    }
}

int main(int argc, char const *argv[]) {
    bool unknown = false;
    int update_next = 0; // 1 suite, 2 epsilon, 3 agenda file
    const char *suite = "ldmat";
    const char **param = argv + 1;
    while (*param && !unknown) {
        if (update_next) {
            if (update_next == 1) {
                suite = *param;
            } else if (update_next == 2) {
                double eps = atof(*param);
                if (eps > 0) {
                    epsilon = eps;
                } else {
                    cerr << "invalid epsilon value ignored" << endl;
                }
            } else {
                agenda_file_name = *param;
            }

            update_next = 0;
        } else {
            if (!strcmp(*param, "-a")) {
                update_next = 3;
            } else if (!strcmp(*param, "-c")) {
                update_next = 1;
            } else if (!strcmp(*param, "-e")) {
                update_next = 2;
            } else if (!strcmp(*param, "-f")) {
                fast = true;
            } else {
                unknown = true;
            }
        }

        ++param;
    }

    if (unknown || update_next) {
        cerr << "usage: " << argv[0] << " [ -f ] [ -e epsilon ] [ -a agenda_input.cfg ] [ -c class ]" << endl;
    } else {
        UnitTest::TestReporterStdout reporter;
        UnitTest::TestRunner runner(reporter);
        return runner.RunTestsIf(UnitTest::Test::GetTestList(),
            suite,
            UnitTest::True(),
            0);
    }
}