diff --git a/docs/projectq.ops.rst b/docs/projectq.ops.rst index 8fbe2e287..6b9ec5936 100755 --- a/docs/projectq.ops.rst +++ b/docs/projectq.ops.rst @@ -53,7 +53,9 @@ The operations collection consists of various default gates and is a work-in-pro projectq.ops.UniformlyControlledRy projectq.ops.UniformlyControlledRz projectq.ops.StatePreparation + projectq.ops.QPE projectq.ops.FlipBits + projectq.ops.QAA Module contents diff --git a/docs/projectq.setups.decompositions.rst b/docs/projectq.setups.decompositions.rst index 883395db6..26cd392fc 100755 --- a/docs/projectq.setups.decompositions.rst +++ b/docs/projectq.setups.decompositions.rst @@ -26,6 +26,8 @@ The decomposition package is a collection of gate decomposition / replacement ru projectq.setups.decompositions.time_evolution projectq.setups.decompositions.toffoli2cnotandtgate projectq.setups.decompositions.uniformlycontrolledr2cnot + projectq.setups.decompositions.phaseestimation + projectq.setups.decompositions.amplitudeamplification Submodules @@ -172,6 +174,20 @@ projectq.setups.decompositions.uniformlycontrolledr2cnot module :members: :undoc-members: +projectq.setups.decompositions.phaseestimation module +--------------------------------------------------------------- + +.. automodule:: projectq.setups.decompositions.phaseestimation + :members: + :undoc-members: + +projectq.setups.decompositions.amplitudeamplification module +--------------------------------------------------------------- + +.. automodule:: projectq.setups.decompositions.amplitudeamplification + :members: + :undoc-members: + Module contents --------------- diff --git a/examples/bernstein_vazirani_tutorial.ipynb b/examples/bernstein_vazirani_tutorial.ipynb new file mode 100644 index 000000000..8d93c0339 --- /dev/null +++ b/examples/bernstein_vazirani_tutorial.ipynb @@ -0,0 +1,302 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Bernstein-Varzirani Algorithm\n", + "The Bernstein-Varzirani algorithm solves the following problem:\n", + "1. Given an n-bit secret number x\n", + "2. There is an oracle holding the secret key that one can query with a number y, and the oracle will return $x \\bullet y$, the bit-wise dot product modulo 2 of the two numbers. For example, given a 4-bit secret number $x=13$ with binary representation 1101 and $y=15$ with binary representation 1111, then $x \\bullet y =$ 3 mod 2, which is 1\n", + "3. The question is: how many oracle queries does one need to figure out x?\n", + "\n", + "With the classical approach, if we query the oracle with 1, 2(10), 4(100) and 8(1000), then we can figure out each bit of the secret key and thus we need four queries. For a general n-bit number, we need n queries. Surprisingly, in the quantum computing paradigm, the Bernstein-Varzirani algorithm needs only one measurement, regardless of the size of the secret number" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import projectq\n", + "from projectq.backends import CircuitDrawer\n", + "from projectq.ops import H, Z, All, Measure, Barrier\n", + "\n", + "from pathlib import Path" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "n = 3 #No. of qubits" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## The classical approach" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "def query_oracle(y):\n", + " \"\"\"\n", + " The oracle returns the bit-wise dot product modulo 2 \n", + " of the secret key and y\n", + " \"\"\"\n", + " x = 6 #secret key\n", + " s = bin(x & y).count('1') % 2\n", + " \n", + " return s" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Secret key is 6 after 3 tries\n" + ] + } + ], + "source": [ + "b = ''.join([str(query_oracle(2**(k))) for k in reversed(range(n))])\n", + "print(\"Secret key is {} after {} tries\".format(int(b, 2), n))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## The quantum computing approach\n", + "There are five steps in the Bernstein-Varzirani circuit construction:\n", + "1. Initialize all qubits to 0 (This is the default value anyway)\n", + "2. Apply Hadamard transform to all qubits\n", + "3. Implement the oracle \n", + "4. Apply Hadamard transform to all qubits\n", + "5. Measure the circuit\n", + "\n", + "The workings of the Bernstein-Varzirani algorithm can only be understood by working through the mathematics. After step 2, the application of the Hadamard transform to all qubits results in\n", + "$$\\begin{align} H^{\\otimes n} |0\\rangle = \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} | y \\rangle \\end{align}$$\n", + "where $N = 2^n$\n", + "\n", + "The oracle seeks to implement the following operation in step 3:\n", + "$$\\begin{align} \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} | y \\rangle \\rightarrow \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} (-1)^{x \\bullet y} | y \\rangle \\end{align}$$\n", + "The oracle effectively multiplies each superposition state $|y \\rangle$ with -1 raised to the power of $x \\bullet y$, the bit-wise inner product modulo 2. But why does the oracle implement this operation?\n", + "\n", + "Using the identity (proof given below)\n", + "$$\\begin{align} H^{\\otimes n} |x\\rangle = \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} (-1)^{x \\bullet y} | y \\rangle \\end{align}$$\n", + "If we further apply another Hadamard transform to all qubits in step 4, then $$ H^{\\otimes n} H^{\\otimes n} |x\\rangle = H^{\\otimes n} \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} (-1)^{x \\bullet y} | y \\rangle $$\n", + "i.e. $$ |x\\rangle = H^{\\otimes n} \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} (-1)^{x \\bullet y} | y \\rangle $$\n", + "since the Hadamard transform is symmetric and unitary (meaning that the transform is its own inverse). Thus measuring the circuit in step 5 will give the secret key.\n", + "\n", + "Summarizing,\n", + "$$ H^{\\otimes n} |0\\rangle ^{\\otimes n} \\rightarrow \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} | y \\rangle \\rightarrow \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} (-1)^{x \\bullet y} | y \\rangle \\rightarrow H^{\\otimes n} \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} (-1)^{x \\bullet y} | y \\rangle = |x\\rangle$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### The oracle implementation\n", + "First note that in implementing\n", + "$$ \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} (-1)^{x \\bullet y} | y \\rangle = \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} (-1)^{(x_{n-1} y_{n-1} \\oplus \\ldots \\oplus x_0 y_0)} | y_{n-1}\\ldots y_0\\rangle $$\n", + "$$ = \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} (-1)^{x_{n-1} y_{n-1}} (-1)^{x_{n-2} y_{n-2}} \\ldots (-1)^{x_0 y_0}) | y_{n-1}\\ldots y_0\\rangle $$\n", + "the individual product ${x_{n-1} y_{n-1}, \\ldots, x_0 y_0}$ only needs to be evaluated for bits in x that are 1. \n", + "\n", + "For the bits in x that are 1, we only introduce a minus sign for the bits in y that are 1 i.e.\n", + "$ |0 \\rangle \\rightarrow |0 \\rangle$ and $ |1 \\rangle \\rightarrow -|1 \\rangle$. This is precisely what the Z gate does. Thus the oracle can be implemented as: if a bit in x is 1, apply the Z gate to the corresponding bit in y." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "def make_bv_circuit(engine, n):\n", + " circuit = engine.allocate_qureg(n)\n", + " All(H) | circuit\n", + " Barrier | circuit\n", + " \n", + " # Oracle\n", + " x = 6 # secret key\n", + " for i in range(n):\n", + " if x & 1: # Only apply Z if the current bit of the secret key x is 1\n", + " Z | circuit[i]\n", + " x >>= 1 # Move the next bit to the 1 position\n", + " Barrier | circuit\n", + " \n", + " All(H) | circuit\n", + " All(Measure) | circuit\n", + " engine.flush()\n", + "\n", + " return circuit" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Draw circuit diagram" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "drawing_engine = CircuitDrawer()\n", + "main_engine = projectq.MainEngine(drawing_engine)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "circuit = make_bv_circuit(main_engine, n)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAR0AAACBCAQAAABXsKZ8AAAAAmJLR0QA/4ePzL8AAAAJcEhZcwAAAEgAAABIAEbJaz4AAAAHdElNRQfjChIULw91KGppAAAKh0lEQVR42u2dX6gcVx3HP3vN9ZKKKXNTYyAS60QF9aVhEi1iHwJzQSSIPuxFKLUUYbcIYlOqOyhUFJFdERRE5a4++OLL3Uf/vNxRqgVBegdpDa1B3BjaUDXJHdqYFk2T9WFmZ2dnZ2Z3zszendn7+zwkOztz5vzOme+cOX/uzrc2oIYiK3yMD2ZO9Tp/YY+9DLk8xAOsZs7lIle4qlq0AljjQaW4X+QlbsTuW+WTrGeO42Uuc52bwBE+o36xJzmUK/X7Bz/PmqT2JP/mtUxJTg5+kDmXr+erlkJQibuVsnOFdYX6vsBt/0ZdZX3w/aIKV3typahTzcxqkcoXppKvcUjmyP5LR1gOVjzpmIkHmLOeSVhmmjHfTZNO2rNXOCDYGDHfRp+Epi+jHo6fqk5vtgz69P1PGkZkSyie/avvHlsx3473dUxaWFhYbPn5d6nPnoWLxYavOXBphraE4tmf+k6SwHir0wpaGJsWm4AL6IGgU9HRcXD8Z5y31ZEn3tzYn/p2cWjE7gm3OhpmIBI7kFovS7sjLBvJYgxLJ/qI1AHoySjr4NJnKINJ0iaMhg8qB0O6LMtPH5s+OjouDhomBt2UB+Asc41dWrED+wQ2QsEI86eI+u7ioGMGvZo6fWx+yKWc0umjEzf3U+MjcYfvBJ8sOpO7T3CGd2dYxUrIZQonOMtJrqtUY0GsKsZ9hmO8nnDG++O+nlLfH+AO69wCjnA0Lr1Nj/pEV1ingcMvsDASuslh6biRfU5wlj5gT6Rd4bhC5Vxll0tcm/l41VyeW/jK+X1Kce8mrpyvcUzhjH9jlxd4DTjKyehOFwsjdtbGmwjU2cKmSSNmrigsHQcXzf9s0A+k1Jh1UlCoFi4W7eCSRxlOBJqYWEyOosanBLvBfj0kFyOmxREqT7pwxicC23QnHkrj0ulgYOKtXFn+d3URznLSSRGOixPp3LaxIuJZiaTYwKWBHRpRmbM+rmwsukCT7sSWUDz56tuiniicuIlAjXbQmnhMjrCcsTkcDW3WMZ+JSTthSyiePPVtpy6SeqtPUTQaWKE8pv2pV0MeV8uHi5261pW0+GBAqF2ZJh1TnjfLR/oiaR8SH2XtkBw86SS3LCKcJcRI6eUQrMXH0wj6L15fJ1k6MqOzhNRz7B31kOTP2gVFRDqCIvl+pXOx9kTm3zcucl2p6rxcu5D5iv019Plm7anCfgV3M590rvFnBplTHeJuQeEfLAb0ua2Q7u3+NbrDZaX0CeSRzl32uKmU8s3iCnCAuM2NiYWkWRjwBgD/Ybc8vzl/o8CKEaYxULxRh7yV8HdAikg3WVBEpCMoItIRFBHpCIqIdARFRDqCIiIdQRGRjqCISEdQRKQjKCLSERQR6QiKiHQERUQ6giIiHUERkY6giEhHUESkIygi0hEUEekIiuT7s/ZDvE0p3f8UfoIzr/JXqQSrSrf6gNsMgFqxXmR5pFPj7D4YOM6TapUgr4HjO8tj4FjbJwPH+VGtEoiBY+UNHKtVAjFwFEqGGDgKsyAGjoIS5TZwjB7kKP0wvwD6kRe2Tol+gXGLgeMQjRY7/kvudHbYSn3hXSId1nM54nTZoOefwWUj8q7gecWtFrUYOHrl7mGEXu/dY5ftkF/PvuGwHVTVJlqC9UYZ4hYDx6SIe5iJBnAptNhTSRagB8W16dHIei7FuPNGPT+qaOCo4y6ivzOsBJcmhsr75hcU93yIf2+7R1kNHA3qNP1L8I6O/6GFhk0Xh21semxhRLY1LPqjC+590Oii0cBgE5M6zbTSDBvnDn1288Stcd4KoiYSpxsbtX+PPs2bftQOZijuL823yvs4OGgYuDjoGBipD8CyGTia/r8mZ4ILvOnSBjpssoOJzilsWnRwMCLbDbZGyzwdXvGNnLY5RZ+/Y9KhNcuN4PgH5oj729wzipqJOBOj/gr/CEVth+I+za3C69un64tleNPVcXD4Kc8XYuA42XbX+HDc4coGjvcHn2yu8Aw/49dAjcNem2P4ox0d77EyXH8c39aGuXyWh/m8f0kP8y3+AAww2WQ0j5VoDtlET660aAni4z7Ce8JRR+OMjRpO82m+GIl6GPcqH1Wo70/xANf4L3A4/m2z8QaOBgZNflyIgePkKHWFhzlPVpINHE1G7423ucMWf6THCsfbD3YZv6fMmKSRXD4EvOpv/ZJXsIHvTpQg1hyyg8O2fzm7kxUXLUF83M/yWPcn0ZbATNkC4Djwr0jUw7jVDBx/w58CA8dHojunGTie9g0c45yzymvg2GfoxXXuDFuY2HH3VDp25P8ZM+5QD8ZZmR8Cw7i/wDd0tah/73t/7oMnUIdWYjd4ZOBoxPr1ldfAUWPYwX+sHbo/Z76UvyM8OMgw8m1C6D7MLB0vbp0fsa0Q9QvAe1WiVsGinphFeCJQi/HgK6+Bo0YLlw7wvmFPwAt95mAu8WzwrGnMPsPbw/bHRF4FZhxnD+PWwevUZoz6Krs8mj1qFSzqiQOBqIFjnHgORVJsYNDADo1BzFlbXBvbtxD0ThHemopOgzqjzriJwwYu8Ft+1TnvXQKTJvewCWwGdu5eLsNtGxvowime4XG+TBsbgz4a28C2f3gKTTRcLLw1ImfapFZS3Dbf47wdRG1MxBkbNcDTfC6I2vEPiI07V32nCifZwHGsVzQAUma+vCqPZ4VHBpnhq3yCd81QtlAue4OdwV62XL7GxzmBN0OupZ17biVY4/GsUfvnPpoQ9RqPKkT7FA9xLwBHeWL47c6gnZJmb9CK/X439D3frISBY/r1T8XFXtjc7q2yRu3Sm2LgGN9uiYHjgaeTusDSJ7l/LgaOBxwttTVMb5HqQbsjBo4HkFaOvaNhg/xZu6CISEdQ5KAbOFarBEtk4Hi98gaOVSrBQNGAcWTg2Oet4sLJZ+B4Q9HXrSwGjtUqwW2uK/3SfWTg6JTlN+fLYOBYpRKIgaOwHIh0BEVEOoIiIh1BEZGOoIhIR1BEpCMoItIRFBHpCIqIdARFRDqCIiIdQRGRjqCISEdQRKQjKCLSERQR6QiKiHQERUQ6giIiHUERMXCsUgnEwJEbJXkpdbVKIAaOYuCoiBg4Vsr+sPolEANHoWSIgaMwC3GuYGLgKEzFjn1BXFkMHEtDRgPHUkR60A0cS0JmA8cFIgaOpSKjgeMCEQPHkpHTwHHpqKKB44LIbeC4VFTRwHFh5DRwrDAudsjA0cBAT3xvO5TPwLEkZDZwXCBF1HcPGw0zZOBo0+Myz6e0u2UycIwn0WQxlROc5STXVc+dycAxnlXFuM9wLOHFbashq8gQRRg4mhPDAROTJt9Jcbopk4FjPAkmi1NzeY4rU19TW4yBYzxr3KcU9y4v+S58k2cs3MARrMRRpI3BOc5V0MBxYeQ0cKwUaa5Y1TVwXBg5DRwrRJpwqmvguDByGjhWiHQfvuoYOJaGjAaOC2SeBo7WxDkmDRwnR1jO2ByOhjZrm22GBnfRreqg8jr0xZCnvu3URVIXYvZq1Mfmlith4CgUi6qBo4krBo4Hm/RFUjFwFBLRUhd1xcBRSEQMHIWFItIRFKmCgePF2gXu5XCmNP8sQQku1i6gsZbhzHd5NXV/qQwc/w9+31fX9CH7ZQAAACV0RVh0ZGF0ZTpjcmVhdGUAMjAxOS0xMC0xOFQyMDo0NzoxNSswODowMK2N2cUAAAAldEVYdGRhdGU6bW9kaWZ5ADIwMTktMTAtMThUMjA6NDc6MTUrMDg6MDDc0GF5AAAAFHRFWHRwZGY6VmVyc2lvbgBQREYtMS41IAVcCzkAAABKdEVYdHNpZ25hdHVyZQA2YTczNWVhMDZhYjBhZWQxNGFkNzdlYzI5MDM3MDg0OWE1ZGM5MDFkYmIwZjc0ODE2MThmYzIyYzJlOTE5OGM1OseklwAAAABJRU5ErkJggg==\n", + "text/plain": [ + "" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "p = Path('diagram')\n", + "if not p.exists(): #if the diagram directory doesn't exist, create it\n", + " p.mkdir()\n", + "with open('diagram/bv.tex', 'w') as f:\n", + " latex = drawing_engine.get_latex() #get circuit diagram as latex\n", + " f.write(latex) \n", + "#Change the pdf scale to 1.8 from 0.8 to have better visual effect\n", + "!sed -i 's@tikzpicture\\}\\[scale=0.8@tikzpicture\\}\\[scale=1.8@g' diagram/bv.tex\n", + "!cd diagram; pdflatex bv.tex > /dev/null #convert tex to latex, piping to /dev/null to silent output \n", + "\n", + "#Wand package needed to convert pdf to image\n", + "from wand.image import Image as WImage\n", + "img = WImage(filename='diagram/bv.pdf')\n", + "img" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Run circuit" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Measurement outcome: 110\n", + "Secret key: 6\n" + ] + } + ], + "source": [ + "main_engine = projectq.MainEngine()\n", + "circuit = make_bv_circuit(main_engine, n)\n", + "measurement = ''.join(str(int(c)) for c in reversed(circuit))\n", + "print('Measurement outcome:', measurement)\n", + "print('Secret key:', int(measurement, 2))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Proof of the key identity" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\\begin{align} H^{\\otimes n} |x\\rangle = \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} (-1)^{x \\bullet y} | y \\rangle \\end{align}$$ where $ N = 2^n $ and the $ \\bullet$ operator denotes bitwise dot product modulo 2 i.e. $ x \\bullet y = x_0y_0 \\oplus x_1y_1 \\oplus \\ldots x_{N-1}y_{N-1}$ where $\\oplus$ denotes the sum modulo 2 operation or the XOR operation. We prove this identity by induction.\n", + "\n", + "Base case, n = 1:\n", + "$$ H |x_0\\rangle = \\frac{1}{\\sqrt{2}} (|0\\rangle + (-1)^{x_0}|1\\rangle) = \\frac{1}{\\sqrt{2}} \\sum_{y=0}^{1} (-1)^{x_0y} | y \\rangle$$\n", + "\n", + "Assuming n-1 (n > 1) is true i.e. $$ H^{\\otimes (n-1)} |x_{n-2} \\ldots x_0\\rangle = \\frac{1}{\\sqrt{2^{n-1}}} \\sum_{y=0}^{2^{n-1} - 1} (-1)^{(x_{n-2} \\ldots x_0) \\bullet (y_{n-2} \\ldots y_0)} | y \\rangle$$ \n", + "\n", + "$$ = \\frac{1}{\\sqrt{2^{n-1}}} \\sum_{y=0}^{2^{n-1} - 1} (-1)^{(x_{n-2} y_{n-2} \\oplus \\ldots \\oplus x_0 y_0)} | y \\rangle$$ \n", + "\n", + "where $y = y_{n-2} \\ldots y_0$ is the binary representation\n", + "\n", + "\n", + "Then $$ H^{\\otimes n} |x_{n-1} \\ldots x_0 \\rangle = \\frac{1}{\\sqrt{2}} (|0\\rangle + (-1)^{x_{n-1}}|1\\rangle) \\times \\frac{1}{\\sqrt{2^{n-1}}} \\sum_{y=0}^{2^{n-1} - 1} (-1)^{(x_{n-2} y_{n-2} \\oplus \\ldots \\oplus x_0 y_0)} | y \\rangle $$\n", + "\n", + "$$ \\begin{align} = \\frac{1}{\\sqrt{2^n}} \\sum_{y=0}^{2^{n-1} - 1} (-1)^{(x_{n-2} y_{n-2} \\oplus \\ldots \\oplus x_0 y_0)} |0\\rangle| y \\rangle + \\frac{1}{\\sqrt{2^n}} \\sum_{y=0}^{2^{n-1} - 1} (-1)^{x_{n-1}} (-1)^{(x_{n-2} y_{n-2} \\oplus \\ldots \\oplus x_0 y_0)} |1\\rangle| y \\rangle \\end{align} $$\n", + "\n", + "$$ \\begin{align} = \\frac{1}{\\sqrt{2^n}} \\sum_{y=0}^{2^{n-1} - 1} (-1)^{(x_{n-1}0 \\oplus x_{n-2}y_{n-2} \\oplus \\ldots \\oplus x_0y_0)} | y \\rangle + \\frac{1}{\\sqrt{2^n}} \\sum_{y=2^{n-1}}^{2^{n} - 1} (-1)^{(x_{n-1}1 \\oplus x_{n-2}y_{n-2}\\oplus \\ldots \\oplus x_0y_0)} | y \\rangle \\end{align}$$\n", + "\n", + "$$ = \\frac{1}{\\sqrt{2^n}} \\sum_{y=0}^{2^{n-1} - 1} (-1)^{(x_{n-1}x_{n-2} \\ldots x_0) \\bullet (0y_{n-2} \\ldots y_0)} | y \\rangle + \\frac{1}{\\sqrt{2^n}} \\sum_{y=2^{n-1}}^{2^{n} - 1} (-1)^{(x_{n-1} \\ldots x_0) \\bullet (1y_{{n}-2} \\ldots y_0)} | y \\rangle $$\n", + "\n", + "$$ = \\frac{1}{\\sqrt{2^{n}}} \\sum_{y=0}^{2^{n} - 1} (-1)^{(x_{n-1} \\ldots x_0) \\bullet (y_{n-1} \\ldots y_0)} | y \\rangle $$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "projectq", + "language": "python", + "name": "projectq" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.8" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/projectq/backends/_resource_test.py b/projectq/backends/_resource_test.py index 664b687ad..cf7122c01 100755 --- a/projectq/backends/_resource_test.py +++ b/projectq/backends/_resource_test.py @@ -20,7 +20,7 @@ from projectq.cengines import DummyEngine, MainEngine, NotYetMeasuredError from projectq.meta import LogicalQubitIDTag -from projectq.ops import All, Allocate, CNOT, Command, H, Measure, QFT, Rz, X +from projectq.ops import All, Allocate, CNOT, Command, H, Measure, QFT, Rz, Rzz, X from projectq.types import WeakQubitRef from projectq.backends import ResourceCounter @@ -74,6 +74,7 @@ def test_resource_counter(): CNOT | (qubit1, qubit3) Rz(0.1) | qubit1 Rz(0.3) | qubit1 + Rzz(0.5) | qubit1 All(Measure) | qubit1 + qubit3 @@ -81,7 +82,7 @@ def test_resource_counter(): int(qubit1) assert resource_counter.max_width == 2 - assert resource_counter.depth_of_dag == 5 + assert resource_counter.depth_of_dag == 6 str_repr = str(resource_counter) assert str_repr.count(" HGate : 1") == 1 diff --git a/projectq/backends/_sim/_cppkernels/simulator.hpp b/projectq/backends/_sim/_cppkernels/simulator.hpp index c4e611eba..d248ed038 100755 --- a/projectq/backends/_sim/_cppkernels/simulator.hpp +++ b/projectq/backends/_sim/_cppkernels/simulator.hpp @@ -38,7 +38,7 @@ class Simulator{ public: using calc_type = double; using complex_type = std::complex; - using StateVector = std::vector>; + using StateVector = std::vector>; using Map = std::map; using RndEngine = std::mt19937; using Term = std::vector>; @@ -55,11 +55,18 @@ class Simulator{ void allocate_qubit(unsigned id){ if (map_.count(id) == 0){ map_[id] = N_++; - auto newvec = StateVector(1UL << N_); - #pragma omp parallel for schedule(static) + StateVector newvec; // avoid large memory allocations + if( tmpBuff1_.capacity() >= (1UL << N_) ) + std::swap(newvec, tmpBuff1_); + newvec.resize(1UL << N_); +#pragma omp parallel for schedule(static) for (std::size_t i = 0; i < newvec.size(); ++i) newvec[i] = (i < vec_.size())?vec_[i]:0.; - vec_ = std::move(newvec); + std::swap(vec_, newvec); + // recycle large memory + std::swap(tmpBuff1_, newvec); + if( tmpBuff1_.capacity() < tmpBuff2_.capacity() ) + std::swap(tmpBuff1_, tmpBuff2_); } else throw(std::runtime_error( @@ -113,12 +120,18 @@ class Simulator{ } } else{ - StateVector newvec((1UL << (N_-1))); - #pragma omp parallel for schedule(static) + StateVector newvec; // avoid costly memory reallocations + if( tmpBuff1_.capacity() >= (1UL << (N_-1)) ) + std::swap(tmpBuff1_, newvec); + newvec.resize((1UL << (N_-1))); + #pragma omp parallel for schedule(static) if(0) for (std::size_t i = 0; i < vec_.size(); i += 2*delta) std::copy_n(&vec_[i + static_cast(value)*delta], delta, &newvec[i/2]); - vec_ = std::move(newvec); + std::swap(vec_, newvec); + std::swap(tmpBuff1_, newvec); + if( tmpBuff1_.capacity() < tmpBuff2_.capacity() ) + std::swap(tmpBuff1_, tmpBuff2_); for (auto& p : map_){ if (p.second > pos) @@ -189,8 +202,8 @@ class Simulator{ } template - void apply_controlled_gate(M const& m, std::vector ids, - std::vector ctrl){ + void apply_controlled_gate(M const& m, const std::vector& ids, + const std::vector& ctrl){ auto fused_gates = fused_gates_; fused_gates.insert(m, ids, ctrl); @@ -209,8 +222,8 @@ class Simulator{ } template - void emulate_math(F const& f, QuReg quregs, std::vector ctrl, - unsigned num_threads=1){ + void emulate_math(F const& f, QuReg quregs, const std::vector& ctrl, + bool parallelize = false){ run(); auto ctrlmask = get_control_mask(ctrl); @@ -218,37 +231,76 @@ class Simulator{ for (unsigned j = 0; j < quregs[i].size(); ++j) quregs[i][j] = map_[quregs[i][j]]; - StateVector newvec(vec_.size(), 0.); - std::vector res(quregs.size()); - - #pragma omp parallel for schedule(static) firstprivate(res) num_threads(num_threads) - for (std::size_t i = 0; i < vec_.size(); ++i){ - if ((ctrlmask&i) == ctrlmask){ - for (unsigned qr_i = 0; qr_i < quregs.size(); ++qr_i){ - res[qr_i] = 0; - for (unsigned qb_i = 0; qb_i < quregs[qr_i].size(); ++qb_i) - res[qr_i] |= ((i >> quregs[qr_i][qb_i])&1) << qb_i; - } - f(res); - auto new_i = i; - for (unsigned qr_i = 0; qr_i < quregs.size(); ++qr_i){ - for (unsigned qb_i = 0; qb_i < quregs[qr_i].size(); ++qb_i){ - if (!(((new_i >> quregs[qr_i][qb_i])&1) == ((res[qr_i] >> qb_i)&1))) - new_i ^= (1UL << quregs[qr_i][qb_i]); - } - } - newvec[new_i] += vec_[i]; - } - else - newvec[i] += vec_[i]; + StateVector newvec; // avoid costly memory reallocations + if( tmpBuff1_.capacity() >= vec_.size() ) + std::swap(newvec, tmpBuff1_); + newvec.resize(vec_.size()); +#pragma omp parallel for schedule(static) + for (std::size_t i = 0; i < vec_.size(); i++) + newvec[i] = 0; + +//#pragma omp parallel reduction(+:newvec[:newvec.size()]) if(parallelize) // requires OpenMP 4.5 + { + std::vector res(quregs.size()); + //#pragma omp for schedule(static) + for (std::size_t i = 0; i < vec_.size(); ++i){ + if ((ctrlmask&i) == ctrlmask){ + for (unsigned qr_i = 0; qr_i < quregs.size(); ++qr_i){ + res[qr_i] = 0; + for (unsigned qb_i = 0; qb_i < quregs[qr_i].size(); ++qb_i) + res[qr_i] |= ((i >> quregs[qr_i][qb_i])&1) << qb_i; + } + f(res); + auto new_i = i; + for (unsigned qr_i = 0; qr_i < quregs.size(); ++qr_i){ + for (unsigned qb_i = 0; qb_i < quregs[qr_i].size(); ++qb_i){ + if (!(((new_i >> quregs[qr_i][qb_i])&1) == ((res[qr_i] >> qb_i)&1))) + new_i ^= (1UL << quregs[qr_i][qb_i]); + } + } + newvec[new_i] += vec_[i]; + } + else + newvec[i] += vec_[i]; + } } - vec_ = std::move(newvec); + std::swap(vec_, newvec); + std::swap(tmpBuff1_, newvec); + } + + // faster version without calling python + template + inline void emulate_math_addConstant(int a, const QuReg& quregs, const std::vector& ctrl) + { + emulate_math([a](std::vector &res){for(auto& x: res) x = x + a;}, quregs, ctrl, true); + } + + // faster version without calling python + template + inline void emulate_math_addConstantModN(int a, int N, const QuReg& quregs, const std::vector& ctrl) + { + emulate_math([a,N](std::vector &res){for(auto& x: res) x = (x + a) % N;}, quregs, ctrl, true); + } + + // faster version without calling python + template + inline void emulate_math_multiplyByConstantModN(int a, int N, const QuReg& quregs, const std::vector& ctrl) + { + emulate_math([a,N](std::vector &res){for(auto& x: res) x = (x * a) % N;}, quregs, ctrl, true); } calc_type get_expectation_value(TermsDict const& td, std::vector const& ids){ run(); calc_type expectation = 0.; - auto current_state = vec_; + + StateVector current_state; // avoid costly memory reallocations + if( tmpBuff1_.capacity() >= vec_.size() ) + std::swap(tmpBuff1_, current_state); + current_state.resize(vec_.size()); +#pragma omp parallel for schedule(static) + for (std::size_t i = 0; i < vec_.size(); ++i) + current_state[i] = vec_[i]; + for (auto const& term : td){ auto const& coefficient = term.second; apply_term(term.first, ids, {}); @@ -260,17 +312,29 @@ class Simulator{ auto const a2 = std::real(vec_[i]); auto const b2 = std::imag(vec_[i]); delta += a1 * a2 - b1 * b2; + // reset vec_ + vec_[i] = current_state[i]; } expectation += coefficient * delta; - vec_ = current_state; } + std::swap(current_state, tmpBuff1_); return expectation; } void apply_qubit_operator(ComplexTermsDict const& td, std::vector const& ids){ run(); - auto new_state = StateVector(vec_.size(), 0.); - auto current_state = vec_; + StateVector new_state, current_state; // avoid costly memory reallocations + if( tmpBuff1_.capacity() >= vec_.size() ) + std::swap(tmpBuff1_, new_state); + if( tmpBuff2_.capacity() >= vec_.size() ) + std::swap(tmpBuff2_, current_state); + new_state.resize(vec_.size()); + current_state.resize(vec_.size()); +#pragma omp parallel for schedule(static) + for (std::size_t i = 0; i < vec_.size(); ++i){ + new_state[i] = 0; + current_state[i] = vec_[i]; + } for (auto const& term : td){ auto const& coefficient = term.second; apply_term(term.first, ids, {}); @@ -280,7 +344,9 @@ class Simulator{ vec_[i] = current_state[i]; } } - vec_ = std::move(new_state); + std::swap(vec_, new_state); + std::swap(tmpBuff1_, new_state); + std::swap(tmpBuff2_, current_state); } calc_type get_probability(std::vector const& bit_string, @@ -452,6 +518,8 @@ class Simulator{ #pragma omp parallel kernel(vec_, ids[4], ids[3], ids[2], ids[1], ids[0], m, ctrlmask); break; + default: + throw std::invalid_argument("Gates with more than 5 qubits are not supported!"); } fused_gates_ = Fusion(); @@ -500,6 +568,12 @@ class Simulator{ unsigned fusion_qubits_min_, fusion_qubits_max_; RndEngine rnd_eng_; std::function rng_; + + // large array buffers to avoid costly reallocations + static StateVector tmpBuff1_, tmpBuff2_; }; +Simulator::StateVector Simulator::tmpBuff1_; +Simulator::StateVector Simulator::tmpBuff2_; + #endif diff --git a/projectq/backends/_sim/_cppsim.cpp b/projectq/backends/_sim/_cppsim.cpp index 74498d4e2..cab68d0ee 100755 --- a/projectq/backends/_sim/_cppsim.cpp +++ b/projectq/backends/_sim/_cppsim.cpp @@ -50,6 +50,9 @@ PYBIND11_PLUGIN(_cppsim) { .def("measure_qubits", &Simulator::measure_qubits_return) .def("apply_controlled_gate", &Simulator::apply_controlled_gate) .def("emulate_math", &emulate_math_wrapper) + .def("emulate_math_addConstant", &Simulator::emulate_math_addConstant) + .def("emulate_math_addConstantModN", &Simulator::emulate_math_addConstantModN) + .def("emulate_math_multiplyByConstantModN", &Simulator::emulate_math_multiplyByConstantModN) .def("get_expectation_value", &Simulator::get_expectation_value) .def("apply_qubit_operator", &Simulator::apply_qubit_operator) .def("emulate_time_evolution", &Simulator::emulate_time_evolution) diff --git a/projectq/backends/_sim/_simulator.py b/projectq/backends/_sim/_simulator.py index 2218c3471..19e884d6d 100755 --- a/projectq/backends/_sim/_simulator.py +++ b/projectq/backends/_sim/_simulator.py @@ -33,10 +33,12 @@ TimeEvolution) from projectq.types import WeakQubitRef +FALLBACK_TO_PYSIM = False try: from ._cppsim import Simulator as SimulatorBackend except ImportError: from ._pysim import Simulator as SimulatorBackend + FALLBACK_TO_PYSIM = True class Simulator(BasicEngine): @@ -384,14 +386,34 @@ def _handle(self, cmd): ID = cmd.qubits[0][0].id self._simulator.deallocate_qubit(ID) elif isinstance(cmd.gate, BasicMathGate): + # improve performance by using C++ code for some commomn gates + from projectq.libs.math import (AddConstant, + AddConstantModN, + MultiplyByConstantModN) qubitids = [] for qr in cmd.qubits: qubitids.append([]) for qb in qr: qubitids[-1].append(qb.id) - math_fun = cmd.gate.get_math_function(cmd.qubits) - self._simulator.emulate_math(math_fun, qubitids, - [qb.id for qb in cmd.control_qubits]) + if FALLBACK_TO_PYSIM: + math_fun = cmd.gate.get_math_function(cmd.qubits) + self._simulator.emulate_math(math_fun, qubitids, + [qb.id for qb in cmd.control_qubits]) + else: + # individual code for different standard gates to make it faster! + if isinstance(cmd.gate, AddConstant): + self._simulator.emulate_math_addConstant(cmd.gate.a, qubitids, + [qb.id for qb in cmd.control_qubits]) + elif isinstance(cmd.gate, AddConstantModN): + self._simulator.emulate_math_addConstantModN(cmd.gate.a, cmd.gate.N, qubitids, + [qb.id for qb in cmd.control_qubits]) + elif isinstance(cmd.gate, MultiplyByConstantModN): + self._simulator.emulate_math_multiplyByConstantModN(cmd.gate.a, cmd.gate.N, qubitids, + [qb.id for qb in cmd.control_qubits]) + else: + math_fun = cmd.gate.get_math_function(cmd.qubits) + self._simulator.emulate_math(math_fun, qubitids, + [qb.id for qb in cmd.control_qubits]) elif isinstance(cmd.gate, TimeEvolution): op = [(list(term), coeff) for (term, coeff) in cmd.gate.hamiltonian.terms.items()] diff --git a/projectq/backends/_sim/_simulator_test.py b/projectq/backends/_sim/_simulator_test.py index 9e301be0d..9f7d298cb 100755 --- a/projectq/backends/_sim/_simulator_test.py +++ b/projectq/backends/_sim/_simulator_test.py @@ -683,3 +683,55 @@ def receive(command_list): qubit1[0].id: qubit0[0].id} assert (sim._convert_logical_to_mapped_qureg(qubit0 + qubit1) == qubit1 + qubit0) + + +def test_simulator_constant_math_emulation(): + if "cpp_simulator" not in get_available_simulators(): + pytest.skip("No C++ simulator") + return + + results = [[[1, 1, 0, 0, 0]], [[0, 1, 0, 0, 0]], [[0, 1, 1, 1, 0]]] + + import projectq.backends._sim._simulator as _sim + from projectq.backends._sim._pysim import Simulator as PySim + from projectq.backends._sim._cppsim import Simulator as CppSim + from projectq.libs.math import (AddConstant, AddConstantModN, + MultiplyByConstantModN) + + def gate_filter(eng, cmd): + g = cmd.gate + if isinstance(g, BasicMathGate): + return False + return eng.next_engine.is_available(cmd) + + def run_simulation(sim): + eng = MainEngine(sim, []) + quint = eng.allocate_qureg(5) + AddConstant(3) | quint + All(Measure) | quint + eng.flush() + results[0].append([int(qb) for qb in quint]) + + AddConstantModN(4, 5) | quint + All(Measure) | quint + eng.flush() + results[1].append([int(qb) for qb in quint]) + + MultiplyByConstantModN(15, 16) | quint + All(Measure) | quint + eng.flush() + results[2].append([int(qb) for qb in quint]) + + cppsim = Simulator(gate_fusion=False) + cppsim._simulator = CppSim(1) + run_simulation(cppsim) + + _sim.FALLBACK_TO_PYSIM = True + pysim = Simulator() + pysim._simulator = PySim(1) + # run_simulation(pysim) + + for result in results: + ref = result[0] + for res in result[1:]: + assert ref == res diff --git a/projectq/libs/math/_gates.py b/projectq/libs/math/_gates.py index fe1df6784..ef5cade99 100755 --- a/projectq/libs/math/_gates.py +++ b/projectq/libs/math/_gates.py @@ -90,6 +90,14 @@ class AddConstantModN(BasicMathGate): qunum = eng.allocate_qureg(5) # 5-qubit number X | qunum[1] # qunum is now equal to 2 AddConstantModN(3, 4) | qunum # qunum is now equal to 1 + + .. note:: + + Pre-conditions: + + * c < N + * c >= 0 + * The value stored in the quantum register must be lower than N """ def __init__(self, a, N): """ @@ -145,6 +153,14 @@ def SubConstantModN(a, N): qunum = eng.allocate_qureg(3) # 3-qubit number X | qunum[1] # qunum is now equal to 2 SubConstantModN(4,5) | qunum # qunum is now -2 = 6 = 1 (mod 5) + + .. note:: + + Pre-conditions: + + * c < N + * c >= 0 + * The value stored in the quantum register must be lower than N """ return AddConstantModN(N - a, N) @@ -162,6 +178,15 @@ class MultiplyByConstantModN(BasicMathGate): qunum = eng.allocate_qureg(5) # 5-qubit number X | qunum[2] # qunum is now equal to 4 MultiplyByConstantModN(3,5) | qunum # qunum is now 2. + + .. note:: + + Pre-conditions: + + * c < N + * c >= 0 + * gcd(c, N) == 1 + * The value stored in the quantum register must be lower than N """ def __init__(self, a, N): """ diff --git a/projectq/ops/__init__.py b/projectq/ops/__init__.py index db4a38b79..cac384d9e 100755 --- a/projectq/ops/__init__.py +++ b/projectq/ops/__init__.py @@ -37,3 +37,5 @@ from ._uniformly_controlled_rotation import (UniformlyControlledRy, UniformlyControlledRz) from ._state_prep import StatePreparation +from ._qpegate import QPE +from ._qaagate import QAA diff --git a/projectq/ops/_gates.py b/projectq/ops/_gates.py index d524e4b3b..be2240d00 100755 --- a/projectq/ops/_gates.py +++ b/projectq/ops/_gates.py @@ -16,17 +16,29 @@ Contains definitions of standard gates such as * Hadamard (H) * Pauli-X (X / NOT) +* Pauli-Y (Y) * Pauli-Z (Z) +* S and its inverse (S / Sdagger) * T and its inverse (T / Tdagger) +* SqrtX gate (SqrtX) * Swap gate (Swap) +* SqrtSwap gate (SqrtSwap) +* Entangle (Entangle) * Phase gate (Ph) +* Rotation-X (Rx) +* Rotation-Y (Ry) * Rotation-Z (Rz) +* Rotation-XX on two qubits (Rxx) +* Rotation-YY on two qubits (Ryy) +* Rotation-ZZ on two qubits (Rzz) * Phase-shift (R) * Measurement (Measure) and meta gates, i.e., * Allocate / Deallocate qubits * Flush gate (end of circuit) +* Barrier +* FlipBits """ import math @@ -110,7 +122,7 @@ def __str__(self): #: Shortcut (instance of) :class:`projectq.ops.SGate` S = SGate() -#: Shortcut (instance of) :class:`projectq.ops.SGate` +#: Inverse (and shortcut) of :class:`projectq.ops.SGate` Sdag = Sdagger = get_inverse(S) @@ -125,7 +137,7 @@ def __str__(self): #: Shortcut (instance of) :class:`projectq.ops.TGate` T = TGate() -#: Shortcut (instance of) :class:`projectq.ops.TGate` +#: Inverse (and shortcut) of :class:`projectq.ops.TGate` Tdag = Tdagger = get_inverse(T) @@ -145,10 +157,9 @@ def __str__(self): SqrtX = SqrtXGate() -class SwapGate(SelfInverseGate, BasicMathGate): +class SwapGate(SelfInverseGate): """ Swap gate class (swaps 2 qubits) """ def __init__(self): - BasicMathGate.__init__(self, lambda x, y: (y, x)) SelfInverseGate.__init__(self) self.interchangeable_qubit_indices = [[0, 1]] @@ -217,7 +228,7 @@ def matrix(self): class Ry(BasicRotationGate): - """ RotationX gate class """ + """ RotationY gate class """ @property def matrix(self): return np.matrix([[math.cos(0.5 * self.angle), @@ -234,6 +245,36 @@ def matrix(self): [0, cmath.exp(.5 * 1j * self.angle)]]) +class Rxx(BasicRotationGate): + """ RotationXX gate class """ + @property + def matrix(self): + return np.matrix([[cmath.cos(.5 * self.angle), 0, 0, -1j*cmath.sin(.5 * self.angle)], + [0, cmath.cos( .5 * self.angle), -1j*cmath.sin(.5 * self.angle), 0], + [0, -1j*cmath.sin(.5 * self.angle), cmath.cos( .5 * self.angle), 0], + [-1j*cmath.sin(.5 * self.angle), 0, 0, cmath.cos( .5 * self.angle)]]) + + +class Ryy(BasicRotationGate): + """ RotationYY gate class """ + @property + def matrix(self): + return np.matrix([[cmath.cos(.5 * self.angle), 0, 0, 1j*cmath.sin(.5 * self.angle)], + [0, cmath.cos( .5 * self.angle), -1j*cmath.sin(.5 * self.angle), 0], + [0, -1j*cmath.sin(.5 * self.angle), cmath.cos( .5 * self.angle), 0], + [1j*cmath.sin(.5 * self.angle), 0, 0, cmath.cos( .5 * self.angle)]]) + + +class Rzz(BasicRotationGate): + """ RotationZZ gate class """ + @property + def matrix(self): + return np.matrix([[cmath.exp(-.5 * 1j * self.angle), 0, 0, 0], + [0, cmath.exp( .5 * 1j * self.angle), 0, 0], + [0, 0, cmath.exp( .5 * 1j * self.angle), 0], + [0, 0, 0, cmath.exp(-.5 * 1j * self.angle)]]) + + class R(BasicPhaseGate): """ Phase-shift gate (equivalent to Rz up to a global phase) """ @property diff --git a/projectq/ops/_gates_test.py b/projectq/ops/_gates_test.py index efcd63b0a..88efa3a19 100755 --- a/projectq/ops/_gates_test.py +++ b/projectq/ops/_gates_test.py @@ -156,6 +156,42 @@ def test_rz(angle): assert np.allclose(gate.matrix, expected_matrix) +@pytest.mark.parametrize("angle", [0, 0.2, 2.1, 4.1, 2 * math.pi, + 4 * math.pi]) +def test_rxx(angle): + gate = _gates.Rxx(angle) + expected_matrix = np.matrix([[cmath.cos(.5 * angle), 0, 0, -1j * cmath.sin(.5 * angle)], + [0, cmath.cos(.5 * angle), -1j * cmath.sin(.5 * angle), 0], + [0, -1j * cmath.sin(.5 * angle), cmath.cos(.5 * angle), 0], + [-1j * cmath.sin(.5 * angle), 0, 0, cmath.cos(.5 * angle)]]) + assert gate.matrix.shape == expected_matrix.shape + assert np.allclose(gate.matrix, expected_matrix) + + +@pytest.mark.parametrize("angle", [0, 0.2, 2.1, 4.1, 2 * math.pi, + 4 * math.pi]) +def test_ryy(angle): + gate = _gates.Ryy(angle) + expected_matrix = np.matrix([[cmath.cos(.5 * angle), 0, 0, 1j * cmath.sin(.5 * angle)], + [0, cmath.cos(.5 * angle), -1j * cmath.sin(.5 * angle), 0], + [0, -1j * cmath.sin(.5 * angle), cmath.cos(.5 * angle), 0], + [ 1j * cmath.sin(.5 * angle), 0, 0, cmath.cos(.5 * angle)]]) + assert gate.matrix.shape == expected_matrix.shape + assert np.allclose(gate.matrix, expected_matrix) + + +@pytest.mark.parametrize("angle", [0, 0.2, 2.1, 4.1, 2 * math.pi, + 4 * math.pi]) +def test_rzz(angle): + gate = _gates.Rzz(angle) + expected_matrix = np.matrix([[cmath.exp(-.5 * 1j * angle), 0, 0, 0], + [0, cmath.exp( .5 * 1j * angle), 0, 0], + [0, 0, cmath.exp( .5 * 1j * angle), 0], + [0, 0, 0, cmath.exp(-.5 * 1j * angle)]]) + assert gate.matrix.shape == expected_matrix.shape + assert np.allclose(gate.matrix, expected_matrix) + + @pytest.mark.parametrize("angle", [0, 0.2, 2.1, 4.1, 2 * math.pi]) def test_ph(angle): gate = _gates.Ph(angle) diff --git a/projectq/ops/_qaagate.py b/projectq/ops/_qaagate.py new file mode 100755 index 000000000..751b9dc60 --- /dev/null +++ b/projectq/ops/_qaagate.py @@ -0,0 +1,81 @@ +# Copyright 2019 ProjectQ-Framework (www.projectq.ch) +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +from ._basics import BasicGate + + +class QAA(BasicGate): + """ + Quantum Aplitude Aplification gate. + + (Quick reference https://en.wikipedia.org/wiki/Amplitude_amplification. + Complete reference G. Brassard, P. Hoyer, M. Mosca, A. Tapp (2000) + Quantum Amplitude Amplification and Estimation + https://arxiv.org/abs/quant-ph/0005055) + + Quantum Amplitude Amplification (QAA) executes the algorithm, but not + the final measurement required to obtain the marked state(s) with high + probability. The starting state on wich the QAA algorithm is executed + is the one resulting of aplying the Algorithm on the |0> state. + + Example: + .. code-block:: python + + def func_algorithm(eng,system_qubits): + All(H) | system_qubits + + def func_oracle(eng,system_qubits,qaa_ancilla): + # This oracle selects the state |010> as the one marked + with Compute(eng): + All(X) | system_qubits[0::2] + with Control(eng, system_qubits): + X | qaa_ancilla + Uncompute(eng) + + system_qubits = eng.allocate_qureg(3) + # Prepare the qaa_ancilla qubit in the |-> state + qaa_ancilla = eng.allocate_qubit() + X | qaa_ancilla + H | qaa_ancilla + + # Creates the initial state form the Algorithm + func_algorithm(eng, system_qubits) + # Apply Quantum Amplitude Amplification the correct number of times + num_it = int(math.pi/4.*math.sqrt(1 << 3)) + with Loop(eng, num_it): + QAA(func_algorithm, func_oracle) | (system_qubits, qaa_ancilla) + + All(Measure) | system_qubits + + Warning: + No qubit allocation/deallocation may take place during the call + to the defined Algorithm :code:`func_algorithm` + + Attributes: + func_algorithm: Algorithm that initialite the state and to be used + in the QAA algorithm + func_oracle: The Oracle that marks the state(s) as "good" + system_qubits: the system we are interested on + qaa_ancilla: auxiliary qubit that helps to invert the amplitude of the + "good" states + + """ + def __init__(self, algorithm, oracle): + BasicGate.__init__(self) + self.algorithm = algorithm + self.oracle = oracle + + def __str__(self): + return 'QAA(Algorithm = {0}, Oracle = {1})'.format( + str(self.algorithm.__name__), str(self.oracle.__name__)) diff --git a/projectq/ops/_qaagate_test.py b/projectq/ops/_qaagate_test.py new file mode 100755 index 000000000..3e15e6801 --- /dev/null +++ b/projectq/ops/_qaagate_test.py @@ -0,0 +1,27 @@ +# Copyright 2019 ProjectQ-Framework (www.projectq.ch) +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Tests for projectq.ops._qaagate.""" + +from projectq.ops import _qaagate, All, H, X + + +def test_qaa_str(): + + def func_algorithm(): All(H) + + def func_oracle(): All(X) + + gate = _qaagate.QAA(func_algorithm, func_oracle) + assert str(gate) == "QAA(Algorithm = func_algorithm, Oracle = func_oracle)" diff --git a/projectq/ops/_qpegate.py b/projectq/ops/_qpegate.py new file mode 100755 index 000000000..08beee743 --- /dev/null +++ b/projectq/ops/_qpegate.py @@ -0,0 +1,29 @@ +# Copyright 2017 ProjectQ-Framework (www.projectq.ch) +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +from ._basics import BasicGate + + +class QPE(BasicGate): + """ + Quantum Phase Estimation gate. + + See setups.decompositions for the complete implementation + """ + def __init__(self, unitary): + BasicGate.__init__(self) + self.unitary = unitary + + def __str__(self): + return 'QPE({})'.format(str(self.unitary)) diff --git a/projectq/ops/_qpegate_test.py b/projectq/ops/_qpegate_test.py new file mode 100755 index 000000000..5ffcbf185 --- /dev/null +++ b/projectq/ops/_qpegate_test.py @@ -0,0 +1,23 @@ +# Copyright 2017 ProjectQ-Framework (www.projectq.ch) +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Tests for projectq.ops._qpegate.""" + +from projectq.ops import _qpegate, X + + +def test_qpe_str(): + unitary = X + gate = _qpegate.QPE(unitary) + assert str(gate) == "QPE(X)" diff --git a/projectq/setups/decompositions/__init__.py b/projectq/setups/decompositions/__init__.py index aab71b28c..883f04581 100755 --- a/projectq/setups/decompositions/__init__.py +++ b/projectq/setups/decompositions/__init__.py @@ -31,7 +31,9 @@ swap2cnot, toffoli2cnotandtgate, time_evolution, - uniformlycontrolledr2cnot) + uniformlycontrolledr2cnot, + phaseestimation, + amplitudeamplification) all_defined_decomposition_rules = [ rule @@ -54,6 +56,8 @@ swap2cnot, toffoli2cnotandtgate, time_evolution, - uniformlycontrolledr2cnot] + uniformlycontrolledr2cnot, + phaseestimation, + amplitudeamplification] for rule in module.all_defined_decomposition_rules ] diff --git a/projectq/setups/decompositions/amplitudeamplification.py b/projectq/setups/decompositions/amplitudeamplification.py new file mode 100644 index 000000000..517aadeed --- /dev/null +++ b/projectq/setups/decompositions/amplitudeamplification.py @@ -0,0 +1,111 @@ +# Copyright 2019 ProjectQ-Framework (www.projectq.ch) +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +""" +Registers a decomposition for quantum amplitude amplification. + +(Quick reference https://en.wikipedia.org/wiki/Amplitude_amplification. +Complete reference G. Brassard, P. Hoyer, M. Mosca, A. Tapp (2000) +Quantum Amplitude Amplification and Estimation +https://arxiv.org/abs/quant-ph/0005055) + +Quantum Amplitude Amplification (QAA) executes the algorithm, but not +the final measurement required to obtain the marked state(s) with high +probability. The starting state on wich the QAA algorithm is executed +is the one resulting of aplying the Algorithm on the |0> state. + +Example: + .. code-block:: python + + def func_algorithm(eng,system_qubits): + All(H) | system_qubits + + def func_oracle(eng,system_qubits,qaa_ancilla): + # This oracle selects the state |010> as the one marked + with Compute(eng): + All(X) | system_qubits[0::2] + with Control(eng, system_qubits): + X | qaa_ancilla + Uncompute(eng) + + system_qubits = eng.allocate_qureg(3) + # Prepare the qaa_ancilla qubit in the |-> state + qaa_ancilla = eng.allocate_qubit() + X | qaa_ancilla + H | qaa_ancilla + + # Creates the initial state form the Algorithm + func_algorithm(eng, system_qubits) + # Apply Quantum Amplitude Amplification the correct number of times + num_it = int(math.pi/4.*math.sqrt(1 << 3)) + with Loop(eng, num_it): + QAA(func_algorithm, func_oracle) | (system_qubits, qaa_ancilla) + + All(Measure) | system_qubits + +Warning: + No qubit allocation/deallocation may take place during the call + to the defined Algorithm :code:`func_algorithm` + +Attributes: + func_algorithm: Algorithm that initialite the state and to be used + in the QAA algorithm + func_oracle: The Oracle that marks the state(s) as "good" + system_qubits: the system we are interested on + qaa_ancilla: auxiliary qubit that helps to invert the amplitude of the + "good" states + +""" + +import math +import numpy as np + +from projectq.cengines import DecompositionRule +from projectq.meta import Control, Compute, Uncompute, CustomUncompute, Dagger +from projectq.ops import X, Z, Ph, All + +from projectq.ops import QAA + + +def _decompose_QAA(cmd): + """ Decompose the Quantum Amplitude Apmplification algorithm as a gate. """ + eng = cmd.engine + + # System-qubit is the first qubit/qureg. Ancilla qubit is the second qubit + system_qubits = cmd.qubits[0] + qaa_ancilla = cmd.qubits[1] + + # The Oracle and the Algorithm + Oracle = cmd.gate.oracle + A = cmd.gate.algorithm + + # Apply the oracle to invert the amplitude of the good states, S_Chi + Oracle(eng, system_qubits, qaa_ancilla) + + # Apply the inversion of the Algorithm, + # the inversion of the aplitude of |0> and the Algorithm + + with Compute(eng): + with Dagger(eng): + A(eng, system_qubits) + All(X) | system_qubits + with Control(eng, system_qubits[0:-1]): + Z | system_qubits[-1] + with CustomUncompute(eng): + All(X) | system_qubits + A(eng, system_qubits) + Ph(math.pi) | system_qubits[0] + + +#: Decomposition rules +all_defined_decomposition_rules = [DecompositionRule(QAA, _decompose_QAA)] diff --git a/projectq/setups/decompositions/amplitudeamplification_test.py b/projectq/setups/decompositions/amplitudeamplification_test.py new file mode 100644 index 000000000..f99681713 --- /dev/null +++ b/projectq/setups/decompositions/amplitudeamplification_test.py @@ -0,0 +1,182 @@ +# Copyright 2019 ProjectQ-Framework (www.projectq.ch) +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"Tests for projectq.setups.decompositions.amplitudeamplification.py." + +import math +import pytest + +from projectq import MainEngine +from projectq.backends import Simulator +from projectq.cengines import (AutoReplacer, DecompositionRuleSet, MainEngine) + +from projectq.ops import (X, H, Ry, All, Measure) +from projectq.meta import Loop, Control, Compute, Uncompute + +from projectq.ops import QAA +from projectq.setups.decompositions import amplitudeamplification as aa + + +def hache_algorithm(eng, qreg): + All(H) | qreg + + +def simple_oracle(eng, system_q, control): + # This oracle selects the state |1010101> as the one marked + with Compute(eng): + All(X) | system_q[1::2] + with Control(eng, system_q): + X | control + Uncompute(eng) + + +def test_simple_grover(): + rule_set = DecompositionRuleSet(modules=[aa]) + + eng = MainEngine(backend=Simulator(), + engine_list=[ + AutoReplacer(rule_set), + ]) + + system_qubits = eng.allocate_qureg(7) + + # Prepare the control qubit in the |-> state + control = eng.allocate_qubit() + X | control + H | control + + # Creates the initial state form the Algorithm + hache_algorithm(eng, system_qubits) + + # Get the amplitude of the marked state before the AA + # to calculate the number of iterations + eng.flush() + prob1010101 = eng.backend.get_probability('1010101', system_qubits) + + total_amp_before = math.sqrt(prob1010101) + theta_before = math.asin(total_amp_before) + + # Apply Quantum Amplitude Amplification the correct number of times + # Theta is calculated previously using get_probability + # We calculate also the theoretical final probability + # of getting the good state + num_it = int(math.pi / (4. * theta_before) + 1) + theoretical_prob = math.sin((2 * num_it + 1.) * theta_before)**2 + with Loop(eng, num_it): + QAA(hache_algorithm, simple_oracle) | (system_qubits, control) + + # Get the probabilty of getting the marked state after the AA + # to compare with the theoretical probability after teh AA + eng.flush() + prob1010101 = eng.backend.get_probability('1010101', system_qubits) + total_prob_after = prob1010101 + + All(Measure) | system_qubits + H | control + Measure | control + result = [int(q) for q in system_qubits] + control_result = int(control) + + eng.flush() + + assert total_prob_after == pytest.approx(theoretical_prob, abs=1e-6), ( + "The obtained probability is less than expected %f vs. %f" % + (total_prob_after, theoretical_prob)) + + +def complex_algorithm(eng, qreg): + All(H) | qreg + with Control(eng, qreg[0]): + All(X) | qreg[1:] + All(Ry(math.pi / 4)) | qreg[1:] + with Control(eng, qreg[-1]): + All(X) | qreg[1:-1] + + +def complex_oracle(eng, system_q, control): + # This oracle selects the subspace |000000>+|111111> as the good one + with Compute(eng): + with Control(eng, system_q[0]): + All(X) | system_q[1:] + H | system_q[0] + All(X) | system_q + + with Control(eng, system_q): + X | control + + Uncompute(eng) + + +def test_complex_aa(): + rule_set = DecompositionRuleSet(modules=[aa]) + + eng = MainEngine(backend=Simulator(), + engine_list=[ + AutoReplacer(rule_set), + ]) + + system_qubits = eng.allocate_qureg(6) + + # Prepare the control qubit in the |-> state + control = eng.allocate_qubit() + X | control + H | control + + # Creates the initial state form the Algorithm + complex_algorithm(eng, system_qubits) + + # Get the probabilty of getting the marked state before the AA + # to calculate the number of iterations + eng.flush() + prob000000 = eng.backend.get_probability('000000', system_qubits) + prob111111 = eng.backend.get_probability('111111', system_qubits) + + total_amp_before = math.sqrt(prob000000 + prob111111) + theta_before = math.asin(total_amp_before) + + # Apply Quantum Amplitude Amplification the correct number of times + # Theta is calculated previously using get_probability + # We calculate also the theoretical final probability + # of getting the good state + num_it = int(math.pi / (4. * theta_before) + 1) + theoretical_prob = math.sin((2 * num_it + 1.) * theta_before)**2 + with Loop(eng, num_it): + QAA(complex_algorithm, complex_oracle) | (system_qubits, control) + + # Get the probabilty of getting the marked state after the AA + # to compare with the theoretical probability after the AA + eng.flush() + prob000000 = eng.backend.get_probability('000000', system_qubits) + prob111111 = eng.backend.get_probability('111111', system_qubits) + total_prob_after = prob000000 + prob111111 + + All(Measure) | system_qubits + H | control + Measure | control + result = [int(q) for q in system_qubits] + control_result = int(control) + + eng.flush() + + assert total_prob_after == pytest.approx(theoretical_prob, abs=1e-2), ( + "The obtained probability is less than expected %f vs. %f" % + (total_prob_after, theoretical_prob)) + + +def test_string_functions(): + algorithm = hache_algorithm + oracle = simple_oracle + gate = QAA(algorithm, oracle) + assert (str(gate) == + "QAA(Algorithm = hache_algorithm, Oracle = simple_oracle)") diff --git a/projectq/setups/decompositions/phaseestimation.py b/projectq/setups/decompositions/phaseestimation.py new file mode 100644 index 000000000..faf7523cf --- /dev/null +++ b/projectq/setups/decompositions/phaseestimation.py @@ -0,0 +1,129 @@ +# Copyright 2018 ProjectQ-Framework (www.projectq.ch) +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +""" +Registers a decomposition for phase estimation. + +(reference https://en.wikipedia.org/wiki/Quantum_phase_estimation_algorithm) + +The Quantum Phase Estimation (QPE) executes the algorithm up to the inverse +QFT included. The following steps measuring the ancillas and computing the +phase should be executed outside of the QPE. + +The decomposition uses as ancillas (qpe_ancillas) the first qubit/qureg in +the Command and as system qubits teh second qubit/qureg in the Command. + +The unitary operator for which the phase estimation is estimated (unitary) +is the gate in Command + +Example: + .. code-block:: python + + # Example using a ProjectQ gate + + n_qpe_ancillas = 3 + qpe_ancillas = eng.allocate_qureg(n_qpe_ancillas) + system_qubits = eng.allocate_qureg(1) + angle = cmath.pi*2.*0.125 + U = Ph(angle) # unitary_specfic_to_the_problem() + + # Apply Quantum Phase Estimation + QPE(U) | (qpe_ancillas, system_qubits) + + All(Measure) | qpe_ancillas + # Compute the phase from the ancilla measurement + #(https://en.wikipedia.org/wiki/Quantum_phase_estimation_algorithm) + phasebinlist = [int(q) for q in qpe_ancillas] + phase_in_bin = ''.join(str(j) for j in phasebinlist) + phase_int = int(phase_in_bin,2) + phase = phase_int / (2 ** n_qpe_ancillas) + print (phase) + + # Example using a function (two_qubit_gate). + # Instead of applying QPE on a gate U one could provide a function + + def two_qubit_gate(system_q, time): + CNOT | (system_q[0], system_q[1]) + Ph(2.0*cmath.pi*(time * 0.125)) | system_q[1] + CNOT | (system_q[0], system_q[1]) + + n_qpe_ancillas = 3 + qpe_ancillas = eng.allocate_qureg(n_qpe_ancillas) + system_qubits = eng.allocate_qureg(2) + X | system_qubits[0] + + # Apply Quantum Phase Estimation + QPE(two_qubit_gate) | (qpe_ancillas, system_qubits) + + All(Measure) | qpe_ancillas + # Compute the phase from the ancilla measurement + #(https://en.wikipedia.org/wiki/Quantum_phase_estimation_algorithm) + phasebinlist = [int(q) for q in qpe_ancillas] + phase_in_bin = ''.join(str(j) for j in phasebinlist) + phase_int = int(phase_in_bin,2) + phase = phase_int / (2 ** n_qpe_ancillas) + print (phase) + +Attributes: + unitary (BasicGate): Unitary Operation either a ProjectQ gate or a function f. + Calling the function with the parameters system_qubits(Qureg) and time (integer), + i.e. f(system_qubits, time), applies to the system qubits a unitary defined in f + with parameter time. + + +""" + +import numpy as np + +from projectq.cengines import DecompositionRule +from projectq.meta import Control, Loop, get_control_count +from projectq.ops import H, Tensor, get_inverse, QFT + +from projectq.ops import QPE + + +def _decompose_QPE(cmd): + """ Decompose the Quantum Phase Estimation gate. """ + eng = cmd.engine + + # Ancillas is the first qubit/qureg. System-qubit is the second qubit/qureg + qpe_ancillas = cmd.qubits[0] + system_qubits = cmd.qubits[1] + + # Hadamard on the ancillas + Tensor(H) | qpe_ancillas + + # The Unitary Operator + U = cmd.gate.unitary + + # Control U on the system_qubits + if (callable(U)): + # If U is a function + for i in range(len(qpe_ancillas)): + with Control(eng, qpe_ancillas[i]): + U(system_qubits, time=2**i) + else: + for i in range(len(qpe_ancillas)): + ipower = int(2**i) + with Loop(eng, ipower): + with Control(eng, qpe_ancillas[i]): + U | system_qubits + + # Inverse QFT on the ancillas + get_inverse(QFT) | qpe_ancillas + +#: Decomposition rules +all_defined_decomposition_rules = [ + DecompositionRule(QPE, _decompose_QPE) +] diff --git a/projectq/setups/decompositions/phaseestimation_test.py b/projectq/setups/decompositions/phaseestimation_test.py new file mode 100644 index 000000000..1b8f8e63c --- /dev/null +++ b/projectq/setups/decompositions/phaseestimation_test.py @@ -0,0 +1,162 @@ +# Copyright 2017 ProjectQ-Framework (www.projectq.ch) +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"Tests for projectq.setups.decompositions.phaseestimation.py." + +import cmath +import numpy as np +import pytest + +from projectq import MainEngine +from projectq.backends import Simulator +from projectq.cengines import (AutoReplacer, DecompositionRuleSet, + DummyEngine, InstructionFilter, MainEngine) + +from projectq.ops import X, H, All, Measure, Tensor, Ph, CNOT, StatePreparation + +from projectq.ops import (BasicGate) + +from projectq.ops import QPE +from projectq.setups.decompositions import phaseestimation as pe +from projectq.setups.decompositions import qft2crandhadamard as dqft +import projectq.setups.decompositions.stateprep2cnot as stateprep2cnot +import projectq.setups.decompositions.uniformlycontrolledr2cnot as ucr2cnot + + +def test_simple_test_X_eigenvectors(): + rule_set = DecompositionRuleSet(modules=[pe, dqft]) + eng = MainEngine(backend=Simulator(), + engine_list=[AutoReplacer(rule_set), + ]) + results = np.array([]) + for i in range(100): + autovector = eng.allocate_qureg(1) + X | autovector + H | autovector + unit = X + ancillas = eng.allocate_qureg(1) + QPE(unit) | (ancillas, autovector) + All(Measure) | ancillas + fasebinlist = [int(q) for q in ancillas] + fasebin = ''.join(str(j) for j in fasebinlist) + faseint = int(fasebin, 2) + phase = faseint / (2. ** (len(ancillas))) + results = np.append(results, phase) + All(Measure) | autovector + eng.flush() + + num_phase = (results == 0.5).sum() + assert num_phase/100. >= 0.35, "Statistics phase calculation are not correct (%f vs. %f)" % (num_phase/100., 0.35) + + +def test_Ph_eigenvectors(): + rule_set = DecompositionRuleSet(modules=[pe, dqft]) + eng = MainEngine(backend=Simulator(), + engine_list=[AutoReplacer(rule_set), + ]) + results = np.array([]) + for i in range(100): + autovector = eng.allocate_qureg(1) + theta = cmath.pi*2.*0.125 + unit = Ph(theta) + ancillas = eng.allocate_qureg(3) + QPE(unit) | (ancillas, autovector) + All(Measure) | ancillas + fasebinlist = [int(q) for q in ancillas] + fasebin = ''.join(str(j) for j in fasebinlist) + faseint = int(fasebin, 2) + phase = faseint / (2. ** (len(ancillas))) + results = np.append(results, phase) + All(Measure) | autovector + eng.flush() + + num_phase = (results == 0.125).sum() + assert num_phase/100. >= 0.35, "Statistics phase calculation are not correct (%f vs. %f)" % (num_phase/100., 0.35) + + +def two_qubit_gate(system_q, time): + CNOT | (system_q[0], system_q[1]) + Ph(2.0*cmath.pi*(time * 0.125)) | system_q[1] + CNOT | (system_q[0], system_q[1]) + + +def test_2qubitsPh_andfunction_eigenvectors(): + rule_set = DecompositionRuleSet(modules=[pe, dqft]) + eng = MainEngine(backend=Simulator(), + engine_list=[AutoReplacer(rule_set), + ]) + results = np.array([]) + for i in range(100): + autovector = eng.allocate_qureg(2) + X | autovector[0] + ancillas = eng.allocate_qureg(3) + QPE(two_qubit_gate) | (ancillas, autovector) + All(Measure) | ancillas + fasebinlist = [int(q) for q in ancillas] + fasebin = ''.join(str(j) for j in fasebinlist) + faseint = int(fasebin, 2) + phase = faseint / (2. ** (len(ancillas))) + results = np.append(results, phase) + All(Measure) | autovector + eng.flush() + + num_phase = (results == 0.125).sum() + assert num_phase/100. >= 0.34, "Statistics phase calculation are not correct (%f vs. %f)" % (num_phase/100., 0.34) + + +def test_X_no_eigenvectors(): + rule_set = DecompositionRuleSet(modules=[pe, dqft, stateprep2cnot, ucr2cnot]) + eng = MainEngine(backend=Simulator(), + engine_list=[AutoReplacer(rule_set), + ]) + results = np.array([]) + results_plus = np.array([]) + results_minus = np.array([]) + for i in range(100): + autovector = eng.allocate_qureg(1) + amplitude0 = (np.sqrt(2) + np.sqrt(6))/4. + amplitude1 = (np.sqrt(2) - np.sqrt(6))/4. + StatePreparation([amplitude0, amplitude1]) | autovector + unit = X + ancillas = eng.allocate_qureg(1) + QPE(unit) | (ancillas, autovector) + All(Measure) | ancillas + fasebinlist = [int(q) for q in ancillas] + fasebin = ''.join(str(j) for j in fasebinlist) + faseint = int(fasebin, 2) + phase = faseint / (2. ** (len(ancillas))) + results = np.append(results, phase) + Tensor(H) | autovector + if np.allclose(phase, .0, rtol=1e-1): + results_plus = np.append(results_plus, phase) + All(Measure) | autovector + autovector_result = int(autovector) + assert autovector_result == 0 + elif np.allclose(phase, .5, rtol=1e-1): + results_minus = np.append(results_minus, phase) + All(Measure) | autovector + autovector_result = int(autovector) + assert autovector_result == 1 + eng.flush() + + total = len(results_plus) + len(results_minus) + plus_probability = len(results_plus)/100. + assert total == pytest.approx(100, abs=5) + assert plus_probability == pytest.approx(1./4., abs = 1e-1), "Statistics on |+> probability are not correct (%f vs. %f)" % (plus_probability, 1./4.) + + +def test_string(): + unit = X + gate = QPE(unit) + assert (str(gate) == "QPE(X)") diff --git a/setup.py b/setup.py index 604f006af..a0254b0e1 100755 --- a/setup.py +++ b/setup.py @@ -124,6 +124,7 @@ def build_extensions(self): opts.append('/arch:AVX') else: opts.append('-march=native') + opts.append('-ffast-math') opts.append(openmp) if ct == 'unix':