wave_plate_tsp.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Polarization tomography TSP: six-state scheme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from itertools import product\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## One-qubit case"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The tomography scheme"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We assume a common tomography scheme using wave plates and a single detector:\n",
    "\n",
    "![aaa](images/tomo_scheme_6s.png)\n",
    "\n",
    "* **HWP**: a half-wave plate\n",
    "* **QWP**: a quarter-wave plate\n",
    "* **PBS**: a polarizing beam splitter\n",
    "* **SPD**: a single-photon detector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Wave plate angle specification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The wave plate angles are stored in the `ANGLES_6s` array. The angles for our six projections are as follows:\n",
    "\n",
    "| projection | HWP angle | QWP angle |\n",
    "|------------|-----------|-----------|\n",
    "|      H     |     0.0   |     0     |\n",
    "|      V     |    45.0   |     0     |\n",
    "|      D     |    22.5   |     0     |\n",
    "|      A     |   -22.5   |     0     |\n",
    "|      R     |     0.0   |    45     |\n",
    "|      L     |     0.0   |   -45     |\n",
    "\n",
    "`ANGLES_6s` is defined in correspondence with this table:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "ANGLES_6s = np.array([[0, 0],\n",
    "                      [45, 0],\n",
    "                      [22.5, 0],\n",
    "                      [-22.5, 0],\n",
    "                      [0, 45],\n",
    "                      [0, -45]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The TSP adjacency matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the one-qubit tomography case, we can directly proceed to the construction of the _TSP adjacency matrix_."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "def adj_mat(angles):\n",
    "    \"\"\"Construct an adjacency matrix from the waveplate proj. angles.\n",
    "    \n",
    "    Args:\n",
    "       * angles (numpy.array): A list of wave plate angles\n",
    "                               for our tomographic projections.\n",
    "    \n",
    "    Returns:\n",
    "       * adjacency_matrix (numpy.array): The appropriate\n",
    "                                         TSP adjacency matrix.\n",
    "    \n",
    "    Note:\n",
    "       * The 'distance' between any two measurement settings is taken as\n",
    "         the maximum of the absolute-valued element-wise difference between\n",
    "         the waveplate angles of the first and the second measurement.\n",
    "    \"\"\"\n",
    "    # Compute the wave plate angle differences for the i -> j transition.\n",
    "    # The result is an i,j-indexed matrix, element i,j is a vector with\n",
    "    # two elements that indicate HWP angle change and QWP angle change.\n",
    "    angular_diff = angles - angles[:, None]\n",
    "    # Apply absolute value to this matrix, then take the highest element\n",
    "    # of each of the size-2 i,j vectors (need to apply max along axis 2\n",
    "    # of the numpy array).\n",
    "    adjacency_matrix = np.max(np.abs(angular_diff), axis=2)\n",
    "    return adjacency_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now look at the adjacency matrix for one-qubit, six-state tomography:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0. , 45. , 22.5, 22.5, 45. , 45. ],\n",
       "       [45. ,  0. , 22.5, 67.5, 45. , 45. ],\n",
       "       [22.5, 22.5,  0. , 45. , 45. , 45. ],\n",
       "       [22.5, 67.5, 45. ,  0. , 45. , 45. ],\n",
       "       [45. , 45. , 45. , 45. ,  0. , 90. ],\n",
       "       [45. , 45. , 45. , 45. , 90. ,  0. ]])"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "adj_mat(ANGLES_6s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7f31da1f9198>"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASQAAAD+CAYAAACTK3gRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAE2tJREFUeJzt3V2sXNV5xvH/gw01OFA+7CCKobYKJUGRIMmRm4g0IlAiBxBYVYQwDYoqVN8kqWkiBehNFKmViFTl4yIXtQINVTAEAU4RrQCLmEZIqYMPOAm2oTh8BFsONgkUSEpS4OnF3ocMls/MHp89M+t4Pz9pNLP3zFnzcoxfr73WXuuVbSIiSnDEpAOIiJiRhBQRxUhCiohiJCFFRDGSkCKiGElIEVGMIhOSpFWSnpS0S9L1BcRzs6R9kh6fdCwAkk6TtFnSDknbJa2bcDyLJP1I0o/reL48yXhmSFog6TFJ9046FgBJz0r6qaRtkrZOOp4SqbT7kCQtAP4buAjYDTwCrLG9Y4IxfRR4DfhX2++bVBw98ZwCnGL7UUnHAtPA6kn9jiQJWGz7NUlHAg8D62z/1yTi6Ynr88AUcJztSycZSx3Ps8CU7RcnHUupSuwhrQR22X7a9u+A24HLJxmQ7R8Av5pkDL1s77X9aP36VWAncOoE47Ht1+rDI+vHRP+lk7QMuAT41iTjiOGUmJBOBZ7vOd7NBP+ylU7ScuD9wJYJx7FA0jZgH7DJ9kTjAb4OfBF4a8Jx9DLwgKRpSWsnHUyJSkxI0ZCkdwF3AdfafmWSsdh+0/a5wDJgpaSJXdpKuhTYZ3t6UjHM4iO2PwB8AvhMPRQQPUpMSHuA03qOl9Xnokc9VnMXcKvtuycdzwzbLwObgVUTDOM84LJ6zOZ24AJJ35lgPADY3lM/7wM2Ug1PRI8SE9IjwJmSVkg6CrgSuGfCMRWlHkS+Cdhp+6sFxLNU0vH166OpJiSemFQ8tm+wvcz2cqr/f75v+1OTigdA0uJ6AgJJi4GPA0XM2pakuIRk+w3gs8D9VIO1d9jePsmYJN0G/BA4S9JuSddMMh6qHsDVVP/yb6sfF08wnlOAzZJ+QvUPyibbRUy1F+Rk4GFJPwZ+BPy77fsmHFNxipv2j4juKq6HFBHdlYQUEcVIQoqIYiQhRUQxik5Ipd3Nmnj6Ky0eKC+m0uJpi6R1kh6vF1dfW587UdImSU/VzycMaqfohASU9oeXePorLR4oL6bS4pmz+q78v6G60fMc4FJJZwDXAw/aPhN4sD7uq/SEFBHley+wxfZv6vsI/xP4S6pF8bfUn7kFWD2ooZHch7RkyRIvX758zu3s37+fpUuXzj2gZ9pZ0rT/dVi6qJWmYMXcm9i/H9r49bSl1XieaaeZVv/MWtBWPM++Bi/+1ppLGwulxn/734LtwOs9p9bbXg8g6b3AvwEfBv6Xqje0Fbja9swd/AJemjmeNabh/hOaWb58OVu3FrT/1FVz+nMbjQ3pnPZ1VUmL9Mszdf/c2zBwdMPP/hpetz110HbsnZK+AjxQfZRtwJsHfMaSBua//K2I6LAjGj4GsX2T7Q/a/ijwEtUmiy/UmwnObCq4r0k8EdFB4ve76Q16DGxLenf9fDrV+NEGqkXxn64/8mmqy7q+RnLJFhHlE7CgvebuknQS8H/AZ2y/LOlG4I56MfpzwBWDGklCiuiwthKS7T8/yLlfAhcO004SUkRHifLGbJKQIjqsxUu2ViQhRXRUy2NIrUhCiuiomVm2kiQhRXRYxpAiogi5ZIuIYiQhRURRcskWEUUQcNSkgzhAowQpaZWkJyXtkjRwk6WIKN/MjZFtLK5ty8DvkrQA+CZVPfKzgTWSzh51YBExegsaPsalSfJbCeyy/bTt31HVSr98tGFFxKjNyx4ScCrwfM/x7vrcO0haK2mrpK379+9vK76IGKH52ENqxPZ621O2p1rZdjYiRmpm2r+khNRklm0PcFrP8bL6XETMYyUuHWnSQ3oEOFPSCklHAVdS7QQXEfNcWz0kSX9X12R7XNJtkhbVOWNLPTv/3Tp/9DUwIdVlTT4L3A/sBO6wvb1BjBFRsLYGtSWdCvwtMGX7fVQ57ErgK8DXbJ9Btc/2NYNiajSGZPs/bP+p7T+x/Y9NfiYiytfiGNJC4GhJC4FjgL3ABcCd9fuN6rKVdud4RIxJWz0k23uAfwJ+TpWI/geYBl6ur7Bgltn5AyUhRXTUzNKRJg9gycxtPfXj7ZLgkk6gujdxBfBHwGJg1aHElLVsER015J7aL85WKBL4C+AZ2/sBJN0NnAccL2lh3UtqNDufHlJEh7U0hvRz4EOSjqlLZl8I7AA2A5+sP9OoLlsSUkRHtXVjpO0tVIPXjwI/pcor64HrgM9L2gWcBNw0KKZcskV0WFs9EttfAr50wOmnqdbCNpaEFNFR2TEyIopR4tKRJKSIDksPKSKKkFLaEVGU9JAiogjdGdR+Zhqu0kiaPiQbPOkIDqKs/xW+p7cmHcI7rF4z6Qi6IZdsEVGEzLJFRDG6c8kWEcVLQoqIomQMKSKKkB5SRBQjg9oRUZT0kCKiCFk6EhFFKa2HVFqCjIgxaWvHSElnSdrW83hF0rWSTpS0SdJT9fMJg2JKQorosJbKID1p+1zb5wIfBH4DbASuBx60fSbwYH08MJ6I6KAhyyA1dSHwM9vPUZVGuqU+36hQZMaQIjpqyEHtJZK29hyvt73+IJ+7Eritfn2y7b31618AJw/6kiSkiI4a8sbIfnXZqvako4DLgBsOfM+2JQ3cdiOXbBEd1sYYUo9PAI/afqE+fkHSKQD1874m8fQl6WZJ+yQ93jyuiChdW7NsPdbw+8s1gHuoCkRCi4Uiv80h1umOiLK1lZAkLQYuAu7uOX0jcJGkp6jKbd84qJ2BY0i2fyBpeYOYImIeaXMtm+1fU1Wn7T33S6pZt8ZaG9SWtBZYC3D6MW21GhGjUuJq/9YGtW2vtz1le2rporZajYhRanlQe84y7R/RUSX2kJKQIjqstPt+mkz73wb8EDhL0m5J14w+rIgYtREtHZmTJrNsqZAVcRjKfkgRUZSMIUVEETKoHRFFySVbRBQhPaSIKEbKIEVEMdJDioiiZAwpIoqQHlJEFKW0hFRajy0ixmTmTu02VvtLOl7SnZKekLRT0odTly0iGpuZZWvyaOAbwH223wOcA+wkddkiYhgtVa79Q+CjwE0Atn9n+2WKqcu2AthQUq4r7UoZ4M1JB/AOq9do0iHEmA05qN2vLtsKYD/wL5LOAaaBdaQuW0QMY4huQ7+6bAuBDwCfs71F0jc44PIsddkioq8WyyDtBnbb3lIf30mVoNqvyxYRh6e2BrVt/wJ4XtJZ9akLgR0cQl22XLJFdFTLN0Z+Dri1Lqf9NPDXVB2eO+pdZp8DrhjUSBJSRIe1dYlkextwsDGmydRli4j5JUtHIqIoSUgRUY6m12xvjTSKtyUhRXTVTB2kJl4fZSC/l4QU0VUF1kFKQorossIGkZKQIrqqwGm2JKSILsslW0QUocCyIwPzo6TTJG2WtEPSdknrxhFYRIxYi6tr29Kkh/QG8AXbj0o6FpiWtMn2jhHHFhGjNt/GkOoNlvbWr1+VtBM4lWo1b0TMV/N92l/ScuD9wJaDvLcWWAtw+uktRBYRo1dYD6lxfpT0LuAu4Frbrxz4vu31tqdsTy1d2maIETES83QMCUlHUiWjW23fPdqQImIsCpxlG5iQJImqmsBO218dfUgRMTYtjSFJehZ4lap6xRu2pySdCHwXWA48C1xh+6W5hnMecDVwgaRt9ePiOcQeESVo/5LtY7bP7SkGMHRdtiazbA/XoUfE4WT0S0cuB86vX98CPARc1+8HCpv0i4ixal5Le4mkrT2PtQe0ZOABSdM976UuW0Q0NFwPqV9dNoCP2N4j6d3AJklP9L7ZtC5bElJEV7U4y2Z7T/28T9JGYCV1XTbbe1OXLSIGa2FQW9LielkZkhYDHwceJ3XZIqKx9paOnAxsrO4QYiGwwfZ9kh4hddkiorEWZtlsPw2cc5DzvyR12SKikewYGRHFmI9LRyLiMJYeUkQUYb7vhxQRh5GMIUVEUdJDGr/vaUyFyYewek1h65U3DLyrf7yuKuz3czhKDykiipFZtogoSnpIEVGEXLJFRFEyqB0RRUgPKSKKkh5SRBRBwFGTDuKdkpAiuqrApSOFhRMRY9ViGSRJCyQ9June+niFpC2Sdkn6rqSB/bEkpIiuar8u2zpgZ8/xV4Cv2T4DeAm4ZlADSUgRXda8DFJfkpYBlwDfqo8FXADcWX/kFmD1oHYyhhTRVcMtHVkiaWvP8Xrb63uOvw58ETi2Pj4JeNn2G/XxbuDUQV+ShBTRVS3VZZN0KbDP9rSk8+cSUhJSRFe1d2PkecBlki4GFgHHAd8Ajpe0sO4lLQP2DGpo4NWhpEWSfiTpx5K2S/ryHIOPiFK0MIZk+wbby2wvB64Evm/7r4DNwCfrjzWqy9ZkUPu3wAW2zwHOBVZJ+lCDn4uIkrU/y3ag64DPS9pFNaZ006AfGHjJZtvAa/XhkfWjsN28IuKQtLyWzfZDwEP166epSmo31mjav77haRtVbe5Ntrcc5DNrJW2VtHX//mFCiIiJmJlla/IYk0YJyfabts+lGphaKel9B/nMettTtqeWLm07zIho3czSkRbuQ2rLUF9l+2WqgapVowknIsZqtGNIQ2syy7ZU0vH166OBi4AnRh1YRIzY6Ae1h9bkPqRTgFskLaBKYHfYvne0YUXEWBS2eKzJLNtPgPePIZaIGKfsGBkRxUgZpIgoRnpIEVGU+TaGFBGHqfSQIqIoSUgRUYQMakdEMQqsOpKEFNFluWSLiCIUOKhdWIctIsaqhdX+s+0qm7psEdFce4trZ9tVNnXZIqKhljZoc+Vgu8oWUpftGeCqt0bS9KFYvWbSEcwDV2nSEbzThvJ2SV6scn5Hr7fVUPMxpL512erdQKaBM4BvAj8jddkiorHhpv1nrcsG1a6ywLn13mkbgfccSkhJSBFd1v4m/y9L2gx8mFHUZYuIw1RLg9qz7Cq7k0Ooy5YeUkRXtbd05KC7ykraAdwu6R+Ax2ijLltEHKZaujFytl1lD6UuWxJSRJcVNmiThBTRVQUuHUlCiuiy9JAiogjpIUVEMbJBW0QUJT2kiChCLtkioigZ1I6IIhTYQ2qcHyUtkPSYpHtHGVBEjFELO0a2aZge0jqqBXPHjSiWiBinI4CBm8qOV6PcJ2kZcAnwrdGGExFjNU97SF8HvggcO9sHJK0F1gKcfszcA4uIMVDTQaQ3RxrGjIG5T9KlwD7b0/0+Z3u97SnbU0sXtRZfRIyMqPokTR7j0eSbzgMuk3QxsAg4TtJ3bH9qtKFFxGjNJKQmfjvKQN42sIdk+wbby2wvB64Evp9kFHE4EFUfo8mjTyvSaZI2S9pR12VbV58/UdImSU/VzycMiqiw26IiYnxau2R7A/iC7bOBDwGfkXQ2cD3woO0zgQfr476Guji0/RDw0DA/ExGlGuaSbXa29wJ769evStpJVfLocuD8+mO3UOWO6/q1lTu1IzprqITUty7b2y1Ky6m2s90CnFwnK4BfACcP+pIkpIhOa7x2pG9dNgBJ7wLuAq61/Yp6CmvatqSB1T+TkCI6q51LNgBJR1Ilo1tt312ffkHSKbb3SjoF2DeonQxqR3TWEbQ0yyaqEkc7bX+15617qOqxQeqyRUR/rfWQzgOuBn4qaVt97u+BG4E7JF0DPAdcMaihJKSITmtllu1hqux2MBeON5qImKfaG0NqS1nRRMQYJSFFRDGSkCKiGDNr2cqRhBTRWekhRUQxkpAiohhJSBGNLNZst7VMzq89cCnW2ExN9V1W1lASUkQUQ8AfTDqId0hCiuis9JAiohhJSBFRjCSkiChKWSmgrGgiYozSQ4qIYsxs0FaO7BgZ0VntVa6VdLOkfZIe7zmXumwRMYzWSml/G1h1wLmh67IlIUV0Vns9JNs/AH51wOnLqeqxUT+vHtROxpAiOmvkg9qpyxYRTbVfKHI2qcsWEQMMtUHbwEKRB5G6bBHRVHtjSLMYTV02Sc8CrwJvAm8cQqaMiOKIIUpp929Jug04n+rSbjfwJUZcl+1jtl88hFgjokjtDWrbXjPLW6nLFhFNzN+lIwYeqEfJ//lgo+uS1gJrAU4/pr0AI2JU5m/VkY/Y3iPp3cAmSU/UN0K9rU5S6wGmTho8vRcRk1ZeD6nRLJvtPfXzPmAjsHKUQUXEOIx8lm1oAxOSpMWSjp15DXwceLz/T0VE+cpLSE2+6WRgo6oqEAuBDbbvG2lUETEmZV2yDYzG9tPAOWOIJSLGqrwxpLKiiYgxKm+DtiSkiM5KDykiilJWCigrmogYo/SQIqIYSUgRUYz2Vvu3JQkporPm71q2iDjs5JItIopRXkLKFrYRndVqochVkp6UtEvSwPprsykrPUbEGLXTQ5K0APgmcBGwG3hE0j22dwzbVhJSRGe1Nqi9EthVr3tF0u1URSLLSEjTv+JF3cZzLTS1BChpH+/E019p8UCLMdU7XsxVW/H88VwbmJ6evl/SkoYfX9SnLtupwPM97+0G/uxQYhpJQrK9tI12JG0tqcJJ4umvtHigvJhKisf2qknHcKAMakfEXO0BTus5XlafG1oSUkTM1SPAmZJWSDoKuJKqSOTQSh/Ublw7fEwST3+lxQPlxVRaPHNm+w1JnwXup1qLcrPt7YfSluwUCImIMuSSLSKKkYQUEcVIQoqIYiQhRUQxkpAiohhJSBFRjCSkiCjG/wMRT1dtQ4fDMgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "c = ax.matshow(adj_mat(ANGLES_6s), cmap=\"hot_r\")\n",
    "fig.colorbar(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Formulating the TSPLIB compliant TSP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use [Concorde](http://www.math.uwaterloo.ca/tsp/concorde.html), a [TSPLIB95](https://www.iwr.uni-heidelberg.de/groups/comopt/software/TSPLIB95/) compliant TSP solver.\n",
    "\n",
    "TSPLIB95 Documentation, which specifies how the TSP file should look like is available [here](https://www.iwr.uni-heidelberg.de/groups/comopt/software/TSPLIB95/tsp95.pdf).\n",
    "\n",
    "The `.tsp` file is a plain-text specification of the TSP. It starts with a header, where the problem name, type, dimension, are given, along with additional information.\n",
    "\n",
    "The header always has the same basic form, and only the name, dimension, and optional comment fields change."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The basic form of the header, with {} placeholders for\n",
    "# name, dimension, and comment.\n",
    "HEADER = \"\"\"\\\n",
    "NAME : {n}\n",
    "TYPE : TSP\n",
    "COMMENT : {c}\n",
    "DIMENSION : {d}\n",
    "EDGE_WEIGHT_TYPE : EXPLICIT\n",
    "EDGE_WEIGHT_FORMAT : FULL_MATRIX\n",
    "EDGE_WEIGHT_SECTION\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `{}` placeholder substitution is achieved using the `format()` function on the `HEADER` string. This is shown below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NAME : 1q_6s.tsp\n",
      "TYPE : TSP\n",
      "COMMENT : Wave plate tomography, 1 qubit, six-state\n",
      "DIMENSION : 6\n",
      "EDGE_WEIGHT_TYPE : EXPLICIT\n",
      "EDGE_WEIGHT_FORMAT : FULL_MATRIX\n",
      "EDGE_WEIGHT_SECTION\n",
      "\n"
     ]
    }
   ],
   "source": [
    "name = \"1q_6s.tsp\"\n",
    "dimension = 6\n",
    "comment = \"Wave plate tomography, 1 qubit, six-state\"\n",
    "\n",
    "print(HEADER.format(n=name, c=comment, d=dimension))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The adjacency matrix follows the header. The file ends with the string `\"EOF\"`.\n",
    "\n",
    "All that remains is to save the .tsp file, with the adjacency matrix of the problem. It is useful to define a function for this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "def save_tsp(filename, adjacency_matrix, comment=\"\"):\n",
    "    \"\"\"Save a given adjacency matrix in TSPLIB95 form.\n",
    "    \n",
    "    Args:\n",
    "       * filename (str): The name of the .tsp file.\n",
    "       * adjacency_matrix (numpy.array): The adjacency matrix.\n",
    "    \n",
    "    Kwargs:\n",
    "       * comment (str): The comment that appears in the header\n",
    "                        of the .tsp file (optional).\n",
    "                        Defaults to an empty string.\n",
    "    \n",
    "    Notes:\n",
    "       We ran into trouble when solving the TSP with float matrices.\n",
    "       As we will encounter matrix elements like 67.5, we multiply\n",
    "       the adjacency matrix by 2 for the solver.\n",
    "       \n",
    "       As each adjacency matrix element has been multiplied by the\n",
    "       same factor, this will not affect the solution. When we evaluate\n",
    "       the speedup, we use adjacency matrices which **were not** multiplied\n",
    "       to provide the correct results.\n",
    "    \"\"\"\n",
    "    dimension = len(adjacency_matrix)  # determine the dim of TSP\n",
    "    adjacency_matrix *= 2  # multiply by 2 to avoid floats\n",
    "    np.savetxt(filename, adjacency_matrix,\n",
    "               header=HEADER.format(n=filename, d=dimension, c=comment),\n",
    "               footer=\"EOF\", comments=\"\", fmt=\"%d\")  # converted to int"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save the one-qubit, six-state tomography TSP:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = \"1q_6s.tsp\"\n",
    "comment = \"Wave plate tomography, 1 qubit, six-state\"\n",
    "\n",
    "save_tsp(filename, adj_mat(ANGLES_6s), comment)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solving the TSP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Call Concorde. For this to work, make sure that you can run concorde from the terminal (its location must be in the `PATH` variable). Alternatively, placing the `concorde` executable in this directory should work as well.\n",
    "\n",
    "The `-x` parameter allows concorde to delete temporary files after solving the TSP.\n",
    "\n",
    "_Also note that this way of calling external commands is specific to Jupyter notebooks. In ordinary Python, functions like `call()` from the `subprocess` module would have to be used._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "concorde -x 1q_6s.tsp\n",
      "Host: nadezdan  Current process id: 6360\n",
      "Using random seed 1542024784\n",
      "Problem Name: 1q_6s.tsp\n",
      "Problem Type: TSP\n",
      "Wave plate tomography, 1 qubit, six-state\n",
      "Number of Nodes: 6\n",
      "Explicit Lengths (CC_MATRIXNORM)\n",
      "Optimal Solution: 450.00\n",
      "Total Running Time: 0.00 (seconds)\n"
     ]
    }
   ],
   "source": [
    "!concorde -x 1q_6s.tsp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another file, `1q_6s.sol` should appear in the working directory. It contains, in plain text, a sequence of node indexes that form the solution of the TSP. The numbers 0 through 5 can be identified with the six projection measurements `H`, `V`, `D`, `A`, `R`, `L`, in this order.\n",
    "\n",
    "Opening the file, we can see that in our case, the sequence (or _path_)is:\n",
    "\n",
    "1. 0 (H)\n",
    "2. 5 (L)\n",
    "3. 3 (A)\n",
    "4. 4 (R)\n",
    "5. 1 (V)\n",
    "6. 2 (D)\n",
    "\n",
    "We define a function that extracts this sequence from the `.sol` file. The function is more complicated than seems necessary at this point. However, it will work with the `.sol` files of higher-dimension TSPs, where the solution is split onto multiple lines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_projs(filename):\n",
    "    \"\"\"Retrieves the order of projections from the .sol file.\n",
    "    \n",
    "    Args:\n",
    "       * filename (str): The name of the .sol file.\n",
    "    \n",
    "    Returns:\n",
    "       * path (list): The TSP-optimal sequence of projection indices.\n",
    "    \"\"\"\n",
    "    with open(filename, \"r\") as f:  # open the file\n",
    "        # collect the lines of the file\n",
    "        # skip the first line, strip endline characters\n",
    "        lines = map(lambda x: x.rstrip(\" \\n\"), f.readlines()[1:])\n",
    "    # join the lines into space-separated string of indices\n",
    "    joined_lines = \" \".join(list(lines))\n",
    "    # convert to a list of integers\n",
    "    path = list(map(int, joined_lines.split(\" \")))\n",
    "    return path"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can test this function on our solution to the one-qubit, six-state tomography TSP:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0, 5, 3, 4, 1, 2]\n"
     ]
    }
   ],
   "source": [
    "tsp_path = get_projs(\"1q_6s.sol\")\n",
    "print(tsp_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluating speedup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now use the TSP-optimal path to determine the total (angular) duration of optimized one-qubit, six-state tomography. To do this, we convert the path into a list of _edges_, or transitions between the individual measurements. Then, we will be able to sum the adjacency matrix elements that correspond to these edges to obtain the total duration.\n",
    "\n",
    "We are interested in the _cycle length_, meaning that we consider the transition from the last to the first measurement (so that our tomography ends where it started off)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cycle_length(adj_matrix, path):\n",
    "    \"\"\"Compute the length of a cycle for a given adj. matrix and path.\n",
    "    \n",
    "    Args:\n",
    "       * adj_matrix (numpy.array): The appropriate adjacency matrix.\n",
    "       * path (list): A list of projection indices.\n",
    "    \n",
    "    Returns:\n",
    "       * cycle_length (float): The cycle length.\n",
    "    \"\"\"\n",
    "    # convert the path to a list of edges, or transitions\n",
    "    # create pairs of nth and (n+1)th path indices\n",
    "    edges = [[path[n], path[n+1]] for n in range(len(path) - 1)]\n",
    "    # add the final transition that completes the cycle\n",
    "    edges += [[path[-1], path[0]]]\n",
    "    cycle_length = sum([adj_matrix[edge[0]][edge[1]] for edge in edges])\n",
    "    return cycle_length"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see the TSP-optimal cycle length:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "225.0\n"
     ]
    }
   ],
   "source": [
    "tsp_cycle = cycle_length(adj_mat(ANGLES_6s), tsp_path)\n",
    "print(tsp_cycle)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The conventional ordering of tomography goes through the polarization projection in the order `H`, `V`, `D`, `A`, `R`, `L`. When we defined `ANGLES_6s`, we gave the wave the wave plate angles for the projections in this order. We can then obtain the conventional cycle length by giving `[0, 1, 2, 3, 4, 5]`, or `range(6)` as the `path` argument to the function `cycle_length()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "292.5\n"
     ]
    }
   ],
   "source": [
    "conv_cycle = cycle_length(adj_mat(ANGLES_6s), range(6))\n",
    "print(conv_cycle)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define the TSP speedup as the ratio between the conventional and the TSP-optimal cycle length:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.3\n"
     ]
    }
   ],
   "source": [
    "speedup = conv_cycle / tsp_cycle\n",
    "print(speedup)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The _reduction_ of total duration is also of importance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "67.5\n"
     ]
    }
   ],
   "source": [
    "reduction = conv_cycle - tsp_cycle\n",
    "print(reduction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiple qubits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The procedure is almost identical for multiple qubits. The only exception is that `ANGLES_6s` cannot be used anymore. For an $n$-qubit system, $6^n$ projections need to be measured. Furthermore, $2\\times n$ wave plates are required. We need to construct a wave plate angle list, similar to `ANGLES_6s`, for all the projections. This can be understood as a 'Cartesian power' of `ANGLES_6s`\n",
    "\n",
    "We will use the `product()` function from the `itertools` module to generate said Cartesian power."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [],
   "source": [
    "def nqubit_angles(one_qubit_angles, n):\n",
    "    cartesian = product(one_qubit_angles, repeat=n)\n",
    "    # We have just generated a list of n-tuples of 1x2 arrays.\n",
    "    # To convert it to a list of 1x2n arrays, stack the n 1x2 arrays\n",
    "    # horizontally for each n-tuple of the list.\n",
    "    array_list = map(np.hstack, cartesian)\n",
    "    # Finally, convert the list of 1x2n arrays to one 2D array\n",
    "    angles = np.vstack(array_list)\n",
    "    return angles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To check, see how the 2-qubit wave plate angles array turns out:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  0.    0.    0.    0. ]\n",
      " [  0.    0.   45.    0. ]\n",
      " [  0.    0.   22.5   0. ]\n",
      " [  0.    0.  -22.5   0. ]\n",
      " [  0.    0.    0.   45. ]\n",
      " [  0.    0.    0.  -45. ]\n",
      " [ 45.    0.    0.    0. ]\n",
      " [ 45.    0.   45.    0. ]\n",
      " [ 45.    0.   22.5   0. ]\n",
      " [ 45.    0.  -22.5   0. ]\n",
      " [ 45.    0.    0.   45. ]\n",
      " [ 45.    0.    0.  -45. ]\n",
      " [ 22.5   0.    0.    0. ]\n",
      " [ 22.5   0.   45.    0. ]\n",
      " [ 22.5   0.   22.5   0. ]\n",
      " [ 22.5   0.  -22.5   0. ]\n",
      " [ 22.5   0.    0.   45. ]\n",
      " [ 22.5   0.    0.  -45. ]\n",
      " [-22.5   0.    0.    0. ]\n",
      " [-22.5   0.   45.    0. ]\n",
      " [-22.5   0.   22.5   0. ]\n",
      " [-22.5   0.  -22.5   0. ]\n",
      " [-22.5   0.    0.   45. ]\n",
      " [-22.5   0.    0.  -45. ]\n",
      " [  0.   45.    0.    0. ]\n",
      " [  0.   45.   45.    0. ]\n",
      " [  0.   45.   22.5   0. ]\n",
      " [  0.   45.  -22.5   0. ]\n",
      " [  0.   45.    0.   45. ]\n",
      " [  0.   45.    0.  -45. ]\n",
      " [  0.  -45.    0.    0. ]\n",
      " [  0.  -45.   45.    0. ]\n",
      " [  0.  -45.   22.5   0. ]\n",
      " [  0.  -45.  -22.5   0. ]\n",
      " [  0.  -45.    0.   45. ]\n",
      " [  0.  -45.    0.  -45. ]]\n"
     ]
    }
   ],
   "source": [
    "ANGLES_6s_2qubit = nqubit_angles(ANGLES_6s, 2)\n",
    "print(ANGLES_6s_2qubit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a list of $6^2 = 36$ polarization projections `HH`, `HV`, ..., `HR`, `HL`, `VH`, `VV`, ..., `LR`, `LL`. The four angles now correspond to $\\mathrm{HWP}_1$, $\\mathrm{QWP}_1$, $\\mathrm{HWP}_2$, $\\mathrm{QWP}_2$, respectively, with the index corresponding to qubit 1 or qubit 2 of the two-qubit system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the rest is identical to the one-qubit case, we will now proceed faster and generate, solve, and evaluate the tomograpy TSPs for up to 4-qubit system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "metadata": {},
   "outputs": [],
   "source": [
    "for nqubits in range(1, 5):\n",
    "    nqu_angles = nqubit_angles(ANGLES_6s, nqubits)\n",
    "    adjacency_matrix = adj_mat(nqu_angles)\n",
    "    \n",
    "    tsp_filename = \"{}q_6s.tsp\".format(nqubits)\n",
    "    tsp_comment = \"Wave plate tomography, {} qubits, six-state\".format(nqubits)\n",
    "\n",
    "    save_tsp(tsp_filename, adjacency_matrix, tsp_comment)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With all the TSP files generated, call Concorde for each (except one-qubit, it has been solved already). In the summary Concorde gives upon completion, notice how execution time increases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "concorde -x 2q_6s.tsp\n",
      "Host: nadezdan  Current process id: 7801\n",
      "Using random seed 1542028785\n",
      "Problem Name: 2q_6s.tsp\n",
      "Problem Type: TSP\n",
      "Wave plate tomography, 2 qubits, six-state\n",
      "Number of Nodes: 36\n",
      "Explicit Lengths (CC_MATRIXNORM)\n",
      "Set initial upperbound to 2025 (from tour)\n",
      "  LP Value  1: 2025.000000  (0.01 seconds)\n",
      "New lower bound: 2025.000000\n",
      "Final lower bound 2025.000000, upper bound 2025.000000\n",
      "Exact lower bound: 2025.000000\n",
      "DIFF: 0.000000\n",
      "Final LP has 45 rows, 77 columns, 221 nonzeros\n",
      "Optimal Solution: 2025.00\n",
      "Number of bbnodes: 1\n",
      "Total Running Time: 0.09 (seconds)\n"
     ]
    }
   ],
   "source": [
    "!concorde -x 2q_6s.tsp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 141,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "concorde -x 3q_6s.tsp\n",
      "Host: nadezdan  Current process id: 7806\n",
      "Using random seed 1542028799\n",
      "Problem Name: 3q_6s.tsp\n",
      "Problem Type: TSP\n",
      "Wave plate tomography, 3 qubits, six-state\n",
      "Number of Nodes: 216\n",
      "Explicit Lengths (CC_MATRIXNORM)\n",
      "Set initial upperbound to 10935 (from tour)\n",
      "  LP Value  1: 10800.000000  (0.29 seconds)\n",
      "  LP Value  2: 10890.000000  (0.46 seconds)\n",
      "  LP Value  3: 10935.000000  (0.93 seconds)\n",
      "  LP Value  4: 10935.000000  (1.55 seconds)\n",
      "New lower bound: 10935.000000\n",
      "Final lower bound 10935.000000, upper bound 10935.000000\n",
      "Exact lower bound: 10935.000000\n",
      "DIFF: 0.000000\n",
      "Final LP has 479 rows, 569 columns, 7429 nonzeros\n",
      "Optimal Solution: 10935.00\n",
      "Number of bbnodes: 1\n",
      "Total Running Time: 5.05 (seconds)\n"
     ]
    }
   ],
   "source": [
    "!concorde -x 3q_6s.tsp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "concorde -x 4q_6s.tsp\n",
      "Host: nadezdan  Current process id: 7814\n",
      "Using random seed 1542028810\n",
      "Problem Name: 4q_6s.tsp\n",
      "Problem Type: TSP\n",
      "Wave plate tomography, 4 qubits, six-state\n",
      "Number of Nodes: 1296\n",
      "Explicit Lengths (CC_MATRIXNORM)\n",
      "Set initial upperbound to 61965 (from tour)\n",
      "  LP Value  1: 61065.000000  (0.67 seconds)\n",
      "  LP Value  2: 61875.000000  (2.08 seconds)\n",
      "  LP Value  3: 61920.000000  (6.20 seconds)\n",
      "  LP Value  4: 61920.000000  (51.25 seconds)\n",
      "New lower bound: 61919.999999\n",
      "  LP Value  1: 61965.000000  (81.99 seconds)\n",
      "New lower bound: 61965.000000\n",
      "Final lower bound 61965.000000, upper bound 61965.000000\n",
      "Exact lower bound: 61965.000000\n",
      "DIFF: 0.000000\n",
      "Final LP has 2309 rows, 3162 columns, 76413 nonzeros\n",
      "Optimal Solution: 61965.00\n",
      "Number of bbnodes: 1\n",
      "Total Running Time: 121.09 (seconds)\n"
     ]
    }
   ],
   "source": [
    "!concorde -x 4q_6s.tsp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can evaluate how the TSP speedup and total duration reduction scales with the number of qubits in the system:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 qubits:\n",
      "225.0\n",
      "292.5\n",
      "1.3\n",
      "67.5\n",
      "2 qubits:\n",
      "1012.5\n",
      "1800.0\n",
      "1.7777777777777777\n",
      "787.5\n",
      "3 qubits:\n",
      "5467.5\n",
      "10845.0\n",
      "1.9835390946502058\n",
      "5377.5\n",
      "4 qubits:\n",
      "30982.5\n",
      "65115.0\n",
      "2.101670297748729\n",
      "34132.5\n"
     ]
    }
   ],
   "source": [
    "speedups, reductions = [], []\n",
    "for nqubits in range(1, 5):\n",
    "    print(\"{} qubits:\".format(nqubits))\n",
    "    nqu_angles = nqubit_angles(ANGLES_6s, nqubits)\n",
    "    adjacency_matrix = adj_mat(nqu_angles)\n",
    "    sol_filename = \"{}q_6s.sol\".format(nqubits)\n",
    "    solution = get_projs(sol_filename)\n",
    "    \n",
    "    optimal_cycle = cycle_length(adjacency_matrix, solution)\n",
    "    conv_cycle = cycle_length(adjacency_matrix, range(len(adjacency_matrix)))\n",
    "    speedup = conv_cycle / optimal_cycle\n",
    "    reduction = conv_cycle - optimal_cycle\n",
    "    print(optimal_cycle)\n",
    "    print(conv_cycle)\n",
    "    print(speedup)\n",
    "    print(reduction)"
   ]
  }
 ],
 "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}