#! /usr/bin/env python

import openturns as ot

ot.TESTPREAMBLE()


def quadM(m, n):
    res = ot.Matrix(m, n)
    for i in range(m):
        for j in range(n):
            res[i, j] = (i + 1.0) ** (j + 1.0)
    return res


def testQR(m, n, full, keep):
    matrix1 = quadM(m, n)
    print("M=", matrix1)
    if keep:
        Q, R = matrix1.computeQR(full)
    else:
        Q, R = matrix1.computeQRInPlace(full)
    print("full=", full, "keep=", keep)
    print("Q= ", Q)
    print("R=", R)
    print("Q*R=", Q * R)
    if keep:
        print("M2=", matrix1)


# Square case
matrix1 = quadM(3, 3)
matrix1.setName("matrix1")
print("matrix1 = ", repr(matrix1))

result1 = matrix1.computeSingularValues()
print("svd (svd only)= ", repr(result1))

result1, u, v = matrix1.computeSVD(True)
print("svd (svd + U, V full)= ", repr(result1))
result1, u, v = matrix1.computeSVD(False)
print("svd (svd + U, V small)= ", repr(result1), ", U=", repr(u), ", v=", repr(v))

# Rectangular case, m < n
matrix1 = quadM(3, 5)
matrix1.setName("matrix1")
print("matrix1 = ", repr(matrix1))

result1 = matrix1.computeSingularValues()
print("svd (svd only)= ", repr(result1))

result1, u, v = matrix1.computeSVD(True)
print("svd (svd + U, V full)= ", repr(result1))
result1, u, v = matrix1.computeSVD(False)
print("svd (svd + U, V small)= ", repr(result1), ", U=", repr(u), ", v=", repr(v))

# Rectangular case, m > n
matrix1 = quadM(5, 3)
matrix1.setName("matrix1")
print("matrix1 = ", repr(matrix1))

result1 = matrix1.computeSingularValues()
print("svd (svd only)= ", repr(result1))

result1, u, v = matrix1.computeSVD(True)
print("svd (svd + U, V full)= ", repr(result1))
# result1, u, v = matrix1.computeSVD(False)
# print "svd (svd + U, V small)= ", repr(result1), ", U=", repr(u), ",
# v=", repr(v)

for iFull in range(2):
    for iKeep in range(2):
        testQR(3, 3, iFull == 1, iKeep == 1)
        testQR(3, 5, iFull == 1, iKeep == 1)
        testQR(5, 3, iFull == 1, iKeep == 1)