Coverage for /home/runner/work/tket/tket/pytket/pytket/utils/symbolic.py: 69%

1# Copyright Quantinuum 

2# 

3# Licensed under the Apache License, Version 2.0 (the "License"); 

4# you may not use this file except in compliance with the License. 

5# You may obtain a copy of the License at 

6# 

7# http://www.apache.org/licenses/LICENSE-2.0 

8# 

9# Unless required by applicable law or agreed to in writing, software 

10# distributed under the License is distributed on an "AS IS" BASIS, 

11# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 

12# See the License for the specific language governing permissions and 

13# limitations under the License. 

14 

15"""Collection of methods to calculate symbolic statevectors and unitaries, 

16for symbolic circuits. This uses the sympy.physics.quantum module and produces 

17sympy objects. The implementations are slow and scale poorly, so this is 

18only suitable for very small (up to 5 qubit) circuits.""" 

19from collections.abc import Callable 

20from typing import cast 

21 

22import numpy as np 

23import sympy 

24from sympy import ( 

25 BlockDiagMatrix, 

26 BlockMatrix, 

27 Expr, 

28 I, 

29 Identity, 

30 ImmutableMatrix, 

31 Matrix, 

32 Mul, 

33 diag, 

34 eye, 

35 zeros, 

36) 

37from sympy.physics.quantum import gate as symgate 

38from sympy.physics.quantum import represent 

39from sympy.physics.quantum.qapply import qapply 

40from sympy.physics.quantum.qubit import Qubit, matrix_to_qubit 

41from sympy.physics.quantum.tensorproduct import matrix_tensor_product 

42 

43from pytket.circuit import Circuit, Op, OpType 

44 

45# gates that have an existing definition in sympy 

46_FIXED_GATE_MAP: dict[OpType, type[symgate.Gate]] = { 

47 OpType.H: symgate.HadamardGate, 

48 OpType.S: symgate.PhaseGate, 

49 OpType.CX: symgate.CNotGate, 

50 OpType.SWAP: symgate.SwapGate, 

51 OpType.T: symgate.TGate, 

52 OpType.X: symgate.XGate, 

53 OpType.Y: symgate.YGate, 

54 OpType.Z: symgate.ZGate, 

55} 

56 

57ParamsType = list[Expr | float] 

58# Make sure the return matrix is Immutable https://github.com/sympy/sympy/issues/18733 

59SymGateFunc = Callable[[ParamsType], ImmutableMatrix] 

60SymGateMap = dict[OpType, SymGateFunc] 

61 

62# Begin matrix definitions for symbolic OpTypes 

63# matches internal TKET definitions 

64# see OpType documentation 

65 

66 

67def symb_controlled(target: SymGateFunc) -> SymGateFunc: 

68 return lambda x: ImmutableMatrix(BlockDiagMatrix(Identity(2), target(x))) 

69 

70 

71def symb_rz(params: ParamsType) -> ImmutableMatrix: 

72 return ImmutableMatrix( 

73 [ 

74 [sympy.exp(-I * (sympy.pi / 2) * params[0]), 0], 

75 [0, sympy.exp(I * (sympy.pi / 2) * params[0])], 

76 ] 

77 ) 

78 

79 

80def symb_rx(params: ParamsType) -> ImmutableMatrix: 

81 costerm = sympy.cos((sympy.pi / 2) * params[0]) 

82 sinterm = -I * sympy.sin((sympy.pi / 2) * params[0]) 

83 return ImmutableMatrix( 

84 [ 

85 [costerm, sinterm], 

86 [sinterm, costerm], 

87 ] 

88 ) 

89 

90 

91def symb_ry(params: ParamsType) -> ImmutableMatrix: 

92 costerm = sympy.cos((sympy.pi / 2) * params[0]) 

93 sinterm = sympy.sin((sympy.pi / 2) * params[0]) 

94 return ImmutableMatrix( 

95 [ 

96 [costerm, -sinterm], 

97 [sinterm, costerm], 

98 ] 

99 ) 

