Coverage for /usr/lib/python3/dist-packages/sympy/combinatorics/tensor_can.py: 6%

1from sympy.combinatorics.permutations import Permutation, _af_rmul, \ 

2 _af_invert, _af_new 

3from sympy.combinatorics.perm_groups import PermutationGroup, _orbit, \ 

4 _orbit_transversal 

5from sympy.combinatorics.util import _distribute_gens_by_base, \ 

6 _orbits_transversals_from_bsgs 

7 

8""" 

9 References for tensor canonicalization: 

10 

11 [1] R. Portugal "Algorithmic simplification of tensor expressions", 

12 J. Phys. A 32 (1999) 7779-7789 

13 

14 [2] R. Portugal, B.F. Svaiter "Group-theoretic Approach for Symbolic 

15 Tensor Manipulation: I. Free Indices" 

16 arXiv:math-ph/0107031v1 

17 

18 [3] L.R.U. Manssur, R. Portugal "Group-theoretic Approach for Symbolic 

19 Tensor Manipulation: II. Dummy Indices" 

20 arXiv:math-ph/0107032v1 

21 

22 [4] xperm.c part of XPerm written by J. M. Martin-Garcia 

23 http://www.xact.es/index.html 

24""" 

25 

26 

27def dummy_sgs(dummies, sym, n): 

28 """ 

29 Return the strong generators for dummy indices. 

30 

31 Parameters 

32 ========== 

33 

34 dummies : List of dummy indices. 

35 `dummies[2k], dummies[2k+1]` are paired indices. 

36 In base form, the dummy indices are always in 

37 consecutive positions. 

38 sym : symmetry under interchange of contracted dummies:: 

39 * None no symmetry 

40 * 0 commuting 

41 * 1 anticommuting 

42 

43 n : number of indices 

44 

45 Examples 

46 ======== 

47 

48 >>> from sympy.combinatorics.tensor_can import dummy_sgs 

49 >>> dummy_sgs(list(range(2, 8)), 0, 8) 

50 [[0, 1, 3, 2, 4, 5, 6, 7, 8, 9], [0, 1, 2, 3, 5, 4, 6, 7, 8, 9], 

51 [0, 1, 2, 3, 4, 5, 7, 6, 8, 9], [0, 1, 4, 5, 2, 3, 6, 7, 8, 9], 

52 [0, 1, 2, 3, 6, 7, 4, 5, 8, 9]] 

53 """ 

54 if len(dummies) > n: 

55 raise ValueError("List too large") 

56 res = [] 

57 # exchange of contravariant and covariant indices 

58 if sym is not None: 

59 for j in dummies[::2]: 

60 a = list(range(n + 2)) 

61 if sym == 1: 

62 a[n] = n + 1 

63 a[n + 1] = n 

64 a[j], a[j + 1] = a[j + 1], a[j] 

65 res.append(a) 

66 # rename dummy indices 

67 for j in dummies[:-3:2]: 

68 a = list(range(n + 2)) 

69 a[j:j + 4] = a[j + 2], a[j + 3], a[j], a[j + 1] 

70 res.append(a) 

71 return res 

72 

73 

74def _min_dummies(dummies, sym, indices): 

75 """ 

76 Return list of minima of the orbits of indices in group of dummies. 

77 See ``double_coset_can_rep`` for the description of ``dummies`` and ``sym``. 

78 ``indices`` is the initial list of dummy indices. 

79 

80 Examples 

81 ======== 

82 

83 >>> from sympy.combinatorics.tensor_can import _min_dummies 

84 >>> _min_dummies([list(range(2, 8))], [0], list(range(10))) 

85 [0, 1, 2, 2, 2, 2, 2, 2, 8, 9] 

86 """ 

87 num_types = len(sym) 

88 m = [min(dx) if dx else None for dx in dummies] 

89 res = indices[:] 

90 for i in range(num_types): 

91 for c, i in enumerate(indices): 

92 for j in range(num_types): 

93 if i in dummies[j]: 

94 res[c] = m[j] 

95 break 

96 return res 

97 

98 

99def _trace_S(s, j, b, S_cosets): 

100 """ 

101 Return the representative h satisfying s[h[b]] == j 

102 

103 If there is not such a representative return None 

104 """ 

105 for h in S_cosets[b]: 

106 if s[h[b]] == j: 

107 return h 

108 return None 

109 

110 

111def _trace_D(gj, p_i, Dxtrav): 

112 """ 

113 Return the representative h satisfying h[gj] == p_i 

114 

115 If there is not such a representative return None 

116 """ 

117 for h in Dxtrav: 

118 if h[gj] == p_i: 

119 return h 

120 return None 

121 

122 

123def _dumx_remove(dumx, dumx_flat, p0): 

124 """ 

125 remove p0 from dumx 

126 """ 

127 res = [] 

128 for dx in dumx: 

129 if p0 not in dx: 

130 res.append(dx) 

131 continue 

132 k = dx.index(p0) 

133 if k % 2 == 0: 

134 p0_paired = dx[k + 1] 

135 else: 

136 p0_paired = dx[k - 1] 

137 dx.remove(p0) 

138 dx.remove(p0_paired) 

139 dumx_flat.remove(p0) 

140 dumx_flat.remove(p0_paired) 

141 res.append(dx) 

142 

143 

144def transversal2coset(size, base, transversal): 

145 a = [] 

146 j = 0 

