Anders and Briegel in Python
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

147 lines
4.7KB

  1. #!/usr/bin/python
  2. # -*- coding: utf-8 -*-
  3. """
  4. Exposes a few basic QI operators
  5. And a circuit-model simulator
  6. """
  7. import numpy as np
  8. import itertools as it
  9. from fractions import Fraction
  10. def hermitian_conjugate(u):
  11. """ Shortcut to the Hermitian conjugate """
  12. return np.conjugate(np.transpose(u))
  13. # Constants
  14. ir2 = 1 / np.sqrt(2)
  15. # Operators
  16. id = np.array(np.eye(2, dtype=complex))
  17. px = np.array([[0, 1], [1, 0]], dtype=complex)
  18. py = np.array([[0, -1j], [1j, 0]], dtype=complex)
  19. pz = np.array([[1, 0], [0, -1]], dtype=complex)
  20. ha = hadamard = np.array([[1, 1], [1, -1]], dtype=complex) * ir2
  21. ph = np.array([[1, 0], [0, 1j]], dtype=complex)
  22. t = np.array([[1, 0], [0, np.exp(1j * np.pi / 4)]], dtype=complex)
  23. sqx = np.array(
  24. [[1. + 0.j, -0. + 1.j], [-0. + 1.j, 1. - 0.j]], dtype=complex) * ir2
  25. msqx = np.array(
  26. [[1. + 0.j, 0. - 1.j], [0. - 1.j, 1. - 0.j]], dtype=complex) * ir2
  27. sqy = np.array(
  28. [[1. + 0.j, 1. + 0.j], [-1. - 0.j, 1. - 0.j]], dtype=complex) * ir2
  29. msqy = np.array(
  30. [[1. + 0.j, -1. - 0.j], [1. + 0.j, 1. - 0.j]], dtype=complex) * ir2
  31. sqz = np.array(
  32. [[1. + 1.j, 0. + 0.j], [0. + 0.j, 1. - 1.j]], dtype=complex) * ir2
  33. msqz = np.array(
  34. [[1. - 1.j, 0. + 0.j], [0. + 0.j, 1. + 1.j]], dtype=complex) * ir2
  35. # CZ gate
  36. cz = np.array(np.eye(4), dtype=complex)
  37. cz[3, 3] = -1
  38. # States
  39. zero = np.array([[1], [0]], dtype=complex)
  40. one = np.array([[0], [1]], dtype=complex)
  41. plus = np.array([[1], [1]], dtype=complex) / np.sqrt(2)
  42. bond = cz.dot(np.kron(plus, plus))
  43. nobond = np.kron(plus, plus)
  44. # Labelling stuff
  45. common_us = id, px, py, pz, ha, ph, sqz, msqz, sqy, msqy, sqx, msqx
  46. names = "identity", "px", "py", "pz", "hadamard", "phase", "sqz", "msqz", "sqy", "msqy", "sqx", "msqx"
  47. by_name = dict(zip(names, common_us))
  48. paulis = px, py, pz
  49. operators = id, px, py, pz
  50. def normalize_global_phase(m):
  51. """ Normalize the global phase of a matrix """
  52. v = (x for x in m.flatten() if np.abs(x) > 0.001).next()
  53. phase = np.arctan2(v.imag, v.real) % np.pi
  54. rot = np.exp(-1j * phase)
  55. return rot * m if rot * v > 0 else -rot * m
  56. class CircuitModel(object):
  57. def __init__(self, nqubits):
  58. self.nqubits = nqubits
  59. self.d = 2 ** nqubits
  60. self.state = np.zeros((self.d, 1), dtype=complex)
  61. self.state[0, 0] = 1
  62. def act_cz(self, control, target):
  63. """ Act a CU somewhere. """
  64. control = 1 << control
  65. target = 1 << target
  66. for i in xrange(self.d):
  67. if (i & control) and (i & target):
  68. self.state[i, 0] *= -1
  69. def act_cnot(self, control, target):
  70. """ Act a CNOT. """
  71. self.act_hadamard(target)
  72. self.act_cz(control, target)
  73. self.act_hadamard(target)
  74. def act_hadamard(self, qubit):
  75. """ Act a hadamard somewhere. """
  76. where = 1 << qubit
  77. output = np.zeros((self.d, 1), dtype=complex)
  78. for i, v in enumerate(self.state):
  79. q = int(i & where > 0)
  80. output[i] += v * ha[q, q]
  81. output[i ^ where] += v * ha[int(not q), q]
  82. self.state = output
  83. def act_local_rotation(self, qubit, u):
  84. """ Act a local unitary somwhere. """
  85. where = 1 << qubit
  86. output = np.zeros((self.d, 1), dtype=complex)
  87. for i, v in enumerate(self.state):
  88. q = int(i & where > 0)
  89. output[i] += v * u[q, q]
  90. output[i ^ where] += v * u[int(not q), q]
  91. self.state = output
  92. def act_circuit(self, circuit):
  93. """ Act a sequence of gates. """
  94. for node, operation in circuit:
  95. if operation == "cz":
  96. self.act_cz(*node)
  97. else:
  98. self.act_local_rotation(node, operation)
  99. def __eq__(self, other):
  100. """ Check whether two quantum states are the same or not,
  101. up to a global phase. """
  102. a = normalize_global_phase(self.state)
  103. b = normalize_global_phase(other.state)
  104. return np.allclose(a, b)
  105. def __setitem__(self, key, value):
  106. """ Set a matrix element """
  107. self.state[key] = value
  108. def __getitem__(self, key):
  109. """ Get a matrix element """
  110. return self.state[key]
  111. def __str__(self):
  112. s = ""
  113. for i in range(self.d):
  114. label = bin(i)[2:].rjust(self.nqubits, "0")[::-1]
  115. if abs(self.state[i, 0]) > 0.00001:
  116. term = self.state[i, 0]
  117. real_sign = " " if term.real>=0 else "-"
  118. real_frac = Fraction(str(term.real**2)).limit_denominator()
  119. imag_sign = "+" if term.imag>=0 else "-"
  120. imag_frac = Fraction(str(term.imag**2)).limit_denominator()
  121. s += "|{}❭: \t{}√{}\t{} i √{}\n".format(
  122. label, real_sign, real_frac, imag_sign, imag_frac)
  123. return s