Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

""" 

Diamond cutting implementation 

 

AUTHORS: 

 

- Jan Poeschko (2012-07-02): initial version 

""" 

 

#***************************************************************************** 

# Copyright (C) 2012 Jan Poeschko <jan@poeschko.com> 

# 

# Distributed under the terms of the GNU General Public License (GPL) 

# as published by the Free Software Foundation; either version 2 of 

# the License, or (at your option) any later version. 

# http://www.gnu.org/licenses/ 

#***************************************************************************** 

from __future__ import print_function 

 

from sage.geometry.polyhedron.constructor import Polyhedron 

from sage.matrix.constructor import matrix, identity_matrix 

from sage.modules.free_module_element import vector 

from sage.rings.rational_field import QQ 

 

from math import sqrt, floor, ceil 

 

def plane_inequality(v): 

""" 

Return the inequality for points on the same side as the origin 

with respect to the plane through ``v`` normal to ``v``. 

 

EXAMPLES:: 

 

sage: from sage.modules.diamond_cutting import plane_inequality 

sage: ieq = plane_inequality([1, -1]); ieq 

[2, -1, 1] 

sage: ieq[0] + vector(ieq[1:]) * vector([1, -1]) 

0 

""" 

v = vector(v) 

c = -v * v 

if c < 0: 

c, v = -c, -v 

return [c] + list(v) 

 

def jacobi(M): 

r""" 

Compute the upper-triangular part of the Cholesky/Jacobi 

decomposition of the symmetric matrix ``M``. 

 

Let `M` be a symmetric `n \times n`-matrix over a field `F`. 

Let `m_{i,j}` denote the `(i,j)`-th entry of `M` for any 

`1 \leq i \leq n` and `1 \leq j \leq n`. Then, the 

upper-triangular part computed by this method is the 

upper-triangular `n \times n`-matrix `Q` whose 

`(i,j)`-th entry `q_{i,j}` satisfies 

 

.. MATH:: 

 

q_{i,j} = 

\begin{cases} 

\frac{1}{q_{i,i}} \left( m_{i,j} - \sum_{r<i} q_{r,r} q_{r,i} q_{r,j} \right) & i < j, \\ 

a_{i,j} - \sum_{r<i} q_{r,r} q_{r,i}^2 & i = j, \\ 

0 & i > j, 

\end{cases} 

 

for all `1 \leq i \leq n` and `1 \leq j \leq n`. (These 

equalities determine the entries of `Q` uniquely by 

recursion.) This matrix `Q` is defined for all `M` in a 

certain Zariski-dense open subset of the set of all 

`n \times n`-matrices. 

 

EXAMPLES:: 

 

sage: from sage.modules.diamond_cutting import jacobi 

sage: jacobi(identity_matrix(3) * 4) 

[4 0 0] 

[0 4 0] 

[0 0 4] 

 

sage: def testall(M): 

....: Q = jacobi(M) 

....: for j in range(3): 

....: for i in range(j): 

....: if Q[i,j] * Q[i,i] != M[i,j] - sum(Q[r,i] * Q[r,j] * Q[r,r] for r in range(i)): 

....: return False 

....: for i in range(3): 

....: if Q[i,i] != M[i,i] - sum(Q[r,i] ** 2 * Q[r,r] for r in range(i)): 

....: return False 

....: for j in range(i): 

....: if Q[i,j] != 0: 

....: return False 

....: return True 

 

sage: M = Matrix(QQ, [[8,1,5], [1,6,0], [5,0,3]]) 

sage: Q = jacobi(M); Q 

[ 8 1/8 5/8] 

[ 0 47/8 -5/47] 

[ 0 0 -9/47] 

sage: testall(M) 

True 

 

sage: M = Matrix(QQ, [[3,6,-1,7],[6,9,8,5],[-1,8,2,4],[7,5,4,0]]) 

sage: testall(M) 

True 

""" 

dim = M.dimensions() 

if dim[0] != dim[1]: 

raise ValueError("the matrix must be square") 

dim = dim[0] 

q = [list(row) for row in M] 

for i in range(dim - 1): 

for j in range(i + 1, dim): 

q[j][i] = q[i][j] 

q[i][j] = q[i][j] / q[i][i] 

for k in range(i + 1, dim): 

for l in range(k, dim): 

q[k][l] -= q[k][i] * q[i][l] 

for i in range(1, dim): 

for j in range(i): 

q[i][j] = 0 

return matrix(q) 

 

def diamond_cut(V, GM, C, verbose=False): 

r""" 

Perform diamond cutting on polyhedron ``V`` with basis matrix ``GM`` 

and radius ``C``. 

 

INPUT: 

 

- ``V`` -- polyhedron to cut from 

 

- ``GM`` -- half of the basis matrix of the lattice 

 

- ``C`` -- radius to use in cutting algorithm 

 

- ``verbose`` -- (default: ``False``) whether to print debug information 

 

OUTPUT: 

 

A :class:`Polyhedron` instance. 

 

EXAMPLES:: 

 

sage: from sage.modules.diamond_cutting import diamond_cut 

sage: V = Polyhedron([[0], [2]]) 

sage: GM = matrix([2]) 

sage: V = diamond_cut(V, GM, 4) 

sage: V.vertices() 

(A vertex at (2), A vertex at (0)) 

""" 

