1# This code is part of Qiskit. 2# 3# (C) Copyright IBM 2018, 2019. 4# 5# This code is licensed under the Apache License, Version 2.0. You may 6# obtain a copy of this license in the LICENSE.txt file in the root directory 7# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. 8# 9# Any modifications or derivative works of this code must retain this 10# copyright notice, and modified files need to carry a notice indicating 11# that they have been altered from the originals. 12""" 13Helper functions for noise model creation. 14""" 15 16import numpy as np 17 18from qiskit.quantum_info.operators.operator import Operator 19from qiskit.quantum_info.operators.channel.kraus import Kraus 20from qiskit.quantum_info.operators.channel.superop import SuperOp 21from qiskit.quantum_info.operators.predicates import is_identity_matrix 22from qiskit.quantum_info.operators.predicates import is_unitary_matrix 23from qiskit.quantum_info.operators.predicates import matrix_equal 24from qiskit.quantum_info.operators.predicates import ATOL_DEFAULT 25 26from ..noiseerror import NoiseError 27 28 29def standard_gates_instructions(instructions): 30 """Convert a list with unitary matrix instructions into standard gates. 31 32 Args: 33 instructions (list): A list of qobj instructions. 34 35 Returns: 36 list: a list of qobj instructions equivalent to in input instruction. 37 """ 38 output_instructions = [] 39 for instruction in instructions: 40 output_instructions += standard_gate_instruction(instruction) 41 return output_instructions 42 43 44# pylint: disable=too-many-return-statements 45def standard_gate_instruction(instruction, ignore_phase=True): 46 """Convert a unitary matrix instruction into a standard gate instruction. 47 48 Args: 49 instruction (dict): A qobj instruction. 50 ignore_phase (bool): Ignore global phase on unitary matrix in 51 comparison to canonical unitary. 52 53 Returns: 54 list: a list of qobj instructions equivalent to in input instruction. 55 """ 56 57 name = instruction.get("name", None) 58 if name not in ["mat", "unitary", "kraus"]: 59 return [instruction] 60 qubits = instruction["qubits"] 61 params = instruction["params"] 62 63 # Check for single-qubit reset Kraus 64 if name == "kraus": 65 if len(qubits) == 1: 66 superop = SuperOp(Kraus(params)) 67 # Check if reset to |0> 68 reset0 = reset_superop(1) 69 if superop == reset0: 70 return [{"name": "reset", "qubits": qubits}] 71 # Check if reset to |1> 72 reset1 = reset0.compose(Operator(standard_gate_unitary('x'))) 73 if superop == reset1: 74 return [{"name": "reset", "qubits": qubits}, {"name": "x", "qubits": qubits}] 75 # otherwise just return the kraus instruction 76 return [instruction] 77 78 # Check single qubit gates 79 mat = params[0] 80 if len(qubits) == 1: 81 # Check clifford gates 82 for j in range(24): 83 if matrix_equal( 84 mat, 85 single_qubit_clifford_matrix(j), 86 ignore_phase=ignore_phase): 87 return single_qubit_clifford_instructions(j, qubit=qubits[0]) 88 # Check t gates 89 for name in ["t", "tdg"]: 90 if matrix_equal( 91 mat, 92 standard_gate_unitary(name), 93 ignore_phase=ignore_phase): 94 return [{"name": name, "qubits": qubits}] 95 # TODO: u1,u2,u3 decomposition 96 # Check two qubit gates 97 if len(qubits) == 2: 98 for name in ["cx", "cz", "swap"]: 99 if matrix_equal( 100 mat, 101 standard_gate_unitary(name), 102 ignore_phase=ignore_phase): 103 return [{"name": name, "qubits": qubits}] 104 # Check reversed CX 105 if matrix_equal( 106 mat, 107 standard_gate_unitary("cx_10"), 108 ignore_phase=ignore_phase): 109 return [{"name": "cx", "qubits": [qubits[1], qubits[0]]}] 110 # Check 2-qubit Pauli's 111 paulis = ["id", "x", "y", "z"] 112 for pauli0 in paulis: 113 for pauli1 in paulis: 114 pmat = np.kron( 115 standard_gate_unitary(pauli1), 116 standard_gate_unitary(pauli0)) 117 if matrix_equal(mat, pmat, ignore_phase=ignore_phase): 118 if pauli0 == "id": 119 return [{"name": pauli1, "qubits": [qubits[1]]}] 120 elif pauli1 == "id": 121 return [{"name": pauli0, "qubits": [qubits[0]]}] 122 else: 123 return [{ 124 "name": pauli0, 125 "qubits": [qubits[0]] 126 }, { 127 "name": pauli1, 128 "qubits": [qubits[1]] 129 }] 130 # Check three qubit toffoli 131 if len(qubits) == 3: 132 if matrix_equal( 133 mat, 134 standard_gate_unitary("ccx_012"), 135 ignore_phase=ignore_phase): 136 return [{"name": "ccx", "qubits": qubits}] 137 if matrix_equal( 138 mat, 139 standard_gate_unitary("ccx_021"), 140 ignore_phase=ignore_phase): 141 return [{ 142 "name": "ccx", 143 "qubits": [qubits[0], qubits[2], qubits[1]] 144 }] 145 if matrix_equal( 146 mat, 147 standard_gate_unitary("ccx_120"), 148 ignore_phase=ignore_phase): 149 return [{ 150 "name": "ccx", 151 "qubits": [qubits[1], qubits[2], qubits[0]] 152 }] 153 154 # Else return input in 155 return [instruction] 156 157 158def single_qubit_clifford_gates(j): 159 """Return a QASM gate names for a single qubit Clifford. 160 161 The labels are returned in a basis set consisting of 162 ('id', 's', 'sdg', 'z', 'h', x', 'y') gates decomposed to 163 use the minimum number of X-90 pulses in a (u1, u2, u3) 164 decomposition. 165 166 Args: 167 j (int): Clifford index 0, ..., 23. 168 169 Returns: 170 tuple(str): The tuple of basis gates. 171 172 Raises: 173 NoiseError: If index is out of range [0, 23]. 174 """ 175 176 if not isinstance(j, int) or j < 0 or j > 23: 177 raise NoiseError( 178 "Index {} must be in the range [0, ..., 23]".format(j)) 179 180 labels = [ 181 ('id', ), 182 ('s', ), 183 ('sdg', ), 184 ('z', ), 185 # u2 gates 186 ( 187 'h', ), 188 ('h', 'z'), 189 ('z', 'h'), 190 ('h', 's'), 191 ('s', 'h'), 192 ('h', 'sdg'), 193 ('sdg', 'h'), 194 ('s', 'h', 's'), 195 ('sdg', 'h', 's'), 196 ('z', 'h', 's'), 197 ('s', 'h', 'sdg'), 198 ('sdg', 'h', 'sdg'), 199 ('z', 'h', 'sdg'), 200 ('s', 'h', 'z'), 201 ('sdg', 'h', 'z'), 202 ('z', 'h', 'z'), 203 # u3 gates 204 ( 205 'x', ), 206 ('y', ), 207 ('s', 'x'), 208 ('sdg', 'x') 209 ] 210 return labels[j] 211 212 213def single_qubit_clifford_matrix(j): 214 """Return Numpy array for a single qubit Clifford. 215 216 Args: 217 j (int): Clifford index 0, ..., 23. 218 219 Returns: 220 np.array: The matrix for the indexed clifford. 221 222 Raises: 223 NoiseError: If index is out of range [0, 23]. 224 """ 225 226 if not isinstance(j, int) or j < 0 or j > 23: 227 raise NoiseError( 228 "Index {} must be in the range [0, ..., 23]".format(j)) 229 230 basis_dict = { 231 'id': np.eye(2), 232 'x': np.array([[0, 1], [1, 0]], dtype=complex), 233 'y': np.array([[0, -1j], [1j, 0]], dtype=complex), 234 'z': np.array([[1, 0], [0, -1]], dtype=complex), 235 'h': np.array([[1, 1], [1, -1]], dtype=complex) / np.sqrt(2), 236 's': np.array([[1, 0], [0, 1j]], dtype=complex), 237 'sdg': np.array([[1, 0], [0, -1j]], dtype=complex) 238 } 239 mat = np.eye(2) 240 for gate in single_qubit_clifford_gates(j): 241 mat = np.dot(basis_dict[gate], mat) 242 return mat 243 244 245# pylint: disable=invalid-name 246def single_qubit_clifford_instructions(index, qubit=0): 247 """Return a list of qobj instructions for a single qubit Cliffords. 248 249 The instructions are returned in a basis set consisting of 250 ('id', 's', 'sdg', 'z', 'h', x', 'y') gates decomposed to 251 use the minimum number of X-90 pulses in a (u1, u2, u3) 252 decomposition. 253 254 Args: 255 index (int): Clifford index 0, ..., 23. 256 qubit (int): the qubit to apply the Clifford to. 257 258 Returns: 259 list(dict): The list of instructions. 260 261 Raises: 262 NoiseError: If index is out of range [0, 23] or qubit invalid. 263 """ 264 265 if not isinstance(index, int) or index < 0 or index > 23: 266 raise NoiseError( 267 "Index {} must be in the range [0, ..., 23]".format(index)) 268 if not isinstance(qubit, int) or qubit < 0: 269 raise NoiseError("qubit position must be positive integer.") 270 271 instructions = [] 272 for gate in single_qubit_clifford_gates(index): 273 instructions.append({"name": gate, "qubits": [qubit]}) 274 return instructions 275 276 277def standard_gate_unitary(name): 278 """Return the unitary matrix for a standard gate.""" 279 280 unitary_matrices = { 281 ("id", "I"): 282 np.eye(2, dtype=complex), 283 ("x", "X"): 284 np.array([[0, 1], [1, 0]], dtype=complex), 285 ("y", "Y"): 286 np.array([[0, -1j], [1j, 0]], dtype=complex), 287 ("z", "Z"): 288 np.array([[1, 0], [0, -1]], dtype=complex), 289 ("h", "H"): 290 np.array([[1, 1], [1, -1]], dtype=complex) / np.sqrt(2), 291 ("s", "S"): 292 np.array([[1, 0], [0, 1j]], dtype=complex), 293 ("sdg", "Sdg"): 294 np.array([[1, 0], [0, -1j]], dtype=complex), 295 ("t", "T"): 296 np.array([[1, 0], [0, np.exp(1j * np.pi / 4)]], dtype=complex), 297 ("tdg", "Tdg"): 298 np.array([[1, 0], [0, np.exp(-1j * np.pi / 4)]], dtype=complex), 299 ("cx", "CX", "cx_01"): 300 np.array([[1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0]], dtype=complex), 301 ("cx_10",): 302 np.array([[1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0]], dtype=complex), 303 ("cz", "CZ"): 304 np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, -1]], dtype=complex), 305 ("swap", "SWAP"): 306 np.array([[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]], dtype=complex), 307 ("ccx", "CCX", "ccx_012", "ccx_102"): 308 np.array([[1, 0, 0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0, 0], 309 [0, 0, 1, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 1], 310 [0, 0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 0, 0, 1, 0, 0], 311 [0, 0, 0, 0, 0, 0, 1, 0], [0, 0, 0, 1, 0, 0, 0, 0]], 312 dtype=complex), 313 ("ccx_021", "ccx_201"): 314 np.array([[1, 0, 0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0, 0], 315 [0, 0, 1, 0, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0, 0], 316 [0, 0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 1], 317 [0, 0, 0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 0, 1, 0, 0]], 318 dtype=complex), 319 ("ccx_120", "ccx_210"): 320 np.array([[1, 0, 0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0, 0], 321 [0, 0, 1, 0, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0, 0], 322 [0, 0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 0, 0, 1, 0, 0], 323 [0, 0, 0, 0, 0, 0, 0, 1], [0, 0, 0, 0, 0, 0, 1, 0]], 324 dtype=complex) 325 } 326 327 return next((value for key, value in unitary_matrices.items() if name in key), None) 328 329 330def reset_superop(num_qubits): 331 """Return a N-qubit reset SuperOp.""" 332 reset = SuperOp( 333 np.array([[1, 0, 0, 1], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]])) 334 if num_qubits == 1: 335 return reset 336 reset_n = reset 337 for _ in range(num_qubits - 1): 338 reset_n.tensor(reset) 339 return reset_n 340 341 342def standard_instruction_operator(instr): 343 """Return the Operator for a standard gate instruction.""" 344 # Convert to dict (for QobjInstruction types) 345 if hasattr(instr, 'as_dict'): 346 instr = instr.as_dict() 347 # Get name and parameters 348 name = instr.get('name', "") 349 params = instr.get('params', []) 350 # Check if standard unitary gate name 351 mat = standard_gate_unitary(name) 352 if isinstance(mat, np.ndarray): 353 return Operator(mat) 354 355 # Check if standard parameterized waltz gates 356 if name == 'u1': 357 lam = params[0] 358 mat = np.diag([1, np.exp(1j * lam)]) 359 return Operator(mat) 360 if name == 'u2': 361 phi = params[0] 362 lam = params[1] 363 mat = np.array([[1, -np.exp(1j * lam)], 364 [np.exp(1j * phi), 365 np.exp(1j * (phi + lam))]]) / np.sqrt(2) 366 return Operator(mat) 367 if name == 'u3': 368 theta = params[0] 369 phi = params[1] 370 lam = params[2] 371 mat = np.array( 372 [[np.cos(theta / 2), -np.exp(1j * lam) * np.sin(theta / 2)], 373 [np.exp(1j * phi) * np.sin(theta / 2), np.exp(1j * (phi + lam)) * np.cos(theta / 2)]]) 374 return Operator(mat) 375 376 # Check if unitary instruction 377 if name == 'unitary': 378 return Operator(params[0]) 379 380 # Otherwise return None if we cannot convert instruction 381 return None 382 383 384def standard_instruction_channel(instr): 385 """Return the SuperOp channel for a standard instruction.""" 386 # Check if standard operator 387 oper = standard_instruction_operator(instr) 388 if oper is not None: 389 return SuperOp(oper) 390 391 # Convert to dict (for QobjInstruction types) 392 if hasattr(instr, 'as_dict'): 393 instr = instr.as_dict() 394 # Get name and parameters 395 name = instr.get('name', "") 396 397 # Check if reset instruction 398 if name == 'reset': 399 # params should be the number of qubits being reset 400 num_qubits = len(instr['qubits']) 401 return reset_superop(num_qubits) 402 # Check if Kraus instruction 403 if name == 'kraus': 404 params = instr['params'] 405 return SuperOp(Kraus(params)) 406 return None 407 408 409def circuit2superop(circuit, min_qubits=1): 410 """Return the SuperOp for a standard instruction.""" 411 # Get number of qubits 412 max_qubits = 1 413 for instr in circuit: 414 qubits = [] 415 if hasattr(instr, 'qubits'): 416 qubits = instr.qubits 417 elif isinstance(instr, dict): 418 qubits = instr.get('qubits', []) 419 max_qubits = max(max_qubits, 1 + max(qubits)) 420 421 num_qubits = max(max_qubits, min_qubits) 422 423 # Initialize N-qubit identity superoperator 424 superop = SuperOp(np.eye(4**num_qubits)) 425 # compose each circuit element with the superoperator 426 for instr in circuit: 427 instr_op = standard_instruction_channel(instr) 428 if instr_op is None: 429 raise NoiseError('Cannot convert instruction {} to SuperOp'.format(instr)) 430 if hasattr(instr, 'qubits'): 431 qubits = instr.qubits 432 else: 433 qubits = instr['qubits'] 434 superop = superop.compose(instr_op, qubits) 435 return superop 436 437 438def make_unitary_instruction(mat, qubits, standard_gates=True): 439 """Return a qobj instruction for a unitary matrix gate. 440 441 Args: 442 mat (matrix): A square or diagonal unitary matrix. 443 qubits (list[int]): The qubits the matrix is applied to. 444 standard_gates (bool): Check if the matrix instruction is a 445 standard instruction. 446 447 Returns: 448 dict: The qobj instruction object. 449 450 Raises: 451 NoiseError: if the input is not a unitary matrix. 452 """ 453 if not is_unitary_matrix(mat): 454 raise NoiseError("Input matrix is not unitary.") 455 456 if isinstance(qubits, int): 457 qubits = [qubits] 458 459 instruction = {"name": "unitary", "qubits": qubits, "params": [mat]} 460 if standard_gates: 461 return standard_gate_instruction(instruction) 462 else: 463 return [instruction] 464 465 466def make_kraus_instruction(mats, qubits): 467 """Return a qobj instruction for a Kraus error. 468 469 Args: 470 mats (list[matrix]): A list of square or diagonal Kraus matrices. 471 qubits (list[int]): The qubits the matrix is applied to. 472 Returns: 473 dict: The qobj instruction object. 474 475 Raises: 476 NoiseError: if the input is not a CPTP Kraus channel. 477 """ 478 kraus = Kraus(mats) 479 if not kraus.is_cptp() or kraus._input_dim != kraus._output_dim: 480 raise NoiseError("Input Kraus matrices are not a CPTP channel.") 481 if isinstance(qubits, int): 482 qubits = [qubits] 483 return [{"name": "kraus", "qubits": qubits, "params": kraus.data}] 484 485 486def qubits_from_mat(mat): 487 """Return the number of qubits for a multi-qubit matrix.""" 488 arr = np.array(mat) 489 shape = arr.shape 490 num_qubits = int(np.log2(shape[1])) 491 if shape[1] != 2**num_qubits: 492 raise NoiseError("Input Kraus channel is not a multi-qubit channel.") 493 return num_qubits 494 495 496def is_matrix_diagonal(mat): 497 """Test if row-vector representation of diagonal matrix.""" 498 mat = np.array(mat) 499 shape = mat.shape 500 return len(shape) == 2 and shape[0] == 1 501 502 503def kraus2instructions(kraus_ops, standard_gates, atol=ATOL_DEFAULT): 504 """ 505 Convert a list of Kraus matrices into qobj circuits. 506 507 If any Kraus operators are a unitary matrix they will be converted 508 into unitary qobj instructions. Identity unitary matrices will also be 509 converted into identity qobj instructions. 510 511 Args: 512 kraus_ops (list[matrix]): A list of Kraus matrices for a CPTP map. 513 standard_gates (bool): Check if the matrix instruction is a 514 standard instruction (default: True). 515 atol (double): Threshold for testing if probabilities are zero. 516 517 518 Returns: 519 list: A list of pairs (p, circuit) where `circuit` is a list of qobj 520 instructions, and `p` is the probability of that circuit for the 521 given error. 522 523 Raises: 524 NoiseError: If the input Kraus channel is not CPTP. 525 """ 526 # Check threshold 527 if atol < 0: 528 raise NoiseError("atol cannot be negative") 529 if atol > 1e-5: 530 raise NoiseError( 531 "atol value is too large. It should be close to zero.") 532 533 # Check CPTP 534 if not Kraus(kraus_ops).is_cptp(atol=atol): 535 raise NoiseError("Input Kraus channel is not CPTP.") 536 537 # Get number of qubits 538 num_qubits = int(np.log2(len(kraus_ops[0]))) 539 if len(kraus_ops[0]) != 2**num_qubits: 540 raise NoiseError("Input Kraus channel is not a multi-qubit channel.") 541 542 # Check if each matrix is a: 543 # 1. scaled identity matrix 544 # 2. scaled non-identity unitary matrix 545 # 3. a non-unitary Kraus operator 546 547 # Probabilities 548 prob_identity = 0 549 prob_unitary = 0 # total probability of all unitary ops (including id) 550 prob_kraus = 0 # total probability of non-unitary ops 551 probabilities = [] # initialize with probability of Identity 552 553 # Matrices 554 unitaries = [] # non-identity unitaries 555 non_unitaries = [] # non-unitary Kraus matrices 556 557 for mat in kraus_ops: 558 # Get the value of the maximum diagonal element 559 # of op.H * op for rescaling 560 # pylint: disable=no-member 561 prob = abs(max(np.diag(np.conj(np.transpose(mat)).dot(mat)))) 562 if prob > 0.0: 563 if abs(prob - 1) > 0.0: 564 # Rescale the operator by square root of prob 565 rescaled_mat = np.array(mat) / np.sqrt(prob) 566 else: 567 rescaled_mat = mat 568 # Check if identity operator 569 if is_identity_matrix(rescaled_mat, ignore_phase=True): 570 prob_identity += prob 571 prob_unitary += prob 572 # Check if unitary 573 elif is_unitary_matrix(rescaled_mat): 574 probabilities.append(prob) 575 prob_unitary += prob 576 unitaries.append(rescaled_mat) 577 # Non-unitary op 578 else: 579 non_unitaries.append(mat) 580 581 # Check probabilities 582 prob_kraus = 1 - prob_unitary 583 if prob_unitary - 1 > atol: 584 raise NoiseError("Invalid kraus matrices: unitary probability " 585 "{} > 1".format(prob_unitary)) 586 if prob_unitary < -atol: 587 raise NoiseError("Invalid kraus matrices: unitary probability " 588 "{} < 1".format(prob_unitary)) 589 if prob_identity - 1 > atol: 590 raise NoiseError("Invalid kraus matrices: identity probability " 591 "{} > 1".format(prob_identity)) 592 if prob_identity < -atol: 593 raise NoiseError("Invalid kraus matrices: identity probability " 594 "{} < 1".format(prob_identity)) 595 if prob_kraus - 1 > atol: 596 raise NoiseError("Invalid kraus matrices: non-unitary probability " 597 "{} > 1".format(prob_kraus)) 598 if prob_kraus < -atol: 599 raise NoiseError("Invalid kraus matrices: non-unitary probability " 600 "{} < 1".format(prob_kraus)) 601 602 # Build qobj instructions 603 instructions = [] 604 qubits = list(range(num_qubits)) 605 606 # Add unitary instructions 607 for unitary in unitaries: 608 instructions.append( 609 make_unitary_instruction( 610 unitary, qubits, standard_gates=standard_gates)) 611 612 # Add identity instruction 613 if prob_identity > atol: 614 if abs(prob_identity - 1) < atol: 615 probabilities.append(1) 616 else: 617 probabilities.append(prob_identity) 618 instructions.append([{"name": "id", "qubits": [0]}]) 619 620 # Add Kraus 621 if prob_kraus < atol: 622 # No Kraus operators 623 return zip(instructions, probabilities) 624 if prob_kraus < 1: 625 # Rescale kraus operators by probabilities 626 non_unitaries = [ 627 np.array(op) / np.sqrt(prob_kraus) for op in non_unitaries 628 ] 629 instructions.append(make_kraus_instruction(non_unitaries, qubits)) 630 probabilities.append(prob_kraus) 631 # Normalize probabilities to account for any rounding errors 632 probabilities = list(np.array(probabilities) / np.sum(probabilities)) 633 return zip(instructions, probabilities) 634