Coverage for src\pqlattice\linalg\_linalg.py: 91%

1from functools import reduce 

2 

3import numpy as np 

4 

5from ..typing import Matrix, SquareMatrix, Vector, validate_aliases 

6from ._utils import row_add, row_scale, row_swap 

7 

8 

9@validate_aliases 

10def hnf(A: Matrix) -> tuple[Matrix, SquareMatrix]: 

11 """_summary_ 

12 

13 Parameters 

14 ---------- 

15 A : Matrix 

16 _description_ 

17 

18 Returns 

19 ------- 

20 tuple[Matrix, SquareMatrix] 

21 _description_ 

22 """ 

23 H, U, _ = _hnf(A) 

24 return H, U 

25 

26 

27@validate_aliases 

28def _hnf(a: Matrix) -> tuple[Matrix, SquareMatrix, int]: 

29 H = np.array(a, dtype=object) 

30 m, n = H.shape 

31 U = np.eye(m, dtype=object) # The transformation matrix 

32 pivot_row = 0 

33 pivot_col = 0 

34 det_U = 1 

35 

36 while pivot_row < m and pivot_col < n: 

37 # pivot selection 

38 if np.all(H[pivot_row:, pivot_col] == 0): 

39 pivot_col += 1 

40 continue 

41 

42 candidates = [(abs(H[i, pivot_col]), i) for i in range(pivot_row, m) if H[i, pivot_col] != 0] 

43 _, best_row = min(candidates) 

44 

45 row_swap(H, pivot_row, best_row) 

46 row_swap(U, pivot_row, best_row) 

47 det_U *= -1 

48 

49 # clear below pivot 

50 for i in range(pivot_row + 1, m): 

51 while H[i, pivot_col] != 0: 

52 factor = H[i, pivot_col] // H[pivot_row, pivot_col] 

53 

54 row_add(H, i, pivot_row, -factor) 

55 row_add(U, i, pivot_row, -factor) 

56 

57 if H[i, pivot_col] != 0: 

58 row_swap(H, pivot_row, i) 

59 row_swap(U, pivot_row, i) 

60 det_U *= -1 

61 

62 if H[pivot_row, pivot_col] < 0: 

63 row_scale(H, pivot_row, -1) 

64 row_scale(U, pivot_row, -1) 

65 det_U *= -1 

66 

67 pivot_val = H[pivot_row, pivot_col] 

68 

69 # Reduce rows above pivot 

70 for i in range(pivot_row): 

71 factor = H[i, pivot_col] // pivot_val 

72 row_add(H, i, pivot_row, -factor) 

73 row_add(U, i, pivot_row, -factor) 

74 

75 pivot_row += 1 

76 pivot_col += 1 

77 

78 return H, U, det_U 

79 

80 

81@validate_aliases 

82def left_kernel(a: Matrix): 

83 """_summary_ 

84 

85 Parameters 

86 ---------- 

87 a : Matrix 

88 _description_ 

89 

90 Returns 

91 ------- 

92 _type_ 

93 _description_ 

94 """ 

95 return right_kernel(a.T) 

96 

97 

98@validate_aliases 

99def right_kernel(a: Matrix) -> Matrix: 

100 """_summary_ 

101 

102 Parameters 

103 ---------- 

104 a : Matrix 

105 _description_ 

106 

107 Returns 

108 ------- 

109 Matrix 

110 _description_ 

111 """ 

112 H, U = hnf(a.T) 

113 kernel_basis: list[Vector] = [] 

114 

115 m, _ = H.shape 

116 for i in range(m): 

117 if np.all(H[i] == 0): 

118 kernel_basis.append(U[i]) 

119 

120 return np.array(kernel_basis, dtype=object) 

121 

122 

123@validate_aliases 

124def left_nullity(a: Matrix) -> int: 

125 """_summary_ 

126 

127 Parameters 

128 ---------- 

129 a : Matrix 

130 _description_ 

131 

132 Returns 

133 ------- 

134 int 

135 _description_ 

136 """ 

137 kernel = left_kernel(a) 

138 return kernel.shape[0] 

139 

140 

141@validate_aliases 

142def right_nullity(a: Matrix) -> int: 

143 """_summary_ 

144 

145 Parameters 

146 ---------- 

147 a : Matrix 

148 _description_ 

149 

150 Returns 

151 ------- 

152 int 

153 _description_ 

154 """ 

155 kernel = right_kernel(a) 

156 return kernel.shape[0] 

157 

158 

159def rank(a: Matrix) -> int: 

160 """_summary_ 

161 

162 Parameters 

163 ---------- 

164 a : Matrix 

165 _description_ 

166 

167 Returns 

168 ------- 

169 int 

170 _description_ 

171 """ 

172 m, n = a.shape 

173 l_rank = m - left_nullity(a) 

174 r_rank = n - right_nullity(a) 

175 assert l_rank == r_rank 

176 return l_rank 

177 

178 

179@validate_aliases 

180def det(A: SquareMatrix) -> int: 

181 """_summary_ 

182 

183 Parameters 

184 ---------- 

185 A : SquareMatrix 

186 _description_ 

187 

188 Returns 

189 ------- 

190 int 

191 _description_ 

192 """ 

193 H, _, det_U = _hnf(A) 

194 

195 return reduce(lambda a, b: a * b, np.diagonal(H), 1) * det_U 

196 

197 

198@validate_aliases 

199def minor(A: SquareMatrix, i: int, j: int) -> int: 

200 """_summary_ 

201 

202 Parameters 

203 ---------- 

204 A : SquareMatrix 

205 _description_ 

206 i : int 

207 _description_ 

208 j : int 

209 _description_ 

210 

211 Returns 

212 ------- 

213 int 

214 _description_ 

215 """ 

216 return det(np.delete(np.delete(A, i, axis=0), j, axis=1)) 

217 

218 

219@validate_aliases 

220def cofactor(A: SquareMatrix, i: int, j: int) -> int: 

221 """_summary_ 

222 

223 Parameters 

224 ---------- 

225 A : SquareMatrix 

226 _description_ 

227 i : int 

228 _description_ 

229 j : int 

230 _description_ 

231 

232 Returns 

233 ------- 

234 int 

235 _description_ 

236 """ 

237 return minor(A, i, j) * ((-1) ** (i + 1 + j + 1)) 

238 

239 

240@validate_aliases 

241def cofactor_matrix(A: SquareMatrix) -> SquareMatrix: 

242 """_summary_ 

243 

244 Parameters 

245 ---------- 

246 A : SquareMatrix 

247 _description_ 

248 

249 Returns 

250 ------- 

251 SquareMatrix 

252 _description_ 

253 """ 

254 n = A.shape[0] 

255 C = np.zeros((n, n), dtype=object) 

256 for i in range(n): 

257 for j in range(n): 

258 C[i, j] = cofactor(A, i, j) 

259 return C