mario-ambrosino/von_neumann_entropy_pure_state.ipynb

## von_neumann_entropy_pure_state.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reply to \"Entanglement lost\" Question in QML Course\n",
    "\n",
    "[Check the question here](https://courses.edx.org/courses/course-v1:University_of_TorontoX+UTQML101x+2T2019/courseware/5833ff13fd1f4507a4016b4a60fa75bb/8022840dfd7a438981ea6aa15501a0e5/?child=first)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First let's import the necessary code from the Jupyter Notebook, in this we define the density mixed matrix obtained from a pure state density matrix added with a maximally_mixed_state "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def mixed_state(pure_state, visibility):\n",
    "    density_matrix = pure_state * pure_state.T.conj()\n",
    "    maximally_mixed_state = np.eye(4)/2**2\n",
    "    return visibility*density_matrix + (1-visibility)*maximally_mixed_state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, let's define the Von Neumann Entropy. [I borrowed this code forme Computer Science Stack Exchange](https://cs.stackexchange.com/questions/56261/computing-von-neumann-entropy-efficiently), it works very well (avoids the infinity problem that the log_2 implies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "def von_neumann_entropy(density_matrix, cutoff=10):\n",
    "    \n",
    "    # assign as numpy matrix the input\n",
    "    x = np.mat(density_matrix)\n",
    "    # define a identity matrix with the same dimension \n",
    "    one = np.identity(x.shape[0])\n",
    "    # performs the polynomial expansion\n",
    "    base = one - x\n",
    "    power = base*base\n",
    "    result = np.trace(base)\n",
    "    for k in range(2, cutoff):\n",
    "        result -= np.trace(power) / (k*k - k)\n",
    "        power = power.dot(base)\n",
    "\n",
    "    # Twiddly hacky magic.\n",
    "    a = cutoff\n",
    "    for k in range(3):\n",
    "        d = (a+1) / (4*a*(a-1))\n",
    "        result -= np.trace(power) * d\n",
    "        power = power.dot(power)\n",
    "        result -= np.trace(power) * d\n",
    "        a *= 2\n",
    "    result -= np.trace(power) / (a-1) * 0.75\n",
    "\n",
    "    return result / math.log(2)  # convert from nats to bits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another way to test the purity of a quantum state in Von Neumann representation (density operator representation) is the [**purity test**](https://en.wikipedia.org/wiki/Purity_(quantum_mechanics)):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def purity_test (density_matrix):\n",
    "    return np.trace(np.matmul(density_matrix,density_matrix))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, let's define the state as a Bell's state:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "# State definition\n",
    "zero_ket = [[0],[1]]\n",
    "one_ket = [[1],[0]]\n",
    "\n",
    "state_00 = np.kron(zero_ket,zero_ket)\n",
    "state_11 = np.kron(one_ket,one_ket)\n",
    "\n",
    "ϕ = (state_00 + state_11)\n",
    "ϕ = ϕ / np.linalg.norm(ϕ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we define a data structure which contains the Von Neumann entropy values and purity test of the the mixed state in function of the visibility parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "visibilities = np.linspace(0, 1, 100)\n",
    "\n",
    "entropy = []\n",
    "ptest = []\n",
    "for v in visibilities:\n",
    "    m_state = mixed_state(ϕ, v)\n",
    "    entropy.append(von_neumann_entropy(m_state))\n",
    "    ptest.append(purity_test(m_state))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we plot this function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f697833c780>"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEKCAYAAAAMzhLIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl4lNXZ+PHvyWTfVwIBsgCBsIdFQARURAUVcaOiiKKiL1q1LrVqf6+Vt7VvfVsttdVqUakbiooLiFpXUJB933dCEhLIRvZ1Muf3xzMkgQyTIclsyf25rueaedY5ea5k7jxnuY/SWiOEEEKci4+7CyCEEMKzSaAQQghhlwQKIYQQdkmgEEIIYZcECiGEEHZJoBBCCGGXBAohhBB2SaAQQghhlwQKIYQQdvm6uwDnKzY2VicnJ7u7GEII4VU2b95coLWOa825XhcokpOT2bRpk7uLIYQQXkUpday150rVkxBCCLskUAghhLBLAoUQQgi7vK6NQgjhverq6sjOzqa6utrdRemwAgMD6dGjB35+fu12TQkUQgiXyc7OJiwsjOTkZJRS7i5Oh6O1prCwkOzsbFJSUtrtulL1JIRwmerqamJiYiRIOIlSipiYmHZ/YpNAIYRwKQkSzuWM+yuBQgghhF3SRiGE6FRMJhODBw9uWJ8xYwZPPvnkOY9fuXIl/v7+jB071hXF80gSKIQQnUpQUBDbtm1z+PiVK1cSGhpqM1CYzWZ8fTv+16hUPQkhBEZ6oGeeeYbhw4czePBg9u3bR0ZGBq+++irz588nPT2dVatWMXv2bB599FEuvfRSnnjiCYqKirjuuusYMmQIY8aMYceOHQDMmzePWbNmMXHiRFJTU3nttdcAmDVrFkuXLm343JkzZ7Js2TK3/MyO6vihUAjhkZKf/MJp18547upz7quqqiI9Pb1h/amnnuLmm28GIDY2li1btvDPf/6T559/ntdff525c+cSGhrKr3/9awDeeOMNDhw4wHfffYfJZOLBBx9k2LBhfPbZZ/zwww/cfvvtDU8sO3bsYN26dVRUVDBs2DCuvvpq5syZw/z585k2bRolJSWsWbOGt956y2n3oj1IoBBCdCr2qp5uuOEGAEaMGMEnn3xyzmtMnz4dk8kEwOrVq/n4448BmDhxIoWFhZSUlAAwbdo0goKCCAoK4tJLL2XDhg1cd911/PKXvyQvL49PPvmEG2+80eOrrzy7dEII4UIBAQGA0eBtNpvPeVxISEjDe611s/2nu6ie3VX19PqsWbNYtGgRixcvZuHChW0ut7NJoBBCuIW96iFPEhYWRmlp6Tn3T5gwgUWLFvH000+zcuVKYmNjCQ8PB2Dp0qU89dRTVFRUsHLlSp577jkAZs+ezahRo+jatSsDBw50yc/RFhIohBCdytltFJMnT274Ardl6tSp3HTTTSxdupR//OMfzfbPmzePO++8kyFDhhAcHHxGe8OoUaO4+uqryczM5OmnnyYhIQGA+Ph4+vfvz3XXXdeOP5nzSKAQQnQq9fX1NrdnZGQ0vB85ciQrV64EoG/fvg09mQDGjx9/xnnR0dFn9GJqqm/fvixYsKDZ9srKSg4ePMgtt9xynqV3D+keK4QQLvTdd9+RlpbGgw8+SEREhLuL4xB5ohBCCCeYN2+eze2TJk0iMzPTtYVpI3miEEIIYZcECiGEEHZJoBBCCGGXBAohhBB2SaAQQnQqJpOJ9PR0Bg0axPTp06msrDyv81999VXefvttAN58801ycnIcPvf6668nPT2dPn36EBERQXp6Ounp6axZs+a8yvDDDz+wbt268zqnLSRQCCE6ldO5nnbt2oW/vz+vvvqqw+eazWbmzp3L7bffDpx/oPj000/Ztm0br7/+OuPHj2fbtm1s27btvOe6kEAhhBAuMn78eA4dOkRGRgaDBg1q2P788883dG+95JJL+O1vf8vFF1/Miy++yLx583j++edZsmQJmzZtYubMmaSnp/PFF19w/fXXN1zj22+/bUgy6IiNGzdy8cUXM2LECKZMmcLJkycBmD9/PgMGDGDo0KHcdtttHD58mNdff52//OUvrXoaaQ0ZRyGEcI95ThxsNq+kxUPMZjNfffUVkydPbvHY4uJifvzxR+PS1gBy00038dJLL/H8888zcuRItNY89thj5OfnExcXx7///W/uvPNOh4pbU1PDr371K5YtW0ZsbGxD7qgFCxbw5z//mWPHjuHv709xcTGRkZHMmTOH2NhYHn74YYeu31byRCGE6FRO53oaOXIkiYmJ3H333S2ec3q+CnuUUsyaNYt3332X4uJi1q5dy5QpUxwq0969e9m9ezeTJk0iPT2d5557jqysLAAGDhzIbbfdxqJFi/Dz83Poeu1NniiEEJ2KrfkofH19sVgsDevV1dVn7G+aVtyeO++8k6lTpxIYGMj06dMdnmdCa82QIUNYtWpVs31ff/01P/74I0uXLuXZZ59l165dDl2zPUmgEEK4hwPVQ64SHx9PXl4ehYWFhIaGsnz5coeqpMLCwigrK2tYT0hIICEhgWeffZZvv/3W4c8fMGAAx48fZ8OGDYwaNYra2loOHjxIWloa2dnZTJw4kXHjxrFo0SIqKyubfa6zSdWTEKLT8/Pz43e/+x2jR4/mmmuuIS0tzaHzZs+ezdy5c0lPT6eqqgow5sDu2bMnAwYMcPjzAwICWLJkCY8++ihDhw5l2LBhrF+/HrPZzK233sqQIUMYPnw4TzzxBGFhYUybNo0PP/yQYcOGuaQxW9manaldLqxUT+BtoCtgARZorV886xgFvAhcBVQCs7XWW+xdd+TIkXrTpk1OKbMQwrn27t1L//793V0Mp3rggQcYNmyYQ20fzmLrPiulNmutR7bmes6sejIDj2mttyilwoDNSqlvtdZ7mhwzBUi1LqOBV6yvQgjhdUaMGEFISAgvvPCCu4vSrpwWKLTWuUCu9X2ZUmov0B1oGiimAW9r47FmnVIqUinVzXquTXtzSxn9v9/hoxQ+SmHyUfj6KHxNCj+TD/6+PvibfAj0MxHkZyLI30Sgn4kQfxMhAb6EBfoSGuBLeJAfYYG+hAf6ERnsR2SwP+GBvs3muBVCCEdt3rzZ3UVwCpc0ZiulkoFhwPqzdnUHspqsZ1u3nTNQmC2ak6U17VxCg8lHERXsR3SIP9Eh/sSEBhAXGkBcmPHaJTyA+PBAuoYHEhnsJ0FFiFbQWsvfjhM5oznB6YFCKRUKfAw8rLU+e4ZyW78tzX5KpdS9wL0A/l37tHsZT6u3aArKaykor23x2ABfH7pFBNItIoiEyCC6RwXRIyqIHpFB9IwOpltEIL4m6SsgRFOBgYEUFhYSExMjwcIJtNYUFhYSGBjYrtd1WmM2gFLKD1gOfK21/quN/f8CVmqt37eu7wcusVf1NHTYCP31yp+xaE29RWPRmrp6jdlioc6sqa2vp8ZsoabOQlVdPVW19VTW1VNZY6aixkx5TT3lNXWUVpkpra6jpKqO4so6iitrqai1PZdua/j6KBIig0iKCSYpJpjkmBCSY0LoFRdCz+hg/CSIiE6orq6O7OzsZuMURPsJDAykR48ezQbneWRjtrVH0xvAXltBwmoZ8IBSajFGI3aJvSAB4GdSdI1o32h5Wq3ZwqnKWgrLaymqqKWgvIaC8hryy2rIK6vhZGm1damhvMZs91pmiyazqJLMokpWHTxzn6+PIjEmmD5xofTpEkpqfCipXcLo0yWUQD+TU342ITyBn58fKSkp7i6GOE/OrHq6CJgF7FRKnR4G+VsgEUBr/SrwJUbX2EMY3WMdS4ziJP6+PsSHBxIf3nIgKquuI7ekmpziKnKKqzleXMnxU1Vknaoiq6iSvLJzt6OYLZoj+RUcya/gmz0nG7YrBUnRwfSNDyOtWzj9uxqvSdHB+PjIY7oQwj2cWvXkDN4yjqKqtp6sU5UcK6wko6CCjMIKjhYYweFE6fk9dgf7m0jrGsbAhAgGJoQzqHsEfePD8PeV6ishhGM8suqpswvyN9E3Poy+8WHN9lXUmDmSX8Hh/HIO5pVxKK+cg3nlZBRUYLERtytr69mSWcyWzOKGbf4mH9K6hTG4ewRDe0YytEckfbqEYpInDyFEO5MnCg9SXVfP4fxy9p8oY9+JMvbmlrI3t4yCcse6A4f4mxjcI4JhiVEM6xnJsMQo4sICnFxqIYQ3kCeKDiLQz2StXjozT39eWTW7c0rZk1PK7pwSdh4vIauoqtn5FbX1rDtSxLojRQ3bkmKCGZEYxYjkKC5IjqZPXKi0dwghzos8UXipUxW17MopYUd2CduzitmeXezQQMTIYD9GJkUxKiWa0SkxDEwIl/EeQnQCbXmikEDRgeSWVLE1s5htWcVsPnaKncdLqDVb7J4TGuDLyOQoLuwVw4W9YxiYECHtHEJ0QBIohE015np255Sy5dgpNmYUsTHjFEUV9kedhwX6cmGvGC7qE8tFfWLoHRcqI2iF6AAkUAiHaK05nF/Bxowi1h8pZP3RInJL7HfV7RoeyPjUWMalxjKuTywxodI4LoQ3kkAhWkVrTVZRFWuPFLD2cCFrjxTabedQCoZ0j2BC3zgu7htHes9Iad8QwktIoBDt4vQTx5rDBaw+WMDaI4WUVZ87VUlEkB8T+sZxaT8jcMjThhCeSwKFcApzvYUdx0tYdaCAVQfz2ZpVTL2tEYEYTxvDE6OYmNaFy/p3oV98mLRtCOFBJFAIlyipqmPNoQJ+PJDPyv35dlOR9IwO4vL+XZk0oAujkqOlikoIN5NAIVxOa83e3DJW7M9jxb48tmSespl+BCAq2I/L+sdzxYB4JvSNkwy5QriBBArhdkUVtfx4II/v9ubx4/78c6ZhD/Y3cWm/LkwZ3JVL+3UhJECSAwjhChIohEepNVtYd6SQb/ec5Ns9J89ZRRXg68Ml/eK4ekgCl6VJ0BDCmSRQCI9lsWh2HC/h690n+HrXCY4UVNg8LsDXh4lpXbhmSAIT07oQ5C/VU0K0JwkUwitorTlwspyvduXyn10n2HeizOZxwf4mLh8Qz7T0BManxsm0sUK0AwkUwisdyivny525fLkz95xBIyrYj6sGd+P6Yd0ZkRQlXW6FaCUJFMLrHcor4/PtuXy+I4cj+barp3pGB3FdeneuG9ad3nGhLi6hEN5NAoXoMLTW7MktZdn2HD7flkPOOXJRpfeM5MYRPZg6pBuRwf4uLqUQ3kcCheiQLBbNxowiPtuWwxc7cii1kU7E3+TD5QPiuWlED8anxsrAPiHOQQKF6PBqzPWs2JfPJ1uyWbE/j7r65r+38eEB3DC8B78Y2ZOU2BA3lFIIzyWBQnQqRRW1LN+Rw8ebs9meXWLzmFEp0cy4oCdXDe4mI8GFQAKF6MQOnCxjyeZsPtlynILy5inSwwN9uWF4D2aM6kla13A3lFAIzyCBQnR6dfUWVu7P54ONWazYn2czy+2IpChuHZXI1UPkKUN0PhIohGgir7SajzZn88HGLDKLKpvtjwjy46YRPZg5OpFe0s1WdBISKISwwWLRrD1SyHsbMvlm9wmbDeDjU2O5bUwSl6V1kR5TokOTQCFECwrKa/hoUzbvbThGVlFVs/3dI4OYOSaRGRckEh0i4zJExyOBQggHWSyaVYcKeGftMX7Yd7LZHBr+vj5cOzSBOy9KZmBChHsKKYQTSKAQohWyT1Xy3vpMFm/Moqiittn+USnR3HVRMpcP6IrJR3JMCe8mgUKINqiuq+eLHbm8tTaDHTbGZfSMDmL22BR+MbIHYYF+ri+gEO1AAoUQ7UBrzZbMYt5ck8FXO3Mxn1UvFRrgy4wLenLnuBS6Rwa5qZRCtI4ECiHa2YmSat5Zl8Gi9ZkUV9adsc/ko7h6cDfundCLQd2lHUN4BwkUQjhJVW09n2zN5o3VR22mP7+oTwz/NaE341NjZa4M4dEkUAjhZBaLZuWBPBb8dIR1R4qa7e/fLZz7LunNVYO6yngM4ZEkUAjhQruOl7DgpyMs35HTrHttYnQw907oxU0jekiaEOFRJFAI4QZZRZW8sfooizdmUl1nOWNfXFgA94xPYeboJEICfN1UQiEaSaAQwo2KKmp5c00Gb63JoKTqzIbvyGA/7hybwuyLkokIkq61wn0kUAjhASpqzLy/IZPXVh3hZOmZKc/DAny5Y2wyd41LkRQhwi0kUAjhQWrM9Xy8+Tiv/ni4WfbaYH8Tt1+YzD3jU4gJDXBTCUVnJIFCCA9krrewbHsOL6041KxrbbC/iVkXJnHv+F4SMIRLSKAQwoPVWzRf7crlH98fYv/JsjP2BfubuGNsMveO70WUVEkJJ5JAIYQXsFg0X+8+wYvfH2TfiTMDRoi/ibvGpTBnfC9p9BZOIYFCCC9isWi+2XOCv33XPGCEB/py74Re3HlRinSrFe1KAoUQXshi0Xy16wR/++4AB/PKz9gXE+LP/Zf2YeboRBm4J9qFBAohvFi9RbN8Rw7zvz1ARuGZvaQSIgJ5eFJfbhjeXVKDiDaRQCFEB2Cut/DJluO8+P1BjhefOV1rny6hPH5lP64YEC/JB0WrSKAQogOpMdezaF0mL684ROFZM+8NT4zkySn9GZUS7abSCW8lgUKIDqiixszC1UdZ8NMRymrMZ+yb1D+eJ6f0o0+XMDeVTngbCRRCdGBFFbW8vOIQ76w9Rm19Y/JBHwU3X5DII5en0iUs0I0lFN5AAoUQnUBWUSV//fYAn249fsb2YH8Tcy/uzT3jexHkLz2khG1tCRQOdaNQSl2jlJIuF0K4Uc/oYObfnM7yB8cxrk9sw/bK2nr++u0BLnl+BUs2Z2M5e5IMIdrI0S//GcBBpdSflVL9nVkgIYR9g7pH8O6c0bx11yj6xTe2UZwsreHXH21n2ss/szGj+Sx8QrSWw1VPSqlw4BbgTkAD/wbe11qX2T2xnUnVkxCN6i2aDzdl8cI3BygoPzO1+dWDu/HUVWn0iAp2U+mEJ3F61ROA1roU+BhYDHQDrge2KKUebM0HCyHazuSjuGVUIisfv4QHLu1DgG/jn/QXO3O57IUf+es3+6msNdu5ihD2OdpGMVUp9SnwA+AHjNJaTwGGAr92YvmEEA4IDfDl11f24/vHLmbq0ISG7TVmC3//4RCXvfAjS7cdx9s6rwjP4FDVk1LqbeB1rfVPNvZdprX+3hmFs0WqnoRo2eZjRcxbtoedx0vO2D4qJZr/uXYg/buFu6lkwl1c0j1WKdUVGIXRPrFRa32iNR/YVhIohHCMxaJZsjmbP3+9j4LyxhHePgpuG5PEY5f3IyJYUpp3Fq7oHns3sAG4AbgJWKeUuqs1HyiEcA0fH8UvLujJil9fwpxxKfj6GDmiLBreXnuMiS+s5KNNWdKdVrTI0aqn/cBYrXWhdT0GWKO17ufk8jUjTxRCtM6hvDLmLdvD6kMFZ2wfkRTFH6YNYkCCVEd1ZK7o9ZQNNO0GWwZkteYDhRDu0adLGO/cPYpXZg4nIaIx5cfmY6eY+tJq/rB8D+U10jtKNHc+jdmDgaUYbRTTMKqiDgBorf/qxDKeQZ4ohGi7ylozL/1wiNdWHaGuvvE7ID48gGemDmTKoK6SzryDccUTxWHgM4wgAUbAyAXCrIsQwosE+/vym8lp/OfhCYztHdOw/WRpDfcv2sLdb20iq6jSzhVEZ3JeSQGVUmGA1lqXt3iwk8gThRDtS2vNsu05/GH53jNGdwf5mXh4Uip3j0uR2fU6AFf0ehqklNoK7AJ2K6U2K6UGtuYDhRCeRSnFtPTufP/Yxcwak8TpGqequnr+9NU+rn3pZ3ZkF7u3kMKtHP03YQHwqNY6SWudBDwGvOa8YgkhXC0iyI8/XDeIj+8bS1rXxhrlPbmlXPfyz/z+8z1USGN3p+RooAjRWq84vaK1XgmEOKVEQgi3Gp4YxecPjuPJKWkNuaMsGhb+fJQr5v/ETwfy3VxC4WqOBoojSqmnlVLJ1uW/gaPOLJgQwn38TD7Mvbg33zwy4Yy5L44XV3H7wg38+qPtlFTWubGEwpUcDRR3AXHAJ9YlFiPduBCiA0uKCeGdu0fxwvShRDZJ97FkczaT5v/I17vdkslHuFiLvZ6UUibgOa31464pkn3S60kI98gvq2He57v5YkfuGdunDk1g3tQBxIQGuKlkwhFO7fWkta4HRrTm4kKIjiMuLICXbx3Ov2aNIC6sMSh8vj2Hy+f/xJc7c+2cLbyZo1VPW5VSy5RSs5RSN5xenFoyIYRHunJgV7575GJuHN6jYVtRRS33L9rCA+9t4VRFrZ2zhTdyNFBEA4XARGCqdbnGWYUSQni2iGA/XvjFUP49+wK6hjfmjVq+I5fL5//EN9J20aE4muvpIq31zy1tcwVpoxDCs5RU1fHs8j18tDn7jO03Du/BM9cOIDxQ5rzwBK7I9fQPB7cJITqZiCA//jJ9KAtnj6RLk7aLj7dkM3n+T6w5K6258D6+9nYqpS4ExgJxSqlHm+wKB0zOLJgQwrtMTIvnm0ei+N3S3SzbngNATkk1t76+nrvHpfD4lf0I9JOvDW/U0hOFPxCKEVDCmiylGDPdCSFEg8hgf/5+yzBeunXYGeMu3lh9lGtfWs2enFI3lk60lqNtFEla62MuKE+LpI1CCO+QV1rNEx/vYMX+xpQffibF41f2Y864Xvj4yHwXruSKNooApdQCpdQ3SqkfTi+t+UAhROfQJTyQhbMv4NnrBhHoZ3zV1NVr/vfLfcxauJ4TJdVuLqFwlKNPFNuBV4HNQP3p7Vrrzc4rmm3yRCGE9zmSX84jH2xje3ZJw7bIYD+eu2Ewkwd1c2PJOg9XPFGYtdavaK03aK03n15a84FCiM6nV1woS+4bywOX9mmY76K4so65727ht5/upKq23v4FhFs5Gig+V0rdr5TqppSKPr04tWRCiA7Fz+TDr6/sx+J7xpAQ0ThI7731mUx7eTX7T5S5sXTCHkcDxR3A48AajOqnzYDd+h+l1AGllEUpZbMiUhm2KaVqlVJVSqlbz6fgQgjvNLpXDF/9agJXD26scjpwspypL63m3XXHOJ/pmYVrOBQotNYpNpZeLZz2d+A2O/ufBroDAcADGLPoCSE6gYhgP166dRjP3TC4oaG71mzhvz/bxQPvbaW0Wua68CR2A4VS6jdN3k8/a9//2jtXa/0SkGnnkBnAYm14A/BXSg1tuchCiI5AKcWMUYksf3DcGVOvfrEzl2v+vlrm6W4v9WbYu7xNl2jpiWJGk/dPnbVvcps+GWKAPU3Wy4EhbbymEMLL9OkSxme/vIiZoxMbtmUWVXLTK2t5a02GVEW1VnEW/PBH+Nsg+GBmmy5lN4UHoM7x3tb6+bJ1vsXmgUq9A9wAEB0tbehCdDSBfib+eP1gxvaO5cmPd1BWY6a23sIzy3az4WgRf7pxsCQXdES9GQ59C5vfhIPfgLb5lXreWnqi0Od4b2v9fBUAA5qshwK7bBZC61la6xCtdUhKSkobP1YI4amuHtKN5Q+NY1D38IZtX+zMZeo/JP2HXcVZsOJP8OIQeH8GHPjPmUEipEubLt9SoBiqlCpVSpUBQ6zvT68PbtMnw4fADGvvp7uBWq319jZeUwjh5ZJiQlgydyyzxiQ1bDtWWMn1//yZjzZlubFkHqa+zmh7ePcm+Ntg+PE5KD1+5jEpF8P0t+DRPbav4SCHRma36sJKHQN6YASjeuAdjCSDaK1nKqUUsANIs+6/R2v9TkvXlZHZQnQen2/P4cmPd1DRZEDejAt6Mu/agZ03E23hYdj6Dmx7D8pPNt8fHAvDZsLwOyCmd8PmtozMbqmNotW01kkt7Ne0/alECNGBTR2aQP9u4dz37mYO5pUDsHhjFntyS3nlthF0jwxycwldpK4K9iwzAkTGKtvH9LrECA5p14Cvf7t+vNOeKJxFniiE6Hwqasw89cnOhnkuAKJD/PnHLcO4qE+sG0vmRFpD7jbY+i7s+AhqSpofE9oV0m+F4bdDtP3227Y8UUigEEJ4Ba01b63J4Nkv9mK2GN9bPgp+MzmN/5rQC6U6SNryyiLY8aERIE7ubL5f+UDqlTDiDuhzOZgcqxjyyKonIYRoT0opZl+UwsDuEdy/aAv5ZTVYNDz31T52HS/hzzcNIdjfS7/S6s1w+HsjOOz/Ciw2RqZHpcDwWTD0Vgh3bcZdeaIQQnidvNJq7l+0hU3HTjVsS+saxmu3j6RndLAbS3ae8vbB9vdg+wdQfqL5ft8gGDANht0GSReBj6Pp+ZqTqichRKdTa7bw++W7eXddY6agyGA//jlzOGN7e3C7RWUR7PoYtr8Px88xW0P3kUbPpUE3QmBEu3ysBAohRKf1/oZMfrd0F3X1xneZr4/imWsHnjEOw+3q6+Dgt0ZwOPAfqK9tfkxIFxg6A9JnQpe0di+CtFEIITqtW0Yl0jc+jLnvbia/rAazRfP0Z7vYl1vKvGsH4mdqfXVNm2gNOVth+2LYtQQqC5sfY/KHflOM4ND7Mocbpl1NniiEEB1CbkkV9769mZ3HG7uRXtgrhlduG05kcPuOK7CrONPotbTjAyg4YPuY7iNg6C1G1VKwa/LXSdWTEEIA1XX1PL5kB583GW+REhvCG3eMpFdcqPM+uOoU7FlqBIhjP9s+Jrw7DPmFESDi+jmvLOcgVU9CCIGRhfbvM9LpFx/K898Y/80fLajgupd/5pXbRrTv4Ly6ajj4tREcDn5ju93BLwT6T4X0WyB5PPh4Z9oRCRRCiA5FKcUDE1PpHRfKIx9uo7rOQmm1mTsWbuCP1w/i5gsSW77IuVjqjRQaOz+CPZ/bHi2tfIx0GkNmQP9rwD+k9Z/nISRQCCE6pCmDu9EjKpg5b2/kZKnRyP3ExzvJKKzk8Sv64ePj4EhureH4FqNBetfHthPxAXRLhyE3G+0OYfHt94N4AAkUQogOa3CPCJb+chx3vbmRPbnGfBavrDxMZlElL0wfaj8D7ck9jcHhVIbtY6KSYfAvYPB0iOvb7uX3FBIohBAdWteIQD6ceyEPvb+VH/blAfDFjlxOllTz2u0jiQpp0iOq8LARGHZ9Avl7bV8wJM54ahh0E/QYCR0lx5Qd0utJCNEpmOst/GH5Ht5ae6xhW6+4EN69oSscOGVlAAAW70lEQVQJ2V/C7k/hxA7bJweEG43Sg240JgPy0PEO9kivJyGEaIGvyYf/mTaIntHBvPnlT1zls46rStaT8NaRc5wQBP0mG08OfSaBX6BrC+xBJFAIITqHUxmwZylz9n7GnIAtto8x+RupuwfdAH0nQ4ATx154EQkUQoiOq/CwMRBuz1JjEiAb6rSJVZbBfKkv5JKrZnPNqPbPs+TtJFAIIToOrSFvL+z9HPYug5O7bB/n4wu9LuVk4hRuXx3L/lLjq3DJJ4fJrfbnngm9XFhozyeBQgjh3U6Pc9j3uREgCg/ZPs7HD3pPhIHXGYn4gqKIB94cWsUdCzdw4KQxJ/cfv9xLYUUtT0zu13FmzWsjCRRCCO9TXwfH1sC+5bDvCyg9bvs4U4DRED3gWqPNISiy2SHdIoL46L/GMuftjWzMMCZCevXHwxRX1vLH6wdjcnRgXgcmgUII4R1qK+DwD0Zg2P8VVBfbPs4vBPpeAWnXQN8rISCsxUtHBPvxzt2jeeC9LXy31xhrsXhjFiVVdfxtRjoBvt6Zo6m9yDgKIYTnKs8zJvrZ9yUcWQHmatvHBUVB3ylGbqXeE8EvqFUfV1dv4YklO/hka+MTyvjUWBbMGkmQv3cHCxlHIYToGLSG/H3GE8P+ryB7I3COf2bDu0Pa1caTQ9JF7TIIzs/kw/PThxIZ7M/Cn48CsOpgAbcvXM8bsy8gPNCvzZ/hjSRQCCHcy1wLmWtg/3/gwFfnzqsEENcf0q4ygkPCMKekz/DxUTx9TX/Cg3z523cHAdiYcYqZr63n7btGnZnyo5OQQCGEcL2KAmMO6QP/MdodakptH6d8IPFCo5dS2tUQ7Zpuq0opHp7Ul9AAX579wsj5tPN4CTcvWMuiOWOICwtwSTk8hQQKIYTzWSxGHqWD38CBr+H4Zs5ZpeQfarQzpF0NqVe4bKpQW+aM70Wwvy//77OdaA0HTpYzY8Fa3rtnDPHhnSelhwQKIYRzVBXDkZXGk8Ohb889jwNARKLRQ6nfFEgeB76e8x/7raMTCfY38eiH27BoOJxfwc3/MoJFQmTrGs29jQQKIUT70BpO7DSCwsHvIGs96Hrbxyof6DnaCA59J0Ncmken675uWHf8TD48tHgr9RZNRmElNy9Yy+J7L6R7JwgWEiiEEK1XUWh0Wz30PRz+3v5TQ1A09LnMCAy9J7q1Sqk1rh7SDV+T4oH3tlBXr8kqqmJGJwkWMo5CCOE4c63RZfXwD0ZgyNnGOdsawJgeNPVySL0Sug8HH+8eiwDww76TzH1nC7X1FgASo4NZfK/nV0O1ZRyFBAohxLlpDQUHjaeGwysgYxXUlp/7+KBo6H2pkaq7z2UQ2sV1ZXUhbwwWEiiEEO2n7CQc/cloiD6y4tx5lMBoa+hxAfS+DFInGU8QHeCpwRHf7z3J3Hc3U1dvfIcmxwTzwX9d6LG9oSRQCCFar7rESLB35Ec4+iPk7bF/fEQi9JlotDOkXGwz0V5ncXaw6NMllMX3jiE21HN6bZ0mgUII4bjaSqNH0tGfjCVn67l7J4ExX3TKBOh1iREcont5dA8lV/tm9wnuX7QFs8X4Lk3rGsb794zxuBHcEiiEEOdWV2U0QB9dBRmr4fgmqK899/E+fkbX1V4XG8EhYXi75FHqyJbvyOGh97dijRUM7h7BontGe1RuKEkKKIRoVFthBIaMnx0LDCjoNtR4aki5GJIuBP8QlxW3I7hmSAK1ZguPfbQdrY10H3Pe3MRbd43y+qyz4ORAoZSaDLwImIDXtdbPnbV/NvAX4HRr2Uta69edWSYhOpzqEshcbyTWy/gZcraAxWz/nLg0a2CYYGRe9bIxDZ7ohuE9qDFbeOqTnQBsyCji/kWb+deskfj7+ri5dG3jtEChlDIBLwOXA9nARqXUMq312S1lH2itH3BWOYTocEpzIXMtZK4zgsOJXdgdywAQkwop4yF5vJEio4N2W3W3W0YlUlZdx/9+uQ+AFfvzefTDbbw4Y5hXz5TnzCeKUcAhrfURAKXUYmAa0EKXCiFEA4sFCg4YgSFrvfFqLw33aXH9IWmsERSSLoKweKcXVRjundCb0iozL60w5u5eviOXyGA//jBtkNfOwe3MQNEdyGqyng2MtnHcjUqpCcAB4BGtddbZByil7gXuBUhMTHRCUYXwEDXlRmbV7A2QZV3ONeXnacoHug4xAkPSWEgcCyExrimvsOmxK/pSWl3H22uPAfDuuky6hAXy0GWpbi5Z6zgzUNgKnWc/H38OvK+1rlFKzQXeAiY2O0nrBcACMHo9tXdBhXALraHoiNHwnLXBCA4nd4O22D/PNxC6j4TEMcZcDT1HQWC4a8osHKKUYt7UgRRX1rFsew4Af/32ALGhAdw62vv+2XVmoMgGejZZ7wHkND1Aa13YZPU14P+cWB4h3KuyCI5vMXohZW8yXqtOtXxecKwRFHqONgJDt6Hg61l99EVzPj6K56cPpaiiltWHCgD47892Eh3iz+RBXd1cuvPjzECxEUhVSqVg9GqaAdza9AClVDetda519VpgrxPLI4Tr1FUZKbePbzGqko5vhqLDDpyooEt/Iy3G6eAgA9y8lr+vD6/OGsEtC9ax83gJFg2/WryV9+4Zw4ikKHcXz2FOCxRaa7NS6gHga4zusQu11ruVUr8HNmmtlwEPKaWuBcxAETDbWeURwmnMtZC/1wgKOVuN7ql5e1vuogoQFGUEhR4XQI+R0H0EBEY4v8zCZUIDfPn3nRdw0ytryCispMZs4Z63N/HJfWNJjvWO8SoyMluI82GuMYJA7jYjxXbuNqNdwe6ANisfX4gf1CQojISY3vK00ElkFFRwwytrKKowfleSY4L55P6LiHZRqg8ZmS2EM9SUGWMUTuyEE9shdzvk7QNLnWPnx/Qx0l90H2EsXQeDn2dmFhXOlxwbwmu3j+TW19ZRY7aQUVjJPW9vYtGc0QT6efbobQkUQmhtpNI+sQtO7rQGhp1GjyRHRSYaQSFhmHVJlyok0cyIpChenJHOfYu2oDVsPnaK3yzZwYsz0j16jIUECtG51JRD/j6juqhh2dXyWIWmolKg2xBj7oWEdONVUmAIB00e1I3/vnoAf1hujD1etj2H3nGh/GqS546xkEAhOiZzLRQeNNoT8vYacyzk7XFsVPNpygRx/YzBbN2GGK9dB3fq+RdE+7jromSO5JezaH0mAPO/O0DvLiFcMyTBzSWzTQKF8G51VVB4CPL3W5e9xmvhYftzLJwtIAK6DjIam7sONpa4NGlTEE6hlGLetQPJKKzg50PGcLLHPtxOj6hg0nt63j8i0utJeIeKQuMJoeCAseQfgIL9cOoYLSbEa0qZIDYVugyA+IGNS0RP6X0kXK6kso7r//kzRwoqAIgPD+DzB8bRxQnTqUqvJ9Ex1FYaDchFh42nhELra8FBqCo6/+tFJhmD17r0N5LkxQ+A2L7g63nTVIrOKSLYj4WzL2Dayz9TUlXHydIa5r67mffvHUOAr+f0hJJAIVyrphxOHbUGhNOv1qX0eMvnN6MgKtmoJorrawSELmlGQJDJd4QXSI4N4aVbh3HHwg1YNGzJLOaZpbv50w2DPaYnlAQK0b4s9cYX/qljUHzMeD2VYV2OQkV+667rF2yMS4hNNYJAbCrE9jMGrPkFtedPIITLjU+N46kp/fnjl0YWo8UbsxjYPYJZY5LcXDKDBApxfsy1RiAoyYbiTCjJguIsKMm0rmc7lrrCFmWCqCSI7m0EgNOBIaYPhCWAj3fPEiaEPXPGp7Ant5RPtxpP1v+zbDcDuoV7RE4oCRSiUX0dlJ2A0hwjGDS8HocSa3AoP8l5NR6fzcfPCAZRKRCdYgSF6F7GEpUEJs+ZjF4IV1JK8acbBnMwr4xdx0sxWzS/XLSF5Q+NIzbUve1qEig6g3qzUeVTfgLK84xgUHYCynKtrznG9JoV+bQpCJwWEmc0JEclGa/RKUY7QlQyhHcHH89ppBPCkwT6mXhl5gimvrSa4so6TpRW89D7W3nn7tFunUpVAoU30hrqKqGiACoLjNeKfOtSYASDijwozzeeACoLaZcAAMZsamHdjC/8yJ5Gt9KIHkYQiEw01v2D2+ezhOiEekYH87eb07nzzY1oDWsOF/LXb/fz+JVpbiuTBAp30hrM1VBdAlXFRhqJqlPG+6oi431lkfG+6WtloXFeu1PG00BEd6NNINwaECJ6NHlNkOohIZzskn5deGhiKi9+fxCAl1ccZkRSFBPT3DP3uQSK82WxgLnK6PNfV2F9rYTacqitMLp/1pZZX8uN15pS61JmBIVq63p1iWPpqdtMQUgshMY3LmFdjSeDsHgjKIR1NRYJAkJ4hIcuS2VrVjE/HTB6Cj724Xa++tUEuka4PluA943MHthLb1r0R+uaNv4r15YmS72xzVJv9L7R1td6s5Eeur7Oul5nfEnX10F9jTHPQH2t8WquMf5jN1cbKSLMNUZwqKty0n/yrWAKML78Q2KNqTJD4qzrcY1LmDUoBMeCSf4nEMLbFFXUctWLqzhRanzvjE6J5r17xrSqvaItI7O9L1AkmPSme0PdXYz2Y/KHwEgjJXVQpDHjWWCkkY00KKpxCY6G4BgIijYCgl+wpJwQohNYd6SQW19bh8X6Vf2ry1J55PK+530dSeHhar6Bxhe1f4ixnH4fEGa8DwgF/1Bj3T8UAsON9wFhRvK5wAhjW2CEcS35whdCnMOYXjE8dFkqf/vOaK/4+w8HGdMrhgt7x7isDN4XKIKiYcj1gPXLVflYF4wBW8rH6H6pfIypJ31MxnaTn9GH3+RrbDcFGNtM/kbuH5Ofsc339BJofQ0yMoj6BhkjgH0DZeCXEMKlHpyYyrojhaw7UoTW8MgH2/j64QlEBLumTdH7qp4ke6wQohM6WVrN5L/9xKlKYyrea4Z04x+3DHM4H1Rbqp7kX2MhhPAC8eGBPHfjkIb15TtyWbotxyWfLYFCCCG8xJUDu3LzyJ4N608v3UX2qUqnf64ECiGE8CJPTx1AYrSR/aCs2sxjH27HYnFuE4IECiGE8CKhAb7Mvzmd00Mp1h8t4q21GU79TAkUQgjhZUYkRfHLS/s0rP/ff/aRYZ1O1RkkUAghhBd6cGIqaV3DAKius/Cbj3c4rQpKAoUQQnghf18fnp8+tCGdx4ajRbyz7phTPksChRBCeKlB3SO47+LeDevPfbWPzML27wUlgUIIIbzYg5f1oW+8kf+uqq6e3366k/YeSC2BQgghvFiAr4nnpw9t6AW1+lBBuw/Ek0AhhBBebkiPSGaPTWlY/8PyPRRXtt9cNxIohBCiA3j0ir50s05qVFhRy3Nf7Wu3a0ugEEKIDiA0wJf/uXZgw/rijVlsOFrULteWQCGEEB3EFQO7cuXAxnm1f/vpTurqLW2+rgQKIYToQOZdO5AQfxMAh/LKeWtNRpuvKYFCCCE6kG4RQfxqUmrD+ovfHSS/rKZN13RqoFBKTVZK7VdKHVJKPWljf4BS6gPr/vVKqWRnlkcIITqD2WNT6BUXAkBZjZm/fN22hm2nBQqllAl4GZgCDABuUUoNOOuwu4FTWus+wHzg/5xVHiGE6Cz8fX343TWNX7cfbspu0/Wc+UQxCjiktT6ita4FFgPTzjpmGvCW9f0S4DLl6Lx+QgghzumSfl2Y1D++5QMd4MxA0R3IarKebd1m8xittRkoAWKcWCYhhOg0nr6mP/6+bf+a922HspyLrSeDsxOQOHIMSql7gXutqzVKqV1tLFtHEQsUuLsQHkLuRSO5F43kXjTq19oTnRkosoGeTdZ7AGcnIDl9TLZSyheIAJqNENFaLwAWACilNmmtRzqlxF5G7kUjuReN5F40knvRSCm1qbXnOrPqaSOQqpRKUUr5AzOAZWcdswy4w/r+JuAH3d5pD4UQQrSJ054otNZmpdQDwNeACViotd6tlPo9sElrvQx4A3hHKXUI40lihrPKI4QQonWcWfWE1vpL4Muztv2uyftqYPp5XnZBOxSto5B70UjuRSO5F43kXjRq9b1QUtMjhBDCHknhIYQQwi6PDRSS/qORA/fiUaXUHqXUDqXU90qpJHeU0xVauhdNjrtJKaWVUh22x4sj90Ip9Qvr78ZupdR7ri6jqzjwN5KolFqhlNpq/Tu5yh3ldDal1EKlVN65hhAow9+t92mHUmq4QxfWWnvcgtH4fRjoBfgD24EBZx1zP/Cq9f0M4AN3l9uN9+JSINj6/r7OfC+sx4UBPwHrgJHuLrcbfy9Sga1AlHW9i7vL7cZ7sQC4z/p+AJDh7nI76V5MAIYDu86x/yrgK4wxbGOA9Y5c11OfKCT9R6MW74XWeoXWutK6ug5jzEpH5MjvBcAfgD8D1a4snIs5ci/uAV7WWp8C0FrnubiMruLIvdBAuPV9BM3HdHUIWuufsDEWrYlpwNvasA6IVEp1a+m6nhooJP1HI0fuRVN3Y/zH0BG1eC+UUsOAnlrr5a4smBs48nvRF+irlPpZKbVOKTXZZaVzLUfuxTzgNqVUNkZPzAddUzSPc77fJ4CTu8e2Qbul/+gAHP45lVK3ASOBi51aIvexey+UUj4YWYhnu6pAbuTI74UvRvXTJRhPmauUUoO01sVOLpurOXIvbgHe1Fq/oJS6EGP81iCtddunf/Murfre9NQnivNJ/4G99B8dgCP3AqXUJOD/Addqrds2S4nnaulehAGDgJVKqQyMOthlHbRB29G/kaVa6zqt9VFgP0bg6GgcuRd3Ax8CaK3XAoEYeaA6G4e+T87mqYFC0n80avFeWKtb/oURJDpqPTS0cC+01iVa61itdbLWOhmjveZarXWrc9x4MEf+Rj7D6OiAUioWoyrqiEtL6RqO3ItM4DIApVR/jECR79JSeoZlwO3W3k9jgBKtdW5LJ3lk1ZOW9B8NHLwXfwFCgY+s7fmZWutr3VZoJ3HwXnQKDt6Lr4ErlFJ7gHrgca11oftK7RwO3ovHgNeUUo9gVLXM7oj/WCql3seoaoy1tsc8A/gBaK1fxWifuQo4BFQCdzp03Q54r4QQQrQjT616EkII4SEkUAghhLBLAoUQQgi7JFAIIYSwSwKFEEIIuyRQiE5JKbVSKXXlWdsetmbfXNLCuWusr5copWymClFKfamUirS+L7e+Jpy+tlIqvaNmMBUdjwQK0Vm9T/OxNzOAf2utb7J3otZ6bEsX11pfdXaqDK11TpNrp2P0ZxfC40mgEJ3VEuAapVQAgHU+kwQg+3Quf6XUQKXUBqXUNmvu/lTr9vIm1wlXSn1qnfPhVWu+KZRSGdbR0A2UUslKqV3W0cO/B262XvtmpdRBpVSc9Tgf63wBnTHFhPBAEihEp2QdobwBOJ1RdQbwAWcmSJsLvKi1TsdItpht41KjMEb9DgZ6Azc48Nm1wO8w5g1J11p/ALwLzLQeMgnYrrUuON+fSwhnkEAhOrOm1U8zrOtNrQV+q5R6AkjSWlfZuMYG6zwI9dbzx7WyLAuB263v7wL+3crrCNHuJFCIzuwzjAmvhgNBWustTXdqrd8DrgWqgK+VUhNtXOPsHDityomjtc4CTlo/YzQdd04R4YUkUIhOS2tdDqzE+G/+7KcJlFK9gCNa679jZN0cYuMyo6xZS32Am4HVDn58GUZa9KZex6iC+tD6hCKER5BAITq794GhGNNnnu1mYJdSahuQBrxt45i1wHPALuAo8KmDn7sCGHC6Mdu6bRlGFmCpdhIeRbLHCuEhrBMszddaj3d3WYRoyiPnoxCis1FKPQncR2PPJyE8hjxRCCGEsEvaKIQQQtglgUIIIYRdEiiEEELYJYFCCCGEXRIohBBC2CWBQgghhF3/H3LA0hQg0CxJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(visibilities, entropy, linewidth=3, label = \"Entropy\")\n",
    "ax.plot(visibilities, ptest, linewidth=3, label = \"Purity Test\")\n",
    "ax.set_xlim(0, 1)\n",
    "ax.set_ylim(0, 1.2*max(entropy))\n",
    "ax.set_xticks([0, 0.2, 0.4, 0.6, 0.8, 1])\n",
    "ax.set_yticks([2,1,0.5,0,1])\n",
    "ax.set_xlabel('Visibility')\n",
    "ax.set_ylabel('Entropy')\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the plot, we can see that the entropy is maximum when the visibility is 0. That parameter has the theoretical minimum purity value: 1 / D where D = 4 is the dimension of the Hilbert space. (cfr. wikipedia article). "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}