<commit_before># TODO:
# *implement high-level qr() and lu()
# *cache LU decompositon
# *test unitvector
# *iterative solving
# *iterative improving of solution

from __future__ import division

from mptypes import *
from matrices import matrix, swap_row, extend, mnorm_1, norm_p
from copy import copy

def LU_decomp(A):
    """
    LU-factorization of a n*n matrix using the Gauss algorithm.
    Returns L and U in one matrix and the pivot indices.
    """
    if not A.rows == A.cols:
        raise ValueError('need n*n matrix')
    tol = mnorm_1(A) * eps # each pivot element has to be bigger than this
    n = A.rows
    p = [None]*(n - 1)
    for j in xrange(n - 1):
        # pivoting, choose max(abs(reciprocal row sum)*abs(pivot element))
        biggest = 0
        for k in xrange(j, n):
            current = 1/sum([abs(A[k,l]) for l in xrange(j, n)]) * abs(A[k,j])
            if current > biggest: # TODO: what if equal?
                biggest = current
                p[j] = k
        # swap rows according to p
        swap_row(A, j, p[j])
        if abs(A[j,j]) < tol:
            raise ZeroDivisionError('matrix is numerically singular')
        # calculate elimination factors and add rows
        for i in xrange(j + 1, n):
            A[i,j] /= A[j,j]
            for k in xrange(j + 1, n):
                A[i,k] -= A[i,j]*A[j,k]
    return A, p

def L_solve(L, b, p=None):
    """
    Solve the lower part of a LU factorized matrix for y.
    """
    L.rows == L.cols, 'need n*n matrix'
    n = L.rows
    assert len(b) == n
    if p: # swap b according to p
        for k in xrange(0, len(p)):
            swap_row(b, k, p[k])
    # solve
    for i in xrange(1, n):
        for j in xrange(i):
            b[i] -= L[i,j] * b[j]
    return b

def U_solve(U, y):
    """
    Solve the upper part of a LU factorized matrix for x.
    """
    assert U.rows == U.cols, 'need n*n matrix'
    n = U.rows
    assert len(y) == n
    x = copy(y)
    for i in xrange(n - 1, -1, -1):
        for j in xrange(i + 1, n):
            x[i] -= U[i,j] * x[j]
        x[i] /= U[i,i]
    return x

@extraprec(10)
def lu_solve(A, b, **kwargs):
    """
    Ax = b => x

    Solve a determined or overdetermined linear equations system.
    Fast LU decomposition is used, which is less accurate than QR decomposition
    (especially for overdetermined systems), but it's twice as efficient.
    Use qr_solve if you want more precison or have to solve a very ill-
    conditioned system.
    """
    # do not overwrite A nor b
    A, b = matrix(A, **kwargs).copy(), matrix(b, **kwargs).copy()
    if A.rows < A.cols:
        raise ValueError('cannot solve underdetermined system')
    if A.rows > A.cols:
        # use least-squares method if overdetermined
        # (this increases errors)
        AT = A.T
        A = AT * A
        b = AT * b
        return cholesky_solve(A, b)
    else:
        # LU factorization
        A, p = LU_decomp(A)
        b = L_solve(A, b, p)
        x = U_solve(A, b)
        return x

def unitvector(n, i):
    """
    Return the i-th n-dimensional unit vector.
    """
    assert 0 < i <= n, 'this unit vector does not exist'
    return [0]*(i-1) + [1] + [0]*(n-i)

@extraprec(10)
def inverse(A, **kwargs):
    """
    Calculate the inverse of a matrix.

    If you want to solve an equation system Ax = b, it's recommended to use
    solve(A, b) instead, it's about 3 times more efficient.
    """
    # do not overwrite A
    A = matrix(A, **kwargs).copy()
    n = A.rows
    # get LU factorisation
    A, p = LU_decomp(A)
    cols = []
    # calculate unit vectors and solve corresponding system to get columns
    for i in xrange(1, n + 1):
        e = unitvector(n, i)
        y = L_solve(A, e, p)
        cols.append(U_solve(A, y))
    # convert columns to matrix
    inv = []
    for i in xrange(n):
        row = []
        for j in xrange(n):
            row.append(cols[j][i])
        inv.append(row)
    return matrix(inv, **kwargs)

def householder(A):
    """
    (A|b) -> H, p, x, res
    """
    assert isinstance(A, matrix)
    m = A.rows
    n = A.cols
    assert m >= n - 1
    # calculate Householder matrix
    p = []
    for j in xrange(0, n - 1):
        s = 0.
        for i in xrange(j, m):
            s += (A[i,j])**2<commit_msg>improved householder() docstring<commit_after># TODO:
# *implement high-level qr() and lu()
# *cache LU decompositon
# *test unitvector
# *iterative solving
# *iterative improving of solution

from __future__ import division