147 for i in range(size): 

148 if i in base: 

149 a.append(sorted(transversal[j].values())) 

150 j += 1 

151 else: 

152 a.append([list(range(size))]) 

153 j = len(a) - 1 

154 while a[j] == [list(range(size))]: 

155 j -= 1 

156 return a[:j + 1] 

157 

158 

159def double_coset_can_rep(dummies, sym, b_S, sgens, S_transversals, g): 

160 r""" 

161 Butler-Portugal algorithm for tensor canonicalization with dummy indices. 

162 

163 Parameters 

164 ========== 

165 

166 dummies 

167 list of lists of dummy indices, 

168 one list for each type of index; 

169 the dummy indices are put in order contravariant, covariant 

170 [d0, -d0, d1, -d1, ...]. 

171 

172 sym 

173 list of the symmetries of the index metric for each type. 

174 

175 possible symmetries of the metrics 

176 * 0 symmetric 

177 * 1 antisymmetric 

178 * None no symmetry 

179 

180 b_S 

181 base of a minimal slot symmetry BSGS. 

182 

183 sgens 

184 generators of the slot symmetry BSGS. 

185 

186 S_transversals 

187 transversals for the slot BSGS. 

188 

189 g 

190 permutation representing the tensor. 

191 

192 Returns 

193 ======= 

194 

195 Return 0 if the tensor is zero, else return the array form of 

196 the permutation representing the canonical form of the tensor. 

197 

198 Notes 

199 ===== 

200 

201 A tensor with dummy indices can be represented in a number 

202 of equivalent ways which typically grows exponentially with 

203 the number of indices. To be able to establish if two tensors 

204 with many indices are equal becomes computationally very slow 

205 in absence of an efficient algorithm. 

206 

207 The Butler-Portugal algorithm [3] is an efficient algorithm to 

208 put tensors in canonical form, solving the above problem. 

209 

210 Portugal observed that a tensor can be represented by a permutation, 

211 and that the class of tensors equivalent to it under slot and dummy 

212 symmetries is equivalent to the double coset `D*g*S` 

213 (Note: in this documentation we use the conventions for multiplication 

214 of permutations p, q with (p*q)(i) = p[q[i]] which is opposite 

215 to the one used in the Permutation class) 

216 

217 Using the algorithm by Butler to find a representative of the 

218 double coset one can find a canonical form for the tensor. 

219 

220 To see this correspondence, 

221 let `g` be a permutation in array form; a tensor with indices `ind` 

222 (the indices including both the contravariant and the covariant ones) 

223 can be written as 

224 

225 `t = T(ind[g[0]], \dots, ind[g[n-1]])`, 

226 

227 where `n = len(ind)`; 

228 `g` has size `n + 2`, the last two indices for the sign of the tensor 

229 (trick introduced in [4]). 

230 

231 A slot symmetry transformation `s` is a permutation acting on the slots 

232 `t \rightarrow T(ind[(g*s)[0]], \dots, ind[(g*s)[n-1]])` 

233 

234 A dummy symmetry transformation acts on `ind` 

235 `t \rightarrow T(ind[(d*g)[0]], \dots, ind[(d*g)[n-1]])` 

236 

237 Being interested only in the transformations of the tensor under 

238 these symmetries, one can represent the tensor by `g`, which transforms 

239 as 

240 

241 `g -> d*g*s`, so it belongs to the coset `D*g*S`, or in other words 

242 to the set of all permutations allowed by the slot and dummy symmetries. 

243 

244 Let us explain the conventions by an example. 

245 

246 Given a tensor `T^{d3 d2 d1}{}_{d1 d2 d3}` with the slot symmetries 

247 `T^{a0 a1 a2 a3 a4 a5} = -T^{a2 a1 a0 a3 a4 a5}` 

248 

249 `T^{a0 a1 a2 a3 a4 a5} = -T^{a4 a1 a2 a3 a0 a5}` 

250 

251 and symmetric metric, find the tensor equivalent to it which 

252 is the lowest under the ordering of indices: 

253 lexicographic ordering `d1, d2, d3` and then contravariant 

254 before covariant index; that is the canonical form of the tensor. 

255 

256 The canonical form is `-T^{d1 d2 d3}{}_{d1 d2 d3}` 

257 obtained using `T^{a0 a1 a2 a3 a4 a5} = -T^{a2 a1 a0 a3 a4 a5}`. 

258 

259 To convert this problem in the input for this function, 

260 use the following ordering of the index names 

261 (- for covariant for short) `d1, -d1, d2, -d2, d3, -d3` 

262 

263 `T^{d3 d2 d1}{}_{d1 d2 d3}` corresponds to `g = [4, 2, 0, 1, 3, 5, 6, 7]` 

264 where the last two indices are for the sign 

265 

266 `sgens = [Permutation(0, 2)(6, 7), Permutation(0, 4)(6, 7)]` 

267 

268 sgens[0] is the slot symmetry `-(0, 2)` 

269 `T^{a0 a1 a2 a3 a4 a5} = -T^{a2 a1 a0 a3 a4 a5}` 

270 

271 sgens[1] is the slot symmetry `-(0, 4)` 

272 `T^{a0 a1 a2 a3 a4 a5} = -T^{a4 a1 a2 a3 a0 a5}` 

273 

274 The dummy symmetry group D is generated by the strong base generators 

275 `[(0, 1), (2, 3), (4, 5), (0, 2)(1, 3), (0, 4)(1, 5)]` 

276 where the first three interchange covariant and contravariant 

277 positions of the same index (d1 <-> -d1) and the last two interchange 

278 the dummy indices themselves (d1 <-> d2). 

279 

280 The dummy symmetry acts from the left 

281 `d = [1, 0, 2, 3, 4, 5, 6, 7]` exchange `d1 \leftrightarrow -d1` 

282 `T^{d3 d2 d1}{}_{d1 d2 d3} == T^{d3 d2}{}_{d1}{}^{d1}{}_{d2 d3}` 

283 

284 `g=[4, 2, 0, 1, 3, 5, 6, 7] -> [4, 2, 1, 0, 3, 5, 6, 7] = _af_rmul(d, g)` 

285 which differs from `_af_rmul(g, d)`. 

286 

287 The slot symmetry acts from the right 

288 `s = [2, 1, 0, 3, 4, 5, 7, 6]` exchanges slots 0 and 2 and changes sign 

289 `T^{d3 d2 d1}{}_{d1 d2 d3} == -T^{d1 d2 d3}{}_{d1 d2 d3}` 

290 

291 `g=[4,2,0,1,3,5,6,7] -> [0, 2, 4, 1, 3, 5, 7, 6] = _af_rmul(g, s)` 

292 

293 Example in which the tensor is zero, same slot symmetries as above: 

294 `T^{d2}{}_{d1 d3}{}^{d1 d3}{}_{d2}` 

295 

296 `= -T^{d3}{}_{d1 d3}{}^{d1 d2}{}_{d2}` under slot symmetry `-(0,4)`; 

297 

298 `= T_{d3 d1}{}^{d3}{}^{d1 d2}{}_{d2}` under slot symmetry `-(0,2)`; 

299 

300 `= T^{d3}{}_{d1 d3}{}^{d1 d2}{}_{d2}` symmetric metric; 

301 

302 `= 0` since two of these lines have tensors differ only for the sign. 

303 

304 The double coset D*g*S consists of permutations `h = d*g*s` corresponding 

305 to equivalent tensors; if there are two `h` which are the same apart 

306 from the sign, return zero; otherwise 

307 choose as representative the tensor with indices 

308 ordered lexicographically according to `[d1, -d1, d2, -d2, d3, -d3]` 

309 that is ``rep = min(D*g*S) = min([d*g*s for d in D for s in S])`` 

310 

311 The indices are fixed one by one; first choose the lowest index 

312 for slot 0, then the lowest remaining index for slot 1, etc. 

313 Doing this one obtains a chain of stabilizers 

314 

315 `S \rightarrow S_{b0} \rightarrow S_{b0,b1} \rightarrow \dots` and 

316 `D \rightarrow D_{p0} \rightarrow D_{p0,p1} \rightarrow \dots` 

317 

318 where ``[b0, b1, ...] = range(b)`` is a base of the symmetric group; 

319 the strong base `b_S` of S is an ordered sublist of it; 

320 therefore it is sufficient to compute once the 

321 strong base generators of S using the Schreier-Sims algorithm; 

322 the stabilizers of the strong base generators are the 

323 strong base generators of the stabilizer subgroup. 

324 

325 ``dbase = [p0, p1, ...]`` is not in general in lexicographic order, 

326 so that one must recompute the strong base generators each time; 

327 however this is trivial, there is no need to use the Schreier-Sims 

328 algorithm for D. 

329 

330 The algorithm keeps a TAB of elements `(s_i, d_i, h_i)` 

331 where `h_i = d_i \times g \times s_i` satisfying `h_i[j] = p_j` for `0 \le j < i` 

332 starting from `s_0 = id, d_0 = id, h_0 = g`. 

333 

334 The equations `h_0[0] = p_0, h_1[1] = p_1, \dots` are solved in this order, 

335 choosing each time the lowest possible value of p_i 

336 

337 For `j < i` 

338 `d_i*g*s_i*S_{b_0, \dots, b_{i-1}}*b_j = D_{p_0, \dots, p_{i-1}}*p_j` 

339 so that for dx in `D_{p_0,\dots,p_{i-1}}` and sx in 

340 `S_{base[0], \dots, base[i-1]}` one has `dx*d_i*g*s_i*sx*b_j = p_j` 

341 

342 Search for dx, sx such that this equation holds for `j = i`; 

343 it can be written as `s_i*sx*b_j = J, dx*d_i*g*J = p_j` 

344 `sx*b_j = s_i**-1*J; sx = trace(s_i**-1, S_{b_0,...,b_{i-1}})` 

345 `dx**-1*p_j = d_i*g*J; dx = trace(d_i*g*J, D_{p_0,...,p_{i-1}})` 

346 

347 `s_{i+1} = s_i*trace(s_i**-1*J, S_{b_0,...,b_{i-1}})` 

348 `d_{i+1} = trace(d_i*g*J, D_{p_0,...,p_{i-1}})**-1*d_i` 

349 `h_{i+1}*b_i = d_{i+1}*g*s_{i+1}*b_i = p_i` 

350 

351 `h_n*b_j = p_j` for all j, so that `h_n` is the solution. 

352 

353 Add the found `(s, d, h)` to TAB1. 

354 

355 At the end of the iteration sort TAB1 with respect to the `h`; 

356 if there are two consecutive `h` in TAB1 which differ only for the 

357 sign, the tensor is zero, so return 0; 

358 if there are two consecutive `h` which are equal, keep only one. 

359 

360 Then stabilize the slot generators under `i` and the dummy generators 

361 under `p_i`. 

362 

363 Assign `TAB = TAB1` at the end of the iteration step. 

364 

365 At the end `TAB` contains a unique `(s, d, h)`, since all the slots 

366 of the tensor `h` have been fixed to have the minimum value according 

367 to the symmetries. The algorithm returns `h`. 

368 

369 It is important that the slot BSGS has lexicographic minimal base, 

370 otherwise there is an `i` which does not belong to the slot base 

371 for which `p_i` is fixed by the dummy symmetry only, while `i` 

372 is not invariant from the slot stabilizer, so `p_i` is not in 

373 general the minimal value. 

374 

375 This algorithm differs slightly from the original algorithm [3]: 

376 the canonical form is minimal lexicographically, and 

377 the BSGS has minimal base under lexicographic order. 

378 Equal tensors `h` are eliminated from TAB. 

379 

380 

381 Examples 

382 ======== 

383 

384 >>> from sympy.combinatorics.permutations import Permutation 

385 >>> from sympy.combinatorics.tensor_can import double_coset_can_rep, get_transversals 

386 >>> gens = [Permutation(x) for x in [[2, 1, 0, 3, 4, 5, 7, 6], [4, 1, 2, 3, 0, 5, 7, 6]]] 

387 >>> base = [0, 2] 

388 >>> g = Permutation([4, 2, 0, 1, 3, 5, 6, 7]) 

389 >>> transversals = get_transversals(base, gens) 

390 >>> double_coset_can_rep([list(range(6))], [0], base, gens, transversals, g) 

391 [0, 1, 2, 3, 4, 5, 7, 6] 

392 

393 >>> g = Permutation([4, 1, 3, 0, 5, 2, 6, 7]) 

394 >>> double_coset_can_rep([list(range(6))], [0], base, gens, transversals, g) 

395 0 

396 """ 

