From a3793d762569e630287f26479802856b5f677a2c Mon Sep 17 00:00:00 2001 From: Oscar Higgott Date: Mon, 26 Oct 2020 21:35:24 +0000 Subject: [PATCH] More documentation added --- .gitignore | 3 + docs/conf.py | 3 +- docs/index.rst | 3 +- .../{getting_started.rst => installation.rst} | 11 +- docs/usage.ipynb | 608 ++++++++++++++++++ src/pymatching/matching.py | 59 +- src/pymatching/mwpm.cpp | 15 - 7 files changed, 654 insertions(+), 48 deletions(-) rename docs/{getting_started.rst => installation.rst} (60%) create mode 100644 docs/usage.ipynb diff --git a/.gitignore b/.gitignore index 96b3789c..139dfdd2 100644 --- a/.gitignore +++ b/.gitignore @@ -82,3 +82,6 @@ htmlcov ### Timing ### timing + +### Jupyter ### +.ipynb_checkpoints diff --git a/docs/conf.py b/docs/conf.py index 11ef201a..80d16c62 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -34,7 +34,8 @@ "sphinx.ext.autodoc", 'sphinx.ext.mathjax', 'sphinx.ext.napoleon', - 'sphinx_rtd_theme' + 'sphinx_rtd_theme', + 'nbsphinx' ] # Add any paths that contain templates here, relative to this directory. diff --git a/docs/index.rst b/docs/index.rst index 8ca8e36d..1a0f5766 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -10,7 +10,8 @@ Welcome to PyMatching's documentation! :maxdepth: 2 :caption: Contents: - getting_started + installation + usage.ipynb api diff --git a/docs/getting_started.rst b/docs/installation.rst similarity index 60% rename from docs/getting_started.rst rename to docs/installation.rst index e06a80ae..7a161fee 100644 --- a/docs/getting_started.rst +++ b/docs/installation.rst @@ -1,4 +1,4 @@ -Getting Started +Installation =============== Installing PyMatching @@ -12,12 +12,3 @@ first clone the repository (you must use the ``--recursive`` flag):: and then install using ``pip``:: pip install -e ./PyMatching - - -Usage ------ - -Much of the functionality of PyMatching is available through the -:py:class:`pymatching.matching.Matching` class, which can be imported in Python with:: - - from pymatching import Matching \ No newline at end of file diff --git a/docs/usage.ipynb b/docs/usage.ipynb new file mode 100644 index 00000000..4b381547 --- /dev/null +++ b/docs/usage.ipynb @@ -0,0 +1,608 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Usage\n", + "\n", + "## Getting Started" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Much of the functionality of PyMatching is available through the `pymatching.matching.Matching` class, which can be imported in Python with:" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from pymatching import Matching" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The simplest way to construct a `Matching` object is from the X or Z check matrix of the code, which can be given as a numpy or a scipy array. For example, we can construct the matching graph for a five-qubit quantum bit-flip repetition code (which has stabilisers $ZZIII$, $IZZII$, $IIZZI$ and $IIIZZ$) from the $Z$ check matrix using:" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import numpy as np\n", + "\n", + "Hz = np.array([\n", + " [1,1,0,0,0],\n", + " [0,1,1,0,0],\n", + " [0,0,1,1,0],\n", + " [0,0,0,1,1]\n", + "])\n", + "\n", + "m = Matching(Hz)\n", + "m" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that, since two qubits (0 and 4) are incident to only a single stabiliser, a boundary node has automatically been created in the matching graph, and is connected to the stabilisers acting non-trivially on qubits 0 and 4.\n", + "\n", + "If $X$ errors occur on the third and fourth qubits we have a binary noise vector:" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "noise = np.array([0,0,1,1,0])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "and a the resulting syndrome vector is:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[0 1 0 1]\n" + ] + } + ], + "source": [ + "z = Hz@noise % 2\n", + "print(z)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This syndrome vector `z` can then be decoded simply using:" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "c: [0 0 1 1 0], of type \n" + ] + } + ], + "source": [ + "c = m.decode(z)\n", + "print(f\"c: {c}, of type {type(c)}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where `c` is the $X$ correction operator (i.e. $IIXXI$)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that for larger check matrices you may instead prefer to use a scipy sparse matrix to represent the check matrix:" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import scipy\n", + "\n", + "Hz = scipy.sparse.csr_matrix(Hz)\n", + "m = Matching(Hz)\n", + "m" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Noisy Syndromes" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Spacetime matching graph\n", + "\n", + "If stabiliser measurements are instead noisy, then each stabiliser measurement must be repeated, with each defect in the matching graph corresponding to a change in the syndrome (see IV B of [this paper](https://arxiv.org/abs/quant-ph/0110143)). We will repeat each stabiliser measurement 5 times, with each qubit suffering an $X$ error with probability `p`, and each stabiliser will be measured incorrectly with probability `q`. Spacelike edges will be weighted with $\\log((1-p)/p)$ and timelike edges will be weighted with $\\log((1-q)/1)$. The Matching object representing this spacetime matching graph can be constructed using:" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "repetitions=5\n", + "p = 0.05\n", + "q = 0.05\n", + "m2d = Matching(Hz, \n", + " spacelike_weights=np.log((1-p)/p),\n", + " repetitions=repetitions,\n", + " timelike_weights=np.log((1-q)/q)\n", + " )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Simulate noisy syndromes\n", + "\n", + "Now if each qubit suffers an $X$ error with probability `p` in each round of stabiliser measurements, the errors on the data qubits can be given as a 2D numpy array:" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[0, 0, 1, 0, 0],\n", + " [0, 0, 0, 0, 0],\n", + " [0, 0, 0, 0, 1],\n", + " [0, 0, 0, 0, 0],\n", + " [0, 0, 0, 0, 0]])" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "num_stabilisers, num_qubits = Hz.shape\n", + "np.random.seed(1) # Keep RNG deterministic\n", + "noise_new = (np.random.rand(num_qubits, repetitions) < p).astype(int)\n", + "noise_new # New errors in each time step" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[0, 0, 1, 1, 1],\n", + " [0, 0, 0, 0, 0],\n", + " [0, 0, 0, 0, 1],\n", + " [0, 0, 0, 0, 0],\n", + " [0, 0, 0, 0, 0]])" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "noise_cumulative = np.cumsum(noise_new, 1) % 2\n", + "noise_total = noise_cumulative[:,-1] # Total cumulative noise at the last round\n", + "noise_cumulative # Cumulative errors in each time step" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The corresponding noiseless syndrome would be:" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[0, 0, 1, 1, 1],\n", + " [0, 0, 0, 0, 1],\n", + " [0, 0, 0, 0, 1],\n", + " [0, 0, 0, 0, 0]])" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "z_noiseless = Hz@noise_cumulative % 2\n", + "z_noiseless # Noiseless syndrome" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We assume each syndrome measurement is incorrect with probability `q`, but that the perfect the last round of measurements is perfect to ensure an even number of defects (a simple approximation - the overlapping recovery method that would be used in practice):" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[0, 0, 1, 0, 0],\n", + " [0, 0, 0, 0, 0],\n", + " [0, 0, 0, 1, 0],\n", + " [0, 0, 0, 0, 0]])" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "z_err = (np.random.rand(num_stabilisers, repetitions) < q).astype(int)\n", + "z_err[:,-1] = 0\n", + "z_err # Syndrome errors" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[0, 0, 0, 1, 1],\n", + " [0, 0, 0, 0, 1],\n", + " [0, 0, 0, 1, 1],\n", + " [0, 0, 0, 0, 0]])" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "z_noisy = (z_noiseless + z_err) % 2\n", + "z_noisy # Noisy syndromes" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[0, 0, 0, 1, 0],\n", + " [0, 0, 0, 0, 1],\n", + " [0, 0, 0, 1, 0],\n", + " [0, 0, 0, 0, 0]])" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "z_noisy[:,1:] = (z_noisy[:,1:] - z_noisy[:,0:-1]) % 2\n", + "z_noisy # Convert to difference syndrome" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Decode\n", + "\n", + "Decoding can now be done just by inputting this 2D syndrome vector to the `Matching.decode` method:" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([1, 0, 1, 0, 0], dtype=uint8)" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "correction = m2d.decode(z_noisy)\n", + "correction" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And we see that this correction operator successfully corrects the cumulative total noise:" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([0, 0, 0, 0, 0])" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "(noise_total + correction) % 2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Loading from NetworkX graphs\n", + "\n", + "While it can be convenient to decode directly from the check matrices, especially when simulating under a standard independent or phenomenological noise model, it is sometimes necessary to construct the matching graph nodes, edges, weights and boundaries explicitly. This is useful for decoding under more complicated (e.g. circuit-level) noise models, for which matching graph edges can be between nodes separated in both space and time (\"diagonal edges\"). Furthermore, there can be so called \"hook errors\" which are single faults (matching graph edges) corresponding to errors on two or more qubits.\n", + "\n", + "To provide the functionality to handle these use cases, PyMatching allows Matching objects to be constructed explicitly from [NetworkX](https://networkx.org/documentation/stable/index.html) graphs.\n", + "\n", + "Each node in the matching graph with $n$ nodes, represented by the `pymatching.Matching` object, should be uniquely identified by an integer between $0$ and $n-1$ (inclusive). Edges are then added between these integer nodes, with optional attributes `weight`, `qubit_id` and `error_probability`. \n", + "\n", + "We will again use the quantum repetition code as an example. Let's create a quantum repetition code `Matching` object on five qubits from a NetworkX graph:" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "import networkx as nx\n", + "\n", + "p = 0.2\n", + "g = nx.Graph()\n", + "g.add_edge(0, 1, qubit_id=0, weight=np.log((1-p)/p), error_probability=p)\n", + "g.add_edge(1, 2, qubit_id=1, weight=np.log((1-p)/p), error_probability=p)\n", + "g.add_edge(2, 3, qubit_id=2, weight=np.log((1-p)/p), error_probability=p)\n", + "g.add_edge(3, 4, qubit_id=3, weight=np.log((1-p)/p), error_probability=p)\n", + "g.add_edge(4, 5, qubit_id=4, weight=np.log((1-p)/p), error_probability=p)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Just for the purpose of demonstration, we'll assume that there is also an error process that gives a single hook error on qubits $2$ and $3$, corresponding to a single edge between node $2$ and node $4$. This error will occur with probability `p2`. This can be added using:" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [], + "source": [ + "p2 = 0.12\n", + "g.add_edge(2, 4, qubit_id={2, 3}, weight=np.log((1-p2)/p2), error_probability=p2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Since nodes 0 and 5 are incident to only a single edge, they must be specified as being boundary nodes, which can be done by setting their optional `is_boundary` attribute to `True`:" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [], + "source": [ + "g.node[0]['is_boundary'] = True\n", + "g.node[5]['is_boundary'] = True" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, we will connect these boundary nodes with an edge of `weight` zero, and with a `qubit_id` either unspecified or set to `set()` or `-1` (since edges between boundaries do not correspond to qubit faults):" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [], + "source": [ + "g.add_edge(0, 5, weight=0.0, qubit_id=-1, error_probability=0.0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, we can now use this NetworkX graph to construct the `Matching` object:" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [], + "source": [ + "m = Matching(g)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "While the noise and syndrome can be generated separately without PyMatching, if the optional `error_probability` attribute is given to every edge, then the edges can be flipped independently using the `add_noise` method:" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[0 1 0 0 0]\n", + "[0 1 1 0 0 0]\n" + ] + } + ], + "source": [ + "from pymatching import set_seed\n", + "set_seed(1) # Keep RNG deterministic\n", + "\n", + "noise, syndrome = m.add_noise()\n", + "print(noise)\n", + "print(syndrome)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can now decode as before using the `decode` method:" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[0 0 0 0 0]\n" + ] + } + ], + "source": [ + "correction = m.decode(syndrome)\n", + "print((correction+noise)%2)" + ] + } + ], + "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.7.4" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/src/pymatching/matching.py b/src/pymatching/matching.py index 33809217..6b9ad12f 100644 --- a/src/pymatching/matching.py +++ b/src/pymatching/matching.py @@ -31,27 +31,44 @@ def _find_boundary_nodes(G): class Matching: """A class for constructing matching graphs and decoding using the minimum-weight perfect matching decoder - [extended_summary] + The Matching class provides most of the core functionality of PyMatching. + A PyMatching object can be constructed from the :math:`Z` or + :math:`X` check matrix of the quantum code, given as a `scipy.sparse` + matrix or `numpy.ndarray`, along with additional argument specifying the + edge weights, error probabilities and number of repetitions. + Alternatively, a Matching object can be constructed from a NetworkX + graph, with node and edge attributes used to specify edge weights, + qubit ids, boundaries and error probabilities. """ - def __init__(self, H, weights=None, + def __init__(self, H, spacelike_weights=None, error_probabilities=None, repetitions=None, - timelike_weight=None, + timelike_weights=None, measurement_error_probability=None, precompute_shortest_paths=False): """Constructor for the Matching class - [extended_summary] - Parameters ---------- H : `scipy.spmatrix` or `numpy.ndarray` or `networkx.Graph` object The quantum code to be decoded with minimum-weight perfect matching, given either as a binary check matrix (scipy sparse - matrix or numpy.ndarray), or as a matching grap (NetworkX graph). - weights : float or numpy.ndarray, optional - If `H` is given as a scipy or numpy array, `weights` gives the weights - of edges in the matching graph. by default None + matrix or numpy.ndarray), or as a matching graph (NetworkX graph). + If `H` is given as a NetworkX graph with `M` nodes, each node + `m` in `H` should be an integer :math:`01`, - `timelike_weight` gives the weight of timelike edges, by default None + `timelike_weights` gives the weight of timelike edges. By default + None, in which case all weights are set to 1.0 measurement_error_probability : float, optional If `H` is given as a scipy or numpy array and `repetitions>1`, gives the probability of a measurement error to be used for @@ -94,7 +112,7 @@ def __init__(self, H, weights=None, else: num_edges = nx.number_of_edges(H) - weights = 1.0 if weights is None else weights + weights = 1.0 if spacelike_weights is None else spacelike_weights if isinstance(weights, (int, float)): weights = np.array([weights]*num_edges).astype(float) weights = np.asarray(weights) @@ -119,7 +137,7 @@ def __init__(self, H, weights=None, if np.any(weights < 0.): raise ValueError("All weights must be non-negative.") - timelike_weight = 1.0 if timelike_weight is None else timelike_weight + timelike_weights = 1.0 if timelike_weights is None else timelike_weights repetitions = 1 if repetitions is None else repetitions p_meas = measurement_error_probability if measurement_error_probability is not None else -1 boundary = [H.shape[0]*repetitions] if 1 in unique_column_weights else [] @@ -134,7 +152,7 @@ def __init__(self, H, weights=None, for t in range(repetitions-1): for i in range(H.shape[0]): self.stabiliser_graph.add_edge(i+t*H.shape[0], i+(t+1)*H.shape[0], - set(), timelike_weight, p_meas, p_meas >= 0) + set(), timelike_weights, p_meas, p_meas >= 0) else: boundary = _find_boundary_nodes(H) num_nodes = H.number_of_nodes() @@ -179,8 +197,6 @@ def boundary(self): def decode(self, z, num_neighbours=20): """Decode the syndrome `z` using minimum-weight perfect matching - [extended_summary] - Parameters ---------- z : numpy.ndarray @@ -199,8 +215,8 @@ def decode(self, z, num_neighbours=20): than 10. If `num_neighbours=None`, then instead full matching is performed, with the all-pairs shortest paths precomputed and cached the first time it is used. Since full matching is more - memory intensive, it is only recommended for matching graphs - with less than around 10,000 nodes, and is only faster than + memory intensive, it is not recommended to be used for matching graphs + with more than around 10,000 nodes, and is only faster than local matching for matching graphs with less than around 1,000 nodes. By default 20 @@ -242,9 +258,7 @@ def decode(self, z, num_neighbours=20): return decode_match_neighbourhood(self.stabiliser_graph, defects, num_neighbours) def add_noise(self): - """Add noise by flipping edges in the stabiliser graph - - Add noise by flipping edges in the stabiliser graph with + """Add noise by flipping edges in the stabiliser graph with a probability given by the error_probility edge attribute. This is currently only supported for weighted matching graphs initialised using a NetworkX graph. @@ -261,3 +275,6 @@ def add_noise(self): if not self.stabiliser_graph.all_edges_have_error_probabilities: return None return self.stabiliser_graph.add_noise() + + def __repr__(self): + return f"" diff --git a/src/pymatching/mwpm.cpp b/src/pymatching/mwpm.cpp index 6adae6ce..1e0423e5 100644 --- a/src/pymatching/mwpm.cpp +++ b/src/pymatching/mwpm.cpp @@ -18,12 +18,6 @@ py::array_t Decode(IStabiliserGraph& sg, const py::array_t& d int num_nodes = d.shape(0); int num_edges = num_nodes*(num_nodes-1)/2; - // std::clock_t start, before_solve, after_solve, end; - // double create_matching, solve_matching, make_correction; - - // start = std::clock(); - - PerfectMatching *pm = new PerfectMatching(num_nodes, num_edges); for (py::size_t i = 0; i Decode(IStabiliserGraph& sg, const py::array_t& d } }; pm->options.verbose = false; - // before_solve = std::clock(); pm->Solve(); - // after_solve = std::clock(); int N = sg.GetNumQubits(); auto correction = new std::vector(N, 0); std::set qids; @@ -52,13 +44,6 @@ py::array_t Decode(IStabiliserGraph& sg, const py::array_t& d } } delete pm; - // end = std::clock(); - // create_matching = ( before_solve - start ) / (double) CLOCKS_PER_SEC; - // solve_matching = ( after_solve - before_solve ) / (double) CLOCKS_PER_SEC; - // make_correction = ( end - after_solve ) / (double) CLOCKS_PER_SEC; - // std::cout<<"Create matching: "<*>(correction); }); return py::array_t(correction->size(), correction->data(), capsule);