100 

101 

102def symb_u3(params: ParamsType) -> ImmutableMatrix: 

103 theta, phi, lam = params 

104 costerm = sympy.cos((sympy.pi / 2) * theta) 

105 sinterm = sympy.sin((sympy.pi / 2) * theta) 

106 return ImmutableMatrix( 

107 [ 

108 [costerm, -sinterm * sympy.exp(I * sympy.pi * lam)], 

109 [ 

110 sinterm * sympy.exp(I * sympy.pi * phi), 

111 costerm * sympy.exp(I * sympy.pi * (phi + lam)), 

112 ], 

113 ] 

114 ) 

115 

116 

117def symb_u2(params: ParamsType) -> ImmutableMatrix: 

118 return symb_u3([0.5] + params) 

119 

120 

121def symb_u1(params: ParamsType) -> ImmutableMatrix: 

122 return symb_u3([0.0, 0.0] + params) 

123 

124 

125def symb_tk1(params: ParamsType) -> ImmutableMatrix: 

126 return symb_rz([params[0]]) * symb_rx([params[1]]) * symb_rz([params[2]]) 

127 

128 

129def symb_tk2(params: ParamsType) -> ImmutableMatrix: 

130 return ( 

131 symb_xxphase([params[0]]) 

132 * symb_yyphase([params[1]]) 

133 * symb_zzphase([params[2]]) 

134 ) 

135 

136 

137def symb_iswap(params: ParamsType) -> ImmutableMatrix: 

138 alpha = params[0] 

139 costerm = sympy.cos((sympy.pi / 2) * alpha) 

140 sinterm = sympy.sin((sympy.pi / 2) * alpha) 

141 return ImmutableMatrix( 

142 [ 

143 [1, 0, 0, 0], 

144 [0, costerm, I * sinterm, 0], 

145 [0, I * sinterm, costerm, 0], 

146 [0, 0, 0, 1], 

147 ] 

148 ) 

149 

150 

151def symb_phasediswap(params: ParamsType) -> ImmutableMatrix: 

152 p, alpha = params 

153 costerm = sympy.cos((sympy.pi / 2) * alpha) 

154 sinterm = I * sympy.sin((sympy.pi / 2) * alpha) 

155 phase = sympy.exp(2 * I * sympy.pi * p) 

156 return ImmutableMatrix( 

157 [ 

158 [1, 0, 0, 0], 

159 [0, costerm, sinterm * phase, 0], 

160 [0, sinterm / phase, costerm, 0], 

161 [0, 0, 0, 1], 

162 ] 

163 ) 

164 

165 

166def symb_xxphase(params: ParamsType) -> ImmutableMatrix: 

167 alpha = params[0] 

168 c = sympy.cos((sympy.pi / 2) * alpha) 

169 s = -I * sympy.sin((sympy.pi / 2) * alpha) 

170 return ImmutableMatrix( 

171 [ 

172 [c, 0, 0, s], 

173 [0, c, s, 0], 

174 [0, s, c, 0], 

175 [s, 0, 0, c], 

176 ] 

177 ) 

178 

179 

180def symb_yyphase(params: ParamsType) -> ImmutableMatrix: 

181 alpha = params[0] 

182 c = sympy.cos((sympy.pi / 2) * alpha) 

183 s = I * sympy.sin((sympy.pi / 2) * alpha) 

184 return ImmutableMatrix( 

185 [ 

186 [c, 0, 0, s], 

187 [0, c, -s, 0], 

188 [0, -s, c, 0], 

189 [s, 0, 0, c], 

190 ] 

191 ) 

192 

193 

194def symb_zzphase(params: ParamsType) -> ImmutableMatrix: 

195 alpha = params[0] 

196 t = sympy.exp(I * (sympy.pi / 2) * alpha) 

197 return ImmutableMatrix(diag(1 / t, t, t, 1 / t)) 

198 

199 

200def symb_xxphase3(params: ParamsType) -> ImmutableMatrix: 

201 xxphase2 = symb_xxphase(params) 