397 size = g.size 

398 g = g.array_form 

399 num_dummies = size - 2 

400 indices = list(range(num_dummies)) 

401 all_metrics_with_sym = not any(_ is None for _ in sym) 

402 num_types = len(sym) 

403 dumx = dummies[:] 

404 dumx_flat = [] 

405 for dx in dumx: 

406 dumx_flat.extend(dx) 

407 b_S = b_S[:] 

408 sgensx = [h._array_form for h in sgens] 

409 if b_S: 

410 S_transversals = transversal2coset(size, b_S, S_transversals) 

411 # strong generating set for D 

412 dsgsx = [] 

413 for i in range(num_types): 

414 dsgsx.extend(dummy_sgs(dumx[i], sym[i], num_dummies)) 

415 idn = list(range(size)) 

416 # TAB = list of entries (s, d, h) where h = _af_rmuln(d,g,s) 

417 # for short, in the following d*g*s means _af_rmuln(d,g,s) 

418 TAB = [(idn, idn, g)] 

419 for i in range(size - 2): 

420 b = i 

421 testb = b in b_S and sgensx 

422 if testb: 

423 sgensx1 = [_af_new(_) for _ in sgensx] 

424 deltab = _orbit(size, sgensx1, b) 

425 else: 

426 deltab = {b} 

