1# Copyright 2011-2013 Kwant authors. 

2# 

3# This file is part of Kwant. It is subject to the license terms in the file 

4# LICENSE.rst found in the top-level directory of this distribution and at 

5# https://kwant-project.org/license. A list of Kwant authors can be found in 

6# the file AUTHORS.rst at the top-level directory of this distribution and at 

7# https://kwant-project.org/authors. 

8 

9__all__ = ['lll', 'cvp', 'voronoi'] 

10 

11import numpy as np 

12from scipy import linalg as la 

13from itertools import product 

14 

15 

16def gs_coefficient(a, b): 

17 """Gram-Schmidt coefficient.""" 

18 return np.dot(a, b) / np.linalg.norm(b)**2 

19 

20 

21def gs(mat): 

22 """Compute Gram-Schmidt decomposition on a matrix.""" 

23 mat = np.copy(mat) 

24 for i in range(len(mat)): 

25 for j in range(i): 

26 mat[i] -= gs_coefficient(mat[i], mat[j]) * mat[j] 

27 return mat 

28 

29 

30def is_c_reduced(vecs, c): 

31 """Check if a basis is c-reduced.""" 

32 vecs = gs(vecs) 

33 r = np.apply_along_axis(lambda x: np.linalg.norm(x)**2, 1, vecs) 

34 return np.all((r[: -1] / r[1:]) < c) 

35 

36 

37def lll(basis, c=1.34): 

38 """ 

39 Calculate a reduced lattice basis using LLL algorithm. 

40 

41 Reduce a basis of a lattice to an almost orthonormal form. For details see 

42 e.g. https://en.wikipedia.org/wiki/Lenstra–Lenstra–Lovász_lattice_basis_reduction_algorithm. 

43 

44 Parameters 

45 ---------- 

46 basis : 2d array-like of floats 

47 The lattice basis vectors to be reduced. 

48 c : float 

49 Reduction parameter for the algorithm. Must be larger than 1 1/3, 

50 since otherwise a solution is not guaranteed to exist. 

51 

52 Returns 

53 ------- 

54 reduced_basis : numpy array 

55 The basis vectors of the LLL-reduced basis. 

56 transformation : numpy integer array 

57 Coefficient matrix for tranforming from the reduced basis to the 

58 original one. 

59 """ 

60 vecs_orig = np.asarray(basis, float) 

61 if vecs_orig.ndim != 2: 61 ↛ 62line 61 didn't jump to line 62, because the condition on line 61 was never true

62 raise ValueError('"basis" must be a 2d array-like object.') 

63 if vecs_orig.shape[0] > vecs_orig.shape[1]: 63 ↛ 64line 63 didn't jump to line 64, because the condition on line 63 was never true

64 raise ValueError('The number of basis vectors exceeds the ' 

65 'space dimensionality.') 

66 vecs = vecs_orig.copy() 

67 vecsstar = vecs_orig.copy() 

68 m = vecs.shape[0] 

69 u = np.identity(m) 

70 def ll_reduce(i): 

71 for j in reversed(range(i)): 

72 vecs[i] -= np.round(u[i, j]) * vecs[j] 

73 u[i] -= np.round(u[i, j]) * u[j] 

74 

75 # Initialize values. 

76 for i in range(m): 

77 for j in range(i): 

78 u[i, j] = gs_coefficient(vecs[i], vecsstar[j]) 

79 vecsstar[i] -= u[i, j] * vecsstar[j] 

80 ll_reduce(i) 

81 

82 # Main part of LLL algorithm. 

83 i = 0 

84 while i < m-1: 

85 if (np.linalg.norm(vecsstar[i]) ** 2 < 

86 c * np.linalg.norm(vecsstar[i+1]) ** 2): 

87 i += 1 

88 else: 

89 vecsstar[i+1] += u[i+1, i] * vecsstar[i] 

90 

91 u[i, i] = gs_coefficient(vecs[i], vecsstar[i+1]) 

92 u[i, i+1] = u[i+1, i] = 1 

93 u[i+1, i+1] = 0 

94 vecsstar[i] -= u[i, i] * vecsstar[i+1] 

95 vecs[[i, i+1]] = vecs[[i+1, i]] 

96 vecsstar[[i, i+1]] = vecsstar[[i+1, i]] 

97 u[[i, i+1]] = u[[i+1, i]] 

98 for j in range(i+2, m): 

99 u[j, i] = gs_coefficient(vecs[j], vecsstar[i]) 

100 u[j, i+1] = gs_coefficient(vecs[j], vecsstar[i+1]) 

