Repository: rigetticomputing/pyquil
Pull request #251 (2018-01-02 19:24:48)

Port kraus code from docs and add a high-level interface
cc #228 #243 

mpharrigan: Alright, I think I've addressed everything. Ready for re-review

Patch:
diff --git a/examples/NoiseInvestigation.ipynb b/examples/NoiseInvestigation.ipynb
new file mode 100644
index 0000000..8f397f1
--- /dev/null
+++ b/examples/NoiseInvestigation.ipynb
@@ -0,0 +1,234 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Investigating Noise\n",
+    "\n",
+    "In this example, we investigate how a program might behave on a near-term device that is subject to noise using the convenience function `add_noise_to_program`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from pyquil.quil import Program\n",
+    "from pyquil.paulis import PauliSum, PauliTerm, exponentiate, exponential_map, trotterize\n",
+    "from pyquil.gates import MEASURE, H, Z, RX, RZ, CZ\n",
+    "import numpy as np"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The task\n",
+    "We want to prepare $e^{i \\theta XY}$ and measure it in the $Z$ basis."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from numpy import pi\n",
+    "theta = pi/3\n",
+    "xy = PauliTerm('X', 0) * PauliTerm('Y', 1)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### The idiomatic Pyquil program"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "prog = exponential_map(xy)(theta)\n",
+    "print(prog)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### The compiled program\n",
+    "To run on a real device, we must compile each program to the native gate set for the device. The high-level noise model is similarly constrained to use a small, native gate set. In particular, we can use\n",
+    "\n",
+    " - $I$\n",
+    " - $RZ(\\theta)$\n",
+    " - $RX(\\pm \\pi/2)$\n",
+    " - $CZ$\n",
+    "\n",
+    "For simplicity, the compiled program is given below but generally you will want to use a compiler to do this step for you."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def get_compiled_prog(theta):\n",
+    "    return Program([\n",
+    "        RZ(-pi/2, 0),\n",
+    "        RX(-pi/2, 0),\n",
+    "        RZ(-pi/2, 1),\n",
+    "        RX( pi/2, 1),\n",
+    "        CZ(1, 0),\n",
+    "        RZ(-pi/2, 1),\n",
+    "        RX(-pi/2, 1),\n",
+    "        RZ(theta, 1),\n",
+    "        RX( pi/2, 1),\n",
+    "        CZ(1, 0),\n",
+    "        RX( pi/2, 0),\n",
+    "        RZ( pi/2, 0),\n",
+    "        RZ(-pi/2, 1),\n",
+    "        RX( pi/2, 1),\n",
+    "        RZ(-pi/2, 1),\n",
+    "    ])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Scan over noise parameters\n",
+    "We perform a scan over three levels of noise each at 20 theta points.\n",
+    "\n",
+    "Specifically, we investigate T1 values of 1, 3, and 10 us. By default, T2 = T1 / 2, 1 qubit gates take 50 ns, and 2 qubit gates take 150 ns. \n",
+    "\n",
+    "In alignment with the device, $I$ and parametric $RZ$ are noiseless while $RX$ and $CZ$ gates experience 1q and 2q gate noise, respectively."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from pyquil.api import QVMConnection\n",
+    "cxn = QVMConnection()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "t1s = np.logspace(-6, -5, num=3)\n",
+    "thetas = np.linspace(-pi, pi, num=20)\n",
+    "t1s * 1e6 # us"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from pyquil.kraus import add_noise_to_program\n",
+    "records = []\n",
+    "for theta in thetas:\n",
+    "    for t1 in t1s:\n",
+    "        prog = get_compiled_prog(theta)\n",
+    "        noisy = add_noise_to_program(prog, T1=t1).inst([\n",
+    "            MEASURE(0, 0),\n",
+    "            MEASURE(1, 1),\n",
+    "        ])\n",
+    "        bitstrings = np.array(cxn.run(noisy, [0,1], 1000))\n",
+    "        \n",
+    "        # Expectation of Z0 and Z1\n",
+    "        z0, z1 = 1 - 2*np.mean(bitstrings, axis=0)\n",
+    "        \n",
+    "        # Expectation of ZZ by computing the parity of each pair\n",
+    "        zz = 1 - (np.sum(bitstrings, axis=1) % 2).mean() * 2 \n",
+    "        \n",
+    "        record = {\n",
+    "            'z0': z0,\n",
+    "            'z1': z1,\n",
+    "            'zz': zz,\n",
+    "            'theta': theta,\n",
+    "            't1': t1,\n",
+    "        }\n",
+    "        records += [record]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Plot the results"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "from matplotlib import pyplot as plt\n",
+    "import seaborn as sns\n",
+    "sns.set(style='ticks', palette='colorblind')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import pandas as pd\n",
+    "df_all = pd.DataFrame(records)\n",
+    "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12,4))\n",
+    "\n",
+    "for t1 in t1s:\n",
+    "    df = df_all.query('t1 == @t1')\n",
+    "    \n",
+    "    ax1.plot(df['theta'], df['z0'], 'o-')    \n",
+    "    ax2.plot(df['theta'], df['z1'], 'o-')    \n",
+    "    ax3.plot(df['theta'], df['zz'], 'o-', label=f'T1 = {t1*1e6:.0f} us')\n",
+    "    \n",
+    "ax3.legend(loc='best')\n",
+    "\n",
+    "ax1.set_ylabel('Z0')\n",
+    "ax2.set_ylabel('Z1')\n",
+    "ax3.set_ylabel('ZZ')\n",
+    "ax2.set_xlabel(r'$\\theta$')\n",
+    "fig.tight_layout()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "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.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/pyquil/kraus.py b/pyquil/kraus.py
index 4dbe0d5..1d18046 100644
--- a/pyquil/kraus.py
+++ b/pyquil/kraus.py
@@ -1,5 +1,5 @@
 ##############################################################################