427 # p1 = min(IMAGES) = min(Union D_p*h*deltab for h in TAB) 

428 if all_metrics_with_sym: 

429 md = _min_dummies(dumx, sym, indices) 

430 else: 

431 md = [min(_orbit(size, [_af_new( 

432 ddx) for ddx in dsgsx], ii)) for ii in range(size - 2)] 

433 

434 p_i = min([min([md[h[x]] for x in deltab]) for s, d, h in TAB]) 

435 dsgsx1 = [_af_new(_) for _ in dsgsx] 

436 Dxtrav = _orbit_transversal(size, dsgsx1, p_i, False, af=True) \ 

437 if dsgsx else None 

438 if Dxtrav: 

439 Dxtrav = [_af_invert(x) for x in Dxtrav] 

440 # compute the orbit of p_i 

441 for ii in range(num_types): 

442 if p_i in dumx[ii]: 

443 # the orbit is made by all the indices in dum[ii] 

444 if sym[ii] is not None: 

445 deltap = dumx[ii] 

446 else: 

447 # the orbit is made by all the even indices if p_i 

448 # is even, by all the odd indices if p_i is odd 

449 p_i_index = dumx[ii].index(p_i) % 2 

450 deltap = dumx[ii][p_i_index::2] 

451 break 

452 else: 

453 deltap = [p_i] 

454 TAB1 = [] 

455 while TAB: 

456 s, d, h = TAB.pop() 

457 if min([md[h[x]] for x in deltab]) != p_i: 

458 continue 

459 deltab1 = [x for x in deltab if md[h[x]] == p_i] 

460 # NEXT = s*deltab1 intersection (d*g)**-1*deltap 

461 dg = _af_rmul(d, g) 

462 dginv = _af_invert(dg) 

463 sdeltab = [s[x] for x in deltab1] 

464 gdeltap = [dginv[x] for x in deltap] 

465 NEXT = [x for x in sdeltab if x in gdeltap] 

466 # d, s satisfy 

467 # d*g*s*base[i-1] = p_{i-1}; using the stabilizers 

468 # d*g*s*S_{base[0],...,base[i-1]}*base[i-1] = 