from mptypes import *
from matrices import matrix, swap_row, extend, mnorm_1, norm_p
from copy import copy

def LU_decomp(A):
    """
    LU-factorization of a n*n matrix using the Gauss algorithm.
    Returns L and U in one matrix and the pivot indices.
    """
    if not A.rows == A.cols:
        raise ValueError('need n*n matrix')
    tol = mnorm_1(A) * eps # each pivot element has to be bigger than this
    n = A.rows
    p = [None]*(n - 1)
    for j in xrange(n - 1):
        # pivoting, choose max(abs(reciprocal row sum)*abs(pivot element))
        biggest = 0
        for k in xrange(j, n):
            current = 1/sum([abs(A[k,l]) for l in xrange(j, n)]) * abs(A[k,j])
            if current > biggest: # TODO: what if equal?
                biggest = current
                p[j] = k
        # swap rows according to p
        swap_row(A, j, p[j])
        if abs(A[j,j]) < tol:
            raise ZeroDivisionError('matrix is numerically singular')
        # calculate elimination factors and add rows
        for i in xrange(j + 1, n):
            A[i,j] /= A[j,j]
            for k in xrange(j + 1, n):
                A[i,k] -= A[i,j]*A[j,k]
    return A, p

def L_solve(L, b, p=None):
    """
    Solve the lower part of a LU factorized matrix for y.
    """
    L.rows == L.cols, 'need n*n matrix'
    n = L.rows
    assert len(b) == n
    if p: # swap b according to p
        for k in xrange(0, len(p)):
            swap_row(b, k, p[k])
    # solve
    for i in xrange(1, n):
        for j in xrange(i):
            b[i] -= L[i,j] * b[j]
    return b

def U_solve(U, y):
    """
    Solve the upper part of a LU factorized matrix for x.
    """
    assert U.rows == U.cols, 'need n*n matrix'
    n = U.rows
    assert len(y) == n
    x = copy(y)
    for i in xrange(n - 1, -1, -1):
        for j in xrange(i + 1, n):
            x[i] -= U[i,j] * x[j]
        x[i] /= U[i,i]
    return x

@extraprec(10)
def lu_solve(A, b, **kwargs):
    """
    Ax = b => x

    Solve a determined or overdetermined linear equations system.
    Fast LU decomposition is used, which is less accurate than QR decomposition
    (especially for overdetermined systems), but it's twice as efficient.
    Use qr_solve if you want more precison or have to solve a very ill-
    conditioned system.
    """
    # do not overwrite A nor b
    A, b = matrix(A, **kwargs).copy(), matrix(b, **kwargs).copy()
    if A.rows < A.cols:
        raise ValueError('cannot solve underdetermined system')
    if A.rows > A.cols:
        # use least-squares method if overdetermined
        # (this increases errors)
        AT = A.T
        A = AT * A
        b = AT * b
        return cholesky_solve(A, b)
    else:
        # LU factorization
        A, p = LU_decomp(A)
        b = L_solve(A, b, p)
        x = U_solve(A, b)
        return x

def unitvector(n, i):
    """
    Return the i-th n-dimensional unit vector.
    """
    assert 0 < i <= n, 'this unit vector does not exist'
    return [0]*(i-1) + [1] + [0]*(n-i)

@extraprec(10)
def inverse(A, **kwargs):
    """
    Calculate the inverse of a matrix.

    If you want to solve an equation system Ax = b, it's recommended to use
    solve(A, b) instead, it's about 3 times more efficient.
    """
    # do not overwrite A
    A = matrix(A, **kwargs).copy()
    n = A.rows
    # get LU factorisation
    A, p = LU_decomp(A)
    cols = []
    # calculate unit vectors and solve corresponding system to get columns
    for i in xrange(1, n + 1):
        e = unitvector(n, i)
        y = L_solve(A, e, p)
        cols.append(U_solve(A, y))
    # convert columns to matrix
    inv = []
    for i in xrange(n):
        row = []
        for j in xrange(n):
            row.append(cols[j][i])
        inv.append(row)
    return matrix(inv, **kwargs)

def householder(A):
    """
    (A|b) -> H, p, x, res

    (A|b) is the coefficient matrix with left hand side of an optionally
    overdetermined linear equation system.
    H and p contain all information about the transformation matrices.
    x is the solution, res the residual.
    """
    assert isinstance(A, matrix)
    m = A.rows
    n = A.cols
    assert m >= n - 1
    # calculate Householder matrix
    p = []
    for j in xrange(0, n - 1):
        s = 0.
        for i in xrange(j, m):
            s += (A[i,j])**2
{'fraction_non_alphanumeric': 0.12440570522979398, 'fraction_numerical': 0.0070183382386235, 'mean_word_length': 2.0319148936170213, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '716661', 'n_tokens_mistral': 3143, 'n_tokens_neox': 3020, 'n_words': 1088}