/** @file test_suite.cpp
 *  @brief Test suite for project 2
 *
 *  This is where all the tests for the project are implemented.
 *  The tests test the creation of symmetrical tridiagonals,
 *  finding the largest off-diagonal absolute value of a symmetric matrix,
 *  and the Jacobi rotation method.
 *
 *  @author Cory Alexander Balaton (coryab)
 *  @author Janita Ovidie Sandtrøen Willumsen (janitaws)
 *  @bug No known bugs
 */
#include <cassert>
#include <cmath>
#include <iostream>

#include "utils.hpp"
#include "matrix.hpp"
#include "jacobi.hpp"

void test_create_symmetric_tridiagonal() {
    double tol = 10e-8;
    double diff;
    int N = 6;
    double h = 1. / (double) (N+1);
    double a = -1. / (h*h), d = 2. / (h*h);

    arma::mat A = create_symmetric_tridiagonal(N, a, d);

    arma::vec eigval;
    arma::mat eigvec;

    arma::eig_sym(eigval, eigvec, A);

    // Test eigenvalues
    for (int i=0; i < N; i++) {
        diff = eigval(i) - (d + 2.*a*std::cos(((i+1.)*M_PI)/(N+1.)));
        assert(std::abs(diff) < tol); 
    }

    // Test eigenvectors
    arma::vec v, res, analytic_vec = arma::vec(N);
    for (int i=0; i < N; i++) {
        v = eigvec.col(i);

        // Build the analytic vector
        for (int j=0; j < N; j++) {
            analytic_vec(j) = std::sin(((j+1.)*(i+1.)*M_PI) / (N+1.));
        }

        analytic_vec = arma::normalise(analytic_vec);

        // Flip the sign of the analytic vector if they are different
        if (analytic_vec(0)*v(0) < 0.) {
            analytic_vec *= -1;
        }

        res = v - analytic_vec;
        for (int j=0; j < N; j++) {
            assert(std::abs(res(j)) < tol);
        }
    }
}

void test_max_off_diag_symmetric() {
    arma::mat A = arma::mat(4,4,arma::fill::eye);

    double tol = 10e-8;
    A(0, 3) = .5;
    A(3, 0) = .5;
    A(1, 2) = -.7;
    A(2, 1) = -.7;
    
    int k,l;
    assert(max_offdiag_symmetric(A, k, l) - .7 < tol);
    assert(k == 1 && l == 2);
}

void test_jacobi_eigensolver()
{
    double tol = 10e-8;
    double diff;
    int N = 6;
    double h = 1. / (double) (N+1);
    double a = -1. / (h*h), d = 2. / (h*h);

    arma::mat A = create_symmetric_tridiagonal(N, a, d);

    arma::vec eigval;
    arma::mat eigvec;

    int iters;
    bool converged;
    jacobi_eigensolver(A, 10e-14, eigval, eigvec, 10000, iters, converged);

    // Test eigenvalues
    for (int i=0; i < N; i++) {
        diff = eigval(i) - (d + 2.*a*std::cos(((i+1.)*M_PI)/(N+1.)));
        assert(std::abs(diff) < tol); 
    }

    // Test eigenvectors
    arma::vec v, res, analytic_vec = arma::vec(N);
    for (int i=0; i < N; i++) {
        v = eigvec.col(i);

        // Build the analytic vector
        for (int j=0; j < N; j++) {
            analytic_vec(j) = std::sin(((j+1.)*(i+1.)*M_PI) / (N+1.));
        }

        analytic_vec = arma::normalise(analytic_vec);

        // Flip the sign of the analytic vector if they are different
        if (analytic_vec(0)*v(0) < 0.) {
            analytic_vec *= -1;
        }

        res = v - analytic_vec;
        for (int j=0; j < N; j++) {
            assert(std::abs(res(j)) < tol);
        }
    }
}

int main() {

    test_create_symmetric_tridiagonal();
    test_max_off_diag_symmetric();
    test_jacobi_eigensolver();

    return 0;
}