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

1# Copyright Quantinuum

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

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

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.

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."""

20from collections.abc import Callable

21from typing import cast

23import numpy as np

24import sympy

25from sympy import (

26 BlockDiagMatrix,

27 BlockMatrix,

28 Expr,

29 I,

30 Identity,

31 ImmutableMatrix,

32 Matrix,

33 Mul,

34 diag,

35 eye,

36 zeros,

37)

38from sympy.physics.quantum import gate as symgate

39from sympy.physics.quantum import represent

40from sympy.physics.quantum.qapply import qapply

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

42from sympy.physics.quantum.tensorproduct import matrix_tensor_product

44from pytket.circuit import Circuit, Op, OpType

46# gates that have an existing definition in sympy

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

48 OpType.H: symgate.HadamardGate,

49 OpType.S: symgate.PhaseGate,

50 OpType.CX: symgate.CNotGate,

51 OpType.SWAP: symgate.SwapGate,

52 OpType.T: symgate.TGate,

53 OpType.X: symgate.XGate,

54 OpType.Y: symgate.YGate,

55 OpType.Z: symgate.ZGate,

56}

58ParamsType = list[Expr | float]

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

60SymGateFunc = Callable[[ParamsType], ImmutableMatrix]

61SymGateMap = dict[OpType, SymGateFunc]

63# Begin matrix definitions for symbolic OpTypes

64# matches internal TKET definitions

65# see OpType documentation

68def _symb_controlled(target: SymGateFunc) -> SymGateFunc:

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

72def _symb_rz(params: ParamsType) -> ImmutableMatrix:

73 return ImmutableMatrix(

74 [

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

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

77 ]

78 )

81def _symb_rx(params: ParamsType) -> ImmutableMatrix:

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

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

84 return ImmutableMatrix(

85 [

86 [costerm, sinterm],

87 [sinterm, costerm],

88 ]

89 )

92def _symb_ry(params: ParamsType) -> ImmutableMatrix:

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

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

95 return ImmutableMatrix(

96 [

97 [costerm, -sinterm],

98 [sinterm, costerm],

99 ]

100 )

101

102

103def _symb_u3(params: ParamsType) -> ImmutableMatrix:

104 theta, phi, lam = params

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

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

107 return ImmutableMatrix(

108 [

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

110 [

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

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

113 ],

114 ]

115 )

116

117

118def _symb_u2(params: ParamsType) -> ImmutableMatrix:

119 return _symb_u3([0.5, *params])

120

121

122def _symb_u1(params: ParamsType) -> ImmutableMatrix:

123 return _symb_u3([0.0, 0.0, *params])

124

125

126def _symb_tk1(params: ParamsType) -> ImmutableMatrix:

127 return _symb_rz([params[0]]) * _symb_rx([params[1]]) * _symb_rz([params[2]])

128

129

130def _symb_tk2(params: ParamsType) -> ImmutableMatrix:

131 return (

132 _symb_xxphase([params[0]])

133 * _symb_yyphase([params[1]])

134 * _symb_zzphase([params[2]])

135 )

136

137

138def _symb_iswap(params: ParamsType) -> ImmutableMatrix:

139 alpha = params[0]

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

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

142 return ImmutableMatrix(

143 [

144 [1, 0, 0, 0],

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

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

147 [0, 0, 0, 1],

148 ]

149 )

150

151

152def _symb_phasediswap(params: ParamsType) -> ImmutableMatrix:

153 p, alpha = params

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

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

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

157 return ImmutableMatrix(

158 [

159 [1, 0, 0, 0],

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

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

162 [0, 0, 0, 1],

163 ]

164 )

165

166

167def _symb_xxphase(params: ParamsType) -> ImmutableMatrix:

168 alpha = params[0]

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

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

171 return ImmutableMatrix(

172 [

173 [c, 0, 0, s],

174 [0, c, s, 0],

175 [0, s, c, 0],

176 [s, 0, 0, c],

177 ]

178 )

179

180

181def _symb_yyphase(params: ParamsType) -> ImmutableMatrix:

182 alpha = params[0]

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

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

185 return ImmutableMatrix(

186 [

187 [c, 0, 0, s],

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

189 [0, -s, c, 0],

190 [s, 0, 0, c],

191 ]

192 )

193

194

195def _symb_zzphase(params: ParamsType) -> ImmutableMatrix:

196 alpha = params[0]

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

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

199

200

201def _symb_xxphase3(params: ParamsType) -> ImmutableMatrix:

202 xxphase2 = _symb_xxphase(params)

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

204 res2 = Matrix(

205 BlockMatrix(

206 [

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

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

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

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

211 ]

212 )

213 )

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

215 return ImmutableMatrix(res1 * res2 * res3)

216

217

218def _symb_phasedx(params: ParamsType) -> ImmutableMatrix:

219 alpha, beta = params

220

221 return _symb_rz([beta]) * _symb_rx([alpha]) * _symb_rz([-beta])

222

223

224def _symb_eswap(params: ParamsType) -> ImmutableMatrix:

225 alpha = params[0]

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

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

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

229

230 return ImmutableMatrix(

231 [

232 [t, 0, 0, 0],

233 [0, c, s, 0],

234 [0, s, c, 0],

235 [0, 0, 0, t],

236 ]

237 )

238

239

240def _symb_fsim(params: ParamsType) -> ImmutableMatrix:

241 alpha, beta = params

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

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

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

245

246 return ImmutableMatrix(

247 [

248 [1, 0, 0, 0],

249 [0, c, s, 0],

250 [0, s, c, 0],

251 [0, 0, 0, t],

252 ]

253 )

254

255

256def _symb_gpi(params: ParamsType) -> ImmutableMatrix:

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

258

259 return ImmutableMatrix(

260 [

261 [0, 1 / t],

262 [t, 0],

263 ]

264 )

265

266

267def _symb_gpi2(params: ParamsType) -> ImmutableMatrix:

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

269 c = 1 / sympy.sqrt(2)

270

271 return c * ImmutableMatrix(

272 [

273 [1, -I / t],

274 [-I * t, 1],

275 ]

276 )

277

278

279def _symb_aams(params: ParamsType) -> ImmutableMatrix:

280 alpha, beta, gamma = params

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

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

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

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

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

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

287

288 return ImmutableMatrix(

289 [

290 [c, 0, 0, s1],

291 [0, c, s2, 0],

292 [0, s3, c, 0],

293 [s4, 0, 0, c],

294 ]

295 )

296

297

298# end symbolic matrix definitions

299

300

301class SymGateRegister:

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

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

304

305 _g_map: SymGateMap = { # noqa: RUF012

306 OpType.Rx: _symb_rx,

307 OpType.Ry: _symb_ry,

308 OpType.Rz: _symb_rz,

309 OpType.TK1: _symb_tk1,

310 OpType.TK2: _symb_tk2,

311 OpType.U1: _symb_u1,

312 OpType.U2: _symb_u2,

313 OpType.U3: _symb_u3,

314 OpType.CRx: _symb_controlled(_symb_rx),

315 OpType.CRy: _symb_controlled(_symb_ry),

316 OpType.CRz: _symb_controlled(_symb_rz),

317 OpType.CU1: _symb_controlled(_symb_u1),

318 OpType.CU3: _symb_controlled(_symb_u3),

319 OpType.ISWAP: _symb_iswap,

320 OpType.PhasedISWAP: _symb_phasediswap,

321 OpType.XXPhase: _symb_xxphase,

322 OpType.YYPhase: _symb_yyphase,

323 OpType.ZZPhase: _symb_zzphase,

324 OpType.XXPhase3: _symb_xxphase3,

325 OpType.PhasedX: _symb_phasedx,

326 OpType.ESWAP: _symb_eswap,

327 OpType.FSim: _symb_fsim,

328 OpType.GPI: _symb_gpi,

329 OpType.GPI2: _symb_gpi2,

330 OpType.AAMS: _symb_aams,

331 }

332

333 @classmethod

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

335 """Register a callable for an optype.

