121 lines
2.8 KiB
Python
121 lines
2.8 KiB
Python
from typing import List, Tuple
|
|
import numpy as np
|
|
|
|
|
|
def det(matrix: np.matrix):
|
|
if matrix.shape[0] != matrix.shape[1]:
|
|
return None
|
|
if matrix.shape[0] == 1:
|
|
return matrix[0, 0]
|
|
|
|
return _det(matrix)
|
|
|
|
|
|
def _det(matrix: np.matrix):
|
|
matrix, swaps = triangle(np.copy(matrix))
|
|
if matrix is None:
|
|
return 0.0
|
|
|
|
n = matrix.shape[0]
|
|
determinant = (-1) ** swaps
|
|
for i in range(n):
|
|
determinant *= matrix[i, i]
|
|
|
|
return determinant
|
|
|
|
|
|
def triangle(matrix: np.matrix):
|
|
swaps = 0
|
|
n = matrix.shape[0]
|
|
|
|
for i in range(n):
|
|
if matrix[i, i] == 0:
|
|
if not swap_first_non_zero(matrix, i)[0]:
|
|
return None, swaps
|
|
swaps += 1
|
|
|
|
for j in range(i + 1, n):
|
|
matrix[j, :] = matrix[j, :] - matrix[i, :] * (matrix[j, i] / matrix[i, i])
|
|
|
|
return matrix, swaps
|
|
|
|
|
|
def swap_first_non_zero(matrix, i):
|
|
n = matrix.shape[0]
|
|
for j in range(i, n):
|
|
if matrix[j, i] != 0:
|
|
if i != j:
|
|
matrix[[i, j], :] = matrix[[j, i], :]
|
|
return True, j
|
|
return False, -1
|
|
|
|
|
|
def swap_max_row(matrix, j):
|
|
n = matrix.shape[0]
|
|
|
|
current = float("inf")
|
|
if matrix[j, j] != 0:
|
|
current = np.sum(np.abs(matrix[j, :])) / abs(matrix[j, j])
|
|
max_i = j
|
|
for i in range(j + 1, n):
|
|
if matrix[i, j] == 0 or matrix[j, i] == 0:
|
|
continue
|
|
|
|
target_current = np.sum(np.abs(matrix[i, :])) / abs(matrix[i, i])
|
|
new = np.sum(np.abs(matrix[i, :])) / abs(matrix[i, j])
|
|
target_new = np.sum(np.abs(matrix[j, :])) / abs(matrix[j, i])
|
|
|
|
if current - new > target_new - target_current:
|
|
current = new
|
|
max_i = i
|
|
|
|
if j != max_i:
|
|
matrix[[j, max_i], :] = matrix[[max_i, j], :]
|
|
|
|
return max_i
|
|
|
|
|
|
def swap_max_column(matrix, i):
|
|
n = matrix.shape[0]
|
|
|
|
current = float("inf")
|
|
if matrix[i, i] != 0:
|
|
current = np.sum(np.abs(matrix[i, :])) / abs(matrix[i, i])
|
|
best_j = i
|
|
for j in range(i + 1, n):
|
|
if matrix[i, j] == 0 or matrix[j, i] == 0:
|
|
continue
|
|
|
|
target_current = np.sum(np.abs(matrix[j, :])) / abs(matrix[j, j])
|
|
new = np.sum(np.abs(matrix[i, :])) / abs(matrix[i, j])
|
|
target_new = np.sum(np.abs(matrix[j, :])) / abs(matrix[j, i])
|
|
|
|
if current - new > target_new - target_current:
|
|
current = new
|
|
best_j = j
|
|
|
|
if i != best_j:
|
|
matrix[:, [i, best_j]] = matrix[:, [best_j, i]]
|
|
|
|
return best_j
|
|
|
|
|
|
def norm(matrix: np.matrix):
|
|
n = matrix.shape[0]
|
|
m = matrix.shape[1]
|
|
|
|
m_norm = float("-inf")
|
|
for i in range(n):
|
|
s = 0.0
|
|
for j in range(m):
|
|
s += abs(matrix[i, j])
|
|
|
|
if s > m_norm:
|
|
m_norm = s
|
|
|
|
return m_norm
|
|
|
|
|
|
def cond(matrix: np.matrix):
|
|
return norm(matrix) * norm(np.linalg.inv(matrix))
|