202 res1 = matrix_tensor_product(xxphase2, eye(2)) 

203 res2 = Matrix( 

204 BlockMatrix( 

205 [ 

206 [xxphase2[:2, :2], zeros(2), xxphase2[:2, 2:], zeros(2)], 

207 [zeros(2), xxphase2[:2, :2], zeros(2), xxphase2[:2, 2:]], 

208 [xxphase2[2:, :2], zeros(2), xxphase2[2:, 2:], zeros(2)], 

209 [zeros(2), xxphase2[2:, :2], zeros(2), xxphase2[2:, 2:]], 

210 ] 

211 ) 

212 ) 

213 res3 = matrix_tensor_product(eye(2), xxphase2) 

214 return ImmutableMatrix(res1 * res2 * res3) 

215 

216 

217def symb_phasedx(params: ParamsType) -> ImmutableMatrix: 

218 alpha, beta = params 

219 

220 return symb_rz([beta]) * symb_rx([alpha]) * symb_rz([-beta]) 

221 

222 

223def symb_eswap(params: ParamsType) -> ImmutableMatrix: 

224 alpha = params[0] 

225 c = sympy.cos((sympy.pi / 2) * alpha) 

226 s = -I * sympy.sin((sympy.pi / 2) * alpha) 

227 t = sympy.exp(-I * (sympy.pi / 2) * alpha) 

228 

229 return ImmutableMatrix( 

230 [ 

231 [t, 0, 0, 0], 

232 [0, c, s, 0], 

233 [0, s, c, 0], 

234 [0, 0, 0, t], 

235 ] 

236 ) 

237 

238 

239def symb_fsim(params: ParamsType) -> ImmutableMatrix: 

240 alpha, beta = params 

241 c = sympy.cos(sympy.pi * alpha) 

242 s = -I * sympy.sin(sympy.pi * alpha) 

243 t = sympy.exp(-I * sympy.pi * beta) 

244 

245 return ImmutableMatrix( 

246 [ 

247 [1, 0, 0, 0], 

248 [0, c, s, 0], 

249 [0, s, c, 0], 

250 [0, 0, 0, t], 

251 ] 

252 ) 

253 

254 

255def symb_gpi(params: ParamsType) -> ImmutableMatrix: 

256 t = sympy.exp(I * sympy.pi * params[0]) 

257 

258 return ImmutableMatrix( 

259 [ 

260 [0, 1 / t], 

261 [t, 0], 

262 ] 

263 ) 

264 

265 

266def symb_gpi2(params: ParamsType) -> ImmutableMatrix: 

267 t = sympy.exp(I * sympy.pi * params[0]) 

268 c = 1 / sympy.sqrt(2) 

269 

270 return c * ImmutableMatrix( 

271 [ 

272 [1, -I / t], 

273 [-I * t, 1], 

274 ] 

275 ) 

276 

277 

278def symb_aams(params: ParamsType) -> ImmutableMatrix: 

279 alpha, beta, gamma = params 

280 c = sympy.cos(sympy.pi / 2 * alpha) 

281 s = sympy.sin(sympy.pi / 2 * alpha) 

282 s1 = -I * sympy.exp(I * sympy.pi * (-beta - gamma)) * s 

283 s2 = -I * sympy.exp(I * sympy.pi * (-beta + gamma)) * s 

284 s3 = -I * sympy.exp(I * sympy.pi * (beta - gamma)) * s 

285 s4 = -I * sympy.exp(I * sympy.pi * (beta + gamma)) * s 

286 

287 return ImmutableMatrix( 

288 [ 

289 [c, 0, 0, s1], 

290 [0, c, s2, 0], 

291 [0, s3, c, 0], 

292 [s4, 0, 0, c], 

293 ] 

294 ) 

295 

296 

297# end symbolic matrix definitions 

298 

299 

300class SymGateRegister: 

301 """Static class holding mapping from OpType to callable generating symbolic matrix. 

302 Allows users to add their own definitions, or override existing definitions.""" 

303 