336

337 :param typ: OpType to register

338 :param f: Callable for generating symbolic matrix.

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

340 """

341 if typ not in cls._g_map or replace:

342 cls._g_map[typ] = f

343

344 @classmethod

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

346 """Get registered callable."""

347 return cls._g_map[typ]

348

349 @classmethod

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

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

352 return typ in cls._g_map

353

354

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

356 # convert Op to sympy gate

357 if op.type in _FIXED_GATE_MAP:

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

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

360 # check if symbolic definition is needed

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

362 else:

363 raise ValueError(

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

365 )

366

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

368 targets.reverse()

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

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

371 else:

372 try:

373 # use internal tket unitary definition

374 u_mat = ImmutableMatrix(op.get_unitary())

375 except RuntimeError as e:

376 # to catch tket failure to get Op unitary

377 # most likely due to symbolic parameters.

378 raise ValueError(

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

380 " Try registering your own symbolic matrix representation"

381 " with SymGateRegister.func."

382 ) from e

383 return symgate.UGate(targets, u_mat)

384

385

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

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

388

389 :param circ: Input circuit

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

391 :return: Symbolic gate multiplication expression.

392 """

393 outmat = symgate.IdentityGate(0)

394 nqb = circ.n_qubits

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

396 for com in circ:

397 op = com.op

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

399 continue

400 args = com.args

401 try:

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

403 except KeyError as e:

404 raise ValueError(

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

406 ) from e

407 gate = _op_to_sympy_gate(op, targs)

408

409 outmat = gate * outmat

410

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

412 outmat = symgate.IdentityGate(i) * outmat

413

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

415

416

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

418 """Generate a symbolic unitary from Circuit.

419

420 Unitary matches pytket default ILO BasisOrder.

421

422 :param circ: Input circuit

423 :return: Symbolic unitary.

424 """

425 gates = circuit_to_symbolic_gates(circ)

426 nqb = circ.n_qubits

427 try:

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

429 except NotImplementedError:

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

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

432 # for all possible input basis states

433 matrix_dim = 1 << nqb

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

435 outmat = Matrix([])

436 for col, input_state in enumerate(input_states):

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

438

439 return ImmutableMatrix(outmat)

440

441

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

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

444

445 :param circ: Input Circuit.

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

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

448 """

449 gates = circuit_to_symbolic_gates(circ)

450

451 return cast("Qubit", qapply(gates * input_qb))

452

453

454def circuit_apply_symbolic_statevector(

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

456) -> ImmutableMatrix:

457 """Apply circuit to an optional input statevector

458 to calculate output symbolic statevector.

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

460 Statevector follows pytket default ILO BasisOrder.

461

462 :param circ: Input Circuit.

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

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

465 """

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

467 input_qb = matrix_to_qubit(input_state)

468 else:

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

470 return cast(

471 "ImmutableMatrix",

472 represent(circuit_apply_symbolic_qubit(circ, cast("Qubit", input_qb))),

473 )