Coverage for src/pqlattice/linalg/_linalg.py: 84%

1from functools import reduce 

2 

3import numpy as np 

4 

5from .._utils import as_integer 

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

7from ._utils import row_add, row_scale, row_swap 

8 

9 

10@validate_aliases 

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

12 """_summary_ 

13 

14 Parameters 

15 ---------- 

16 A : Matrix 

17 _description_ 

18 

19 Returns 

20 ------- 

21 tuple[Matrix, SquareMatrix] 

22 _description_ 

23 """ 

24 H, U, _ = _hnf(A) 

25 return H, U 

26 

27 

28@validate_aliases 

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

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

31 m, n = H.shape 

32 U = np.eye(m, dtype=object) 

33 pivot_row = 0 

34 pivot_col = 0 

35 det_U = 1 

36 

37 while pivot_row < m and pivot_col < n: 

38 # pivot selection 

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

40 pivot_col += 1 

41 continue 

42 

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

44 _, best_row = min(candidates) 

45 

46 row_swap(H, pivot_row, best_row) 

47 row_swap(U, pivot_row, best_row) 

48 det_U *= -1 

49 

50 # clear below pivot 

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

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

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

54 

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

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

57 

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

59 row_swap(H, pivot_row, i) 

60 row_swap(U, pivot_row, i) 

61 det_U *= -1 

62 

63 if H[pivot_row, pivot_col] < 0: 

64 row_scale(H, pivot_row, -1) 

65 row_scale(U, pivot_row, -1) 

66 det_U *= -1 

67 

68 pivot_val = H[pivot_row, pivot_col] 

69 

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 right_image(A: Matrix) -> Matrix: 

83 """ 

84 _summary_ 

85 

86 Parameters 

87 ---------- 

88 A : Matrix 

89 _description_ 

90 

91 Returns 

92 ------- 

93 Matrix 

94 _description_ 

95 """ 

96 return left_image(A.T).T 

97 

98 

99@validate_aliases 

100def left_image(A: Matrix) -> Matrix: 

101 """ 

102 _summary_ 

103 

104 Parameters 

105 ---------- 

106 A : Matrix 

107 _description_ 

108 

109 Returns 

110 ------- 

111 Matrix 

112 _description_ 

113 """ 

114 H, _ = hnf(A) 

115 

116 m, _ = H.shape 

117 k = 0 

118 for i in range(m): 

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

120 k = i 

121 break 

122 

123 return as_integer(H[:k]) 

124 

125 

126@validate_aliases 

127def left_kernel(A: Matrix): 

128 """ 

129 {x : xA = 0} 

130 

131 Parameters 

132 ---------- 

133 A : Matrix 

134 _description_ 

135 

136 Returns 

137 ------- 

138 _type_ 

139 _description_ 

140 """ 

141 H, U = hnf(A) 

142 kernel_basis: list[Vector] = [] 

143 

144 m, _ = H.shape 

145 for i in range(m): 

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

147 kernel_basis.append(U[i]) 

148 

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

150 

151 

152@validate_aliases 

153def right_kernel(A: Matrix) -> Matrix: 

154 """ 

155 {x : Ax = 0} 

156 

157 Parameters 

158 ---------- 

159 A : Matrix 

160 _description_ 

161 

162 Returns 

163 ------- 

164 Matrix 

165 _description_ 

166 """ 

167 return left_kernel(A.T) 

168 

169 

170@validate_aliases 

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

172 """_summary_ 

173 

174 Parameters 

175 ---------- 

176 a : Matrix 

177 _description_ 

178 

179 Returns 

180 ------- 

181 int 

182 _description_ 

183 """ 

184 kernel = left_kernel(a) 

185 return kernel.shape[0] 

186 

187 

188@validate_aliases 

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

190 """_summary_ 

191 

192 Parameters 

193 ---------- 

194 a : Matrix 

195 _description_ 

196 

197 Returns 

198 ------- 

199 int 

200 _description_ 

201 """ 

202 kernel = right_kernel(a) 

203 return kernel.shape[0] 

204 

205 

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

207 """_summary_ 

208 

209 Parameters 

210 ---------- 

211 a : Matrix 

212 _description_ 

213 

214 Returns 

215 ------- 

216 int 

217 _description_ 

218 """ 

219 m, n = a.shape 

220 l_rank = m - left_nullity(a) 

221 r_rank = n - right_nullity(a) 

222 assert l_rank == r_rank 

223 return l_rank 

224 

225 

226@validate_aliases 

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

228 """_summary_ 

229 

230 Parameters 

231 ---------- 

232 A : SquareMatrix 

233 _description_ 

234 

235 Returns 

236 ------- 

237 int 

238 _description_ 

239 """ 

240 H, _, det_U = _hnf(A) 

241 

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

243 

244 

245@validate_aliases 

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

247 """_summary_ 

248 

249 Parameters 

250 ---------- 

251 A : SquareMatrix 

252 _description_ 

253 i : int 

254 _description_ 

255 j : int 

256 _description_ 

257 

258 Returns 

259 ------- 

260 int 

261 _description_ 

262 """ 

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

264 

265 

266@validate_aliases 

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

268 """_summary_ 

269 

270 Parameters 

271 ---------- 

272 A : SquareMatrix 

273 _description_ 

274 i : int 

275 _description_ 

276 j : int 

277 _description_ 

278 

279 Returns 

280 ------- 

281 int 

282 _description_ 

283 """ 

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

285 

286 

287@validate_aliases 

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

289 """_summary_ 

290 

291 Parameters 

292 ---------- 

293 A : SquareMatrix 

294 _description_ 

295 

296 Returns 

297 ------- 

298 SquareMatrix 

299 _description_ 

300 """ 

301 n = A.shape[0] 

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

303 for i in range(n): 

304 for j in range(n): 

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

306 return C