304 _g_map: SymGateMap = { 

305 OpType.Rx: symb_rx, 

306 OpType.Ry: symb_ry, 

307 OpType.Rz: symb_rz, 

308 OpType.TK1: symb_tk1, 

309 OpType.TK2: symb_tk2, 

310 OpType.U1: symb_u1, 

311 OpType.U2: symb_u2, 

312 OpType.U3: symb_u3, 

313 OpType.CRx: symb_controlled(symb_rx), 

314 OpType.CRy: symb_controlled(symb_ry), 

315 OpType.CRz: symb_controlled(symb_rz), 

316 OpType.CU1: symb_controlled(symb_u1), 

317 OpType.CU3: symb_controlled(symb_u3), 

318 OpType.ISWAP: symb_iswap, 

319 OpType.PhasedISWAP: symb_phasediswap, 

320 OpType.XXPhase: symb_xxphase, 

321 OpType.YYPhase: symb_yyphase, 

322 OpType.ZZPhase: symb_zzphase, 

323 OpType.XXPhase3: symb_xxphase3, 

324 OpType.PhasedX: symb_phasedx, 

325 OpType.ESWAP: symb_eswap, 

326 OpType.FSim: symb_fsim, 

327 OpType.GPI: symb_gpi, 

328 OpType.GPI2: symb_gpi2, 

329 OpType.AAMS: symb_aams, 

330 } 

331 

332 @classmethod 

333 def register_func(cls, typ: OpType, f: SymGateFunc, replace: bool = False) -> None: 

334 """Register a callable for an optype. 

335 

336 :param typ: OpType to register 

337 :type typ: OpType 

338 :param f: Callable for generating symbolic matrix. 

339 :type f: SymGateFunc 

340 :param replace: Whether to replace existing entry, defaults to False 

341 :type replace: bool, optional 

342 """ 

343 if typ not in cls._g_map or replace: 

344 cls._g_map[typ] = f 

345 

346 @classmethod 

347 def get_func(cls, typ: OpType) -> SymGateFunc: 

348 """Get registered callable.""" 

349 return cls._g_map[typ] 

350 

351 @classmethod 

352 def is_registered(cls, typ: OpType) -> bool: 

353 """Check if type has a callable registered.""" 

354 return typ in cls._g_map 

355 

356 

357def _op_to_sympy_gate(op: Op, targets: list[int]) -> symgate.Gate: 

358 # convert Op to sympy gate 

359 if op.type in _FIXED_GATE_MAP: 

360 return _FIXED_GATE_MAP[op.type](*targets) 

361 if op.is_gate(): 361 ↛ 365line 361 didn't jump to line 365 because the condition on line 361 was always true

362 # check if symbolic definition is needed 

363 float_params = all(isinstance(p, float) for p in op.params) 

364 else: 

365 raise ValueError( 

366 f"Circuit can only contain unitary gates, operation {op} not valid." 

367 ) 

368 

369 # pytket matrix basis indexing is in opposite order to sympy 

370 targets.reverse() 

371 if (not float_params) and SymGateRegister.is_registered(op.type): 

372 u_mat = SymGateRegister.get_func(op.type)(op.params) 

373 else: 

374 try: 

375 # use internal tket unitary definition 

376 u_mat = ImmutableMatrix(op.get_unitary()) 

377 except RuntimeError as e: 

378 # to catch tket failure to get Op unitary 

379 # most likely due to symbolic parameters. 

380 raise ValueError( 

381 f"{op.type} is not supported for symbolic conversion." 

382 " Try registering your own symbolic matrix representation" 

383 " with SymGateRegister.func." 

384 ) from e 

385 return symgate.UGate(targets, u_mat) 

386 

387 

388def circuit_to_symbolic_gates(circ: Circuit) -> Mul: 

389 """Generate a multiplication expression of sympy gates from Circuit 

390 

391 :param circ: Input circuit 

392 :type circ: Circuit 

393 :raises ValueError: If circ does not match a unitary operation. 

394 :return: Symbolic gate multiplication expression. 

395 :rtype: Mul 

396 """ 