-# Copyright 2016-2017 Rigetti Computing
+# Copyright 2018 Rigetti Computing
 #
 #    Licensed under the Apache License, Version 2.0 (the "License");
 #    you may not use this file except in compliance with the License.
@@ -16,10 +16,11 @@
 """
 Module for creating and verifying noisy gate definitions in terms of Kraus maps.
 """
+
 import numpy as np
 
-from pyquil.quilbase import Pragma
 from pyquil.parameters import format_parameter
+from pyquil.quilbase import Pragma, Gate
 
 
 def _check_kraus_ops(n, kraus_ops):
@@ -35,7 +36,7 @@ def _check_kraus_ops(n, kraus_ops):
                 "Kraus operators for {0} qubits must have shape {1}x{1}: {2}".format(n, 2 ** n, k))
 
     kdk_sum = sum(np.transpose(k).conjugate().dot(k) for k in kraus_ops)
-    if not np.allclose(kdk_sum, np.eye(2**n), atol=1e-5):
+    if not np.allclose(kdk_sum, np.eye(2 ** n), atol=1e-5):
         raise ValueError(
             "Kraus operator not correctly normalized: sum_j K_j^*K_j == {}".format(kdk_sum))
 
@@ -56,3 +57,223 @@ def _create_kraus_pragmas(name, qubit_indices, kraus_ops):
                       "({})".format(" ".join(map(format_parameter, np.ravel(k)))))
                for k in kraus_ops]
     return pragmas
