Quantum_physics/01_GHZ_circuit.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulate preparing a GHZ state and compute covariance\n",
    "\n",
    "In this notebook, we simulate preparation of the $n$-qubit GHZ state\n",
    "$$\n",
    "\\vert\\text{GHZ}\\rangle = \\frac1{\\sqrt2} \\Bigl( \\vert0\\dots0\\rangle + \\vert1\\dots1\\rangle \\Bigr)\n",
    "$$\n",
    "via the unitary circuit\n",
    "$$\n",
    "\\vert\\text{GHZ}\\rangle = \\text{CNOT}_{n-1,n} \\dots \\text{CNOT}_{1,2} H_1 \\vert 0 \\dots 0 \\rangle\n",
    "$$\n",
    "where the Hadamard and CNOT gates are given by\n",
    "\\begin{align}\n",
    "H &= \\frac1{\\sqrt2} \\begin{pmatrix}1 & 1 \\\\ 1 & -1\\end{pmatrix}, \n",
    "&\n",
    "\\text{CNOT} &= \\vert0\\rangle\\langle0\\vert \\otimes \\mathbb 1 + \\vert1\\rangle\\langle1\\vert \\otimes \\sigma_X, &\n",
    "&\n",
    "\\sigma_X = \\begin{pmatrix}0 & 1\\\\1 & 0\\end{pmatrix}.\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "Applicable keywords:\n",
    "\n",
    "- Efficient classical simulation of a unitary circuit\n",
    "\n",
    "- Quantum computation\n",
    "\n",
    "- Quantum information\n",
    "\n",
    "- Matrix product states, tensor trains\n",
    "\n",
    "- Matrix product operators (MPO), matrix product density operators (MPDO)\n",
    "\n",
    "- Tensor network methods\n",
    "\n",
    "The code below uses Python and `mpnum` (http://mpnum.readthedocs.io/). This notebook is part of https://github.com/milan-hl/mpnum-examples/. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Definitions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Open ``init.ipynb`` to see which MPSs and MPOs it defines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%run init.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Let's do 6 qubits\n",
    "n_sites = 6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Hadamard and CNOT gates are constructed as follows:\n",
    "\\begin{align}\n",
    "  H &= \\frac1{\\sqrt 2}(\\sigma_X + \\sigma_Z) \\\\\n",
    "  \\text{CNOT} &= \\vert0\\rangle\\langle0\\vert \\otimes \\mathbb 1 + \\vert1\\rangle\\langle1\\vert \\otimes \\sigma_X\n",
    "\\end{align}\n",
    "where $\\sigma_{X,Y,Z}$ are Pauli matrices.\n",
    "\n",
    "With `mpnum` and the definitions from `init.ipynb`, we have e.g.\n",
    "\\begin{align}\n",
    "  \\vert1\\rangle\\langle1\\vert \\otimes \\sigma_X &= \\texttt{mp.chain((mdd, mpx))}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Define Hadamard gate (single qubit)\n",
    "mh = (mpx + mpz) / 2**0.5\n",
    "\n",
    "# Define CNOT gate (two qubits)\n",
    "mcnot = mp.chain((muu, mid)) + mp.chain((mdd, mpx))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we compute\n",
    "$$\n",
    "\\Vert H H^\\dagger - \\mathbb 1 \\Vert_F\n",
    "$$\n",
    "where the Norm is the Frobenius norm $\\Vert A \\Vert_F^2 = \\operatorname{tr}(A^\\dagger A)$ and $A^\\dagger$ denotes the Hermitian adjoint (cojugate transpose). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.1401849173675503e-16"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Verify that Hadamard is unitary\n",
    "mp.normdist(mp.dot(mh, mh.adj()), mp.factory.eye(1, ldim))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Verify that CNOT is unitary\n",
    "mp.normdist(mp.dot(mcnot, mcnot.adj()), mp.factory.eye(2, ldim))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prepare GHZ as MPS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a pure state, containing spin up on `n_sites` sites:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "psi = mp.chain([mup] * n_sites)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Number of sites of the pure state:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(psi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Number of sites on which the two gates act:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, 2)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(mh), len(mcnot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1, 1, 1, 1, 1)\n",
      "(2, 1, 1, 1, 1)\n",
      "(2, 2, 1, 1, 1)\n",
      "(2, 2, 2, 1, 1)\n",
      "(2, 2, 2, 2, 1)\n",
      "(2, 2, 2, 2, 2)\n"
     ]
    }
   ],
   "source": [
    "psi = mp.chain([mup] * n_sites)\n",
    "# Apply Hadamard at site 0 (the first site)\n",
    "psi = mp.partialdot(mh, psi, 0)\n",
    "\n",
    "print(psi.ranks)\n",
    "\n",
    "for pos in range(n_sites - 1):\n",
    "    # Apply CNOT at sites (pos, pos + 1)\n",
    "    psi = mp.partialdot(mcnot, psi, pos)\n",
    "    print(psi.ranks)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output shows that each application of a CNOT gate increases the bond dimension (=rank) by a factor of two (=the bond dimension/rank of the CNOT gate)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function computes\n",
    "$$\n",
    "\\text{cov}_\\rho(A_i, B_j) = \\text{tr}(A_i B_j \\rho) - \\text{tr}(A_i \\rho) \\text{tr}(B_j \\rho)\n",
    "$$\n",
    "for a pure state $\\rho = \\vert\\psi\\rangle\\langle\\psi\\vert$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cov_ab(psi, a, b, pos1, pos2):\n",
    "    if pos1 == pos2:\n",
    "        return np.nan  # Return not a number\n",
    "    # Todo: Verify that imaginary parts are small\n",
    "    # Apply a on site pos1\n",
    "    ept_a = mp.inner(psi, mp.partialdot(a, psi, pos1)).real\n",
    "    # Apply b on site pos2\n",
    "    ept_b = mp.inner(psi, mp.partialdot(b, psi, pos2)).real\n",
    "    ept_ab = mp.inner(\n",
    "        psi,\n",
    "        mp.partialdot(a, mp.partialdot(b, psi, pos2), pos1)\n",
    "    ).real\n",
    "    cov = ept_ab - ept_a * ept_b\n",
    "    return cov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "covariances = [\n",
    "    [\n",
    "        cov_ab(psi, mpz, mpz, pos1, pos2)\n",
    "        for pos2 in range(n_sites)\n",
    "    ]\n",
    "    for pos1 in range(n_sites)\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATMAAAEKCAYAAAB+LbI7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAF3NJREFUeJzt3X20XXV95/H3597cCEOQ1EYgEBSnYgqigxJjLRapQBuLBbsGV7F2lC4wM8tBy6iVjM5Cpa1SO2p11bam1AXaVnyuqUlFRGOHESFBQXlohogoMYGAykMEhOR+5o+9b3o4nHPvOTnn3HPu735ea+11z977d377e05yv/f3sB9km4iIuW5s2AFERPRDkllEFCHJLCKKkGQWEUVIMouIIiSZRUQRkswioghJZl2Q9B5J5w/p2NdJenYP718u6duSHpT0xg7fc4ekU/b1mBGzKcmsQ5KeCrwG+Ei9fqGkXU3LI5Is6Xc7rLObOv43cFEPH+GtwEbbB9r+UItYBp64JJ0l6VpJP5O0s379eklqF4OksyVdXb9+dYvva1f9fV3YZSxdfd4k9tGXZNa5s4ENth8GsH2R7UVTC3AocAOwAfhcJxV2Wcc64NclLd3H+J8O3LyP7+2ZpDcDHwT+nOpzHgL8N+AEYGEnddj+h8bvq/7OzgfuBv52MJHHnGG7yAV4FXA9cD/wPeAkQMAFwA+A+4BPAQfVi4ElDe8/huqX5Mn1+leB329zrP2BjcBVwH77GO+MdQBXAq+dpo6j6zruo0pcpzfEvgd4BNgFPKvpfR8HJoGH6/1vrbffAbwF+E79PX5yKjbgMOCzwD3A94E3ThPXQcDPgP88w3dwB3BK07azgavblH8e8CBw0jR1XgD8qC63BTh5ms+7pv6/8iBwC/A7M3w/HX8HWQa/DD2AgXwoeHP9y/x8qtbnc4AjgT+pf9mXAvsBXwDeV7/nTuDEhjo+D7ypYf0e4AUtjrUQ+BfgG8CifYy3ozqADwHvb7NvAtgKvK2u76X1L+Xyev9G4Nxp6m6VSO4Arqt/aZ8C3ErVmhqj+kNxYX2s/wjcDvxmm7pXAbuBBTN8Dx0nM2BxnXgumKa+5fW/62H1+pHAL01zrFfWn3UM+F2qBLy0Vfluv4Msg1+K62bWY1vvAH7P9rdsT9r+LtVf1TfU23fYfgT4DLCifutNVC0bJK2kSoQfbqh6MVVyaDzWAuByqi7Ty2zvatj3eknP6iDetnXU+8+VdGy9+mAdRyu/AiwCLrb9qO2vAl+kaqH24kO2t9v+CfDPwHHAC4CnuuomP2r7dqpu3llt6lgC3Gt7d8Pn+oak+yQ9LOnEhrL/VG+/T9J9wF81V1aPsV1G9W/23mli3wM8CThG0oTtO2x/r11h25+uP+uk7U8CtwEr2xTv9juIAVsw7AAG4BTgu7ZvbNr+a/X27Q3blgA76tc3UXUtAd4DvNP2zxvK/hQ4cGpF0hhwKdVf/5fYvr/xYLaf8EvYbKY66nouaVg9kKoL2cphwJ22Jxu2/QA4fKY4ZnBXw+uH6uM8HTisTjZTxoH/06aOHwNLJC2YSmi2fxVA0jYeP3b7CttfmVqRdDZwblN9FwDHAsfbbnvbF9tb69nndwLPlnQFVWt7e6vykl4DvImqBQfVH4clbarv9juIASuuZUbVHWr1C/9UqnGfRmcAV9evb6L6C34KVTf0Y01lvwM0trT+mqo1dIrte5sPJqmT/9TT1lHXs7Fh9WigOUlP2Q4cUSfIKU+jGi/qRDf3groT+L7txQ3LgbZ/q035a4CfU33fPZF0EvB24Ezb7RL7Xrb/0faLqZKPgT+b2tVU79OpWlbnAb9oezHV/wm1Kk/330EMWInJ7NvAiyX9J1WOknQ0sAl4kaRfkrRI0kVUXbuP1u+bapm9G3ib7T1N9W4AXgIg6f3Ay4CTbe9oKoekJcDOhvVLJV3aVGbaOuoyB1J3bSU9CTieahKglWupxnjeKmmi/qX/baoubCfuphr36cR1wAOSLpC0v6RxScdKekGrwnXSeRfwV5LOrL//MUnHAQd0eEzqmdzLgfNtf7uD8sslvbT+7h6hGmqY+ndt/rwHUCWse+r3/gFV64825bv6DmIWDHvQbhAL8EdUfzl3USWp59fb30zVrbwP+DT14G69b3+q/+jXtKlzCbCNajLBwKN1/Y3L3VR/IF4KvKPhvVcBr2tYP3amOupyLwLeU79+JfC5GT73s4GvU7VA987G1fs2Mv0EwBnAD+vv5i31tjt4/KD3O4G/r18fBnyCqhv6U+CbNA2otzjGq6mSwENUSeNaYDWwsNXx6m1nU08AUA22u8V3tgv4mxbHe259vAeBn1CNIR42zef907rcvcD76+/y3GnKd/0dZBncovofJTog6d3ATtt/MUO584Ef2P68pIVUXcPn2n6sy+OtBnbZ/kdJ1wLn2L5pX+OPKFmJEwADY/ttHRZ9DtXMH7YfpZ4l3QfPAdbW9bxwH+uImBfSMhth9eD/qd226CLmoxInAIog6Z+orqVMIovoQFpmEVGEtMwiogizMQEwUk2/9RPLhx1CRN+d9tgWzVyqvfUTyzv+Pe31WIOSlllEFCGnZkQEmhjJxlZXkswigvH9x4cdQs+SzCKCsQVpmUVEAdLNjIgipGUWEUVIyywiipCWWUQUYXzh3D/lNMksItBYWmYRUQCNp2UWEQUYG0/LLCIKkG5mRBQhEwARUQSNzYNkJumXqR6zdTjVvcm2A+ts3zrg2CJilpTQzZw2HUu6gOqhq6J6/uCm+vUnJK0ZfHgRMRvGxtXxMqpmapmdAzy7+aEa9dO4bwYubvWm+nmPqwE+8pGPsHr16j6EGhGDUkLLbKZkNkn11OYfNG1fWu9ryfZa6uc9MmK3zY6IJ+r3mJmkVcAHgXHgEtsXN+1/GnAZsLgus8b2hl6OOVMyOx+4StJtwJ31tqcBzwTO6+XAETE6xif6l8wkjQMfBk4FtgGbJK2zfUtDsf8FfMr2X0s6BtgAHNnLcadNZra/JOlZwEqqCQBNBWd7Ty8HjojR0edu5kpgq+3bASRdTjWJ2JjMDDy5fn0Q1cRiT2aczbQ9CXyz1wNFxOjqppvZOCZeW1sPLU05nH/vyUHVAHphUzXvBL4s6Q3AAcAp3cTbSs4zi4iuWmZNY+Itq2v1tqb1VwGX2n6fpBcBH5d0bN142idJZhHR727mNuCIhvVlPLEbeQ6wCsD2NZL2A5YAO/f1oHP/tN+I6NnYgvGOlw5sAo6S9AxJC4GzgHVNZX4InAwg6WhgP+CeXj5DWmYR0deTYW3vlnQecAXVaRcftX2zpIuAzbbXAW8G/lbS/6Dqgp5tu6fTuJLMIqLvJ83W54xtaNp2YcPrW4AT+nnMJLOImB8XmkdE+ebD5UwRMQ8kmUVEETqcpRxpSWYRkTGziCiE0s2MiAJkzCwiipBuZgfWTywf9CEiokdpmc1Rpz22Zdgh7JVkH6Mgs5kRUYS0zCKiDBkzi4gSKKdmREQJMpsZEUVQJgAiogSZAIiIIkjpZkZECdIyi4gSZAIgIoqQMbOIKILGM5sZESVINzMiSpArACKiDGmZRUQJMgEQEWUo4KTZff4Ekv6gn4FExPBofLzjZVT1ko7f1W6HpNWSNkva/KXJ+3o4RETMijF1voyoabuZkr7TbhdwSLv32V4LrAVYP7Hc+xxdRMyK+XAFwCHAbwI/bdou4BsDiSgiZt88ODXji8Ai2zc075C0cSARRcTsK71lZvucafb9Xv/DiYhhGOWB/U7l1IyIKOLUjCSziBjpWcpOzf10HBE9k8Y6XjqrT6skbZG0VdKaacqdKcmSVvT6GdIyi4i+tswkjQMfBk4FtgGbJK2zfUtTuQOBNwLX9uO4aZlFRDVm1ukys5XAVtu3234UuBw4o0W5PwbeCzzSj4+QZBYRMD7e+TKzw4E7G9a31dv2kvQ84AjbX+zXR0g3MyK6ms2UtBpY3bBpbX3Vz94iLd6290ogVQNvHwDO7i7I6SWZRURXY2aNlyu2sQ04omF9GbC9Yf1A4FhgY31TyEOBdZJOt72540CaJJlFRL/PM9sEHCXpGcCPgLOAvSfZ274fWLL30NXVRG/pJZFBxswiAqprMztdZmB7N3AecAVwK/Ap2zdLukjS6YP6CGmZRUTfr820vQHY0LTtwjZlT+rHMZPMIgLGcm3mnLR+YvmwQ4gYLaXfNSNmx2mPbRl2CHsl0c9T8+B+ZhExH+SuGRFRhLTMIqIEzs0ZI6II6WZGRBGSzCKiBM6YWUQUIS2ziChCWmYRUYLMZkZEGdLNjIgSOMksIoqQMbOIKEFaZhFRhtzPLCJKUMJJszO2LSX9sqSTJS1q2r5qcGFFxKzq70OAh2LayCS9EfgC8AbgJkmNTyV+9yADi4jZY9TxMqpm6ma+Djje9i5JRwKfkXSk7Q/S+kGfwOMfEnre2MGsGlvcp3AjYhDmwwTAuO1dALbvkHQSVUJ7OtMks8aHhK6fWO525SJiRBSQzGb6BHdJOm5qpU5sL6d6gOdzBhlYRMyeybHxjpdRNVMyew1wV+MG27ttvwY4cWBRRcTs6uNDgIdl2m6m7W3T7Pu//Q8nIoZhPoyZRcQ8MMqzlJ1KMouItMwiohAjPBbWqSSziGBSoztL2akks4hINzMiypAJgIgoQlpmEVGEEm4BlGQWEZkAiIgylDBmNvc7yhHRM2us46UTklZJ2iJpq6Q1LfY/SdIn6/3X1rcY60mSWUT09eaMksaBDwMvA44BXiXpmKZi5wA/tf1M4APAn/X6GZLMIqLfLbOVwFbbt9t+FLgcOKOpzBnAZfXrzwAnS73NQmTMbASsn1g+7BBinutmzKzxTtK1tfUNWaccDtzZsL4NeGFTNXvL2N4t6X7gF4F7uwj7cZLM4glOe2zLsEN4nCT7wZvsopPWeCfpNlplxuY7TndSpitJZhGB+zvitA04omF9GbC9TZltkhYABwE/6eWgGTOLiH4/nWkTcJSkZ0haCJwFrGsqsw54bf36TOCrttMyi4je9PM8s3oM7DzgCmAc+KjtmyVdBGy2vQ74O+DjkrZStcjO6vW4SWYR0feTZm1vADY0bbuw4fUjwCv7ecwks4go4gqAJLOIYNJzf/g8ySwi0jKLiDIkmUVEEewks4gowGRaZhFRgkwAREQRMmYWEUXImFlEFGFetMwkrQRse1N9t8hVwL/VlytERAGKb5lJegfVrW8XSLqS6gZrG4E1kp5n+08HH2JEDNrksAPog5laZmcCxwFPAu4Cltl+QNKfA9cCLZNZ450ozxs7mFVji/sXcUT03XyYzdxtew/wkKTv2X4AwPbDktom88Y7Ua6fWN7TPYoiYvCK72YCj0r6D7YfAo6f2ijpIMpomUYE82MC4ETbPwew3Zi8Jvj3u0RGxBw3WUD/adpkNpXIWmy/lx6eohIRo2U+tMwiYh6YD2NmETEP7Ekyi4gSpGUWEUXo7SFvoyHJLCIyARARZSj+1IyImB8mJ9Myi4gC5LbZEVGETABERBFyakZEFCETABFRhHQzI6IIuZwpirR+YvmwQ4hZlpZZxCw57bEtww5hrxKTfZJZRBRhMt3MiChBCS2zuf9Ilojo2Z7JzpdeSHqKpCsl3Vb//IVpyj5Z0o8k/WUndSeZRQS2Ol56tAa4yvZRwFX1ejt/DHy904qTzCICu/OlR2cAl9WvLwNe0aqQpOOBQ4Avd1pxkllEMOnOlx4dYnsHQP3z4OYCksaA9wF/1E3FmQCIiK5aXJJWA6sbNq2tH/w9tf8rwKEt3vr2Dg/xemCD7Tulzru1SWYR0VUyqxPX2mn2n9Jun6S7JS21vUPSUmBni2IvAn5N0uuBRcBCSbtsTze+lmQWEb3PUnZhHdUDxC+uf36huYDtV0+9lnQ2sGKmRAYZM4sIYHKy86VHFwOnSroNOLVeR9IKSZf0UnFaZhExayfN2v4xcHKL7ZuBc1tsvxS4tJO6k8wioogrAJLMIiI3Z4yIMrirptloXpTe9QSApI8NIpCIGJ49ezpfRtW0LTNJ65o3Ab8uaTGA7dMHFVhEzJ75MGa2DLgFuAQwVTJbQXWpQVuNZwifN3Ywq8YW9x5pRAxMCWNmM3UzVwDXU12GcL/tjcDDtr9uu+3V7LbX2l5he0USWcTom8ULzQdm2paZ7UngA5I+Xf+8e6b3RMTc466aZqM5AdBRYrK9DXilpNOABwYbUkTMtlm8nGlgumpl2V4PrB9QLBExJJMFDJqlyxgRIz0W1qkks4hIMouIMkwWkM2SzCICz7cJgIgo0549aZlFRAG6u9B8NCWZRUQRlzMlmUVEl1cAjKYks4jIqRkRUYZcARARRZjMbGZElCAnzUZEEXJqRsQsWT+xfNghFC1jZhHz1GmPbRl2CH1VQMMsySwiYE8Bd2dMMouInDQbEWVIMouIIhSQy5LMIiIts4goRM4zi4giZDYzIoqQbmZEFCHJLCKKUMKF5mPDDiAihs+T7njphaSnSLpS0m31z19oU+69km6WdKukD0nSTHUnmUUEtjteerQGuMr2UcBV9frjSPpV4ATgucCxwAuAl8xUcbqZEcGe3bM2m3kGcFL9+jJgI3BBUxkD+wELAQETwN0zVZyWWUR01TKTtFrS5oZldReHOsT2jvqYO4CDW8RyDfA1YEe9XGH71pkqTsssIvBk5y0z22uBte32S/oKcGiLXW/vpH5JzwSOBpbVm66UdKLtf53ufV0lM0kvBlYCN9n+cjfvjYjR1c+bM9o+pd0+SXdLWmp7h6SlwM4WxX4H+KbtXfV7/gX4FWDaZDZtN1PSdQ2vXwf8JXAg8A5JTxi4i4i5aRYnANYBr61fvxb4QosyPwReImmBpAmqwf8Zu5kzjZlNNLxeDZxq+13AbwCvbvemxj71lybvmymGiBiyyd2THS89uhg4VdJtwKn1OpJWSLqkLvMZ4HvAd4EbgRtt//NMFc/UzRyrzwMZA2T7HgDbP5O0u92bGvvU6yeWz/2z8SIKN+nZmc20/WPg5BbbNwPn1q/3AP+127pnSmYHAddTTY9a0qG275K0qN4WEQUo/nIm20e22TVJNUgXEQUoPpm1Y/sh4Pt9jiUihiT3M4uIIkx2cZ7ZqEoyiwgm9+wZdgg9SzKLiPk7ZhYRZUkyi4gizNZ5ZoOUZBYRaZlFRBm6uWvGqEoyi4jMZkZEGfp5C6BhSTKLiHQzI6IMmQCIiCI4p2ZERAkmd8/9CQDNlavlJa2ub/o4MkYtpsQzvVGLB0YzprlqLj1qrpvHWc2WUYsp8Uxv1OKB0YxpTppLySwioq0ks4gowlxKZqM4rjBqMSWe6Y1aPDCaMc1Jc2YCICJiOnOpZRYR0VaSWUQUYU4kM0mrJG2RtFXSmiHH8lFJOyXdNMw4pkg6QtLXJN0q6WZJfzgCMe0n6TpJN9YxvWvYMQFIGpf0bUlfHIFY7pD0XUk3SNo87HhKMPJjZpLGgf9H9Sj3bcAm4FW2bxlSPCcCu4CP2T52GDE0xbMUWGr7W5IOpHpo8yuG9f3UMQk4wPYuSRPA1cAf2v7msGKq43oTsAJ4su2XDzmWO4AVtu8dZhwlmQsts5XAVtu3234UuBw4Y1jB2P5X4CfDOn4z2ztsf6t+/SBwK3D4kGOy7V316kS9DPWvpqRlwGnAJcOMIwZnLiSzw4E7G9a3MeRf1lEl6UjgecC1w41kb5fuBmAncKXtYcf0F8BbgVG5otrAlyVdLylXAfTBXEhmarFttPvGQyBpEfBZ4HzbDww7Htt7bB8HLANWShpal1zSy4Gdtq8fVgwtnGD7+cDLgP9eD19ED+ZCMtsGHNGwvgzYPqRYRlI9LvVZ4B9sf27Y8TSyfR+wEVg1xDBOAE6vx6kuB14q6e+HGA+2t9c/dwKfpxpOiR7MhWS2CThK0jMkLQTOAtYNOaaRUQ+2/x1wq+33DzseAElPlbS4fr0/cArwb8OKx/b/tL3M9pFU/3++avv3hxWPpAPqyRokHQD8BjASs+Nz2cgnM9u7gfOAK6gGtz9l++ZhxSPpE8A1wHJJ2ySdM6xYaicA/4WqtXFDvfzWkGNaCnxN0neo/hhdaXvop0OMkEOAqyXdCFwHrLf9pSHHNOeN/KkZERGdGPmWWUREJ5LMIqIISWYRUYQks4goQpJZRBQhySwiipBkFhFF+P8S+2zBLXy8/QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc9123cc8d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.heatmap(covariances, vmin=-1, vmax=1, cmap='coolwarm', square=True) \\\n",
    "    .set(title='$cov(Z_i, Z_j)$ of the GHZ state')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prepare GHZ as MPO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "psi = mp.chain([mup] * n_sites)\n",
    "rho = mp.mps_to_mpo(psi)\n",
    "# Apply Hadamard at site 0 (the first site)\n",
    "rho = mp.partialdot(mh, rho, 0)\n",
    "rho = mp.partialdot(rho, mh.adj(), 0)\n",
    "\n",
    "for pos in range(n_sites - 1):\n",
    "    # Apply CNOT at sites (pos, pos + 1)\n",
    "    rho = mp.partialdot(mcnot, rho, pos)\n",
    "    rho = mp.partialdot(rho, mcnot.adj(), pos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#rho.to_array_global().reshape([ldim**n_sites] * 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function computes the same covariance as `cov_ab()`, but for a mixed state provided as MPO."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cov_ab_mpo(rho, a, b, pos1, pos2):\n",
    "    # Todo: Verify that imaginary parts are small\n",
    "    if pos1 == pos2:\n",
    "        return np.nan\n",
    "    ept_a = mp.trace(mp.partialdot(rho, a, pos1)).real\n",
    "    ept_b = mp.trace(mp.partialdot(rho, b, pos2)).real\n",
    "    ept_ab = mp.trace(mp.partialdot(a, mp.partialdot(b, rho, pos2), pos1)).real\n",
    "    cov = ept_ab - ept_a * ept_b\n",
    "    return cov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "covariances_mpo = [\n",
    "    [\n",
    "        cov_ab_mpo(rho, mpz, mpz, pos1, pos2)\n",
    "        for pos2 in range(n_sites)\n",
    "    ]\n",
    "    for pos1 in range(n_sites)\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATMAAAEKCAYAAAB+LbI7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAF3NJREFUeJzt3X20XXV95/H3597cCEOQ1EYgEBSnYgqigxJjLRapQBuLBbsGV7F2lC4wM8tBy6iVjM5Cpa1SO2p11bam1AXaVnyuqUlFRGOHESFBQXlohogoMYGAykMEhOR+5o+9b3o4nHPvOTnn3HPu735ea+11z977d377e05yv/f3sB9km4iIuW5s2AFERPRDkllEFCHJLCKKkGQWEUVIMouIIiSZRUQRkswioghJZl2Q9B5J5w/p2NdJenYP718u6duSHpT0xg7fc4ekU/b1mBGzKcmsQ5KeCrwG+Ei9fqGkXU3LI5Is6Xc7rLObOv43cFEPH+GtwEbbB9r+UItYBp64JJ0l6VpJP5O0s379eklqF4OksyVdXb9+dYvva1f9fV3YZSxdfd4k9tGXZNa5s4ENth8GsH2R7UVTC3AocAOwAfhcJxV2Wcc64NclLd3H+J8O3LyP7+2ZpDcDHwT+nOpzHgL8N+AEYGEnddj+h8bvq/7OzgfuBv52MJHHnGG7yAV4FXA9cD/wPeAkQMAFwA+A+4BPAQfVi4ElDe8/huqX5Mn1+leB329zrP2BjcBVwH77GO+MdQBXAq+dpo6j6zruo0pcpzfEvgd4BNgFPKvpfR8HJoGH6/1vrbffAbwF+E79PX5yKjbgMOCzwD3A94E3ThPXQcDPgP88w3dwB3BK07azgavblH8e8CBw0jR1XgD8qC63BTh5ms+7pv6/8iBwC/A7M3w/HX8HWQa/DD2AgXwoeHP9y/x8qtbnc4AjgT+pf9mXAvsBXwDeV7/nTuDEhjo+D7ypYf0e4AUtjrUQ+BfgG8CifYy3ozqADwHvb7NvAtgKvK2u76X1L+Xyev9G4Nxp6m6VSO4Arqt/aZ8C3ErVmhqj+kNxYX2s/wjcDvxmm7pXAbuBBTN8Dx0nM2BxnXgumKa+5fW/62H1+pHAL01zrFfWn3UM+F2qBLy0Vfluv4Msg1+K62bWY1vvAH7P9rdsT9r+LtVf1TfU23fYfgT4DLCifutNVC0bJK2kSoQfbqh6MVVyaDzWAuByqi7Ty2zvatj3eknP6iDetnXU+8+VdGy9+mAdRyu/AiwCLrb9qO2vAl+kaqH24kO2t9v+CfDPwHHAC4CnuuomP2r7dqpu3llt6lgC3Gt7d8Pn+oak+yQ9LOnEhrL/VG+/T9J9wF81V1aPsV1G9W/23mli3wM8CThG0oTtO2x/r11h25+uP+uk7U8CtwEr2xTv9juIAVsw7AAG4BTgu7ZvbNr+a/X27Q3blgA76tc3UXUtAd4DvNP2zxvK/hQ4cGpF0hhwKdVf/5fYvr/xYLaf8EvYbKY66nouaVg9kKoL2cphwJ22Jxu2/QA4fKY4ZnBXw+uH6uM8HTisTjZTxoH/06aOHwNLJC2YSmi2fxVA0jYeP3b7CttfmVqRdDZwblN9FwDHAsfbbnvbF9tb69nndwLPlnQFVWt7e6vykl4DvImqBQfVH4clbarv9juIASuuZUbVHWr1C/9UqnGfRmcAV9evb6L6C34KVTf0Y01lvwM0trT+mqo1dIrte5sPJqmT/9TT1lHXs7Fh9WigOUlP2Q4cUSfIKU+jGi/qRDf3groT+L7txQ3LgbZ/q035a4CfU33fPZF0EvB24Ezb7RL7Xrb/0faLqZKPgT+b2tVU79OpWlbnAb9oezHV/wm1Kk/330EMWInJ7NvAiyX9J1WOknQ0sAl4kaRfkrRI0kVUXbuP1u+bapm9G3ib7T1N9W4AXgIg6f3Ay4CTbe9oKoekJcDOhvVLJV3aVGbaOuoyB1J3bSU9CTieahKglWupxnjeKmmi/qX/baoubCfuphr36cR1wAOSLpC0v6RxScdKekGrwnXSeRfwV5LOrL//MUnHAQd0eEzqmdzLgfNtf7uD8sslvbT+7h6hGmqY+ndt/rwHUCWse+r3/gFV64825bv6DmIWDHvQbhAL8EdUfzl3USWp59fb30zVrbwP+DT14G69b3+q/+jXtKlzCbCNajLBwKN1/Y3L3VR/IF4KvKPhvVcBr2tYP3amOupyLwLeU79+JfC5GT73s4GvU7VA987G1fs2Mv0EwBnAD+vv5i31tjt4/KD3O4G/r18fBnyCqhv6U+CbNA2otzjGq6mSwENUSeNaYDWwsNXx6m1nU08AUA22u8V3tgv4mxbHe259vAeBn1CNIR42zef907rcvcD76+/y3GnKd/0dZBncovofJTog6d3ATtt/MUO584Ef2P68pIVUXcPn2n6sy+OtBnbZ/kdJ1wLn2L5pX+OPKFmJEwADY/ttHRZ9DtXMH7YfpZ4l3QfPAdbW9bxwH+uImBfSMhth9eD/qd226CLmoxInAIog6Z+orqVMIovoQFpmEVGEtMwiogizMQEwUk2/9RPLhx1CRN+d9tgWzVyqvfUTyzv+Pe31WIOSlllEFCGnZkQEmhjJxlZXkswigvH9x4cdQs+SzCKCsQVpmUVEAdLNjIgipGUWEUVIyywiipCWWUQUYXzh3D/lNMksItBYWmYRUQCNp2UWEQUYG0/LLCIKkG5mRBQhEwARUQSNzYNkJumXqR6zdTjVvcm2A+ts3zrg2CJilpTQzZw2HUu6gOqhq6J6/uCm+vUnJK0ZfHgRMRvGxtXxMqpmapmdAzy7+aEa9dO4bwYubvWm+nmPqwE+8pGPsHr16j6EGhGDUkLLbKZkNkn11OYfNG1fWu9ryfZa6uc9MmK3zY6IJ+r3mJmkVcAHgXHgEtsXN+1/GnAZsLgus8b2hl6OOVMyOx+4StJtwJ31tqcBzwTO6+XAETE6xif6l8wkjQMfBk4FtgGbJK2zfUtDsf8FfMr2X0s6BtgAHNnLcadNZra/JOlZwEqqCQBNBWd7Ty8HjojR0edu5kpgq+3bASRdTjWJ2JjMDDy5fn0Q1cRiT2aczbQ9CXyz1wNFxOjqppvZOCZeW1sPLU05nH/vyUHVAHphUzXvBL4s6Q3AAcAp3cTbSs4zi4iuWmZNY+Itq2v1tqb1VwGX2n6fpBcBH5d0bN142idJZhHR727mNuCIhvVlPLEbeQ6wCsD2NZL2A5YAO/f1oHP/tN+I6NnYgvGOlw5sAo6S9AxJC4GzgHVNZX4InAwg6WhgP+CeXj5DWmYR0deTYW3vlnQecAXVaRcftX2zpIuAzbbXAW8G/lbS/6Dqgp5tu6fTuJLMIqLvJ83W54xtaNp2YcPrW4AT+nnMJLOImB8XmkdE+ebD5UwRMQ8kmUVEETqcpRxpSWYRkTGziCiE0s2MiAJkzCwiipBuZgfWTywf9CEiokdpmc1Rpz22Zdgh7JVkH6Mgs5kRUYS0zCKiDBkzi4gSKKdmREQJMpsZEUVQJgAiogSZAIiIIkjpZkZECdIyi4gSZAIgIoqQMbOIKILGM5sZESVINzMiSpArACKiDGmZRUQJMgEQEWUo4KTZff4Ekv6gn4FExPBofLzjZVT1ko7f1W6HpNWSNkva/KXJ+3o4RETMijF1voyoabuZkr7TbhdwSLv32V4LrAVYP7Hc+xxdRMyK+XAFwCHAbwI/bdou4BsDiSgiZt88ODXji8Ai2zc075C0cSARRcTsK71lZvucafb9Xv/DiYhhGOWB/U7l1IyIKOLUjCSziBjpWcpOzf10HBE9k8Y6XjqrT6skbZG0VdKaacqdKcmSVvT6GdIyi4i+tswkjQMfBk4FtgGbJK2zfUtTuQOBNwLX9uO4aZlFRDVm1ukys5XAVtu3234UuBw4o0W5PwbeCzzSj4+QZBYRMD7e+TKzw4E7G9a31dv2kvQ84AjbX+zXR0g3MyK6ms2UtBpY3bBpbX3Vz94iLd6290ogVQNvHwDO7i7I6SWZRURXY2aNlyu2sQ04omF9GbC9Yf1A4FhgY31TyEOBdZJOt72540CaJJlFRL/PM9sEHCXpGcCPgLOAvSfZ274fWLL30NXVRG/pJZFBxswiAqprMztdZmB7N3AecAVwK/Ap2zdLukjS6YP6CGmZRUTfr820vQHY0LTtwjZlT+rHMZPMIgLGcm3mnLR+YvmwQ4gYLaXfNSNmx2mPbRl2CHsl0c9T8+B+ZhExH+SuGRFRhLTMIqIEzs0ZI6II6WZGRBGSzCKiBM6YWUQUIS2ziChCWmYRUYLMZkZEGdLNjIgSOMksIoqQMbOIKEFaZhFRhtzPLCJKUMJJszO2LSX9sqSTJS1q2r5qcGFFxKzq70OAh2LayCS9EfgC8AbgJkmNTyV+9yADi4jZY9TxMqpm6ma+Djje9i5JRwKfkXSk7Q/S+kGfwOMfEnre2MGsGlvcp3AjYhDmwwTAuO1dALbvkHQSVUJ7OtMks8aHhK6fWO525SJiRBSQzGb6BHdJOm5qpU5sL6d6gOdzBhlYRMyeybHxjpdRNVMyew1wV+MG27ttvwY4cWBRRcTs6uNDgIdl2m6m7W3T7Pu//Q8nIoZhPoyZRcQ8MMqzlJ1KMouItMwiohAjPBbWqSSziGBSoztL2akks4hINzMiypAJgIgoQlpmEVGEEm4BlGQWEZkAiIgylDBmNvc7yhHRM2us46UTklZJ2iJpq6Q1LfY/SdIn6/3X1rcY60mSWUT09eaMksaBDwMvA44BXiXpmKZi5wA/tf1M4APAn/X6GZLMIqLfLbOVwFbbt9t+FLgcOKOpzBnAZfXrzwAnS73NQmTMbASsn1g+7BBinutmzKzxTtK1tfUNWaccDtzZsL4NeGFTNXvL2N4t6X7gF4F7uwj7cZLM4glOe2zLsEN4nCT7wZvsopPWeCfpNlplxuY7TndSpitJZhGB+zvitA04omF9GbC9TZltkhYABwE/6eWgGTOLiH4/nWkTcJSkZ0haCJwFrGsqsw54bf36TOCrttMyi4je9PM8s3oM7DzgCmAc+KjtmyVdBGy2vQ74O+DjkrZStcjO6vW4SWYR0feTZm1vADY0bbuw4fUjwCv7ecwks4go4gqAJLOIYNJzf/g8ySwi0jKLiDIkmUVEEewks4gowGRaZhFRgkwAREQRMmYWEUXImFlEFGFetMwkrQRse1N9t8hVwL/VlytERAGKb5lJegfVrW8XSLqS6gZrG4E1kp5n+08HH2JEDNrksAPog5laZmcCxwFPAu4Cltl+QNKfA9cCLZNZ450ozxs7mFVji/sXcUT03XyYzdxtew/wkKTv2X4AwPbDktom88Y7Ua6fWN7TPYoiYvCK72YCj0r6D7YfAo6f2ijpIMpomUYE82MC4ETbPwew3Zi8Jvj3u0RGxBw3WUD/adpkNpXIWmy/lx6eohIRo2U+tMwiYh6YD2NmETEP7Ekyi4gSpGUWEUXo7SFvoyHJLCIyARARZSj+1IyImB8mJ9Myi4gC5LbZEVGETABERBFyakZEFCETABFRhHQzI6IIuZwpirR+YvmwQ4hZlpZZxCw57bEtww5hrxKTfZJZRBRhMt3MiChBCS2zuf9Ilojo2Z7JzpdeSHqKpCsl3Vb//IVpyj5Z0o8k/WUndSeZRQS2Ol56tAa4yvZRwFX1ejt/DHy904qTzCICu/OlR2cAl9WvLwNe0aqQpOOBQ4Avd1pxkllEMOnOlx4dYnsHQP3z4OYCksaA9wF/1E3FmQCIiK5aXJJWA6sbNq2tH/w9tf8rwKEt3vr2Dg/xemCD7Tulzru1SWYR0VUyqxPX2mn2n9Jun6S7JS21vUPSUmBni2IvAn5N0uuBRcBCSbtsTze+lmQWEb3PUnZhHdUDxC+uf36huYDtV0+9lnQ2sGKmRAYZM4sIYHKy86VHFwOnSroNOLVeR9IKSZf0UnFaZhExayfN2v4xcHKL7ZuBc1tsvxS4tJO6k8wioogrAJLMIiI3Z4yIMrirptloXpTe9QSApI8NIpCIGJ49ezpfRtW0LTNJ65o3Ab8uaTGA7dMHFVhEzJ75MGa2DLgFuAQwVTJbQXWpQVuNZwifN3Ywq8YW9x5pRAxMCWNmM3UzVwDXU12GcL/tjcDDtr9uu+3V7LbX2l5he0USWcTom8ULzQdm2paZ7UngA5I+Xf+8e6b3RMTc466aZqM5AdBRYrK9DXilpNOABwYbUkTMtlm8nGlgumpl2V4PrB9QLBExJJMFDJqlyxgRIz0W1qkks4hIMouIMkwWkM2SzCICz7cJgIgo0549aZlFRAG6u9B8NCWZRUQRlzMlmUVEl1cAjKYks4jIqRkRUYZcARARRZjMbGZElCAnzUZEEXJqRsQsWT+xfNghFC1jZhHz1GmPbRl2CH1VQMMsySwiYE8Bd2dMMouInDQbEWVIMouIIhSQy5LMIiIts4goRM4zi4giZDYzIoqQbmZEFCHJLCKKUMKF5mPDDiAihs+T7njphaSnSLpS0m31z19oU+69km6WdKukD0nSTHUnmUUEtjteerQGuMr2UcBV9frjSPpV4ATgucCxwAuAl8xUcbqZEcGe3bM2m3kGcFL9+jJgI3BBUxkD+wELAQETwN0zVZyWWUR01TKTtFrS5oZldReHOsT2jvqYO4CDW8RyDfA1YEe9XGH71pkqTsssIvBk5y0z22uBte32S/oKcGiLXW/vpH5JzwSOBpbVm66UdKLtf53ufV0lM0kvBlYCN9n+cjfvjYjR1c+bM9o+pd0+SXdLWmp7h6SlwM4WxX4H+KbtXfV7/gX4FWDaZDZtN1PSdQ2vXwf8JXAg8A5JTxi4i4i5aRYnANYBr61fvxb4QosyPwReImmBpAmqwf8Zu5kzjZlNNLxeDZxq+13AbwCvbvemxj71lybvmymGiBiyyd2THS89uhg4VdJtwKn1OpJWSLqkLvMZ4HvAd4EbgRtt//NMFc/UzRyrzwMZA2T7HgDbP5O0u92bGvvU6yeWz/2z8SIKN+nZmc20/WPg5BbbNwPn1q/3AP+127pnSmYHAddTTY9a0qG275K0qN4WEQUo/nIm20e22TVJNUgXEQUoPpm1Y/sh4Pt9jiUihiT3M4uIIkx2cZ7ZqEoyiwgm9+wZdgg9SzKLiPk7ZhYRZUkyi4gizNZ5ZoOUZBYRaZlFRBm6uWvGqEoyi4jMZkZEGfp5C6BhSTKLiHQzI6IMmQCIiCI4p2ZERAkmd8/9CQDNlavlJa2ub/o4MkYtpsQzvVGLB0YzprlqLj1qrpvHWc2WUYsp8Uxv1OKB0YxpTppLySwioq0ks4gowlxKZqM4rjBqMSWe6Y1aPDCaMc1Jc2YCICJiOnOpZRYR0VaSWUQUYU4kM0mrJG2RtFXSmiHH8lFJOyXdNMw4pkg6QtLXJN0q6WZJfzgCMe0n6TpJN9YxvWvYMQFIGpf0bUlfHIFY7pD0XUk3SNo87HhKMPJjZpLGgf9H9Sj3bcAm4FW2bxlSPCcCu4CP2T52GDE0xbMUWGr7W5IOpHpo8yuG9f3UMQk4wPYuSRPA1cAf2v7msGKq43oTsAJ4su2XDzmWO4AVtu8dZhwlmQsts5XAVtu3234UuBw4Y1jB2P5X4CfDOn4z2ztsf6t+/SBwK3D4kGOy7V316kS9DPWvpqRlwGnAJcOMIwZnLiSzw4E7G9a3MeRf1lEl6UjgecC1w41kb5fuBmAncKXtYcf0F8BbgVG5otrAlyVdLylXAfTBXEhmarFttPvGQyBpEfBZ4HzbDww7Htt7bB8HLANWShpal1zSy4Gdtq8fVgwtnGD7+cDLgP9eD19ED+ZCMtsGHNGwvgzYPqRYRlI9LvVZ4B9sf27Y8TSyfR+wEVg1xDBOAE6vx6kuB14q6e+HGA+2t9c/dwKfpxpOiR7MhWS2CThK0jMkLQTOAtYNOaaRUQ+2/x1wq+33DzseAElPlbS4fr0/cArwb8OKx/b/tL3M9pFU/3++avv3hxWPpAPqyRokHQD8BjASs+Nz2cgnM9u7gfOAK6gGtz9l++ZhxSPpE8A1wHJJ2ySdM6xYaicA/4WqtXFDvfzWkGNaCnxN0neo/hhdaXvop0OMkEOAqyXdCFwHrLf9pSHHNOeN/KkZERGdGPmWWUREJ5LMIqIISWYRUYQks4goQpJZRBQhySwiipBkFhFF+P8S+2zBLXy8/QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc90ff453c8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.heatmap(covariances_mpo, vmin=-1, vmax=1, cmap='coolwarm', square=True) \\\n",
    "    .set(title='$cov(Z_i, Z_j)$ of the GHZ state')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prepare GHZ as MPO, with errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n_sites = 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Prepare GHZ as MPO with some error probability eta\n",
    "\n",
    "eta = 0.1\n",
    "\n",
    "psi = mp.chain([mup] * n_sites)\n",
    "rho = mp.mps_to_mpo(psi)\n",
    "# Apply Hadamard at site 0 (the first site)\n",
    "rho = mp.partialdot(mh, rho, 0)\n",
    "rho = mp.partialdot(rho, mh.adj(), 0)\n",
    "\n",
    "for pos in range(n_sites - 1):\n",
    "    rho_last = rho\n",
    "    # Apply CNOT at sites (pos, pos + 1)\n",
    "    rho = mp.partialdot(mcnot, rho, pos)\n",
    "    rho = mp.partialdot(rho, mcnot.adj(), pos)\n",
    "    # Apply CNOT with probability (1 - eta) only.\n",
    "    rho = eta * rho_last + (1 - eta) * rho\n",
    "    rho.compress(method='svd', relerr=1e-10)\n",
    "    #print(rho.ranks)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "covariances_mpo_eta = [\n",
    "    [\n",
    "        cov_ab_mpo(rho, mpz, mpz, pos1, pos2)\n",
    "        for pos2 in range(n_sites)\n",
    "    ]\n",
    "    for pos1 in range(n_sites)\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATYAAAENCAYAAACFAUnHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3XmcHFW99/HPt2cSsrGGHcIDjyyibGoMCCpIAAERxBW8KiIaH1FAvC744Etc7lXE3Yv4GCEEvBBEFonsEQxclC0sYQsIIksgEvY9kkz/nj+qJtQ0PTPndFV3V3f/3q9XvdJd/ZtTp3smp09Vnd85MjOcc66bVNpdAeecK5o3bM65ruMNm3Ou63jD5pzrOt6wOee6jjdszrmu4w2bc67reMPmAJD0fUlfbNOxb5D0xnYc23Unb9gcktYBPgH8On3+TUkv1GzLJJmkjwSWGVPGj4DvFPuuXC+TZx44SV8BtjSzzwzz+iTgT8CTwPvMbHkDxxi2DEnjgEeBN5rZkgbegnNDeI+tw0iaIuk8SY9LelLSien+rSXNl/SMpDsl7Z/uP0bSOTVl/FzSLzK79gGuGuZ444ELgReBDzTYqI1YhpktA24C9oot27l6vGHrIJL6SBqIB4FNgY2AsySNAf4IXA6sCxwBnCFpK2AOsK+k1TJlfBg4M1P0tsA9dY43FjgPGAsckDZAsXUOLWMRsH1s+c7V4w1bZ5kGbAh8xcxeNLNlZnYNsBMwCTjezF4xsytJGsCDzexB4GbgfWkZuwMvmdl1mXLXAJ7PHkhSP3AWsB6wj5m9kHntcElbjlbZkcpIX/+0pG3Sp8+n9XAut/52V8BFmQI8aGYravZvCDxsZtXMvgdJenSQ9M4OBk4HPsrQ3hrA08Cqg08kVYDZwFbArmb2bDbYzE4araKjlZGWc3Lm6arAM6OV61wI77F1loeBTdKeUNajwJS0MRm0CfBI+vj3wG6SNgYO5LUN221Atgf2K5Je4B5m9kRtJST9T0BdRywjLWd+5unWwMKAcp0blTdsneUGYAlwvKSJksZJ2gW4nuTC/FcljZG0G/BektNAzOxxYD5wKvAPM1tUU+7FwK4Akn5CcjNher07lJLWBpZmns+WNLsmZsQy0phVSU9/Ja0CvAWYF/xJODcCb9g6iJkNkDRYmwMPAYuBj5jZK8D+JI3JE8BJwCfM7O7Mj58J7MFre2uQnKLuK2lb4GhgfeDOmjFoj6U9wu1IeniDpgB/GXySXjMbrQyAbYA70sf7A/PN7NEGPxrnhvBxbA4ASd8DlprZz0aJ+yLJdb7z0zueC4HtYoeBSJoBvGBmZ0q6HjjMzO4Y7eecC+E3DxwAZvZ/A0O3JRlaQtpT3LrBQ24LzEzL2bHBMpyry3tsri3SGwd7NjLg17nR+DU213KS/kByTc0bNdcU3mNzznUd77E557qON2zOua7T6ruiQee9f5361uAC15iyWnDsxHVWHT0oNX5yeLljVpsUHNs/aWJwbN/ECUFxlfFhcQCMHx8eO3ZceOyYscGhNmaViNjwcqt9Y8LL7Qv707eIMqvqCz9+pTmxa237dgUHD+OiMVsFX596z/J7ch+vGXy4h3NuCI0pZVsVJdepqKS9Jd0j6T5JxxRVKedc+1T6FbyVVcM9tnRer18Ce5Kk9twoaa6Z3VVU5ZxzracxnX/pPc+p6DTgPjO7H0DSWcABgDdsznWwvvGd37DleQcbkUyjM2gxr87/tZKkGZIWSFowc+bMHIdzzrVCT5+KAvXe1WvuppjZTNKcwHqvO+fKpRtuHuRp2BaTTFkzaGOSCQ+dcx2szD2xUHkathuBLSRtRjJT60Ek00475zqY+nq4YTOzFZK+AFwG9AGzzOzOwmrmnGuLvh6/K4qZXUwyrXSQmIyCUM89+vzoQSlVwn9hMbEx4uoQ+M0ZUWYltMwmivkMTOH1jfmNVSPKDT5+eIIA1eroMe2ivh5v2Jpp5wU3BsVdt+O0JtfEud5S6YJT0byZB7MkLZXkUzo71yVUUfBWVnn7nLOBvQuoh3OuJCp9Ct7KKu81tqslbVpMVZxzZdDTd0Wdc92p0h9xF6Skmn77I5tSdcHjjzf7cM65nHr+VDRENqXqr1Pf6ilVzpVcmW8KhOr8ASvOuUKpUgnegsobZd5GSZtI+rOkWyTdJmnfvO8h73CPOcC1wFaSFks6LG+FnHPtVeRwj8y8jfsAbwAOlvSGmrBvAGeb2ZtIUjNPyvse8t4VPThvBZxz5VJwSlXIvI0GDC4ysjoFTKbR0ruioQuvPPfo803JKHjxiReCY/vGht8ZqowJ/xgrY8IXB1FfYB1C4wBiUpQU8Qcek08Uk1YWUV+LqK/UhLWaY44fcx3LWnvFKCrtT5oBzMjsmpleVx9Ub97GHWuK+RZwuaQjgInAHjH1racrhnvsdP0NQXELdn1bk2viXOeLaXRr5lusW1y9H6t5fjAw28x+LOltwG8lbWNmDWfUNvxVIGlKesFvkaQ7JR3VaFnOufIoOKUqZN7Gw4CzAczsWmAcsHae95Cnj7sC+Hcz2xrYCfh8nYuCzrkOU3DDtnLeRkljSW4OzK2JeQiYDiBpa5KGLdeg1zzzsS0BlqSPn5e0iOR82hdzca6DFTll13DzNkr6DrDAzOYC/w78RtLRJKepnzSzXGNeC7nGluaLvgm4vojynHPtU/REk/XmbTSzb2Ye3wXsUuQxc78DSZOAc4EvmtlzdV5fmVJ19gOP5D2cc67Jih6g2w65emySxpA0ameY2Xn1YrJ3Te46cLqnVDlXct2QUpVnJXgBpwCLzOwnxVXJOddO3dCw5elL7gJ8HNhd0q3pljvHyznXXj19Kmpm11B/8J1zroN1Q4+tpZkHE9dZNSgu5pvgxSdeaEpGwbJnlwXH9q8S/jH2jY1IvwqMVcTEgMFpWoDFpH8tj0ipiqiDVsSki0X0IEJTtSLKjEvpihnd3+KUqpgUvZLqipSqGFOvujYo7rZ9d2tuRZwrqZ7usUkaB1wNrJKWc46ZHVdUxZxz7VHma2eh8vTY/gXsbmYvpMM+rpF0iZldV1DdnHNt0NM9tjTlYXAeoDHp5uPUnOtwvd5jG5wd8yZgc+CXZuYpVc51uEp/5zdsud6BmQ2Y2Q4kU5FMk7RNbUw2perMRQ/kOZxzrhUqlfCtpAqpmZk9A8ynzqrwZjbTzKaa2dSPbr1pEYdzzjWRpOCtrPJMNLmOpDXSx+NJpvO9u6iKOefao6czD4ANgNPS62wVklVmLiymWs65dun1u6K3kczB5pzrJiXuiYVqaebB+Mlhq1TFdHFjVpNa9uyypmQUvPLiK8GxYyb8Kzi2b9zYoLiYla8qY8PKBLDl4Ss5qT+8DqxYER7bF/4nqupAcKyFxlpE+lfM2iMRE8RWLPx9FaHiKVXdbbuL5wfF3f2hvZpbEedaqZdPRQel19gWAI+Y2X75q+Sca6cy3xQIVUSP7ShgEa+u5Oyc62DdcPMgV9MsaWPgPcDJxVTHOdd2qoRvJZW3x/Yz4KtA2ERrzrnS6+kem6T9gKVmdtMocStTqk6/9W+NHs451yLq6wveyirvmgf7S3oAOItk7YP/rg3KplR9YoctcxzOOdcSvZwramZfN7ONzWxTkmXrrzSzjxVWM+dcW6ii4K2sfBybc26oEt8UCFVIw2Zm80lm9xjRmNUmFXG4ISpjwt9CzKIrr7z4SlMG3q5YFj6av39ZWEbDy8ueYpU1w+7fxGQTxGQI2IomZSkMRGQpVCKyBCpho/n1r5ex/sBsjSYt/BK88ExRStwTC+U9toK8/veXB8X941P7F37s0EbNxQtu1LpImW8KhPKGzTk3RM9nHqR3RJ8HBoAVZja1iEo559qoxBNIhiqix/YuM3uigHKcc2XQ6z0251wX6oIeW96m2YDLJd0kaUa9gGzmwewb7sp5OOdcs/X61OAAu5jZo5LWBeZJutvMrs4GmNlMYCbA098/3Ncdda7suuCuaN7l9x5N/10KnA9MK6JSzrn2kSrBW1h52lvSPZLuk3TMCHEflGSSct+EzJMEP1HSqoOPgb2AO/JWyDnXZhWFb6NIJ6L9JbAP8AbgYElvqBO3KnAkUMii63l6bOsB10haCNwAXGRmlxZRKedcGxU7H9s04D4zu9/MXiGZMOOAOnHfBU4AlhXxFvKsUnU/sH3UwSZNDIqLuSgZs5BJ39jwtxuz6MqKZcubklEw8EpYOlE1MA6gujw8VhHpV30Ri8QQkX5Ff8SfaF/EoiehC69ELLoSXCbNW/ilEMXeFd0IeDjzfDGw49DD6U3AFDO7UNKXizioD/dog81mzQ2KW3zEh5tcE+fqiOhYpKMhsiMiZqY3DFeG1PmxlS21kgt1PwU+GVfJkXnD5pwbKuKuaHbUwzAWA1MyzzcGHs08XxXYBpivpKe4PjBX0v5mtiC4IjXyrnmwhqRzJN0taZGkt+UpzzlXAsVeY7sR2ELSZpLGkszduPKUxcyeNbO1zWzTdG7H64BcjRrk77H9HLjUzD6YVnpCzvKcc+1W4DU2M1sh6QvAZUAfMMvM7pT0HWCBmYVdl4nUcMMmaTXgnaTnxukdj/Al0Z1z5VRwRoGZXQxcXLPvm8PE7lbEMfO8g/8NPA6cKukWSSen49mGyKZUnfqXhTkO55xrCSl8K6k8DVs/8GbgV2b2JuBF4DWjirOLuRy6S9ToEOdcO/T1h28lladhWwwsNrPBkcLnkDR0zrlO1ss9NjP7J/CwpK3SXdMBn77DuU7nK8FzBHBGekf0fuDQ/FVyzrVViXtioXI1bGZ2KxCcid83MWw0SMx6hTELT1QiUqr6xoWnCIWuJgVJmlTRGQXViNSnmNi+gfAUpZjVr6JWqapGpElVI1KaQt9b4GpWAFjUwNaIciPSr4pQ4nnWQpX36p8DYOP/OnvUmKVf/0QLauJ6hfV6j80514Uqnd8s5JmPbStJt2a25yR9scjKOedaz6TgrazyTFt0D7ADrJxM7hGSWXSdc52sxHc7QxXV55wO/N3MHiyoPOdcu5S4JxaqqKb5IGBOvReyKVWzrrqpoMM555qmUgnfSip3jy0dw7Y/8PV6r2fna3pp1nG+SpVzJWeVzl+lqohT0X2Am83ssQLKcs61mfk1NgAOZpjTUOdcB+qCa2y5GjZJE4A9gc8WUx3nXLv1fI/NzF4CJofGV8YHTrAbc1EyIqVK/RHpVxGrX1XHR6Q0Ra0otbzwrAKLSJOKWtFqbET6VcxqThH1jUo9Ck1piikzJqWrEh4blX5VhF7vsbnyWPf7pwfFPfUf3rl2o+j1Hptzrvt0w13RvKtUHS3pTkl3SJojaVxRFXPOtYeh4K2s8uSKbgQcCUw1s21IVqA5qKiKOefaw1QJ3soq76loPzBe0nKSpfceHSXeOVd2JW6wQuWZGvwR4EfAQ8AS4Fkzu7w2LptSdcoV19e+7JwrmW6Y3SPPqeiawAHAZsCGwERJH6uNy65Sddj0HRuvqXOuJbrhVDRPzfYA/mFmj5vZcuA8YOdiquWcaxer9AVvZZXnGttDwE5p9sHLJFMXLSikVs65tinz3c5QeSaavF7SOcDNwArgFtJZPJxznavMp5ih8qZUHQccF/wD48cHhVUiVqmKSf+IW9EqfJWqmBWaYtKUQleUsoGB5mQURKQIUY1I+4lJk4qpQ0y5/YHlNiNNq+xKfFMglGce9KC1vvHroLjnfuJLWPQiK2z+2fbJm3lwVJp1cKcv5OJcd6hW+oK3ssoz3GMb4DPANGB7YD9JWxRVMedce/R0ShWwNXCdmb1kZiuAq4ADi6mWc65den0c2x3AOyVNTod87AtMKaZazrl26enMAzNbBPwAmAdcCiwkGfYxxJCUqsv/2nBFnXOt0eunopjZKWb2ZjN7J/AUcG+dmFdTqvbyxATnyq7XT0WRtG767ybA+/FFXZzreFX1BW8hJO0t6R5J90k6ps7rq0j6Xfr69ZI2zfse8o5jO1fSZGA58HkzezpvhZxz7VXkKaakPuCXJIs+LQZulDTXzO7KhB0GPG1mm0s6iOQS10fyHDdv5sE78vy8c658Cj7FnAbcZ2b3A0g6i2RWoGzDdgDwrfTxOcCJkmQ5VrFpbebB2OJnDq9E/BIsYuWpmDQpVkSs5hRRbl9gilBMmhbValMyCmJWv7KY1Zxi/rabEBtz/JjVt2JStRRTbgFiemySZgAzMrtmmlk2Z3wj4OHM88VA7fxlK2PMbIWkZ0lWv3siotpDeEqVG9FqX/pZUNwLJ73m0onrUDHDONJGbKTJL+oVVvttERITZdTujqRZkpZKuiOzby1J8yTdm/67Zp5KOOfKw0zBW4DFDB3fujGvXUJgZYykfmB1klEWDQs5j5sN7F2z7xjgCjPbArgife6c6wJV+oK3ADcCW0jaTNJYkgWf5tbEzAUOSR9/ELgyz/U1CGjYzOxqXtt6HgCclj4+DXhfnko458qjyAG6abrlF4DLgEXA2WZ2p6TvSNo/DTsFmCzpPuBLFNBRavQa23pmtgTAzJYMjmdzznW+ojMKzOxi4OKafd/MPF4GfKjIYzZ96PCQlKpLrm724ZxzOXVDSlWjPbbHJG2Q9tY2AJYOF5i9a7Ls4pldMsWoc90r8KZAqTXaY8te7DsEuKCY6jjn2q1KJXgrq1F7bJLmALsBa0taTLLGwfHA2ZIOI1mtqtDzY+dc+5T5FDPUqA2bmR08zEvTC66Lc64EuuFUtLWZB2PCV34KFjHvupZHxPZHpF+tiEiTasLqVxobs+pTRIrQwEBzMgqiUp+ak6YUtfpVE8SlarX20nS1F3pszoWadPjxQXEvzf52k2vi8uiGU9FGU6o+lK5MVZU0tblVdM61UsEpVW3RaErVHSQTS/rANOe6TNUqwVtZhdw8uLp2Rst0vQNU4sUcnHON6YZTUb/G5pwbosynmKFam1J10fxmH845l1M1YiurpvfYhqRUzZvtKVXOlVw39Nj8VNQ5N0RPXGMbJqXqKeC/gHWAiyTdambvbmZFnXOtUea7naHypFSdX3BdnHMlEJGcUlotPRW1MasExakS8Y0RE9sXnlIVtfJURPoVEelXoeVGrY7UpNWkMGtORkHM/7KY2NDPrEkrapW57eiJU1HnmmHCJ48Lils25wdNromr1Q03DxpNqfqhpLsl3SbpfElrNLeazrlWMQvfyqrRlKp5wDZmth3wN+DrBdfLOdcmA1YJ3sqqoVWqzOzydPUZgOtI1gp0znWBbuixFXGN7VPA7wooxzlXAt1w8yBXX1LSscAK4IwRYl5NqfrjFXkO55xrgaqFb2XVcI9N0iHAfsD0kVZtzqZUvTx/Tok/CuccdMdd0YYaNkl7A18DdjWzl4qtknOuncp87SxUoylVXwdWAealc7JdZ2b/p4n1dM61yEAv9NiGSak6pZGDWeBiLhYxgWXMZJdaEZF50BfRmR0Iz1KgP6LcaliWgCKyCZo1kj5uIRVrzsDbqIVfAt9b0z6DiN9DpbUTBPXsqahzrTTu4K8FxS0796dNrklvKPNNgVDesDnnhuiGa2yNplR9N02nulXS5ZI2bG41nXOtYih4K6tGU6p+aGbbmdkOwIXAN4uumHOuPXpiHNswq1Q9l3k6kXLPwuKcixAzW1VZNZx5IOk/JT0M/Bsj9NiymQezLpjX6OGccy1SNQVvZdVww2Zmx5rZFJJ0qi+MEDfTzKaa2dRPHbBno4dzzrVINyTBFzHvyJnABwooxzlXAq1q2CStJWmepHvTf9ccIXY1SY9IOjGk7IYaNklbZJ7uD9zdSDnOufJp4c2DY4ArzGwL4Ir0+XC+C1wVWnCjKVX7StqKZM3UBwFPp3KuS1SrLbt2dgBJ2wJwGjCfJAd9CElvAdYDLgWmhhTc0pSqal/Y4iQx3UhTRHRErALTmQCoxKRqNSH9qUmLuTQvnSiu3LZmFMS8r6hyy3uBqoXDONYzsyUAZrZE0rq1AZIqwI+BjwPTQwv2zAPXVcZ94OiguGUX/qrJNelcMW2upBnAjMyumelUZYOv/wlYv86PHht4iMOBi83s4Zi8cG/YnHNDRHXUM/MtDvP6HsO9JukxSRukvbUNgKV1wt4GvEPS4cAkYKykF8xspOtxjaVUZV77siSTtPZo5TjnOkMLbx7MBQ5JHx8CXFAbYGb/ZmabmNmmwJeB00dr1KDxlCokTQH2BB4KKMM51yFaOI7teGBPSfeStCXHA0iaKunkPAU3lFKV+inwVeq0ss65zhVzfykPM3uSOjcEzGwB8Ok6+2eTdLRG1eg4tv2BR8xsYUDsqylVf7iskcM551qoGzIPom8eSJpAckdjr5D47MXFF66bW+KPwjkH5Z61I1Qjd0VfB2wGLExvv24M3Cxpmpn9s8jKOedab4RF5+ooZyJ8dMNmZrcDKwfSSXoAmGpmTxRYL+dcm5T5FDNUyHCPOcC1wFaSFks6rPnVcs61y8BA+FZWjaZUZV/fNPRgFrjyUzVqlarlwbFElGsRKVWqxKQpRaxOFPqXE/MV29+ElZwgbnbCZqwmlWprRkHT0q9aO/Njr15jc64rjNvvc6PGLLusobTojtYNp6LesDnnhrCoLls5bx40ukrVt9JJ325Nt32bW03nXKt0w2IuDadUAT81sx3S7eJiq+Wca5eeGKA7QkqVc64LDQyUuMUKlGfNgy+kiybPGmWu8pUpVaeef0mOwznnWqEbemyNNmy/IslA2AFYQjLDZV3ZVaoOPXCfBg/nnGuVqlnwVlYN3RU1s8cGH0v6Dclq8M65LtDiYXNN0VDDNjjrZfr0QOA1k1A65zpTXK5oOTW6StVuknYADHgA+GwT6+ica6GYJJKyaukqVRa4SlXTxKxoZRErTzVr5afQVK2oFKWIlK5mfXNHpV81aaUsejOrIMRAmQeoBfLMA+dGMe7dYfM+vHzlb5tck9aIyzwop4YXc5F0hKR7JN0p6YTmVdE510rdMNwjpMc2GzgROH1wh6R3kazivJ2Z/aveQqfOuc5U7YIeW6OZB58Djjezf6Ux9dYDdM51oG64K9roAN0tSRYxvV7SVZLeWmSlnHPtMzBgwVtZNdqw9QNrAjsBXwHO1jDrzw9JqTrvogYP55xrFTML3sqq0buii4HzLHlnN0iqAmsDj9cGZlepeu7meeX9JJxzQHdcY2u0x/YHYHcASVsCYwFfzMW5LtATd0WHyTyYBcxKh4C8AhxiZe6XOueCdcM4tjyLuXys4Lo450qgzLN2hGpp5kFVYWlKlYhsppg0KYuIVUx6TlSqVkw6UeAH0aQUJWvSNA9RqVolWPmpKRkFJZ5Co7qivHUL5SlVzhVo/O4fD4p76aqzmlyTxnXBmWjQNbZZwH7AUjPbJt33O2CrNGQN4Bkz26FptXTOtUxPXGOjTkqVmX1k8LGkHwPPFl4z51xbdMN9wFyLuaSDcj9MOvTDOdf5umEcW95rbO8AHjOze4uojHOu/bqhx5ZnlSqAg4E5IwVkU6pmn+tLIzhXdtUV1eCtrBrusUnqB94PvGWkuGxK1TO3XNn5XwXOdbluGMeWp8e2B3C3mS0uqjLOufazqgVveUhaS9I8Sfem/9Zdn1jSCemEtosk/WK4CTeyQmbQnQNcC2wlabGkwXmSD2KU01DnXOdp4ewexwBXmNkWwBXp8yEk7QzsAmwHbAO8Fdh1tIIbTqkys0+O9rPOuc7TwruiB5DkoQOcBswHvlYTY8A4kok2BIwBHmMULc08sMBcqZgMIVVG7ZW+Gjt6D/ZVMalHTUrVCv1GVKX9aUcxaVJR/22ataJVM1i11BkFoaoDEb93aQYwI7NrZnpdPcR6g+sTm9mSeksMmNm1kv4MLCFp2E40s0WjFewpVc61yYRdDwqKe/Gv5zW5JkNZxJdJ9uZgPZL+BKxf56VjQ8qXtDmwNbBxumuepHea2dUj/VyjKVU7AP+PpIu4AjjczG4IqahzrtyKPBU1sz2Ge03SY5I2SHtrGwD11k45ELjOzF5If+YSkpm7R2zYQs6hZgN71+w7Afh2mh/6zfS5c64LtPDmwVzgkPTxIcAFdWIeAnaV1C9pDMmNg1FPRUdt2NIu31O1u4HV0serA4+OVo5zrjO0argHcDywp6R7gT3T50iaKunkNOYc4O/A7cBCYKGZ/XG0ghu9xvZF4DJJPyJpHHdusBznXMm0anYPM3sSmF5n/wLg0+njAeCzsWU3OkD3c8DRZjYFOBo4ZbjAISlV54za0Drn2mxgYCB4K6tGe2yHAEelj38PnDxcYPauydMLr+r8XA3nulw3zMfWaI/tUV4d/bs74LN7ONclemJd0WFWqfoM8PM0EX4ZQwfoOec6WDVmUHRJ5VmlasRZPZxznakbTkVLmVIVWWh4aERsxSIujDYpVSs0palZpwTNSP+KrkMJ0sWapdUZBaGatTpZK3lKlXMdYOLO72/Zsco8gWSokGmLZklamq76Prhve0nXSrpd0h8lrTZSGc65zlG1avBWVo2mVJ0MHGNm2wLnA18puF7OuTZpYeZB0zSaUrUVryahzgM+UHC9nHNtYtVq8FZWjY5juwPYP338IWBKMdVxzrVbT/TYhvEp4POSbgJWBV4ZLjCbUnXaOXMbPJxzrlV6NqXKzO4G9gKQtCXwnhFiV6ZUPXX7NeVt4p1zQNxEk2XVUMMmaV0zWyqpAnyDZNJJ51wXKPMpZqhGU6omSfp8GnIecGrTauica6meGKA7QkrVzwuui3OuBLqhxxaVyd+MDZjhseGx7T5+N8e2+/jNjO21rf0VgAUeGx7b7uN3c2y7j9/M2F7bGh3u4ZxzpeUNm3Ou65ShYQtdNdpjy3H8bo5t9/GbGdtTlJ6rO+dc1yhDj8055wrlDZtzrut4w+ac6zotnxpc0uuBA4CNACNZym+umS0qoNyNgOvN7IXM/r3N7NKa2GmAmdmNkt5AMpHm3WZ28SjHON3MPhFYn7cD04A7zOzyzP4dgUVm9pyk8cAxwJuBu4DvmdmzmdgjgfPN7OGA440FDgIeNbM/SfoosDOwCJhpZstr4l8HHEgy5dQKkiUU52SP71ynaunNA0lfAw4GzgIWp7s3JvkPeZaZHR9YzqFmdmrm+ZHA50ljhomOAAAF3UlEQVT+E+8AHGVmF6Sv3Wxmb87EHgfsQ9KozwN2BOYDewCXmdl/pnG1cywJeBdwJYCZ7T/kRekGM5uWPv5MWp/zSWZB+ePge5N0J7C9ma2QNBN4CTgHmJ7uf3+mzGeBF4G/A3OA35vZ48N8Jmek72kC8AwwiSSPdzrJ7/mQms/rvcBVwL7ArcDTJA3d4WY2v94xOtXgpA0FlznZzJ4sskxXoFaOBgb+Boyps38scG9EOQ/VPL8dmJQ+3hRYQNK4AdxSJ7aPpAF4Dlgt3T8euC0TdzPw3yQTAOya/rskfbxrnTrdknl8I7BO+ngicHvmtUXZY9SUcWttmSSXC/YCTgEeBy4FDgFWrYm9Lf23H3gM6EufK/u+sp9B+ngCMD99vEmdz2t14HjgbuDJdFuU7lsj4nd2Sc3z1YDvA78FPlrz2kk1z9cHfgX8EpgMfCt9D2cDG9TErlWzTQYeANYE1qqJ3bvmfZ4C3AacCayXee14YO308VTgfuA+4MHav4X07+YbwOsCPpOpwJ/Tv7MpJF+0z6Z/P2+qiZ0EfAe4M415HLgO+GQr/w93ytbqa2xVYMM6+zdIX1tJ0m3DbLcD69X8fJ+lp59m9gBJI7SPpJ+Q/MfOWmFmA2b2EvB3M3su/bmXa+owFbgJOBZ41pJezMtmdpWZXVXnPVQkrSlpMkkP6fG03BdJTvUG3SHp0PTxQklT0/e7JTDkdDH5caua2eVmdhjJZ3cSyanz/XWOP5Zk4s8JJP9RAVYBxtSpb3/m9VXTgz1UJ/Zskt7cbmY22cwmk/RcnwZ+nw2U9OZhtreQ9KSzTiX53ZwLHCTpXEmrpK/tVBM7m+RU/WGShuBlkjkA/4fXTpn1BMnvbXBbQHKJ4ub0cdb3Mo9/TPLF9V6ShuXXmdfeY2ZPpI9/CHzEzDYH9kx/LmtNYA3gz5JukHS0pHp/85D8Lk8ALgL+CvzazFYnuTxxUk3sGSS/83cD3wZ+AXwceJek7+GGamUrSvIf8j7gEpLBhTNJeiD3kfn2TGMfI/nP8L9qtk1JriNlY68EdqjZ1w+cDgzU7L8emJA+rmT2r05NDyrdvzHJf+ATqekp1sQ9QPKH94/03/Xt1W/aW2uOM5vk9PJ6ksbsfpLTwu1ryrxlhOONr3l+dFrOg8CRwBXAb0h6NsfVxB5F0jOZSdITOzTdvw5wdU3sPSPU4Z6a5wPp7+LPdbaXa2Jre6fHAn8h6WHV9mSzveHa3nptOV9O/6a2zez7xzD1v3mEcrK/s7uB/vTxdTVxt49Q5jtIGqh/pp/BjJrYkd5Xbc95Yc3zGwf/hkmuD7f0/3LZt9YfMPlF7ESyAMwH08d9deJOAd4+TBln1jzfeLAhqRO7S83zVYaJWzv7n6HO6+8hubgf+34nAJvV2b8qsD3wFjKnPTUxW0Yea0Ngw/TxGunnO22Y2Demr79+lDIvB77K0FOz9YCvAX+qib0D2GKYch6ueb6IzBdLuu8QklOtB2v2L8w8/o+a126vc6zBL6OfpJ/z/cPUaTHwJeDfSb4UlHkte1niiPRz2J3kNPhnwDtJek6/rSmz3pdjH8mX+qk1+68luczwIZIvpPel+3elJsGdpEf39vTxe0muBw++NuyXT69uba+Ab+XeSE6tfkDSa3kq3Ral+9asif0gsNUw5byv5vkJwB514vam5norybWlSXViNwfOGaHu7yW5DvXPYV4/rmYbvC66PnB6TexuwO9IrnveDlwMzCDtyWXizor4bLcHLiM5g3k9yRyHz6SN+841sdsBN6SvX0P6pUfSyz6y3X8nZdvaXgHfOncjPYUtcyzJTaFtmlWHTvgMenHzXFHXMEkPmdkmvRzb7uPHxvaKlg/QdZ1F0m3DvUTN3elujW338WNjnTdsbnTrkQwxeLpmv0guaPdCbLuPHxvb87xhc6O5kOTC/a21L0ia3yOx7T5+bGzP82tszrmu47N7OOe6jjdszrmu4w2bc67reMPmnOs63rA557rO/wdfRZUXsgTyxgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc9100aee10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.heatmap(covariances_mpo_eta, vmin=-1, vmax=1, cmap='coolwarm', square=True) \\\n",
    "    .set(title='$\\mathrm{cov}(Z_i, Z_j)$')\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}