Skip to content

Commit 3a1a3c9

Browse files
[Andrew Polyakov] Added code samples for numpy testing
1 parent 1a34b8d commit 3a1a3c9

File tree

6 files changed

+238
-0
lines changed

6 files changed

+238
-0
lines changed
Lines changed: 71 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,71 @@
1+
import numpy as np
2+
3+
4+
def solve_linear_system(matrix: np.ndarray) -> np.ndarray:
5+
"""
6+
Solve a linear system of equations using Gaussian elimination with partial pivoting
7+
8+
Args:
9+
- matrix: Coefficient matrix with the last column representing the constants.
10+
11+
Returns:
12+
- Solution vector.
13+
14+
Raises:
15+
- ValueError: If the matrix is not correct (i.e., singular).
16+
17+
https://courses.engr.illinois.edu/cs357/su2013/lect.htm Lecture 7
18+
19+
Example:
20+
>>> A = np.array([[2, 1, -1], [-3, -1, 2], [-2, 1, 2]], dtype=float)
21+
>>> B = np.array([8, -11, -3], dtype=float)
22+
>>> solution = solve_linear_system(np.column_stack((A, B)))
23+
>>> np.allclose(solution, np.array([2., 3., -1.]))
24+
True
25+
>>> solve_linear_system(np.array([[0, 0], [0, 0]], dtype=float))
26+
array([nan, nan])
27+
"""
28+
ab = np.copy(matrix)
29+
num_of_rows = ab.shape[0]
30+
num_of_columns = ab.shape[1] - 1
31+
x_lst: list[float] = []
32+
33+
for column_num in range(num_of_rows):
34+
for i in range(column_num, num_of_columns):
35+
if abs(ab[i][column_num]) > abs(ab[column_num][column_num]):
36+
ab[[column_num, i]] = ab[[i, column_num]]
37+
if ab[column_num, column_num] == 0.0:
38+
raise ValueError("Matrix is not correct")
39+
else:
40+
pass
41+
if column_num != 0:
42+
for i in range(column_num, num_of_rows):
43+
ab[i, :] -= (
44+
ab[i, column_num - 1]
45+
/ ab[column_num - 1, column_num - 1]
46+
* ab[column_num - 1, :]
47+
)
48+
49+
for column_num in range(num_of_rows):
50+
for i in range(column_num, num_of_columns):
51+
if abs(ab[i][column_num]) > abs(ab[column_num][column_num]):
52+
ab[[column_num, i]] = ab[[i, column_num]]
53+
if ab[column_num, column_num] == 0.0:
54+
raise ValueError("Matrix is not correct")
55+
else:
56+
pass
57+
if column_num != 0:
58+
for i in range(column_num, num_of_rows):
59+
ab[i, :] -= (
60+
ab[i, column_num - 1]
61+
/ ab[column_num - 1, column_num - 1]
62+
* ab[column_num - 1, :]
63+
)
64+
65+
for column_num in range(num_of_rows - 1, -1, -1):
66+
x = ab[column_num, -1] / ab[column_num, column_num]
67+
x_lst.insert(0, x)
68+
for i in range(column_num - 1, -1, -1):
69+
ab[i, -1] -= ab[i, column_num] * x
70+
71+
return np.asarray(x_lst)
Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,24 @@
1+
import numpy as np
2+
3+
4+
def lower_upper_decomposition(table: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
5+
rows, columns = np.shape(table)
6+
if rows != columns:
7+
msg = (
8+
"'table' has to be of square shaped array but got a "
9+
f"{rows}x{columns} array:\n{table}"
10+
)
11+
raise ValueError(msg)
12+
lower = np.zeros((rows, columns))
13+
upper = np.zeros((rows, columns))
14+
for i in range(columns):
15+
for j in range(i):
16+
total = np.sum(lower[i, :i] * upper[:i, j])
17+
if upper[j][j] == 0:
18+
raise ArithmeticError("No LU decomposition exists")
19+
lower[i][j] = (table[i][j] - total) / upper[j][j]
20+
lower[i][i] = 1
21+
for j in range(i, columns):
22+
total = np.sum(lower[i, :i] * upper[:i, j])
23+
upper[i][j] = table[i][j] - total
24+
return lower, upper
Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,37 @@
1+
import numpy as np
2+
3+
4+
def _is_matrix_spd(matrix: np.ndarray) -> bool:
5+
"""
6+
Returns True if input matrix is symmetric positive definite.
7+
Returns False otherwise.
8+
9+
For a matrix to be SPD, all eigenvalues must be positive.
10+
11+
>>> import numpy as np
12+
>>> matrix = np.array([
13+
... [4.12401784, -5.01453636, -0.63865857],
14+
... [-5.01453636, 12.33347422, -3.40493586],
15+
... [-0.63865857, -3.40493586, 5.78591885]])
16+
>>> _is_matrix_spd(matrix)
17+
True
18+
>>> matrix = np.array([
19+
... [0.34634879, 1.96165514, 2.18277744],
20+
... [0.74074469, -1.19648894, -1.34223498],
21+
... [-0.7687067 , 0.06018373, -1.16315631]])
22+
>>> _is_matrix_spd(matrix)
23+
False
24+
"""
25+
# Ensure matrix is square.
26+
assert np.shape(matrix)[0] == np.shape(matrix)[1]
27+
28+
# If matrix not symmetric, exit right away.
29+
if np.allclose(matrix, matrix.T) is False:
30+
return False
31+
32+
# Get eigenvalues and eignevectors for a symmetric matrix.
33+
eigen_values, _ = np.linalg.eigh(matrix)
34+
35+
# Check sign of all eigenvalues.
36+
# np.all returns a value of type np.bool_
37+
return bool(np.all(eigen_values > 0))
Lines changed: 76 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,76 @@
1+
def rank_of_matrix(matrix: list[list[int | float]]) -> int:
2+
"""
3+
Finds the rank of a matrix.
4+
Args:
5+
matrix: The matrix as a list of lists.
6+
Returns:
7+
The rank of the matrix.
8+
Example:
9+
>>> matrix1 = [[1, 2, 3],
10+
... [4, 5, 6],
11+
... [7, 8, 9]]
12+
>>> rank_of_matrix(matrix1)
13+
2
14+
>>> matrix2 = [[1, 0, 0],
15+
... [0, 1, 0],
16+
... [0, 0, 0]]
17+
>>> rank_of_matrix(matrix2)
18+
2
19+
>>> matrix3 = [[1, 2, 3, 4],
20+
... [5, 6, 7, 8],
21+
... [9, 10, 11, 12]]
22+
>>> rank_of_matrix(matrix3)
23+
2
24+
>>> rank_of_matrix([[2,3,-1,-1],
25+
... [1,-1,-2,4],
26+
... [3,1,3,-2],
27+
... [6,3,0,-7]])
28+
4
29+
>>> rank_of_matrix([[2,1,-3,-6],
30+
... [3,-3,1,2],
31+
... [1,1,1,2]])
32+
3
33+
>>> rank_of_matrix([[2,-1,0],
34+
... [1,3,4],
35+
... [4,1,-3]])
36+
3
37+
>>> rank_of_matrix([[3,2,1],
38+
... [-6,-4,-2]])
39+
1
40+
>>> rank_of_matrix([[],[]])
41+
0
42+
>>> rank_of_matrix([[1]])
43+
1
44+
>>> rank_of_matrix([[]])
45+
0
46+
"""
47+
48+
rows = len(matrix)
49+
columns = len(matrix[0])
50+
rank = min(rows, columns)
51+
52+
for row in range(rank):
53+
# Check if diagonal element is not zero
54+
if matrix[row][row] != 0:
55+
# Eliminate all the elements below the diagonal
56+
for col in range(row + 1, rows):
57+
multiplier = matrix[col][row] / matrix[row][row]
58+
for i in range(row, columns):
59+
matrix[col][i] -= multiplier * matrix[row][i]
60+
else:
61+
# Find a non-zero diagonal element to swap rows
62+
reduce = True
63+
for i in range(row + 1, rows):
64+
if matrix[i][row] != 0:
65+
matrix[row], matrix[i] = matrix[i], matrix[row]
66+
reduce = False
67+
break
68+
if reduce:
69+
rank -= 1
70+
for i in range(rows):
71+
matrix[i][row] = matrix[i][rank]
72+
73+
# Reduce the row pointer by one to stay on the same row
74+
row -= 1
75+
76+
return rank
Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,23 @@
1+
import numpy as np
2+
from numpy import ndarray
3+
4+
5+
def norm_squared(vector: ndarray) -> float:
6+
"""
7+
Return the squared second norm of vector
8+
norm_squared(v) = sum(x * x for x in v)
9+
10+
Args:
11+
vector (ndarray): input vector
12+
13+
Returns:
14+
float: squared second norm of vector
15+
16+
>>> norm_squared([1, 2])
17+
5
18+
>>> norm_squared(np.asarray([1, 2]))
19+
5
20+
>>> norm_squared([0, 0])
21+
0
22+
"""
23+
return np.dot(vector, vector)
Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
1+
import numpy as np
2+
3+
4+
def function_to_test_default(a: np.ndarray):
5+
if a.shape == (2, 1,):
6+
return 2
7+
return a.shape

0 commit comments

Comments
 (0)