+
+
+def append_kraus_to_gate(kraus_ops, gate_matrix):
+    """
+    Follow a gate ``gate_matrix`` by a Kraus map described by ``kraus_ops``.
+
+    :param list kraus_ops: The Kraus operators.
+    :param numpy.ndarray gate_matrix: The unitary gate.
+    :return: A list of transformed Kraus operators.
+    """
+    return [kj.dot(gate_matrix) for kj in kraus_ops]
+
+
+def damping_kraus_map(p=0.10):
+    """
+    Generate the Kraus operators corresponding to an amplitude damping
+    noise channel.
+
+    :param float p: The one-step damping probability.
+    :return: A list [k1, k2] of the Kraus operators that parametrize the map.
+    :rtype: list
+    """
+    damping_op = np.sqrt(p) * np.array([[0, 1],
+                                        [0, 0]])
+
+    residual_kraus = np.diag([1, np.sqrt(1 - p)])
+    return [residual_kraus, damping_op]
+
+
+def dephasing_kraus_map(p=0.10):
+    """
+    Generate the Kraus operators corresponding to a dephasing channel.
+
+    :params float p: The one-step dephasing probability.
+    :return: A list [k1, k2] of the Kraus operators that parametrize the map.
+    :rtype: list
+    """
+    return [np.sqrt(1 - p) * np.eye(2), np.sqrt(p) * np.diag([1, -1])]
+
+
+def tensor_kraus_maps(k1, k2):
+    """
+    Generate the Kraus map corresponding to the composition
+    of two maps on different qubits.
+
+    :param list k1: The Kraus operators for the first qubit.
+    :param list k2: The Kraus operators for the second qubit.
+    :return: A list of tensored Kraus operators.
+    """
+    return [np.kron(k1j, k2l) for k1j in k1 for k2l in k2]
+
+
+def combine_kraus_maps(k1, k2):
+    """
+    Generate the Kraus map corresponding to the composition
+    of two maps on the same qubits with k1 being applied to the state
+    after k2.
+
+    :param list k1: The list of Kraus operators that are applied second.
+    :param list k2: The list of Kraus operators that are applied first.
+    :return: A combinatorially generated list of composed Kraus operators.
+    """
+    return [np.dot(k1j, k2l) for k1j in k1 for k2l in k2]
+
+
+def damping_after_dephasing(T1, T2, gate_time):
+    """
+    Generate the Kraus map corresponding to the composition
+    of a dephasing channel followed by an amplitude damping channel.
+
+    :param float T1: The amplitude damping time
+    :param float T2: The dephasing time
+    :param float gate_time: The gate duration.
+    :return: A list of Kraus operators.
+    """
+    damping = damping_kraus_map(p=gate_time / float(T1))
+    dephasing = dephasing_kraus_map(p=gate_time / float(T2))
+    return combine_kraus_maps(damping, dephasing)
+
+
+# You can only apply gate-noise to non-parametrized gates,
+# so we need to define placeholders for RX(+/- pi/2)
+SINGLE_Q = {
+    "noisy-x-plus90": (np.array([[1, -1j],
+                                 [-1j, 1]]) / np.sqrt(2)),
+    "noisy-x-minus90": (np.array([[1, 1j],
+                                  [1j, 1]]) / np.sqrt(2)),
+}
+
+TWO_Q = {
+    "noisy-cz": np.diag([1, 1, 1, -1]),
+}
+
+
+def _noisy_instruction(instruction):
+    """
+    Translate an ordinary gate instruction into its noisy version, where applicable.
+
+    In an attempt to closely model the QPU, noisy versions of RX(+-pi/2) and CZ are supported;
+    I and parametric RZ are returned as-is (noiseless), and other gates are not allowed.
+    Pragmas are returned as-is.
+
+    Note: this function doesn't actually define the noisy gates. Please see
+    :py:func`add_noise_to_program` for a full solution.
+
+    :param instruction: The instruction
+    :return: A noisy version of the instruction
+    """
+    if not isinstance(instruction, Gate):
+        return instruction
+
+    if instruction.name == 'RZ':
+        return instruction
+
+    if instruction.name == 'I':
+        return instruction
+
+    if instruction.name == 'RX':
+        assert len(instruction.params) == 1
+        assert len(instruction.qubits) == 1
+        if instruction.params[0] == np.pi / 2.0:
+            return Gate('noisy-x-plus90', [], instruction.qubits)
+        if instruction.params[0] == -np.pi / 2.0:
+            return Gate('noisy-x-minus90', [], instruction.qubits)
+        raise ValueError("Can't add noise to a parametric gate. "
+                         "Try compiling to RX(pi/2) or RX(-pi/2)")
+
+    if instruction.name == 'CZ':
+        return Gate('noisy-cz', [], instruction.qubits)
+
+    raise ValueError('Gate {} is not in the native instruction set'.format(instruction))
+
+
+def _get_program_topology(prog):
+    """
+    Get the graph of two-qubit gates used in a program.
+
+    :param prog: The program
+    :return: (qubits, edges). ``qubits`` is a set of all integer qubits.
+        ``edges`` is a set of tuples of two-qubit pairings.
+    """
+    qubits = set()
+    edges = set()
+    for instruction in prog:
+        if not isinstance(instruction, Gate):
+            continue
+        qb = [q.index for q in instruction.qubits]
+        qubits.update(qb)
+        if len(qb) == 2:
+            edges.add(tuple(sorted(qb)))
+        if len(qb) > 2:
+            raise ValueError("Encountered a >2 qubit instruction")
+
+    return qubits, edges
+
+
+def add_noise_to_program(prog, T1=30e-6, T2=None, gate_time_1q=50e-9, gate_time_2q=150e-09):
+    """
+    Add generic damping and dephasing noise to a program.
+
+    This high-level function is provided as a convenience to investigate the effects of a
+    generic noise model on a program. For more fine-grained control, please investigate
+    the other methods available in the ``pyquil.kraus`` module.
+
+    In an attempt to closely model the QPU, noisy versions of RX(+-pi/2) and CZ are provided;
+    I and parametric RZ are noiseless, and other gates are not allowed. To use this function,
+    you need to compile your program to this native gate set.
+
+    The default noise parameters
+
+        T1 = 30 us
+        T2 = T1 / 2
+        1q gate time = 50 ns
+        2q gate time = 150 ns
+
+    are currently typical for near-term devices.
+
+    This function will define new gates and add Kraus noise to these gates. It will translate
+    the input program to use the noisy version of the gates.
+
+    :param prog: A pyquil program consisting of I, RZ, CZ, and RX(+-pi/2) instructions
+    :param T1: The T1 amplitude damping time. By default, this is 30 us
+    :param T2: The T2 dephasing time. By default, this is one-half of the T1 time.
+    :param gate_time_1q: The duration of the one-qubit gates, namely RX(+pi/2) and RX(-pi/2).
+        By default, this is 50 ns.
+    :param gate_time_2q: The duration of the two-qubit gates, namely CZ.
+        By default, this is 150 ns.
+    :return: A new program with noisy operators.
+    """
+
+    if T2 is None:
+        T2 = T1 / 2.0
+
+    from pyquil.quil import Program  # Avoid circular dependency
+    qubits, edges = _get_program_topology(prog)
+    new_prog = Program()
+
+    # Define noisy 1q gates
+    for name, matrix in SINGLE_Q.items():
+        new_prog.defgate(name, matrix)
+        k_ops = append_kraus_to_gate(damping_after_dephasing(T1=T1, T2=T2, gate_time=gate_time_1q),
+                                     gate_matrix=matrix)
+        for qubit in qubits:
+            new_prog.define_noisy_gate(name, (qubit,), k_ops)
+
+    # Define noisy 2q gates
+    for name, matrix in TWO_Q.items():
+        new_prog.defgate(name, matrix)
+        k1q = damping_after_dephasing(T1=T1, T2=T2, gate_time=gate_time_2q)
+        k_total = append_kraus_to_gate(tensor_kraus_maps(k1q, k1q), gate_matrix=matrix)
+        for q1, q2 in edges:
+            # Note! q1 must be less than q2. This is done in _get_program_topology
+            assert q1 < q2
+            new_prog.define_noisy_gate(name, (q1, q2), k_total)
+
+    # Translate noiseless gates to noisy gates.
+    new_instrs = [_noisy_instruction(instruction) for instruction in prog]
+    new_prog.inst(new_instrs)
+
+    return new_prog