# coerce to floats 

GM = GM.n() 

C = float(C) 

if verbose: 

print("Cut\n{}\nwith radius {}".format(GM, C)) 

 

dim = GM.dimensions() 

if dim[0] != dim[1]: 

raise ValueError("the matrix must be square") 

dim = dim[0] 

T = [0] * dim 

U = [0] * dim 

x = [0] * dim 

L = [0] * dim 

 

# calculate the Gram matrix 

q = matrix([[sum(GM[i][k] * GM[j][k] for k in range(dim)) for j in range(dim)] for i in range(dim)]) 

if verbose: 

print( "q:\n{}".format(q.n()) ) 

# apply Cholesky/Jacobi decomposition 

q = jacobi(q) 

if verbose: 

print( "q:\n{}".format(q.n()) ) 

 

i = dim - 1 

T[i] = C 

U[i] = 0 

 

new_dimension = True 

cut_count = 0 

inequalities = [] 

while True: 

if verbose: 

print("Dimension: {}".format(i)) 

if new_dimension: 

Z = sqrt(T[i] / q[i][i]) 

if verbose: 

print("Z: {}".format(Z)) 

L[i] = int(floor(Z - U[i])) 

if verbose: 

print("L: {}".format(L)) 

x[i] = int(ceil(-Z - U[i]) - 1) 

new_dimension = False 

 

x[i] += 1 

if verbose: 

print("x: {}".format(x)) 

if x[i] > L[i]: 

i += 1 

elif i > 0: 

T[i - 1] = T[i] - q[i][i] * (x[i] + U[i]) ** 2 

i -= 1 

U[i] = 0 

for j in range(i + 1, dim): 

U[i] += q[i][j] * x[j] 

new_dimension = True 

else: 

if all(elmt == 0 for elmt in x): 

break 

hv = [0] * dim 

for k in range(dim): 

for j in range(dim): 

hv[k] += x[j] * GM[j][k] 

hv = vector(hv) 

 

for hv in [hv, -hv]: 

cut_count += 1 

if verbose: 

print("\n%d) Cut using normal vector %s" % (cut_count, hv)) 

hv = [QQ(round(elmt, 6)) for elmt in hv] 

inequalities.append(plane_inequality(hv)) 

#cut = Polyhedron(ieqs=[plane_inequality(hv)]) 

#V = V.intersection(cut) 

 

if verbose: 

print("Final cut") 

cut = Polyhedron(ieqs=inequalities) 

V = V.intersection(cut) 

 

if verbose: 

print("End") 

 

return V 

 

def calculate_voronoi_cell(basis, radius=None, verbose=False): 

""" 

Calculate the Voronoi cell of the lattice defined by basis 

 

INPUT: 

 

- ``basis`` -- embedded basis matrix of the lattice 

 

- ``radius`` -- radius of basis vectors to consider 

 

- ``verbose`` -- whether to print debug information 

 

OUTPUT: 

 

A :class:`Polyhedron` instance. 

 

EXAMPLES:: 

 

sage: from sage.modules.diamond_cutting import calculate_voronoi_cell 

sage: V = calculate_voronoi_cell(matrix([[1, 0], [0, 1]])) 

sage: V.volume() 

1 

""" 

 

dim = basis.dimensions() 

artificial_length = None 

if dim[0] < dim[1]: 

# introduce "artificial" basis points (representing infinity) 

artificial_length = max(abs(v) for v in basis).ceil() * 2 

additional_vectors = identity_matrix(dim[1]) * artificial_length 

basis = basis.stack(additional_vectors) 

# LLL-reduce to get quadratic matrix 

basis = basis.LLL() 

basis = matrix([v for v in basis if v]) 

dim = basis.dimensions() 

if dim[0] != dim[1]: 

raise ValueError("invalid matrix") 

basis = basis / 2 

 

ieqs = [] 

for v in basis: 

ieqs.append(plane_inequality(v)) 

ieqs.append(plane_inequality(-v)) 

Q = Polyhedron(ieqs=ieqs) 

 

# twice the length of longest vertex in Q is a safe choice 

if radius is None: 

radius = 2 * max(abs(v) ** 2 for v in basis) 

 

V = diamond_cut(Q, basis, radius, verbose=verbose) 

 

if artificial_length is not None: 

# remove inequalities introduced by artificial basis points 

H = V.Hrepresentation() 

H = [v for v in H if all(not V._is_zero(v.A() * w / 2 - v.b() and 

not V._is_zero(v.A() * (-w) / 2 - v.b())) for w in additional_vectors)] 

V = Polyhedron(ieqs=H) 

 

return V