397 outmat = symgate.IdentityGate(0) 

398 nqb = circ.n_qubits 

399 qubit_map = {qb: nqb - 1 - i for i, qb in enumerate(circ.qubits)} 

400 for com in circ: 

401 op = com.op 

402 if op.type == OpType.Barrier: 402 ↛ 403line 402 didn't jump to line 403 because the condition on line 402 was never true

403 continue 

404 args = com.args 

405 try: 

406 targs = [qubit_map[q] for q in args] # type: ignore 

407 except KeyError as e: 

408 raise ValueError( 

409 f"Gates can only act on qubits. Operation {com} not valid." 

410 ) from e 

411 gate = _op_to_sympy_gate(op, targs) 

412 

413 outmat = gate * outmat 

414 

415 for i in range(len(qubit_map)): 

416 outmat = symgate.IdentityGate(i) * outmat 

417 

418 return outmat * sympy.exp(circ.phase * sympy.pi * I) 

419 

420 

421def circuit_to_symbolic_unitary(circ: Circuit) -> ImmutableMatrix: 

422 """Generate a symbolic unitary from Circuit. 

423 

424 Unitary matches pytket default ILO BasisOrder. 

425 

426 :param circ: Input circuit 

427 :type circ: Circuit 

428 :return: Symbolic unitary. 

429 :rtype: ImmutableMatrix 

430 """ 

431 gates = circuit_to_symbolic_gates(circ) 

432 nqb = circ.n_qubits 

433 try: 

434 return cast(ImmutableMatrix, represent(gates, nqubits=circ.n_qubits)) 

435 except NotImplementedError: 

436 # sympy can't represent n>1 qubit unitaries very well 

437 # so if it fails we will just calculate columns using the statevectors 

438 # for all possible input basis states 

439 matrix_dim = 1 << nqb 

440 input_states = (Qubit(f"{i:0{nqb}b}") for i in range(matrix_dim)) 

441 outmat = Matrix([]) 

442 for col, input_state in enumerate(input_states): 

443 outmat = outmat.col_insert(col, represent(qapply(gates * input_state))) 

444 

445 return ImmutableMatrix(outmat) 

446 

447 

448def circuit_apply_symbolic_qubit(circ: Circuit, input_qb: Expr) -> Qubit: 

449 """Apply circuit to an input state to calculate output symbolic state. 

450 

451 :param circ: Input Circuit. 

452 :type circ: Circuit 

453 :param input_qb: Sympy Qubit expression corresponding to a state. 

454 :type input_qb: Expr 

455 :return: Output state after circuit acts on input_qb. 

456 :rtype: Qubit 

457 """ 

458 gates = circuit_to_symbolic_gates(circ) 

459 

460 return cast(Qubit, qapply(gates * input_qb)) 

461 

462 

463def circuit_apply_symbolic_statevector( 

464 circ: Circuit, input_state: np.ndarray | ImmutableMatrix | None = None 

465) -> ImmutableMatrix: 

466 """Apply circuit to an optional input statevector 

467 to calculate output symbolic statevector. 

468 If no input statevector given, the all zero state is assumed. 

469 Statevector follows pytket default ILO BasisOrder. 

470 

471 :param circ: Input Circuit. 

472 :type circ: Circuit 

473 :param input_state: Input statevector as a column vector, defaults to None. 

474 :type input_state: Optional[Union[np.ndarray, ImmutableMatrix]], optional 

475 :return: Symbolic state after circ acts on input_state. 

476 :rtype: ImmutableMatrix 

477 """ 

478 if input_state: 478 ↛ 479line 478 didn't jump to line 479 because the condition on line 478 was never true

479 input_qb = matrix_to_qubit(input_state) 

480 else: 

481 input_qb = Qubit("0" * circ.n_qubits) 

482 return cast( 

483 ImmutableMatrix, 

484 represent(circuit_apply_symbolic_qubit(circ, cast(Qubit, input_qb))), 

485 )