Test patch:
diff --git a/pyquil/tests/test_kraus.py b/pyquil/tests/test_kraus.py
new file mode 100644
index 0000000..6974910
--- /dev/null
+++ b/pyquil/tests/test_kraus.py
@@ -0,0 +1,69 @@
+import numpy as np
+import pytest
+
+from pyquil.gates import CZ, RZ, RX, I, H
+from pyquil.kraus import damping_kraus_map, dephasing_kraus_map, tensor_kraus_maps, \
+    combine_kraus_maps, damping_after_dephasing, _noisy_instruction
+from pyquil.quil import Pragma
+
+
+def test_damping_kraus_map():
+    p = 0.05
+    k1, k2 = damping_kraus_map(p=p)
+    assert k1[1, 1] == np.sqrt(1 - p)
+    assert k2[0, 1] == np.sqrt(p)
+
+
+def test_dephasing_kraus_map():
+    p = 0.05
+    k1, k2 = dephasing_kraus_map(p=p)
+    np.testing.assert_allclose(np.diag(k1), [np.sqrt(1 - p)] * 2)
+    np.testing.assert_allclose(np.abs(np.diag(k2)), [np.sqrt(p)] * 2)
+
+
+def test_tensor_kraus_maps():
+    damping = damping_kraus_map()
+    k1, k2, k3, k4 = tensor_kraus_maps(damping, damping)
+    assert k1.shape == (4, 4)
+    assert k2.shape == (4, 4)
+    assert k3.shape == (4, 4)
+    assert k4.shape == (4, 4)
+    np.testing.assert_allclose(k1[-1, -1], 1 - 0.1)
+
+
+def test_combine_kraus_maps():
+    damping = damping_kraus_map()
+    dephasing = dephasing_kraus_map()
+    k1, k2, k3, k4 = combine_kraus_maps(damping, dephasing)
+    assert k1.shape == (2, 2)
+    assert k2.shape == (2, 2)
+    assert k3.shape == (2, 2)
+    assert k4.shape == (2, 2)
+
+
+def test_damping_after_dephasing():
+    damping = damping_kraus_map()
+    dephasing = dephasing_kraus_map()
+    ks_ref = combine_kraus_maps(damping, dephasing)
+
+    ks_actual = damping_after_dephasing(10, 10, 1)
+    np.testing.assert_allclose(ks_actual, ks_ref)
+
+
+def test_noisy_instruction():
+    # Unaffected
+    assert I(0) == _noisy_instruction(I(0))
+    assert RZ(.234)(0) == _noisy_instruction(RZ(.234)(0))
+    assert Pragma('lalala') == _noisy_instruction(Pragma('lalala'))
+
+    # Noisified
+    assert _noisy_instruction(CZ(0, 1)).out() == 'noisy-cz 0 1'
+    assert _noisy_instruction(RX(np.pi / 2)(0)).out() == 'noisy-x-plus90 0'
+    assert _noisy_instruction(RX(-np.pi / 2)(23)).out() == 'noisy-x-minus90 23'
+
+    # Unsupported
+    with pytest.raises(ValueError):
+        _noisy_instruction(H(0))
+
+    with pytest.raises(ValueError):
+        _noisy_instruction(RX(2 * np.pi / 3)(0))