Qutrit Decay Dynamics.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import scipy as sc\n",
    "import inspect\n",
    "import copy\n",
    "from scipy.linalg import expm\n",
    "import qutip as qt\n",
    "from functools import reduce\n",
    "import itertools\n",
    "import seaborn\n",
    "from numpy import exp\n",
    "\n",
    "## Defining qutrit operators\n",
    "\n",
    "P0 = qt.basis(3,0)*qt.basis(3,0).dag()\n",
    "P1 = qt.basis(3,1)*qt.basis(3,1).dag()\n",
    "P2 = qt.basis(3,2)*qt.basis(3,2).dag()\n",
    "\n",
    "Id = qt.Qobj(np.eye(3))\n",
    "\n",
    "import sys\n",
    "sys.path.append('/home/qnl')\n",
    "\n",
    "import time\n",
    "\n",
    "from importlib import reload\n",
    "\n",
    "from qnl_analysis import SimTools as ST"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $T_1$ decay of a qutrit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For $T_1$ decay, the operators we need are $C_0$ and $C_1$ below, essentially the qutrit analogues of the qubit $\\sigma_-$ lowering operator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "gamma1GE = 1./(10e-6) # 1/s\n",
    "gamma1EF = 1./(5e-6) # 1/s\n",
    "\n",
    "C0 = np.sqrt(gamma1GE)*qt.Qobj(np.array([[0,1,0],\n",
    "                                         [0,0,0],\n",
    "                                         [0,0,0]]))\n",
    "C1 = np.sqrt(gamma1EF)*qt.Qobj(np.array([[0,0,0],\n",
    "                                         [0,0,1],\n",
    "                                         [0,0,0]]))\n",
    "\n",
    "times = np.linspace(0,30e-6,50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = qt.Qobj(0*np.eye(3))\n",
    "\n",
    "psi0 = qt.Qobj(np.array([[0, 0, 0],\n",
    "                        [0, 0.5, 0.5],\n",
    "                        [0, 0.5, 0.5]]))\n",
    "\n",
    "result = qt.mesolve(H, psi0, times, [C0, C1],[])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x120645e80>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl8AAAE/CAYAAAB4o6baAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd0VNXax/Hvk15pSWgJNYQOUkIXQQQBG00RLICKXuy93+u1N16vFetVsSEiCCKCSFekgwhSAgEChBJC6Anp+/3jBG6MIRlgMmfK81krK5k57RcWbJ7ZZ5+9xRiDUkoppZRyDT+7AyillFJK+RItvpRSSimlXEiLL6WUUkopF9LiSymllFLKhbT4UkoppZRyIS2+lFJKKaVcSIsvVWFEpKeIpJ7H8e+LyL+cmUkp5VvOtx1SqiJo8eVDRCRFRE6KyAkRSRORT0Ukwu5cACIySkQWF3/PGDPGGPOcXZmUUu5FRK4TkVVFbdg+EZklIhfanUups6XFl++50hgTAbQDOgD/tDmPUkqVS0QeAN4AXgRqAHWBd4EBLswQ4KprKe+mxZePMsbsAWYBLUWktohMF5FDIpIsIree2k9EnhaRySLyjYgcF5E1InJBse1GRBoVez1eRJ4v7Zoi8piIbCs6z0YRGVT0fjPgfaBL0SfaI6WdS0RuLcp3qChv7RI5xojIVhE5LCLjRESKtjUSkUUiclREDorIN876c1RKVTwRqQw8C9xpjPnOGJNpjMkzxvxgjHlYRIJF5A0R2Vv09YaIBJc4x4MicqCox+ymYu8Hi8j/iciuojsC74tIaNG2niKSKiKPish+4NOi968QkbUickRElohI62LnSxGRh0RkXVGb842IhBTbPqDo2GNF7WG/U7+jiHxclG+PiDwvIv5F27QN8zJafPkoEakDXAb8DnwNpAK1gauBF0XkkmK7DwC+BaoBE4BpIhJ4DpfdBnQHKgPPAF+KSC1jzCZgDLDUGBNhjKlSSt5ewEvAUKAWsBOYWGK3K7B68y4o2q9v0fvPAT8DVYE44O1zyK6Usk8XIASYeobtTwKdgTZY//478tde/ZpY7U4scAswTkSqFm17BWhcdGyjon2eKnFsNaAecJuItAM+Af4BRAEfANNLFHtDgX5AA6A1MApARDoCnwMPA1WAi4CUomM+A/KLMrQFLgVGF23TNszLaPHle6YV9SwtBhYBHwIXAo8aY7KNMWuB/wI3FjtmtTFmsjEmD/gPViPY+WwvbIz51hiz1xhTaIz5BtiK1Ug64nrgE2PMGmNMDvA4Vk9Z/WL7vGyMOWKM2QUswGpMAfKwGs7aRb/jX8aWKaXcXhRw0BiTf4bt1wPPGmMOGGPSsT7cFW/D8oq25xljZgIngCZFveO3AvcbYw4ZY45j3dYcVuzYQuDfxpgcY8zJov0/MMYsN8YUGGM+A3L4a5v4VlFbdwj4gf+1RbdgtWNzitrBPcaYzSJSA+gP3FfUq3cAeL1YDm3DvIwWX75noDGmijGmnjHmDqzerlONzik7sT79nbL71A/GmEL+10t2VkRkRLGu+iNASyDawcNrF+U6leMEkFEi5/5iP2cBpx4meAQQYIWIbBCRm882u1LKVhlAdBljrv7SPhT9XLyNyihRuJ1qH2KAMGB1sXbpp6L3T0k3xmQXe10PePDU/kXH1ClxvTO1RXWw7gCUVA8IBPYVO+cHQPWi7dqGeRkdPKj2AtVEJLJYAVYX2FNsnzqnfhARP6xu771Fb2VhNV6n1MQqzv5CROoBHwGXYN1eLBCRtVgNCoBxIGe9YucLx/o0vOeMR5w6sTH7sT6tItaTUXNF5BdjTHJ5xyql3MJSIBsYCEwuZfup9mFD0eu6/K+NKstB4CTQomgcbGlKtk27gReMMS84cP6SdgPxZ3g/B4gurXdP2zDvoz1fPs4YsxtYArwkIiFFA0dvAb4qtlt7ERlc9KnzPqxGYlnRtrXAdSLiXzRwtMcZLhWO1YilAxQNeG1ZbHsaECciQWc4fgJwk4i0KRpb8SKw3BiTUt7vKCLXiEhc0cvDRTkKyjtOKeUejDFHscZhjRORgSISJiKBItJfRF7FGrf6TxGJEZHoon2/dOC8hVgfCl8XkeoAIhIrIn3LOOwjYIyIdBJLuIhcLiKRDvwqH2O1Y5eIiF/RtZoaY/Zhjel6TUQqFW2LF5EeRZm0DfMyWnwpgOFAfaxPilOxxjfMKbb9e+BarH/0NwKDi8Z/AdwLXAkcwRp3Ma20CxhjNgKvYX2CTQNaAb8V22U+1qfW/SJysJTj5wH/AqYA+7A+PQ4rud8ZdACWi8gJYDpwrzFmh4PHKqXcgDHmP8ADWAPp07F6i+7CanOeB1YB64D1wJqi9xzxKJAMLBORY8BcoEkZOVZh9UK9g9UmJlM0oN6B32EFcBPWeK6jWONuT/XojwCCgI1F552M9XARaBvmdcSY8u72KF8mIk8DjYwxN9idRSmllPIG2vOllFJKKeVCWnwppZRSSrmQ3nZUSimllHIh7flSSimllHIhLb6UUkoppVzItklWo6OjTf369e26vFLKBqtXrz5ojIkpf0/3pu2XUr7Hme2XbcVX/fr1WbVqlV2XV0rZQER2lr+X+9P2Synf48z2S287KqWUUkq5kBZfSimllFIupMWXUkoppZQL2TbmqzR5eXmkpqaSnZ1td5TzFhISQlxcHIGBgXZHUUq5gDe0X9puKeUablV8paamEhkZSf369RERu+OcM2MMGRkZpKam0qBBA7vjKKVcwNPbL223lHIdt7rtmJ2dTVRUlEc2XMWJCFFRUR79CVgpdXY8vf3Sdksp13Gr4gvw2IarJG/5PZRSjvP0f/eenl8pT+F2xZc7SEtL47rrrqNhw4a0b9+eLl26MHXqVABWrFhBz549SUhIoF27dlx++eWsX7/e5sRKqTMRkX4ikiQiySLyWCnbR4lIuoisLfoabUdOZzlT+7Vw4UIqV65MmzZtTn/NnTvX7rhK+SS3GvPlDowxDBw4kJEjRzJhwgQAdu7cyfTp00lLS2Po0KFMmDCBrl27ArB48WK2bdtGq1at7IytlCqFiPgD44A+QCqwUkSmG2M2ltj1G2PMXS4P6GRltV9Vq1ale/fuzJgxw+aUSint+Sph/vz5BAUFMWbMmNPv1atXj7vvvpt33nmHkSNHni68AC688EIGDhxoR1SlbJOTX8Cy7RlMWrXb7ijl6QgkG2O2G2NygYnAAJcmyD5qfblAWe2XUsp9aM9XCRs2bKBdu3Zn3DZy5EgXJ1LKfgWFhnWpR1iyLYOl2zJYmXKInPxCwoL8GdQ2lkB/t/0cFwsUrxBTgU6l7DdERC4CtgD3G2P+VlWKyG3AbQB169Z1PMHxNOt7SGXHjzlHZbVfAL/++itt2rQ5/XrKlCnEx8dXeC6l1F+5bfH1zA8b2Lj3mFPP2bx2Jf59ZYuzOubOO+9k8eLFBAUFUadOnb9s69SpE8eOHePSSy/lzTffdGZUpWyXejiLX7ce5Net6SzeepBj2fkANK0ZyXWd6tI1PpqODaq5c+EFUNoIclPi9Q/A18aYHBEZA3wG9PrbQcZ8CHwIkJiYWPIcf/GX9qsgBwryICjj7NMXc77t19ixY/W2o1Juwm2LL7u0aNGCKVOmnH49btw4Dh48SGJiIv369WPNmjUMGGDdtVi+fDmTJ0/Wxkx5hZz8AlbuOMz8zQdYmHSA7QczAahVOYT+LWtxYUI0XeKjiI4ItjnpWUkFin9qigP2Ft/BGFO8KvoIeMWpCcQfyANTUPRzxSmr/VJKuQ+3Lb7O9hOes/Tq1YsnnniC9957j9tvvx2ArKwswPoU2alTJ/r27Xt63NepbUp5ovTjOczblMa8zQf4LfkgWbkFBAX40blhFDd0rsdFjaOJj4nw5CkIVgIJItIA2AMMA64rvoOI1DLG7Ct6eRWw6Xwv+pf2Kz8HDmyEynEQHnO+py5TWe2XUsp9uG3xZRcRYdq0adx///28+uqrxMTEEB4eziuvvELNmjX55ptvePTRR9mzZw/Vq1cnOjqap556yu7YSjks+cBxft6YxpyNaazdfQRjILZKKIPbxXJxk+p0jY8mNKhie2hcxRiTLyJ3AbMBf+ATY8wGEXkWWGWMmQ7cIyJXAfnAIWCUU0P4B4FfAORmVnjxVVb7BX8f8/XPf/6Tq6++ukIzKaX+TouvUtSqVYuJEyeWuq1z584sWrTIxYmUOnfGGDbuO8as9fuZ9ec+tqVbtxNbxVbm/t6N6dO8Bk1rRnpy71aZjDEzgZkl3nuq2M+PA49XWAARCAyHXNf0QJXVfh096pqnLpVSZdPiSykvZIxhw95jzFi3j1l/7mNnRhZ+Ap0bRjGya336NK9Brcqhdsf0HUFhkHMUCvLBX5tdpXydtgJKeZFt6SeYvnYvP6zby/b0TAL8hK6Norm9Rzx9mtcgyrMGy3uPoHDre14W+FeyN4tSynZafCnl4dKP5zD9j71M/T2VP/ccQwQ6NajG6Asb0q9lTaqFB9kdUQWGWd9zMyFEiy+lfJ0WX0p5oOy8AuZuSuO7NXtYtCWdgkJDq9jK/OuK5lzRuhY1KoXYHVEV5+cPASGQl2l3EqWUG9DiSykPsmnfMb5ZuZupv+/h6Mk8alYK4dbuDRncLpbGNSLtjqfKEhQOJ4+AMdYgfKWUz9LiSyk3dyInn+/X7uGblbtZl3qUIH8/Lm1Rg6GJdejWKBp/P/2P3CMEhkNWhjXvV6D2TCrly7T4KsHf359WrVqdfj1s2DAee+wxevbsyb59+wgNtZ4Qa9SoEZMnT7YrpvIBm/cf48tlO5m6Zg+ZuQU0rRnJv69szsA2sVTVcVyeJ6ho3FdeZoUVX462Xzq/l1L20uKrhNDQUNauXVvqtq+++kqX6VAVKje/kFl/7uPLZTtZmXKYoAA/rmhdixs616NtnSpeOxeXTwgIsZYXys2CsKgKuYS2X0p5Bi2+lHIDhzJz+XrFLj5fmkLasRzqRYXxxGVNuaZ9He3l8hYiEBhqPfGolPJpWnyVcPLkyb8sv/H4449z7bXXAnD99def7rbv06cPY8eOtSWj8h5b047z8eIdTP19Dzn5hXRPiOblwa3p0TgGPx3L5X2CwuFEGhQWWE9AOpmj7de8efOIiqqY3jelVPnct/ia9RjsX+/cc9ZsBf1fLnMX7bZXFc0Yw8qUw3ywaBvzNh8gOMCPwe3iuKlbfX1i0Vucqf0qzIf8k9a8X3KWxZe2X0p5DfctvpTyMoWFhjmb0vhg0TbW7DpC1bBA7uudwIgu9XUiVF/h52d9LywAf+9YvFwpdfbct/gq5xOeUp6ioNAwY91exi1IZkvaCepUC+XZAS24pn0dQoP0P2CvVFb7lbbBGvtVraHr8iil3Ir7Fl82KTlmol+/frz8stWQFh8zER0dzdy5c23JqDxDfkEh36+1iq7tBzNJqB7Bm8PacHmrWgT4+9kdT9klMBxyT1TIqctqv5RS7kOLrxIKCgpKfX/hwoWuDaI8VkGh4fu1e3hz3lZ2ZmTRrFYl3ru+HX1b1NRB9Mqa7yv7MOTnQoBzbzdr+6WUZ9DiSyknKSw0zPpzP/+Zk8S29Exa1K7ERyMS6d2sus7Ppf4nKNz6npfp9OJLKeUZtPhS6jwZY1iQdID/m72FjfuOkVA9Qnu61JkFhgJizfcVWtXuNEopG2jxpdR5WJd6hBdnbmLZ9kPUiwrj9Wsv4KoLYnW9RXVm4mf1fulkq0r5LC2+lDoHuw9lMXZ2EtP/2EtUeBDPDmjB8I51CdSB9MoRFTzZqlLKvWnxpdRZOJadxzvzk/n0tx34+wl392rEbRc1JDIk0O5oypMERQBpVu9XSCW70yilXEyLL6UcUFho+Hb1bsbOTiIjM5er28Xx4KVNqFk5xO5oyhOdGnSvxZdSPkmLrxL8/f1p1arV6dfTpk0jJSWFAQMG0KBBA0Dn+PI1q1IO8cwPG1m/5yjt61Xlk1EdaB1Xxe5YypP5+Rctsu38+b7S0tK4//77WbZsGVWrViUoKIhHHnmEqlWrajumlJvQ4quE0tZGS0lJoXv37syYMcOmVMoOGSdyeHHmZqasSaVW5RDeHNaGqy6ordNGKOcIioDMDDCF1iB8JzDGMHDgQEaOHMmECRMA2LlzJ9OnT6dq1arajinlJrT4UqqEwkLD1yt38epPSWTl5nNHz3ju6tWIsCD956KcKCgCMtMh7+T/bkOep/nz5xMUFMSYMWNOv1evXj3uvvtunWhVKTfi0MctEeknIkkikiwij5Wyva6ILBCR30VknYhc5vyornFqeY42bdowaNCg0+//+uuvp99/4YUXbEyoKtKGvUcZ/N4Snpz6J81qRTLr3u480q+pFl7K+U6P+3LerccNGzbQrl27M27Xdkwp91Du/ygi4g+MA/oAqcBKEZlujNlYbLd/ApOMMe+JSHNgJlD/fIK9suIVNh/afD6n+Jum1ZryaMdHy9yntNuOgHbXe7nsvALemLuVj37dTtWwQF6/9gIGtonVW4zqnDjcfuVmgUjRxKtlc6T9KunOO+9k8eLFBAUFMXbsWG3HlHITjnyc7wgkG2O2A4jIRGAAULz4MsCpR3YqA3udGVKpirR8ewaPfbeeHQczuTaxDk9c1ozKYTp1hHIBPz9rri8nadGiBVOmTDn9ety4cRw8eJDExESnXUMpdf4cKb5igd3FXqcCnUrs8zTws4jcDYQDvc832Nl+wlPqbB3PzuPlWZv5avku6lQL5avRnejWKNruWMoLONx+ZWXAkV0Q09Sh3q/y9OrViyeeeIL33nuP22+/3bpEVtZ5n1cp5VyOFF+l3XcxJV4PB8YbY14TkS7AFyLS0hhT+JcTidwG3AZQt27dc8mrlFMs2XaQh79dx76jJxl9YQMeuLSxjutSrhcUYX3PPeGU4ktEmDZtGvfffz+vvvoqMTExhIeH88orr5z3uZVSzuPI/zapQJ1ir+P4+23FW4B+AMaYpSISAkQDB4rvZIz5EPgQIDExsWQB5xZOnPj74NeePXvSs2dP14dRTncyt4BXftrM+CUpNIgO59sxXWlfTxc3VjbxDwK/QMjJhPAYp5yyVq1aTJw4sdRt2o4p5R4cKb5WAgki0gDYAwwDriuxzy7gEmC8iDQDQoB0ZwZV6nz9vuswD076g+0HMxnVtT6P9mtKaJCuq6dsJFK0yLbzJ1tVSrmvcosvY0y+iNwFzAb8gU+MMRtE5FlglTFmOvAg8JGI3I91S3KUMcYte7aU78kvKGTcgm28NX8rNSKDdWyXci9BEZB9BPJzICDY7jRKKRdwaJCLMWYm1vQRxd97qtjPG4Fuzo2m1Pnbc+Qk9038nZUphxnUNpZnBrSgki6CrdxJcLF1HrX4UsonuN0IY2OMV8ytpB1/9pu5fh+PTVlHoYE3rm3DwLaxdkdSXu6c2q+AUBB/69ZjWLWKCeYgbbeUcg23Kr5CQkLIyMggKirKowswYwwZGRmEhITYHcUnncwt4JkfNjBx5W4uqFOFt4a1oV6Uc5ZvUepMzrn9Oj3uK7PiwjlA2y2lXMetiq+4uDhSU1NJT/f8sfohISHExcXZHcPnbE8/wR1frSEp7Th39Izn/j6NCfR3zqLFSpXlvNqv7GPWuK/0AvCz7yEQbbeUcg23Kr4CAwNp0KCB3TGUh5qxbi+PTl5HUIAfn47qQM8m1e2OpHzIebVfu1fCx0PhmvHQbFC5uyulPJtbFV9KnYvc/EJenLmJ8UtSaFe3Cu9c147aVc5/wkqlXKZ2WwiKhO2LoIUWX0p5Oy2+lEdLO5bNmC9X8/uuI9xyYQMe7deUoAC9zag8jH8A1O8GO36xO4lSygW0+FIea/XOw4z5cjWZOfm8e307LmtVy+5ISp27Bj1gy09wNBUq67grpbyZdhEoj/T1il0M+3ApYUH+TL2jmxZe6oxEpJ+IJIlIsog8VsZ+V4uIEZFEV+Y7rWEP6/v2RbZcXinlOlp8KY+Sm1/Ik1PX8/h36+kSH830Oy+kSc1Iu2MpNyUi/sA4oD/QHBguIs1L2S8SuAdY7tqExcQ0g7Bo2KHFl1LeTosv5TEOZ+Zy48fL+Wr5Lm7vGc+nozpQOUxnq1dl6ggkG2O2G2NygYnAgFL2ew54Fch2Zbi/8PODBhdZPV862alSXk2LL+URtqefYPB7S/h91xHeuLYNj/Zrir+f507Eq1wmFthd7HVq0XuniUhboI4xZoYrg5WqYQ84sR8ObrU7iVKqAumAe+X2lmw7yO1friHAT/j6tk60r2fvEizKo5RWoZ/uVhIRP+B1YFS5JxK5DbgNoG7duk6KV0KDonFfOxZBTOOKuYZSynba86Xc2jcrdzHi4xVUjwxm2p3dtPBSZysVqFPsdRywt9jrSKAlsFBEUoDOwPTSBt0bYz40xiQaYxJjYmIqJm3V+lC5LmxfWDHnV0q5Be35Um7JGMN/5mzh7fnJdE+IZtz17agUouO71FlbCSSISANgDzAMuO7URmPMUSD61GsRWQg8ZIxZ5eKcpwJAw4tg0w9QaO9SQ0qpiqM9X8rt5BUU8uiUdbw9P5mhiXF8OqqDFl7qnBhj8oG7gNnAJmCSMWaDiDwrIlfZm+4MGvSE7KOwf53dSZRSFUR7vpRbycrN586v1rAgKZ17Lkng/t4JiOjAenXujDEzgZkl3nvqDPv2dEWmMjXobn3fvshadkgp5XW050u5jYMnchj+4TIWbUnnxUGteKBPYy28lO+JrAkxTXW+L6W8mBZfyi2kHs7imveXkpR2nA9uTOS6ThX0NJlSnqBBD9i5FPJz7E6ilKoAWnwp221LP8HQ95dy8EQOX43uRJ/mNeyOpJS9GvaA/JOQutLuJEqpCqDFl7LVhr1HGfr+UnILCpl4W2edSkIpgHrdQPxgxy92J1FKVQAtvpRtVu88xLAPlxEc4Mekf3ShRe3KdkdSyj2EVoFabXS+L6W8lBZfyha/JR/khv+uIDoimG9v70rDmAi7IynlXuJ7QeoqOHnY7iRKKSfT4ku53C9b0rl5/ErqRYUx6R9diK0SanckpdxPQh8wBbBtgd1JlFJOpsWXcqmFSQcY/fkq4mMi+PrWzsREBtsdSSn3FJsIIZUhea7dSZRSTqbFl3KZ+ZvTuO3z1SRUj2DCrZ2oGh5kdySl3Jd/gHXrMXkuFBbanUYp5URafCmXmLsxjX98sZomNSOZMLozVcK08FKqXI36wIk0SFtvdxKllBNp8aUq3PzNadz+1Wqa16rEl6M7UTlM12lUyiGNelvft86xN4dSyqm0+FIVavHWg4z5cg1Na1bi81s6UTlUCy+lHBZZA2pdoOO+lPIyWnypCrNixyFu/XwVDaPD+fzmjlp4KXUuGvWB3ct1ygmlvIgWX6pCrN19hJvHr6R2lRC+HK2D65U6Zwl9wBTqlBNKeREtvpTTbdh7lBEfL6daeBBfje5MdIROJ6HUOdMpJ5TyOlp8Kafann6CER+vICI4gAm3dqJm5RC7Iynl2fwDIP4SnXJCKS+ixZdymv1Hs7nx4xUAfHVrZ+KqhtmcSCkvkaBTTijlTbT4Uk5xJCuXEZ8s5+jJPD67uSMNosPtjqSU99ApJ5TyKlp8qfOWlZvPzeNXknIwiw9HtKdlbGW7IynlXSKqW1NOaPGllFfQ4kudl7yCQu74ag1rdx/hreFt6BofbXckpbxToz6QukKnnFDKC2jxpc6ZMYZHp6xjYVI6LwxqRb+WteyOpJT3SrhUp5xQykto8aXO2etztvDdmj080KcxwzvWtTuOUt4tLhFCq8KW2XYnUUqdJy2+1DmZtHI3b81PZmhiHHf3amR3HKW8n58/NO4HW36Cgjy70yilzoMWX+qs/bIlncenrqd7QjQvDGqFiNgdSSnf0PRyyD4CO5fYnUQpdR60+FJnZePeY9zx1RoSqkfw7vXtCPTXv0JKuUz8JRAQCpt/tDuJUuo86P+cymH7j2Zz8/iVRAQH8OlNHYgM0YWylXKpoDCI72UVX8bYnUYpdY60+FIOOZlbwOjPV3I8O49PRnWgVuVQuyMp5ZuaXg7HUmHfWruTKKXOkRZfqlyFhYaHvv2DDXuP8dbwtjSvXcnuSEr5rib9Qfz01qNSHkyLL1WuN+dt5cf1+3isX1MuaVbD7jhK+bawalCvmxZfSnkwh4ovEeknIkkikiwij51hn6EislFENojIBOfGVHaZsW4vb87bypB2cdx2UUO74yilwLr1eGAjZGyzO4lS6hyUW3yJiD8wDugPNAeGi0jzEvskAI8D3YwxLYD7KiCrcrF1qUd4cNIfJNaryouDW+qUEkq5iyaXWd+190spj+RIz1dHINkYs90YkwtMBAaU2OdWYJwx5jCAMeaAc2MqVztwLJtbP19FdEQw79/YnuAAf7sjKaVOqVoParbW4kspD+VI8RUL7C72OrXoveIaA41F5DcRWSYi/Uo7kYjcJiKrRGRVenr6uSVWFS4331os+9jJfD4akUh0RLDdkZRSJTW9AnYvhxP6WVcpT+NI8VXavaaSE8wEAAlAT2A48F8RqfK3g4z50BiTaIxJjImJOdusykWe/3Ejq3Ye5pWrW+uTjUq5q6aXAwaSZtmdRCl1lhwpvlKBOsVexwF7S9nne2NMnjFmB5CEVYwpD/Ptqt18vnQnt3ZvwFUX1LY7jlLqTGq0gCr19NajUh7IkeJrJZAgIg1EJAgYBkwvsc804GIAEYnGug253ZlBVcVbl3qEJ6f9Sdf4KB7t19TuOEqpsohAsyth+wLIOW53GqXUWSi3+DLG5AN3AbOBTcAkY8wGEXlWRK4q2m02kCEiG4EFwMPGmIyKCq2c7+CJHMZ8sZqYiGDeHt6WAF2zUSn31/QKKMiFLbPtTqKUOgsBjuxkjJkJzCzx3lPFfjbAA0VfysMUFBru+fp3DmbmMmVMV6J0gL1SnqFOJ4isDX9+B62utjuNUspB2r2heGPuFpZsy+D5AS1pFVfZ7jhKKUf5+UGLQZA8B04esTuNUspBWnz5uIVJB3h7fjJXt49jaIc65R+glHIvLYdYtx514L1SHkOLLx+298hJ7v9mLU1rRvLcgJZ2x1GqQpS3PJqIjBGR9SKyVkTD5UkeAAAgAElEQVQWl1zBw+3FtrOeevxzit1JlFIO8r3iKz8HMg/CoR2wfz3sWgZ7VkNupt3JXCqvoJC7JqwhN7+Qcde3IzRIZ7BX3seR5dGACcaYVsaYNsCrwH9cHPP8iFi9X9sXWm2bUsrtOTTg3uMd3gmbfrC+di/n73PEAghUrW/NnVOjBTToAXW7WGMqvNArszazZtcR3h7elviYCLvjKFVRTi+PBiAip5ZH23hqB2PMsWL7h1N6A+HeWg6Bxf+Bjd9Dh1vsTqOUKof3Fl/5ObDiQ1j/Lez7w3qvRivo/gBE1ICgCAiOsL7nZUHaRjiwwfqeNBMWvQKV60LroXDBMIj2njljf96wn/8u3sGILvW4UidSVd6ttOXROpXcSUTuxHpaOwjo5ZpoTlSjBUQ3gQ1TtfhSygN4Z/G1axlMvxsOboHYROjzrDUZYbWGZz6m2ZX/+znnhDV4dd1E69Pkr/8Hse2h233WflLaikueYe+Rkzw8eR0tYyvx5OXN7I6jVEVzZHk0jDHjgHEich3wT2Dk304kchtwG0DdunWdHPM8nbr1uPAlOLYPKtWyO5FSqgzedU8t5zjMfBg+6Qd5J+H6KXDrPOh2b9mFV0nBEXDBtXDjVHhgE1z6PGQfhUk3woc9YMvPYDzvzkRBoeG+b9aSV1DIW8PaEhyg47yU13NkebTiJgIDS9vg9mvTthwMGNg4ze4kSqlyeE/xteMXeLcLrPgIOt4GdyyFhN7nf97ImtD1brhjOQx835pLZ8I18ElfSFl8/ud3oXELklmx4xDPDWhJQx3npXxDucujiUjxMQWXA1tdmM95ohOgZmt96lEpD+AdxVfyPPhyCASGwi0/w2WvQnCkc6/hHwBthsNdq+CK1+HIbhh/OUy9HbIOOfdaFWBVyiHemLuFgW1qM7hdrN1xlHIJB5dHu0tENojIWqxxX3+75egxWg6B1JVwOMXuJEqpMnh+8bXjF5h4HcQ0gZtnQ52OFXu9gCBIvBnuWQPdH4L1k+CdDrB+stveijyalce9E9cSVzWM5wa2RDx4zJpSZ8sYM9MY09gYE2+MeaHovaeMMdOLfr7XGNPCGNPGGHOxMWaDvYnPQ4tB1vcNU+3NoZQqk2cXXzuXwoRroWoDuHEahFVz3bUDQ+GSf8Fti6BKXZhyi5XlaKrrMjjAGMPjU9eRdiybt4a3JTIk0O5ISqmKUrUexHWA9XrrUSl35rnFV+oq+OoaqFQbRnwP4dH25KjZEkbPhb4vQcqv8F5X2DTDniyl+HZ1KjPX7+fBS5vQpk4Vu+MopSpaq2sgbT3s/9PuJEqpM/DM4ittA3w5GMKjYOQPEFnD3jx+/tDlDrj9N6sX7pvrYeYj1lxjNtp9KItnpm+gc8Nq/OOis3jaUynluVpeDX6BsHaC3UmUUmfgecVXfg5MGQ0BIVbhVcmNJgmt1tAa8N/5DljxAXzcBzK22RKloNDwwKS1+Inw2tA2+PnpOC+lfEJ4FDTpD+u+gYI8u9MopUrhecXXwpfhwEa46m1rrJW7CQiGfi/B8IlwZBd8cJG15IeLffDLNlamHObZgS2IrRLq8usrpWzU5nrIOghbf7Y7iVKqFJ5VfKWuht/egDY3QOO+dqcpW5P+MGYxVG8Gk0bAghehsNAll/5zz1Fen7OFy1vVYmAbnVZCKZ/TqDeEV9dbj0q5Kc8pvvJOwrQxEFkb+r1odxrHVI6DUT9an0IXvWLNkJ9zokIvmZ1XwP3frKVqWBDP67QSSvkm/wBrlY4tP0HmQbvTKKVK8Jzia/7z1lqNA96GkMp2p3FcQDAMGGc9DZk0Ez6+tEInQBw7O4mtB04w9poLqBoeVGHXUUq5uTbXQ2E+rJtkdxKlVAmeUXztXApLx1mTm8b3sjvN2ROxnoa8YQoc2wMfXmwt/u1ky7dn8PHiHYzoUo8ejd1w7TmllOtUbwa12+mtR6XckNsXX3sObeWRn//B5qg60OdZu+Ocn/hecOt8CK0Kn10FG5y3AG5Wbj4PT15H3WphPNa/qdPOq5TyYG2us+b82veH3UmUUsW4ffEVlnWI2YGFzG3Z3/nrNdohKh5umQO128C3o2DJO05ZlujVn5LYdSiLsVe3Jiwo4PxzKqU8X6urwT9Ie7+UcjNuX3xVjetE2xrtWXh8u91RnCc8ypqVv9mV8POTMOtRKCw459Mt3ZbB+CUpjOpan04No5wYVCnl0UKrQtPLrXFf+bl2p1FKFXH74gvg4rq9SDqcxJ4Te+yO4jyBoXDNZ9DlLmtC1kkjIC/7rE+TmZPPI1P+oH5UGI/0a1IBQZVSHq3NDXDykPXko1LKLXhE8dWzTk8AFu5eaGsOp/Pzg74vQL+XYfMM+HIIZB89q1O8PGszqYdPMvaaC/R2o1Lq7+Ivhsha8PsXdidRShXxiOKrXqV6NKzc0PuKr1M63w6DP4Ldy2D8FXAi3aHDlmw7yBfLdnJT1wZ0qF+tgkMqpTySnz+0GwFb51ToNDdKKcd5RPEFVu/Xqv2rOJZ7zO4oFaP1UGtJooNb4ZNL4fDOMnc/mVvAY1PWUz8qjIf76u1GpVQZ2o0E8YNVn9qdRCmFBxVfF9e5mHyTz297frM7SsVJ6GMNxM/KgE/6QtrGM+76+twt7DqUxctDWhMa5O/CkEopj1M51lry7PcvID/H7jRK+TyPKb5aRbeiWkg1FuxaYHeUilW3E9z0kzX9xPjLYM+av+2yLvUI//11O8M71qWzPt2olHJEh9HWB7uN39udRCmf5zHFl7+fPz3ierB4z2LyCvLsjlOxajSHm3+y5jX77CrYueT0pryCQh6ZvI7oiGCdTFUp5bgGPSCqEaz8r91JlPJ5HlN8gXXr8XjecVYfWG13lIpXrYHVA1apFnwxGJLnAvDhL9vZvP84zw1sSeXQQJtDKqU8hp+ftUTb7uWwf73daZTyaR5VfHWu3Zlg/2Dvv/V4SuVYGDUTohvBhGHsXzaJN+dt5bJWNenboqbd6ZRSnqbNdRAQCis/tjuJUj7No4qv0IBQutTqwsLdCzFOWJLHI0TEwMgZmNptifnpHwwJWMLTV7WwO5VSyhOFVoVWQ6wZ789yTkGllPN4VPEFcHHdi9mbuZcth7fYHcV1QqvwbbO3WV7QlBd5m+pbJ9mdSCnlqRJvgbxM+OMbu5Mo5bM8rvi6KO4iBGHBbh+59QgcPJHD83N28m7sSxDfC6bfDSs+sjuWUsoTxbaD2u1g1cfWU9VKKZfzuOIrOjSa1jGtfar4euHHTZzMK+DpwYnI8K+hyWUw8yFY8o7d0ZRSnqjDaEjfDCmL7U6ilE/yuOILrNnuN2ZsZO+JvXZHqXC/JR9k6u97GNMjnkbVIyAgGIZ+Ds0Hws9Pwi9j7Y6olPI0LQdDWBQse9fuJEr5JI8svvrW6wvA7JTZNiepWNl5Bfxz2p/Uiwrjzosb/W+DfyAM+RhaD4P5z8OCF/X2gVLKcYGh0OFWSJppLWmmlHIpjyy+6lSqQ6voVszaMcvuKBXq/UXb2HEwk+cGtCQksMQSQv4BMPBdaHsjLHoF5j6tBZhSynEdRkNACCzV4QtKuZpHFl8A/Rv0Z9OhTew4usPuKBVie/oJ3l2wjSsvqM1FjWNK38nPH658y3p66bc3YPYTWoAppRwTEQMXDIe1X8OJdLvTKOVTPLb46lu/L4Lw046f7I7idMYY/vX9nwQH+vGvK5qVvbOfH1z+GnS63Rq/MfMhKCx0TVCllGfrcicU5MBKfXpaKVfy2OKrelh1EmsmMnPHTK+bcHXGun38lpzBw32bUD0ypPwDRKDfS9D1Hmvdthn3agGmlCpfdIL19PSKjyA3y+40SvkMjy2+wLr1mHIshaTDSXZHcZrMnHxe+HETLWpX4vpO9Rw/UAT6PAvdH4I1n8P3d0JhQcUFVUp5h653w8lD8McEu5Mo5TM8uvjqU7cPARLAzB0z7Y7iNG/N38r+Y9k8O6Al/n5ydgeLwCX/gouftBrS726DgvyKCaqU8g51u0Bse1j6rn5gU8pFHCq+RKSfiCSJSLKIPFbGfleLiBGRROdFPLMqIVXoUrsLP+34yStuPSYfOM7Hv+7gmvZxtK9X9dxP1OMR6P00/DkZJt8EBXnOiqiU8jYiVu/XoW2Q5N1PkCvlLsotvkTEHxgH9AeaA8NFpHkp+0UC9wDLnR2yLP0b9Gdf5j7+SP/DlZd1OmMM/56+gbAgfx7t3/T8T3jh/dD3Rdg0HSaNgPyc8z+nUso7Nb0SqtSFJW/bnUQpn+BIz1dHINkYs90YkwtMBAaUst9zwKtAthPzlatX3V4E+wd7/JxfM9fv57fkDB7q24ToiGDnnLTLnXDZ/1kTKU68DvJOOue8Sinv4h8AXe6C3csg5Te70yjl9RwpvmKB3cVepxa9d5qItAXqGGNmODGbQ8IDw7ko7iJmp8wmv9Azxzdl5uTz3IyNNK91loPsHdHxVmsusOR5MGEo5GY69/xKKe/QbgSEV4dFL9udRCmv50jxVdqo79MDrETED3gdeLDcE4ncJiKrRGRVerrzJvXr36A/GdkZrNy/0mnndKV3FiSz/1g2zw1scfaD7B3RfiQMet9aRPfLIZB9zPnXUEp5tsBQ6HYv7PgFdi61O41SXs2R4isVqFPsdRxQfEXrSKAlsFBEUoDOwPTSBt0bYz40xiQaYxJjYs4wa/s56B7bnfDAcI+89ZhyMJOPf93B4HaxtK9XreIudMEwaz3I1JXwxSA4eaTirqWU8kyJN0N4jPZ+KVXBHCm+VgIJItJARIKAYcD0UxuNMUeNMdHGmPrGmPrAMuAqY8yqCklcipCAEHrX7c3PO38mK8+zJgp8/sdNBPoLj/VzwiD78rQcDEM/h31/wGdXQmZGxV9TKeU5gsKsyZq3L4RdLn12SimfUm7xZYzJB+4CZgObgEnGmA0i8qyIXFXRAR01pPEQMvMymZ0y2+4oDvtlSzpzN6VxV68EqldyYCZ7Z2h6OQz/Gg5ugfGXw/H9rrmuUsozdLgFwqK190upCuTQPF/GmJnGmMbGmHhjzAtF7z1ljJleyr49XdnrdUqbmDY0qNyA77Z+5+pLn5O8gkKenbGRelFh3HxhfddePKEPXP8tHNkFn/aHI7vLP0YpD1XePIUi8oCIbBSRdSIyT0Sc/NSLhwkKt+b92jYfdnvmOFql3J1Hz3BfnIgwJGEIa9PXsu3INrvjlOvLZTtJPnCCJy9rRnCAv+sDNLgIRkyzbj1+2h8y3P/PTKmz5eA8hb8DicaY1sBkrClzfFuH0RAWpb1fSlUQrym+AK5oeAUBfgFu3/t1KDOX1+dsoXtCNH2a17AvSJ2OMHK6Nf3Ep5fBgc32ZVGqYpQ7T6ExZoEx5tRg0WVYDxX5tuAIa96v5LmQutruNEp5Ha8qvqJCo7i4zsX8sO0Hcgty7Y5zRq/9nERmbgFPXdEckQqYWuJs1G4DN80EDIy/DPautTePUs5V7jyFJdwCeN5j0xWh460QWhUWvmh3EqW8jlcVXwBDEoZwOOcwC3YvsDtKqTbtO8bXK3ZxY+d6JNSItDuOpXozuGkWBIZbT0HuXGJ3IqWcpcx5Cv+yo8gNQCIw9gzbK2SeQrcVHGktU5Y815r7SynlNF5XfHWu1Zla4bWYsmWK3VH+xhjDczM2Uik0kPt6J9gd56+i4uHmnyCyJnwxGLbOtTuRUs5Q3jyFAIhIb+BJrGlySl0ItaLmKXRrHW+DSnEw5ykoLLQ7jVJew+uKL38/fwY1GsTSfUtJPZ5qd5y/mL/5AEu2ZXDfJQlUCQuyO87fVY61esCiE+DrYbBhqt2JlDpfZc5TCKeXR/sAq/A6YENG9xUYCr2ehL2/w0ZtD5RyFq8rvgAGJQxCEKYlT7M7yml5BYW8MHMTDaPDub6zGz/JHh4No2ZAXCJMvhlWf2Z3IqXOmYPzFI4FIoBvRWStiPxtCh2f1vpaqN4C5j0L+e47llYpT+KVxVfN8Jp0i+3G1OSpFBQW2B0HgAnLd7E9PZMnLmtGoL+b/7GHVIYbvoP4XvDDPbD4dTClDpNRyu2VN0+hMaa3MaaGMaZN0ZfbTB7tFvz8oc8zcDgFVn9qdxqlvIKbVwHnbkjCEA5kHeC3vb/ZHYWjWXm8MXcLXeOjuKRZdbvjOCYoDIZ9DS2vhrlPw8//1AJMKV/VqDfU7w6LXoHsY3anUcrjeW3x1SOuB9Gh0UzYNMHuKLyzYCtHTubx5OXN7J9a4mwEBMHgj6xBt0vfge/vhIJ8u1MppVxNBPo8C1kZ8NubdqdRyuN5bfEV6B/IsCbD+G3vb7bOeL8zI5PxS1K4pn0cLWpXti3HOfPzg/6vQs/HYe1XMOlGyDtpdyqllKvFtoOWQ2DpODi2z+40Snk0ry2+AIY2GUqwfzBfbPzCtgwvz9pMoL8fD13axLYM500Eej4Gl/0fJM2CzwdC1iG7UymlXK3Xv6Aw3xp8r5Q6Z15dfFUNqcoVDa9gxvYZHMp2fbGwMuUQs/7cz5ge8VSvFOLy6ztdx1vhmvGwdw180k8X5FbK11RrAF3uhD8mwO4VdqdRymN5dfEFcGPzG8kpyGFS0iSXXtcYw4szN1GjUjC3dm/o0mtXqBYDrSchj++Hjy+FtA12J1JKudJFD0NkbfjxQXCTp8mV8jReX3zFV4mnW2w3Jm6e6NL1Hn/6cz+/7zrCA30aExrk77LrukSD7nDzLMDAJ/0hZbHdiZRSrhIcAX2fh/3rdOoJpc6R1xdfACOajSAjO4NZO1yzXm5eQSGv/LSZxjUiuLp9nfIP8EQ1WsAtc4qWIxoE6yfbnUgp5SotBltTT8x7DjIz7E6jlMfxieKrS+0uNKrSiC82foFxwVxVE5bvIiUji8f7N8Pfz4OmljhbVepY60HGdYApt8Cvr+lcYEr5AhG4bCzknoB5z9idRimP4xPFl4hwQ7MbSDqcxMr9Kyv0Wsez83hz3la6NIyiZxMfWHw3rBrcOBVaXWM9AfXDvToXmFK+oHoz6DQG1nwOe1bbnUYpj+ITxRfA5Q0vp2pwVT7f+HmFXueDRds5lJnLE5d52ISq5yMgGAZ9CN0fhDWfwdfXQs5xu1MppSpaj0chojr8+BAUFtqdRimP4TPFV0hACEObDGVR6iK2H91eIdfYfzSb/y7ezoA2tWkV54ETqp4PPz+45Cm48k3YtkCnolDKF4RUgkuft6afWfWx3WmU8hg+U3wBXNfsOkIDQvlw3YcVcv7X52yhsBDPnlD1fLUfBdd/C0d2wUe9IFVvRyjl1VpdA/GXwJx/w+GddqdRyiP4VPFVLaQaw5oMY9aOWU7v/dqadpxvV+/mxi71qFMtzKnn9jiNLrGehAwMhfGXwZ/f2Z1IKVVRRKwebxFrzKc+dKNUuXyq+AIY2WIkwf7BTu/9enV2EuFBAdx1cSOnntdjVW8Kt86HWm1g8k2w6FVtlJXyVlXqQO+nYfsCaw1YpVSZfK74igqNcnrv1+qdh5izMY1/9GhI1fAgp5zTK4RHw8jp0PpaWPCCNR1FbpbdqZRSFSHxFqjXDX56QhfeVqocPld8gXN7v4wxvDIriZjIYG6+sIET0nmZgGAY9AFc8m/r9uMnfXUgvlLeyM8PrnobCnLgxwe0p1upMvhk8eXM3q8FSQdYkXKIey5JICwowEkJvYwIdH8ArvsGDqfARxfDzqV2p1JKOVtUPFz8JCTNhD+n2J1GKbflk8UXOKf3q6DQ8OpPSdSPCmNYBy9dRsiZGveF0XMhOBI+uxJWf2Z3IqWUs3W5E2Lbw8yH9fajUmfgs8WXM3q/vl+7h837j/PgpU0I9PfZP8qzE9PEGojfoDv8cI/1dFR+jt2plFLO4ucPA9+HvJMwbYxOvqpUKXy6YjjV+/XBHx+c9bE5+QW89vMWWsZW4vJWtSognRcLrQrXT4YL74fV4+HT/nA01e5USilniWkM/V6C7Qth2Ti70yjldny6+IoKjWJ40+HM2jGLjRkbz+rYr5btYs+Rkzzaryl+3rx4dkXx87ceTb/2S0jfAh/0gB2/2J1KKeUs7UdB0ytg7jOw7w+70yjlVny6+AIY3Wo0VYKrMHblWIyDT+ecyMln3IJkusZH0T3BBxbPrkjNrrRuQ4ZVg88HwOI39CkppbyBiPX0Y3g0TBmt08woVYzPF1+RQZHc2eZOVqWtYv6u+Q4d88niHWRk5vJwXx9eRsiZYhpbBVizK2Huv+Hr4ZB1yO5USqnzFVbNmmrm4FaY/YTdaZRyGz5ffAEMaTyE+MrxvLb6NfIK8src93BmLh/9sp1Lm9egbd2qLkroA4Ij4ZrPoP+rkDzXug2p60Iq5fka9oBu98DqT2HTD3anUcotaPEFBPgF8FCHh9h9fDcTNk8oc9/3F23jRG4+D2mvl/OJQKd/wM2zrdef9IVl7+ttSKU83cX/hNrtYNodkLHN7jRK2U6LryIXxl5It9hufPDHBxzOPlzqPmnHshm/JIVBbWJpXCPSxQl9SFx7+MciaNQbfnoUvrlBb0Mq5ckCgmDoZ+AXABOvh5wTdidSylZafBXzcOLDZOVn8d4f75W6/a15Wyk0hvv7NHZxMh8UVg2Gfw2XvgBbZsP7F0LKb3anUkqdqyp14epP4GASTL9be7SVT9Piq5j4KvFc3fhqJiVNYvuRv068ujMjk29W7mZYh7rUqRZmU0IfIwJd74LRc6w1Ij+7Aha8BAX5didTSp2L+Iuh179gw3ew7F270yhlGy2+SrijzR2EBYTx4ooX/zL1xOtzthDgL9zdq5GN6XxU7bbwj1+g9bWw6GWrCDu80+5USqlzceH91vxfP/8LUhbbnUYpW2jxVUK1kGrc2+5elu9bzvRt0wFI2n+c7//Yy6iuDaheKcTmhD4qOBIGvQ+DPoT9f8J73WDt13rrQilPIwID34NqDeHbUXB0j92JlHI5Lb5KcU2Ta2hbvS1jV40l42QGr/2cRERQAGN6NLQ7mrrgWrj9N6jZylo3btIIyMywO5VS6myEVLJWt8jLhgnXQs5xuxMp5VJafJXCT/x4usvTZOVl8fii5/l5Yxq3XtSQKmFBdkdTAFXrwagZ0PsZSJoF73WBrXPsTqWUOhvVm8LQ8XBgI0y+WcdyKp+ixdcZNKzSkFtb3crStLlUrpbMTd3q2x1JFefnDxfeZ82MH1oNvroavr8Lso/anUwp5ahGveHy12Drz9a0MjqMQPkILb7K0DpyEAU51QmrPQ0//1y746jS1GoNty20BvGu/Qre7WLNkK+U8gyJN0HXe2Dlf/UJSOUztPg6A2MMb87ZQejRYWQWZPD272/bHUmdSWAI9H4abpkLQRHw5RBrHiHtBVPKM/R+BppdBbOfhE0z7E6jVIVzqPgSkX4ikiQiySLyWCnbHxCRjSKyTkTmiUg950d1rV+3HmRFyiHuvbAP1za5lgmbJvD7gd/tjqXKEtfempKi233w+5cwrhNs/tHuVEqp8vj5weAPIbY9TBkNO5fYnUipClVu8SUi/sA4oD/QHBguIs1L7PY7kGiMaQ1MBl51dlBXMsbwfz8nEVsllGs71OHedvdSO6I2j//6OMdz9akctxYYAn2egdFzISwKJl4H39wIx/fbnUwpVZbAULjuG6hSx3oCcq9+2FXey5Ger45AsjFmuzEmF5gIDCi+gzFmgTEmq+jlMiDOuTFda87GNNalHuXeSxIIDvAnIiiCVy56hf2Z+3lu6XN/mXxVuanY9tZYsEuespYneqcjrPoUCgvtTqaUOpPwaLhxGoRUgS8Gw4HNdidSqkI4UnzFAruLvU4teu9MbgFmnU8oOxUWGv4zZwsNo8MZ3O5/v+YFMRdwR5s7mJUy6/Tkq8rN+QdC9wfh9iXWwPwZ98Gn/axJWpVS7qlyLIyYZv37/WIgHNphdyKlnM6R4ktKea/Urh8RuQFIBMaeYfttIrJKRFalp6c7ntKFfly/j837j3Nv7wQC/P/6x3NLy1tIrJHIC8tfYOcxXd7GY0Q3gpE/wIBxkJEMH1wEPz2hEzv6CAfGrF4kImtEJF9ErrYjoyohKt7qAcvPhs8HwLG9didSyqkcKb5SgTrFXscBf/uXICK9gSeBq4wxOaWdyBjzoTEm0RiTGBMTcy55K1R+QSGvz91CkxqRXNm69t+2+/v581L3lwj0C+SRXx4hryDPhpTqnIhA2xvgrlXQ7kbrkfZ3OsCf3+ncQl7MwTGru4BRwATXplNlqtEcbpgCWYfgsyt1GSLlVRwpvlYCCSLSQESCgGHAX+67iUhb4AOswuuA82O6xvdr97I9PZP7+yTg51dahx/UDK/JM12fYWPGRp1+whOFVYMr34Rb5kB4DEy+CcZfobcivZcjY1ZTjDHrAB0Q6G5i28MNk+F4GnzaHw7rHQflHcotvowx+cBdwGxgEzDJGLNBRJ4VkauKdhsLRADfishaEfG4QVF5BYW8OW8rLWpXom+LmmXu27teb65pfA2fbviUuTt1Qk+PVKeDNSD/8v/AgQ3wQXeY8YD1KVt5k7Mds6rcTd3OMOJ7yD4Cn14GGdvsTqTUeXNoni9jzExjTGNjTLwx5oWi954yxkwv+rm3MaaGMaZN0ddVZZ/R/UxencquQ1k8eGljRErv9Sru0Y6P0jq6NU8sfoIth7e4IKFyOj9/6HAL3L0GOoyG1ePhrbaw7H3I1xUNvITDY1bLPZEHjFn1WnHtYeQMyD9pFWDpSXYnUuq86Az3QHZeAW/N20qbOlW4uEl1h44J9g/m9YtfJyIwgnvm38Ph7MMVnFJVmLBqcNlYGLMYal1grTH37v+3d+fRUVT5Ase/1ekknT0hIYSQQBZ2QhBlGQRBBQEdQB0VEBXXQXEfdVSO85z3ZnvzcHTkqKMwooyOjoKODiCIioCIggHZd0gIhARIIGTfOl3vj1uBtiUQpNNVnf59zqnTlerqW7++Sd/8+tatWz+DnYtkPJj/a0f8/E0AABy9SURBVNGY1Zaw+pjVNq9jNtzxCegulYAVbTE7IiF+Mkm+gPe+O0hRWS1PjO7Rol6vJonhicy6YhbF1cU8vupxGlwyAN+vdeitTm9MWaAuc3//VnhjLBSsNzsy8dOdc8yq8COJveDOpWB3qDFg+5abHZEQP0nAJ1819Y28snI/g9LbMbRr/Hm/vm/7vvz20t+ScySHmd/59cT+AtRVkd1Hw31rYNyLcCIXXh+pZskvltPL/qYlY1Y1TRuoaVoBcBMwW9O07eZFLM4poSvc8znEpcG7E2GTXKQq/I/d7ADM9vbaAxRX1PHyzf3Pq9fL3YTMCew+sZu3drxFt7huTOwx0ctRCp8LssOAO6HvjfDNy/Dty7BrMfSbApc/rW6BIvyCrutLgCUe2551W8/Bz+/KEXCik+HOJepL0cfT1TQUw59QX56E8AMB3fNVWefktVW5XNYtgcEZ59/r5e5Xl/yKoZ2G8sd1f2T5QekKbzNCo+CKGfDIZhg8HbYugJcuhqVPq8vfhRDmcMTALR9A9iRY8QdY9AjI3IvCTwR08vWPbw5woqqex67qfsFl2W12XhjxAlnxWTy56klyjuR4IUJhGREJMPZP8PD3qrH/bg7M6gfLnoFKv53aTgj/Zg+B62er24h9/w81G36lXIkqrC9gk6+ymgZmr9rPyJ6J9O8c55Uyw4PDeWXkK6RGpfLQlw+x4/gOr5QrLCQmBa59GR7MgT7Xq5nyX8w2kjBp9IXwOU2Dkc/CL/4OhzfAnMuhcKPZUQlxVgGbfM39Oo/yWie/8kKvl7tYRyyzr5pNTEgM07+YzoGyA14tX1hEfCZc/6q6XVGf64wkrC8sfQrKCsyOTojAkz0R7lqm1t8YC5vfNzceIc4iIJOv0qp63vg6j7F9ksjqFOP18jtEdGD2VbMBmPb5NI5UHfH6MYRFxGfC9a/BAznQ9wbIeV2djvzPA1Cyz+zohAgsyRepO1d0GgAfTVNjM2XCZGFBAZl8zVmdS1W993u93KXFpPHqqFepqK/gjk/voKBCekPatISucO0r8PAmGHA3bP0AXh6grsY6JOP/hPCZyPYw9WN1gcy6V+GN0XJLImE5AZd8FVfUMW/NAcZnJ9MjKapVj9U7vjevj36dyoZKbv/0dvLK8lr1eMICYlPhmpnw6Da47DHIWwVzR8HcMbBzMbjk3s1CtLqgYLj6zzDpn2quvtkj1BciISwi4JKv11btp87ZyCOjuvnkeH0S+jB39FycLid3fHqH3AcyUES2V4OAf7UDxv4fVBTC+7eo3rB1c6CuwuwIhWj7eo1Xtw3r0Bs+vBs+fgDqq8yOSojASr6Oltfyz7X5XN8/hcz2kT47bo92PXhz7JvYNTt3LbuL7SUygXbACI2En90HD22EG9+EsDhY+mt4oTd8OkNOhwjR2mI7wx1L4LInYNM78OpQyP/G7KhEgAuo5OvlL/fR6NJ5ZKRver3cZcRkMO/qeUQGR3LPZ/fwTaF8+ANKkB2yfgG/XA73LIfuY9RcYS9dAu9MhD3LwNVodpRCtE1Bdhj5X3DH4tM35v50BtRXmx2ZCFABk3wVlFbzXs5BJg5MpXN8uCkxpEalMm/sPDpGduT+L+5n/u75psQhTJYyAG54XY0LG/5rNSfRuxNh1kXw1V9k0lYhWkvaMJj+DQy8R00P89owOLjO7KhEAAqY5Oul5fvQNI2HruxqahxJEUm8NfYthiQP4fdrf89zOc/RKD0egSm6I1z5DDy2A26aB+3S4Mvfq1OS86fCvi+kN0wIbwuNhJ//BaYuBFcDvDEGFj8GNaVmRyYCSEAkX3klVXzwfQG3DO5Mx5gws8MhMiSSl658iSk9p/DWjrd4dOWjVDdI93fACgpWs+XfvkhN2jpoGuSthn/eoGbPX/EnKM03O0oh2paMEaoXbPB9sOFNeHkgbH4PdN3syEQACIjka9YXewgJsjH98kyzQznFbrMzY/AMZgyawVcFX3Hb0ttkNnwBCd3UPSQf36V6w9r3gFUzYVY2zBsHG9+RKyWF8JbQKDUlxbSVENsFPrpXfc6O7TQ7MtHGtfnka8/RCv6zuZDbL00jMcphdjg/MqXXFP428m8cqz7GpMWTWJK7xOyQhBXYQ1Vv2G3/hke3whXPQPlh+M/98Fw3+PCX6rRko9PsSIXwfx37wd2fw7gX4eg2dUXkJ49DVYnZkYk2qs0nX3/9fA8RIXbuHZ5hdijNGtppKAvGL6BHux48tfopfvft76h11podlrCK2FQY8SQ89L36B3HRzbB3mTot+XwP9U8i/1uZwFWIC2GzwYA74aEN6nH9m+oimNUvQIO0x8K72nTyte1wGUu3HeHuYenERYSYHc5ZJUUkMXfMXO7KuosFexZw65JbyS3LNTssYSWaBqmDYNxf4Ym9MOkdSL9MnYp8c6w6NfnZb6BgvYxbEeKnikiAnz8P93+rro5c/j9qcuRN/5ILYITXaLpJjfSAAQP09evXt+ox7pqXw4b8UlY/dQXRjuBWPZY3fVXwFc98/Qw1zhoevOhBbut9G0G2ILPDElZVVwG7l8LWBbB/hbqCKzoFel+rlpSB6lu9BWiatkHX9QFmx3GhfNF+CYvIXaW+1BzZAvHd4PKn1ZAAaZMDjjfbL2u0yK1gQ/4Jvtx1jGnDM/wq8QIYnjKcf0/4N5cmX8rzG55n6qdTpRdMNC80CrInwi0L4Nd74brXIKkv5Pxd3VT4hZ6w6BHY+zk468yOVgj/kjECpq2CiW+Bza5uU/TqUNj+sZzqFz9Zm+z50nWdyXPWsr+4kq+evILwEHurHKe16brO0ryl/Om7P1HTUMMD/R9gau+p2G3++X6Ej9WWqZnzd32iBufXV0JIJHQdCd2vhm5XqVMsPiQ9X8KvuVyw4yNY+Wco2QMJ3eHSh9WXH3uo2dGJVubN9qtNJl+r9xZz29zv+O/xvbljaHqrHMOXSmpK+MPaP7D84HK6xnZlxqAZDOo4yOywhD9pqIW8r2DXYtjzKVQeBTQ12363MSoRS8pu9dOTknyJNsHVCNs/gq9fhKNbIaoj/Gw6XHInOKLNjk60Ekm+zkLXda59ZQ3HK+v58okRhNrbxnl5XddZcWgFM3NmcrjyMKO7jOaJAU/QMbKj2aEJf+NywZHNqldszzIo/F5tj2gPmVdC11GQcQVEtvf6oSX5Em2KrsP+L2HNi+rLTWg0XHQLDPolxFtnXknhHZJ8ncWn245w3z83MPPGbCYOSPV6+WarddYyb/s85m6dC8BdWXcxtc9UIoIjTI5M+K2Ko7B/Oexbrv6R1JxQ25P6QsblkH45dBkCIRf+NybJl2izDn+v7he5/WN10UvXUTDoXvVokQtexIWR5KsZjS6dsS9+RaOu89mjw7EHtd0/+MLKQp5f/zyf5X9GXGgcd/e9m0k9JuGwW28iWeFHXC4o2qQSsbxVcGgdNNaDLVhNc5F2mbr8PmUgBJ//35okX6LNqzgKG+bB+jeg8oiaOb//rXDRFIhJMTs6cQEk+WrGv78v4LH5m3llysX8PDswTsdtLd7KSxtf4tuib0kMS+TefvdyfdfrCQ7yrys8hUXVV8PBbyF3pUrGjmwF3QVBISoBSxsGw3+t7k/ZApJ8iYDhrIddi2DDP9RnB02d1r/4NuhxjQzQ90OSfJ1BvdPFyBdWEu0IZtGDw7DZNK+V7Q9yjuTw0saX2HhsI0kRSdzW6zZu6H6DnI4U3lVzEg6uhQOrIX8NVJ+AR7e0+OWSfImAVHpATYa86R11m7DQGOg9AfrepL7AyJxhfkGSrzN4e20+//XxNt68cyBX9Ej0Wrn+RNd11hSuYe7Wuaw/up6okCgm95jMlF5TSAjz7ZQCIkA468He8rtHSPIlApqrUfUib10AOxep6V8ikyDrF9D7OktNiCx+TJIvD9X1TkY8t5K0+HDm3zsETQusXq8z2VK8hXnb5/FF/hcE24IZkzaGiT0m0q99P6kfYRpJvoQwNNSoaV+2fgB7P1NjKyM7QM9x0Gu86hGT4SOW4s32q03M1vnG13kUV9Tx6i0XS2JhyG6fzQuXv0B+eT5v73ibRfsXsSh3Ed3jujOx+0TGZY6TU5JCCGGW4DB1m6I+16sJkfd+DjsXwuZ/wfq54IiBzJHQfQx0vQoi4s2OWHiR3/d8naiqZ/jMFQzJjOfvU/3+C3WrqWqoYkneEubvns+uE7sIs4cxsvNIxmWMY3DHwTJrvvAJ6fkS4hwaatTVxruXqh6xqmOcmhC5aQ6+ThdLr5gJpOfLzSsr9lFd7+TJMT3MDsXSIoIjuKn7TdzY7Ua2lWzjw70f8ln+ZyzOXUy8I56r06/mmvRryErIkt5DIYQwS3AY9BqnlqapX/Z+piZEXvlnWPm/EBIF6ZepefjShkH7XjJWzM/4dc/XoRPVjHx+Fdf1T2bmjf28FFngqGusY3XBahbnLmZVwSqcLicdwjswsvNIRnUZRf/E/tIjJrxKer6EuADVJ9RM+rkrIXeFuooSICwOOl8KXYwlqa/0jLUC6fky/PXzPWgaPDqqu9mh+KXQoFBGdRnFqC6jKKsrY1XBKr7I/4IP937Iu7veJS40jmGdhjGs0zCGJA8hzhFndshCCBG4wttBn+vUAlCar6Z8yV8DB9bA7k/UdnsYJPeH1IGQMkidsoxKMi9u8SN+m3ztLCrno02HmTY8g+TYMLPD8XsxoTFMyJzAhMwJVDdUs6ZwDcsPLmf14dUsyl2EhkZWQhaXJl/K4I6DyW6fTWiQTBIohBCmieuiloumqJ/LC9WkyAXr4dB38O3fwDVLPReVDMkXqaQsuT8kZUNUB/NiD3B+e9rxzje/Y0N+KaufvJKYcOlebS2NrkZ2ntjJ6sOrWXN4DVtLtuLSXYTYQujbvi8DkwZySYdL6JvQV66eFOckpx2F8KGGWijaDIXfQ+FGtZTsBYz/+xHtoUOWOk2Z1Bfa94SE7j/p1mGBIOBPO67NPc6K3cU8fXVPSbxaWZAtiKyELLISspjebzrl9eVsPLqRnCM55BzNYc6WObh0FzbNRtfYrmS3z6Zf+35kxWeRFpMmY8aEEMIswQ7oPFgtTWrL4cgWdauwI9vU+rrX1DxjAJoN2mWoRKx9T0joppb4buCINud9tEF+95/R5dL54yc76Rjj4I5L08wOJ+BEh0QzInUEI1JHAFBeX87W4q1sLt7MluItLMtbxgd7PgDAEeSge7vu9GrXi17tetE9rjuZsZmEB4eb+RaEECJwOaLVFZJpw05va2yA4/vg2E61FO+EY7vUdBd64+n9IpMgPhPapasErWmJ7QJhsb5/L37M75KvjzcdZuvhMv46qR+OYLkfltmiQ6IZ2mkoQzsNBcClu8gry2PH8R3sPLGTncd3sjh3Me/vfv/Ua1IiU+ga15WusV1Jj0knLTqNtJg0okPkW5UQQvhcUDAk9lKLO2c9lOapU5Ule1SCdiJXTQhbefSH+zpiVBIW10U9xqRCTCeISVHr4fEg0xid4lfJV019I88t2012SgzX9utkdjjiDGyajczYTDJjMxmfOR5QCVlBRQF7T+5lX+m+U49fF3yNU3eeem07Rzu6RHchNSqVlMgUUqJSSI1KJTkymYSwBGyazGMjhBA+Yw+B9j3U4qmuQk11cXw/nDwIJ/PV1ZfFu1Vy5qz94f5BoRDdUQ38j05W65FJ6irMyA6nH0OjAiJJ86vk6/XVuRSV1TJrcn9strb/y2krbJqNztGd6RzdmZGdR57a3uBqoKCigANlBzhQfoC8sjwOVRxiXdE6FlYv/EEZdpudpPAkkiOTSYpIokN4BzqEdyAxPJHEiEQSwxKJc8TJGDMhhPCF0KjTA/U96bqak6zsEJQVqMfyw1BeBBVFcHg97CyCxrofv9buUBcC/GCJVz1n4QkQkaDWw+LU1BuhMX45wazf/Kc6Vl7Lq6v2M7ZPEoPS25kdjvCCYFsw6THppMek/+i5usY6DlcepqCigMLKQoqqiiiqLKKwqpC1RWspqSnBpbt+8BoNjThHHAlhCSSEJdDO0Y44RxztHO3UemgcsY5YYkJjiA2NJTokWpI1IYTwNk1TCVNEvJre4kx0HWpPQsVRqDwClceg4ghUFaul8hhUFKoLA6pLTl8Q8KNj2cARq5KxsFhj3Xh0xLgt0SpRC40y1qPUEhJlSvLWov88mqaNBWYBQcDruq7/2eP5UOAt4BLgODBJ1/UD3gz0hc/30NDo4umre3qzWGFRoUGhZMRkkBGTccbnG12NHK89zrHqYxytOkpxTTElNSWU1JRwvOY4JTUlHCg7QGldKTXOmmaPExkcSXRINNGh0USFRBEdEk1kcCRRIVFEhkQSGayWiOAIwoPDiQiOUIs9grDgMMLt4TjsDjklamFWaL+EEB40zUiY4iDxHP/XdV2d5qwugarjUHNC9azVnICaUrVeexJqTqqfS/PUel05uJxnL3vyv6DnNd57Xy10zuRL07Qg4BXgKqAAyNE0baGu6zvcdrsbKNV1vaumaZOB/wMmeSvInUXlvL/+EHcNTSctQeaSEmoKjMTwRBLDE8lKyDrrvjXOGkprSymtLeVk3UlO1p2krK5MLfVlVNRXUF5XTnl9Ofnl+VTUV1DZUElVQ1WL4wmzh51aHEEOHHZjCXIQGhRKqD301LrD7iAkKIQQWwihQaEEBwUTEqTWQ2whBAcFE2xT24Jtwae22TX7qefsNvvpRbOf2ib35fwhK7RfQogLpGmqt8oRra6ubCldh4ZqNb1G7UmVwNWVq5+b1jv0br24z6IlPV+DgH26rucCaJr2HnAt4N54XQv8t7H+AfCypmma7oUZXHVdTS0RExbMw1d2u9DiRAAKs4cRFhlGcmTyeb3Opbuoaqiisl4lYlXOKqoaqqhuqKaqoYoaZw01zhqqndVUN1RT66yltrGWGmcNtU71WF5fTm1jLbXOWuqcddQ21tLgaqDWWYuO9yc4DtKC1GILwq6p5CzIFoRNs2HX1HrTPjab2mbTbOpnzabWjf3dt9k0GzbUcxraqW0RwRE8O+RZr78PLzK1/QK4/r0nKazZ742ihBBelnzwEB9Nnunz47Yk+eoEHHL7uQAY3Nw+uq47NU0rA+KBEvedNE2bBkwD6Ny5c4sC3FJQxtf7Svjt+N4yoarwKZtmIyokiqiQKK+Xres6Tt1JfWM9dY111DfW09DYQL2rXq27GmhwNZxeb2ygQVePTpeTBpd6dLqcOHXn6XWXk0a98UePja5G9Whsc+muH2xz6a7Tj65GGvQGXLrr1HYdXT3qpx9dussf7mpgavslhBBn0pLk60znMTy/EbZkH3RdnwPMAXV7jhYcm36psSy4bwj9UmQCN9F2aJpGsKZOIfpBAuPPTG2/AFO+VQshrK0lo4QLgFS3n1OAwub20TTNDsQAJ7wRIMDAtHaE2GVAsxDivJnefgkhhKeWZDQ5QDdN09I1TQsBJgMLPfZZCNxurN8IfOmt8RJCCHEBpP0SQljOOU87GmMgHgSWoS7VfkPX9e2apv0OWK/r+kJgLvC2pmn7UN8YJ7dm0EII0RLSfgkhrKhF83zpur4EWOKx7Vm39VrgJu+GJoQQF07aLyGE1chAKiGEEEIIH5LkSwghhBDChyT5EkIIIYTwIUm+hBBCCCF8SJIvIYQQQggfkuRLCCGEEMKHJPkSQgghhPAhzayJnDVNKwbyz+MlCXjc6NYiJK7zZ9XYrBoXWDe2842ri67r7VsrGF+R9ssnrBqbVeMC68Zm1bjg/GLzWvtlWvJ1vjRNW6/r+gCz4/AkcZ0/q8Zm1bjAurFZNS6rsWo9WTUusG5sVo0LrBubVeMC82KT045CCCGEED4kyZcQQgghhA/5U/I1x+wAmiFxnT+rxmbVuMC6sVk1Lquxaj1ZNS6wbmxWjQusG5tV4wKTYvObMV9CCCGEEG2BP/V8CSGEEEL4P13XvboAY4HdwD7g6TM8Hwq8bzy/Dkhze26GsX03MOZcZQLpRhl7jTJDmjuGWxnHUZeVeh7DlLiM7bcDLqAOKAReM6HOhgPfA07gRo8yjhp1the43UJx6UZ9bQIW+jiux4AdwBZgOeoS5KZyioB6o85O1ZeZsVmkzu4DthrH/hrofa5j+Hpp7j2Z3U4g7Vebab/cytHd6myhCXV2tnbCzDo7W9tqdp15rQ3zdsMVBOwHMoAQYLN7cMY+92N8OIHJwPvGem9j/1Djje83ymu2TGA+MNlYfw2Y3swx5htljDZ+oVuAq5qOYWJc7xvl5Bu/MDPrLA3IBt4CbnQroz+QC2wDfmasx1kgrgyg0sT6ugIIN9anu/0u84CDQAejzg411ZeJsc3HGnUW7Xa8CcCnZzuGN9smab+k/cIC7ZfH31mVyXXWXDthdp0117Zaoc681oZ5u/EaAixz+3kGMMNjn2XAEGPdjsquNc99m/ZrrkzjNSWA3fPYZzjGSWPbDI/FfT8z4ioxnl8FbDOzztz2nYdqvIYY5d0MzHYrYzZws9lxGdsqza4vY3t/YI3x/GZgtls565rqy8TYtliwzm4Glp7pb949Fl8uzb0nab+k/fJWXO7PA5VWqDPPdsIqdeYWV1PbarU6u6A2zNtjvjqhvuk3KTC2nXEfXdedQBkQf5bXNrc9HjhplOF5LM9j1ALH3LY37Xu21/girjKgJ+o0VbqmaRuBqYDnhG++iM1TUxlm11lzcQE4gF8Cj2iadt2Z9vFRXHcDS43nat3KKkB1k7u/xozYtmOROtM07QFN0/YDM4GHPY99jvfS2qT9kvarteNyP65D07T1qM/osGb2MaOdsFKdubetlqgzb7Vh3k6+tDNs01u4j7e2N3cMz+26x6OZcVUDnXVd7w+8C4zWNC26mbhbKzZPmsej+76+rLPm4gLoDPwW9S3jRU3TMn0dl6Zpt6L+2TzXzP6erzEjtoVum02tM13XX9F1PRN4CvjNOY7ta9J+SfvV2nG5H7ezrmZW/xtwmdtn0aexnaGdsESdNdO2ml5n3mrDvJ18FQCpbj+noAbGnXEfTdPsQAxw4iyvbW57CRBrlOF5LM9jOIBEt+1N+57tNb6IKwbYBXTSdf24sY8TKAa6+7jOPDWVYXadNRcXuq43xbMLWInqovZZXJqmjQKeASboul5nlONwKysF9aF0fy8+jw04YJU6c/Me0NTz1pJ2wxek/ZL2q7XjOnVc47MIapzQXk5/Fn0WWzPthOl11kzbaok6c3NhbdjZzkme74I6z5qLGnDWNICtj8c+D+AxmNRY78MPB6zlogbENVsmsIAfDoq7v5ljLDDKcB+wOrrpGCbGNd8o5wCQaZSzE3WlSTtf1pnbseahxkw0lXERahD5NtS577ym2EyOqx8QaZQxFPVB7O3D32V/1KDKbh5//wdQ3c+JRp0V+Pp36RmbherMva7GA+vPdgxvtk3Sfkn7hQXaL492oodRzjbURQu9fVlnNN9OmFpnnnFZrM681oa1RgN2DbDHqLxnjG2/Q2WwoL7FLUBdkvkdkOH22meM1+0Grj5bmcb2DKOMfUaZoc0dw62M48ayG3jH7LiM7X9EXT5bh8qWx5tQZwNRiUKVUT/b3co4ZmzbB/zHInEdNOqrCHXp70Ifx/UF6p/MJtymbTDKOQI0GPHeacLv8kexWaTOZhm/v03ACtwSm+aO4evlTO/JhN+ftF9tuP0ytj/uVmdHUGOYrNJOmF1nzbWtVqgzr7VhMsO9EEIIIYQPyQz3QgghhBA+JMmXEEIIIYQPSfIlhBBCCOFDknwJIYQQQviQJF9CCCGEED4kyZcQQgghhA9J8iWEEEII4UOSfAkhhBBC+ND/AwC0w9KHS+/NAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x360 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize = (10,5))\n",
    "plt.subplot(121)\n",
    "plt.plot(result.times, np.real([a[0,0] for a in result.states]), label = 'GG')\n",
    "plt.plot(result.times, np.real([a[1,1] for a in result.states]), label = 'EE')\n",
    "plt.plot(result.times, np.real([a[2,2] for a in result.states]), label = 'FF')\n",
    "plt.title(\"Populations\")\n",
    "plt.legend(loc = 2)\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.plot(result.times, np.abs([a[0,1] for a in result.states]), label = 'GE')\n",
    "plt.plot(result.times, np.abs([a[1,2] for a in result.states]), label = 'EF')\n",
    "plt.plot(result.times, np.real([a[0,2] for a in result.states]), label = 'GF')\n",
    "plt.title(\"Coherences\")\n",
    "plt.legend(loc = 2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $T_1$ and $T_2$ decay of a qutrit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For qutrit $T_2$ decay, we introduce the operators $D_0$ and $D_1$ below, which are analogues of the $\\sigma_z$ operator in a qubit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "gamma1GE = 1./(100000e-6) # 1/s\n",
    "gamma1EF = 1./(5e-6) # 1/s\n",
    "\n",
    "C0 = np.sqrt(gamma1GE)*qt.Qobj(np.array([[0,1,0],[0,0,0],[0,0,0]]))\n",
    "C1 = np.sqrt(gamma1EF)*qt.Qobj(np.array([[0,0,0],[0,0,1],[0,0,0]]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "gammaphiGE = 1/(100000e-6) # 1/s\n",
    "gammaphiEF = 1/(100000e-6) # 1/s\n",
    "\n",
    "D0 = np.sqrt(gammaphiGE) * qt.Qobj(np.diag([1,-1,0]))\n",
    "D1 = np.sqrt(gammaphiEF) * qt.Qobj(np.diag([0,1,-1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "psi0 = qt.Qobj(np.array([[0.5, 0, 0.5],\n",
    "                        [0, 0, 0.5],\n",
    "                        [0.5, 0, 0.5]]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = qt.mesolve(H, psi0, times, [C0, C1, D0, D1], [])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1208232e8>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl8AAAE/CAYAAAB4o6baAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd4lFXexvHvmUmlBULooYMoAgIGFBVFFAUbYAUb+u6KKPa1sBbWde1lLWvfXV0VEVQUUbGgogKKEhApIlKkhBJKICF9MnPePybBiIEMyUyeKffnunINM0+74y6H35znPOcYay0iIiIiUjdcTgcQERERiSUqvkRERETqkIovERERkTqk4ktERESkDqn4EhEREalDKr5ERERE6pCKLwkZY8wgY0xWLY5/3hhzVzAziUhsqW07JBIKKr5iiDFmnTGmyBiTb4zJNsa8bIxp4HQuAGPMZcaYuZU/s9aOs9b+w6lMIhJejDEXGmMyy9uwLcaYj4wxxzmdS+RgqfiKPWdaaxsAfYF+wJ0O5xERqZYx5ibgCeB+oAXQDngWGF6HGeLq6loS3VR8xShr7SbgI6CHMaa1MWaGMSbHGLPaGHNFxX7GmLuNMW8bY6YaY/YYYxYZY46otN0aY7pUev8/Y8y9VV3TGDPBGLOm/Dw/GWNGln9+GPA8MKD8G+3uqs5ljLmiPF9Oed7W++QYZ4xZZYzZZYx5xhhjyrd1McZ8ZYzJNcbsMMZMDdZ/RxEJPWNMCnAPMN5a+461tsBa67HWvm+tvcUYk2iMecIYs7n85wljTOI+5/iLMWZbeY/Z5ZU+TzTGPGqM2VB+R+B5Y0xy+bZBxpgsY8xtxpitwMvln59hjFlsjNltjPnGGNOr0vnWGWNuNsYsKW9zphpjkiptH15+bF55ezi04nc0xvy3PN8mY8y9xhh3+Ta1YVFGxVeMMsa0BU4DfgDeALKA1sC5wP3GmJMq7T4ceAtIBSYD040x8TW47BpgIJAC/B2YZIxpZa1dAYwDvrXWNrDWNq4i72DgAeB8oBWwHpiyz25n4O/NO6J8v1PLP/8H8CnQBEgH/lWD7CLinAFAEvDufrbfARwN9Mb/978/v+/Vb4m/3WkD/Al4xhjTpHzbQ8Ah5cd2Kd9n4j7HpgLtgbHGmL7AS8CVQFPgBWDGPsXe+cBQoCPQC7gMwBjTH3gVuAVoDBwPrCs/5hWgrDxDH+AU4M/l29SGRRkVX7FnennP0lzgK+BF4DjgNmttsbV2MfAf4JJKxyy01r5trfUA/8TfCB59sBe21r5lrd1srfVZa6cCq/A3koG4CHjJWrvIWlsC/BV/T1mHSvs8aK3dba3dAMzG35gCePA3nK3Lf8ffjS0TkbDXFNhhrS3bz/aLgHustdustdvxf7mr3IZ5yrd7rLUzgXygW3nv+BXAjdbaHGvtHvy3NUdVOtYH/M1aW2KtLSrf/wVr7XfWWq+19hWghN+3iU+Vt3U5wPv81hb9CX87Nqu8Hdxkrf3ZGNMCGAbcUN6rtw14vFIOtWFRRsVX7BlhrW1srW1vrb0af29XRaNTYT3+b38VNlb8wVrr47desoNijLm0Ulf9bqAHkBbg4a3Lc1XkyAd27pNza6U/FwIVDxPcChjge2PMcmPM/x1sdhFx1E4g7QBjrn7XPpT/uXIbtXOfwq2ifWgG1AMWVmqXPi7/vMJ2a21xpfftgb9U7F9+TNt9rre/tqgt/jsA+2oPxANbKp3zBaB5+Xa1YVFGgwdlM5BqjGlYqQBrB2yqtE/bij8YY1z4u703l39UiL/xqtASf3H2O8aY9sC/gZPw3170GmMW429QAGwAOdtXOl99/N+GN+33iIoTW7sV/7dVjP/JqM+MMV9ba1dXd6yIhIVvgWJgBPB2Fdsr2ofl5e/b8VsbdSA7gCLg8PJxsFXZt23aCNxnrb0vgPPvayPQeT+flwBpVfXuqQ2LPur5inHW2o3AN8ADxpik8oGjfwJer7TbkcaYs8u/dd6Av5GYX75tMXChMcZdPnD0hP1cqj7+Rmw7QPmA1x6VtmcD6caYhP0cPxm43BjTu3xsxf3Ad9baddX9jsaY84wx6eVvd5Xn8FZ3nIiEB2ttLv5xWM8YY0YYY+oZY+KNMcOMMQ/jH7d6pzGmmTEmrXzfSQGc14f/S+HjxpjmAMaYNsaYUw9w2L+BccaYo4xffWPM6caYhgH8Kv/F346dZIxxlV/rUGvtFvxjuh4zxjQq39bZGHNCeSa1YVFGxZcAjAY64P+m+C7+8Q2zKm1/D7gA/1/6S4Czy8d/AVwPnAnsxj/uYnpVF7DW/gQ8hv8bbDbQE5hXaZcv8H9r3WqM2VHF8Z8DdwHTgC34vz2O2ne//egHfGeMyQdmANdba38N8FgRCQPW2n8CN+EfSL8df2/RNfjbnHuBTGAJsBRYVP5ZIG4DVgPzjTF5wGdAtwPkyMTfC/U0/jZxNeUD6gP4Hb4HLsc/nisX/7jbih79S4EE4Kfy876N/+EiUBsWdYy11d3tkVhmjLkb6GKtvdjpLCIiItFAPV8iIiIidUjFl4iIiEgd0m1HERERkTqkni8RERGROqTiS0RERKQOOTbJalpamu3QoYNTlxcRByxcuHCHtbZZ9XuGN7VfIrEnmO2XY8VXhw4dyMzMdOryIuIAY8z66vcKf2q/RGJPMNsv3XYUERERqUMqvkRERETqkIovERERkTrk2Jivqng8HrKysiguLnY6Sq0lJSWRnp5OfHy801FEpA5EQ/uldkukboRV8ZWVlUXDhg3p0KEDxhin49SYtZadO3eSlZVFx44dnY4jInUg0tsvtVsidSesbjsWFxfTtGnTiGy4KjPG0LRp04j+BiwiByfS2y+1WyJ1J6yKLyBiG659RcvvISKBi/S/95GeXyRSBFR8GWOGGmNWGmNWG2MmVLH9MmPMdmPM4vKfPwc/at3Jzs7mwgsvpFOnThx55JEMGDCAd999F4Dvv/+eQYMG0bVrV/r27cvpp5/O0qVLHU4sIvuj9svffn355ZekpKTQu3fvvT+fffaZ03FFYlK1Y76MMW7gGWAIkAUsMMbMsNb+tM+uU62114QgY52y1jJixAjGjBnD5MmTAVi/fj0zZswgOzub888/n8mTJ3PMMccAMHfuXNasWUPPnj2djC0iVVD79Vv71aRJEwYOHMgHH3zgcEoRCWTAfX9gtbV2LYAxZgowHNi38QoJj7eMHYV5JLiScJvQPx/w1ewvcMXFM+rS/2N3YSkAKc1accmfruS+e+7mggsvpnvvjL3bevTtD7D3fWWFpV7e/3HzHz63NcxmbU2PrOn1anhcDX/DQK9X81w1PC7AC9b4f50ADgz1f9M/Xi8w8W4X5x6ZXrOL1A1H2y+APQXZYKFhgxYhv9YXX3xBQkIC48aN2/tZ+/btufbaa/nyyy9Dfn0RCUwg1UwbYGOl91nAUVXsd44x5njgF+BGa+3GfXcwxowFxgK0a9cuoIAlZR5ySrbg8zTBeusFdExtfLPwRzp168GGnMI/bFv041LOOnd0lduqklNQyrUzfgh2RJH9sMTjJZFSEvH4f4yHBDwkUOZ/NWUk4iGesr0/CaaMeLzEU0YcZSTgJY4y4k0ZcfiI27vNixsvcfhwGy/xePG6E+HI95z+xQ/E0fYLYHvRTiy2Toqv5cuX07dv3/1unzNnDr179977ftq0aXTu3DnkuUTk9wIpvqoagbnvF+P3gTestSXGmHHAK8DgPxxk7YvAiwAZGRkH/HL99/eX89PmPAAKywqJMzuIc9Vu7pnDWjXi9tMOO+A+zRsmUlg/nkNaNATghuuu5Zt584hPSCA9PZ3WjZP3bht47AD25O3hpCFDeOyfj//xZLsT+eym4/dzpZoNbK3peNiaDqOt6QDcml8v0PPX7X+/oJ7f+jCl+bhK8vyvpXtweQownnxcpfmY0gKMpwCXpxDjKcCUv7rKijCeIkzl17JiKCve+2dT8/63qqO64rGuOHDHQ8WfXXFY4/a/JqUE9Xoh4Hj75SkrxmO9JMftwn8XtGa6t27E3848/KCOGT9+PHPnziUhIYFHHnlEtx1FwkQgxVcW0LbS+3Tgd/fSrLU7K739N/BQ7aP9xoXB4sNVy3853S5DUvyBG7/evXoyY/q7e/d7/rln2bFjBxkZGZw2bCjLlizmvHNGArDg++95++23+eCDD6o8b7zbRZfmDWuVWcKYzwtFu6BgBxTugMKdUJjj/6woBwp3QfFuKNoNxbn+PxfnQkleYOc3bkio7/+JrwcJ9fyvyQ2gUQuIT4K45PLXJIhP9r/GJf726k4sf58I7oTfXit+4hL3Flb+z+IrvY/DUPNCOkw43n7FuePxlHkp83mId9e8+ArE4YcfzrRp0/a+f+aZZ/a2XyISPgIpvhYAXY0xHYFNwCjgwso7GGNaWWu3lL89C1hR22CVv+Fl7cmisKyQQ5ocUtvTVmvw4MHcfvvtPPfcc1x11VUAFBb6bzOOHz+eo446ilNPPXXvgPuKbRJFPMWwZ4v/J2+z/zU/G/K3VXrd5i+wrK/qc8QlQXIqJDeGpMbQuC0k9YCkFEhsBEmN/K+JDctfG0BCg/LXhv5Xd0Lou+qin+PtF8DqHctxY+iY1r22pz6gA7VfIhI+qi2+rLVlxphrgE8AN/CStXa5MeYeINNaOwO4zhhzFlAG5ACXBTNkojuR3JJcvD4vbldovzkaY5g+fTo33ngjDz/8MM2aNaN+/fo89NBDtGzZkqlTp3LbbbexadMmmjdvTlpaGhMnTgxpJgmy0kLYtc7/s3s97N4IuRvKX7P8vVj7ikuCBs2hQQtI7QRt+0P9Zv6fek2hfpr/NTkV6qX6e6HEceHQfgE0ciez3VuEx1NEfAj/v3Gg9gv+OObrzjvv5Nxzzw1ZHhGpmqnrJ+gqZGRk2MzMzN99tmLFCg477I9jsvJK8ti4ZyOdUjqRHEH/qO3v95E64PX4i6sdq2DHL7BzFexcAzm/Qv7W3+8bl+zvmUppCynp/p9GraFhq99ek1LUCxUExpiF1tqIvwd2MO0XQHFxLmvys2gV35DUlMAH6ztB7ZZI1YLZfoXV2o77k+hOBKDEWxJRxZfUAWv9vVXZy2DbT5D9k/91xyrweX7br35zaNoFupwMqR2gScfynw7+nioVVhJCiYmNSMiHPE8BqU6HERHHRUTxleBOwBhDsVdrjsU0a2H3Bti0ELb8+NtPUc5v+6S0gxbdoesp0OxQSOvqL7qSGzuXW2KeMYZG7kR2eksoKyshLi7R6Ugi4qCIKL6MMSS4Eyj1/nEiU4liZSWwaRFkfQ8bv4esBf4B7+B/Gq9FdzjsDGh1BLTs5S+2kho5m1lkPxompbKjYAv5hdtp3CisJ6YVkRCLiOIL/Lcei8qKnI4hoVRW6u/VWjcX1s3xF1wV/5s36QidBkF6P0jPgOaHQ1yCk2lFDkpyUhPiCraQV5qP+mFFYltEFV95JXn4rA+XCWg9cIkEu9bB6s9g1Wfw69fgKQAMtOgBR14GHY6FtkdDg2YOBxWpHWMMjVwJ7PKV4vV6cLtrN2m0iESuiCq+oHzQfZwG3Ucsnw82ZcKK92HlR/6nEAEat4MjRkHnE6H9sf5B8CJRplFSY3IKt5FfuJ2Uhq2djiMiDomY4ivJnQRASVloiy+3203Pnj33vh81ahQTJkxg0KBBbNmyheRk/7W7dOnC22+/HbIcUcXn9d9GXPE+rPjAP9WDKx46DoR+f/I/gdi0i544lKhXL6kp7sJt7CnNI4XgF1+Btl+a30vEWRFTfMW74zEYSrwlIb1OcnIyixcvrnLb66+/rmU6AmUtbF0KS6bC0rf9BVdcMnQ9GQ47y/80op5AlBhjXC4amnjyrAefz4sryJNGq/0SiQwRU3y5jIsEd0LIiy+ppYKdsPh1+PEN/3xbrnh/odXrPOh6qn99QpEY1igxhd3FOygo3EHDBi2cjiMiDoiY4gsgMS6R4rLQzvVVVFT0u+U3/vrXv3LBBRcAcNFFF+3tth8yZAiPPPJISLNEDGthw3zIfAl+mg7eUv9Tiac9Cj3O0fgtkUrq10vDVbyDvJLdQS++Am2/Pv/8c5o2bRrUa4tI4MK3+Ppogv+2VSUtfKV4vB5sfD0MNRgf1LInDHvwgLuo2/4geIr9PVzfv+jv5UpsBEdeDhmXQ3MtTyIxrIr2q4IL6FRWSJm1B9eWqf0SiRrhW3xVoWKKCZ/14TahXWBbDqBoN2T+F+Y/DwXb/BOcnvkU9DwXEuo7nU4k7LldcZR5PXi9HuLcmq9OJNaEb/FVxTe8srJi1u1eQ5uGbWicqMHadW5PNnzzFCz8H5TmQ+eT4LgboMNAPakoUlk1PVQun5eNOT/TwLhJb3poHYUSkXARvsVXFRLKvyGWlJVAiJZG23fMxNChQ3nwQX9DWnnMRFpaGp999lloQoSbwhx/0fXdC/4lf3qcA8de578NIiIHzeVy08jEk2s9+LxluNzBaYoP1H6JSPiIqOKrLp549Hq9VX7+5ZdfhuyaYatkj//W4jdP+f/c8zwYNAGadnY6mUjEa5TUhF1F29hTuC1oE66q/RKJDBFVfIF/stVib2ifeIx5Pp9/uojP/w4F26Hb6TD4DmhxuNPJRKJG/XpNiSvaRl6IJlwVkfAVccVXYlwieaVa4zFkshbCR7f4F7huexSMnuJfyFpEgsoYl9Z6FIlRkVd8la/xWOotJSkuyeE0UaRgB3x2N/zwGjRoASNfgF4XaCC9SAg1SmpCTmE2ewq30bhhG6fjiEgdidjiq8RbouIrGKyFZdNg5i1QkgfHXAvH3wpJjZxOJhL16iU3Ja4wm7ySPBVfIjEk4oqvvU88apmh2svfBh/e5F/wus2RMPxZaK7H3kXqijGGFFciOb4SvGWluOM055dILIi44ktrPAaBtbD8HfjwZigtgJP/DgOugSA97i4igWuUnMrOgi3kFW6jSaN0p+OISB2IyH9tE92J/rm+QsDtdtOz52/zV02fPp1169YxfPhwOnbsCET4HF8l+fDBjbD0TfV2iYSB5KQmxBdsIa90D02CcL7s7GxuvPFG5s+fT5MmTUhISODWW2+lSZMm0dOOiUS4iC2+8kvzQ/LEY1Vro61bt46BAwfywQcfBPVadS77J3hrDOxcDSfeAcfdpN4uEYcZY0hxJ7LTW0JZWTFxtRjLaq1lxIgRjBkzhsmTJwOwfv16ZsyYQZMmTaKjHROJAhE5V0NSXBIWS6m31OkokeOHSfDvwVCcC5e+ByfcqsJLJEykJKdhgdyCbbU6zxdffEFCQgLjxo3b+1n79u259tpra5lQRIIpMosvt/+bYXFZ8CdbrVieo3fv3owcOXLv53PmzNn7+X333Rf064aMpwimXw3vjffP13XlHOh4vNOpRKSSpKTGJAK5nvxanWf58uX07dt3v9sjth0TiTJh2/Xx0PcP8XPOz/vdXugpJM4Vt/fpx0Acmnoot/W/7YD7VHXbEYjM7vr87TBlNGQtgBNu8/+43E6nEol61bVfVfGUFVNqvSS7k3BV8fc0kPZrX+PHj2fu3LkkJCTwyCOPRGY7JhKFwrb4qo7LuPBZn9Mxwte2n2Hyef4C7PzXoPtZTicSkQOIcydQWlaE1+epsvgKxOGHH860adP2vn/mmWfYsWMHGRlapUIknIRt8VXdN7zN+ZvJK8mjW2o3jGZh/701X8CbYyA+GS7/0P9Uo4jUmYPtoaqwbsdPlGLp2rR7jdq1wYMHc/vtt/Pcc89x1VVXAVBYWFijLCISOhE55gv8g+691ovH53E6SnhZ+ApMOhdS2sKfP1fhJRJBGic0wgMUFeXU6HhjDNOnT+err76iY8eO9O/fnzFjxvDQQw8FN6iI1ErY9nxVp/Kg+4MZ91Wd/Pw/DngdNGgQgwYNCto1Quabp+HTO6DLEDjvZUhs6HQiETkIDeu3wFWay+7indSr17RG52jVqhVTpkypcltEtGMiMSBie74q1ngs9gb/iceI9NUj/sKr+wgY/YYKL5EI5HbH09DEkefz4PN5nY4jIiESscWX2+UmwZ0QkukmIoq18Pk9MPte6DUKzvkvuOOdTiUiNZSS1BgvkF9Yuzm/RCR8RWzxBf5xXzHd82UtfHIHzHkMjrwMRjyniVNFIlyD5GbEAbtLcp2OIiIhEnbFl7U24H2T3El4vB68Ydg9fzC/Rw0vAB//FeY/A0eNgzOeAFfY/c8pElOC8ffeuFykuBLIt17KQrSG7f6EvN0SESDMiq+kpCR27twZcAOQVL4GWrj1fllr2blzJ0lJNV+jrVpfPwrfPQdHXw1DHwRNtyHiqINtvw6kYrmhvILs2gcLUJ20WyIChNnTjunp6WRlZbF9+/aA9vf6vGQXZlOUWET9+PohTndwkpKSSE9PD83JF/z3tzFep9ynwkskDBxs+3VA1rKjYCs5ZJPWoKD25wtQSNstEdkrrIqv+Ph4OnbsGPD+1lrGvzmeE9JP4J5j7wlhsjCy/F348C/Q9VQY/rRuNYqEiYNtv6rz7fuP8njOAt474Sk6dTgxaOcVEedF9L/cxhi6Nel20GuoRaw1s2HaFdDuaDjvf3qqUSSKnXX0LbitZXrmv5yOIiJBFtHFF/gXm129e3X0z3S/5UeYchGkHQKjp0BCPacTiUgIpTU7jONdjZiR9wsej5YIEokmEV98dUvthsfnYV3uOqejhE7+dn/hldwYLnnH/yoiUW9k17PZ6TbMWfCM01FEJIgCKr6MMUONMSuNMauNMRMOsN+5xhhrjMkIXsQDOzT1UIDovfXo9cBbY6BgO4x6HRq2dDqRSEQJ5/arOgMzxpPmtby7+l2no4hIEFVbfBlj3MAzwDCgOzDaGNO9iv0aAtcB3wU75IG0b9SeRHciK3NW1uVl687HE2D9PDjraWjdx+k0IhEl3Nuv6sTFJ3NWyiHM8eWxfdtPTscRkSAJpOerP7DaWrvWWlsKTAGGV7HfP4CHgTqddCvOFUfXxl35eVcU9nxlvgwL/gPHXg+9znM6jUgkCuv2KxAj+l6L1xhmfPeI01FEJEgCKb7aABsrvc8q/2wvY0wfoK219oMgZgtYt9RurMxZGV2zM6//FmbeAl1OhpP+5nQakUgV9u1XdTp2PJG+NoHp2xdifT6n44hIEARSfFU1g+feKscY4wIeB/5S7YmMGWuMyTTGZAZlIsJy3VK7sbtkN9mFdTcbdEjlb4M3L4XG7fwLZbvcTicSiVRh334FYmTbk1nntvyw7PU6va6IhEYgxVcW0LbS+3Rgc6X3DYEewJfGmHXA0cCMqgatWmtftNZmWGszmjVrVvPU+6gYdP/Lrl+Cdk7H+Hww/SooyYMLJunJRpHaCfv2KxCnHPUX6vks7yx/pU6vKyKhEUjxtQDoaozpaIxJAEYBMyo2WmtzrbVp1toO1toOwHzgLGttZkgSV+GQJocAUfLE4/cvwOrP4JR7ocUfxgWLyMEJ+/YrEPUaNGdYUis+Ld5KQf5Wp+OISC1VW3xZa8uAa4BPgBXAm9ba5caYe4wxZ4U6YCDqx9enXcN2kV98bV0KsybCIcOg35+dTiMS8SKh/QrUiJ6XU+QyfDz/UaejiEgtBbS2o7V2JjBzn88m7mffQbWPdfAOa3oYy3Ysc+LSwVFaCG//CZJTYfgzWixbJEgiof0KxBHdR9Hl+wd5a+PnnON0GBGplYif4b5Cz7SebMrfRE5xjtNRauaT22HHShj5PNRv6nQaEQkzxuXivNYDWe4qY/mKaU7HEZFaiJriq0daD4DI7P1a8QEsfBmOuQ46n+h0GhEJU2ceczvJPsvUxc87HUVEaiFqiq/DUg/DZVws3bHU6SgHpzAHPrgBWvaCwXc5nUZEwljDRm04Lak1HxVvITd3Y/UHiEhYipriq158Pbo07hJ5xdesu/wF2PCnIS7B6TQiEubO73MVxS7D+9884HQUEamhqCm+wD/ua9mOZZEz0/2vX8MPk2DAeGh1hNNpRCQCdD90JD19cby5Za5mvBeJUFFVfPVI60FuSS5Ze7KcjlI9TxG8fwM06QCD/up0GhGJIOe3O4Vf3ZbMJf9zOoqI1EBUFV8903oCRMatx68fgZw1cMbjkFDP6TQiEkGGDriNRj4fU5f9z+koIlIDUVV8dW7cmSR3UvgXX9nLYd6TcMRo6DzY6TQiEmGS6qUyvEFnPi/LYceOCJ9cWiQGRVXxFeeKo3vT7uE93YTPCzOug6QUOOU+p9OISIQ6r99NlBnDO9/c73QUETlIUVV8gX/c14qcFXh8HqejVG3h/2BTJpz6gCZTFZEa69hhEEeRxNs7FuEtK3U6jogchKgrvnqm9aTEW8LqXaudjvJHxbkw+35ofyz0Ot/pNCIS4c7vNJwtbsNX3z/pdBQROQhRV3xVzHQfluO+5vwTCnfAqfdp7UYRqbXBR/+Fll7L66vedDqKiByEqCu+2jRoQ5PEJuE37mvXOpj/rH+Qfes+TqcRkSgQF5/MqGYZfE8xK1fPrP4AEQkLUVd8GWPokdYj/Hq+Pvs7GLeWEBKRoDp34N0k+SyvL3jc6SgiEqCoK77AP+5rze41FHgKnI7it/F7WP4OHHsdpLRxOo2IRJGUxh04M7kNH5ZsISdnjdNxRCQAUVl89UjrgcXy086fnI4C1sLHf4UGLeGY65xOIyJR6KJ+N1NqDG/NudvpKCISgKgtviBMBt0vm+afWuKkuyCxgdNpRCQKde48hGNIZurOH/CUFjodR0SqEZXFV5OkJqQ3SHd+0L2n2D/Wq2VP/0B7EZEQufjQC9nuNnzy7YNORxGRakRl8QX+cV+O93z98BrkboAh94DL7WwWEYlqx2ZcQwcvTFr7PtZap+OIyAFEbfHVI60HWwu2sr1wuzMBykpg7uPQ9mjodKIzGUQkZrjccVzU+gSWu8r4cfkUp+OIyAFEbfHVs1lPAOduPf7wGuRtgkETNKGqiNSJs46bSEOfZdLi55yOIiIHELXF16GphxJn4vhx+491f/GyEv9s9m2Pgk6D6v76IhKT6jVozjkNuvBZWQ5Zm753Oo6I7EfUFl/Jccl0b9qdRdsW1f3Ff5ikXi8RccRFA/+GAV6dd4/TUURAZjqcAAAgAElEQVRkP6K2+ALo26Ivy3Yso8RbUncXLSv193ql99dYLxGpcy1b9uH0hBa8W7COXbvWOh1HRKoQ3cVX8754fB6Wbq/Dpx4XT4K8LBh0m3q9RMQRlx91G8Uuwxtf3el0FBGpQlQXX32a+xewrrNbj3t7vfpB55Pq5poiIvvo3PkUBpkGvLFrCYWFO52OIyL7iOriq3FSY7o07sKi7Doqvn6cDLkb4QSN9RIRZ11+xFXsdhne/fpvTkcRkX1EdfEF/luPi7cvxuvzhvZCPq9/Xq82R0IX9XqJiLP69rqE3r54Xt38FWWeYqfjiEgl0V98tehLgaeAlbtWhvZCK2fCrnVw7PXq9RIR5xnD5d1GsdkNn3zzgNNpRKSSqC++jmxxJEDobz3Ofw4at4NDzwjtdUREAjToqJvo5DW8vHY61udzOo6IlIv64qtl/Za0adAmtIPuNy+G9fOg/5Vaw1FEwobLHcdl7U5lpcvHN4s0671IuIj64gv8474WZi8M3WKz85+FhAbQ95LQnF9EpIZOHziR5l7Lv5e97HQUESkXG8VXi77kFOewPm998E+etwWWTYM+l0BSSvDPLyJSCwmJDbm85XEsNCUs+PF/TscREWKo+IIQzfe14N/+Jx2PGhv8c4uIBME5g+6nqdfywuJnnY4iIsRI8dWxUUdSk1JZmL0wuCcuLYTMl+HQ0yG1U3DPLSISJMn1Urm8+dF8RxE/LH3d6TgiMS8mii9jDH2a9wl+8bVkKhTlwNFXB/e8IiJBdt6JD9DEZ3lh0VNORxGJeTFRfIF/0P2m/E1kF2QH54TW+qeXaHUEtD8mOOcUEQmRevWbMSYtg3kUsuSnN52OIxLTYqb42jvfV7DGfa3+HHas9Pd6aVJVEYkAo058kBSf5YXMx52OIhLTYqb46pbajXpx9YJ36zHzJajfHA4/OzjnExEJsfoNWjImtS9f23yWr5zudByRmBUzxVecK44jmh0RnJ6vPVvhl4+h92iIS6j9+URE6sjowQ/S0Gd58ftHnI4iErMCKr6MMUONMSuNMauNMROq2D7OGLPUGLPYGDPXGNM9+FFr78gWR7J612pyS3Jrd6LFk8F6oc+lwQkmIiETLe1XsDRo2JpLmvTiC18eP//ygdNxRGJStcWXMcYNPAMMA7oDo6tonCZba3taa3sDDwP/DHrSIMhomYHFkrk1s+YnsRZ+eA3aHwtpXYIXTkSCLprar2C66ER/79cz8+93OopITAqk56s/sNpau9ZaWwpMAYZX3sFam1fpbX0gROv41E6vtF7Ui6vHt1u+rflJ1s+DnLX+Ge1FJNxFTfsVTI1S2vF/Tfvypd3D4uVTnI4jEnMCKb7aABsrvc8q/+x3jDHjjTFr8H9zvC448YIr3h1Pv5b9+HZzLYqvRa9CYiPoPrz6fUXEaVHTfgXbhSc9RqrX8lTmY6Fb91ZEqhRI8VXVPAp/+JtqrX3GWtsZuA24s8oTGTPWGJNpjMncvn37wSUNkgGtB7BhzwY25W86+IOLdsNP70HP8yChXvDDiUiwRVX7FUz16jdjbKuBLKCY+YtecDqOSEwJpPjKAtpWep8ObD7A/lOAEVVtsNa+aK3NsNZmNGvWLPCUQTSg1QCAmvV+LX0Lyoqhrwbai0SIqGq/gu28Ex+ildfy1JIXsD6f03FEYkYgxdcCoKsxpqMxJgEYBcyovIMxpmult6cDq4IXMbg6pnSkeb3mNSu+Fr0KLXtC697BDyYioRBV7VewJSQ14qp2w1jmKuOL+Y85HUckZlRbfFlry4BrgE+AFcCb1trlxph7jDFnle92jTFmuTFmMXATMCZkiWvJGMOAVgP4but3eH3ewA/cvBi2LoG+Yfuricg+oq39CoUzT/gHHbzw9MpJeMs8TscRiQlxgexkrZ0JzNzns4mV/nx9kHOF1IDWA3hvzXv8nPMzh6cdHthBP7wGcUnQ89zQhhORoIq29ivY4uKTuKbredy89i1mzr2HMwf9w+lIIlEvZma4r+zoVkcDBD7lhKcIlrwFh50FyU1CmExEpO4NOeZ2DvO5eGbtdDylhU7HEYl6MVl8NU1uSrcm3QIf97XyIyjJhT4XhzaYiIgDXO44ruv+Jza5YernNzsdRyTqxWTxBf5bjz9s+4GisqLqd142DRq2gg7HhT6YiIgDju13DQNI4vmtX5Obu7H6A0SkxmK3+Go1AI/Pw8LshQfesTgXVn0Kh48El7tuwomI1DHjcvGXARPJM/DvWTExz6yIY2K2+Orboi8JroTqbz3+/CF4S6HHOXUTTETEId0OOZPhCS2YnL+KjVnfOR1HJGrFbPGVFJdEnxZ9qh90v2waNG4HbY6sm2AiIg66ZvA/ibPw5Je3Oh1FJGrFbPEF/luPq3atYnvhfpYKKdgJa2b7e71MVauUiIhElxYtj2BM4x584s1h8fKpTscRiUqxXXy19i81NH/L/Kp3+Gk6WC/00NxeIhI7Lh/yJGley2MLHtayQyIhENPF16Gph9Ikscn+x30tewfSukGLACdiFRGJAvUatOCatqew2JQy69sHnY4jEnViuvhyGRdHtzqabzZ/g8/u8+0ubzOsn6dbjiISk0YMeoAuXsM/V75BcdFup+OIRJWYLr4Ajm97PDuLd7J8x/Lfb1g+HbDQ42xHcomIOMkdn8iEI8azyQ3/+/Rap+OIRJWYL74GthmI27iZvXH27zcsmwYte0FaV2eCiYg47Kgjr+QUVwr/zfmBzZsznY4jEjVivvhKSUyhb4u+fJn15W8f5vwKmzK1iLaIxLybT3oCgEc/v9HhJCLRI+aLL4BB6YNYtWsVWXuy/B8sf8f/evhI50KJiISBVq0zuCK1N7N8u/l24QtOxxGJCiq+gBPbngjAV1lf+T9Y9g60Pco/uaqISIwbc+rTpHvhwSXP4PEUOh1HJOKp+ALaNmpL55TO/nFfOb9C9jLoPtzpWCIiYSExqTG3HXYpa12WybN0+1GktlR8lRvUdhALty4kb8V7/g+6DXM2kIhIGDnhqL8wkHo8lz2P7dtXOB1HJKKp+Co3qO0gymwZ81a/B80Og9ROTkcSEQkbxuXitkEPU2rg0U+vdjqOSERT8VWuZ1pPUhObMDt/g3q9RESq0L79Cfw5pQczy3YwL/MZp+OIRCwVX+XcLjcnNGjP3OQkPF1PcTqOiEhY+vOwF+nghX8seZ6iol1OxxGJSCq+KhmUn88et4uFcU4nEREJTwlJjZjY+3o2ueG5mX92Oo5IRFLxVaGslAHrFpKIiy83feV0GhGRsNWv758ZEd+cV/esZOWqmU7HEYk4Kr4qrJ9HcnEuAxp348uNX2KtdTqRiEjY+svQf9PIwt/n3om3zON0HJGIouKrwsqPIC6JQYeMZFP+JlbtXuV0IhGRsNU4tRO3dDiLpS4PUz+/yek4IhFFxReAtf7iq9OJnNBhCACzN8yu5iARkdh2xgn3MoAknto8m63ZPzodRyRiqPgCyF4Ouf4pJtKS0+jdrDefrv/U6VQiImHNuFzcdeLj+IC7P7kS6/M5HUkkIqj4An+vFwYOGQrA0I5D+WXXL6zZvcbZXCIiYa5tu+O4vsVA5tkCps/+q9NxRCKCii+AlTMhPQMatgDg1A6n4jIuPl73scPBRETC3+hT/0WGTeDhDR+ydetip+OIhD0VX3lbYPOi381qn5acRr+W/fj414/11KOISDVc7jjuGfwUXuDuT8fp9qNINVR8/VLeu9XttN99PKzDMNblrWNFjhaQFRGpTtt2x3JDy+OZZwt4d/YEp+OIhDUVX6s+hcbtoNmhv/v45PYnE2fi+PhX3XoUEQnEqFOeop9N5JENM9m6RbcfRfYntosvrwd+nQOdTwJjfrcpJTGFY9ocw8frPsZn1YUuIlIdlzuOe07+F15g4qdX4vN5nY4kEpZiu/jatBBK90DnE6vcPLTDULYUbGHJ9iV1HExEJDKlpw/g5taD+ZZC3vj0OqfjiISl2C6+1swG44KOx1e5eXC7wSS6E/no14/qOJiISOQ6b8gTnGDq88+tX/HLmk+cjiMSdmK7+Fo7G1r3geQmVW6uH1+f49OP55N1n1DmK6vjcCIikcm4XPz9tFdo6IPbvr6V4uJcpyOJhJXYLb6KcyErEzpVfcuxwrCOw9hZvJPM7Mw6CiYiEvmapnXj3sOvYLXLxxPvX+J0HJGwErvF169zwHr3O96rwsA2A6kXV09PPYqIHKTjjrqei5La8Xrhr8xZ8LTTcUTCRuwWX2tnQ3x9SO9/wN2S4pIY3G4ws9bPwuP11FE4EZHocOPw1+niNdy17Hl27lztdByRsBC7xdea2dDhWIhLqHbXYR2HkVeax9xNc+sgmIhI9EhMasxDx93PHgN3fngJPq/Gz4rEZvG1ewPkrKl2vFeFAa0HkJacxjur3wlxMBGR6HPIIWdwS8tBzLX5vPTRWKfjiDguNouvNbP9r9WM96oQ74pneOfhzMmaw/bC7SEMJiISnS445SmGuhvzrx3fk7nkVafjiDgqoOLLGDPUGLPSGLPaGPOHRbuMMTcZY34yxiwxxnxujGkf/KhBtHY2NGz1hyWFDmRk15F4rZf31rwXwmAiEmxR135FKONy8bfhU2nrM9ya+TA7d65yOpKIY6otvowxbuAZYBjQHRhtjOm+z24/ABnW2l7A28DDwQ4aND4frP0KOg36w5JCB9K+UXsyWmTwzqp3sNaGLJ6IBE/UtV8RrkHD1jx23H3kGZjwwUV4y/QQk8SmQHq++gOrrbVrrbWlwBRgeOUdrLWzrbWF5W/nA+nBjRlEW3+EopyAx3tVdnbXs9m4Z6Pm/BKJHNHVfkWBboecxV/ThzCfIl788HKn44g4IpDiqw2wsdL7rPLP9udPQPiux1Mx3qvToIM+dEj7ITSMb8i0VdOCGklEQia62q8ocfZJj3GGuynP7VrMNwufdzqOSJ0LpPiq6t5clffdjDEXAxnAI/vZPtYYk2mMydy+3aGB62tnQ4se0LDFQR+aFJfEaZ1O47P1n5FbouUyRCJAdLVfUcK4XNw1YiqdfS5uXfI0G7O+czqSSJ0KpPjKAtpWep8ObN53J2PMycAdwFnW2pKqTmStfdFam2GtzWjWrFlN8tZOaSFsmF+jXq8K53Q9hxJvCTN/nRm0WCISMtHTfkWZeg1a8NTJz+IDrp81lsLCnU5HEqkzgRRfC4CuxpiOxpgEYBQwo/IOxpg+wAv4G65twY8ZJBu+BW9pwFNMVOWwpodxWOphvLNKc36JRIDoab+iUNt2x/HI4VeyxniZ+M5IrM/ndCSROlFt8WWtLQOuAT4BVgBvWmuXG2PuMcacVb7bI0AD4C1jzGJjzIz9nM5Z6+eBcUPbo2t1mrO7ns3POT/z086fghRMREIhqtqvKHVs/2u5vmkGn3h38dLMPzsdR6ROxAWyk7V2JjBzn88mVvrzyUHOFRrrv4XWvSGxQa1Oc1qn03g081HeWfUO3Zvu+9S6iISTqGm/otjlp7/Eiskn8OSO7+m24GmO63eN05FEQip2ZrgvK4FNC6HdgFqfqlFCI4a0H8KHaz+k0FNY/QEiIrJfxuXi7yPf4RCfm1uXPc/adV86HUkkpGKn+Nq0CLwl0P6YoJxu1KGjyPfkM3319KCcT0QkltWr34wnT3mBeAvjv7iWnJw1TkcSCZnYKb42fON/DULPF8ARzY6gV1ovXl/xOj6rQaIiIrXVJv1o/tX/TrYby/UzzqekJM/pSCIhETvF1/pv/Ws51ksN2ikv6X4JG/Zs4Ousr4N2ThGRWNarxyju73I+i00pd719lp6AlKgUG8WXzwsbvwtar1eFk9ufTMv6LXntp9eCel4RkVh2ysCJXN+4Dx+V7eTZGRc7HUck6GKj+MpeBiV5QRvvVSHOFceFh17I91u/5+ecn4N6bhGRWPanM//HyPgWPJ+7lBmz73A6jkhQxUbxtf5b/2uQe77AP+dXclwyk36aFPRzi4jEKuNycdc573GUTeJv699jbuYzTkcSCZrYKL42fAMpbaFx2+r3PUgpiSkM7zycmb/OZEfRjqCfX0QkVsUn1ufxke/Sxbq5aelzLPnpLacjiQRF9Bdf1vp7vkLQ61Xh4u4XU+YrY+rKqSG7hohILGqYks5zp0+mqTWM/+7vrF3/ldORRGot+ouvnLVQsA3ah674at+oPSekn8CbK9+kxFvlmrwiIlJDac0P54WTnsVlYdzn17A1e6nTkURqJfqLr/UV83sFd7D9vi7ufjE5xTnMXDuz+p1FROSgtGs/kOcG3EOusVw182Jyd29wOpJIjUV/8bXhW0hOhWbdQnqZ/i37061JN15e/jJenzek1xIRiUXdDzubJ3uNZ73xctW7I8jP3+p0JJEaif7ia/03/vFexoT0MsYYxvYay6+5v/Lp+k9Dei0RkVh19JFX8WjXi1lhShn/9hkUFu50OpLIQYvu4itvC+z6Nejze+3Pye1PpkvjLjz/4/Pq/RIRCZHBx07ggY7nsphirntrGMVFu52OJHJQorv4qljPMYSD7StzGRdXHnEla3PXMmv9rDq5pohILBp6wt38o90ZfG8LuemtYXhKCpyOJBKw6C6+1n8L8fWh5RF1dslT2p9C55TOPP/j81pwW0QkhM4a/CB3tTqJOTafW948FY+nyOlIIgGJ7uJrw7fQth+44+rski7jYtwR41iTu0Zjv0REQuy8U59kQrNj+dyXy81ThlBaku90JJFqRW/xVZwL2ctDOrnq/gxpP4ROKZ144ccX1PslIhJiF532PH9tdhxf+HK5YeoQSkrynI4kckDRW3xtXgxYSM+o80u7XW7GHTGO1btXa+yXiEgduPC055jYcjBzbD7XThlCUdEupyOJ7Ff0Fl+bFvpfW/d15PKntD+FjikdNfZLRKSOnHfqk9zT5lTm2wLGv3kKhQVab1fCU3QXX6mdoV6qI5d3u9xc2etKVu9ezafrNPZLRKQujDz5Ue5vP5yFtogr3zqV3LyNTkcS+YMoLr4WQZsjHY0wtMNQujbpyhOLnqDUW+poFhGRWHHGiffxSOfzWUYJl087g23bljsdSeR3orP4ytsMezY7Xny5XW5uzriZTfmbmLxisqNZRERiySkDJ/Jsj6vJwsulH4xi/YZ5TkcS2Ss6i6+K8V4OF18Ax7Q+hoFtBvLCkhfIKc5xOo6ISMwYkHE1Lx01kULj49LPr2TFL+87HUkEiObiyxUHLXs6nQSAv2T8haKyIp5b/JzTUUREYkqP7ufzyqB/kWjh8nl/5bsf/uN0JJEoLr5a9ID4JKeTANC5cWfOPeRc3vrlLdbuXut0HBGRmNKx42BeHfYarayLcT8+wXuz73A6ksS46Cu+fD7Y9ENY3HKs7OreV5Mcl8xjCx9zOoqISMxp2aoPr5zzARkkc+eGGTw9/UKsT9MAiTOir/jauQpK94Rd8ZWalMrYXmP5Outrvtn8jdNxRERiTqOUdjx74ZeMjG/OC7lLmfDGyVqOSBwRfcVXGA2239eFh11ImwZteDTzUcp8ZU7HERGJOfEJ9fn7qFlc37g3M8u2c8Ubg9i1S8NBpG5FZ/GV0BDSujqd5A8S3YncdORNrNq1SlNPiIg4xLhc/Hn4azzcfiTLKGb0u8NZufojp2NJDInO4qt1b3C5nU5SpSHth3B8+vE8vfhpNuVvcjqOiEjMGjboHl7pNxEPlkvm3MIn8x5wOpLEiOgqvjzFsHVZWN5yrGCM4c6j7gTgH/P/gbXW4UQiIrGrx+HnM/WMKRxCPDevnsy/po/G59WwEAmt6Cq+speBzxPWxRdAqwatuK7PdczbNI+PflVXt4iIk9Ka9+Cl0V9xdnwLXsxdxnWTj9eakBJS0VV8VQy2T89wNkcARh86mh5Ne/DQgofYXbzb6TgiIjEtIakRd4/6lNubD2SeN48L3j6N5SvfczqWRKnoK74atoJGrZ1OUi23y83dx9xNbkkuj2Y+6nQcEZGYZ1wuRg97lleOvB0vlku+vYOpn96g+cAk6KKv+ArzW46VdUvtxmWHX8Z7a97juy3fOR1HRESAXj0v5K0R79Hf1OPeLZ8z4Y2TKCzY7nQsiSLRU3wV7YKdq6FNX6eTHJRxR4yjbcO2/O2bv5Ffqsn+RETCQePUzjx78TyubXwEH3u2c8GbJ+k2pARN9BRfm3/wv0ZQzxdAUlwS9x13H1sKtnD/d/c7HUdERMq53PGMHT6J//S8jkLr4+Jv7+DlD6/Q05BSa9FTfFUMtm/dx9kcNdCneR/G9RrH+2vf54O1HzgdR0REKul35FimjZjBIFcK/9wxn7GTjmHbtuVOx5IIFkXF1yJo2hWSUpxOUiNX9LqCPs37cO/8e9m4R484i4iEk8apnfjnxXP4W6uTWOIr5JwPL2DWNw86HUsiVEDFlzFmqDFmpTFmtTFmQhXbjzfGLDLGlBljzg1+zABsWRKRvV4V4lxxPDDwAQyGCXMmaO1HkSCJiPZLIoJxuTj3lCeYcsKTtLZx3LTqdW59/UR2717ndDSJMNUWX8YYN/AMMAzoDow2xnTfZ7cNwGWAMwsWFuZAXha07OHI5YOlTYM2TBwwkSXbl/D8j887HUck4kVE+yURp1Onk5h08TzGp/Rklmc7I989gy+/e8LpWBJBAun56g+sttautdaWAlOA4ZV3sNaus9YuAZyZDCV7mf+1ZU9HLh9MwzoO46zOZ/Hvpf9mwdYFTscRiXTh335JRIpPqM+4EZN5Y8C9pFo31/78X+6YfDK5uzc4HU0iQCDFVxug8iCkrPLPwsfWpf7XFpFffAHcftTttG3Yllu+uoWtBVudjiMSycK//ZKIdmi3EUy5aB5XNDyUD0u3cta7p/HRnHs0MascUCDFl6nisxqtBm2MGWuMyTTGZG7fHsQJ67Yu9c9s36BZ8M7poPrx9XnyxCcp9hZz/ezrKS4rdjqSSKQK//ZLIl58YgOuO/stpg64j9bEcevat7h60rFs2rzQ6WgSpgIpvrKAtpXepwOba3Ixa+2L1toMa21Gs2ZBLJS2Lo2KW46VdW7cmQcHPsiKnSu4+9u7sbZG/16IxLrwb78kanTrNpxJF3/HbWkDWOjdw8hPxvC/mWPxeIqcjiZhJpDiawHQ1RjT0RiTAIwCZoQ21kEoK4HtP0OLyB5sX5VBbQdxTZ9r+HDth7yy/BWn44hEovBuvyTquOMTufj0F3nvlJfo56rPY9u/5bxJR/P9Dy85HU3CSLXFl7W2DLgG+ARYAbxprV1ujLnHGHMWgDGmnzEmCzgPeMEYU3ezz23/GXxlUdfzVeGKnlcwpP0QHl/0OPM2zXM6jkhECfv2S6JWqzb9efqSb3mqy8UUWy9/WvI4t7w+iOxty5yOJmHAOHU7KyMjw2ZmZtb+RD9MgvfGwzULIa1L7c8Xhgo9hVzy0SVsKdjCpGGT6NS4k9ORRGrEGLPQWpvhdI7aClr7JTGhuDCHlz65mv/mLsNtYWxaPy4++XGSkhs7HU0OQjDbr8if4X7rMoivD6kdnU4SMvXi6/HkiU+S4Epg7KyxbMnf4nQkEREJUFK9VK4eOYXpg19ggLsRT+ZkctYbA/nw67vx+bxOxxMHREHxtRRaHA4ut9NJQiq9YTrPD3meAk8BY2eNZVfxLqcjiYjIQWjb7lievPQbXupxLY2Niwm/TuPiV/rxw1LN7xtrIrv4sjYqn3Tcn0NTD+Vfg//FloItXP3Z1RR4CpyOJCIiB6nfkWOZcmkm96afRrYt5dJFD3Dta8eyau0sp6NJHYns4mv3BijJjfhlhQ5GRssMHjn+EVbkrOCG2TdQ6i11OpKIiBwklzue4Sc9xPsXzObaxkeQWZbLOV/fyO2TTyJrk1Y3iXaRXXxVzGzfspezOerYie1O5O5j7mb+lvlMmDMBj8/jdCQREamBevWbMXb4JD4e8T6XNejKpyXZnDnrcu5983S2Zi91Op6ESOQXX8YFzfddJzf6jegygpszbmbW+lnc8tUt6gETEYlgKU06ctO57/LhsEmMTGzNtML1nPbRaO6dejpbty52Op4EWWQXX9nLoGkXSKjndBJHjDl8DBP6T+DzDZ9zw+wbKPGWOB1JRERqoUXL3kwc/SkfnPxvhie2YlrReoZ9fDH/mHoam7douaJoEdnF19YlMTPYfn8uOuwiJg6YyNxNcxn/+XgKPYVORxIRkVpqkz6Av42excwhL3F2YhveKdrAaZ+M4fbJJ2lgfhSI3OKraLd/wH2MF18A5x1yHvcddx8Lti7gqs+uIr803+lIIiISBK3a9Oeu0Z/w0amvclG9TnxWks3Zc27imlePYdHS152OJzUUucVXdvkSDS1UfAGc2flMHj7+YZZsX8L/ffJ/ZBdkOx1JRESCpGWrvtxy/gxmjXif8Sk9WVKWy5hFD3LR//ry0Zx7tHh3hInc4mvvk44qviqc2uFUnhr8FOvz1nPhhxeyYucKpyOJiEgQpTTpyLgRk/nkgq+4o/nx5Po83Lr2LYa91o//fvB/5O7e4HRECUBkF1/1m0PDFk4nCSsD0wfy6rBXcblcjPl4DF9u/NLpSCIiEmTJ9dMYNewZZoxZxNOHjKGDK4kndi7g5HdPY+Ibp7B85XtOR5QDiODiS4Pt96dbajcmnzaZzimdue6L63jtp9dwagF1EREJHZc7nhMG3Mx/Lsvk7QH3c3pSKz4u3syo+Xdy4ct9eW/27RQX7XY6puwjMouvslLYvlLF1wE0q9eMl4a+xMntT+bhBQ9z57w79SSkiEgU63bImdw9ehafjfyQCc2OJd96uHPD+wyechz3v3kGK1fNdDqilIvM4mvHL+AtVfFVjeS4ZB494VGuPuJq3l/zPhd+eCFrdq9xOpaIiIRQo8btuei053lvzGL+2+MaBsY3ZVrhOs795jYueLkPUz+9gby8LKdjxrTILL402D5gLuPiqt5X8cKQF1uQJkIAABNLSURBVNhVsovRH47m/TXvOx1LRERCzLjd9D/ySh66+Cu+GP4+E5odR5n1ce+WzznxnaHcNOl4vpz/uJ6UdEBkFl/bfgJ3AqR2djpJxBjQegBvnfkWhzc9nNvn3s7EeRMp8BQ4HUtEROpASpOOXHTac7x92Q9M6Xc35yW3Z2HpTq5d+RInTerH/W+eweLlU7E+n9NRY4JxaiB2RkaGzczMrNnBr58HeZvhqnnBDRUDynxlPLv4Wf6z9D+0qt+Ke469h6NaHeV0LIkRxpiF1toMp3PUVq3aL5Ew4Skt4JuFzzFj9XS+9O6m1Bhae+HUlEM4rdf/0e3/27v76KjqO4/j719mJs8kIYEEQgiQEGICCAjSUpTqUWnVKm6rxWrPth7Odrd27fZ0j6fr2q7bnrN7tvWf0rU91bXr0y5btban2OragqVSUQuLgCQE88BTCBgIhJDnebj7x71JhzEDoUzm3mE+r3PumTt3fvO73/lN5pvv3KepvhmTkZrbaCZCIvNXahZf318IFcvhzp8kNqg08m7nu3zrzW9xqOcQa2vX8vWlXyc3kJ6/kSnJo+JLxJt6e47yux2P8cqR13nb6iNkDLMihhsLa7mp/h7q561J+0IsvYuv4T7413K4/pvw8QcTH1gaGQgN8Ni7j/Fc43OU55fzyIpHWFG+wu2w5DKm4kvE+053tbDp//6d33RsYzsDhI1hehhumFTFDbWfYXH9Wvz+LLfDTLr0Lr6O7oT/uB4++xzU3574wNJQ9Fawm2bdxIPLHmR6/nS3w5LLkIovkdTSfaqVLTt/zKajW9lm9RI0hsKIxaqsMq6bdQMrF60jLz89LnaeyPzlT0QnSXViv3079Qp347iMLCldwku3v8TTe5/myfeeZGv7VtYtXMd9C+4jy5d+325ERMRWVFzNHTc+yh1A79kOtu16ki2Hf8cbQ8d5ueV/8DdvYKnJ5dopi7mm/m6qZl2X9rsnxyP1tnz99hF464fw8HHwpV7t6HXHeo/x6I5H+e2h3zIjfwZfW/o1Vs9aTYbRh0kunbZ8iVweQsEBdjf8lC3Nv+QPZ9to8dm1RHkEVuZWsKLi4yxfeC+FBTNdjjRx0nu344a10H0Y7n8r8UHJqLePvc13//hdWrpbqCuu44ElD3DNjGswxrgdmqQwFV8il6djR7fzh4b/Yuvx7bwT7qE/w5BhWcy3Mvlo0TxWzPkEV9Z9hqysArdD/bOld/G1fhGUL4G7nk54THKucCTMKwde4Ue7fkR7bztXlV7FA0seYNm0lP/fKS5R8SVy+QsO9/Hevhd5q/VV3jrzPntNkLAxZFoWi0wOVxfVsnzOJ1hY+xdkZuW7He64pW/xNdxvn+l43UNw3TcmJjD5kGA4yC9afsHjux+nc6CTq0qv4r4F97GqYpV2R8pFUfElkn7OnjnCzsbn+eORN9jed4gmE8ZyirEFZLO0oJqrZl7L4ro7yc+f5na4caVv8XVsNzy+Cu56BubfMTGBSVyDoUFean6JZxuepaOvg+rCar644IvcOudWAr6A2+FJClDxJSJnTh9gR+Pz7Ox4i529h9hnQoSNwVgWNZaPRTnTWVR2FYuqb2FWxcc8cwB/+hZfu5+HX3wJ7n8HSnW2o1uCkSC/Ofgbntr7FPtP72dKzhQ+XfNp7pp3F9PyvPutRdyn4ktEYvX3fsDuppd498gb7O5pY0+kn94M+/jiwojFAl8+CwuqWDj9Iyyo+RTFLv20YPoWX5u+Ddt+4JzpqC0tbrMsi20d29jQtIGt7VsxxrCqYhVra9fysfKPaZekfIiKLxG5kEg4SGvbJnYdeI29Jxt4b/ADWjMiRJwTvqaHoT5QSF3BHOrLl1Nf9UlKSmomPK70vc7Xif1QMleFl0cYY1g5YyUrZ6zkaO9Rfvb+z/h588/ZcmQLZbll3Fp1K7dV3cbcyXPdDlVERFJEhi9ATc3N1NTczF3Osv7eD2hs/jV7j75JY3czjcOn2dy9C7p3QeMTTA1b1PrzuSJ/JrVTr2TezGuprFjh2Svxp9aWrx8sgWkL4bPPTkxQcsmC4SCbj2xmY8tGtnVsI2yFqSuu49aqW1k9a7WunJ/mtOVLRBKlt+co+1peZd+xd9jf3cr+oS5aM8KEnC1kmZZFteWnJquEmsIq5pYuorpiJdPKFv1Zx5Gl527H4IB9puOqB+H6f5y4wCRhTg6c5LWDr/Fy68s0dDUAsKBkATfMuoEbK29kduFsdwOUpFPxJSITaXiwh7ZDr/P+0bdp7mqiub+D5nA/nb4/XaMyN2JRTYCqrBI+v/SrXFE7vp8qTM/djl0tYEX0s0IpZErOFO6tu5d76+7l4JmDbD68mc2HN7N+53rW71xPVWEV1864lpUzVrK0bCmZvky3QxYRkRSWmV3AFbV3cEXtuVdE6D7VSuvhN2jt3EVbdyutAx+wbfA4awZOuhJn6hRfnU32rYqvlDS7cDbrFq5j3cJ1HO87zuuHX2fLkS1saNrAM43PkOPPYfm05awoX8HV065mbtFcHbAvIiIJUVRczdLiapbGLLciEVfiSZ3i60QTGB+UuHOKqSTOtLxp3FN3D/fU3UN/sJ8dH+xga/tW3ux4k9+3/x6AoqwilpUtY9m0ZSwpXcK8yfPwZ6TOn6uIiHifW9cQS53/Ziea7MLLo2cuyJ8nN5DLqopVrKpYBUBHbwfbj28fnTYd3gRAjj+H+SXzWTR1EVdOvZL6knrKcsv0W5MiIpJyUqj42q8Lq6aB8vxy1sxdw5q5awC7GNvVuYs9J/ewu3M3zzQ8Q8gKAVCcXUxdSR31xfXUFtdSM7mGykmV2kImIiKelhr/pUJDcKpNPymUhsrzyynPL+eWqlsA+yeOmk41se/UPhq7GtnXtY+nOp4aLcgyMzKpKqpibtFcqgqrmFM4h9kFs6ksqNQB/SIi4gmpUXx1tYAV1sH2QrY/m8Wli1lcunh02VB4iLbuNlq6W2g+3UxzdzPbj2/nV22/Gm2TYTIozyunsqCSmZNmUpFfYd9OqmB6/nQmBSZpF6aIiCRFahRfJ3Smo8SX5cuirqSOupK6c5b3Bfs42HOQg2cOcuDMAQ71HOLI2SPsPbmXnuGec9rmBfKYnjed8vxyynLLKMstozS3dPR2au5UCjILVKCJiMglS43iq7MJTIb900Ii45QXyGN+yXzml8z/0GNnhs7Q3ttO+9l2jvcdp6O3g46+Do71HmPPiT10D3V/6DmBjABTcqYwJWcKJdklFOcUMzlrMsXZxUzOnszk7MkUZRVRmFVIYVahtqaJiMiYxlV8GWM+CawHfMCTlmX9W8zjWcCzwFKgC1hrWdbBhEV5ogkmz4FAdsK6lPQ2UiCNVZiBvSuzs79zdDo5cPKcqaOvg4auBk4Pnh493iyWz/iYlDmJgswCe8qyb/MCeUzKnER+IJ/8zHzyAnn25M8jN5A7ej/Hn0NuIJfMjEwVcZfA9fwlIhLjgsWXMcYH/BC4CWgHthtjNlqW1RjVbB1w2rKsucaYu4HvAmsTFuWJ/VBad+F2IgmS5cti5qSZzJw087ztLMvibPAspwZO0T3UzZmhM3QPdY/O9wz30DPUQ0+wh7NDZ+no7aA32EvvcC+D4cFxxeIzPnL8OeT4c8j2Z5Ptz7bnfdlk+bLI9tu3I/OZvszR+5kZmWT67CmQEbDnMzIJ+ALn3Poz/AQyAgQyAva8L4Df+PFl+PBn+O3J+FOuCPRE/hIRiTGeLV/LgRbLstoAjDE/BdYA0clrDfDPzvzPgMeMMcZKxA9HhobhVCvUfeqSuxJJNGPM6JatixWMBOkb7uNs8Cz9wX76gn32FOpjIDhAf6ifgdAA/UH7diA0wGB4kIGgcxsaoDfYy1BoiMHwIEPhIYbCQwyHhxkKD03Aq+WcgsxnfBRlFfHrT/96QtaVIO7mL+DbLzfQ2NFz4YYiknT15QU8ctvYe0Am0niKrxnAkaj77cBH4rWxLCtkjDkDlADn/GiSMeZLwJcAKisrxxfhqVaIhHSwvVx2AhkBirKLKMouSnjflmURjAQZDA8yHB4mGA4yHBlmODzMcMS+H4wER2+HI8OEIiFCkdDo8pAVGl0WioQIWSHCkfA58ylw+Q5385eIyBjGU3yNtZ8h9hvheNpgWdYTwBMAy5YtG9+3yuJq+OutUFgxruYiYm+RG9ndmObczV/gyrdqEfG28fyoUTsQfeBLBdARr40xxg8UAqcSESD+TJh+JeQWJ6Q7EUkr7uYvEZExjKf42g7UGGPmGGMygbuBjTFtNgJfcObvBF5P1PESIiKXQPlLRDzngrsdnWMg/hZ4DftU7f+0LKvBGPMdYIdlWRuBnwDPGWNasL8x3j2RQYuIjIfyl4h40biu82VZ1ivAKzHL/ilqfhC4K7GhiYhcOuUvEfGa8ex2FBEREZEEUfElIiIikkQqvkRERESSSMWXiIiISBKp+BIRERFJIhVfIiIiIkmk4ktEREQkiYxbF3I2xpwADl3EU6YQ80O3HqG4Lp5XY/NqXODd2C42rlmWZU2dqGCSRfkrKbwam1fjAu/G5tW44OJiS1j+cq34uljGmB2WZS1zO45YiuvieTU2r8YF3o3Nq3F5jVfHyatxgXdj82pc4N3YvBoXuBebdjuKiIiIJJGKLxEREZEkSqXi6wm3A4hDcV08r8bm1bjAu7F5NS6v8eo4eTUu8G5sXo0LvBubV+MCl2JLmWO+RERERC4HqbTlS0RERCT1WZaV0An4JLAfaAH+YYzHs4DnncffAWZHPfaQs3w/8IkL9QnMcfpodvrMjLeOqD66sE8rjV2HK3E5y78ARIAhoAP4sQtjtgrYCYSAO2P6+MAZs2bgCx6Ky3LGaxewMclxfR1oBPYAm7FPQR7p5xgw7IzZ6Hi5GZtHxuxvgPecdf8BqL/QOpI9xXtNbucJlL8um/wV1Y8VNWYbXRiz8+UJN8fsfLnV7TFLWA5LdOLyAa1AFZAJ7I4OzmlzP86HE7gbeN6Zr3faZzkvvNXpL26fwAvA3c78j4Evx1nHC04fq503dA9w08g6XIzreaefQ84b5uaYzQauBJ4F7ozqYwnQBuwFPurMT/ZAXFVAr4vjdT2Q68x/Oeq9PAAcBsqcMTsyMl4uxvYC3hizgqj13Q787/nWkcjcpPyl/IUH8lfM31mfy2MWL0+4PWbxcqsXxixhOSzRyWsF8FrU/YeAh2LavAascOb92NW1iW070i5en85zTgL+2HWPsY5uZ9lDMVN0OzfiOuk8/ntgr5tjFtX2aezktcLp73PA41F9PA58zu24nGW9bo+Xs3wJ8Kbz+G7g8ah+3hkZLxdj2+PBMfsc8OpYf/PRsSRzivealL+UvxIVV/TjQK8Xxiw2T3hlzKLiGsmtXhuzS8phiT7mawb2N/0R7c6yMdtYlhUCzgAl53luvOUlQLfTR+y6YtcxCHRGLR9pe77nJCOuM8AV2Lup5hhj3gX+Eoi94FsyYos10ofbYxYvLoBs4K+AvzPG3DFWmyTFtQ541XlsMKqvduzN5NHPcSO2BjwyZsaYrxhjWoHvAV+NXfcFXstEU/5S/prouKLXm22M2YH9Gb0mThs38oSXxiw6t3pizBKVwxJdfJkxllnjbJOo5fHWEbvcirl1M65+oNKyrCXABmC1MaYgTtwTFVssE3Mb3TaZYxYvLoBK4BHsbxnfN8ZUJzsuY8znsf/ZPBqnfexz3IhtY9RiV8fMsqwfWpZVDXwD+OYF1p1syl/KXxMdV/R6Ky37yuo/Aq6N+iwmNbYx8oQnxixObnV9zBKVwxJdfLUDM6PuV2AfGDdmG2OMHygETp3nufGWnwSKnD5i1xW7jmygNGr5SNvzPScZcRUCTcAMy7K6nDYh4AQwL8ljFmukD7fHLF5cWJY1Ek8TsAV7E3XS4jLG3Ag8DNxuWdaQ0092VF8V2B/K6NeS9NiAg14Zsyg/BUa2vI0nbySD8pfy10THNbpe57MI9nFCzfzps5i02OLkCdfHLE5u9cSYRbm0HHa+fZIXO2HvZ23DPuBs5AC2+TFtvkLMwaTO/HzOPWCtDfuAuLh9Ai9y7kFx98dZx4tOH9EHrK4eWYeLcb3g9HMQqHb62Yd9pklxMscsal1PYx8zMdLHYuyDyPdi7/s+MBKby3EtAvKdPlZifxDrk/heLsE+qLIm5u//IPbm51JnzNqT/V7GxuahMYseq9uAHedbRyJzk/KX8hceyF8xeaLW6Wcv9kkL9ckcM+LnCVfHLDYuj41ZwnLYRCSwW4D3ncF72Fn2HewKFuxvcS9in5L5R6Aq6rkPO8/bD9x8vj6d5VVOHy1On1nx1hHVR5cz7Qf+2+24nOX/gn367BB2tXybC2N2NXah0OeMT0NUH53Oshbglx6J67AzXsewT/3dmOS4NmH/k9lF1GUbnH6OA0En3vtceC8/FJtHxmy98/7tAn5HVGETbx3JnsZ6TS68f8pfl3H+cpb/fdSYHcc+hskrecLtMYuXW70wZgnLYbrCvYiIiEgS6Qr3IiIiIkmk4ktEREQkiVR8iYiIiCSRii8RERGRJFLxJSIiIpJEKr5EREREkkjFl4iIiEgSqfgSERERSaL/B/54p9xo/x/NAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x360 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize = (10,5))\n",
    "plt.subplot(121)\n",
    "plt.plot(result.times, np.real([a[0,0] for a in result.states]), label = 'GG')\n",
    "plt.plot(result.times, np.real([a[1,1] for a in result.states]), label = 'EE')\n",
    "plt.plot(result.times, np.real([a[2,2] for a in result.states]), label = 'FF')\n",
    "plt.title(\"Populations\")\n",
    "plt.legend(loc = 2)\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.plot(result.times, np.abs([a[0,1] for a in result.states]), label = 'GE')\n",
    "plt.plot(result.times, np.abs([a[1,2] for a in result.states]), label = 'EF')\n",
    "plt.plot(result.times, np.real([a[0,2] for a in result.states]), label = 'GF')\n",
    "plt.title(\"Coherences\")\n",
    "plt.legend(loc = 2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ZZ-CSUM with the original (4 EF pi pulses) method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first define the entangling ZZ Hamiltonian between our qutrits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha_11 = 2 * np.pi * -0.27935*1e6 # Hz\n",
    "alpha_12 = 2 * np.pi * 0.1599*1e6 # Hz\n",
    "alpha_21 = 2 * np.pi * -0.52793*1e6 # Hz\n",
    "alpha_22 = 2 * np.pi * -0.742967*1e6 # Hz\n",
    "\n",
    "H_ZZ = qt.Qobj(np.diag([0,0,0,0,alpha_11, alpha_12, 0, alpha_21, alpha_22]), dims = [[3,3],[3,3]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We next calculate the required times for these to be turned into a CSUM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "T1 = 0.614 * 1e-6 #s\n",
    "T2 = 0.105 * 1e-6 #s\n",
    "T3 = 0.614 * 1e-6 #s\n",
    "T4 = 0.105 * 1e-6 #s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The CSUM gate uses some single-qutrit unitaries, which we define below. In particular, we use the Hadamard gate and $\\pi$-pulses on the $|E\\rangle\\leftrightarrow |F\\rangle$ subspace of each qutrit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "omega = np.exp(2*np.pi*1j/3)\n",
    "omega_inv = np.exp(-2*np.pi*1j/3)\n",
    "Hadamard = qt.Qobj(1/np.sqrt(3)*np.array([[1,1,1],\n",
    "                                         [1, omega, omega_inv],\n",
    "                                         [1, omega_inv, omega]]))\n",
    "PiEF = qt.Qobj(np.array([[1,0,0],\n",
    "                        [0,0,-1j],\n",
    "                        [0,-1j,0]]))\n",
    "\n",
    "PiEFm = qt.Qobj(np.array([[1,0,0],\n",
    "                        [0,0,1j],\n",
    "                        [0,1j,0]]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "Hadamard_Q1 = qt.tensor(Id, Hadamard)\n",
    "Hadamard_Q1_dag = Hadamard_Q1.dag()\n",
    "\n",
    "PiEF_Q0 = qt.tensor(PiEF, Id)\n",
    "PiEF_Q0_dag = PiEF_Q0.dag()\n",
    "\n",
    "PiEFm_Q0 = qt.tensor(PiEFm, Id)\n",
    "PiEFm_Q0_dag = PiEFm_Q0.dag()\n",
    "\n",
    "PiEF_Q1 = qt.tensor(Id, PiEF)\n",
    "PiEF_Q1_dag = PiEF_Q1.dag()\n",
    "\n",
    "PiEFm_Q1 = qt.tensor(Id,PiEFm)\n",
    "PiEFm_Q1_dag = PiEFm_Q1.dag()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fidelity_of_csum():\n",
    "    \n",
    "    gamma1GE_0 = 1./(40e-6) # 1/s\n",
    "    gamma1EF_0 = 1./(26e-6) # 1/s\n",
    "\n",
    "    C0 = np.sqrt(gamma1GE_0)*qt.tensor(qt.Qobj(np.array([[0,1,0],[0,0,0],[0,0,0]])), Id)\n",
    "    C1 = np.sqrt(gamma1EF_0)*qt.tensor(qt.Qobj(np.array([[0,0,0],[0,0,1],[0,0,0]])), Id)\n",
    "\n",
    "    gammaphiGE_0 = 1/(50e-6) # 1/s\n",
    "    gammaphiEF_0 = 1/(45e-6) # 1/s\n",
    "\n",
    "    D0 = np.sqrt(gammaphiGE_0) * qt.tensor(qt.Qobj(np.diag([1,-1,-1])), Id)\n",
    "    D1 = np.sqrt(gammaphiEF_0) * qt.tensor(qt.Qobj(np.diag([-1,1,-1])), Id)\n",
    "\n",
    "    gamma1GE_1 = 1./(50e-6) # 1/s\n",
    "    gamma1EF_1 = 1./(33e-6) # 1/s\n",
    "\n",
    "    E0 = np.sqrt(gamma1GE_1)*qt.tensor(Id, qt.Qobj(np.array([[0,1,0],[0,0,0],[0,0,0]])))\n",
    "    E1 = np.sqrt(gamma1EF_1)*qt.tensor(Id, qt.Qobj(np.array([[0,0,0],[0,0,1],[0,0,0]])))\n",
    "\n",
    "    gammaphiGE_1 = 1/(60e-6) # 1/s\n",
    "    gammaphiEF_1 = 1/(35e-6) # 1/s\n",
    "\n",
    "    F0 = np.sqrt(gammaphiGE_1) * qt.tensor(Id, qt.Qobj(np.diag([1,-1,-1])))\n",
    "    F1 = np.sqrt(gammaphiEF_1) * qt.tensor(Id, qt.Qobj(np.diag([-1,1,-1])))\n",
    "    \n",
    "    fidelities = []\n",
    "    \n",
    "    for state_0 in range(3):\n",
    "        for state_1 in range(3):\n",
    "            state = [state_0, state_1]\n",
    "            \n",
    "            psi0 = qt.tensor(qt.basis(3,state[0])*qt.basis(3,state[0]).dag(),qt.basis(3,state[1])*qt.basis(3,state[1]).dag())\n",
    "\n",
    "            # now do the CSUM\n",
    "            psi0 = Hadamard_Q1_dag*psi0*Hadamard_Q1\n",
    "\n",
    "            result = qt.mesolve(H_ZZ, psi0, np.linspace(0,T1,100), [C0, C1, D0, D1, E0, E1, F0, F1], [])\n",
    "            psi0 = result.states[-1]\n",
    "\n",
    "            psi0 = PiEF_Q0*psi0*PiEF_Q0_dag\n",
    "\n",
    "            result = qt.mesolve(H_ZZ, psi0, np.linspace(0,T2,100), [C0, C1, D0, D1, E0, E1, F0, F1], [])\n",
    "            psi0 = result.states[-1]\n",
    "\n",
    "            psi0 = PiEF_Q1*psi0*PiEF_Q1_dag\n",
    "\n",
    "            result = qt.mesolve(H_ZZ, psi0, np.linspace(0,T3,100), [C0, C1, D0, D1, E0, E1, F0, F1], [])\n",
    "\n",
    "            psi0 = result.states[-1]\n",
    "            psi0 = PiEFm_Q0*psi0*PiEFm_Q0_dag\n",
    "\n",
    "            result = qt.mesolve(H_ZZ, psi0, np.linspace(0,T4,100), [C0, C1, D0, D1, E0, E1, F0, F1], [])\n",
    "\n",
    "            psi0 = result.states[-1]\n",
    "            psi0 = PiEFm_Q1*psi0*PiEFm_Q1_dag\n",
    "\n",
    "            psi0 = Hadamard_Q1*psi0*Hadamard_Q1_dag\n",
    "\n",
    "            fidelities.append(np.abs(psi0[3* state_0 + ((state_0 + state_1)%3), 3* state_0 + ((state_0 + state_1)%3)]))\n",
    "    \n",
    "    return fidelities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'fidelity_of_csum' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-2-e7b376140e70>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mfidelity_of_csum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'fidelity_of_csum' is not defined"
     ]
    }
   ],
   "source": [
    "fidelity_of_csum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9025060315655272"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(fidelity_of_csum())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ZZ-CSUM with the shortest time "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first define the entangling ZZ Hamiltonian between our qutrits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha_11 = 2 * np.pi * -0.27935 * 1e6 # Hz\n",
    "alpha_12 = 2 * np.pi * 0.1599 * 1e6 # Hz\n",
    "alpha_21 = 2 * np.pi * -0.52793 * 1e6 # Hz\n",
    "alpha_22 = 2 * np.pi * -0.742967 * 1e6 # Hz\n",
    "\n",
    "H_ZZ = qt.Qobj(np.diag([0,0,0,0,alpha_11, alpha_12, 0, alpha_21, alpha_22]), dims = [[3,3],[3,3]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We next calculate the required times for these to be turned into a CSUM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "T1 = 0.3986566 * 1e-6 #s\n",
    "T2 = 0.3986566 * 1e-6  #s\n",
    "T3 = 0.13608932 * 1e-6 #s\n",
    "T4 = 0.13608932 * 1e-6 #s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The CSUM gate uses some single-qutrit unitaries, which we define below. In particular, we use the Hadamard gate and $\\pi$-pulses on the $|E\\rangle\\leftrightarrow |F\\rangle$ subspace of each qutrit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "omega = np.exp(2*np.pi*1j/3)\n",
    "omega_inv = np.exp(-2*np.pi*1j/3)\n",
    "Hadamard = 1/np.sqrt(3)*np.array([[1,1,1],\n",
    "                                         [1, omega, omega_inv],\n",
    "                                         [1, omega_inv, omega]])\n",
    "PiEF = np.array([[1,0,0],\n",
    "                [0,0,-1j],\n",
    "                [0,-1j,0]])\n",
    "\n",
    "PiEFm = np.array([[1,0,0],\n",
    "                [0,0,1j],\n",
    "                [0,1j,0]])\n",
    "\n",
    "PiGE = np.array([[0,-1j,0],\n",
    "                [-1j,0,0],\n",
    "                [0,0,1]])\n",
    "\n",
    "PiGEm = np.array([[0,1j,0],\n",
    "                [1j,0,0],\n",
    "                [0,0,1]])\n",
    "\n",
    "Id = np.eye(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "PiEF_Q0 = np.kron(PiEF, Id)\n",
    "\n",
    "PiEFm_Q0 = np.kron(PiEFm, Id)\n",
    "\n",
    "PiEF_Q1 = np.kron(Id, PiEF)\n",
    "\n",
    "PiEFm_Q1 = np.kron(Id,PiEFm)\n",
    "\n",
    "PiGE_Q0 = np.kron(PiGE, Id)\n",
    "\n",
    "PiGEm_Q0 = np.kron(PiGEm, Id)\n",
    "\n",
    "PiGE_Q1 = np.kron(Id, PiGE)\n",
    "\n",
    "PiGEm_Q1 = np.kron(Id,PiGEm)\n",
    "\n",
    "H_ZZ = np.diag([0,0,0,0,alpha_11, alpha_12, 0, alpha_21, alpha_22])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "## single-qutrit permutation matrices\n",
    "P = {'P0': np.array([[1,0,0],\n",
    "                     [0,1,0],\n",
    "                     [0,0,1]]),\n",
    "     'P1': np.array([[1,0,0],\n",
    "                     [0,0,1],\n",
    "                     [0,1,0]]),\n",
    "     'P2': np.array([[0,1,0],\n",
    "                     [1,0,0],\n",
    "                     [0,0,1]]),\n",
    "     'P3': np.array([[0,1,0],\n",
    "                     [0,0,1],\n",
    "                     [1,0,0]]),\n",
    "     'P4': np.array([[0,0,1],\n",
    "                     [1,0,0],\n",
    "                     [0,1,0]]),\n",
    "     'P5': np.array([[0,0,1],\n",
    "                     [0,1,0],\n",
    "                     [1,0,0]]),\n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 1, 0],\n",
       "       [0, 0, 1],\n",
       "       [1, 0, 0]])"
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(P['P0'].T,P['P3'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.+0.j,  0.+0.j,  0.+0.j],\n",
       "       [ 0.+0.j, -1.+0.j,  0.+0.j],\n",
       "       [ 0.+0.j,  0.+0.j,  1.+0.j]])"
      ]
     },
     "execution_count": 100,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reduce(np.dot, reversed([PiGE, PiEF, PiGE, PiEF, PiEF, PiGE, PiEF, PiGE]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 0, 0],\n",
       "       [0, 0, 1],\n",
       "       [0, 1, 0]])"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "PiEF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "CSUM_permutation_1 = np.dot(np.kron(P['P1'].T,P['P4'].T), np.kron(P['P0'], P['P0']))\n",
    "CSUM_permutation_2 = np.dot(np.kron(P['P4'].T,P['P5'].T), np.kron(P['P1'], P['P4']))\n",
    "CSUM_permutation_3 = np.dot(np.kron(P['P5'].T,P['P3'].T), np.kron(P['P4'], P['P5']))\n",
    "CSUM_permutation_4 = np.dot(np.kron(P['P3'].T,P['P2'].T), np.kron(P['P5'], P['P3']))\n",
    "CSUM_permutation_5 = np.dot(np.kron(P['P0'].T,P['P0'].T), np.kron(P['P3'], P['P2']))\n",
    "\n",
    "Phi_00 = -1*(alpha_22 * T2 + alpha_11 * T3 + alpha_21 * T4)\n",
    "Phi_01 = -1*(alpha_21 * T2 + alpha_22 * T3)\n",
    "Phi_02 = -1*(alpha_12 * T4)\n",
    "Phi_10 = -1*(alpha_22 * T1 + alpha_11 * T3 + alpha_21 * T4)\n",
    "Phi_20 = -1*(alpha_12 * T1 + alpha_12 * T2)\n",
    "\n",
    "Z_corr = np.kron(expm(-1j*np.diag([0,Phi_10 - Phi_00, Phi_20 - Phi_00])),\n",
    "                   expm(-1j*np.diag([0,Phi_01 - Phi_00, Phi_02 - Phi_00])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "Hadamard_Q1 = np.kron(Id, Hadamard)\n",
    "Hadamard_Q1_dag = np.conj(Hadamard_Q1.T)\n",
    "\n",
    "Hadamard_Q0 = np.kron(Hadamard, Id)\n",
    "Hadamard_Q0_dag = np.conj(Hadamard_Q0.T)\n",
    "\n",
    "Hadamard_Q0_map = unitary_op_to_map(Hadamard_Q0)\n",
    "Hadamard_Q0_dag_map = unitary_op_to_map(Hadamard_Q0_dag)\n",
    "\n",
    "Hadamard_Q1_map = unitary_op_to_map(Hadamard_Q1)\n",
    "Hadamard_Q1_dag_map = unitary_op_to_map(Hadamard_Q1_dag)\n",
    "\n",
    "Z_corr_map = unitary_op_to_map(Z_corr)\n",
    "\n",
    "Perm1_map = unitary_op_to_map(CSUM_permutation_1)\n",
    "Perm2_map = unitary_op_to_map(CSUM_permutation_2)\n",
    "Perm3_map = unitary_op_to_map(CSUM_permutation_3)\n",
    "Perm4_map = unitary_op_to_map(CSUM_permutation_4)\n",
    "Perm5_map = unitary_op_to_map(CSUM_permutation_5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'idling_Lindbladian' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-16-882438718ca2>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0midling_map_T1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexpm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0midling_Lindbladian\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mT1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0midling_map_T2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexpm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0midling_Lindbladian\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mT2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0midling_map_T3\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexpm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0midling_Lindbladian\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mT3\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0midling_map_T4\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexpm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0midling_Lindbladian\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mT4\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'idling_Lindbladian' is not defined"
     ]
    }
   ],
   "source": [
    "idling_map_T1 = sc.linalg.expm(idling_Lindbladian*T1)\n",
    "idling_map_T2 = sc.linalg.expm(idling_Lindbladian*T2)\n",
    "idling_map_T3 = sc.linalg.expm(idling_Lindbladian*T3)\n",
    "idling_map_T4 = sc.linalg.expm(idling_Lindbladian*T4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 378,
   "metadata": {},
   "outputs": [],
   "source": [
    "real_CSUMp1_map = reduce(np.dot, [Hadamard_Q1_map, Z_corr_map, Perm5_map, idling_map_T4, Perm4_map, idling_map_T3, Perm3_map, idling_map_T2, Perm2_map, idling_map_T1, Perm1_map, Hadamard_Q1_dag_map])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 421,
   "metadata": {},
   "outputs": [],
   "source": [
    "init_state = (2,2)\n",
    "psi_init = [0,0,0,0,0,0,0,0,0]\n",
    "psi_init[3*init_state[0] + init_state[1]] = 1\n",
    "\n",
    "rho_init = ST.vec_to_dm(psi_init)\n",
    "rho_init_vec = rho_init.T.reshape(rho_init.shape[0]*rho_init.shape[1],-1)\n",
    "\n",
    "rho_final_vec = np.dot(real_CSUMp1_map, rho_init_vec)\n",
    "\n",
    "rho_final = rho_final_vec.reshape(9,9).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 422,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       "       [0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       "       [0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       "       [0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       "       [0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       "       [0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       "       [0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       "       [0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       "       [0, 0, 0, 0, 0, 0, 0, 0, 1]])"
      ]
     },
     "execution_count": 422,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rho_init"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 423,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.001, 0.   , 0.001, 0.   , 0.   , 0.   , 0.   , 0.   , 0.   ],\n",
       "       [0.   , 0.001, 0.001, 0.   , 0.   , 0.   , 0.   , 0.   , 0.   ],\n",
       "       [0.001, 0.001, 0.009, 0.   , 0.   , 0.   , 0.   , 0.   , 0.   ],\n",
       "       [0.   , 0.   , 0.   , 0.002, 0.003, 0.001, 0.   , 0.   , 0.   ],\n",
       "       [0.   , 0.   , 0.   , 0.003, 0.006, 0.002, 0.   , 0.   , 0.   ],\n",
       "       [0.   , 0.   , 0.   , 0.001, 0.002, 0.001, 0.   , 0.   , 0.   ],\n",
       "       [0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.026, 0.002, 0.007],\n",
       "       [0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.002, 0.929, 0.002],\n",
       "       [0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.007, 0.002, 0.026]])"
      ]
     },
     "execution_count": 423,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.round(np.abs(rho_final),3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 429,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9268593749999999"
      ]
     },
     "execution_count": 429,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "0.975*0.975*0.975"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "00: 92.5\n",
    "01: 92.5\n",
    "02: 92.5\n",
    "10: 92.5\n",
    "11: 92.5\n",
    "12: 92.5\n",
    "20: 92.9\n",
    "21: 92.9\n",
    "22: 92.9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_csum(state_0 = 0, state_1 = 0):\n",
    "    \n",
    "    fidelities = []\n",
    "    \n",
    "    psi0 = np.array([0,0,0,0,0,0,0,0,0])\n",
    "    psi0[state_0*3 + state_1] = 1\n",
    "    \n",
    "    psi0 = np.dot(np.kron(np.eye(3), np.conj(Hadamard.T)), psi0)\n",
    "    \n",
    "    # do the first permutation\n",
    "    psi0 = reduce(np.dot, [CSUM_permutation_1,  psi0]) \n",
    "\n",
    "    psi0 = np.dot(expm(-1j*H_ZZ*T1),psi0)\n",
    "\n",
    "    psi0 = reduce(np.dot, [CSUM_permutation_2, psi0])\n",
    "\n",
    "    psi0 = np.dot(expm(-1j*H_ZZ*T2),psi0)\n",
    "\n",
    "    psi0 = reduce(np.dot, [CSUM_permutation_3, psi0])\n",
    " \n",
    "    psi0 = np.dot(expm(-1j*H_ZZ*T3),psi0)            \n",
    "\n",
    "    psi0 = reduce(np.dot, [CSUM_permutation_4, psi0])\n",
    "\n",
    "    psi0 = np.dot(expm(-1j*H_ZZ*T4),psi0)\n",
    "            \n",
    "    psi0 = reduce(np.dot, [CSUM_permutation_5, psi0])\n",
    "   \n",
    "    Phi_00 = -1*(alpha_22 * T2 + alpha_11 * T3 + alpha_21 * T4)\n",
    "    Phi_01 = -1*(alpha_21 * T2 + alpha_22 * T3)\n",
    "    Phi_02 = -1*(alpha_12 * T4)\n",
    "    Phi_10 = -1*(alpha_22 * T1 + alpha_11 * T3 + alpha_21 * T4)\n",
    "    Phi_20 = -1*(alpha_12 * T1 + alpha_12 * T2)\n",
    "\n",
    "    Z_corr = np.kron(expm(-1j*np.diag([0,Phi_10 - Phi_00, Phi_20 - Phi_00])),\n",
    "                       expm(-1j*np.diag([0,Phi_01 - Phi_00, Phi_02 - Phi_00])))\n",
    "    \n",
    "    psi0 = reduce(np.dot, [Z_corr, psi0])\n",
    "    psi0 = psi0*exp(-1j*Phi_00)\n",
    "    \n",
    "    psi0 = np.dot(np.kron(np.eye(3), Hadamard), psi0)\n",
    "    \n",
    "    print(np.abs(psi0))\n",
    "\n",
    "    return psi0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.        +0.j        ,  0.        +0.j        ,\n",
       "         0.        +0.j        ],\n",
       "       [ 0.        +0.j        ,  0.82891884+0.55936889j,\n",
       "         0.        +0.j        ],\n",
       "       [-0.        +0.j        , -0.        +0.j        ,\n",
       "        -0.89888709+0.43818032j]])"
      ]
     },
     "execution_count": 105,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expm(-1j*np.diag([0,Phi_01 - Phi_00, Phi_02 - Phi_00]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "120.00000044742617"
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.angle(-0.89888709+0.43818032j)*180/np.pi - np.angle(0.82891884+0.55936889j)*180/np.pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.00000000e+00 2.49800181e-16 1.94289029e-16 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00]\n"
     ]
    }
   ],
   "source": [
    "test_csum(0,0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Direct calculation of the two-qutrit quantum map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To simulate the performance of the ZZ-based CSUM gate in the scrambling algorithm, we extract here the quantum map by directly (numerically) exponentiating the Lindbladian.  We'll compare this with the results obtained by solving the master equation above.  This calculation does not rely on QuTiP and thus serves as an independent check of the dynamics.\n",
    "\n",
    "Following [this paper](https://arxiv.org/pdf/1510.08634.pdf), if we vectorize the density matrix by simply ordering the columns one below the other, so that the $(a,b)$ entry of $\\rho$ is the $(bn + a)$$^{th}$ entry (zero-indexed) of the vector $\\vec{r}$, and if we write the master equation in standard form as \n",
    "\n",
    "\\begin{equation}\n",
    "\\frac { d } { d t } \\hat { \\boldsymbol { \\rho } } ( t ) = \\mathcal { L } \\hat { \\boldsymbol { \\rho } } = -  i   [ \\hat { \\mathbf { H } } , \\hat { \\boldsymbol { \\rho } } ] + \\sum _ { i } \\gamma _ { i } \\left( \\hat { \\mathbf { A } } _ { i } \\hat { \\boldsymbol { \\rho } } \\hat { \\mathbf { A } } _ { i } ^ { \\dagger } - \\frac { 1 } { 2 } \\left\\{ \\hat { \\mathbf { A } } _ { i } ^ { \\dagger } \\hat { \\mathbf { A } } _ { i } , \\hat { \\boldsymbol { \\rho } } \\right\\} \\right)\n",
    "\\end{equation},\n",
    "\n",
    "then the vectorized versions of the operators are:\n",
    "\n",
    "\\begin{equation}\n",
    "[ \\hat { \\mathbf { H } } , \\hat { \\boldsymbol { \\rho } } ] \\rightarrow  \\left( I \\otimes \\hat { \\mathbf { H } } - \\hat { \\mathbf { H } } ^ { T } \\otimes I \\right) \\vec { r }\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\hat { \\mathbf { A } } _ { i } \\hat { \\boldsymbol { \\rho } } \\hat { \\mathbf { A } } _ { i } ^ { \\dagger } \\rightarrow \\left( \\left( \\hat { \\mathbf { A } } _ { i } ^ { \\dagger } \\right) ^ { T } \\otimes \\hat { \\mathbf { A } } _ { i } \\right) \\vec { r }\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\hat { \\mathbf { A } } _ { i } ^ { \\dagger } \\mathbf { \\hat { A } } _ { i } \\hat { \\boldsymbol { \\rho } } \\rightarrow \\left( I \\otimes \\hat { \\mathbf { A } } _ { i } ^ { \\dagger } \\hat { \\mathbf { A } } _ { i } \\right) \\vec { r }\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\hat { \\boldsymbol { \\rho } } \\hat { \\mathbf { A } } _ { i } ^ { \\dagger } \\hat { \\mathbf { A } } _ { i } \\rightarrow \\left( \\left( \\hat { \\mathbf { A } } _ { i } ^ { \\dagger } \\hat { \\mathbf { A } } _ { i } \\right) ^ { T } \\otimes I \\right) \\vec { r }\n",
    "\\end{equation}\n",
    "\n",
    "That means the whole Lindbladian can be expressed as:\n",
    "\n",
    "\\begin{equation}\n",
    "L = -i \\left( I \\otimes \\hat { \\mathbf { H } } - \\hat { \\mathbf { H } } ^ { T } \\otimes I \\right)+ \\sum _ { i } \\gamma _ { i } \\left( \\left( \\hat { \\mathbf { A } } _ { i } ^ { \\dagger } \\right) ^ { T } \\otimes \\hat { \\mathbf { A } } _ { i } - \\frac { 1 } { 2 } \\left( I \\otimes \\hat { \\mathbf { A } } _ { i } ^ { \\dagger } \\hat { \\mathbf { A } } _ { i } + \\left( \\hat { \\mathbf { A } } _ { i } ^ { \\dagger } \\hat { \\mathbf { A } } _ { i } \\right) ^ { T } \\otimes I \\right) \\right)\n",
    "\\end{equation}\n",
    "\n",
    "In terms of which the Master equation becomes:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac { d } { d t } \\vec { r } = L \\vec { r }\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Take a look at the code here for matrix_Lindbladian, a function which implements the above calculations given a Hamiltonian, collapse rates, and collapse operators:*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 358,
   "metadata": {},
   "outputs": [],
   "source": [
    "def matrix_Lindbladian(Hamiltonian, collapse_rates, collapse_operators):\n",
    "    \"\"\"\n",
    "    Forms a Lindbladian superoperator in 'vectorized' form, \n",
    "    i.e. a matrix which acts on a 'vectorized' density matrix\n",
    "    \n",
    "    Arguments:\n",
    "    Hamiltonian - an n x n matrix (if the Hilbert space dimension is n) \n",
    "                    describing the coherent part of the system's evolution\n",
    "    collapse_rates - a list of numbers, each of which corresponds to a collapse_operator in \n",
    "                    the list collapse_operators.  \n",
    "    collapse_operators - a list of operators which act as collapse operators in the \n",
    "                    master equation describing the evolution\n",
    "                    \n",
    "    Returns:\n",
    "    Lindbladian - an n^2 x n^2 matrix representing the Lindbladian superoperator in vectorized form\n",
    "                The convention is that this acts on a density matrix which has been vectorized by \n",
    "                stacking columns one on top of the other. \n",
    "    \"\"\"\n",
    "    \n",
    "    assert len(collapse_rates) == len(collapse_operators)\n",
    "    \n",
    "    system_dimension = Hamiltonian.shape\n",
    "    \n",
    "    identity_matrix = np.eye(system_dimension[0])\n",
    "    \n",
    "    # start by vectorizing the Hamiltonian part of the Lindbladian\n",
    "    Lindbladian = -1j * (np.kron(identity_matrix, Hamiltonian) - np.kron(Hamiltonian.T, identity_matrix) )\n",
    "    \n",
    "    # then add each collapse term, again in vectorized form\n",
    "    for collapse_index, collapse_rate in enumerate(collapse_rates):\n",
    "        collapse_operator = collapse_operators[collapse_index]\n",
    "        \n",
    "        collapse_term = collapse_rate * (np.kron(np.conj(collapse_operator), collapse_operator) - \\\n",
    "                                        0.5 * np.kron(identity_matrix, np.dot(np.conj(collapse_operator.T), collapse_operator)) - \\\n",
    "                                        0.5 * np.kron(np.dot(np.conj(collapse_operator.T), collapse_operator), identity_matrix))\n",
    "        \n",
    "        Lindbladian = Lindbladian + collapse_term\n",
    "    \n",
    "    return Lindbladian\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 366,
   "metadata": {},
   "outputs": [],
   "source": [
    "Id = np.eye(3)\n",
    "\n",
    "H_ZZ = np.diag([0,0,0,0,alpha_11, alpha_12, 0, alpha_21, alpha_22])\n",
    "\n",
    "gamma1GE_0 = 1./(40e-6) # 1/s\n",
    "gamma1EF_0 = 1./(26e-6) # 1/s\n",
    "\n",
    "C0 = np.kron(np.array([[0,1,0],[0,0,0],[0,0,0]]), Id)\n",
    "C1 = np.kron(np.array([[0,0,0],[0,0,1],[0,0,0]]), Id)\n",
    "\n",
    "gammaphiGE_0 = 1/(50e-6) # 1/s\n",
    "gammaphiEF_0 = 1/(45e-6) # 1/s\n",
    "\n",
    "D0 = np.kron(np.diag([1,-1,-1]), Id)\n",
    "D1 = np.kron(np.diag([-1,1,-1]), Id)\n",
    "\n",
    "gamma1GE_1 = 1./(50e-6) # 1/s\n",
    "gamma1EF_1 = 1./(33e-6) # 1/s\n",
    "\n",
    "E0 = np.kron(Id, np.array([[0,1,0],[0,0,0],[0,0,0]]))\n",
    "E1 = np.kron(Id, np.array([[0,0,0],[0,0,1],[0,0,0]]))\n",
    "\n",
    "gammaphiGE_1 = 1/(60e-6) # 1/s\n",
    "gammaphiEF_1 = 1/(35e-6) # 1/s\n",
    "\n",
    "F0 = np.kron(Id, np.diag([1,-1,-1]))\n",
    "F1 = np.kron(Id, np.diag([-1,1,-1]))\n",
    "    \n",
    "qutrit_collapse_rates = [gamma1GE_0, gamma1EF_0, gammaphiGE_0, gammaphiEF_0, gamma1GE_1, gamma1EF_1, gammaphiGE_1, gammaphiEF_1]\n",
    "qutrit_collapse_ops = [C0, C1, D0, D1, E0, E1, F0, F1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 367,
   "metadata": {},
   "outputs": [],
   "source": [
    "idling_Lindbladian = matrix_Lindbladian(H_ZZ, qutrit_collapse_rates, qutrit_collapse_ops)\n",
    "nodecay_Lindbladian = matrix_Lindbladian(H_ZZ, [0]*8, qutrit_collapse_ops)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 394,
   "metadata": {},
   "outputs": [],
   "source": [
    "idling_map_T1 = sc.linalg.expm(idling_Lindbladian*T1)\n",
    "idling_map_T2 = sc.linalg.expm(idling_Lindbladian*T2)\n",
    "idling_map_T3 = sc.linalg.expm(idling_Lindbladian*T3)\n",
    "idling_map_T4 = sc.linalg.expm(idling_Lindbladian*T4)\n",
    "\n",
    "idling_nodecay_map_T1 = sc.linalg.expm(nodecay_Lindbladian*T1)\n",
    "idling_nodecay_map_T2 = sc.linalg.expm(nodecay_Lindbladian*T2)\n",
    "idling_nodecay_map_T3 = sc.linalg.expm(nodecay_Lindbladian*T3)\n",
    "idling_nodecay_map_T4 = sc.linalg.expm(nodecay_Lindbladian*T4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need the vectorized map representations of the unitary operations representing the $\\pi$ pulses and the Hadamard gate.  From the considerations above, namely the equation\n",
    "\n",
    "\\begin{equation}\n",
    "\\hat { \\mathbf { A } } _ { i } \\hat { \\boldsymbol { \\rho } } \\hat { \\mathbf { A } } _ { i } ^ { \\dagger } \\rightarrow \\left( \\left( \\hat { \\mathbf { A } } _ { i } ^ { \\dagger } \\right) ^ { T } \\otimes \\hat { \\mathbf { A } } _ { i } \\right) \\vec { r }\n",
    "\\end{equation}\n",
    "\n",
    "we gather that a unitary operation $\\hat U$ is represented as $\\hat U^* \\otimes \\hat U$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*This is implemented in ST.unitary_op_to_map(), using the code below*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def unitary_op_to_map(unitary_op):\n",
    "    \"\"\"\n",
    "    Given a unitary operator, converts it into vectorized form.\n",
    "\n",
    "    that is, the n x n matrix unitary_op is defined such that it transforms density matrices via\n",
    "    rho_out = unitary_op * rho_in * unitary_op.dagger\n",
    "    \n",
    "    The vectorized version is defined such that it acts on the vectorized density matrix rho_in_vec by\n",
    "    rho_out_vec = vectorized_unitary_op * rho_in_vec\n",
    "\n",
    "    The vectorized version of a unitary operator is given by np.kron(unitary_op*, unitary_op) \n",
    "    [see https://arxiv.org/pdf/1510.08634.pdf (page 4) for a proof]\n",
    "    \"\"\"\n",
    "    return np.kron(np.conj(unitary_op), unitary_op)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 399,
   "metadata": {},
   "outputs": [],
   "source": [
    "omega = np.exp(2*np.pi*1j/3)\n",
    "omega_inv = np.exp(-2*np.pi*1j/3)\n",
    "Hadamard = 1/np.sqrt(3)*np.array([[1,1,1],\n",
    "                                 [1, omega, omega_inv],\n",
    "                                 [1, omega_inv, omega]])\n",
    "PiEF = np.array([[1,0,0],\n",
    "                [0,0,-1j],\n",
    "                [0,-1j,0]])\n",
    "\n",
    "PiEFm = np.array([[1,0,0],\n",
    "                [0,0,1j],\n",
    "                [0,1j,0]])\n",
    "\n",
    "Hadamard_Q1 = np.kron(Id, Hadamard)\n",
    "Hadamard_Q1_dag = np.conj(Hadamard_Q1.T)\n",
    "\n",
    "Hadamard_Q0 = np.kron(Hadamard, Id)\n",
    "Hadamard_Q0_dag = np.conj(Hadamard_Q0.T)\n",
    "\n",
    "PiEF_Q0 = np.kron(PiEF, Id)\n",
    "PiEFm_Q0 = np.kron(PiEFm, Id)\n",
    "PiEF_Q1 = np.kron(Id, PiEF)\n",
    "PiEFm_Q1 = np.kron(Id,PiEFm)\n",
    "\n",
    "PiEF_Q0_map = ST.unitary_op_to_map(PiEF_Q0)\n",
    "PiEFm_Q0_map = ST.unitary_op_to_map(PiEFm_Q0)\n",
    "PiEF_Q1_map = ST.unitary_op_to_map(PiEF_Q1)\n",
    "PiEFm_Q1_map = ST.unitary_op_to_map(PiEFm_Q1)\n",
    "Hadamard_Q1_map = ST.unitary_op_to_map(Hadamard_Q1)\n",
    "Hadamard_Q1_dag_map = ST.unitary_op_to_map(Hadamard_Q1_dag)\n",
    "\n",
    "Hadamard_Q0_map = ST.unitary_op_to_map(Hadamard_Q0)\n",
    "Hadamard_Q0_dag_map = ST.unitary_op_to_map(Hadamard_Q0_dag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 318,
   "metadata": {},
   "outputs": [],
   "source": [
    "real_CSUMp1_map = reduce(np.dot, [Hadamard_Q1_map, PiEFm_Q1_map, idling_map_T4, PiEFm_Q0_map, idling_map_T3, PiEF_Q1_map, idling_map_T2, PiEF_Q0_map, idling_map_T1, Hadamard_Q1_dag_map])\n",
    "real_CMINp1_map = reduce(np.dot, [Hadamard_Q1_map, PiEFm_Q1_map, idling_map_T3, PiEFm_Q0_map, idling_map_T4, PiEF_Q1_map, idling_map_T1, PiEF_Q0_map, idling_map_T2, Hadamard_Q1_dag_map])\n",
    "\n",
    "ideal_CSUMp1_map = reduce(np.dot, [Hadamard_Q1_map, PiEFm_Q1_map, idling_nodecay_map_T4, PiEFm_Q0_map, idling_nodecay_map_T3, PiEF_Q1_map, idling_nodecay_map_T2, PiEF_Q0_map, idling_nodecay_map_T1, Hadamard_Q1_dag_map])\n",
    "ideal_CMINp1_map = reduce(np.dot, [Hadamard_Q1_map, PiEFm_Q1_map, idling_nodecay_map_T3, PiEFm_Q0_map, idling_nodecay_map_T4, PiEF_Q1_map, idling_nodecay_map_T1, PiEF_Q0_map, idling_nodecay_map_T2, Hadamard_Q1_dag_map])\n",
    "\n",
    "real_CSUMm1_map = reduce(np.dot, [Hadamard_Q0_map, PiEFm_Q1_map, idling_map_T4, PiEFm_Q0_map, idling_map_T3, PiEF_Q1_map, idling_map_T2, PiEF_Q0_map, idling_map_T1, Hadamard_Q0_dag_map])\n",
    "real_CMINm1_map = reduce(np.dot, [Hadamard_Q0_map, PiEFm_Q1_map, idling_map_T3, PiEFm_Q0_map, idling_map_T4, PiEF_Q1_map, idling_map_T1, PiEF_Q0_map, idling_map_T2, Hadamard_Q0_dag_map])\n",
    "\n",
    "ideal_CSUMm1_map = reduce(np.dot, [Hadamard_Q0_map, PiEFm_Q1_map, idling_nodecay_map_T4, PiEFm_Q0_map, idling_nodecay_map_T3, PiEF_Q1_map, idling_nodecay_map_T2, PiEF_Q0_map, idling_nodecay_map_T1, Hadamard_Q0_dag_map])\n",
    "ideal_CMINm1_map = reduce(np.dot, [Hadamard_Q0_map, PiEFm_Q1_map, idling_nodecay_map_T3, PiEFm_Q0_map, idling_nodecay_map_T4, PiEF_Q1_map, idling_nodecay_map_T1, PiEF_Q0_map, idling_nodecay_map_T2, Hadamard_Q0_dag_map])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 401,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Scrambling_map(CSUM_map, CMIN_map, scrambling_param = 1):\n",
    "    \"\"\" \n",
    "    When scrambling_param == 1, returns the scrambler.\n",
    "    When scrambling_param == 0, returns the identity\n",
    "    \n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    TwoQHad_unitary = np.kron(Hadamard_to_identity(scrambling_param), Hadamard_to_identity(scrambling_param))\n",
    "    TwoQHadDag_unitary = np.conj(TwoQHad_unitary).T\n",
    "    \n",
    "    \n",
    "    TwoQHad_map = unitary_op_to_map(TwoQHad_unitary)\n",
    "    TwoQHadDag_map = unitary_op_to_map(TwoQHadDag_unitary)\n",
    "    \n",
    "    return reduce(np.dot, [TwoQHadDag_map, CMIN_map, TwoQHad_map, CSUM_map])\n",
    "\n",
    "def Scrambling_star_map(CSUM_map, CMIN_map, scrambling_param = 1):\n",
    "    \"\"\" When scrambling_param == 1, returns the scrambler conjugate.\n",
    "    When scrambling_param == 0, returns the identity conjugate (which is just the identity)\"\"\"\n",
    "    \n",
    "    TwoQHadStar_unitary = np.kron(HadamardStar_to_identity(scrambling_param), HadamardStar_to_identity(scrambling_param))\n",
    "    TwoQHadDagStar_unitary = np.conj(TwoQHadStar_unitary).T\n",
    "    \n",
    "    TwoQHadStar_map = unitary_op_to_map(TwoQHadStar_unitary)\n",
    "    TwoQHadDagStar_map = unitary_op_to_map(TwoQHadDagStar_unitary)\n",
    "    \n",
    "    return reduce(np.dot, [TwoQHadDagStar_map, CMIN_map, TwoQHadStar_map, CSUM_map])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 402,
   "metadata": {},
   "outputs": [],
   "source": [
    "ideal_scrambling_map = Scrambling_map(ideal_CSUMp1_map, ideal_CMINp1_map, 1)\n",
    "real_scrambling_map = Scrambling_map(real_CSUMp1_map, real_CMINp1_map, 1)\n",
    "\n",
    "ideal_scramblingstar_map = Scrambling_star_map(ideal_CSUMm1_map, ideal_CMINm1_map, 1)\n",
    "real_scramblingstar_map = Scrambling_star_map(real_CSUMm1_map, real_CMINm1_map, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 321,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Hadamard_to_identity(alpha):\n",
    "    \"\"\" returns a single-qutrit unitary matrix which\n",
    "    smoothly interpolates between the Hadamard gate and the identity gate\n",
    "    When alpha is unity, it is the Hadamard\n",
    "    When alpha is zero, it is the identity\n",
    "    \"\"\"\n",
    "    \n",
    "    Hadamard_pulse_sequence = [['EFY-90'],\n",
    "    ['X90'],\n",
    "    ['Z{:.5f}'.format(-109.471*alpha)],\n",
    "    ['EFZ{:.5f}'.format(54.7355*alpha)],\n",
    "    ['X-90'],\n",
    "    ['Z{:.5f}'.format(180*alpha)],\n",
    "    ['EFZ{:.5f}'.format(-90*alpha)],\n",
    "    ['EFY90']]\n",
    "    \n",
    "    return resulting_unitary(Hadamard_pulse_sequence)\n",
    "\n",
    "def HadamardStar_to_identity(alpha):\n",
    "    \"\"\" returns a single-qutrit unitary matrix which\n",
    "    smoothly interpolates between the Hadamard gate and the identity gate\n",
    "    When alpha is unity, it is the Hadamard\n",
    "    When alpha is zero, it is the identity\n",
    "    \"\"\"\n",
    "    \n",
    "    Hadamard_pulse_sequence = [['EFY-90'],\n",
    "    ['X-90'],\n",
    "    ['Z{:.5f}'.format(109.471*alpha)],\n",
    "    ['EFZ{:.5f}'.format(-54.7355*alpha)],\n",
    "    ['X90'],\n",
    "    ['Z{:.5f}'.format(-180*alpha)],\n",
    "    ['EFZ{:.5f}'.format(90*alpha)],\n",
    "    ['EFY90']]\n",
    "    \n",
    "    return resulting_unitary(Hadamard_pulse_sequence)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qnl_ctrl.qnl_analysis import SimTools as ST\n",
    "from qnl_ctrl.qnl_analysis.SimTools import bit_to_trit, rot_x, rot_y, trit_z\n",
    "\n",
    "def resulting_unitary(pulse_sequence):\n",
    "    # starting from the initial state, calculates the unitary after the pulse_sequence is applied                                                                            \n",
    "    # here time flows from left to right, in that the first pulse in the pulse_string is the first pulse applied.\n",
    "    single_qutrit_unitary = np.eye(3)\n",
    "    for x in pulse_sequence:\n",
    "        pulse_string = x[0]\n",
    "        if pulse_string in list(qutrit_mapping.keys()):\n",
    "            single_qutrit_unitary = np.dot(qutrit_mapping[pulse_string], single_qutrit_unitary)\n",
    "        else:\n",
    "            if pulse_string[0] == 'Z':\n",
    "                phase = np.pi*float(pulse_string[1:])/180.\n",
    "                single_qutrit_unitary = np.dot(trit_z(phase,-1)[0], single_qutrit_unitary)\n",
    "\n",
    "            elif pulse_string[0:3] == 'EFZ':\n",
    "                phase = np.pi*float(pulse_string[3:])/180.\n",
    "                single_qutrit_unitary = np.dot(trit_z(phase,1)[0], single_qutrit_unitary)\n",
    "    return single_qutrit_unitary\n",
    "\n",
    "qutrit_mapping = {'I': np.eye(3, dtype='complex'),                                                                                                                               \n",
    "                  'GPrep': np.eye(3, dtype='complex'),\n",
    "                  'X90': bit_to_trit(rot_x(np.pi/2.)),\n",
    "                  'X-90': bit_to_trit(rot_x(-np.pi/2.)),\n",
    "                  'Y90': bit_to_trit(rot_y(np.pi/2.)),\n",
    "                  'Y-90': bit_to_trit(rot_y(-np.pi/2.)),\n",
    "                  'X270': bit_to_trit(rot_x(-np.pi/2.)),\n",
    "                  'Y270': bit_to_trit(rot_y(-np.pi/2.)),\n",
    "                  'X180': bit_to_trit(rot_x(np.pi)),\n",
    "                  'Y180': bit_to_trit(rot_y(np.pi)),\n",
    "                  'EFX90': bit_to_trit(rot_x(np.pi/2.), 1),\n",
    "                  'EFX-90': bit_to_trit(rot_x(-np.pi/2.), 1),\n",
    "                  'EFY90': bit_to_trit(rot_y(np.pi/2.), 1),\n",
    "                  'EFY-90': bit_to_trit(rot_y(-np.pi/2.), 1),\n",
    "                  'EFX270': bit_to_trit(rot_x(-np.pi/2.), 1),\n",
    "                  'EFY270': bit_to_trit(rot_y(-np.pi/2.), 1),\n",
    "                  'EFX180': bit_to_trit(rot_x(np.pi), 1),\n",
    "                  'EFY180': bit_to_trit(rot_y(np.pi), 1),\n",
    "                  'Delay': np.eye(3, dtype='complex'),\n",
    "                  'Hadamard': np.sqrt(1./3)*np.array([[1,1,1],\n",
    "                                        [1,np.exp(1j*2*np.pi/3),np.exp(-1j*2*np.pi/3)],\n",
    "                                        [1,np.exp(-1j*2*np.pi/3),np.exp(1j*2*np.pi/3)]]),\n",
    "                  'R': np.diag([1, np.exp(1j*2*np.pi/3),np.exp(1j*2*np.pi/3)])                   }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulating our algorithm performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 367,
   "metadata": {},
   "outputs": [],
   "source": [
    "EPR_state = np.array([1,0,0,0,1,0,0,0,1])/np.sqrt(3)\n",
    "\n",
    "psi0 = reduce(np.kron, [state_to_teleport, EPR_state, EPR_state])\n",
    "rho0 = ST.vec_to_dm(psi0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 368,
   "metadata": {},
   "outputs": [],
   "source": [
    "EPR_projector = np.outer(EPR_state, np.conj(EPR_state))\n",
    "full_EPR_projector = reduce(np.kron, [np.eye(3), EPR_projector, np.eye(9)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to combine the maps to get the full five-qutrit quantum map.  Based on the equation below, \n",
    "\n",
    "\\begin{equation}\n",
    "\\hat { \\mathbf { A } } _ { i } \\hat { \\boldsymbol { \\rho } } \\hat { \\mathbf { A } } _ { i } ^ { \\dagger } \\rightarrow \\left( \\left( \\hat { \\mathbf { A } } _ { i } ^ { \\dagger } \\right) ^ { T } \\otimes \\hat { \\mathbf { A } } _ { i } \\right) \\vec { r }\n",
    "\\end{equation}\n",
    "\n",
    "we can infer that a general superoperator transforms as \n",
    "\\begin{equation}\n",
    "\\hat { \\mathbf { A } }  \\hat { \\boldsymbol { \\rho } } \\hat { \\mathbf { B } } \\rightarrow \\left( \\left( \\hat { \\mathbf { B } }  \\right) ^ { T } \\otimes \\hat { \\mathbf { A } }\\right) \\vec { r }\n",
    "\\end{equation}\n",
    "\n",
    "Unfortunately, as far as I can tell, there is no general way to combine two vectorized maps $\\mathcal{L}_1(\\rho_1)$ and $\\mathcal{L}_2(\\rho_2)$ to get vectorized form for the general map $\\mathcal{L}_{total}(\\rho_1 \\otimes \\rho_2) \\equiv \\mathcal{L}_1(\\rho_1) \\otimes \\mathcal{L}_2(\\rho_2)$.  But, it is pretty clear that if we have the Kraus representation for both the maps $\\mathcal{L}_1$ and $\\mathcal{L}_2$, then we can easily convert that to a vectorized superoperator.\n",
    "\n",
    "So, we have to now take a detour and figure out how to get the Kraus map from a vectorized superoperator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Getting a Kraus map from a vectorized superoperator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The answer was given in [this paper](https://aip.scitation.org/doi/abs/10.1063/1.1518555) in 2003.  Basically, one converts from a superoperator $\\mathcal{S}$ to the so-called Choi matrix $\\mathcal{T}$ via the transformation \n",
    "\n",
    "\\begin{equation}\n",
    "\\mathcal{T} = \\sum_{i,j=0}^{N-1} \\left(\\mathbb{E}_{ij}\\otimes\\mathbb{I}\\right)\\mathcal{S}\\left( \\mathbb{I}\\otimes\\mathbb{E}_{ij}\\right),\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "where $\\mathbb{I}$ is the $N$ by $N$ identity matrix (where $N$ is the Hilbert space dimension, e.g. $N=9$ for a two-qutrit system), and $\\mathbb{E}_{ij}$ is the $N$ by $N$ matrix with a single element equal to unity (the $i,j$ element).  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 505,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "def ABsuperop_to_matrix(A, B):\n",
      "    \"\"\" given operators A and B, in terms of which the superoperator Sis defined as\n",
      "    S(rho) = A rho B, this method returns the matrix representation of S\"\"\"\n",
      "\n",
      "    return np.kron(B.T, A)\n",
      "\n",
      "def Kraus_ops_to_matrix(left_ops, right_ops):\n",
      "    \"\"\"given a list of left_ops and right_ops such that\n",
      "    S(rho) = sum_i leftops[i] rho rightops[i], this method returns the matrix representation of S\"\"\"\n",
      "    \n",
      "    return np.sum([ABsuperop_to_matrix(left_ops[i], right_ops[i]) for i in range(len(left_ops))],axis = 0)\n",
      "\n",
      "def E(i,j,N):\n",
      "    E_matrix = np.zeros((N,N))\n",
      "    E_matrix[i,j] = 1\n",
      "    return E_matrix\n",
      "\n",
      "def choi_matrix(superoperator_matrix):\n",
      "    \"\"\" Returns the Choi matrix corresopnding to the given superoperator-matrix. \"\"\"\n",
      "    \n",
      "    # make sure the superoperator_matrix is square\n",
      "    assert superoperator_matrix.shape[0] == superoperator_matrix.shape[1]\n",
      "    \n",
      "    Hilbert_space_dimension = int(np.floor(np.sqrt(superoperator_matrix.shape[0])))\n",
      "    assert Hilbert_space_dimension * Hilbert_space_dimension == superoperator_matrix.shape[0]\n",
      "    \n",
      "    choi_matrix_placeholder = np.zeros(superoperator_matrix.shape)\n",
      "    \n",
      "    for i in np.arange(Hilbert_space_dimension):\n",
      "        for j in np.arange(Hilbert_space_dimension):\n",
      "            EkronI = np.kron(E(i,j,Hilbert_space_dimension), np.eye(Hilbert_space_dimension))\n",
      "            IkronE = np.kron(np.eye(Hilbert_space_dimension), E(i,j,Hilbert_space_dimension))\n",
      "            MappedIkronE = np.dot(superoperator_matrix, IkronE)\n",
      "            choi_matrix_placeholder = choi_matrix_placeholder + np.dot(EkronI, MappedIkronE)\n",
      "    \n",
      "    return choi_matrix_placeholder\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(inspect.getsource(ST.ABsuperop_to_matrix))\n",
    "print(inspect.getsource(ST.Kraus_ops_to_matrix))\n",
    "print(inspect.getsource(ST.E))\n",
    "print(inspect.getsource(ST.choi_matrix))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we have the choi matrix $\\mathcal{T}$, we can use the SVD of T, i.e. the decomposition of $\\mathcal{T}$ into\n",
    "\\begin{equation}\n",
    "\\mathcal{T} \\equiv \\mathcal{V}\\mathcal{\\Omega}\\mathcal{W}^\\dagger = \\sum_{k=0}^{N^2-1} \\omega_k \\mathbf{v_k}\\mathbf{w_k}^\\dagger,\n",
    "\\end{equation}\n",
    "we can write that \n",
    "\\begin{equation}\n",
    "\\mathcal{S}\\mathrm{col}\\left( \\mathbf{X} \\right)  = \\mathrm{col} \\left( \\sum_{k=0}^{N^2-1} \\omega_k \\mathbf{V_k}\\mathbf{X}\\mathbf{W_k}^\\dagger \\right),\n",
    "\\end{equation}\n",
    "where $\\mathrm{col}(\\mathbf{V_k}) = \\mathbf{v}_k$ and $\\mathrm{col}(\\mathbf{W_k}) = \\mathbf{w}_k$.\n",
    "\n",
    "This is not exactly the Kraus form, but for our purposes we don't care that much, since we know how to vectorize things in the form above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 506,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "def to_Kraus_ops(superoperator_matrix, svd_threshold = 1e-10):    \n",
      "    \"\"\"\n",
      "    Returns the Kraus operators corresponding to the superoperator_matrix\n",
      "    \n",
      "    Arguments:\n",
      "    superoperator_matrix: the matrix representation of the superoperator, in the vectorized basis\n",
      "    svd_threshold: the threshold value below which singular values are considered to be zero\n",
      "    \n",
      "    \"\"\"\n",
      "    Hilbert_space_dimension = int(np.floor(np.sqrt(superoperator_matrix.shape[0])))\n",
      "    assert Hilbert_space_dimension * Hilbert_space_dimension == superoperator_matrix.shape[0]\n",
      "\n",
      "    \n",
      "    # first convert the superoperator matrix into a Choi matrix\n",
      "    choi = choi_matrix(superoperator_matrix)\n",
      "    \n",
      "    u, s, v = np.linalg.svd(choi)\n",
      "    \n",
      "    #we don't really need u on their own, so we can just take the dot product of u and s\n",
      "    us = np.dot(u,np.diag(s))\n",
      "    \n",
      "    #to get the columns we should take the transpose of v\n",
      "    vT = v.T \n",
      "    \n",
      "    left_ops = []\n",
      "    right_ops = []\n",
      "    for i in range(superoperator_matrix.shape[0]):\n",
      "        if s[i] >= svd_threshold:\n",
      "            left_op = us[:,i].reshape(Hilbert_space_dimension, Hilbert_space_dimension).T\n",
      "            right_op = vT[:,i].reshape(Hilbert_space_dimension, Hilbert_space_dimension)\n",
      "            left_ops.append(left_op)\n",
      "            right_ops.append(right_op)\n",
      "        else:\n",
      "            break\n",
      "        \n",
      "    return left_ops, right_ops\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(inspect.getsource(ST.to_Kraus_ops))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 503,
   "metadata": {},
   "outputs": [],
   "source": [
    "states = {'GS': np.array([1,0,0]), 'ES': np.array([0,1,0]), 'FS': np.array([0,0,1])}\n",
    "\n",
    "scrambling_param = 1\n",
    "\n",
    "teleportation_fidelities_ideal = {}\n",
    "teleportation_fidelities_real = {}\n",
    "\n",
    "EPR_state = np.array([1,0,0,0,1,0,0,0,1])/np.sqrt(3)\n",
    "\n",
    "\n",
    "ideal_scrambling_map = Scrambling_map(ideal_CSUMp1_map, ideal_CMINp1_map, scrambling_param)\n",
    "real_scrambling_map = Scrambling_map(real_CSUMp1_map, real_CMINp1_map, scrambling_param)\n",
    "\n",
    "ideal_scramblingstar_map = Scrambling_star_map(ideal_CSUMm1_map, ideal_CMINm1_map, scrambling_param)\n",
    "real_scramblingstar_map = Scrambling_star_map(real_CSUMm1_map, real_CMINm1_map, scrambling_param)\n",
    "\n",
    "###\n",
    "# calculate real evolution\n",
    "###\n",
    "l1real, r1real = ST.to_Kraus_ops(real_scrambling_map)\n",
    "l2real, r2real = ST.to_Kraus_ops(real_scramblingstar_map)\n",
    "\n",
    "ltotreal, rtotreal = ST.Krauslike_tensor(l1real, r1real, l2real, r2real)\n",
    "ltotreal, rtotreal = ST.Krauslike_tensor(ltotreal, rtotreal, [np.eye(3)], [np.eye(3)])\n",
    "\n",
    "###\n",
    "# calculate ideal evolution\n",
    "###\n",
    "\n",
    "l1ideal, r1ideal = ST.to_Kraus_ops(ideal_scrambling_map)\n",
    "l2ideal, r2ideal = ST.to_Kraus_ops(ideal_scramblingstar_map)\n",
    "\n",
    "ltotideal, rtotideal = ST.Krauslike_tensor(l1ideal, r1ideal, l2ideal, r2ideal)\n",
    "ltotideal, rtotideal = ST.Krauslike_tensor(ltotideal, rtotideal, [np.eye(3)], [np.eye(3)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 504,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.690707263744\n",
      "0.6662395752250001\n",
      "0.696181640625\n",
      "0.7109272172250001\n",
      "0.7552444263040001\n",
      "0.7174784556810001\n",
      "0.6904446648999999\n",
      "0.677635190596\n",
      "0.7319555359359999\n",
      "0.7092166225\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-504-070febcf8fe5>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     15\u001b[0m     \u001b[0mrho_app\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrho0\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     16\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mltotreal\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 17\u001b[0;31m         \u001b[0mrho_app\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mrho_app\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mreduce\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mltotreal\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrho0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrtotreal\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     18\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     19\u001b[0m     \u001b[0mrho_proj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mreduce\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mfull_EPR_projector\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrho_app\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfull_EPR_projector\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "n_randomizations = 100\n",
    "\n",
    "fidelities = []\n",
    "\n",
    "for _ in range(n_randomizations):\n",
    "\n",
    "    input_state = np.random.randn(3) + 1j*np.random.randn(3)\n",
    "    input_state = input_state / np.linalg.norm(input_state)\n",
    "\n",
    "    input_dm = ST.vec_to_dm(input_state)\n",
    "\n",
    "    psi0 = reduce(np.kron, [input_state, EPR_state, EPR_state])\n",
    "    rho0 = ST.vec_to_dm(psi0)\n",
    "\n",
    "    rho_app = np.zeros(rho0.shape)\n",
    "    for i in range(len(ltotreal)):\n",
    "        rho_app = rho_app + reduce(np.dot, [ltotreal[i], rho0, rtotreal[i]])\n",
    "\n",
    "    rho_proj = reduce(np.dot, [full_EPR_projector, rho_app, full_EPR_projector])\n",
    "    rho_proj = rho_proj / np.trace(rho_proj)\n",
    "\n",
    "    final_state = ST.partial_trace(rho_proj, k = 0, dim = [[81,81],[3,3]])\n",
    "    final_state = final_state / np.trace(final_state)\n",
    "\n",
    "    fid = ST.fidelity(final_state, input_dm)\n",
    "    fidelities.append(fid)\n",
    "    print(fid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}