469 # D_{p_0,...,p_{i-1}}*p_{i-1} 

470 # so that to find d1, s1 satisfying d1*g*s1*b = p_i 

471 # one can look for dx in D_{p_0,...,p_{i-1}} and 

472 # sx in S_{base[0],...,base[i-1]} 

473 # d1 = dx*d; s1 = s*sx 

474 # d1*g*s1*b = dx*d*g*s*sx*b = p_i 

475 for j in NEXT: 

476 if testb: 

477 # solve s1*b = j with s1 = s*sx for some element sx 

478 # of the stabilizer of ..., base[i-1] 

479 # sx*b = s**-1*j; sx = _trace_S(s, j,...) 

480 # s1 = s*trace_S(s**-1*j,...) 

481 s1 = _trace_S(s, j, b, S_transversals) 

482 if not s1: 

483 continue 

484 else: 

485 s1 = [s[ix] for ix in s1] 

486 else: 

487 s1 = s 

488 # assert s1[b] == j # invariant 

489 # solve d1*g*j = p_i with d1 = dx*d for some element dg 

490 # of the stabilizer of ..., p_{i-1} 

491 # dx**-1*p_i = d*g*j; dx**-1 = trace_D(d*g*j,...) 

492 # d1 = trace_D(d*g*j,...)**-1*d 

493 # to save an inversion in the inner loop; notice we did 

494 # Dxtrav = [perm_af_invert(x) for x in Dxtrav] out of the loop 

495 if Dxtrav: 

496 d1 = _trace_D(dg[j], p_i, Dxtrav) 

497 if not d1: 

498 continue 

499 else: 

500 if p_i != dg[j]: 

501 continue 

502 d1 = idn 

503 assert d1[dg[j]] == p_i # invariant 

504 d1 = [d1[ix] for ix in d] 

505 h1 = [d1[g[ix]] for ix in s1] 

506 # assert h1[b] == p_i # invariant 

507 TAB1.append((s1, d1, h1)) 

508 

509 # if TAB contains equal permutations, keep only one of them; 

510 # if TAB contains equal permutations up to the sign, return 0 

511 TAB1.sort(key=lambda x: x[-1]) 

512 prev = [0] * size 

513 while TAB1: 

514 s, d, h = TAB1.pop() 

515 if h[:-2] == prev[:-2]: 

516 if h[-1] != prev[-1]: 

517 return 0 

518 else: 

519 TAB.append((s, d, h)) 

520 prev = h 

521 

522 # stabilize the SGS 

523 sgensx = [h for h in sgensx if h[b] == b] 

524 if b in b_S: 

525 b_S.remove(b) 

526 _dumx_remove(dumx, dumx_flat, p_i) 

527 dsgsx = [] 

528 for i in range(num_types): 

529 dsgsx.extend(dummy_sgs(dumx[i], sym[i], num_dummies)) 

530 return TAB[0][-1] 

531 

532 

533def canonical_free(base, gens, g, num_free): 

534 """ 

535 Canonicalization of a tensor with respect to free indices 

536 choosing the minimum with respect to lexicographical ordering 

537 in the free indices. 

538 

539 Explanation 

540 =========== 

541 

542 ``base``, ``gens`` BSGS for slot permutation group 

543 ``g`` permutation representing the tensor 

544 ``num_free`` number of free indices 

545 The indices must be ordered with first the free indices 

546 

547 See explanation in double_coset_can_rep 

548 The algorithm is a variation of the one given in [2]. 

549 

550 Examples 

551 ======== 

552 

553 >>> from sympy.combinatorics import Permutation 

554 >>> from sympy.combinatorics.tensor_can import canonical_free 

555 >>> gens = [[1, 0, 2, 3, 5, 4], [2, 3, 0, 1, 4, 5],[0, 1, 3, 2, 5, 4]] 

556 >>> gens = [Permutation(h) for h in gens] 

557 >>> base = [0, 2] 

558 >>> g = Permutation([2, 1, 0, 3, 4, 5]) 

559 >>> canonical_free(base, gens, g, 4) 

560 [0, 3, 1, 2, 5, 4] 

561 

562 Consider the product of Riemann tensors 

563 ``T = R^{a}_{d0}^{d1,d2}*R_{d2,d1}^{d0,b}`` 

564 The order of the indices is ``[a, b, d0, -d0, d1, -d1, d2, -d2]`` 

565 The permutation corresponding to the tensor is 

566 ``g = [0, 3, 4, 6, 7, 5, 2, 1, 8, 9]`` 

567 

568 In particular ``a`` is position ``0``, ``b`` is in position ``9``. 

569 Use the slot symmetries to get `T` is a form which is the minimal 

570 in lexicographic order in the free indices ``a`` and ``b``, e.g. 

571 ``-R^{a}_{d0}^{d1,d2}*R^{b,d0}_{d2,d1}`` corresponding to 

572 ``[0, 3, 4, 6, 1, 2, 7, 5, 9, 8]`` 

573 

574 >>> from sympy.combinatorics.tensor_can import riemann_bsgs, tensor_gens 

575 >>> base, gens = riemann_bsgs 

576 >>> size, sbase, sgens = tensor_gens(base, gens, [[], []], 0) 

577 >>> g = Permutation([0, 3, 4, 6, 7, 5, 2, 1, 8, 9]) 

578 >>> canonical_free(sbase, [Permutation(h) for h in sgens], g, 2) 

579 [0, 3, 4, 6, 1, 2, 7, 5, 9, 8] 

580 """ 

