Quantum-Error-Decoding-using-NEAT.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0850736d",
   "metadata": {},
   "source": [
    "# Quantum Error Decoding using NEAT"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa395a60",
   "metadata": {},
   "source": [
    "<div style=\"text-align: justify\">This is my account of the QOSF Mentorship project that I realized in the second quarter of 2021 under supervision of Dr. Dario Rosa. The goal was to implement the following paper titled \"A NEAT Quantum Decoder\" by Hugo Théveniaut and Evert van Nieuwenburg [1]: https://arxiv.org/abs/2101.08093 and learn about toric code and neuro-evolution along the way. In order to understand them in depth, I didn't use any additional libraries except for numpy. The paper itself presents a new method to conduct active protection of the toric code against random bit- and phase-flips, based on Neuro-Evolutionary Augmented Topologies algorithm. It performs similarly as good as already known methods, but uses a few orders of magnitude less parameters and thus has a significant advantage over them. This notebook is intended to be a gentle introduction to both toric code and NEAT as well as contain my implementation of the paper. However, the last part is still work in progress as I didn't manage to finish the implementation on time and it is going to be added in the future, so for now the aim of the project was fulfilled only partially.</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aacb5dfe",
   "metadata": {},
   "source": [
    "### Part I. Toric Code"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07db553a",
   "metadata": {},
   "source": [
    "In the era of so-called NISQ devices (Noisy Intermediate-Scale Quantum) researchers have to deal with the noise constantly introducing logical errors to qubits. In classical codes where we operate with zeroes and ones only one type of error can occur, namely a bit-flip which is unintentional swap of $0$ into $1$ or the other way around. But in the quantum realm we can in fact distinguish two types of errors."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "855c2ab6",
   "metadata": {},
   "source": [
    "First one is a <i>bit-flip</i> which transforms an arbitrary quantum state $\\alpha |0> + \\beta |1>$ to $\\beta |0> + \\alpha |1>$, so basically it's equivalent to applying $X$ gate to a qubit. And the second one is a <i>phase-flip</i> which transforms an arbitrary quantum state $\\alpha |0> + \\beta |1>$ to $\\alpha |0> - \\beta |1>$, so basically it's equivalent to applying $Z$ gate to a qubit, sometimes it's called <i>depolarization</i>."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a840b09a",
   "metadata": {},
   "source": [
    "To prevent that it is possible to encode an information stored by one qubit in many qubits - most often the more the better. There are for example three-qubit encodings of one logical qubit allowing to prevent a single bit-flip and analogously similar encoding for preventing phase-flips. There is also a nine-qubit encoding of one logical qubit called Shor code covering all possible errors, namely bit-flips, phase-flips and their linear combinations as well."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea93493a",
   "metadata": {},
   "source": [
    " However, not so long after Shor in 1997 A. Kitaev proposed a new way of such encoding in [2]. Two logical qubits are realized as $2N^2$ physical qubits arranged in a $N \\times N$ lattice whose opposite edges are identified (i.e. glued together), hence it creates a torus on which qubits lie. On the picture from [1] I showed which qubits exactly are identified using colored hexagons. On this picture each corner of the chessboard is a physical qubit whose state of course can be changed by some random Pauli operator. Moreover, the toric code is defined as a stabilizer code, i.e. there are operators called stabilizers which don't change the logical state of the lattice by definition. These are plaquette and star which can be applied to squares of the <i>chessboard</i>. They are respectively a product of $Z$ operators on corners of a black square and a product of $X$ operators on corners of a white square. Again, by definition as long as we can return to the empty lattice (without syndromes) with these operators solely, its logical state hasn't change."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1959c823",
   "metadata": {},
   "source": [
    "<img src=\"images/img11.png\" width=600>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "729d8550",
   "metadata": {},
   "source": [
    "Since they don't change the logical value of the system they can be measured on and on and give us constantly some information on what's going on. The measurement circuit for both types is shown below. And a visualization of a toric code, when the opposite edges are glued together. The following two pictures come from [4]."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b8f9e5c",
   "metadata": {},
   "source": [
    "<img src=\"images/img13.png\" width=600>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9cb5ac6",
   "metadata": {},
   "source": [
    "<img src=\"images/img20.png\" width=600>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c3d043e",
   "metadata": {},
   "source": [
    "As it follows from simple calculations, plaquette is measured $1$ iff there is an even amount of $X$ operators applied to the square we measure and $-1$ otherwise and vice versa - star is measured $1$ iff there is an even amount of $Z$ operators applied to the square we measure and $-1$ otherwise. What is worth repeating here, it is not possible to measure plaquette or star on each square of the <i>chessboard</i> without disturbing the logical state of the lattice. We can only measure plaquettes on black rows and stars on white rows. If stabilizer for a given square is measured $-1$ we say there's a syndrome and they are depicted on the first picture as orange circles."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75a459d2",
   "metadata": {},
   "source": [
    "So, the only knowledge we can get about the system while measuring those stars and plaquettes is where the syndromes currently are. It is not the full knowledge about the system but it will be sufficient to monitor the possible logical errors and correct its state by applying $X$, $Y$, $Z$ to physical qubits before a logical error occurs. Let's observe that one physical error, i.e. a single qubit flipped, generates two syndromes, but if two neighbouring qubits get the same error, then there are again two syndromes and not four for example. We can also add in a third error and so on and there will be still just two syndromes on the ends of what we call an \"error string\". Error strings are important because if they get closed creating a loop as these two endpoints meet each other, all syndromes disappear. That's because on each square's corners there is an even number of errors. If the loop goes even number of times through each of vertical and horizontal boundaries, it's called \"trivial\" and everything is OK, because we can eliminate this loop with stabilizer operators. However, if such a loop goes through boundary as you can see on the right half of the first picture, it is called \"non-trivial\" and can't be eliminated with plaquettes and stars solely. Thus, it introduces a logical error. Not going into details, the key topological concept here is contraction and it's possible to contract a loop (informally squeeze it to a point) only if it doesn't go around one of torus equators. On the other hand it goes through them exactly when it passes through the glued boundaries as you can check imagining the way a torus is assembled out of a square. Lastly, if there's a vertical non-trivial loop made of bit-flips (phase-flips), the logical state of the first logical qubit is bit-flipped (phase-flipped) and if there's a horizontal non-trivial loop, the state of the second qubit is flipped."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e76c225b",
   "metadata": {},
   "source": [
    "<img src=\"images/img30.png\" width=600>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adf1de82",
   "metadata": {},
   "source": [
    "Now for the coding part. Let's import the necessary libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "cafcc3c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math as m\n",
    "from random import random, choice\n",
    "from numpy.random import normal"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ee19ee0",
   "metadata": {},
   "source": [
    "First we need to define a lattice of physical qubits, which is going to be just a $2N \\times N$ matrix consisting of relevant operators applied to different qubits at any moment. These operators can be $X$, $Z$ or $Y$ which is a combination of two former ones. They will be introduced both at random what is intended to simulate the noisy enviroment and by our choice what means applying one of the three aforementioned gates to any of these physical qubits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a370ef42",
   "metadata": {},
   "outputs": [],
   "source": [
    "def lattice(N):\n",
    "    # N is the code distance and the lattice is 2N x N\n",
    "    array = []\n",
    "    for i in range(2*N):\n",
    "        array.append([])\n",
    "    for i in range(2*N**2):\n",
    "        array[i % (2*N)].append(\"╳\")\n",
    "    return array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f994512b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def error(lat, i, j, e):\n",
    "    # i, j are respectively row and column index of a qubit in the lattice lat\n",
    "    # e can be one of \"X\", \"Y\", \"Z\" meaning three Pauli operators\n",
    "    N = len(lat[0])\n",
    "    if (lat[i%(2*N)][j%N] == e): lat[i%(2*N)][j%N] = \"╳\"\n",
    "    elif (lat[i%(2*N)][j%N] == \"Y\"): \n",
    "        if (e == \"X\"): lat[i%(2*N)][j%N] = \"Z\"\n",
    "        else: lat[i%(2*N)][j%N] = \"X\" \n",
    "    elif (lat[i%(2*N)][j%N] != \"╳\"): lat[i%(2*N)][j%N] = \"Y\" \n",
    "    else: lat[i%(2*N)][j%N] = e"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5cbf87d",
   "metadata": {},
   "source": [
    "Another fundamental thing we need to code is the syndrome measurement which will be stored in another matrix. The important thing here is that plaquettes and stars can be measured in different rows of our lattice, so in even rows of the resulting matrix we will keep plaquette measurements and in the odd rows we will keep star measurements. If the measurement is $1$ (there's no syndrome) we write $0$ in our matrix and if it's $-1$ (there's a syndrome) we write $1$, the only purpose of such a conversion is our later convenience."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "18e77f9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def X(lat):\n",
    "    # an auxilliary matrix storing information about bit flips on the lattice\n",
    "    # it keeps -1 in (i, j) iff (i, j) qubit was bit-flipped, otherwise 1\n",
    "    N = len(lat[0])\n",
    "    array = lattice(N)\n",
    "    for i in range(2*N):\n",
    "        for j in range(N):\n",
    "            if (lat[i][j] == \"X\" or lat[i][j] == \"Y\"): array[i][j] = -1\n",
    "            else: array[i][j] = 1\n",
    "    return array\n",
    "\n",
    "def Z(lat):\n",
    "    # an auxilliary matrix storing information about phase flips on the lattice\n",
    "    # it keeps -1 in (i, j) iff (i, j) qubit was phase-flipped, otherwise 1\n",
    "    N = len(lat[0])\n",
    "    array = lattice(N)\n",
    "    for i in range(2*N):\n",
    "        for j in range(N):\n",
    "            if (lat[i][j] == \"Z\" or lat[i][j] == \"Y\"): array[i][j] = -1\n",
    "            else: array[i][j] = 1\n",
    "    return array\n",
    "\n",
    "def syndromes(lat):\n",
    "    N = len(lat[0])\n",
    "    X_lat = X(lat)\n",
    "    Z_lat = Z(lat)\n",
    "    array = lattice(N)\n",
    "    for i in range(2*N):\n",
    "        for j in range(N):\n",
    "            if (i%2 == 0):\n",
    "                array[i][j] = (1-X_lat[i][j]*X_lat[(i+1)%(2*N)][j]*X_lat[(i+2)%(2*N)][j]*X_lat[(i+1)%(2*N)][(j+1)%N])//2\n",
    "            else:\n",
    "                array[i][j] = (1-Z_lat[i][j]*Z_lat[(i+1)%(2*N)][j]*Z_lat[(i+2)%(2*N)][j]*Z_lat[(i+1)%(2*N)][(j-1)%N])//2\n",
    "    return array"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "815108ef",
   "metadata": {},
   "source": [
    "Below is a function that will help us visualise the current state of the lattice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "8eacad11",
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw(lat):\n",
    "    N = len(lat[0])\n",
    "    syn = syndromes(lat)\n",
    "    for i in range(2*N):\n",
    "        line = \"\"\n",
    "        if (i%2 == 0):\n",
    "            if (syn[(i-1)%(2*N)][0] == 1):\n",
    "                line += \"o \"\n",
    "            else:\n",
    "                line += \"  \"\n",
    "            for j in range(N):\n",
    "                line += lat[i][j]\n",
    "                if (syn[(i-1)%(2*N)][(j+1)%N] == 1):\n",
    "                    line += \" o \"\n",
    "                else:\n",
    "                    line += \"   \"\n",
    "            print(line)\n",
    "            line = \"\"\n",
    "            for j in range(N):\n",
    "                line += \" / \\\\\"\n",
    "            print(line)\n",
    "            line = \"\"\n",
    "        else:\n",
    "            for j in range(N):\n",
    "                line += lat[i][j]\n",
    "                if (syn[(i-1)%(2*N)][j%N] == 1):\n",
    "                    line += \" o \"\n",
    "                else:\n",
    "                    line += \"   \"\n",
    "            line += lat[i][0]\n",
    "            print(line)\n",
    "            line = \"\"\n",
    "            for j in range(N):\n",
    "                line += \" \\\\ /\"\n",
    "            print(line)\n",
    "            line = \"\"\n",
    "    if (syn[2*N-1][0] == 1):\n",
    "        line += \"o \"\n",
    "    else:\n",
    "        line += \"  \"\n",
    "    for j in range(N):\n",
    "        line += lat[0][j]\n",
    "        if (syn[2*N-1][(j+1)%N] == 1):\n",
    "            line += \" o \"\n",
    "        else:\n",
    "            line += \"   \"\n",
    "    print(line)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5b00471",
   "metadata": {},
   "source": [
    "For example let's generate a $5$x$10$ lattice and put some errors in there."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "654abb24",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  ╳ o ╳   ╳   X   ╳   \n",
      " / \\ / \\ / \\ / \\ / \\\n",
      "╳   Z   ╳   ╳ o ╳ o ╳\n",
      " \\ / \\ / \\ / \\ / \\ /\n",
      "o ╳ o ╳   ╳   ╳ o Y o \n",
      " / \\ / \\ / \\ / \\ / \\\n",
      "╳   ╳   ╳   ╳   ╳ o ╳\n",
      " \\ / \\ / \\ / \\ / \\ /\n",
      "  ╳   ╳   ╳   ╳   ╳   \n",
      " / \\ / \\ / \\ / \\ / \\\n",
      "╳   ╳   ╳   ╳   ╳   ╳\n",
      " \\ / \\ / \\ / \\ / \\ /\n",
      "  ╳   ╳   ╳   ╳   ╳   \n",
      " / \\ / \\ / \\ / \\ / \\\n",
      "╳   ╳   ╳   ╳   ╳   ╳\n",
      " \\ / \\ / \\ / \\ / \\ /\n",
      "  ╳   ╳   ╳   ╳   ╳   \n",
      " / \\ / \\ / \\ / \\ / \\\n",
      "╳   ╳   ╳   ╳ o ╳   ╳\n",
      " \\ / \\ / \\ / \\ / \\ /\n",
      "  ╳ o ╳   ╳   X   ╳   \n"
     ]
    }
   ],
   "source": [
    "qubits = lattice(5)\n",
    "error(qubits, 0, 3, \"X\")\n",
    "error(qubits, 1, 1, \"Z\")\n",
    "error(qubits, 2, 4, \"Y\")\n",
    "draw(qubits)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f336a9a",
   "metadata": {},
   "source": [
    "As you see syndromes appear inside the squares which have an odd number of qubits bit-flipped or phase-flipped and the first type of syndromes is found in even rows only while the other one in odd rows only. The next thing we want to implement is random bit-flip and depolarizing noise, which correspond to random $X$ and $Z$ operations on physical qubits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "cd21006a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def randombitflips(lat, prob):\n",
    "    # prob is the probability for each physical qubit to be flipped\n",
    "    N = len(lat[0])\n",
    "    for i in range(2*N):\n",
    "        for j in range(N):\n",
    "            r = np.random.rand()\n",
    "            if (r < prob): error(lat, i, j, \"X\")\n",
    "\n",
    "def randomdepolarizing(lat, prob):\n",
    "    # prob is the probability for each physical qubit to be depolarized\n",
    "    N = len(lat[0])\n",
    "    for i in range(2*N):\n",
    "        for j in range(N):\n",
    "            r = np.random.rand()\n",
    "            if (r < prob/3): error(lat, i, j, \"X\")\n",
    "            elif (r < 2*prob/3): error(lat, i, j, \"Z\")\n",
    "            elif (r < prob): error(lat, i, j, \"Y\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e513ccff",
   "metadata": {},
   "source": [
    "Now let's check how does it work:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "89fec8b6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  X   ╳   X   X   X   X   X   X   X   ╳   \n",
      " / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\\n",
      "╳   X o ╳ o ╳   ╳   ╳ o ╳ o X   X   ╳ o ╳\n",
      " \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ /\n",
      "  ╳   ╳   ╳   X   X   ╳   X   X   ╳   X   \n",
      " / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\\n",
      "╳   X o ╳   X o ╳   X   X   ╳ o X o ╳   ╳\n",
      " \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ /\n",
      "  X   ╳   X   X   ╳   ╳   ╳   X   ╳   X   \n",
      " / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\\n",
      "X o ╳ o ╳ o ╳   ╳ o ╳ o X   X   X o ╳   X\n",
      " \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ /\n",
      "  X   X   ╳   X   X   ╳   ╳   X   ╳   ╳   \n",
      " / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\\n",
      "╳   X   ╳   X   ╳   ╳   ╳ o ╳   X o ╳ o ╳\n",
      " \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ /\n",
      "  ╳   ╳   X   ╳   X   ╳   X   ╳   ╳   X   \n",
      " / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\\n",
      "X o ╳ o ╳   ╳   X o X   ╳ o ╳ o ╳   ╳ o X\n",
      " \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ /\n",
      "  ╳   X   X   X   ╳   X   ╳   X   ╳   X   \n",
      " / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\\n",
      "╳ o ╳ o ╳   ╳ o X   ╳ o ╳   ╳ o ╳ o ╳   ╳\n",
      " \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ /\n",
      "  X   ╳   X   X   X   ╳   ╳   ╳   X   X   \n",
      " / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\\n",
      "X   ╳   X   X   X o ╳ o ╳   X o X o ╳   X\n",
      " \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ /\n",
      "  ╳   X   X   X   X   X   X   X   X   ╳   \n",
      " / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\\n",
      "╳   X o ╳ o X   X   X o X o ╳ o ╳ o X   ╳\n",
      " \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ /\n",
      "  X   X   X   X   X   ╳   X   ╳   X   X   \n",
      " / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\\n",
      "X o X   X o X o X   ╳   ╳ o ╳   X   X o X\n",
      " \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ /\n",
      "  ╳   X   ╳   ╳   ╳   ╳   ╳   X   X   ╳   \n",
      " / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\\n",
      "╳   X   ╳   X o X   ╳   X   ╳ o X   X o ╳\n",
      " \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ /\n",
      "  X   ╳   X   X   X   X   X   X   X   ╳   \n"
     ]
    }
   ],
   "source": [
    "qubits10 = lattice(10) \n",
    "randombitflips(qubits10, 0.5)\n",
    "draw(qubits10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "cd9d200a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  Z   ╳   ╳   ╳   ╳   X o Y   \n",
      " / \\ / \\ / \\ / \\ / \\ / \\ / \\\n",
      "╳   Z   ╳   ╳   ╳ o ╳ o ╳ o ╳\n",
      " \\ / \\ / \\ / \\ / \\ / \\ / \\ /\n",
      "o ╳ o ╳   ╳   ╳   Y   ╳ o Z o \n",
      " / \\ / \\ / \\ / \\ / \\ / \\ / \\\n",
      "╳   ╳   ╳ o ╳   Z o Z o ╳ o ╳\n",
      " \\ / \\ / \\ / \\ / \\ / \\ / \\ /\n",
      "o ╳ o Z   Y o ╳   ╳ o X o Y o \n",
      " / \\ / \\ / \\ / \\ / \\ / \\ / \\\n",
      "╳   ╳ o X   ╳ o Y o ╳ o ╳ o ╳\n",
      " \\ / \\ / \\ / \\ / \\ / \\ / \\ /\n",
      "  ╳   ╳   ╳   ╳ o ╳   ╳   ╳   \n",
      " / \\ / \\ / \\ / \\ / \\ / \\ / \\\n",
      "╳   ╳ o ╳   ╳ o ╳   ╳   ╳   ╳\n",
      " \\ / \\ / \\ / \\ / \\ / \\ / \\ /\n",
      "  ╳   X   ╳   X   ╳   ╳ o ╳   \n",
      " / \\ / \\ / \\ / \\ / \\ / \\ / \\\n",
      "╳   ╳ o ╳ o ╳ o ╳   ╳ o Z o ╳\n",
      " \\ / \\ / \\ / \\ / \\ / \\ / \\ /\n",
      "o ╳ o Z o X o Z   ╳ o Y o Y o \n",
      " / \\ / \\ / \\ / \\ / \\ / \\ / \\\n",
      "╳   ╳   ╳ o ╳   Z   ╳ o ╳ o ╳\n",
      " \\ / \\ / \\ / \\ / \\ / \\ / \\ /\n",
      "  ╳   ╳ o Z o ╳ o ╳   ╳   ╳   \n",
      " / \\ / \\ / \\ / \\ / \\ / \\ / \\\n",
      "╳   ╳ o X o ╳   ╳   ╳ o ╳ o ╳\n",
      " \\ / \\ / \\ / \\ / \\ / \\ / \\ /\n",
      "  Z   ╳   ╳   ╳   ╳   X o Y   \n"
     ]
    }
   ],
   "source": [
    "qubits7 = lattice(7)\n",
    "randomdepolarizing(qubits7, 0.25)\n",
    "draw(qubits7)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49c7dd4d",
   "metadata": {},
   "source": [
    "Couple more functions below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "f8a7588e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def hamiltonian(lat):\n",
    "    # computes the hamiltonian of the system\n",
    "    N = len(lat[0])\n",
    "    syn = syndromes(lat)\n",
    "    amount = 0\n",
    "    for i in range(2*N):\n",
    "        for j in range(N):\n",
    "            if (syn[i][j] != 1): amount -= 1\n",
    "    return amount\n",
    "\n",
    "def plaquette(lat, i, j):\n",
    "    # (i,j) points to the square in i-th row and j-th column\n",
    "    N = len(lat[0])\n",
    "    if (i%2 == 0):\n",
    "        error(lat, i, j, \"Z\")\n",
    "        error(lat, (i+1)%(2*N), j, \"Z\")\n",
    "        error(lat, (i+2)%(2*N), j, \"Z\")\n",
    "        error(lat, (i+1)%(2*N), (j+1)%N, \"Z\")\n",
    "    else: print(\"You cannot place a plaquette there.\")\n",
    "    \n",
    "def star(lat, i, j):\n",
    "    # (i,j) points to the square in i-th row and j-th column \n",
    "    N = len(lat[0])\n",
    "    if (i%2 == 1):\n",
    "        error(lat, i, j, \"X\")\n",
    "        error(lat, (i+1)%(2*N), j, \"X\")\n",
    "        error(lat, (i+2)%(2*N), j, \"X\")\n",
    "        error(lat, (i+1)%(2*N), (j-1)%N, \"X\")\n",
    "    else: print(\"You cannot place a star there.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "466e2d1d",
   "metadata": {},
   "source": [
    "The hamiltonian of such a lattice is here just energy of its state. Namely, an empty lattice has the energy of $-2N^2$ and gets bigger by $1$ with each new syndrome, so that a lattice with syndromes inside all its squares has the energy of $0$. We're not gonna use it, but it's a nice physical intepretation what quantum decoding really means - it's minimizing the energy of a given quantum state without introducing any logical errors. We can for example check the hamiltonian of the last drawn lattice and an empty equivalent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "b50527e0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-52\n",
      "-98\n"
     ]
    }
   ],
   "source": [
    "print(hamiltonian(qubits7))\n",
    "print(hamiltonian(lattice(7))) # should be -2*7^2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "193a1bde",
   "metadata": {},
   "source": [
    "### Part II. NeuroEvolution of Augmenting Topologies"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a019fa37",
   "metadata": {},
   "source": [
    "Standard methods to adjust a neural network to perform a specific task are done with some sort of gradient involved. However, in the case of evolutionary algorithms we don't need it. For a task to be trained on we need a game which a neural network can play and get a positive reward if it succeeds (win) or a negative one otherwise (lose). In our case it's going to be a toric code decoding game, described in Part III. Changes to neural networks in some considered population are done by random mutations and they are directed by positive or negative feedback from the game all neural networks play. These which win the most often are passed to the next generation and then combined with others winners and replicated. It is only possible due to the notion of a genome of a neural network. Moreover, NEAT allows to evolve and adjust not only weights and biases but also the topology of a network, i.e. its shape, nodes, connections. The original paper were NEAT algorithm was first introduced is [3]. All images from now on come from [1]."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4894858f",
   "metadata": {},
   "source": [
    "NEAT step-by-step goes like this: \n",
    "1. Initialize a population of trivial neural networks.\n",
    "2. Measure what percentage of games do these neural networks win.\n",
    "3. Mutate them randomly with a certain probability.\n",
    "4. Based on 2. leave the best individuals only.\n",
    "5. Divide them into species, which means briefly group them by similarity.\n",
    "6. Crossover individuals within species to get a new population.\n",
    "7. Repeat steps 2. to 6. how many times you want.\n",
    "8. Choose the best individual basing on its percentage of won games."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4bd6c6e",
   "metadata": {},
   "source": [
    "The crucial feature of NEAT is the genomic encoding of a neural network. In principle any neural network can be described in terms of a genome in which each gene corresponds to a node or connection and its properties."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b30585fa",
   "metadata": {},
   "source": [
    "<img src=\"images/img3.png\" width=500>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed5c67d6",
   "metadata": {},
   "source": [
    "Such an encoding allows us to define easily various types of mutations what I'll touch later. For example a node mutation can occur in which a random edge is divided on two edges and a new hidden neuron between them is created."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19e100f5",
   "metadata": {},
   "source": [
    "<img src=\"images/img4.png\" width=400>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b41cd122",
   "metadata": {},
   "source": [
    "The heart of the algorithm is crossover which given two genomes produces its offspring with parental features inherited randomly."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "586e77a8",
   "metadata": {},
   "source": [
    "<img src=\"images/img1.png\" width=500>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "613adb89",
   "metadata": {},
   "source": [
    "In order to implement NEAT algorithm we need to agree on what does a genome look like and how are we going to define it as part of the code. Well, let's say it's a 2-element list whose first element is a list of all node genes and analogously the second one is a list of all connection genes. Furthermore, a node gene is a list of the node properties as follows: $[\\text{id}, \\text{type}, \\text{bias}]$, where $\\text{id}$ is identification number, $\\text{type}$ is \"in\", \"out\" or \"hid\" which correspond to different types of nodes in a neural network and $\\text{bias}$ is just... well, a bias of that node, i.e. a number that is added to the weighted sum of signals from preceding neurons, but before the neuron activation is computed. Input nodes have always bias equal to 0. Whereas a connection gene of a connection between nodes $\\text{A}$ and $\\text{B}$ is a list of the connection properties as follows: $[\\text{innovation number}, \\text{id of node A}, \\text{id of node B}, \\text{weight}, \\text{status}]$. Innovation number is a unique number assigned to a connection gene whenever a new connection is created via mutation. This way we can easily track the history of the evolution history. Status can be $0$ or $1$ and means just \"turned off\" or \"turned on\" respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa02661e",
   "metadata": {},
   "source": [
    "To give an example, a very simple neural network consisting of two input neurons which are connected to a hid neuron, which is then connected to an output neuron, might have a following genome: "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67afa9c3",
   "metadata": {},
   "source": [
    "<center>$[[[0, \\text{\"in\"}, 0],[1, \\text{\"in\"}, 0],[2, \\text{\"hid\"}, 3.5],[3, \\text{\"out\"}, 0.2]],[[0, 0, 2, 0.1, 1],[1, 1, 2, 0.4, 1],[2, 2, 3, 0.7, 1]]]$</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78e37798",
   "metadata": {},
   "source": [
    "For understanding the further code it's important to remember that $\\text{genome}[0]$ refers here to the list of nodes genes and $\\text{genome}[1]$ refers here to the list of connection genes. Below is defined a function that for given genome and inputs it returns the probabilities for each of output neurons for a neural network described by the genome."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2fdf67f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigmoid(x):\n",
    "    # it will be used as an activation function in hidden and output neurons\n",
    "    return 1 / (1 + m.exp(-x))\n",
    "\n",
    "def compute(genome, nodes, k):\n",
    "    # a subroutine to compute a value of a k^th neuron if all its inputs are known\n",
    "    weighted_sum = 0\n",
    "    for y in genome[1]:\n",
    "        if (y[4] == 1 and y[2] == k):\n",
    "            if (nodes[y[1]] == \"empty\"):\n",
    "                return \"empty\"\n",
    "            weighted_sum += y[3] * nodes[y[1]]\n",
    "    result = sigmoid(weighted_sum + genome[0][k][2])\n",
    "    return result\n",
    "\n",
    "def evaluate(genome, inputs):\n",
    "    # it does the forward propagation given genome and inputs and returns outputs\n",
    "    nodes = []\n",
    "    \n",
    "    for x in genome[0]:\n",
    "            nodes.append(\"empty\")\n",
    "    \n",
    "    i = 0\n",
    "    for k in range(len(nodes)):\n",
    "        if (genome[0][k][1] == \"in\"):\n",
    "            nodes[k] = inputs[i]\n",
    "            i += 1\n",
    "    if (i != len(inputs)):\n",
    "        print(\"Number of input neurons is wrong.\")\n",
    "    \n",
    "    while (\"empty\" in nodes):\n",
    "        for k in range(len(nodes)):\n",
    "            if (nodes[k] == \"empty\"):\n",
    "                nodes[k] = compute(genome, nodes, k)\n",
    "    \n",
    "    outputs = []\n",
    "    for i in range(len(genome[0])):\n",
    "        if (genome[0][i][1] == \"out\"):\n",
    "            outputs.append(nodes[i])\n",
    "    return outputs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ce66a42",
   "metadata": {},
   "source": [
    "For later implementation we'll need a so-called trivial network that takes $9$ inputs and outputs probabilities of $4$ possible choices in a random way, i.e. biases and weights are generated by the normal distribution. The initial population will consist of trivial networks only."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6ea6f59d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def trivial_network():\n",
    "    # the initial values of weights and node biases are taken from gaussian distribution N(0,1)\n",
    "    first = [[0, \"in\"], [1, \"in\"], [2, \"in\"], [3, \"in\"], [4, \"in\"], [5, \"in\"], [6, \"in\"], [7, \"in\"], [8, \"in\"], [9, \"out\"], [10,\"out\"], [11, \"out\"], [12, \"out\"]]\n",
    "    second = []\n",
    "    for i in range(9):\n",
    "        first[i].append(0)\n",
    "        for j in range(4):\n",
    "            weight = normal()\n",
    "            second.append([len(second), i, 9+j, weight, 1])\n",
    "    for i in range(4):\n",
    "        bias = normal()\n",
    "        first[9+i].append(bias)\n",
    "    genome = [first, second]\n",
    "    return genome"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "038145ee",
   "metadata": {},
   "source": [
    "Next thing is a mutation procedure that takes a genome and modifies it in one of $8$ possible ways: either adds a new connection or removes an existing connection or adds a node or removes a node or changes a weight or changes a bias or enables a connection or disables a connection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9946c637",
   "metadata": {},
   "outputs": [],
   "source": [
    "def mutate(genome, max_node, max_conn):\n",
    "    # max_node is the largest node index present in the whole population\n",
    "    # max_conn is the largest innovation number in the whole population\n",
    "    r = random()\n",
    "    if (r < 0.1):\n",
    "        #add connection\n",
    "        pairs = []\n",
    "        for x1 in genome[0]:\n",
    "            for x2 in genome[0]:\n",
    "                if ((\"in\" in x1 and \"in\" in x2) or (\"out\" in x1 and \"out\" in x2)):\n",
    "                    continue\n",
    "                areConnected = 0\n",
    "                for y in genome[1]:\n",
    "                    if (y[1] == x1[0] and y[2] == x2[0]):\n",
    "                        areConnected = 1\n",
    "                if (areConnected == 0):\n",
    "                    pairs.append([x1[0], x2[0]])\n",
    "        if (len(pairs) != 0):\n",
    "            max_conn += 1\n",
    "            p = choice(pairs)\n",
    "            r = normal()\n",
    "            genome[1].append([max_conn, p[0], p[1], r, 1])\n",
    "        \n",
    "    elif (r < 0.2):\n",
    "        #remove connection\n",
    "        y = choice(genome[1])\n",
    "        genome[1].remove(y)\n",
    "        \n",
    "    elif (r < 0.3):\n",
    "        #add node\n",
    "        max_node += 1\n",
    "        max_conn += 2\n",
    "        y = choice(genome[1])\n",
    "        r = normal()\n",
    "        genome[0].append([max_node, \"hid\", r])\n",
    "        r = normal()\n",
    "        genome[1].append([max_conn-1, y[1], max_node, r, 1])\n",
    "        r = normal()\n",
    "        genome[1].append([max_conn, max_node, y[2], r, 1])\n",
    "        y[4] = 0\n",
    "        \n",
    "    elif (r < 0.4):\n",
    "        #remove node\n",
    "        hidden = []\n",
    "        for x in genome[0]:\n",
    "            if (\"hid\" in x):\n",
    "                hidden.append(x)\n",
    "        if (len(hidden) != 0):\n",
    "            x = choice(hidden)\n",
    "            connections = []\n",
    "            for y in genome[1]:\n",
    "                if (y[1] == x[0] or y[2] == x[0]):\n",
    "                    connections.append(y)\n",
    "            for y in connections:\n",
    "                genome[1].remove(y)\n",
    "            genome[0].remove(x)\n",
    "        \n",
    "    elif (r < 0.9):\n",
    "        #weight mutation\n",
    "        y = choice(genome[1])\n",
    "        r = normal()\n",
    "        y[3] = r\n",
    "        \n",
    "    else:\n",
    "        #bias mutation\n",
    "        x = [0, \"in\", 0]\n",
    "        while (x[1] == \"in\"):\n",
    "            x = choice(genome[0])\n",
    "        r = normal()\n",
    "        x[2] = r\n",
    "        \n",
    "    r = random()\n",
    "    if (r < 0.01):\n",
    "        #enable connection\n",
    "        disabled = []\n",
    "        for y in genome[1]:\n",
    "            if (y[4] == 0):\n",
    "                disabled.append(y)\n",
    "        if (len(disabled) != 0):\n",
    "            y = choice(disabled)\n",
    "            y[4] = 1\n",
    "    elif (r < 0.02):\n",
    "        #disable connection\n",
    "        enabled = []\n",
    "        for y in genome[1]:\n",
    "            if (y[4] == 1):\n",
    "                enabled.append(y)\n",
    "        if (len(enabled) != 0):\n",
    "            y = choice(enabled)\n",
    "            y[4] = 0\n",
    "            \n",
    "    result = [genome, max_node, max_conn]\n",
    "    \n",
    "    return result"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "927abc05",
   "metadata": {},
   "source": [
    "And the last on our list in this part is crossover. It's meant to imitate a biological crossover where, in order to create a new unique organism, genes of parents are combined together. Thus, if in both parental genomes there is a specific innovation number, a new genome inherits one of these genes with that number (other parameters might be still different) chosen randomly. The other connection genes, which are present only in one of parental genomes, are inherited randomly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "043ec647",
   "metadata": {},
   "outputs": [],
   "source": [
    "def crossover(genome1, genome2):\n",
    "    first = []\n",
    "    second = []\n",
    "    \n",
    "    for x in genome1[0]:\n",
    "        isyFound = 0\n",
    "        for y in genome2[0]:\n",
    "            if (x[0] == y[0]):\n",
    "                r = random()\n",
    "                if (r < 0.5):\n",
    "                    first.append(x)\n",
    "                else:\n",
    "                    first.append(y)\n",
    "                isyFound = 1\n",
    "        if (isyFound == 0):\n",
    "            first.append(x)\n",
    "        \n",
    "    for x in genome1[1]:\n",
    "        isyFound = 0\n",
    "        for y in genome2[1]:\n",
    "            if (x[0] == y[0]):\n",
    "                r = random()\n",
    "                if (r < 0.5):\n",
    "                    second.append(x)\n",
    "                else:\n",
    "                    second.append(y)\n",
    "                isyFound = 1\n",
    "        if (isyFound == 0):\n",
    "            r = random()\n",
    "            if (r < 0.5):\n",
    "                second.append(x)\n",
    "                \n",
    "    for x in first:\n",
    "        isConnected = 0\n",
    "        for y in second:\n",
    "            if (y[1] == x[0] or y[2] == x[0]):\n",
    "                isConnected = 1\n",
    "                break\n",
    "        if (isConnected == 0):\n",
    "            first.remove(x)\n",
    "                \n",
    "    genome = [first, second]\n",
    "    return genome"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bc030fe",
   "metadata": {},
   "source": [
    "Moreover, in order to protect new innovations that occured in the population, all neural networks are divided into species. A distance function $\\delta = c_1 \\frac{E}{N_\\text{genes}} + c_2 \\frac{D}{N_\\text{genes}} + c_3 + \\bar{W}$ between two neural networks is introduced and basically all representants which are distant from a given neural network for less than some agreed upon threshold $\\delta_c$ belong to the same species. In the formula $E$ stands for excess genes between two genomes and similarly $D$ for disjoint genes, whereas $\\bar{W}$ is the average weight difference in the shared genes. And $N_\\text{genes}$ means the number of genes in the largest genome. The aim is to allow newly-mutated properties to adjust in a less-agressive environment before they are released to the whole population. For instance if there's a new connection it allows it to get better weight to be actually useful and contributing to the whole evolution process, instead of being dismissed before its potential would be unleashed."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "052118df",
   "metadata": {},
   "source": [
    "### Part III. Implementation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a073b438",
   "metadata": {},
   "source": [
    "Now, when we already know the theory behind toric code and NEAT algorithm, we can eventually apply the latter to the former to obtain fault-tolerant encoding of qubits. In order to do so we need first to implement a subroutine checking is there any non-trivial loop in the lattice we work with. Next we can define the toric code decoding game, which can be summed up in two steps: 1. Apply a Pauli gate to one of qubits. 2. Check is there any non-trivial loop. If so, the game is lost. If not, go to the step 1. Repeat these steps until all syndromes disappear."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73b915b9",
   "metadata": {},
   "source": [
    "<img src=\"images/img10.png\" width=750>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6db1bc17",
   "metadata": {},
   "source": [
    "Above you can see how Correction #1 does not introduce a non-trivial loop, while Correction #2 does. Below is a subroutine defined that starts at each of upper squares of the lattice and explores it in the form of a ternary tree along all error strings which go through those squares."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b7835f0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def extend(input_list, index, value):\n",
    "    if (index > len(input_list)-1):\n",
    "        for i in range(index-len(input_list)+1):\n",
    "            input_list.append(0)\n",
    "    input_list[index] = value\n",
    "\n",
    "def newleaf(input_list, parent, value):\n",
    "    if(3*parent+4 > len(input_list)):\n",
    "        extend(input_list, 3*parent+3, 0)\n",
    "    for i in range(3):\n",
    "        if (input_list[3*parent+i+1] == value):\n",
    "            break\n",
    "        if (input_list[3*parent+i+1] == 0):\n",
    "            input_list[3*parent+i+1] = value\n",
    "            break\n",
    "        if (i == 2): print(\"This parent already possesses 3 leaves.\")\n",
    "\n",
    "def Xiteration(lat, tab, n):\n",
    "    # n is index of a vertex in the tree 'tab'\n",
    "    if (n > len(tab)-1 or tab[n] in [0, \"loop\", \"comeback\"]): return None\n",
    "    N = len(lat[0])\n",
    "    i = tab[n][0]               #coords of a given vertex\n",
    "    j = tab[n][1]\n",
    "    g = tab[m.floor((n-1)/3)][0]    #coords of a given vertex' parent\n",
    "    h = tab[m.floor((n-1)/3)][1]\n",
    "    if ((i != g or (j-1)%N != h) and (lat[(i+1)%(2*N)][j] in [\"X\", \"Y\"])):\n",
    "        newleaf(tab, n, [i, (j-1)%N])\n",
    "    if ((i != g or (j+1)%N != h) and (lat[(i+1)%(2*N)][(j+1)%N] in [\"X\", \"Y\"])):\n",
    "        newleaf(tab, n, [i, (j+1)%N])\n",
    "    if (((i+2)%(2*N) != g or j != h) and (lat[(i+2)%(2*N)][j] in [\"X\", \"Y\"])):\n",
    "        newleaf(tab, n, [(i+2)%(2*N), j%N])\n",
    "    if (((i-2)%(2*N) != g or j != h) and (lat[i][j] in [\"X\", \"Y\"])):\n",
    "        newleaf(tab, n, [(i-2)%(2*N), j])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc411881",
   "metadata": {},
   "source": [
    "Another function which if a square is met twice along one of branches of the tree walk tells us is it a trivial or non-trivial loop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "70abdeea",
   "metadata": {},
   "outputs": [],
   "source": [
    "def isLooped(lat, tab, n):\n",
    "    # n is index of a vertex in the tree 'tab'\n",
    "    N = len(lat[0])\n",
    "    i = n\n",
    "    counter = 0     #number of times the horizontal boundary was crossed\n",
    "    while (i > 0):\n",
    "        if (tab[i][0] - tab[m.floor((i-1)/3)][0] == 2*(N-1)):\n",
    "            counter -= 1\n",
    "        if (tab[i][0] - tab[m.floor((i-1)/3)][0] == 2*(-N+1)):\n",
    "            counter += 1\n",
    "        if (i > 0):\n",
    "            i = m.floor((i-1)/3)\n",
    "        if (tab[i] == tab[n]):\n",
    "            if (i == 0 and counter%2 == 1):\n",
    "                return \"comeback\"\n",
    "            else:\n",
    "                return \"loop\"\n",
    "    return tab[n]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdf08d27",
   "metadata": {},
   "source": [
    "And at last a function that returns $0$ iff there is a non-trivial loop going from up to down and after considering also the rotation of the lattice we get the function for both types of non-trivial loops - vertical ones as well as horizontal ones. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "7249d889",
   "metadata": {},
   "outputs": [],
   "source": [
    "def vertical_reward(lat):\n",
    "    N = len(lat[0])\n",
    "    k = 1\n",
    "    walk = []\n",
    "    maxwalk = 2\n",
    "    \n",
    "    for j in range(N):\n",
    "        walk.append([])\n",
    "        walk[j].append([0, j])\n",
    "        if(lat[1][j] in [\"X\",\"Y\"]):\n",
    "            walk[j].append([0,(j-1)%N])\n",
    "        if(lat[1][(j+1)%N] in [\"X\",\"Y\"]):\n",
    "            walk[j].append([0,(j+1)%N])\n",
    "        if(lat[2][j] in [\"X\",\"Y\"]):\n",
    "            walk[j].append([2,j])\n",
    "            \n",
    "    while (k < maxwalk and k < 3**(N**2+1)):\n",
    "        for i in range(N):\n",
    "            Xiteration(lat, walk[i], k)\n",
    "            if (3*k+3 < len(walk[i])):\n",
    "                for l in [3*k+1, 3*k+2, 3*k+3]:\n",
    "                    if (walk[i][l] != 0):\n",
    "                        walk[i][l] = isLooped(lat, walk[i], l)\n",
    "                        #print(k, i, walk[i][l])\n",
    "                        if (walk[i][l] == \"comeback\"): \n",
    "                            return 0\n",
    "            if (len(walk[i]) > maxwalk):\n",
    "                maxwalk = len(walk[i])\n",
    "        k += 1\n",
    "    #print(walk)\n",
    "    return 1\n",
    "\n",
    "def transpose(lat):\n",
    "    N = len(lat[0])\n",
    "    tab = []\n",
    "    for i in range(2*N):\n",
    "        tab.append([])\n",
    "    for j in range(N):\n",
    "        for i in range(N):\n",
    "            tab[2*j].append(lat[2*N-1-2*i][j])\n",
    "            tab[2*j+1].append(lat[2*N-1-2*i-1][j])\n",
    "    return tab\n",
    "\n",
    "def reward(lat):\n",
    "    return vertical_reward(lat) * vertical_reward(transpose(lat))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4cab5cd8",
   "metadata": {},
   "source": [
    "This reward can now be used for implementing the mentioned game as all choices that our networks will undertake can be quantatively measured with rewarding it for good choices."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6900a66c",
   "metadata": {},
   "source": [
    "<i> THE BELOW IMPLEMENTATION OF NEAT IS WORK IN PROGRESS </i>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca3c5284",
   "metadata": {},
   "outputs": [],
   "source": [
    "def empty(syn):\n",
    "    for row in syn:\n",
    "        if (1 in row):\n",
    "            return 0\n",
    "    return 1\n",
    "\n",
    "def vertex(i):\n",
    "    if (i == 0): return [0, 0]\n",
    "    if (i == 1): return [1, 0]\n",
    "    if (i == 2): return [1, 1]\n",
    "    if (i == 3): return [2, 0]\n",
    "    \n",
    "def best(probabilities):\n",
    "    # [probabilities for indices 0..4]\n",
    "    highest = [0, 0]\n",
    "    for i in range(len(probabilities)):\n",
    "        if (probabilities[i] > highest[1]):\n",
    "            highest = [i, probabilities[i]]\n",
    "    return highest \n",
    "\n",
    "def syndrome_neighborhood(syn, i, j):\n",
    "    N = len(syn[0])\n",
    "    synhood = []\n",
    "    synhood.append(syn[(i-2)%(2*N)][(j-1)%N])       #0\n",
    "    synhood.append(syn[(i-2)%(2*N)][j])             #1\n",
    "    synhood.append(syn[(i-2)%(2*N)][(j+1)%N])       #2\n",
    "    synhood.append(syn[i][(j-1)%N])                 #3\n",
    "    synhood.append(syn[i][j])                       #4\n",
    "    synhood.append(syn[i][(j+1)%N])                 #5\n",
    "    synhood.append(syn[(i+2)%(2*N)][(j-1)%N])       #6\n",
    "    synhood.append(syn[(i+2)%(2*N)][j])             #7\n",
    "    synhood.append(syn[(i+2)%(2*N)][(j+1)%N])       #8\n",
    "    return synhood\n",
    "\n",
    "def game(genome, N, prob, training):\n",
    "    # for training purpose use training == 1\n",
    "    # for performance evaluation use training == 0\n",
    "    \n",
    "    state = lattice(N)\n",
    "    randombitflips(state, prob)\n",
    "    syn = syndromes(state)\n",
    "    prob_max = 0\n",
    "    record = [] # for keeping a record of already taken actions\n",
    "    k = 0\n",
    "    \n",
    "    while(not empty(syn)):\n",
    "        \n",
    "        for i in range(N):\n",
    "            for j in range(N):\n",
    "                synhood = syndrome_neighborhood(syn, 2*i, j)\n",
    "                outcome = best(evaluate(genome, synhood))\n",
    "                # [index(0..4), probability]\n",
    "                if (outcome[1] > prob_max):\n",
    "                    prob_max = outcome[1]\n",
    "                    [i_max, j_max] = [2*i + vertex(outcome[0])[0], j + vertex(outcome[0])[1]]\n",
    "                        \n",
    "        if (training == 1 and [i_max, j_max] in record):\n",
    "            return 0\n",
    "        \n",
    "        if (training == 0 and k > 1000):\n",
    "            return 0\n",
    "        \n",
    "        error(state, i_max, j_max, \"X\")\n",
    "        record.append([i_max, j_max])\n",
    "        syn = syndromes(state)\n",
    "        k += 1\n",
    "        \n",
    "    return reward(state)\n",
    "\n",
    "def fitness(genome, N, prob, N_g):\n",
    "    # N_g is number of games which are played to determine fitness\n",
    "    points = 0\n",
    "    for i in range(N_g):\n",
    "        points += game(genome, N, prob, 0)\n",
    "        print(\"Game\", i+1, \"was played.\")\n",
    "    return points/N_g\n",
    "\n",
    "def distance(genome1, genome2):\n",
    "    c_1 = 1.0\n",
    "    c_2 = 1.0\n",
    "    c_3 = 0.4\n",
    "    # to be written\n",
    "    \n",
    "\n",
    "def NEAT(pop_size, num_generations, N, error_prob, N_g):\n",
    "    # hyperparameters from the original paper:\n",
    "    # pop_size from 100 to 300\n",
    "    # num_generations = ?\n",
    "    # number of games N_g = 400\n",
    "    # error_prob in 0.01, 0.05, 0.1, 0.15 (50 per each)\n",
    "    population = []\n",
    "    for i in range(pop_size):\n",
    "        population.append(trivial_network())\n",
    "    parents = []\n",
    "    max_node = len(population[0][0])-1 # the biggest neuron index in the population\n",
    "    max_conn = len(population[0][1])-1 # the biggest innovation number in the population\n",
    "    top = 5 # an arbitrary choice of how many top individuals we take for crossover\n",
    "    threshold_distance = 3.0\n",
    "    \n",
    "    for k in range(num_generations):\n",
    "        for genome in population:\n",
    "            performance = []\n",
    "            performance.append(fitness(genome, N, error_prob, N_g))\n",
    "            mutation = mutate(genome, max_node, max_conn)\n",
    "            max_node = mutation[1]\n",
    "            max_conn = mutation[2]\n",
    "        for j in range(top):\n",
    "            parents.append(population[best(performance)[0]])\n",
    "            performance[best(performance)[0]] = 0\n",
    "        population = []\n",
    "        for i in range(pop_size):\n",
    "            a = m.floor(random()*top)\n",
    "            b = m.floor(random()*top)\n",
    "            while(a == b):\n",
    "                b = m.floor(random()*top)\n",
    "            population.append(crossover(parents[a],parents[b]))\n",
    "            \n",
    "    for network in population:\n",
    "        performance = []\n",
    "        performance.append(fitness(network, N, error_prob, N_g))\n",
    "        \n",
    "    return population[best(performance)[0]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d171ea7",
   "metadata": {},
   "source": [
    "### References"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eba1792b",
   "metadata": {},
   "source": [
    "[1] H. Théveniaut, E. van Nieuwenburg <i>A NEAT Quantum Error Decoder</i>, January 2021, https://arxiv.org/abs/2101.08093"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b830ec07",
   "metadata": {},
   "source": [
    "[2] A. Kitaev <i>Fault-tolerant quantum computation by anyons</i>, July 1997, https://arxiv.org/abs/quant-ph/9707021"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18d4d539",
   "metadata": {},
   "source": [
    "[3] K. Stanley <i>Evolving Neural Networks through Augmenting Topologies</i>, June 2002, http://nn.cs.utexas.edu/downloads/papers/stanley.ec02.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69943048",
   "metadata": {},
   "source": [
    "[4] K. Fujii <i>Quantum Computation with Topological Codes</i>, April 2015, https://arxiv.org/pdf/1504.01444.pdf"
   ]
  }
 ],
 "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}