mail@pastecode.io avatar
2 months ago
12 kB
from typing import List, Callable, TypeVar, Tuple, Optional, Union
from functools import reduce
from math import isinf, isnan

from common.utils import randomize_int, fast_power
from matrix.Vector import Vector
from matrix.Matrix import Matrix

T = TypeVar('T')

class SqMatrix(Matrix):
    def __init__(self, content):
        """Set self.dim and self.content, assert that number of rows equal number of columns
        Also initialize superclass properties on self."""
        self.dim = len(content)
        self.content = tuple(tuple(row) for row in content)
        assert all(len(row) == self.dim for row in content)

    def __repr__(self):
        """Print a square matrix as 'SqMatrix(...)' with one row of content per line,
        and indentation."""
        lines = ",\n        ".join(f"{row}" for row in self.content)
        return f"SqMatrix(({lines}))"

    def random(dim: int, randomizer: Callable[[], int] = randomize_int()) -> 'SqMatrix':
        """Create and return a new square matrix of the given size, being filled
        with randomly selected elements.  Each such element is the return
        value of the given randomizer function.   randomizer will be called
        once for each entry in the resulting matrix."""
        content = [[randomizer() for _ in range(dim)] for _ in range(dim)]
        return SqMatrix(content)

    def identity(dim: int) -> 'SqMatrix':
        """Return a square matrix of the given size with 1 (integer)
        along the main diagonal, and 0 (integer) off the main diagonal."""
        content = [[1 if i == j else 0 for j in range(dim)] for i in range(dim)]
        return SqMatrix(content)

    def zero(dim: int) -> 'SqMatrix':
        """Return a square matrix of the given size consisting entirely
        of 0 (integer)."""
        content = [[0 for _ in range(dim)] for _ in range(dim)]
        return SqMatrix(content)

    def tabulate(dim: int, f: Callable[[int, int], float]) -> 'SqMatrix':
        """Create and return a new square matrix of the given size.
        The SqMatrix is filled with the return value of the given function, f.
        The function, f, is called once per entry in the resulting matrix,
        and the return value of f(i,j) is the value of m[i][j]"""
        content = [[f(i, j) for j in range(dim)] for i in range(dim)]
        return SqMatrix(content)

    def diagonal(entries) -> 'SqMatrix':
        """Given a sequence of values, return a square matrix with those
        values along the main diagonal.   If the sequence has size n,
        then the resulting matrix has dimension n as well.   The diagonal
        matrix has 0 (integer) off the main diagonal."""
        dim = len(entries)
        content = [[entries[i] if i == j else 0 for j in range(dim)] for i in range(dim)]
        return SqMatrix(content)

    def gaussian_elimination_back_substitution(self, v: Vector) -> Optional[Vector]:
        """Adjoin the given column vector to the matrix, self.
        Then perform elementary row operations on the given matrix to reduce
        it to row echelon form.   Then perform back substitution to return
        the vector which solves the system."""
        assert isinstance(v, Vector)
        assert v.dim == self.dim

        # compute an echelon form of the matrix, self.
        # this form has zeros below the main diagonal.
        # and non-zeros on the main diagonal unless the matrix is
        # singular in which case there are zeros on the diagonal.
        ech = self.adjoin_col(v).make_row_echelon()
        if ech is None:
            return None
        assert ech.rows == self.rows
        assert ech.cols == self.cols + 1
        return ech.back_substitution()

    def gauss_jordan_elimination(self, v: Vector) -> Optional[Vector]:
        """Adjoin the given column vector to the matrix, self.
        Then use Gauss Jordan elimination (via elementary row operations)
        to reduce the matrix to diagonal.  This may result in a matrix with
        a zero on the diagonal, in which case None is returned, otherwise
        normalize the diagonal and return the vector which solves the system.
        assert isinstance(v, Vector)
        assert v.dim == self.dim

        diag, _ = self.adjoin_col(v).make_unit_diagonal()
        if diag is None:
            return None
        assert diag.rows == self.rows
        assert diag.cols == self.cols + 1

        return diag.col_vec(v.dim)

    def gauss_jordan_inverse(self) -> Tuple[Optional['SqMatrix'], T]:
        """returns 2-tuple of two values:  (inverse,determinant),
        If the determinant is zero, then inverse=None
        diag, det = self.adjoin_cols(SqMatrix.identity(self.dim)).make_unit_diagonal()
        if diag is None:
            return None, 0
        assert diag.rows == self.rows
        assert diag.cols == 2 * self.cols
        return diag.extract_cols(range(self.dim, 2 * self.dim)), det

    def laplacian_expansion(self, zero: T = 0) -> T:
        """Compute the determinant of the given square matrix using
        Laplacian expansion."""
        if self.dim == 1:
            return self[0][0]
            det = 0
            for j in range(self.dim):
                det += ((-1) ** j) * self[0][j] * self.minor(0, j).laplacian_expansion()
            return det

    def cramers_rule(self, b: Vector) -> Vector:
        """Return a new vector which solves the system of equations Ax=b,
        where A is the matrix self.   If the determinant of A is 0, then None
        is returned.  The solution is found using Cramer's Rule:
        I.e., to compute the k'th component of the returned Vector,
        we replaced the k'th column of a by the column vector b,
        and calculate the determinant of that matrix, then divide by
        the determinant of A.
        assert self.dim == b.dim
        det_A = self.laplacian_expansion()
        if det_A == 0:
            return None

        solutions = []
        for i in range(self.dim):
            mat = self.copy()
            mat.set_col(i, b)
            det_mat = mat.laplacian_expansion()
            solutions.append(det_mat / det_A)

        return Vector(*solutions)

    def power(self, p: int) -> 'SqMatrix':
        """Raise the matrix to the p'th power."""
        assert isinstance(p, int)
        assert p >= 0
        identity_matrix = SqMatrix.identity(self.dim)
        return fast_power(self, lambda x, y: x * y, p, identity_matrix)

    def characteristic_polynomial(self) -> 'Polynomial':
        """Compute the characteristic polynomial of the square matrix.
        This is done by generating the matrix M - λI, which is a SqMatrix
        having a Polynomial in each entry.  Then use laplacian_expansion
        to compute the determinant.
        from matrix.Polynomial import Polynomial
        variable = Polynomial.variable()
        identity_matrix = SqMatrix.identity(self.dim)
        diff_matrix = self - (identity_matrix * variable)
        return diff_matrix.laplacian_expansion().to_polynomial(variable)

    def eigenvalues(self, epsilon: float) -> List[float]:
        """Compute the (some) eigenvalues of the matrix by calling
        characteristic_polynomial to obtain a Polynomial, then using
        the function poly_roots to compute its roots.
        The eigenvalues method returns the roots/eigenvalues in increasing
        Just as poly_roots, this method is allowed to fail to compute
        the eigenvalues in some cases.
        from algebra.roots import poly_roots
        assert isinstance(epsilon, float), f"epsilon={epsilon} is not a float"
        char_poly = self.characteristic_polynomial()
        roots = poly_roots(char_poly, epsilon)
        return sorted(roots)

    def eigenvectors(self, spectrum: List[float], epsilon=0.0000001) -> List[Vector]:
        """Given a list of non-zero eigenvalues, `spectrum`,
        return the corresponding set of eigenvectors.
        It is assumed that spectrum is already sorted in increasing order.
        If the given eigenvalues are not distinct, return None.
        If any eigenvalue is zero (abs(e)<epsilon), return None.
        If there are not exactly self.dim number of eigenvalues given,
          then return None"""
        assert spectrum == sorted(spectrum)
        if len(spectrum) != self.dim:
            return None
        # check for zero eigenvalues
        if any(abs(e) < epsilon for e in spectrum):
            return None
        # check for duplicate eigenvalues
        for k in range(len(spectrum) - 1):
            # we know spectrum[k+1] >= spectrum[k]
            if spectrum[k + 1] - spectrum[k] < epsilon:
                # we found two equal eigenvalues (within epsilon of each other)
                return None

        def condition(m):
            """replace numbers close to zero by zero"""
            return SqMatrix.tabulate(self.dim,
                                     lambda r, c: 0 if abs(m[r][c]) < epsilon else m[r][c])

        id = SqMatrix.identity(self.dim)
        amli = [condition(self - id.scale(spectrum[k]))
                for k in range(self.dim)]  # A - lambda[k] I

        def prod(ms):
            return reduce(lambda a, b: a * b,

        # pamli: product-of-A-minus-lambda-I
        # here we form a list of matrix products.
        # the k-th element in this list is the product of all
        #   elements of amli except the k'th one
        pamli = [prod(ms)
                 for k in range(len(spectrum))
                 for ms in [[amli[j]
                             for j in range(len(spectrum))
                             if j != k]]]

        # each element in pamli is a matrix having a column of zeros.
        # we want to select a column vector from each, but we want to avoid
        # the zero vectors.  thus for each matrix in pamli, we find the column
        # vector with the maximum norm.   This works if all the eigenvalues
        # or distinct; however, if there is a repeated eigenvalue, then this
        # search will result in a repeated eigenvector.  This is a bug which
        # we need to fix.
        vectors = []
        for k in range(len(spectrum)):
            # find all columns of pamli[k],
            # excluding any zero vector in `vectors`
            independent = [cv  # column vector
                           for j in range(pamli[k].cols)
                           for cv in [pamli[k].col_vec(j)]
                           if cv.norm() > epsilon]
            if independent:
                vectors.append(min(independent, key=lambda v: v.norm()))

        return vectors

    def eigen(self, epsilon: float) -> Optional[List[Tuple[float, Vector]]]:
        """Compute (eigenvalue, eigenvector) pairs of the matrix,
        if one (or more) of the eigenvalues is 0, return None.
        By 'e is 0' we mean abs(e) < epsilon.
        If we cannot find self.dim number of unique eigenvalues,
        then return None.
        Else return a list of tuples (e,v) where e is an eigenvalue
        and v is its normalized eigenvector.  Note that if v is a normalized
        eigenvector, then so is -v, so use v or -v as you prefer.
        The list is sorted e1 < e2 < ... < en.

        spectrum = self.eigenvalues(epsilon)
        if len(set(spectrum)) != self.dim:
            return None
        if any(abs(e) < epsilon for e in spectrum):
            return None

        vectors = self.eigenvectors(spectrum)
        return [(spectrum[k], vectors[k])
                for k in range(len(vectors))]

    def minor(self, i, j):
        """Return the minor of the matrix after removing the i-th row and j-th column."""
        return SqMatrix([row[:j] + row[j+1:] for row in (self[:i]+self[i+1:])])
Leave a Comment