581 g = g.array_form 

582 size = len(g) 

583 if not base: 

584 return g[:] 

585 

586 transversals = get_transversals(base, gens) 

587 for x in sorted(g[:-2]): 

588 if x not in base: 

589 base.append(x) 

590 h = g 

591 for i, transv in enumerate(transversals): 

592 h_i = [size]*num_free 

593 # find the element s in transversals[i] such that 

594 # _af_rmul(h, s) has its free elements with the lowest position in h 

595 s = None 

596 for sk in transv.values(): 

597 h1 = _af_rmul(h, sk) 

598 hi = [h1.index(ix) for ix in range(num_free)] 

599 if hi < h_i: 

600 h_i = hi 

601 s = sk 

602 if s: 

603 h = _af_rmul(h, s) 

604 return h 

605 

606 

607def _get_map_slots(size, fixed_slots): 

608 res = list(range(size)) 

609 pos = 0 

610 for i in range(size): 

611 if i in fixed_slots: 

612 continue 

613 res[i] = pos 

614 pos += 1 

615 return res 

616 

617 

618def _lift_sgens(size, fixed_slots, free, s): 

619 a = [] 

620 j = k = 0 

621 fd = list(zip(fixed_slots, free)) 

622 fd = [y for x, y in sorted(fd)] 

623 num_free = len(free) 

624 for i in range(size): 

625 if i in fixed_slots: 

626 a.append(fd[k]) 

627 k += 1 

628 else: 

629 a.append(s[j] + num_free) 

630 j += 1 

631 return a 

632 

633 

634def canonicalize(g, dummies, msym, *v): 

635 """ 

636 canonicalize tensor formed by tensors 

637 

638 Parameters 

639 ========== 

640 

641 g : permutation representing the tensor 

642 

643 dummies : list representing the dummy indices 

644 it can be a list of dummy indices of the same type 

645 or a list of lists of dummy indices, one list for each 

646 type of index; 

647 the dummy indices must come after the free indices, 

648 and put in order contravariant, covariant 

649 [d0, -d0, d1,-d1,...] 

650 

651 msym : symmetry of the metric(s) 

652 it can be an integer or a list; 

653 in the first case it is the symmetry of the dummy index metric; 

654 in the second case it is the list of the symmetries of the 

655 index metric for each type 

656 

657 v : list, (base_i, gens_i, n_i, sym_i) for tensors of type `i` 

658 

659 base_i, gens_i : BSGS for tensors of this type. 

660 The BSGS should have minimal base under lexicographic ordering; 

661 if not, an attempt is made do get the minimal BSGS; 

662 in case of failure, 

663 canonicalize_naive is used, which is much slower. 

664 

665 n_i : number of tensors of type `i`. 

666 

667 sym_i : symmetry under exchange of component tensors of type `i`. 

668 

669 Both for msym and sym_i the cases are 

670 * None no symmetry 

671 * 0 commuting 

672 * 1 anticommuting 

673 

674 Returns 

675 ======= 

676 

677 0 if the tensor is zero, else return the array form of 

678 the permutation representing the canonical form of the tensor. 

679 

680 Algorithm 

681 ========= 

682 

683 First one uses canonical_free to get the minimum tensor under 

684 lexicographic order, using only the slot symmetries. 

685 If the component tensors have not minimal BSGS, it is attempted 

686 to find it; if the attempt fails canonicalize_naive 

687 is used instead. 

688 

689 Compute the residual slot symmetry keeping fixed the free indices 

690 using tensor_gens(base, gens, list_free_indices, sym). 

691 

692 Reduce the problem eliminating the free indices. 

693 

694 Then use double_coset_can_rep and lift back the result reintroducing 

695 the free indices. 

696 

697 Examples 

698 ======== 

699 

700 one type of index with commuting metric; 

701 

702 `A_{a b}` and `B_{a b}` antisymmetric and commuting 

703 

704 `T = A_{d0 d1} * B^{d0}{}_{d2} * B^{d2 d1}` 

705 

706 `ord = [d0,-d0,d1,-d1,d2,-d2]` order of the indices 

707 

708 g = [1, 3, 0, 5, 4, 2, 6, 7] 

709 

710 `T_c = 0` 

711 

712 >>> from sympy.combinatorics.tensor_can import get_symmetric_group_sgs, canonicalize, bsgs_direct_product 

713 >>> from sympy.combinatorics import Permutation 

714 >>> base2a, gens2a = get_symmetric_group_sgs(2, 1) 

715 >>> t0 = (base2a, gens2a, 1, 0) 

716 >>> t1 = (base2a, gens2a, 2, 0) 

717 >>> g = Permutation([1, 3, 0, 5, 4, 2, 6, 7]) 

718 >>> canonicalize(g, range(6), 0, t0, t1) 

719 0 

720 

721 same as above, but with `B_{a b}` anticommuting 

722 

723 `T_c = -A^{d0 d1} * B_{d0}{}^{d2} * B_{d1 d2}` 

724 

725 can = [0,2,1,4,3,5,7,6] 

726 

727 >>> t1 = (base2a, gens2a, 2, 1) 

728 >>> canonicalize(g, range(6), 0, t0, t1) 

729 [0, 2, 1, 4, 3, 5, 7, 6] 

730 

731 two types of indices `[a,b,c,d,e,f]` and `[m,n]`, in this order, 

732 both with commuting metric 

733 

734 `f^{a b c}` antisymmetric, commuting 

735 

736 `A_{m a}` no symmetry, commuting 

737 

738 `T = f^c{}_{d a} * f^f{}_{e b} * A_m{}^d * A^{m b} * A_n{}^a * A^{n e}` 

739 

740 ord = [c,f,a,-a,b,-b,d,-d,e,-e,m,-m,n,-n] 

741 

742 g = [0,7,3, 1,9,5, 11,6, 10,4, 13,2, 12,8, 14,15] 

743 

744 The canonical tensor is 

745 `T_c = -f^{c a b} * f^{f d e} * A^m{}_a * A_{m d} * A^n{}_b * A_{n e}` 

746 

747 can = [0,2,4, 1,6,8, 10,3, 11,7, 12,5, 13,9, 15,14] 

748 

749 >>> base_f, gens_f = get_symmetric_group_sgs(3, 1) 

750 >>> base1, gens1 = get_symmetric_group_sgs(1) 

751 >>> base_A, gens_A = bsgs_direct_product(base1, gens1, base1, gens1) 

752 >>> t0 = (base_f, gens_f, 2, 0) 

753 >>> t1 = (base_A, gens_A, 4, 0) 

754 >>> dummies = [range(2, 10), range(10, 14)] 

755 >>> g = Permutation([0, 7, 3, 1, 9, 5, 11, 6, 10, 4, 13, 2, 12, 8, 14, 15]) 

756 >>> canonicalize(g, dummies, [0, 0], t0, t1) 

757 [0, 2, 4, 1, 6, 8, 10, 3, 11, 7, 12, 5, 13, 9, 15, 14] 

758 """ 