101 if abs(u[i+1, i]) > 0.5: 

102 ll_reduce(i+1) 

103 i = max(i-1, 0) 

104 coefs = np.linalg.lstsq(vecs_orig.T, vecs.T, rcond=None)[0] 

105 if not np.allclose(np.round(coefs), coefs, atol=1e-6): 105 ↛ 106line 105 didn't jump to line 106, because the condition on line 105 was never true

106 raise RuntimeError('LLL algorithm instability.') 

107 if not is_c_reduced(vecs, c): 107 ↛ 108line 107 didn't jump to line 108, because the condition on line 107 was never true

108 raise RuntimeError('LLL algorithm instability.') 

109 return vecs, np.array(np.round(coefs), int) 

110 

111 

112def cvp(vec, basis, n=1, group_by_length=False, rtol=1e-09): 

113 """ 

114 Solve the n-closest vector problem for a vector, given a basis. 

115 

116 This algorithm performs poorly in general, so it should be supplied 

117 with LLL-reduced basis. 

118 

119 Parameters 

120 ---------- 

121 vec : 1d array-like of floats 

122 The lattice vectors closest to this vector are to be found. 

123 basis : 2d array-like of floats 

124 Sequence of basis vectors 

125 n : int 

126 Number of lattice vectors closest to the point that need to be found. If 

127 `group_by_length=True`, all vectors with the `n` shortest distinct distances 

128 are found. 

129 group_by_length : bool 

130 Whether to count points with distances that are within `rtol` as one 

131 towards `n`. Useful for finding `n`-th nearest neighbors of high symmetry points. 

132 rtol : float 

133 If `group_by_length=True`, relative tolerance when deciding whether 

134 distances are equal. 

135 

136 Returns 

137 ------- 

138 coords : numpy array 

139 An array with the coefficients of the `n` closest lattice vectors, or all 

140 the lattice vectors with the `n` shortest distances from the 

141 requested point. 

142 

143 Notes 

144 ----- 

145 This function can also be used to solve the `n` shortest lattice vector 

146 problem if the `vec` is zero, and `n+1` points are requested 

147 (and the first output is ignored). 

148 """ 

149 # Calculate coordinates of the starting point in this basis. 

150 basis = np.asarray(basis) 

151 if basis.ndim != 2: 151 ↛ 152line 151 didn't jump to line 152, because the condition on line 151 was never true

152 raise ValueError('`basis` must be a 2d array-like object.') 

153 vec = np.asarray(vec) 

154 # Project the coordinates on the space spanned by `basis`, in 

155 # units of `basis` vectors. 

156 vec_bas = la.lstsq(basis.T, vec)[0] 

157 # Round the coordinates, this finds the lattice point which is 

158 # the center of the parallelepiped that contains `vec`. 

159 # `vec` is surely in the parallelepiped with edges `basis` 

160 # centered at `center_coords`. 

161 center_coords = np.array(np.round(vec_bas), int) 

162 # Make `vec` the projection of `vec` onto the space spanned by `basis`. 

163 vec = vec_bas @ basis 

164 bbt = basis @ basis.T 

165 # The volume of a parallelepiped spanned by `basis` is `sqrt(det(bbt))`. 

166 # (We use this instead of `det(basis)` because basis does not necessarily 

167 # span the full space, hence not always a square matrix.) 

168 # We calculate the radius `rad` of the largest sphere that can be inscribed 

169 # in the parallelepiped spanned by `basis`. The area of each face spanned 

170 # by one less vector in `basis` is given by the above formula applied to `bbt` 

171 # with one row and column discarded (first minor). We utilizie the relation 

172 # between entries of the inverse matrix and first minors. The perpendicular 

173 # size of the parallelepiped is given by the volume of the parallelepiped 

174 # divided by the area of the face. We choose the smallest perpendicular 

175 # size as the diameter of the largest sphere inside. 

176 rad = 0.5 / np.sqrt(np.max(np.diag(la.inv(bbt)))) 

177 

178 l = 1 

179 while True: 

180 # Make all the lattice points in and on the edges of a parallelepiped 

181 # of `2*l` size around `center_coords`. 

182 points = np.mgrid[tuple(slice(i - l, i + l + 1) for i in center_coords)] 

183 points = points.reshape(basis.shape[0], -1).T 

184 # If there are less than `n` points, make more. 

185 if len(points) < n: 

186 l += 1 

187 continue 

188 point_coords = points @ basis 

189 point_coords = point_coords - vec.T 

190 distances = la.norm(point_coords, axis = 1) 

191 order = np.argsort(distances) 

192 distances = distances[order] 

