root/library/tests/square_mat_stress.cpp @ 467

#include "../bdm/math/square_mat.h"
#include "../bdm/math/chmat.h"
#include "base/user_info.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 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());
    }
}

}

template<typename TMatrix>
void test_matrix(int index, const mat &A) {
    Real_Timer tt;
        
    cout << "agenda[" << index << ']' << endl;
    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 = randu(sz);
    double w = randu();
    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();

    if (!fast) {
        mat 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;
}

template<typename TMatrix>
void test_agenda() {
    UIFile fag(agenda_file_name);
    Array<mat> mag;
    UI::get(mag, fag, "agenda");
    int sz = mag.size();
    CHECK(sz > 0);
    for (int i = 0; i < sz; ++i) {
        test_matrix<TMatrix>(i, mag(i));
    }
}

SUITE(ldmat) {
    TEST(cycle) {
        test_agenda<ldmat>();
    }
}

SUITE(fsqmat) {
    TEST(cycle) {
        test_agenda<fsqmat>();
    }
}

SUITE(chmat) {
    TEST(cycle) {
        test_agenda<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, "-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 ] [ -c class ]" << endl;
    } else {
        UnitTest::TestReporterStdout reporter;
        UnitTest::TestRunner runner(reporter);
        return runner.RunTestsIf(UnitTest::Test::GetTestList(),
            suite,
            UnitTest::True(),
            0);
    }
}