759 from sympy.combinatorics.testutil import canonicalize_naive 

760 if not isinstance(msym, list): 

761 if msym not in (0, 1, None): 

762 raise ValueError('msym must be 0, 1 or None') 

763 num_types = 1 

764 else: 

765 num_types = len(msym) 

766 if not all(msymx in (0, 1, None) for msymx in msym): 

767 raise ValueError('msym entries must be 0, 1 or None') 

768 if len(dummies) != num_types: 

769 raise ValueError( 

770 'dummies and msym must have the same number of elements') 

771 size = g.size 

772 num_tensors = 0 

773 v1 = [] 

774 for base_i, gens_i, n_i, sym_i in v: 

775 # check that the BSGS is minimal; 

776 # this property is used in double_coset_can_rep; 

777 # if it is not minimal use canonicalize_naive 

778 if not _is_minimal_bsgs(base_i, gens_i): 

779 mbsgs = get_minimal_bsgs(base_i, gens_i) 

780 if not mbsgs: 

781 can = canonicalize_naive(g, dummies, msym, *v) 

782 return can 

783 base_i, gens_i = mbsgs 

784 v1.append((base_i, gens_i, [[]] * n_i, sym_i)) 

785 num_tensors += n_i 

786 

787 if num_types == 1 and not isinstance(msym, list): 

788 dummies = [dummies] 

789 msym = [msym] 

790 flat_dummies = [] 

791 for dumx in dummies: 

792 flat_dummies.extend(dumx) 

793 

794 if flat_dummies and flat_dummies != list(range(flat_dummies[0], flat_dummies[-1] + 1)): 

795 raise ValueError('dummies is not valid') 

796 

797 # slot symmetry of the tensor 

798 size1, sbase, sgens = gens_products(*v1) 

799 if size != size1: 

800 raise ValueError( 

801 'g has size %d, generators have size %d' % (size, size1)) 

802 free = [i for i in range(size - 2) if i not in flat_dummies] 

803 num_free = len(free) 

804 

805 # g1 minimal tensor under slot symmetry 

806 g1 = canonical_free(sbase, sgens, g, num_free) 

807 if not flat_dummies: 

808 return g1 

809 # save the sign of g1 

810 sign = 0 if g1[-1] == size - 1 else 1 

811 

812 # the free indices are kept fixed. 

813 # Determine free_i, the list of slots of tensors which are fixed 

814 # since they are occupied by free indices, which are fixed. 

815 start = 0 

816 for i, (base_i, gens_i, n_i, sym_i) in enumerate(v): 

817 free_i = [] 

818 len_tens = gens_i[0].size - 2 

819 # for each component tensor get a list od fixed islots 

820 for j in range(n_i): 

821 # get the elements corresponding to the component tensor 

822 h = g1[start:(start + len_tens)] 

823 fr = [] 

824 # get the positions of the fixed elements in h 

825 for k in free: 

826 if k in h: 

827 fr.append(h.index(k)) 

828 free_i.append(fr) 

829 start += len_tens 

830 v1[i] = (base_i, gens_i, free_i, sym_i) 

831 # BSGS of the tensor with fixed free indices 

832 # if tensor_gens fails in gens_product, use canonicalize_naive 

833 size, sbase, sgens = gens_products(*v1) 

834 

835 # reduce the permutations getting rid of the free indices 

836 pos_free = [g1.index(x) for x in range(num_free)] 

837 size_red = size - num_free 

838 g1_red = [x - num_free for x in g1 if x in flat_dummies] 

839 if sign: 

840 g1_red.extend([size_red - 1, size_red - 2]) 

841 else: 