193 points = points[order] 

194 if group_by_length: 

195 # Group points that are equidistant within `rtol * rad`. 

196 group_boundaries = np.where(np.diff(distances) > rtol * rad)[0] 

197 if len(group_boundaries) + 1 < n: 

198 # Make more points if there are less than `n` groups. 

199 l += 1 

200 continue 

201 elif len(group_boundaries) == n - 1: 

202 # If there are exactly `n` groups, choose largest distance. 

203 dist_n = distances[-1] 

204 points_keep = points 

205 else: 

206 # If there are more than `n` groups, choose the largest 

207 # distance in the `n`-th group. 

208 dist_n = distances[group_boundaries[n-1]] 

209 # Only return points that are in the first `n` groups. 

210 points_keep = points[:group_boundaries[n-1] + 1] 

211 else: 

212 rtol = 0 

213 # Choose the `n`-th distance. 

214 dist_n = distances[n-1] 

215 # Only return the first `n` points. 

216 points_keep = points[:n] 

217 # We covered all points within radius `(2*l-1) * rad` from the 

218 # parallelepiped spanned by `basis` centered at `center_coords`. 

219 # Because `vec` is guaranteed to be in this parallelepiped, if 

220 # the current `dist_n` is smaller than `(2*l-1) * rad`, we surely 

221 # found all the smaller vectors. 

222 if dist_n < (2*l - 1 - rtol) * rad: 

223 return np.array(points_keep, int) 

224 else: 

225 # Otherwise there may be smaller vectors we haven't found, 

226 # increase `l`. 

227 l += 1 

228 continue 

229 

230 

231def voronoi(basis, reduced=False, rtol=1e-09): 

232 """ 

233 Return an array of lattice vectors forming its voronoi cell. 

234 

235 Parameters 

236 ---------- 

237 basis : 2d array-like of floats 

238 Basis vectors for which the Voronoi neighbors have to be found. 

239 reduced : bool 

240 If False, exactly `2 (2^D - 1)` vectors are returned (where `D` 

241 is the number of lattice vectors), these are not always the minimal 

242 set of Voronoi vectors. 

243 If True, only the minimal set of Voronoi vectors are returned. 

244 rtol : float 

245 Tolerance when deciding whether a vector is in the minimal set, 

246 vectors associated with small faces compared to the size of the 

247 unit cell are discarded. 

248 

249 Returns 

250 ------- 

251 voronoi_neighbors : numpy array of ints 

252 All the lattice vectors that (potentially) neighbor the origin. 

253 List of length `2n` where the second half of the list contains 

254 is `-1` times the vectors in the first half. 

255 """ 

256 basis = np.asarray(basis) 

257 if basis.ndim != 2: 257 ↛ 258line 257 didn't jump to line 258, because the condition on line 257 was never true

258 raise ValueError('`basis` must be a 2d array-like object.') 

259 # Find halved lattice points, every face of the VC contains half 

260 # of the lattice vector which is the normal of the face. 

261 # These points are all potentially on a face a VC, 

262 # but not necessarily the VC centered at the origin. 

263 displacements = list(product(*(len(basis) * [[0, .5]])))[1:] 

264 # Find the nearest lattice point, this is the lattice point whose 

265 # VC this face belongs to. 

266 vertices = np.array([cvp(vec @ basis, basis)[0] for vec in 

267 displacements]) 

268 # The lattice vector for a face is exactly twice the vector to the 

269 # closest lattice point from the halved lattice point on the face. 

270 vertices = np.array(np.round((vertices - displacements) * 2), int) 

271 if reduced: 271 ↛ 275line 271 didn't jump to line 275, because the condition on line 271 was never true

272 # Discard vertices that are not associated with a face of the VC. 

273 # This happens if half of the lattice vector is outside the convex 

274 # polytope defined by the rest of the vertices in `keep`. 

275 bbt = basis @ basis.T 

276 products = vertices @ bbt @ vertices.T 

277 keep = np.array([True] * len(vertices)) 

278 for i in range(len(vertices)): 

279 # Relevant if the projection of half of the lattice vector `vertices[i]` 

280 # onto every other lattice vector in `veritces[keep]` is less than `0.5`. 

281 mask = np.array([False if k == i else b for k, b in enumerate(keep)]) 

282 projections = 0.5 * products[i, mask] / np.diag(products)[mask] 

283 if not np.all(np.abs(projections) < 0.5 - rtol): 

284 keep[i] = False 

285 vertices = vertices[keep] 

286 

287 vertices = np.concatenate([vertices, -vertices]) 

288 return vertices