import math
import unittest

import numpy as np


def det(mat):
    """ det(a) --> d """
    assert isinstance(mat, np.ndarray)
    assert len(mat.shape) == 2
    assert mat.shape[0] == mat.shape[1]
    val_a, val_b, val_c, val_d = mat.reshape((-1,))
    return val_a * val_d - val_b * val_c


def det_recursive(mat):
    """ det(a)  = b1 * det (b) + b2 * det(c) + b3* det(d)
    𝒪(n!)
    """
    rows, _ = mat.shape
    if rows == 2:
        return det(mat)
    row = mat[0, :]
    det_mat = 0.0
    sign = 1
    for i, val in enumerate(row):
        # print(i, val)
        mat_i = np.zeros((rows-1, rows-1))
        mat_i[:, 0:i] = mat[1:, 0:i]
        mat_i[:, i:] = mat[1:, i+1:]
        # print(mat_i)
        det_i = det_recursive(mat_i)
        det_mat += sign * val * det_i
        sign *= -1
    return det_mat


def num_ops_det(num):
    if num == 1:
        num_ops = 1
    elif num == 2:
        num_ops = num_ops_det(num - 1) * num + (num - 1) + 1
    elif num == 3:
        num_ops = num_ops_det(num - 1) * num + (num - 1) + 1
    elif num >= 4:
        num_ops = num_ops_det(num - 1) * num + (num - 1) + 1
    else:
        num_ops = -1
    return num_ops


def demo():
    mat_a = np.array([[1, 2], [3, 4]])
    print(det(mat_a))

    mat_b = np.array([[789]])
    print(np.linalg.det(mat_b))

    det_num_ops = [num_ops_det(num) for num in range(2, 16)]
    print(det_num_ops)

    val_fact = [math.factorial(num) for num in range(2, 16)]
    print(val_fact)

    vec_det = np.array(det_num_ops, dtype=float)
    vec_fac = np.array(val_fact, dtype=float)

    print(vec_det/vec_fac)


class TestDeterminant(unittest.TestCase):

    def test_det_2x2(self):
        mat_a = np.random.randint(1, 10, size=(2, 2))
        nrows, ncols = mat_a.shape
        self.assertEqual(nrows, ncols)
        self.assertEqual(nrows, 2)

        val_a, val_b, val_c, val_d = mat_a.reshape((-1,))
        det_a = val_a * val_d - val_b * val_c
        self.assertEqual(det(mat_a), det_a)
        np.testing.assert_allclose(
            np.linalg.det(mat_a), det(mat_a))

    def test_det_lu_nxn(self):
        for num in range(3, 9+1):
            mat_a = np.random.randint(-5, 5, size=(num, num))
            mat_l, mat_u = np.tril(mat_a, k=-1) + np.eye(num, dtype=int), \
                np.triu(mat_a)
            mat_a = mat_l @ mat_u
            det_u = np.prod(np.diag(mat_u))
            np.testing.assert_allclose(np.linalg.det(mat_a), det_u, atol=1e-7)
            self.assertAlmostEqual(np.linalg.det(mat_a), det_u)

    def test_det_nxn(self):
        for num in range(3, 9+1):
            mat_a = np.random.randint(1, 10, size=(num, num))
            self.assertAlmostEqual(np.linalg.det(mat_a), det_recursive(mat_a))


if __name__ == '__main__':
    unittest.main(verbosity=2)