842 g1_red.extend([size_red - 2, size_red - 1]) 

843 map_slots = _get_map_slots(size, pos_free) 

844 sbase_red = [map_slots[i] for i in sbase if i not in pos_free] 

845 sgens_red = [_af_new([map_slots[i] for i in y._array_form if i not in pos_free]) for y in sgens] 

846 dummies_red = [[x - num_free for x in y] for y in dummies] 

847 transv_red = get_transversals(sbase_red, sgens_red) 

848 g1_red = _af_new(g1_red) 

849 g2 = double_coset_can_rep( 

850 dummies_red, msym, sbase_red, sgens_red, transv_red, g1_red) 

851 if g2 == 0: 

852 return 0 

853 # lift to the case with the free indices 

854 g3 = _lift_sgens(size, pos_free, free, g2) 

855 return g3 

856 

857 

858def perm_af_direct_product(gens1, gens2, signed=True): 

859 """ 

860 Direct products of the generators gens1 and gens2. 

861 

862 Examples 

863 ======== 

864 

865 >>> from sympy.combinatorics.tensor_can import perm_af_direct_product 

866 >>> gens1 = [[1, 0, 2, 3], [0, 1, 3, 2]] 

867 >>> gens2 = [[1, 0]] 

868 >>> perm_af_direct_product(gens1, gens2, False) 

869 [[1, 0, 2, 3, 4, 5], [0, 1, 3, 2, 4, 5], [0, 1, 2, 3, 5, 4]] 

870 >>> gens1 = [[1, 0, 2, 3, 5, 4], [0, 1, 3, 2, 4, 5]] 

871 >>> gens2 = [[1, 0, 2, 3]] 

872 >>> perm_af_direct_product(gens1, gens2, True) 

873 [[1, 0, 2, 3, 4, 5, 7, 6], [0, 1, 3, 2, 4, 5, 6, 7], [0, 1, 2, 3, 5, 4, 6, 7]] 

874 """ 

875 gens1 = [list(x) for x in gens1] 

876 gens2 = [list(x) for x in gens2] 

877 s = 2 if signed else 0 

878 n1 = len(gens1[0]) - s 

879 n2 = len(gens2[0]) - s 

880 start = list(range(n1)) 

881 end = list(range(n1, n1 + n2)) 

882 if signed: 

883 gens1 = [gen[:-2] + end + [gen[-2] + n2, gen[-1] + n2] 

884 for gen in gens1] 

885 gens2 = [start + [x + n1 for x in gen] for gen in gens2] 

886 else: 

887 gens1 = [gen + end for gen in gens1] 

888 gens2 = [start + [x + n1 for x in gen] for gen in gens2] 

889 

890 res = gens1 + gens2 

891 

892 return res 

893 

894 

895def bsgs_direct_product(base1, gens1, base2, gens2, signed=True): 

896 """ 

897 Direct product of two BSGS. 

898 

899 Parameters 

900 ========== 

901 

902 base1 : base of the first BSGS. 

903 

904 gens1 : strong generating sequence of the first BSGS. 

905 

906 base2, gens2 : similarly for the second BSGS. 

907 

908 signed : flag for signed permutations. 

909 

910 Examples 

911 ======== 

912 

913 >>> from sympy.combinatorics.tensor_can import (get_symmetric_group_sgs, bsgs_direct_product) 

914 >>> base1, gens1 = get_symmetric_group_sgs(1) 

915 >>> base2, gens2 = get_symmetric_group_sgs(2) 

916 >>> bsgs_direct_product(base1, gens1, base2, gens2) 

917 ([1], [(4)(1 2)]) 

918 """ 

919 s = 2 if signed else 0 

920 n1 = gens1[0].size - s 

921 base = list(base1) 

922 base += [x + n1 for x in base2] 

923 gens1 = [h._array_form for h in gens1] 

924 gens2 = [h._array_form for h in gens2] 

925 gens = perm_af_direct_product(gens1, gens2, signed) 

926 size = len(gens[0]) 

927 id_af = list(range(size)) 

928 gens = [h for h in gens if h != id_af] 

929 if not gens: 

930 gens = [id_af] 

931 return base, [_af_new(h) for h in gens] 

932 

933 

934def get_symmetric_group_sgs(n, antisym=False): 

935 """ 

936 Return base, gens of the minimal BSGS for (anti)symmetric tensor 

937 

938 Parameters 

939 ========== 

940 

941 n : rank of the tensor 

942 antisym : bool 

943 ``antisym = False`` symmetric tensor 

944 ``antisym = True`` antisymmetric tensor 

945 

946 Examples 

947 ======== 

948 

949 >>> from sympy.combinatorics.tensor_can import get_symmetric_group_sgs 

950 >>> get_symmetric_group_sgs(3) 

951 ([0, 1], [(4)(0 1), (4)(1 2)]) 

952 """ 

953 if n == 1: 

954 return [], [_af_new(list(range(3)))] 

955 gens = [Permutation(n - 1)(i, i + 1)._array_form for i in range(n - 1)] 

956 if antisym == 0: 

957 gens = [x + [n, n + 1] for x in gens] 

958 else: 

959 gens = [x + [n + 1, n] for x in gens] 

960 base = list(range(n - 1)) 

961 return base, [_af_new(h) for h in gens] 

962 

963riemann_bsgs = [0, 2], [Permutation(0, 1)(4, 5), Permutation(2, 3)(4, 5),