import numpy as np
from utils import qr_pos

def is_diagonal(A):
    return np.count_nonzero(A - np.diag(np.diagonal(A))) == 0

def test_eig(n=30, tol=0.00001):
    M = np.random.rand(N, N) # Random matrix
    M = 2*M - 1              # Entries from -1 to 1
    Q, R = qr_pos(M)         # QR dec. so that diag(R) is positive
    if is_diagonal(R):
        print 'R is diagonal!'
    RTR = np.transpose(R) * R
    sinvals = np.sqrt(np.linalg.eigvals(RTR))
    diff = np.linalg.norm(sinvals - np.diag(R))

    # print 'Singular values:', sinvals
    # print 'Diagonal of R:', np.diag(R)
    # print 'Difference:', diff
    return diff < tol, M

success = True
N = 300
for i in range(N):
    equal, M = test_eig()
    print 'i',
    if not equal:
        print
        print 'Case found where the hypothesis does not hold:', M
        success = False
        break
if success:
    print
    print 'Hypothesis holds for', N, 'random test cases.'

i i i i ... i i i
Hypothesis holds for 300 random test cases.