docs/source/state_tomography.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "yCpENKsOeRKt"
   },
   "source": [
    "# Quantum state tomography\n",
    "\n",
    "One can open this notebook in Google Colab (is recommended)\n",
    "\n",
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/LuchnikovI/QGOpt/blob/master/docs/source/state_tomography.ipynb)\n",
    "\n",
    "In this tutorial, we perform quantum state tomography via Riemannian optimization. First two blocks of a code (1. Many-qubit, informationally complete, positive operator-valued measure (IC POVM) and 2. Data set generation (measurement outcomes simulation)) are refered to data generation, third bock dedicated to tomography of a state.\n",
    "\n",
    "First, one needs to import all necessary libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "gdotZMMLeOn6"
   },
   "outputs": [],
   "source": [
    "import tensorflow as tf  # tf 2.x\n",
    "from math import sqrt\n",
    "\n",
    "try:\n",
    "    import QGOpt as qgo\n",
    "except ImportError:\n",
    "    !pip install git+https://github.com/LuchnikovI/QGOpt\n",
    "    import QGOpt as qgo\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from tqdm import tqdm\n",
    "\n",
    "# Fix random seed to make results reproducable.\n",
    "tf.random.set_seed(42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "azpIZCm2d8Ou"
   },
   "source": [
    "## 1. Many-qubit, informationally complete, positive operator-valued measure (IC POVM)\n",
    "\n",
    "Before generating measurement outcomes and performing quantum tomography, one needs to introduce POVM describing quantum measurements. For simplicity, we use one-qubit tetrahedral POVM and generalize it on a many-qubit case by taking tensor product between POVM elements, i.e. $\\{M_\\alpha\\}_{\\alpha=1}^4$ is the one-qubit tetrahedral POVM, $\\{M_{\\alpha_1}\\otimes \\dots \\otimes M_{\\alpha_N}\\}_{\\alpha_1=1,\\dots,\\alpha_N=1}^4$ is the many-qubits tetrahedral POVM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "OGPbEyFZd8Ov"
   },
   "outputs": [],
   "source": [
    "# Auxiliary function that returns Kronecker product between two\n",
    "# POVM elements A and B\n",
    "def kron(A, B):\n",
    "    \"\"\"Kronecker product of two POVM elements.\n",
    "\n",
    "    Args:\n",
    "        A: complex valued tensor of shape (q, n, k).\n",
    "        B: complex valued tensor of shape (p, m, l).\n",
    "\n",
    "    Returns:\n",
    "        complex valued tensor of shape (q * p, n * m, k * l)\"\"\"\n",
    "\n",
    "    AB = tf.tensordot(A, B, axes=0)\n",
    "    AB = tf.transpose(AB, (0, 3, 1, 4, 2, 5))\n",
    "    shape = AB.shape\n",
    "    AB = tf.reshape(AB, (shape[0] * shape[1],\n",
    "                         shape[2] * shape[3],\n",
    "                         shape[4] * shape[5]))\n",
    "    return AB\n",
    "\n",
    "# Pauli matrices\n",
    "sigma_x = tf.constant([[0, 1], [1, 0]], dtype=tf.complex128)\n",
    "sigma_y = tf.constant([[0 + 0j, -1j], [1j, 0 + 0j]], dtype=tf.complex128)\n",
    "sigma_z = tf.constant([[1, 0], [0, -1]], dtype=tf.complex128)\n",
    "\n",
    "# All Pauli matrices in one tensor of shape (3, 2, 2)\n",
    "sigma = tf.concat([sigma_x[tf.newaxis],\n",
    "                   sigma_y[tf.newaxis],\n",
    "                   sigma_z[tf.newaxis]], axis=0)\n",
    "\n",
    "# Coordinates of thetrahedron peaks (is needed to build tetrahedral POVM)\n",
    "s0 = tf.constant([0, 0, 1], dtype=tf.complex128)\n",
    "s1 = tf.constant([2 * sqrt(2) / 3, 0, -1/3], dtype=tf.complex128)\n",
    "s2 = tf.constant([-sqrt(2) / 3, sqrt(2 / 3), -1 / 3], dtype=tf.complex128)\n",
    "s3 = tf.constant([-sqrt(2) / 3, -sqrt(2 / 3), -1 / 3], dtype=tf.complex128)\n",
    "\n",
    "# Coordinates of thetrahedron peaks in one tensor of shape (4, 3)\n",
    "s = tf.concat([s0[tf.newaxis],\n",
    "               s1[tf.newaxis],\n",
    "               s2[tf.newaxis],\n",
    "               s3[tf.newaxis]], axis=0)\n",
    "\n",
    "# One qubit thetrahedral POVM\n",
    "M = 0.25 * (tf.eye(2, dtype=tf.complex128) + tf.tensordot(s, sigma, axes=1))\n",
    "\n",
    "n = 2  # number of qubits we experiment with\n",
    "\n",
    "# M for n qubits (Mmq)\n",
    "Mmq = M\n",
    "for _ in range(n - 1):\n",
    "    Mmq = kron(Mmq, M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "s9e5LH3rd8Oy"
   },
   "source": [
    "## 2. Data set generation (measurement outcomes simulation).\n",
    "\n",
    "Here we generate a set of measurement outcomes (training set). First of all, we generate a random density matrix that is a target state we want to reconstruct. Then, we simulate measurement outcomes over the target state driven by many-qubits tetrahedral POVM introduced in the previous cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "jA5UUexod8Oz"
   },
   "outputs": [],
   "source": [
    "#-----------------------------------------------------#\n",
    "num_of_meas = 600000  # number of measurement outcomes\n",
    "#-----------------------------------------------------#\n",
    "\n",
    "# random target density matrix (to be reconstructed)\n",
    "m = qgo.manifolds.DensityMatrix()\n",
    "A = m.random((2 ** n, 2 ** n), dtype=tf.complex128)\n",
    "rho_true = A @ tf.linalg.adjoint(A)\n",
    "\n",
    "# measurements simulation (by using Gumbel trick for sampling from a\n",
    "# discrete distribution)\n",
    "P = tf.cast(tf.tensordot(Mmq, rho_true, [[1, 2], [1, 0]]), dtype=tf.float64)\n",
    "eps = tf.random.uniform((num_of_meas, 2 ** (2 * n)), dtype=tf.float64)\n",
    "eps = -tf.math.log(-tf.math.log(eps))\n",
    "ind_set = tf.math.argmax(eps + tf.math.log(P), axis=-1)\n",
    "\n",
    "# POVM elements came true (data set)\n",
    "data_set = tf.gather_nd(Mmq, ind_set[:, tf.newaxis])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Pvbfxbpqd8O3"
   },
   "source": [
    "## 3. Data processing (tomography)\n",
    "\n",
    "First, we define an example of the density matrices manifold:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = qgo.manifolds.DensityMatrix()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The manifold of density matrices is represneted through the quadratic parametrization $\\varrho = AA^\\dagger$ with an equivalence relation $A\\sim AQ$, where $Q$ is an arbitrary unitary matrix. Thus, we initialize a variable, that represents the parametrization of a density matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# random initial paramterization\n",
    "a = m.random((2 ** n, 2 ** n), dtype=tf.complex128)\n",
    "# in order to make an optimizer works properly\n",
    "# one need to turn a to real representation\n",
    "a = qgo.manifolds.complex_to_real(a)\n",
    "# variable\n",
    "a = tf.Variable(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we initialize Riemannian Adam optimizer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr = 0.07  # optimization step size\n",
    "opt = qgo.optimizers.RAdam(m, lr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we ran part of code that calculate forward pass, gradients, and optimization step several times until convergence is reached:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 68
    },
    "colab_type": "code",
    "id": "6MlYR2h7d8O3",
    "outputId": "4ad95b33-a0a4-461b-830b-4740af82cdcd"
   },
   "outputs": [],
   "source": [
    "# the list will be filled by value of trace distance per iteration\n",
    "trace_distance = []\n",
    "\n",
    "for _ in range(400):\n",
    "    with tf.GradientTape() as tape:\n",
    "        # complex representation of parametrization\n",
    "        # shape=(2**n, 2**n)\n",
    "        ac = qgo.manifolds.real_to_complex(a)\n",
    "\n",
    "        # density matrix\n",
    "        rho_trial = ac @ tf.linalg.adjoint(ac)\n",
    "\n",
    "        # probabilities of measurement outcomes\n",
    "        p = tf.tensordot(rho_trial, data_set, [[0, 1], [2, 1]])\n",
    "        p = tf.math.real(p)\n",
    "\n",
    "        # negative log likelihood (to be minimized)\n",
    "        L = -tf.reduce_mean(tf.math.log(p))\n",
    "\n",
    "    # filling trace_distance list (for further plotting)\n",
    "    trace_distance.append(tf.reduce_sum(tf.math.abs(tf.linalg.eigvalsh(rho_trial - rho_true))))\n",
    "    # gradient\n",
    "    grad = tape.gradient(L, a)\n",
    "    # optimization step\n",
    "    opt.apply_gradients(zip([grad], [a]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "KCToEbFhd8O7"
   },
   "source": [
    "Here we plot trace distance vs number of iteration to validate the result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 298
    },
    "colab_type": "code",
    "id": "xbkSQUK6d8O8",
    "outputId": "a7ec23a6-c806-4c6d-c0ea-4a6e13d85287"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,0,'$iter$')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEMCAYAAADu7jDJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XmUlNW97vHvjx5soghHQFA60hAGm0FbbAFBuOKQo8jgmDhc1EjEY+KQaE40OU4kV5KbpYlg0LtwOA4hZjAeo5gBRY0SlNAIDQ2IICA0GBmMKMrQwL5/7GpooBtqeKt2ddXzWatWdVe/VD2+C/pxv8Pe5pxDREQkWS1CBxARkeZNRSIiIilRkYiISEpUJCIikhIViYiIpERFIiIiKVGRiIhISlQkIiKSEhWJiIikpDB0gExo166dKysrCx1DRKTZmDt37kbnXPt4ts2LIikrK6Oqqip0DBGRZsPMPoh3Wx3aEhGRlKhIREQkJSoSERFJSV6cIxGR7FZXV0dtbS3btm0LHSXvlJSUUFpaSlFRUdLvoSIRkeBqa2tp1aoVZWVlmFnoOHnDOcemTZuora2lS5cuSb9Pszu0ZWaHm9mTZvaImV0ROo+IpG7btm20bdtWJZJhZkbbtm1THglmRZGY2eNmtt7MavZ7/RwzW2pmy83s9tjLFwLPOueuBUZlPKyIpIVKJIwo9ntWFAnwBHBOwxfMrACYDJwL9AIuM7NeQCmwJrbZrnSGmjwZfve7dH6CiEjzlxVF4px7A/h4v5f7A8udcyucczuA3wCjgVp8mUCa8z/2GEyZks5PEBFp/rKiSJrQib0jD/AF0gl4DrjIzB4GXmzqD5vZODOrMrOqDRs2JBVg0CCYPRt2pXXcIyLSvGVzkTR24M455z53zn3DOXe9c25qU3/YOTfFOVfpnKts3z6u6WIOcOqpsGUL1NQcelsREYA1a9YwbNgwysvL6d27NxMnTozsvV9//XXGjBnT6M9mzZrF3XffDfir4H77299G9rmHks1FUgt8ucH3pcC6TAY49VT//NZbmfxUEWnOCgsLuf/++1myZAlvv/02kydPZvHixZG89/z58znppJMa/dmgQYMYP348ADNmzOCdd96J5DPjkc1FMgfobmZdzKwYuBR4IZMBunSBDh1UJCL54oILLuCOO+5gyJAhdOzYkVdeeSXh9zjmmGPo168fAK1ataK8vJy1a9c2uu2yZcs4/fTTqays5Pvf/z7dunUDYODAgaxatQqAtWvXUllZCUB1dTVr165lwIABdO3alddff33Pe11yySXMnDmTmTNncsstt/Dss89SUVHB+PHjOfnkkznhhBMYMmRIwv898ciKGxLN7BngdKCdmdUCdzvnHjOzG4C/AgXA4865RZnN5Ucls2Zl8lNF8tt3vgPz50f7nhUV8MADh96upqaGwYMH8+abb/Lcc88xdepUzjrrrD0/HzJkCJ999tkBf+6+++7bZ7t6q1atYt68eQwYMOCAn+3atYsrr7ySyZMn069fP2688UZ69+6Nc47Vq1fTuXNnABYsWEDfvn0BPyIZPXo0s2fPZvr06dx55528+eabe7L37duX1q1bc8opp3DffffRuXNnBgwYwPz58ykuLuaTTz6Ja38lKiuKxDl3WROv/wn4U7Lva2YjgZH1LZ+MU0+F55+HDRsgyVMtItIMfPHFF2zevJnvfve7AOzcuZM2bdrss039L+14bNmyhYsuuogHHniAI4888oCfP//88/Tq1WvP6KW8vJw2bdqwfPlyunTpsuf+jvoi2blzJ5s2beKHP/whABUVFWzcuBHwN3TW1dXRunVrAJYuXUrPnj2pq6tj69at3HrrrVx11VV7RjZRy4oiSRfn3IvAi5WVldcm+x7150lmz4YRIyIKJiJNimfkkA6LFi3i5JNPpqCgAPC/wPv06bPPNvGOSOrq6rjooou44ooruPDCCxv9vHnz5lFRUbHn++rqas4++2wWLly4ZwQCUFVVxXXXXcfixYvp1q0bxcXFALzzzjuceOKJe7L36tULgE2bNtG6dWuKioooKiqipqaGF198kXHjxvHNb36Tb33rW8nsnoPK6SKJQr9+UFCgIhHJdTU1Nfv8Yl+wYAGjR4/eZ5t4RiTOOcaOHUt5eTm33HJLk9u1bduWd999F4DZs2fz1FNPceutt/LGG2/QsmVLAJYsWcJLL73EL3/5S6ZPn87KlSvZvn07dXV1jB8/nl/84hcALFy4kBNOOAGAlStXcuyxxwL+HEz37t259NJLWbx4cdomxVSRHMLhh0OfPvCPf4ROIiLptHDhwn3OZdTU1BwwIonH3//+d55++mn69u27p5gmTJjA8OHD99luzJgxnHfeefTt25fhw4fTtm1bunXrRsuWLZk0aRJf+9rX6NOnD23btqVDhw5UV1dzxRVXMGjQILZu3cqdd97JwIEDD8h+/PHHs3HjRvr06cOxxx7LBx98wOGHH07v3r155JFHkt09B2XOubS8cTaprKx0qSy1e911fqqUTZugRTZf5ybSTC1ZsoTy8vLQMYJZs2YNF198MbNnzw7y+Y3tfzOb65yL66RKTv9aNLORZjZl8+bNKb1P//7wySewbFlEwUREGqiurt5zaKo5yukicc696JwbV38lQ7LqR7s6vCUi6TBixIi0HXbKhJwukqiUl8MRR/gT7iIisi8VSRwKCqCyUkUiItIYFUmcBgyA6mrQktIi6ZEPF/5koyj2u4okTgMGQF1d9FM3iAiUlJSwadMmlUmG1a/ZXlJSktL75PR9JFFMkVKvf3///PbbELt0W0QiUlpaSm1tLcmuHSTJKykpobS09NAbHoTuI0lAjx7QvTu89FIEoUREspjuI0mTc86B117TeRIRkYZUJAkYPhy2boXp00MnERHJHiqSBJx5pp9K/qmnQicREckeKpIEFBXBFVfACy9AEwueiYjknZwukqjm2mropptg92742c8ie0sRkWYtp4skqrm2GurSBb7xDXj4YXjvvcjeVkSk2crpIkmXH/8YSkrgIGvWiIjkDRVJEjp2hLvu8veT/PnPodOIiISlIknSTTf5mxO/+13YsSN0GhGRcFQkSSouhl/8ApYuhcmTQ6cREQlHRZKC4cP93e733APr14dOIyIShookBWZ+VPLFF3DHHaHTiIiEkdNFko77SPZ3/PFw443w6KOaYl5E8lNOF0k67iNpzF13QatWMGFCWj9GRCQr5XSRZEqbNnD99fCHP8Dy5aHTiIhklookIjffDIWFcN99oZOIiGSWiiQixxwDV18NTzwB//xn6DQiIpmjIonQ977nb06cNCl0EhGRzFGRRKh7d7joInjoIfj009BpREQyQ0USsdtug82b/eXAIiL5QEUSscpKOO00P23Krl2h04iIpJ+KJA1uuglWrIA//Sl0EhGR9MvpIsnEne2NOf98KC3VSXcRyQ85XSSZurN9f0VF/gbFV16BxYsz+tEiIhmX00US0rXXwmGHaYp5Ecl9KpI0ad/eXwr8zDOwfXvoNCIi6aMiSaMxY+Bf/9JJdxHJbSqSNDrrLGjXzk/mKCKSq1QkaVRYCCNHwrRpUFcXOo2ISHqoSNJs9Gh/p/vMmaGTiIikh4okzYYNg4ICmDEjdBIRkfRQkaTZkUdC//4qEhHJXSqSDDjjDJgzBz7/PHQSEZHoqUgyYNAgP4HjnDmhk4iIRC+niyTUXFv7GzDAP7/9dtAYIiJpkdNFEmqurf21bQs9eqhIRCQ35XSRZJPKSpg3L3QKEZHoqUgy5MQTYfVq+Pjj0ElERKKlIsmQigr/XF0dNoeISNRUJBly4on+ef78sDlERKKmIsmQDh38Sfd33w2dREQkWiqSDOrRA957L3QKEZFoqUgyqGdPFYmI5B4VSQb16AHr1sGWLaGTiIhER0WSQT16+Odly8LmEBGJUsJFYt5N6QiT6+qLRIe3RCSXFCb6B5xzzsxOMbPLgM2x17QqeRy6dfPPKhIRySUJF0nMK0Ax0B5w0cXJbS1bwnHHqUhEJLckWyQb6kchZjYqwjw5r0cPWLo0dAoRkegkc47kPOByMxte/3X0sXJXz56+SJzGcSKSI5K5aqsdsBV/WOsoYEKkiXJct27w6aewaVPoJCIi0UjmZPuTwJNpyJIXvvIV/7xiBbRrFzaLiEgUDlkkZnZcnO/1iXPu0xTzRMrMRgIju9VfLpUFunb1zytWQP/+YbOIiEQhnhFJPKMPBzwBPJVSmog5514EXqysrLw2dJZ6Xbr45xUrwuYQEYnKIYvEOTcsE0HyxZe+BB07qkhEJHfk9KGtbNW1q4pERHJHTh/aylZdu8Kbb4ZOISISDR3aCqBrV/j1r2HHDiguDp1GRCQ1CV/+a2bD939Nc20lpmtX2L0bVq/eO/+WiEhzlcwUKe0jT5Fn6u8lef99FYmINH/J3pAoKWh4L4mISHOX1KSNZjYAGAN8CcA5d02UoXJdx45QUgLLl4dOIiKSumRXSLwW+AS4B1gZWZo80aIFlJfDokWhk4iIpC7ZIvkIKAF2Ax2ii5M/+vaFhQtDpxARSV2yRTIVeBj4Pn6RK0lQnz6wbh18/HHoJCIiqUlmPZLBQBnQE/grsCPiTHmhb1//vGBB2BwiIqlKZkRyFNANvy5J/UMS1L+/P1cyY0boJCIiqUm4SGIz6p4C1AEbYw9J0FFHweDBMG1a6CQiIqlJds32V4Bi/GhEi8Ym6fzz4dZbYfx4+OIL/7j+eujVK3QyEZH4mUti8XAzG14/LYqZjXLOvRB5sghVVla6qqqq0DEOsG0bDB0Kc+bsnXOruBhmzdp7DkVEJAQzm+ucq4xn22ROtp8HXG5mw+u/TvQ9xCspgZkzYd482LjRT5lyxBEwdiwk0e8iIkEkVCRm1gI4BvgCf1jrKGBCGnLljeJiqKiAVq2gtBTuvdePUF56KXQyEZH4xF0kZnYD/kbEHwEDgWLn3NPOuQVmVpCugPnmyiuhUyd4+OHQSURE4pPIiORWoK9z7ljgq8BgM7s79rO3I0+WpwoL4eqr4S9/gbVrQ6cRETm0RIpkC7AewDn3T2AscFHsZ8le/SWNuOIKv17JC1l9CYOIiJdIkTwM/N7M6lfQOA74PPa1Dm1F6Pjj/TolKhIRaQ7iLhLn3EP4NdmnmNm/gOXAR2Y2DZibpnx5yQxGjYJXX4XPPgudRkTk4BK6ass590fn3BnA0cCJwEPAXYDWI4nYqFF+Tffp00MnERE5uKRm/3XO1TnnFjnnpjvn3nHJ3NWYJDPramaPmdmzmfrMEAYP9tOo6PCWiGS7ZKeRT4qZPW5m682sZr/XzzGzpWa23MxuP9h7OOdWOOfGpjdpeIWFcPbZ8MorujlRRLJbRosEeAI4p+ELsXtQJgPnAr2Ay8ysl5n1NbNp+z2OznDeoIYN82uWLFsWOomISNMSvmzXzEbGZgBOmHPuDTMr2+/l/sBy59yK2Pv/BhjtnPsJMCKZz8kVw4b559degx49wmYREWlKMiOSeyPO0AlY0+D72thrjTKztmb2/4CTzOwHB9lunJlVmVnVhg0bokubQd27+7vcX301dBIRkaYlcyOhRZyhsfdr8qyAc24T8B+HelPn3BRgCvjZf5NOF5CZH5VMn+7Pk1jUe15EJALJjEii/qVcC3y5wfelwLqIP6PZGjYM1q+HxYtDJxERaVymT7Y3Zg7Q3cy6mFkxcCmgi15jzjjDP+vwlohkq0xf/vsM8BbQ08xqzWysc24ncAPwV2AJ8Dvn3KKIPm+kmU3ZvHlzFG8XRFmZf7z2WugkIiKNS3iFRDN72Tl3dprypEW2rpAYr2uugeef94tftciGMaSI5Ly0rpDY3EokFwwdCv/6FyxZEjqJiMiB9P+3zcDgwf551qywOUREGqMiaQa6dYP27eHvfw+dRETkQEkViZndbGaPxL6+M9pIsj8zGDRIRSIi2SnZEclX2Hs3equIskQuF67aqjd4MCxfDh99FDqJiMi+ki0SB7Q0sz7AsRHmiZRz7kXn3LjWrVuHjpKy+vMkb70VNoeIyP6SLZL78VObjAF+GF0caUq/flBcrMNbIpJ9kplrC+fcauCg64ZItEpKoLJSRSIi2SepIondoe6AI4DWzrn/FWkqadTgwTBxImzb5otFRCQbJLvU7mXOucudc6OA/4k4U2Ry6WQ7+Cu3duyAuXNDJxER2SvZy3+Hxx6jgH4RZ4pMLp1sBzj1VP/89tthc4iINJTUoS2gfex5O3BbRFnkEDp08BM4zp4dOomIyF7JLLVr+PMik9KQRw5hwABNlSIi2SWZSRsdcIqZXVZ/iCsNuaQJAwfCmjWwTkt/iUiWSPY+kreBYvwhrm7RxZFDGTDAP+vwlohki4SLxMzOA04FNgAbgUFRh5KmnXQSFBWpSEQkeyQzImkHbMWPRo4CJkSaKEK5dvkv+PtHKio0VYqIZI+EisTMWgBTnXPXOueedM497ZxbkKZsKcu1y3/rnXaaH5Fs2xY6iYhIAkViZjcAHwEfmFm1mY1t8LOCdISTxg0bBtu3634SEckOiYxIbgX6Ouc6AecAp5nZ3bGf6VdaBg0Z4tdunzEjdBIRkcSKZAuwHsA59yEwFrgo9rNkb2yUJLRpA6efDlOnwu7dodOISL5LpEgeBn5vZvWX+x4HfB77Woe2Muyaa2DlSnjuudBJRCTfmb+/MM6NzUYDNwMn4VdGnIYfjWxwzn0jLQkjUFlZ6aqqqkLHiNT27X424JoaOPNM6NoVevSASy6Bjh1DpxOR5s7M5jrnKuPaNpEiafABRUAPoBP+XpJ5Lpk3ypBcLBLwy+7ec49fo+SDD+DTT6F1a5g2zV/ZJSKSrEiLxMyOi/NzP3HOfRrnthlhZiOBkd26dbt22bJloeOk3ZIlMGqUn2r+3XehZcvQiUSkuYq6SF6L430c8IRz7ql4PjTTcnVE0pjXX/eXB0+cCDfdFDqNiDRXaT+01dzkU5GAX99d06iISCoSKZJDXrbbnA9t5avLL4f//E9YscKfhBcRSad47v94Mo5tHPAEkJWHtvLNiBG+SF59VUUiIul3yCJxzg3LRBCJTs+ecPTR8MYb8M1vhk4jIrku2fVIJIuZwdCh8Le/hU4iIvlARZKjBg6E1athw4bQSUQk16lIclRFhX+urg6bQ0Ryn4okR514on9WkYhIuuV0keTiConxatcOOnWCefNCJxGRXJfTRZKrKyTGq08fP22KiEg65XSR5LuePeG99yAPJi8QkYBUJDmsZ0/YsgXWrQudRERymYokh/Xo4Z/fey9sDhHJbSqSHNazp39eujRsDhHJbSqSHFZaCocd5idvFBFJFxVJDjODzp1h1arQSUQkl6lIclznzn4ZXhGRdFGR5DgViYikm4okx5WVwUcfwdatoZOISK5SkeS4zp398+rVYXOISO7K6SLJ57m26tUXiQ5viUi65HSR5PtcW7C3SHTlloikS04XicCxx0JhoUYkIpI+KpIcV1job0xUkYhIuqhI8oBuShSRdFKR5IGyMo1IRCR9VCR5oHNnWLsWduwInUREcpGKJA+UlfnFrWprQycRkVykIskDupdERNJJRZIHysr8s064i0g6qEjyQGmpn1JeIxIRSYfC0AEk/YqLoVOnvSOSHTtg6lTYtg2+/nU46qig8USkmdOIJE+UlcH778Pu3TBiBFxzDXzrW3453pkzQ6cTkeZMRZInTj4Z5s6FBx+El1/2z++840cjZ50FL7wQOqGINFcqkjwxdKhfk+Q73/Fff/vbcNJJMGsWnHACXHopzJ8fOqWINEcqkjwxZMjer++91598B2jbFqZNgzZt4Oqroa4uSDwRacZUJHmifXt47DF/PuS00/b92dFHw6RJUF0Nv/51mHwi0nyZcy50hrSrrKx0VVVVoWNkNef8oa6tW2HxYigoCJ1IREIys7nOucp4ts3pEYlWSIyfGdxxB7z3Hjz+uC8WEZF45HSRaIXExFx4IVRUwLhxvlREROKR00UiiWnRAt54w590nzAB3nwzdCIRaQ5UJLKPVq3goYegQwe46y4d4hKRQ1ORyAFatoQ774TXX/dTqDz/POzaFTqViGQrFYk06vrrYcwYXyIXXACnnw66ZkFEGqMikUa1aAFPPQWffw7//d8wezZ87Ws61CUiB1KRyEEVFfmT7xMnwvTp8OSToROJSLZRkUhcrrsO+veHu++G7dtDpxGRbKIikbi0aAE/+hGsXg2/+lXoNCKSTVQkErevftXPFDxpks6ViMheKhKJmxncdBMsWAB/+1voNCKSLVQkkpDLL/dTz0+cGDqJiGQLFYkkpGVLf+L9j3+ElStDpxGRbKAikYRdf70/+f7LX4ZOIiLZQEUiCSsthYsv9gtlbdkSOo2IhKYikaTcfLOfMmXSpNBJRCQ0FYkk5dRT/folP/4xzJ0bOo2IhKQikaRNnuzXez/zTD8v1+7doROJSAgqEklax45+qvleveCqq2DQIJgzJ3QqEck0FYmkpEsXmDkTnngCPvjAl8nTT4dOJSKZpCKRlLVo4Uck774LQ4fClVfCz34WOpWIZIqKRCLTujX8+c9+VcXbboMf/EBzconkg8LQASS3FBfD1KnQpg389Kewbh08+CAceeS+2+3eDTNmwLPPwqZNUF4OI0fCKaf4Ob1EpPnQiEQiV1AADz/s1y751a+gd2/4+c/9lCrvv+/n6Sov97MJ//a3sGgRTJgAAwZAv35+RcZt20L/V4hIvMzlwbGHyspKV1VVFTpGXpo9G265BWbN2vf1gQPh29+GSy6Bww6Djz/2o5MHH4SaGj8x5Lnn+pP3ZWVQUuILavdu2LnTPwoK/OstW/rnFi3go4/8KKj+8eGHsH69P8RWUADt2/v3q3906QKdO/sM8dq50xdd/WP79gO/377dj6xatNj7KCg48OtER1/5tr2k7oQT/N+1RJnZXOdcZVzbNrciMbPzgfOAo4HJzrnph/ozKpLwFi3aWyaDBvlRSmOc85cUP/oovPwybNiQ/Gcefjgccwx06OD/Ie3c6Utl9Wqoq9u7nZk/v/OlL/lHQQHs2OG3qX80/F73y0hzsn27P+ScqESKJKPnSMzscWAEsN4516fB6+cAE4EC4FHn3E+beg/n3PPA82b2b8B9wCGLRMLr3bvp8mjIDIYN84/du/2oYvVq/4t81y7/S76w0D/v2rXvSGDnTn+DZKdOvkD2Py9Tb9cu/76rVvnDbStX+hHRF1/4x86d/h9eUZF/NPy6qMiPfkpK/Cim/uv9X6v/h7t7t3/s2nXg17t2JbYPE/1/vua+vUSjMAO/5TM6IjGzocAW4Kn6IjGzAuA94GygFpgDXIYvlZ/s9xbXOOfWx/7c/cBU59w7h/pcjUhERBKTtSMS59wbZla238v9geXOuRUAZvYbYLRz7if40cs+zMyAnwJ/jqdEREQkvbLhqq1OwJoG39fGXmvKjcBZwMVm9h9NbWRm48ysysyqNqRyoF1ERA4qG+4jaew6jiaPtznnJgGHnLzcOTcFmAL+0FbS6URE5KCyYURSC3y5wfelwLpAWUREJEHZUCRzgO5m1sXMioFLgRcCZxIRkThltEjM7BngLaCnmdWa2Vjn3E7gBuCvwBLgd865RRF93kgzm7J58+Yo3k5ERBrR7G5ITIYu/xURSUwil/9mw6EtERFpxvJiRGJmG4APkvzj7YCNEcaJinIlRrkSl63ZlCsxyebq7JxrH8+GeVEkqTCzqniHd5mkXIlRrsRlazblSkwmcunQloiIpERFIiIiKVGRHNqU0AGaoFyJUa7EZWs25UpM2nPpHImIiKREIxIREUmJiqQJZnaOmS01s+VmdnvgLKvMbKGZzTezqthrR5nZy2a2LPb8bxnK8riZrTezmgavNZnFzH4Q24dLzezfM5zrHjNbG9tv881seIBcXzaz18xsiZktMrObY68H3WcHyRV0n5lZiZn9w8yqY7nGx14Pvb+ayhX871jsswrMbJ6ZTYt9n9n95ZzTY78HflGt94GuQDFQDfQKmGcV0G6/134G3B77+nbg/2Yoy1CgH1BzqCxAr9i+OwzoEtunBRnMdQ/wvUa2zWSuY4B+sa9b4Rdx6xV6nx0kV9B9hp8N/IjY10XAbGBgFuyvpnIF/zsW+7xbgF8D02LfZ3R/aUTSuD2LbTnndgC/AUYHzrS/0cCTsa+fBM7PxIc6594APo4zy2jgN8657c65lcBy/L7NVK6mZDLXhy62AJtz7jP8fHKdCLzPDpKrKZnK5ZxzW2LfFsUejvD7q6lcTcnY3zEzKwXOAx7d7/Mztr9UJI1LdLGtdHPAdDOba2bjYq91cM59CP6XAnB0sHRNZ8mG/XiDmS2IHfqqH94HyWV+ddCT8P83mzX7bL9cEHifxQ7TzAfWAy8757JifzWRC8L/HXsA+D6wu8FrGd1fKpLGJbTYVgYMds71A84Fvm1mQwNmSUTo/fgw8BWgAvgQuD/2esZzmdkRwB+A7zjnPj3Ypo28lrZsjeQKvs+cc7uccxX4tYn6m1mfg2weOlfQ/WVmI4D1zrm58f6RRl5LOZeKpHFZtdiWc25d7Hk98D/4oehHZnYMQOx5fah8B8kSdD865z6K/ePfDTzC3iF8RnOZWRH+l/VU59xzsZeD77PGcmXLPotl+QR4HTiHLNhfjeXKgv01GBhlZqvwh+DPMLNfkeH9pSJpXNYstmVmh5tZq/qvga8CNbE8V8U2uwr4Y4h8MU1leQG41MwOM7MuQHfgH5kKVf8PKeYC/H7LaC4zM+AxYIlz7ucNfhR0nzWVK/Q+M7P2ZtYm9nVL4CzgXcLvr0Zzhd5fzrkfOOdKnXNl+N9Trzrn/jeZ3l/puoqguT+A4fgrWd4H/itgjq74qyyqgUX1WYC2wAxgWez5qAzleQY/hK/D/9/N2INlAf4rtg+XAudmONfTwEJgQewf0DEBcp2GP3SwAJgfewwPvc8OkivoPgNOAObFPr8GuOtQf98D5wr+d6zB553O3qu2Mrq/dGe7iIikRIe2REQkJSoSERFJiYpERERSoiIREZGUqEhERCQlKhIREUmJikRERFKiIhFJIzMbZGbjzazUzL4eOo+MgkSJAAABPUlEQVRIOuiGRJEMMLOr8Gva3JbAnylwzu1KYyyRSKhIRNLIzH4PTMTPdfQJ8Bl+Tibw0393wk//PcY5tzS2/Rr8tO4znHP/J/OpRRJTGDqASI7rg5+LaQ5+Jb2a2Ky7fwHGOefejy3PejvwDaAvfiLFYcESiyRIRSKSJmZWAhQ55zabWU/8JHngV6vrDfzBT8JLIfBmbPujgB+FyCuSLBWJSPr0BhabWVtgs3OuLvb6ifhZnB9ruLGZnQzMds7tzHBOkZToqi2R9OmLn168C/suHvQh8O9m1gLAzPrG1gep316kWVGRiKRPfTG8C7QzsxozGwQ8jv+3tyS2Bvhtzl/1oiKRZklXbYmISEo0IhERkZSoSEREJCUqEhERSYmKREREUqIiERGRlKhIREQkJSoSERFJiYpERERS8v8BV8UXTxs1P6YAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(trace_distance, 'b')\n",
    "plt.legend([r'$n=$' + str(n) + r'$\\ qubits$'])\n",
    "plt.yscale('log')\n",
    "plt.ylabel(r'$||\\varrho_{\\rm true} - \\varrho_{\\rm trial}||_{\\rm tr}$')\n",
    "plt.xlabel(r'$iter$')"
   ]
  }
 ],
 "metadata": {
  "accelerator": "GPU",
  "colab": {
   "collapsed_sections": [],
   "name": "state_tomography.ipynb",
   "provenance": []
  },
  "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}