control-pulseoptim-CRAB-2qubitInerac.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculation of control fields for state-to-state transfer of a 2 qubit system using CRAB algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Jonathan Zoller (jonathan.zoller@uni-ulm.de)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example to demonstrate using the control library to determine control\n",
    "pulses using the ctrlpulseoptim.optimize_pulse_unitary function.\n",
    "The CRAB algorithm is used to optimize pulse shapes to minimize the fidelity\n",
    "error, which is equivalent maximising the fidelity to an optimal value of 1.\n",
    "\n",
    "The system in this example are two qubits, where the interaction can be\n",
    "controlled. The target is to perform a pure state transfer from a down-down\n",
    "state to an up-up state.\n",
    "\n",
    "The user can experiment with the timeslicing, by means of changing the\n",
    "number of timeslots and/or total time for the evolution.\n",
    "Different initial (starting) pulse types can be tried as well as\n",
    "boundaries on the control and a smooth ramping of the pulse when\n",
    "switching the control on and off (at the beginning and close to the end).\n",
    "The initial and final pulses are displayed in a plot\n",
    "\n",
    "An in depth discussion of using methods of this type can be found in [1,2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import datetime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qutip import Qobj, identity, sigmax, sigmaz, tensor, mesolve\n",
    "import random\n",
    "import qutip.logging_utils as logging\n",
    "logger = logging.get_logger()\n",
    "#Set this to None or logging.WARN for 'quiet' execution\n",
    "log_level = logging.INFO\n",
    "#QuTiP control modules\n",
    "import qutip.control.pulseoptim as cpo\n",
    "\n",
    "example_name = '2qubitInteract'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining the physics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dynamics of the system are governed by the combined Hamiltonian:\n",
    "H(t) = H_d + sum(u1(t)*Hc1 + u2(t)*Hc2 + ....)\n",
    "That is the time-dependent Hamiltonian has a constant part (called here the drift) and time vary parts, which are the control Hamiltonians scaled by some functions u_j(t) known as control amplitudes\n",
    "In this example we describe an Ising like Hamiltonian, encompassing random coefficients in the drift part and controlling the interaction of the qubits:\n",
    "\n",
    "$ \\hat{H} = \\sum_{i=1}^2 \\alpha_i \\sigma_x^i + \\beta_i \\sigma_z^i + u(t) \\cdot \\sigma_z \\otimes \\sigma_z $\n",
    "\n",
    "Initial $\\newcommand{\\ket}[1]{\\left|{#1}\\right\\rangle} \\ket{\\psi_0} = \\text{U_0}$ and target state $\\ket{\\psi_t} = \\text{U_targ}$ are chosen to be:\n",
    "\n",
    "$ \\ket{\\psi_0} = \\begin{pmatrix} 1 \\\\ 0 \\\\ 0 \\\\ 0 \\end{pmatrix}$\n",
    "\n",
    "$ \\ket{\\psi_t} = \\begin{pmatrix} 0 \\\\ 0 \\\\ 0 \\\\ 1 \\end{pmatrix}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "random.seed(20)\n",
    "alpha = [random.random(),random.random()]\n",
    "beta  = [random.random(),random.random()]\n",
    "\n",
    "Sx = sigmax()\n",
    "Sz = sigmaz()\n",
    "\n",
    "H_d = (alpha[0]*tensor(Sx,identity(2)) + \n",
    "      alpha[1]*tensor(identity(2),Sx) +\n",
    "      beta[0]*tensor(Sz,identity(2)) +\n",
    "      beta[1]*tensor(identity(2),Sz))\n",
    "H_c = [tensor(Sz,Sz)]\n",
    "# Number of ctrls\n",
    "n_ctrls = len(H_c)\n",
    "\n",
    "q1_0 = q2_0 = Qobj([[1], [0]])\n",
    "q1_targ = q2_targ = Qobj([[0], [1]])\n",
    "\n",
    "psi_0 = tensor(q1_0, q2_0)\n",
    "psi_targ = tensor(q1_targ, q2_targ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining the time evolution parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To solve the evolution the control amplitudes are considered constant within piecewise timeslots, hence the evolution during the timeslot can be calculated using U(t_k) = expm(-i*H(t_k)*dt). Combining these for all the timeslots gives the approximation to the evolution from an initial state $\\psi_0$ at t=0 to U(T) at the t=evo_time.\n",
    "The number of timeslots and evo_time have to be chosen such that the timeslot durations (dt) are small compared with the dynamics of the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of time slots\n",
    "n_ts = 100\n",
    "# Time allowed for the evolution\n",
    "evo_time = 18"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set the conditions which will cause the pulse optimisation to terminate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At each iteration the fidelity of the evolution is tested by comparaing the calculated evolution U(T) with the target U_targ. For unitary systems such as this one this is typically:\n",
    "f = normalise(overlap(U(T), U_targ)). The maximum fidelity (for a unitary system) calculated this way would be 1, and hence the error is calculated as fid_err = 1 - fidelity. As such the optimisation is considered completed when the fid_err falls below such a target value.\n",
    "\n",
    "In some cases the optimisation either gets stuck in some local minima, or the fid_err_targ is just not achievable, therefore some limits are set to the time/effort allowed to find a solution.\n",
    "\n",
    "The algorithm uses the CRAB algorithm to determine optimized coefficients that lead to a minimal fidelity error. The underlying optimization procedure is set to be the Nelder-Mead downhill simplex. Therefore, when all vertices shrink together, the algorithm will terminate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fidelity error target\n",
    "fid_err_targ = 1e-3\n",
    "# Maximum iterations for the optisation algorithm\n",
    "max_iter = 600\n",
    "# Maximum (elapsed) time allowed in seconds\n",
    "max_wall_time = 120"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set the initial pulse type"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The control amplitudes must be set to some initial values. Typically these are just random values for each control in each timeslot. These do however result in erratic optimised pulses. For this example, a solution will be found for any initial pulse, and so it can be interesting to look at the other initial pulse alternatives."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pulse type alternatives: RND|ZERO|LIN|SINE|SQUARE|SAW|TRIANGLE|\n",
    "p_type = 'RND'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Give an extension for output files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Set to None to suppress output files\n",
    "f_ext = \"{}_n_ts{}_ptype{}.txt\".format(example_name, n_ts, p_type)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the optimisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this step, the actual optimization is performed. At each iteration the Nelder-Mead algorithm calculates a new set of coefficients that improves the currently worst set among all set of coefficients. For details see [1,2] and a textbook about static search methods. The algorithm continues until one of the termination conditions defined above has been reached. If undesired results are achieved, rerun the algorithm and/or try to change the number of coefficients to be optimized for, as this is a very crucial parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:qutip.control.dynamics:Setting memory optimisations for level 0\n",
      "INFO:qutip.control.dynamics:Internal operator data type choosen to be <class 'numpy.ndarray'>\n",
      "INFO:qutip.control.dynamics:phased dynamics generator caching True\n",
      "INFO:qutip.control.dynamics:propagator gradient caching True\n",
      "INFO:qutip.control.dynamics:eigenvector adjoint caching True\n",
      "INFO:qutip.control.dynamics:use sparse eigen decomp False\n",
      "INFO:qutip.control.pulseoptim:System configuration:\n",
      "Drift Hamiltonian:\n",
      "Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True\n",
      "Qobj data =\n",
      "[[ 1.67112549  0.68625416  0.90563968  0.        ]\n",
      " [ 0.68625416 -0.13810698  0.          0.90563968]\n",
      " [ 0.90563968  0.          0.13810698  0.68625416]\n",
      " [ 0.          0.90563968  0.68625416 -1.67112549]]\n",
      "Control 1 Hamiltonian:\n",
      "Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True\n",
      "Qobj data =\n",
      "[[ 1.  0.  0.  0.]\n",
      " [ 0. -1.  0.  0.]\n",
      " [ 0.  0. -1.  0.]\n",
      " [ 0.  0.  0.  1.]]\n",
      "Initial state / operator:\n",
      "Quantum object: dims = [[2, 2], [1, 1]], shape = (4, 1), type = ket\n",
      "Qobj data =\n",
      "[[ 1.]\n",
      " [ 0.]\n",
      " [ 0.]\n",
      " [ 0.]]\n",
      "Target state / operator:\n",
      "Quantum object: dims = [[2, 2], [1, 1]], shape = (4, 1), type = ket\n",
      "Qobj data =\n",
      "[[ 0.]\n",
      " [ 0.]\n",
      " [ 0.]\n",
      " [ 1.]]\n",
      "INFO:qutip.control.pulseoptim:Initial amplitudes output to file: ctrl_amps_initial_2qubitInteract_n_ts100_ptypeRND.txt\n",
      "INFO:qutip.control.optimizer:Optimising pulse(s) using CRAB with 'fmin' (Nelder-Mead) method\n",
      "INFO:qutip.control.pulseoptim:Final amplitudes output to file: ctrl_amps_final_2qubitInteract_n_ts100_ptypeRND.txt\n"
     ]
    }
   ],
   "source": [
    "result = cpo.opt_pulse_crab_unitary(H_d, H_c, psi_0, psi_targ, n_ts, evo_time, \n",
    "                fid_err_targ=fid_err_targ, \n",
    "                max_iter=max_iter, max_wall_time=max_wall_time, \n",
    "                init_coeff_scaling=5.0, num_coeffs=5, \n",
    "                method_params={'xtol':1e-3},\n",
    "                guess_pulse_type=p_type, guess_pulse_action='modulate',\n",
    "                out_file_ext=f_ext,\n",
    "                log_level=log_level, gen_stats=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Report the results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Firstly the performace statistics are reported, which gives a breakdown of the processing times. In this example it can be seen that the majority of time is spent calculating the propagators, i.e. exponentiating the combined Hamiltonian.\n",
    "\n",
    "The optimised U(T) is reported as the 'final evolution', which is essentially the string representation of the Qobj that holds the full time evolution at the point when the optimisation is terminated.\n",
    "\n",
    "The key information is in the summary (given last). Here the final fidelity is reported and the reason for termination of the algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "------------------------------------\n",
      "---- Control optimisation stats ----\n",
      "**** Timings (HH:MM:SS.US) ****\n",
      "Total wall time elapsed during optimisation: 0:00:20.479439\n",
      "Wall time computing Hamiltonians: 0:00:00.373044 (1.82%)\n",
      "Wall time computing propagators: 0:00:19.826778 (96.81%)\n",
      "Wall time computing forward propagation: 0:00:00.073227 (0.36%)\n",
      "Wall time computing onward propagation: 0:00:00.065351 (0.32%)\n",
      "Wall time computing gradient: 0:00:00 (0.00%)\n",
      "\n",
      "**** Iterations and function calls ****\n",
      "Number of iterations: 500\n",
      "Number of fidelity function calls: 706\n",
      "Number of times fidelity is computed: 706\n",
      "Number of gradient function calls: 0\n",
      "Number of times gradients are computed: 0\n",
      "Number of times timeslot evolution is recomputed: 706\n",
      "\n",
      "**** Control amplitudes ****\n",
      "Number of control amplitude updates: 705\n",
      "Mean number of updates per iteration: 1.41\n",
      "Number of timeslot values changed: 70499\n",
      "Mean number of timeslot changes per update: 99.99858156028368\n",
      "Number of amplitude values changed: 70499\n",
      "Mean number of amplitude changes per update: 99.99858156028368\n",
      "------------------------------------\n",
      "Final evolution\n",
      "Quantum object: dims = [[2, 2], [1, 1]], shape = (4, 1), type = ket\n",
      "Qobj data =\n",
      "[[-0.00479632 -1.54708545e-04j]\n",
      " [ 0.02948317 +1.72419994e-02j]\n",
      " [-0.00998838 +2.15479124e-02j]\n",
      " [-0.82590384 -5.62253680e-01j]]\n",
      "\n",
      "********* Summary *****************\n",
      "Final fidelity error 0.0008772110044362735\n",
      "Final gradient normal 0.0\n",
      "Terminated due to Goal achieved\n",
      "Number of iterations 500\n",
      "Completed in 0:00:20.479439 HH:MM:SS.US\n"
     ]
    }
   ],
   "source": [
    "result.stats.report()\n",
    "print(\"Final evolution\\n{}\\n\".format(result.evo_full_final))\n",
    "print(\"********* Summary *****************\")\n",
    "print(\"Final fidelity error {}\".format(result.fid_err))\n",
    "print(\"Final gradient normal {}\".format(result.grad_norm_final))\n",
    "print(\"Terminated due to {}\".format(result.termination_reason))\n",
    "print(\"Number of iterations {}\".format(result.num_iter))\n",
    "print(\"Completed in {} HH:MM:SS.US\".format(\n",
    "        datetime.timedelta(seconds=result.wall_time)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the initial and final amplitudes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here the (random) starting pulse is plotted along with the pulse (control amplitudes) that was found to produce the target gate evolution to within the specified error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3XncHFWd7/HPN4hGZI1BJEBIQHQmzghyAwFUBgZlAIEoOoAo4HbBGWXG9cqIC4teBxVHREfAFVAwgqKArCrIXIVAwoSwKpBFQgIIkQQGRAK/+0edTiqd7n7qWaq7qvv7fr2e19Nddar619XV/atz6tQpRQRmZmZVM67XAZiZmbXiBGVmZpXkBGVmZpXkBGVmZpXkBGVmZpXkBGVmZpXkBGW1J+kKSUd3mH+mpE8VXNd1kt47dtF1j6S9JC3pdRxmY8UJyipJ0iJJry9SNiL2j4hz0nLvlPT/mua/LyJOGaO4Xi7pQkmPSFohab6kD0tab5TrdXIxa+IEZVaQpO2B2cD9wN9GxCbAPwLTgY268PrPK/s1zKrECcoqr1ErkvQlSX+StFDS/rn510l6r6S/Bs4Edpf0hKTH0vzvSfpseryZpMsk/TGt6zJJWxcM5STgtxHx4YhYBhARv4uIIyKi8VoHS7pD0mMprr/OxblI0kdTrWuFpFmSxkt6EXAFMCnF/YSkSZJOlHSRpO9LWgm8U9ILJH1F0tL09xVJLyi4HU+XdL+klZLmSnpdbt6JqWb4fUmPS7ot1Rb/TdLDabl9m7b55yXdlN7LzyRNSPPGp/U8mrbDzZK2KLiNzVZzgrK6mAH8DpgIfAH4tiTlC0TEXcD7gBsiYsOI2LTFesYB3wW2BSYDTwFfKxjD64GL2s2U9HLgAuCDwObA5cClkp6fK3YosB8wFXgV8M6I+B9gf2BpinvDiFiays9Mr7kp8APgBGA3YCdgR2BX4JMF4785LTcBOB+4UNL43PyDgPOAzYD/Bq4i215bAScDZzWt7yjg3cAkYBXw1TT9aGATYBvgxWSfyVMFYzRbzQnK6mJxRHwzIp4FzgG2BIZ9VB4Rj0bEjyPiyYh4HPgc8HcFF38xsKzD/MOAn0fENRHxDPAl4IXAHrkyX42IpRGxHLiULGF0ckNE/DQinouIp4C3AydHxMMR8UeyWt2RRYKPiO+n978qIk4DXgC8IlfkvyLiqohYBVxIlmT/Pb2XHwJTJOWT/nkRcXtKsJ8CDk3n4p4h21Yvi4hnI2JuRKwsEqNZnhOU1cWDjQcR8WR6uOFwVyJpA0lnSVqcms2uBzYt2MnhUbLE2M4kYHEuzufIzldtlSvzYO7xkwz9Hu7v9Brp8aQh1gGApI9Iuis1yT1GVsuZmCvyUO7xU8Aj6YCg8ZymePOxLQbWT+s7j6z29cPUDPkFSesXidEszwnK+s1Qw/N/hKzWMCMiNgb2TNPVfpHVfgG8pcP8pWRNh9kKsybIbYAHCqy7XdzN09d6DbJmyqUMIZ1v+jhZE+NmqflzBcXedzvbNMXxDFlSeyYiToqIaWS1xwPJmgPNhsUJyvrNQ8DWTed98jYiqw08lk7qf2YY6/4MsIekL0p6KYCkl6UOAZsCPwLeKGmfVGP4CPA08NuCcb9Y0iZDlLsA+KSkzSVNBD4NfL/A+jciO0/0R+B5kj4NbFxguU7eIWmapA3IzlFdFBHPStpb0t+mWulKssT1bMc1mbXgBGX95lfAHcCDkh5pMf8rZOeFHgFuBK4suuKIuA/YHZgC3CFpBfBjYA7weET8DngHcEZa/0HAQRHxlwLrvpss+SxIPd/aNdt9Nr3efOA24JY0bShXkfUU/D1Zc9yfWbf5cLjOA75H1mw5HviXNP2lZB07VgJ3Ab+mWBI1W4t8w0IzGy5J1wHfj4hv9ToW61+uQZmZWSU5QZmZWSW5ic/MzCrJNSgzM6ukSg0+OXHixJgyZUqvwzAzsxLNnTv3kYjYfKhyhRKUpBcCk1M32tJMmTKFOXPmlPkSZmbWY5IWD12qQBOfpIOAeaTrRSTtJOmS0YVnZmbWWZEa1IlkIyZfBxAR8yRNKS2iPnT+7D/ws3lrRruZudNWHDFjcg8jMjOrviIJalVErGi6s4F10JyQZi9cDsCMqRO4c1k2qLMTlJlZZ0US1O2SjgDWk7QD2XAmRcYWG1g/m/cAdy5bybQts6HOZkydsLrWdNhZN3DnspUcdtYNq8u7RmVmtq4iCeo4spukPU02VthVwCllBlVH+VpTIznNOnb3dcrN3GmrtZ67RmVm1tqQCSrde+eE9Gdt5GtN07bceJ1E1HDEjMlrJaN8TcrMzNZom6AkXUqHe+tExMGlRFRj7WpNZmY2fJ1qUF9K/w8hGz6/MVz+24BFJcY0cPLnpHw+ysws0zZBRcSvASSdEhF75mZdKun60iOruOaeevlOEcORbwr0+SgzszWKdJLYXNJ2EbEAQNJUYMghKvpdc0+9TuedOsmfk/L5KDOzNYokqA8B10lakJ5PAY4tLaIa8TknM7PyFOnFd2W6/umv0qS7I+LpcsMyM7NBN2SCknRU06QdJRER55YUk5mZWaEmvl1yj8cD+wC3AE5QJfAoE2ZmmSJNfMfln0vaBDivtIgqrNVoEWPJo0yYma0xkhsWPgnsMNaB1EHR0SJGyqNMmJmtUeQcVH5EiXHANODCMoOqMvfcMzPrjiI1qC/lHq8CFkfEkiIrl7QIeBx4luy2HdOHHaGZmQ2kIgnqgIj4eH6CpFObp3Wwd0Q8MvzQzMxskA15y3fgDS2m7T/WgZiZmeV1Gs38n4B/BraTND83ayPgNwXXH8DVkgI4KyLObvE6xwDHAEyeXL3eamX33BuKB5I1s0HVqYnvfOAK4PPA8bnpj0fE8oLrf01ELJX0EuAaSXdHxFoDzaakdTbA9OnT297eo1fK7rnXiQeSNbNB1ilBRUQskvT+5hmSJhRJUhGxNP1/WNLFwK5A7UZC71XPPQ8ka2aDbKga1IHAXLKmOuXmBbBdpxVLehEwLiIeT4/3BU4eXbhmZjYoOt0P6sD0f+oI170FcLGkxuucHxFXjnBdZmY2YDp1kti504IRccsQ8xcAO44wLjMzG3CdmvhO6zAvgL8f41gqYazulGtmY6v5uwnu2drvOjXx7d3NQKpirO6Ua2ajl09Ksxdm/bJmTJ2w+vnshctXz3ey6j9FxuIbT3Y91GvJak7/BZwZEX8uObae8Xh7Zr3RXEvKJ6UZUyeslYSak1c+WYETVj8oMtTRuWTj6Z2Rnr+N7HYb/1hWUNaaL9q1ftSpltSclPLyl2G0appvlLH6KpKgXhER+c4O10q6tayArDVftGv9Kt+s3ikhdeJb1fSnIgnqvyXtFhE3AkiaQfGhjmyM+KJd6xftOiKNdbO6Wxzqr0iCmgEcJekP6flk4C5Jt5GNNvGq0qLrkl6Pt2c2SLrREcktDv2hSILar/QoeqyX4+2ZDYJWB4FldkRqbnFwbaqehkxQEbFY0mbANvnyQ12oWzfuuWdWHg+6bCNRpJv5KcA7gftYc+v3vr1Q18zGRrdrTe34/G19FWniOxTYPiL+UnYwZlZfna5hctO5jUSRBHU7sCnwcMmxdE0/DGeUb1MHt6tb7zV3fhhpl/Gy+btTH0US1OfJuprfDjzdmBgRB5cWVcnqPpxRc6xuV7deqUozXlH+7tRLkQR1DnAqcBvwXLnhdE/Vv0id+KJEq4q69YD1d6deiiSoRyLiq6VHYmaV162LbLvJXdCrq0iCmivp88AlrN3EV6tu5r4Y12xkOo2VV4daUyfugl5tRRLUq9P/3XLTatfNvG5NEWa9MpwRxevOXdCrrciFun1zX6i6N0WYdUNdeuOVwT38qqVIDQpJbwReCYxvTIuIk8sKaiz0Q1dys26pW2+8MriHX/UUGUniTGADYG/gW8BbgZtKjmvU6t6V3KxMvqh2Xe7hVz1FalB7RMSrJM2PiJMknQb8pOzARuKkS+/gzqXZUc+gHQW6J5INxyA34w2Hv1e9VSRBPZX+PylpEvAoMLW8kMbGIB0FuieSFeFmvOEp+3vVXItt9fqD/j1WRHQuIH2K7Hbv+wBfJ+vB982I+PRYBzN9+vSYM2fOWK92oDSO9vzDY52a8cA/gMPRuGVH/jx2ke3XKQk1fx6d5vXbZyVpbkRMH6pckV58p6SHP5Z0GTA+IlaMNkAzK5eb8cZOc2vM7IXLmb1weccaUKMctE5CnT6P5touDGaryJA1qG5yDWr0XIMabG7G646hmufyRntQMNLaW5WNWQ3KzOrDF6R3R3OPvzINcvd3JyizmnOtqb8Ncvf3tglK0s6dFqzbWHxm/cLXMNmgdH/vVIM6rcO82o3FN0g8XEt/c+eHwTZIl5W0TVBjMQafpP2A04H1gG9FxL+Pdp3W2SC3V/czN+NZwyANcFtkqKP1gX8C9kyTrgPOiohnhlhuPbLrpt4ALAFulnRJRNw5qoito0Furx5rw+mp1clIajduxjMr1kniG8D6wH+m50emae8dYrldgXsjYgGApB8CMwEnKKuF5qa0kSh6vUyr5WDN9TNuxrN2+rlJv0iC2iUidsw9/5WkWwsstxVwf+75EmBGcyFJxwDHAEye3B8b1eprrJvSRloLc0KyIvq9Sb9IgnpW0vYRcR+ApO2AZwsspxbT1rkqOCLOBs6G7ELdAus1GzNlN6V183oZGzz93qRfJEF9DLhW0gKypLMt8K4Cyy0Btsk93xpYOuwIzUrkHnFm1dUxQUkaRzaa+Q7AK8gS1N0R8XSBdd8M7CBpKvAAcDhwxOjCNRs994gzq4eOCSoinpN0WkTsDswfzoojYpWkDwBXkXUz/05E3DHyUG2kBuWivqI8HJD1s376vhdp4rta0luAn8QwR5aNiMuBy0cUmY2JQbqobzhca7J+1G/f9yIJ6sPAi4BVkv5M1swXETHyvrfWNYN0UV87zR0hRtt13Kyq+u37Pm6oAhGxUUSMi4jnR8TG6bm/3VYbjSa9BjfrmdVDkZEkfhkR+ww1zaxK3BHCrP46jWY+HtgAmChpM9Zc17QxMKkLsZmNmDtCmNV/lIlONahjgQ+SJaO5rElQK8nG2LMaqvsO206780yuNdmg6odRJjqNZn46cLqk4yLijC7GZCXphx22neYLbl1rskHXD6NMDHkOKiLOkLQHMCVfPiLOLTEuK0E/7LCduMZk1l+KdJI4D9gemMeaMfgCcIKynmrVEcLM+keR66CmA9OGe5GuWdncEcJseOo2ykSRBHU78FJgWcmxWA/UbYd193GzkanjKBNFEtRE4E5JNwGrB4mNiINLi8q6oo47rGtNZiNTx1EmiiSoE8sOwnqjDjusu4+bDa4ivfh+LWkLYJc06aaIeLjcsKwXqnKNVD4pNd/63LUms7FRh+b9Ir34DgW+CFxHdrHuGZI+FhEXlRybdVGVrpHKN+P5BoJmY68uzftFmvhOAHZp1JokbQ78AnCC6iO9vEbKzXhm3VWH5n0olqDGNTXpPUqBUdCt/spsAnAznll1VKV5v1mRBHWlpKuAC9Lzw4ArygvJqiCfIGYvXM7shcvXquUU2YGba0Z5+aTkZjyz3mk+GGz1fc+bNmljPnPQK7sRGipy/a2kQ4DXkp2Duj4iLi4jmOnTp8ecOXPKWLWNQnOiaa7xtDNUOScls+rpdGAJY5OgJM2NiOlDlmuXoCS9DNgiIn7TNH1P4IGIuG9UEbbgBFUPQ+3AeU5CZtasaILq1MT3FeATLaY/meYdNMLYrOaaO1SYmZWhU2eHKRExv3liRMwhG9nczMysNJ0S1PgO81441oGYmZnldToHdQHwq4j4ZtP09wD7RsRhYx6M9Edg8ShXMxF4ZAzC6RXH31uOv3fqHDs4/uHYNiI2H6pQpwS1BXAx8BeyW75DduuN5wNvjogHxyjQMSVpTpGTb1Xl+HvL8fdOnWMHx1+GTrd8fwjYQ9LewN+kyT+PiF91JTIzMxtoRQaLvRa4tguxmJmZrdaPQxad3esARsnx95bj7506xw6Of8wVGknCzMys2/qxBmVmZn3ACcrMzCqptglK0n6SfifpXknHt5j/Akmz0vzZkqZ0P8p1SdpG0rWS7pJ0h6R/bVFmL0krJM1Lf5/uRaydSFok6bYU3zoDKCrz1bT950vauRdxNpP0itx2nSdppaQPNpWp3PaX9B1JD0u6PTdtgqRrJN2T/m/WZtmjU5l7JB3dvahXv36r2L8o6e60b1wsadM2y3bcz7qhTfwnSnogt48c0GbZjr9T3dAm/lm52BdJmtdm2d5u/4io3R+wHnAfsB3ZdVm3AtOayvwzcGZ6fDgwq9dxp1i2BHZOjzcCft8i9r2Ay3od6xDvYxEwscP8A8huyyJgN2B2r2Nusx89SHbRYKW3P7AnsDNwe27aF4Dj0+PjgVNbLDcBWJD+b5Yeb1aB2PcFnpcen9oq9iL7WQ/jPxH4aIH9q+PvVK/ib5p/GvDpKm7/utagdgXujYgFEfEX4IfAzKYyM4Fz0uOLgH0kqYsxthQRyyLilvT4ceAuoB/vzjcTODcyNwKbStqy10E12Qe4LyJGO3pJ6SLiemB50+T8Pn4O8KYWi/4DcE1ELI+IPwHXAPuVFmgLrWKPiKsjYlV6eiOwdTdjGo42276IIr9TpesUf/pNPJQ19/urlLomqK2A+3PPl7Duj/zqMumLsAJ4cVeiKyg1O74amN1i9u6SbpV0haTu3B1seAK4WtJcSce0mF/kM+q1w2n/xaz69ofsdjjLIDvwAV7SokwdPod30/4mqEPtZ730gdRE+Z02zat12PavAx6KiHvazO/p9q9rgmpVE2ruL1+kTM9I2hD4MfDBiFjZNPsWsmanHYEzgJ92O74CXhMROwP7A+9P9wnLq/r2fz5wMHBhi9l12P5FVf1zOAFYBfygTZGh9rNe+QawPbATsIysmaxZpbd98jY61556uv3rmqCWANvknm8NLG1XRtLzgE0YWTV9zElanyw5/SAiftI8PyJWRsQT6fHlwPqSJnY5zI4iYmn6/zDZmI27NhUp8hn10v7ALZEN6bWWOmz/5KFGs2n6/3CLMpX9HFKHjQOBt0c64dGswH7WExHxUEQ8GxHPAd+kdVyV3faw+nfxEGBWuzK93v51TVA3AztImpqOhA8HLmkqcwnQ6LH0VrKR2Xt+9JLafL8N3BURX25T5qWN82WSdiX7nB7tXpSdSXqRpI0aj8lOeN/eVOwS4KjUm283YEWjOaoi2h45Vn375+T38aOBn7UocxWwr6TNUjPUvmlaT0naD/g4cHBEPNmmTJH9rCeazqe+mdZxFfmd6qXXA3dHxJJWMyux/XvVO2O0f2S9xH5P1kvmhDTtZLIdHrL7WV0I3AvcBGzX65hTXK8lq+bPB+alvwOA9wHvS2U+ANxB1uvnRmCPXsfd9B62S7HdmuJsbP/8exDw9fT53AZM73Xcufg3IEs4m+SmVXr7kyXTZcAzZEfm7yE7p/pL4J70f0IqOx34Vm7Zd6fvwb3AuyoS+71k52ca34FGj9tJwOWd9rOKxH9e2q/nkyWdLZvjT8/X+Z2qQvxp+vca+3yubKW2v4c6MjOzSqprE5+ZmfU5JygzM6skJygzM6skJygzM6skJygzM6skJygzM6skJygzM6skJygzM6skJygzM6skJygzM6skJygzM6skJygzM6skJyjrKkmTJT0hab0RLv+EpO3GOKbvSfrsWK6zWyRNkRTp3j69jmUvSUtyz++QtNcYrn+RpNeP1fqs+pygrCNJ75R0m6QnJT0o6RuSNh3G8mv9qETEHyJiw4h4diTxpGUXjGTZkZK0paRvS1om6XFJd0s6Kd0jZzTr7VpykXSdpD9JekHZr9UQEa+MiOvS658o6fvdem3rD05Q1pakjwCnAh8juyPxbsC2wDXpBmx9T9IE4AbghcDuEbER8AZgU7Jbfpf9+qNOXpKmAK8juw/ZwaNdn1m3OEFZS5I2Bk4CjouIKyPimYhYBBxKlqTekcqdKOkiSbNS7eIWSTumeecBk4FLU9Pc/2muNaQj+89K+m0qc6mkF0v6gaSVkm5OP7CNuELSy9LjAyTdmV73AUkfzZU7UNI8SY+ldb8qN+/VKc7HJc0iu7llOx8GHgfekd4/EXF/RPxrRMxP69sjxbki/d8j91rXSTpF0m/S612tNbePvz79fyy9991TjfU3kv5D0nLgREnjJH1S0mJJD0s6V9Imw/g4jyK78eL3WHMH3kZ835P0n5KuSDH8Rtkdhb+Salx3S3p1rvwiSf+WtvufJH1XUsvt16g9K7t77ieAw9Jr3Jqfnyu/Vi1L0pHpPT8q6YSmdY+TdLyk+9L8H6WDCSSNl/T9NP2x9JlsMYztZRXhBGXt7EH2w/2T/MSIeAK4gqwW0TCT7O7FE4DzgZ9KWj8ijgT+AByUmua+0Oa1DgeOBLYiq5XcAHw3re8u4DNtlvs2cGyq1fwN8CsASTsD3wGOJbvr7FnAJZJekGp+PyW7I+qEFPdbOmyH1wM/iYjnWs1MP4o/B76aXuvLwM8lvThX7AjgXcBLgOcDjUS6Z/q/ado+N6TnM4AFqfzngHemv73J7nK6IfC1DjE3Owr4Qfr7hxY/1ocCnwQmAk+Tbf9b0vOL0nvKezvwD2Sf1cvTsm1FxJXA/wVmpfe541ABS5oGfINsv5hEtm23zhX5F+BNwN+l+X8iu4MzZEl4E2CbtNz7gKeGek2rHicoa2ci8EhErGoxb1ma3zA3Ii6KiGfIfszGkzUHFvXdiLgvIlaQJb/7IuIX6bUvBF7dZrlngGmSNo6IP0XELWn6/wbOiojZEfFsRJxD9sO7W/pbH/hKqhVeBNzcIbYXp/fbzhuBeyLivIhYFREXAHcDBzW9v99HxFPAj4CdOqwPYGlEnJHW9xRZQvhyRCxIBwj/BhxepPlP0mvJarw/ioi5ZLceP6Kp2MURMTci/gxcDPw5Is5N5wlnse72/1qqRS4nS6BvGyqOEXgrcFlEXB8RTwOfAvIHCceS3YJ8SZp/IvDWtE2eIfvcXpY+/7kRsbKEGK1kTlDWziPAxDY/glum+Q33Nx6kmsYSsqPaoh7KPX6qxfMN2yz3FuAAYLGkX0vaPU3fFvhIat55TNJjZEfTk9LfAxERufUs7hDbo2Tvt51JLZZfTFYbbHgw9/jJDu+n4f6m582vsRh4HlCk2epo4OqIaHxe59PUzMfwt38+vsUM77MuahJr71f/Q/ZZNGwLXJz7fO8CniXbJucBVwE/lLRU0hckrV9CjFYyJyhr5wayWsch+YnKeq7tD/wyN3mb3PxxZE0xS9OkfCIYUxFxc0TMJGsK+ylZ7QSyH7bPRcSmub8NUu1mGbCVJOVWNbnDy/wCeHN6X60sJfuxzJsMPFDkLRSc3vwak4FVrJ1I1iHphWTNd3+nrAfmg8CHgB2VzhOO0Da5x5NZ81l30uq9/g+wQe75S3OPl7H2frUBWa2o4X5g/6bPeHxEPJBqxidFxDSypuoDyZo5rWacoKyl1Nx2EnCGpP0krZ86K1xIVkM6L1f8f0k6JNW2PkiW2G5M8x4iO28ypiQ9X9LbJW2SmhZXkh1BA3wTeJ+kGcq8SNIbJW1ElnhXAf8i6XmSDgF27fBSXwY2Bs6RtG167a0kfTl1vLgceLmkI9L6DgOmAZcVeBt/JGu2Gmr7XAB8SNJUSRuy5nxOq+bXvDeRbZNpZM2KOwF/DfwXo/vBfr+krdP5t0+QNQMO5SFgSlOin0fWVLm+pOlkzXoNFwEHSnptOm94Mmv/Xp0JfC73mWwuaWZ6vLekv1V2rd1Ksia/EV3WYL3lBGVtpU4NnwC+RPZFn0125LpPavdv+BlwGNmJ6iOBQ1LSAPg88MnUFPNRxtaRwCJJK8lOhL8jxT2H7DzU11JM95J1MiAi/kJWK3xnmncYTR1B8tJ5lj3IfuRmS3qcrPa4Arg3Ih4lO0L/CFkT1P8BDsw1qbUVEU+SncP5Tdo+7c7bfYfsgOB6YCHwZ+C4odZP1pT33XTt2YONP7Lt8vYi57DaOB+4mqwjxwKgyEXOF6b/j0pqnCv8FFlHiz+RHQyd3ygcEXcA70/TlqUyqy8CBk4HLgGuTp/JjWSdSyCriV1Ets/eBfwa8DVYNaS1m+LNhkfSiWQno9/R61isfJIWAe+NiF/0Ohbrf65BmZlZJTlBmZlZJbmJz8zMKsk1KDMzq6SeD9GfN3HixJgyZUqvwzAzsxLNnTv3kYjYfKhyhRJUuuBvckT8btSRdTBlyhTmzJlT5kuYmVmPSeo0estqQzbxSTqI7IK6K9PznSRdMrrwzMzMOitSgzqR7Er76wAiYp5ytz+w4Tl/9h/42by1R8GZudNWHDGj02g7ZmaDp0iCWhURK9YeusyGI5+UZi9cDsCMqRNWP5+9cPnq+U5WZmaZIgnqdklHAOtJ2oHsPiy/LTes/vKzeQ9w57KVTNtyY2ZMnbBWEsonrzuXZXcEcIIyMyuWoI4DTiAbAPQCsmHsTykzqH40bcuNmXXs7utMP2LG5NUJ6bCzblhnvpnZoBoyQaUBLU9If1ZA83mmRu2piDuXrVwrUbnJz8wGVdsEJelSOtzLJyIOLiWiPpBv0oOs9jRzp62GWIp1yrjJz8wGWaca1JfS/0PIhq9vDFf/NmBRiTH1hXZNep3km/vATX5mNtjaJqiI+DWApFMiYs/crEslXV96ZGZmNtCKjMW3uaTVd/yUNBUYcogKMzOz0SjSi+9DwHWSFqTnU4BjS4uoppq7ixftFDGUfKcJd5gws0FSpBfflen6p79Kk+5uut23sXbHiKKdIoaSX4c7TJjZoBkyQUk6qmnSjpKIiHNLiqm2RtIxohNfI2Vmg6xIE98uucfjgX2AWwAnKDMzK02RJr7j8s8lbQKcV1pEZmZmjOyGhU8COxQpKGkR8DjwLNmgs9NH8HpmZjaAipyDyo8oMQ6YBlw4jNfYOyIeGUFslTaa4YxGyj36zGyQFKlBfSn3eBWwOCKWlBRPbYx0OKORco8+Mxs0RRLUARHx8fwESac2T2sjgKslBXBWRJzdXEDSMcAxAJMn1+sHd6x77XXiHn1mNmiKJKiF6G8zAAAPhUlEQVQ3AM3JaP8W01p5TUQslfQS4BpJd0fEWsMkpaR1NsD06dPbDk5rZv2v1R2nO3FTd39rO9SRpH+SdBvwCknzc38LgflFVh4RS9P/h4GLyW4db2bWUqPpvIg7l60cVjKz+ulUgzofuAL4PHB8bvrjEbF8qBVLehEwLiIeT4/3BU4eTbBm1n9aDRNWpOn8sLNucMehPtcpQUVELJL0/uYZkiYUSFJbABdLarzO+RFx5chDNbN+NNJhwtxxqP8NVYM6EJhL1tlBuXkBbNdqodUFIhYAO442wCopa0DYkfCdd62u2l2iMZr7pzXXpsDfiX7Q6X5QB6b/U7sXTrWVMSDsSPjOu1ZnZVyi4e9Ef+p0y/edOy0YEbeMfTjV182u5e34zrtWNyM9z1SUvxP9qVMT32kd5gXw92Mci5n1qV60PrgDRf11auLbu5uBmFl/62brgztQ9IciY/GNB/4ZeC1Zzem/gDMj4s8lx2ZmNdbLTkUeeaU/tL1QN+dc4JXAGcDXyAaL9e02zKyj/EW3vexUZPVVZKijV0REvrv4tZJuLSsgM6unseo+Xgafj6qnIgnqvyXtFhE3AkiaAfym3LCqo0rXPnXiL6D1WrdH+C/K56Pqq0iCmgEcJekP6flk4K40Tl9ExKtKi64CqnLtUyf+AlpVVKXGlOfzUfVVJEHtV3oUFVfFL12ev4DWK3VpYbB6GjJBRcRiSZsB2+TLD+qFuma2Rh1aGJp5SKT6KNLN/BTgncB9rLn1uy/UNRtAVe4IUYSHRKqXIk18hwLbR8Rfyg7GzKqtqh0hivKQSPVSJEHdDmwKPFxyLGZWA3WqMVm9FUlQnyfran478HRjYkQcXFpUPdSuCaNO3MZuY8kdIaxXiiSoc4BTgduA58oNp/fq3oThNnYba3XsCDEcvoawuookqEci4qulR1IhdW7CcBu7jVbdO0IMh68hrLYiCWqupM8Dl7B2E5+7mZv1obq3IgyHryGstiIJ6tXp/265ae5mbtZHyr6hoNlIFLlQ1/eFMutz/X6eqSh3MKqWIjUoJL2R7JYb4xvTIuLksoIys3IN0nmmotzBqHqKjCRxJrABsDfwLeCtwE0lx9VV/d6N1r2UrNkgnWcqyh2MqqdIDWqPiHiVpPkRcZKk04CflB1YN/Vz84Z7KVmDzzMNnw/ueqtIgnoq/X9S0iTgUWBqeSH1Rr9+Wd1LyRr6+UCsDN08uGtucm0VyyAmxyIJ6jJJmwJfBG4h68H3zVKjMrNR83mm0Sn74C7/+cxeuByAGVMnrFNu9sLlzF64fK3PclASVpFefKekhz+WdBkwPiJWlBuWmY1Epx8915pGZyx6+LX7fGZMndB2fa0ONGAwmuoL9eJriIinyV2sW1f9MN6eWSv5ZrxOP3o2PM2JvblWUzS5FE1Kea06bwzKubFhJajhkrQfcDqwHvCtiPj3Ml+vKPdgsn7izg/la04SzTWh5ia4huZa7FgcNAxSx6fSEpSk9YCvA28AlgA3S7okIu4s6zWHY1C/xL4QsRzNR8plbtdOR+U+2OqOfMLq1MGhjFrsIHV8apugJO3cacECY/HtCtwbEQvS+n4IzAQqkaAGkS9EHJ1OP0T5JNHpiLpZ0R+vTueW3JTXW821q27r54POTjWo0zrMKzIW31bA/bnnS4AZBeMac/1+MW4RvhBxdJqbhvPySWKoLsMNw0lkIzl3Yf2v3w862yaoMRiDT61Wu04h6RjgGIDJk8vbqL4GxIZrpN20ix5RF01k4FqStdbvB51FhjpaH/gnYM806TrgrIh4ZohFlwDb5J5vDSxtLhQRZwNnA0yfPn2dBDYcJ116B3cuXdlynk8e23CV3Zmm101D1p/6qYdfkU4S3wDWB/4zPT8yTXvvEMvdDOwgaSrwAHA4cMQI4xw115pa66edeSy4R5zVWb/18CuSoHaJiB1zz38l6dahFoqIVZI+AFxF1s38OxFxxwjjLOQzB72yzNX3nX7bmceCm4Ktzvqth1+RBPWspO0j4j4ASdsBzxZZeURcDlw+ivisRP22M4+EhwMyq64iCepjwLWSFpB1fNgWeFepUZl1iS/atn5W9y7oHROUpHFko5nvALyCLEHdnYY8MusLrjFZP+qHLugdE1REPCfptIjYHZjfpZisR+p+tFWUr4mzQdAPXdDHFShztaS3SGp1XZP1iZk7bbXWD/Wdy1YWvkanbhrNeuAmPbMqK3IO6sPAi4BVkv5M1swXEeHDzj7SD0dbw+FmPRtEdbuspMj9oDbqRiBmZfHtVczqeVnJkE18kn5ZZJpZVeWb9MDNejaYjpgxmVnH7s6sY3evzQFap9HMxwMbABMlbcaasfU2BiZ1ITbrsbo1B3TiJj2z+unUxHcs8EGyZDSXNQlqJdl9nqyP1bE5wMyKq8MBaKfRzE8HTpd0XESc0cWYrALqPsqEu5KbtVeXA9AinSTOkLQHMCVfPiLOLTEuq5i6XSPlMfXM2qvLAWiR222cB2wPzGPNGHwBOEENiLpeke7zTmbFVPUAtMh1UNOBaRExqns1WX3V4RopdyU3G5kqH4AWSVC3Ay8FlpUci9VI1U6wetBXs5FpdQBale93kQQ1EbhT0k3A6kFiI+Lg0qKySqvqCVY36ZmNXv77PXvhcmYvXL5W68S0SRt37d57RRLUiWUHYfVSlROs7qlnNvby3+/mpvNuK9KL79eStgB2SZNuioiHyw3L6qRXzQHuqWdWrubmv24r0ovvUOCLwHVkF+ueIeljEXFRybFZDXS7ua9VrcnNemb9qUgT3wnALo1ak6TNgV8ATlC2TnNf2d1VXWsyGxxFEtS4pia9Ryl2HykbMM3JovkE60iSVbvu4641mfW/IgnqSklXARek54cBV5QXktVVc3t1Prm06g3ULmE1LwcwY+oEwN3HzQaJilx/K+kQ4LVk56Cuj4iLywhm+vTpMWfOnDJWbT3WXBNqTjx5zfOqcJ2VmY0dSXMjYvqQ5dolKEkvA7aIiN80Td8TeCAi7huTSHOcoAbHUN1XnZTM+lfRBNWpie8rwCdaTH8yzTtohLGZ9bz7qplVX6fODlMiYn7zxIiYQzayuZmZWWk6JajxHea9cKwDMTMzy+t0DuoC4FcR8c2m6e8B9o2Iw8Y8GOmPwOJRrmYi8MgYhNMrjr+3HH/v1Dl2cPzDsW1EbD5UoU4JagvgYuAvZLd8h+zWG88H3hwRD45RoGNK0pwiJ9+qyvH3luPvnTrHDo6/DJ1u+f4QsIekvYG/SZN/HhG/6kpkZmY20IoMFnstcG0XYjEzM1utH4csOrvXAYyS4+8tx987dY4dHP+YKzSShJmZWbf1Yw3KzMz6gBOUmZlVUm0TlKT9JP1O0r2Sjm8x/wWSZqX5syVN6X6U65K0jaRrJd0l6Q5J/9qizF6SVkial/4+3YtYO5G0SNJtKb51BlBU5qtp+8+XtHMv4mwm6RW57TpP0kpJH2wqU7ntL+k7kh6WdHtu2gRJ10i6J/3frM2yR6cy90g6untRr379VrF/UdLdad+4WNKmbZbtuJ91Q5v4T5T0QG4fOaDNsh1/p7qhTfyzcrEvkjSvzbK93f4RUbs/YD3gPmA7suuybgWmNZX5Z+DM9PhwYFav406xbAnsnB5vBPy+Rex7AZf1OtYh3sciYGKH+QeQ3ZZFwG7A7F7H3GY/epDsosFKb39gT2Bn4PbctC8Ax6fHxwOntlhuArAg/d8sPd6sArHvCzwvPT61VexF9rMexn8i8NEC+1fH36lexd80/zTg01Xc/nWtQe0K3BsRCyLiL8APgZlNZWYC56THFwH7SFIXY2wpIpZFxC3p8ePAXUA/3uBoJnBuZG4ENpW0Za+DarIPcF9EjHb0ktJFxPXA8qbJ+X38HOBNLRb9B+CaiFgeEX8CrgH2Ky3QFlrFHhFXR8Sq9PRGYOtuxjQcbbZ9EUV+p0rXKf70m3goa+73Vyl1TVBbAffnni9h3R/51WXSF2EF8OKuRFdQanZ8NTC7xezdJd0q6QpJr+xqYMUEcLWkuZKOaTG/yGfUa4fT/otZ9e0P2e1wlkF24AO8pEWZOnwO76b9TVCH2s966QOpifI7bZpX67DtXwc8FBH3tJnf0+1f1wTVqibU3F++SJmekbQh8GPggxGxsmn2LWTNTjsCZwA/7XZ8BbwmInYG9gfen+4Tllf17f984GDgwhaz67D9i6r653ACsAr4QZsiQ+1nvfINYHtgJ2AZWTNZs0pv++RtdK499XT71zVBLQG2yT3fGljaroyk5wGbMLJq+piTtD5ZcvpBRPykeX5ErIyIJ9Ljy4H1JU3scpgdRcTS9P9hsjEbd20qUuQz6qX9gVsiG9JrLXXY/slDjWbT9P/hFmUq+zmkDhsHAm+PdMKjWYH9rCci4qGIeDYingO+Seu4KrvtYfXv4iHArHZler3965qgbgZ2kDQ1HQkfDlzSVOYSoNFj6a1kI7P3/Ogltfl+G7grIr7cpsxLG+fLJO1K9jk92r0oO5P0IkkbNR6TnfC+vanYJcBRqTffbsCKRnNURbQ9cqz69s/J7+NHAz9rUeYqYF9Jm6VmqH3TtJ6StB/wceDgiHiyTZki+1lPNJ1PfTOt4yryO9VLrwfujoglrWZWYvv3qnfGaP/Ieon9nqyXzAlp2slkOzxk97O6ELgXuAnYrtcxp7heS1bNnw/MS38HAO8D3pfKfAC4g6zXz43AHr2Ou+k9bJdiuzXF2dj++fcg4Ovp87kNmN7ruHPxb0CWcDbJTav09idLpsuAZ8iOzN9Ddk71l8A96f+EVHY68K3csu9O34N7gXdVJPZ7yc7PNL4DjR63k4DLO+1nFYn/vLRfzydLOls2x5+er/M7VYX40/TvNfb5XNlKbX8PdWRmZpVU1yY+MzPrc05QZmZWSU5QZmZWSU5QZmZWSU5QZmZWSUPe8t3MipPU6PoN8FLgWeCP6fmTEbFHTwIzqyF3MzcriaQTgSci4ku9jsWsjtzEZ9Ylkp5I//eS9GtJP5L0e0n/Luntkm5K997ZPpXbXNKPJd2c/l7T23dg1l1OUGa9sSPwr8DfAkcCL4+IXYFvAcelMqcD/xERuwBvSfPMBobPQZn1xs2RxiaUdB9wdZp+G7B3evx6YFruNmYbS9oosvuImfU9Jyiz3ng69/i53PPnWPO9HAfsHhFPdTMws6pwE59ZdV1NNnAtAJJ26mEsZl3nBGVWXf8CTE93bb2TbMR1s4HhbuZmZlZJrkGZmVklOUGZmVklOUGZmVklOUGZmVklOUGZmVklOUGZmVklOUGZmVkl/X/VR9pmV3AJ1QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11e5e2160>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig1 = plt.figure()\n",
    "ax1 = fig1.add_subplot(2, 1, 1)\n",
    "ax1.set_title(\"Initial Control amps\")\n",
    "ax1.set_ylabel(\"Control amplitude\")\n",
    "ax1.step(result.time, \n",
    "         np.hstack((result.initial_amps[:, 0], result.initial_amps[-1, 0])), \n",
    "         where='post')\n",
    "\n",
    "ax2 = fig1.add_subplot(2, 1, 2)\n",
    "ax2.set_title(\"Optimised Control Amplitudes\")\n",
    "ax2.set_xlabel(\"Time\")\n",
    "ax2.set_ylabel(\"Control amplitude\")\n",
    "ax2.step(result.time, \n",
    "         np.hstack((result.final_amps[:, 0], result.final_amps[-1, 0])), \n",
    "         where='post')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Versions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table><tr><th>Software</th><th>Version</th></tr><tr><td>QuTiP</td><td>4.2.0</td></tr><tr><td>Numpy</td><td>1.13.3</td></tr><tr><td>SciPy</td><td>1.0.0</td></tr><tr><td>matplotlib</td><td>2.1.0</td></tr><tr><td>Cython</td><td>0.26.1</td></tr><tr><td>Number of CPUs</td><td>2</td></tr><tr><td>BLAS Info</td><td>INTEL MKL</td></tr><tr><td>IPython</td><td>6.2.1</td></tr><tr><td>Python</td><td>3.6.3 |Anaconda custom (64-bit)| (default, Oct  6 2017, 12:04:38) \n",
       "[GCC 4.2.1 Compatible Clang 4.0.1 (tags/RELEASE_401/final)]</td></tr><tr><td>OS</td><td>posix [darwin]</td></tr><tr><td colspan='2'>Thu Apr 19 19:07:18 2018 PDT</td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from qutip.ipynbtools import version_table\n",
    "\n",
    "version_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[1] Doria, P., Calarco, T. & Montangero, S.: Optimal Control Technique for Many-Body Quantum Dynamics. Phys. Rev. Lett. 106, 1–4 (2011).\n",
    "\n",
    "[2] Caneva, T., Calarco, T. & Montangero, S.: Chopped random-basis quantum optimization. Phys. Rev. A - At. Mol. Opt. Phys. 84, (2011)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "tlist = np.linspace(0,evo_time,n_ts)\n",
    "H = [[H_d, np.ones(n_ts)], [H_c[0],result.final_amps[:,0]]]\n",
    "sim = mesolve(H,psi_0,tlist,e_ops=[psi_targ*psi_targ.dag()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 1)"
      ]
     },
     "execution_count": 91,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xd4nFeZ9/HvmRn13qttSbbc4xbFTiG9B9KoSVgIZTehhIVlYWHh3cBmaS9tgX1DILQQYCEQShJIsJPgFJK4x3YsS7bVLMlqo2JpZqQZTTnvH1M8kkbWSJoije/PdfmK9DzPzBxPxj8d3c8pSmuNEEKIxGKIdwOEEEJEnoS7EEIkIAl3IYRIQBLuQgiRgCTchRAiAUm4CyFEApox3JVSP1VK9SmljkxzXimlvqeUalJKHVZKbYl8M4UQQsxGOD33R4AbznL+RqDW9+ce4KH5N0sIIcR8zBjuWuuXgMGzXHIr8Kj22gXkKqXKItVAIYQQs2eKwHNUAB1B33f6jnVPvlApdQ/e3j0ZGRnnr169OgIvL4QQ5479+/f3a62LZrouEuGuQhwLuaaB1vph4GGAuro6vW/fvgi8vBBCnDuUUifDuS4So2U6gSVB31cCXRF4XiGEEHMUiXB/Enivb9TMhcCw1npKSUYIIUTszFiWUUr9GrgCKFRKdQJfAJIAtNY/AJ4GbgKagFHg/dFqrBBCiPDMGO5a6ztnOK+Bj0asRUIIIeZNZqgKIUQCknAXQogEJOEuhBAJSMJdCCESkIS7EOKcZne6sTvdMXu9pj4LHk/0966WcBdCnNPu/cV+7v3F/og/b2u/bcoxs8XBNd9+iZ/8vTXirzeZhLsQ4pzlcLl5rWWAl06Y6To9FrHnfem4mSu/+QKvNvVPOL6vzbsG4/lVeRF7relIuAshzllHTg0z7vKgNTx1aPpVU4ZHnTy2t53f7u3gDwc6OdRx+qzP+/j+TgBePG6ecHx36yCpSQbWl+fMv/EziMTCYUIIsSjtaxsCoKYwgycOdnHv5cunXKO15sO/2s+rzQMTjj/x0UvYuCR3yvVWh4sdR3sApjxmb9sgm5fkkWyKfr9aeu5CiHPW3rYhqgszeM9FyzjaPcKJXsuUa55v6OPV5gE+e+Nq/v6ZK3n2Xy6jICOZLz/dgHeC/kQ76nuwOz1cWlvIka5hhkedAIzYnTR0j7C1Oj/qfy+QcBdCJLAHdzZx5NRwyHNaaw60D3H+sjzevKEMg4InDk4szTjdHr7yTAM1RRl88E3VVOalU1uSxSeuXcme1kF2HO2d8rx/fP0UlXlp3HflCrSG3a3e3vv+k0N4NBLuQggxH239Nr6x/Rjfee54yPMt/TYGbeNcUJVHcVYql6wo5IlDpyb0xv93dzstZhufu3ENScYzcXnnBUtYXpTB155pxOn2BI73Wey80tTPbZsq2LQ0l9QkQ6A0s7d1EJNBsXnp1FJONEi4CyES0gvH+gDvTU1/aSRYYOTKMm9P+tZNFXQMjnGg3XuzdHjMyXeeO85FNQVcvaZ4wmNNRgOfu2kNrf02frXrzN4ZTx3qxqPhts3lpJiMXFCVz2v+cG8bZH1FDunJsbnVKeEuhEhIO4+ZyUwx4XRrttf3TDm/r22IvPQklhdlAHD9uhJSTAa+9Jej3PPoPt7yPy9zeszJ59+8BqWmbjh31epiLl5ewLefPc6PX25hwOrgiYOnWF+RzYriLAAuWl7AsV4LnUOjHOoYjllJBiTchRAJaGzcO379HXWVVBWk82SIYY77Tw5x/rL8QHBnpSZx++YKjpwapqXfxprSbL759o2srwg9bFEpxZduW09NUSZf+ksD277yPIc7h7ltU0XgmouXFwLwwxdbGHd72FoVu3CXoZBCiITzWks/4y4PV60uJjPFxIM7mzBbHBRlpQDQb3XQ0m/jnRcsmfC4r71tA1+5/TwMhlBbQ09VU5TJnz56Ccd7LfxuXweHOoe5bfOZcF9fnk1WionH9nYAUBeDyUt+0nMXQiScnY1m0pKMbK3O5+aN5Xg0PHPkzO6f+096x7fXLZsatuEGe7CVJVl8/s1r+e29F1GYmRI4bjIa2FaTz7jbw6qSLHLTk+fwt5kbCXchRELRWrPzWB+XrCggxWRkZUkWq0qyJsxA3dc2SLLJwHmV0Z8pemFNARC7IZB+Eu5CiITSbLbSOTTGFavOjHC5eWMZe9uGaOge4RevtfHEwS42VOSQYjJGvT1XrCrCoOCylUVRf61gUnMXQiSUnY3e9VyuWHUmTN+yoZxv7jjOjd99GYA1Zdl8+vpVMWnPiuIsXvv3qynOSpn54giScBdCJJSdx/pYWZJJZV564FhVYQb3XFbDuMvD28+vnHYETLSUZKfG9PVAwl0IESdNfRaazTauX1c6p8drraeMP7c73extG+T9l1RPuf5zN62Z0+ssVlJzF0LExfd3NvPPv34d9xx2JWrrt7H6P/7KG50T141p7bfhdGvOi3HPfCGScBdCxEVLvw2Hy0P74OisH7unbRCHy8OB9qEJx9t8ux9VF2ZEpI2LmYS7ECIu2ga8QRxqmd2ZNHSPANBitk443up7zioJdwl3IUTsDdnGOe1bzOtEn3WGq6cKhPukfUrb+m0UZqaQmSK3EyXchRAx5+9hw+x77lprGrq9j2kxTw73UaoL00M97Jwj4S6EiDl/bbwiN21Kz91id/LUoS4OdpzG6nBNeWz3sJ3hMSfFWSl0DY9hd7oD51oHbFQVSEkGZCikECIO2vptGBRcvaaYx/Z24PZojL41XX72ShvffvbMBhvLizL4/YcvDqzL4i/J3HReGY+82kbbgI3VpdlYHS7MFofU232k5y5EAnG43IzYp25MsdC09NuozEtnbVk2DpeHzqEzI2Zeax5gZUkmD7/nfO69vIZms42dvo034Ey437jeOz7eX5rx/zYgPXcvCXchoszudGOxO7HYnYyNu2e8PtSmy+H69O8Oc9uDr8z58bHSNmCjqjCD2hLvphYner2lGYfLzYH2Id60oojr1pXymetXU5iZzN98SwoANHRbWJqfHphl2uoL9ZMD3h8QVVJzB6QsI0RU7Wsb5J0/fA3/PB2lvD3O+66sZW15Nlprmvqs/K2xjyNdIxzrGaGtf5Tr1pXw1beeR1ZqUtiv1dZv46nDXWgNnUOjE6bfLyRaa9r6R6lbls+K4kzAO2LmmrUlHO4cxuHysK3Gu4KiwaC4fGUxzzX04nJ7MBkNNHSPsKYsi4wUE6XZqTT7hkP6h1ZKz91Lwl2IKPrTwVOkmIx88tqVKAW9I3Z+s6eDp9/o4U0rCukYGg30OCty01hTlsXmJXk8fqCT+q4RHrxrC2vLs8N6rR+93BL4em/b4IIN937rOFaHi6qCdHLSkijNTg2MmNnd4t1vNHjHoqtWF/P7A5283nGadeXZtA7YuGVTOQA1RRmBnntrv43irBQyZBgkIOEuRNR4PJod9b1cvrKIf7qsJnD8vitr+dmrrTy+v5MVxZncc1kNV68uoTTnzOJSb91Swcd+/Tq3f/8VvvXOjbxlQ/lZX6vf6uDx/Z284/xKnjnSw57WQW7fXBm1v9tMRsddGJQiNWnqkrr+MPbf+KwtyQyMmNndOsjq0izyMs5sanHpykJMBsXfGvswGRRae1d1BG+4P3Wo2/fbgE1upgaRcBciSg52nqbP4uD69SUTjuekJ/GJa1byiWtWTvvYbTUFPP3xS/nwL/fz8d8cxGQwcMP66RfYevTVNsbdHu69fDlmi4M9rYMR+3vMltaaOx7eRVVBBt+7c/OU85OXCFhRnMlv9nTgcLnZ1zbEuyZtfZedmkRdVR47G/tY4vttZK0v3KsLMxkeczI06qRtwMbVqye+1+eysG6oKqVuUEodU0o1KaU+G+L8UqXUTqXU60qpw0qpmyLfVCEWl+31PZgMiqtWzS1wCjNT+Nn7t7KhMoeP/frAhBEjwUbHXTy66yTXrilheVEmW6sLaDbb6Lc6Ql7v8Wg65rCeS7iOnBrhcOcwf2/qD3lzuHXARpJRUZGbBni3qBtzuvnrkR7GnG62hdix6MpVxTT2WPhbYx9ZKSYq87yPrSny/oA41Hmafuu49NyDzBjuSikj8CBwI7AWuFMptXbSZf8H+K3WejNwB/D9SDdUiMVEa29J5qLlBeSkh39TdLLMFBOPvH8rK0uy+NAv9vPc0d4pr/OzV9o4Perk3su9pZ+t1d59Qfe1he69f+vZY1z5zRemDf/5+v2BTgAGbeOB+wnB2vptLMlPx2T0xk+t76bqL3edBEJvR3fVau+uSs819LK6LCuw1G+NL8x3Nnp/8Mns1DPC6blvBZq01i1a63HgN8Ctk67RgP+uTw7QhRDnsBN9Vlr7bVw3x7XKg+WkJfGLD26jujCDf3x0H3f/dA+NPSMcOTXMXT/azTe2H+PS2kLOX+YNxfMqckkxGdjTOjTluY71WPjhiy24PDpk8M7XuMvDEwdPsbrUO8Rx8qqN4K25VweNaKkt9l67t22I2uJMCjKn7li0ojgz0Fv319sBKvPSSTIqnm/whvsyGSkTEE64VwAdQd93+o4F+yLwD0qpTuBp4GOhnkgpdY9Sap9Sap/ZbA51iRAJYfuRHgCuWxuZGnB+RjJP3HcJn79pDa+3D3HTd1/m5v/3dxp7RvjPW9bx0/ddELg22WRg89Jc9rQNTHgOj0fzuT++gX9/i+7hsYi0LdjOY30MjTr59PWryEwxTQl3j0cHxrj75aQnBbag828mPZlSiit9e6IGh7vRoFhWkMGp096/iwyDPCOccFchjk0upN0JPKK1rgRuAn6hlJry3Frrh7XWdVrruqKi2G4WK0QsbT/aw+aluRHdXi3FZOSfLqvhxU9fyUeuWMHHrqrlxX+7krsvriLJOPGf29aqfI52jWAJmq362L4O9p8c4t9v9O5I1HU68uH++P5OirJSuHxlEZuW5HLg5OkJ53stduxOz5T11mtLvKUZ//j2UN68oQyjQVG3LG/CcX9ppjQ7lbTk6G94vViEE+6dQPDt60qmll0+CPwWQGv9GpAKFEaigUIsNqdOj3Hk1AjXrZ1/SSaUvIxkPnX9Kj557Uqyp5nktLW6AI+GA+3ecDVbHHztmUa2Vefz/kuqyEox0XXaHtF2DVgd7Gzs47ZN5ZiMBrYszaWxZ2TC4l+t02ym4S/NhKq3+11YU8DB+68NzGr1q/bdVJWZqROFMxRyL1CrlKoGTuG9YXrXpGvagauBR5RSa/CGu9RdxDnJf9PzunXxG5a3eWkuRoNiT+sAY+NuvvDkEcbG3Xz59vNQSlGemxYoZUTKk4e6cHk0bzvfO75+87I8PBoOd5zm4hXevt7kMe5+H3xTNRuX5FCcdfbfdELN2F1e6O31y+5LE80Y7lprl1LqPmA7YAR+qrWuV0o9AOzTWj8J/CvwI6XUv+At2bxPz2eBDCEWsReO9VFVkM7yosy4tSEjxcT6ihx+/HIrDpeHNWXZPPyeusB0//Lc1IiXZf50sIv1FdmsLvXWxLcs8ZZPDrQPBcK9rd9GislA2aRy1ZL8dJbkz63n7R8OKfX2icKaxKS1fhrvjdLgY/cHfX0UuCSyTRNi8bE73bzWMsAdFyyNd1O4enUxx3pG+NxNq/nAJdWBoYcA5blpHOw4fZZHz15zn5W3n39mVmxOehLLizICpSG70832+l7WlGVjMIS6lTc368pzuGZNMVevKY7YcyYCmaEqRATtbh3E7vRw+ar4Dxj46JUruOeympBLAJTnpjE06l2lMhI3IcfG3VgdLoqzJw5j3LI0j+caetFa89ALzbQPjvK1t54379cLlpZs5Md3XzDzhecYWfJXiAh64VgfKSYDF00zpC+WjIbQa7uAtywD0BWh4ZD+CVFFk8aob1mWx9CokxeOmXnoxWZu2VgeKNGI6JJwFyKCXjxm5sKagmlDdaEoz/FOCIpU3b3P4g33wqypPXeAj//mdZKNBv7Pm9dE5PXEzCTchYiQ9oFRWvptXLEASjIzKc+NbLibLaF77rXFmWSlmBixu/jktSspjuC4f3F2UnMXIkJeOO6dAn/FqoV/Y680JxWl4FSExrqbfWWZ4kk9d4NBccmKQrqGx3jvRcsi8loiPBLuQkTIC8fMLCtIXxTjrZOMBkqyIjcc0mxxoJR3mYTJ/ueuzbg9esJoHRF98m4LEQF2p5tXm/sD658sBmURHOtutjgoyEgOGeBJRsOCvweRiCTchYiAhTQEMlzluWl0D0eoLGNxUBhiNUcRPxLuQkTAM290k55s5MLq+A+BDFeFbwmCUJPJ7U43X3yyPnCj1M/p9nDDd17imTe6JxzvtzooypJwX0gk3IWYJ7vTzV8Od3Pj+rJFtSpheU4q4y4PA7bxKecOnBzikVfb+OuRiSHebLbS2GNhV8vE5YTNFseUkTIiviTchZinHUd7sThcvO38ydscLGxnGw7Z7Fvg62i3ZcLxRt/37UHb9GmtMUvPfcGRcBdinv5woJPynNRFVZKBs4d7i9kKQGPPyITjjT1Tw33E7mLc5ZFwX2Ak3IWYh74ROy8dN3P7loqILoYVC/5wDzXWvdns7bkf67Hg8ZypyfvDvmNoLHA8MIFJwn1BkXAXYh6eONiFR8Nbt1TOfPECk5eeRGqSYdqee7LRwOi4e0IvvbHbgsmgGHd56LV4fyhMNztVxJeEuxBzpLXm9wc62bQkN65rt8+Vf9OOyXup2p1uTp0e47KV3gW+Grq9vfXTo+P0jNgD+5y2+zbY9s9OlZ77wiLhLsQcHe0eobHHEth5aDHyDoecWJZp7behNdywvgyDggZfnd1fb/fvMOXv0ff7Fw2TnvuCIuEuxBxtP9KD0aC4eUNZvJsyZ+U5aVPKMs2+m6lry7KpLswI9Nwbff+9anUxBgUdg2d67klGRU5a6P1cRXxIuAsxR01mK8vy08lNn7qeymJRnpuG2eLA4XIHjrWYz2xivbosO3ATtbHHQl56EhW5aZTlpAV67v7ZqYvthnKik3AXYo5azLYpGz0vNmX+TTuCSjPNZisVuWmkJRtZW5ZNx+AYFruTxh4Lq0uzUUqxND99QrhLvX3hkXAXYg601pwcGF30mzJvrMwF4LXmMzNOW8y2wKbTq0uzAGjotnCsx8Iq3/fLCtJpH/SWc2R26sIk4S7EHPSOOBhzuqkuWtzhvrIkkyX5aTzX0At4f2i1mK2B0T9ryrIB2FHfw5jTzZoyb7gvyU+n3+pgdNwls1MXKAl3Ieagpd9707F6kffclVJcvbqEV5r6GR130TviwDbuZrnvh1ZZTirZqSaeONQFwOpSb9gvzU8HoK1/lAEJ9wVJwl2IOWjr99abqwrT49yS+bt2bQkOl4e/n+gPLDvg77krpVhTlh3YjGNlibfn7g/3w52n8WgZBrkQSbgLMQdtAzaSTYbARtOL2dbqfLJSTTzX0BsYBlkTNCnLX5qpKsgIrHrpD/f9J4cAmcC0EEm4i0XnyKlhrvrmC/SNRGajibloMduoKkhPiOF/SUYDV6wq5vmGPk70WclINlKSfSas/XV2/81VgNz0JLJSTOxvl3BfqCTcxaLz0gkzLf02/lrfM+G41pq/HumZMGY7EsbG3QyPOSccaxuwLfqRMsGuWVPMgG2cvxzupqYoE6XO/NDy99z99XbwlmuW5KcHxsTLaJmFR8JdLDoner2lg+ca+iYcf+GYmQ/9cj/f39kc0df71OOHeNcPXwt87/Zo2gdGF8VG2OG6YmUxJoNiwDYeuJnqt6Ysm3+4cCm3bCqfcNxfmgHpuS9EEu5i0TnR513jZFfzAFaHK3D8j6+fAuDHL7cwYHWEfOxsDY85eba+l8YeS2C6fdfpMcbdnoQK95z0JLZW5wMT6+3gLdt86bbzpvx9lxZ4wz092UhGiik2DRVhk3AXi4rHo2nqs7K+Iptxt4eXj5sBsDpc7Djaw5tWFDLmdPPQC5Hpve+o72Hc7QHgRd9rtfp2KVrss1Mnu3qNd0GwmjDH7i/x9dxlpMzCJOEuFpWOoVHsTg93XLCUnLSkQGlmR30PdqeHj19Ty9u2VPLorpMh1ymfracOd7MkP42K3DRe8oV728CZtVcSye2bK3jrlgretKIwrOv9ZRkpySxMEu5iUTnuq7evKcvmylVF7DzWh9uj+dPBLirz0jh/aR4fv6YWNHzv+RPzeq0Bq4NXmvq5eUM5l60s4tXmAZxuDy1mG+nJRooTLNTyM5L59js3hb0QWiDcpee+IEm4i0XFX2+vLcnk6jUlDNrG2VHfw99PmLl1UzkGg6IyL527ti3ld/s7A+O25+KZIz24PZqbN5Zz+cpCrA4Xr7efDoyUCR5Rci6qyE3DoKTnvlBJuItF5USv1TclPonLVxVhMijuf7Iej4bbNlUErrvvqhWkJxm5/4kjaK3P8ozTe+pQFyuKM1ldmsXFKwoxGhQvHTfT1m9LuJLMXCSbDHzl9vN494VL490UEYKEu1hQDnWc5kcvtUwbyMd7Lawo9o7myE5NYltNPmaLg3Xl2dSWnJlkU5iZwmduXM0rTQM8vr9z1u3oGbazp22QmzeUo5QiOzWJzUtyeb6xj46hMQl3nzu2Lp0w/l0sHGGFu1LqBqXUMaVUk1Lqs9Nc806l1FGlVL1S6n8j20xxrvjBi818+ekGttf3Tjnn9o2UWRkU4tf4RngE99r97tq6lAuq8vjSXxoCmziH68+Hu9Aa3rLxzC5Ll60soqF7BLdHJ9xIGZF4Zgx3pZQReBC4EVgL3KmUWjvpmlrg34FLtNbrgE9Eoa0iwWmt2dM6CMB/PlU/YQw7QOfQKA6Xh5UlZ8Zh37apgru2LeUddVP3MTUYFF996wbGxt088OejYbejZ9jOj19uZX1F9oSNry9fWRT4ujoBFgwTiS2cnvtWoElr3aK1Hgd+A9w66Zp/Ah7UWg8BaK37EGKWms02BmzjvKtuCT0jdv772eMTzvtHygSXX/IykvnK7edNO8JjRXEm9121gqcOdfGDF5tx+sasT8did/K+n+3B6nDx9bdtnHBufUUOeenefUKrCzNDPVyIBSOccK8AOoK+7/QdC7YSWKmUekUptUspdUOoJ1JK3aOU2qeU2mc2m+fWYpGw9rZ5e+33Xl7DnVuX8rNXWjlyajhw/nivd6SMv+Yerg9dvpyrVhfztWcaufG7LwcmI03mdHv4yK8O0NRn5fvv3sLa8om1ZKNBcfnKIgoykgMhL8RCFc6c4VDjvSbf7TIBtcAVQCXwslJqvdb69IQHaf0w8DBAXV3d3IYwiIS1p3WQwswUqgsz+Mz1q9lR38Pn//gGv//wxZiMBk70WgIjZWYj2WTgJ3fX8VxDH1/+y1Hu/ukeKnLTWFuezdqybIwGRfewnaPdIxzqOM033r6By4JKMMH+4y1rMVsd5/wwSLHwhRPuncCSoO8rga4Q1+zSWjuBVqXUMbxhvzcirRTnhD2tg2ytzkMpRU56El+4eR0f+/Xr/Pdzx/n09as50WedUJKZDaUU164t4bKVhfx2Xyd7Wgc52jXMcw29aA2FmcmU5qTyn7es4x11S6Z9noLMFApk0o5YBMIJ971ArVKqGjgF3AHcNemaPwF3Ao8opQrxlmlaItlQkdg6h0Y5dXqMf7q0OnDs5o3l/P1EP99/oZm6qnya+qxcVFMwr9dJMRl5z4XLeM+FywCwO90o5T0uRCKZseautXYB9wHbgQbgt1rreqXUA0qpW3yXbQcGlFJHgZ3Ap7XWA6GfUYip/PX2rdUTw/uLt6xjZXEW9/3qAA6Xh9qSyN7ITE0ySrCLhBTWOp1a66eBpycduz/oaw180vdHiFnb0zpIdqqJVaUTyy5pyUYefPdmbv6fVwDmXJYR4lwjM1TFgrC7dZALqvIxhti2bkVxFl9/+wbWlGVP2OpNCDE9WWFfxJ3Z4qDFbOOdZ7mRefPGcm7eWD7teSHERNJzF3G3L1Bvz49zS4RIHBLuIu72tA2SmmRgfXlOvJsiRMKQcBdx19ZvY0VxJskm+TgKESnyr0nEXfewndLstHg3Q4iEIuEu4q53xE5pjsz6FCKSJNxFXNmdboZGnZTlSM9diEiScBdx1TNsB6A0OzXOLREisUi4i7jq9oV7WY6EuxCRJOEu4qpnZAyAUgl3ISJKwl3Elb/nLuEuRGRJuIu46hm2k5OWRHqyrIQhRCRJuIu46h62S71diCiQcBdx1TNsl5KMEFEg4S7iSnruQkSHhLuIm3GXh36rQ5YeECIKJNxF3PSO+EfKyNIDQkSahLuImzPhLj13ISJNwl3EjcxOFSJ6JNxF3PTIBCYhokbCXcRN97CdjGQjWSkygUmISJNwF3HTMzJGaU4qSql4N0WIhCPhLuLGO8ZdbqYKEQ0S7iJuZHaqENEj4S7iwuX20GdxyEgZIaJEwl3ERb91HLdHS89diCiRcBdx0T3s3aRDeu5CRIeEu4gL/xj3Etk7VYiokHAXcXFmdqqMlhEiGiTczzG7Wwb4/gtNaK3j2o7eETvJJgN56UlxbYcQiUqmBp5DXm8f4n0/28uY081FNQVsXpoXt7b413GXCUxCRIf03M8RzWYrH3hkL0VZKWSmmHj0tZNxbU/viJ2SLKm3CxEtEu7ngN4RO+/9yR6MBsWjH9jK28+v5C+HuzFbHHFrU7/VQWFWctxeX4hEJ+F+DvjWjmMM2sb52fu2UlWYwT9cuIxxt4fH9rbHrU391nEKMmSTDiGiJaxwV0rdoJQ6ppRqUkp99izXvV0ppZVSdZFropiv3a2DXFpbyHmVOQCsKM7k0tpCfrmrHZfbE/P2jLs8DI85KcyUcBciWmYMd6WUEXgQuBFYC9yplFob4ros4J+B3ZFupJi7vhE7JwdGuaAqf8Lx915URc+InWeP9sa8TYO2cQApywgRReH03LcCTVrrFq31OPAb4NYQ1/0X8HXAHsH2iXnad3IIgLqqiSNjrlpdTEVuWlxurPZbvbV+6bkLET3hhHsF0BH0fafvWIBSajOwRGv957M9kVLqHqXUPqXUPrPZPOvGitnb0zpIapKBdeU5E44bDYo7LljCay0DMb+xapZwFyLqwgn3UAORAzNglFIG4L+Bf53pibT2C3uZAAAWkElEQVTWD2ut67TWdUVFReG3UszZvpODbFqSS7Jp6v/qOl+ppr5rOKZt6rf4w13KMkJESzjh3gksCfq+EugK+j4LWA+8oJRqAy4EnpSbqvFndbg42jXC1kn1dr+15dkA1HeNxLJZ9Ft9NXfpuQsRNeGE+16gVilVrZRKBu4AnvSf1FoPa60LtdZVWusqYBdwi9Z6X1RaLML2evsQHn2mhz5ZTloSS/PTY95zH7A6SEsykiF7pwoRNTOGu9baBdwHbAcagN9qreuVUg8opW6JdgPF3O1tG8KgYPPS3GmvWVeeHYeeu0xgEiLawuo6aa2fBp6edOz+aa69Yv7NEpGwr22QNWXZZKVOvzjX+oocnjnSw4jdSfZZroukfuu4lGSEiDKZoZqgnG4Pr7efnjK+fTJ/3f1oDHvv/VaHhLsQUSbhnqCOdo0w5nRPGd8+2bo43FT1hruUZYSIJgn3BLW3bRCAumVn77kXZ6VSnJVC/anY3FR1ezSDNinLCBFtEu4J6vX201TmpYW1AXUsb6oOjY7j0TIMUohok3BPUCf6LKwuzQ7r2vUVOTSZrdid7ii3SpYeECJWJNwTkMvtoa1/lOXFGWFdv648G7dHc6zHEuWWQb/FP4FJau5CRJOEewLqHBpj3O1heVFmWNf71505EoPJTP6ee4H03IWIKgn3BNRstgKEHe6VeWlkp5piUnf3h3uRhLsQUSXhnoDOhHt4ZRmlFOvKc2IS7marg2Sjgew0WXpAiGiScE9AzX02CjOTyU0Pv669viKbxu6RqO/M1G8ZpyAzGaVCLTYqhIgUCfcE1Gy2UhNmScZvTVk2DpeHtgFblFrlNWCT2alCxIKEewJqNlvDrrf7rSjO9D02uuHeb3VQICNlhIg6CfcEM2gbZ2jUGXa93a+60Hu9v14fLf0WmZ0qRCxIuCeYwM3U4tn13LNSkyjJTqElij13rbWUZYSIEQn3BNPc5w33FbMsywDUFGZGtec+PObE6dYygUmIGJBwTzDNZispJgPluWmzfuzy4gxazDa01jNfPAeBMe5Z0nMXItok3BNMs9lGdWEGRsPshxrWFGYyPOZkwDYehZad2Tu1IEPCXYhok3BPMM1m66zr7X7+x0Wr7h5YNEy22BMi6iTcE4jd6aZjcHTWwyD9aqI8YqbfIitCChErEu4J5OTAKB4d/rIDk1XkppFiMtASrXC3jmNQkDeLmbNCiLmRcE8gs10wbDKDQVFdmBG1iUz9Vgf5GSlzuh8ghJgdCfcE4h8GWTPHnjt46+7R67nL3qlCxIqEewJpNlupyE0jPXnuKy4uL8ygfXAUhyvyuzKZrTI7VYhYkXBPIK39tsAyAnO1vDgTj4b2gdEIteoM84id4mwJdyFiQcI9gbQNjLKsIH1ez1FT6F9ALLKlGY9H02dxUJo984bdQoj5k3BPEKdHxxkec1JVML+eu79eH+mbqgO2cVweTYmEuxAxIeGeIE76yihL59lzz0gxUZqdGvGee++IHYASKcsIERMS7gni5KA33Ofbc4cza8xEUp/FG+7F0nMXIiYk3BPEyX5vGC/Nn1/PHc6sDhnJBcR6R7yzU6UsI0RsSLgniJODo5Rkp5CWbJz3cy0vysBidwUW+oqEnmFfz11WhBQiJiTcE8TJARvLIlCSgTMLiDX1Ra7u3mexU5iZTJJRPnJCxIL8S0sQJwdGWRaBkgyc2U+1KYI3VXtHHBRnSUlGiFiRcE8Ao+Mu+iwOquY5gcmvNDuVzBQTTb2WiDwfeEfLyEgZIWJHwj0BtPtGykTiZiqAUorlxZkR77mX5kjPXYhYCSvclVI3KKWOKaWalFKfDXH+k0qpo0qpw0qp55VSyyLfVDGdtv7IDYP0qy3O5ERvZMLd6fYwYJOyjBCxNGO4K6WMwIPAjcBa4E6l1NpJl70O1GmtNwCPA1+PdEPF9NoHfcMg5zmBKdiK4kz6LA6Gx5zzfi6zxYHWMgxSiFgKp+e+FWjSWrdorceB3wC3Bl+gtd6ptfavNLULqIxsM8XZtA2MkpeeRE5aUsSeszaCI2ZkdqoQsRdOuFcAHUHfd/qOTeeDwDOhTiil7lFK7VNK7TObzeG3UpxV+8AoSyNYkoGgETN987+pKhOYhIi9cMI91LY5IacuKqX+AagDvhHqvNb6Ya11nda6rqioKPxWirNqG7BRFcGSDEBlXjopJkNEeu7+pQck3IWInXDCvRNYEvR9JdA1+SKl1DXA54FbtNaOyDRPzGTc5aHr9FjExrj7GQ2KmqJMTkQg3HuG7RgNioIM2YVJiFgJJ9z3ArVKqWqlVDJwB/Bk8AVKqc3AD/EGe1/kmymm0znk3RQ7UrNTg9UWZ0ao5u6gOCsFg+ydKkTMzBjuWmsXcB+wHWgAfqu1rldKPaCUusV32TeATOB3SqmDSqknp3k6EWH+1SDnu0lHKLXFmXQOjTE67prX8/RZ7LIapBAxFtZmm1rrp4GnJx27P+jrayLcLhEm/2qQ0ei5+2+qNvfZOK8yZ87P0ztin/f2f0KI2ZEZqovcycFR0pONFGZGvp5dW+JfY2Z+I2Z6RxxyM1WIGJNwX+RODoyyrCADpSJfz15WkIHJoOY1U9XudDM85pRwFyLGJNwXubYBW8RHyvglGQ1UFWbM66aqfwKTrOMuRGxJuC9iTreH9oHRwKbW0bCiaH4jZmQCkxDxIeG+iLUPjuLyaJYXZUbtNWpLMmkbsOFwuef0eH/PXVaEFCK2JNwXMf8m1lHtuRdn4tFnVp6crcC6MrIipBAxJeG+iDX71luviWLPfVVpFgAN3SNzenzviJ0Uk4HstLBG3QohIkTCfRFrMVspzEyJ6GqQk60oyiQ1ycChztNzerx/GGQ0RvMIIaYn4b6INZttLI9iSQbAZDSwvjyHw53Dc3p874idUrmZKkTMSbgvYi1ma1RLMn4bKnOp7xrG6fbM+rHdw3aKZR13IWJOwn2RGrSNMzTqjHrPHWDjkhzsTg/HZ7lhtt3ppmNoNCY/gIQQE0m4L1Itvpup0RwG6behMhdg1qWZpj4rWsPKEgl3IWJNwn2Rao5huFcVpJOdauLwLG+qnvDt4rSyJCsazRJCnIWE+yLVYraRbDJQkZcW9ddSSrFxSS6HOmbXcz/ea8VkUFRFYcVKIcTZSbgvUs1mK9UFGRhjtAHGhsocjvVasDvDn6l6otdCdWEGySb5mAkRa/KvbpFqMduiOjN1sg2Vubg9mvqu8CczneizSklGiDiRcF+Exl0eTg6OxqTe7rcxcFM1vLr72Lib9sHRwJrwQojYknBfhNoHR3F7dEx77qU5qRRnpYQ9YqbZ7B8pIz13IeJBwn0RiuVImWAbKnPDXobAPyZehkEKER8S7otQLFaDDGVjZQ4tZhsjdueM1x7vtZJkVFHZ21UIMTMJ90WoxWylOCuFrNToLRgWyoYl3rr7G2GUZk70WqgpzCTJKB8xIeJB1mGNsCcOnuJLf2kgLclIVqqJqoIM/vtdmyI6HLDZbI15rx1g89JckoyKF4+buWRF4VmvPd5nCcxsFULEnnSrIkhrzfd3NpNiMrBlaS6ZKSb+8kY3r7UMROw13B7N8V4rtcWxv1GZnZrExcsL+euRHrTW0143Ou6iY3CMlXFooxDCS8I9gg60D3Gs18J9V67gO3ds5ucf2Ep6spFnj/ZE7DUaukewOlzUVeVF7Dln44b1pbQPjtLYM/0iYv49V+VmqhDxI+EeQb/a1U5miombN5YDkJpk5PKVRTx7tBePZ/qe7mzs8v0WsLU6PyLPN1vXri1BKfjrkel/YB3v9YZ7rQyDFCJuJNwj5PToOH9+o5vbN1eQkXLmVsZ160roHXFw+NTcNruYbE/rIEvz0ynLif6aMqEUZqZwwbJ8ttdPH+4nei0kGw1UFaTHsGVCiGAS7hHy+P5Oxl0e7tq2dMLxK1cVYzQodpwlDMPl8Wj2tA2yLU69dr/r15fS2GOhtd8W8vyJPu8NX5OMlBEibuRfXwRorfnfPe1sXprLmrLsCedy05PZVp3PjqO9836dE31WTo8641aS8bt+XQlAyN67x6Np6B6RkowQcSbhHgG7WwdpMdt497ZlIc9ft7aEpj5rYIONub+Ot95+YU3BvJ5nvirz0jmvIidkuG+v76F72M5Vq4vi0DIhhJ+E+zxprXlwZxM5aUm8ZUNZyGuuXVcKwLPz7L3vbh2kLCeVyhis4T6TG9aX8nr7aXqG7YFjbo/m288eZ3lRBrdsrIhj64QQEu7z9LfGPl4+0c8/X11LapIx5DUVuWmsK88OWZrpGBzlrd9/hc88fpi/HunGMs3Ufq01u1u89XalYrOG+9lc7/uB9d3njwfGvD91qIsTfVY+ee2qmK0zL4QI7ZwNd601f2vs5ZWm/lltQBFs3OXhv/58lOVFGbz3otAlGb/r1pZyoH2IE0GbTHs8mk/97hAN3RaePtLNh355gPP/67mQPfzWfhv9Vgdbq+NbkvFbUZzJvZfX8Os9HXx9+zGcbg/fee44a8qyuXF9abybJ8Q5LyGXH7A73fzv7naO+UZ09Frs3LV1Kf94aQ1Gg8LhcvPvf3iDPxw4BUCyycD5S/P42FUruHiGafXBHnm1lbaBUR55/wUzrqHyrguW8ItdbXzg53v540cuoTAzhUdebWN36yBff9sGbt9SwYGTQ3zhyXo+/8c3uLAmf8LaMXtaBwHYVhPfm6nBPnvDaqx2Fw+90MzB9tO0DYzy4/fWYZBeuxBxl3A993GXh4/86gAP/Pkozzf2AVCSncpXn2nkbQ+9yu6WAd79o9384cApPnFNLT973wXcfdEyOoZG+cDP97KvbTCs1zFbHHzv+SauWl3MFauKZ7y+NCeVH999AX0jDu55dB8N3SN8fXsjV64q4h11lSQZDWyrKeCrbz0Ps9XBd547MeHxu1sHKcxMoaZw4ayyqJTiv25dz22bynmtZYCNS3K5es3M74UQIvrU2dYIiaa6ujq9b9++iD6ny+3hY79+nWeO9PDl29cHRq9orXnqcDdfeOIIQ6NOUpMMfOsdm3hz0A3QfquDd/7gNfqtDh6796IpQxqD7T85xANP1VPfNcKOf7mMmlmsq/7MG918+FcHSDYZSEsysuNfLqMkO3XCNZ/74xs8treDp+57E2vLs2nqs3DHw7vZVp3Pg+/eMst3Jfqcbg8Pv9TCtWtLZHMOIaJMKbVfa10343XhhLtS6gbgu4AR+LHW+muTzqcAjwLnAwPAu7TWbWd7zkiHu8Pl5jOPH+ZPB7u4/y1r+cCbqqdcY7Y4+PHfW7h5QznrK3KmnO8cGuXtD72GW2seevcWNlTmBlZzHBt3c6RrmJ+83Mpf63soykrhP96yllt8Sw3Mxg9fbOarzzTy3Ts2ceumqaNKhkedXPWtF1hakM7W6nx+8nIr6clGHn5vXdyHQQoh4iti4a6UMgLHgWuBTmAvcKfW+mjQNR8BNmitP6SUugO4XWv9rrM973zDXWuNy6M5OWDjsb0dPL6/k6FRJ/92wyo+csWKOT/viV4L7/jha5wedZJkVKwqzcLl1pzos+L2aDKSjdx7+XL+8dJq0pPnfstiwOqgIDNl2vOP7+/kU787BMA76yr5zA2rz3q9EOLcEG64h5NOW4EmrXWL74l/A9wKHA265lbgi76vHwf+n1JK6SjUfH7y91a+sb2RcZcH/1pcJoPiunUlvHvbshnXGZ9JbUkWz33ycna3DPLGqWHqu4YxKMW1a0s4ryKHC6ryyctInvffY6agftuWCix2Jxsqczh/2cK5iSqEWBzCCfcKoCPo+05g23TXaK1dSqlhoADoD75IKXUPcI/vW6tS6thcGg0UTn7uh3x/FqApbV2gpJ2RtVjaCYunrdJOr7OPu/YJJ9xDjWub3CMP5xq01g8DD4fxmmdvkFL7wvm1ZCFYLG2VdkbWYmknLJ62SjtnJ5yhkJ3AkqDvK4Gu6a5RSpmAHCC8MYVCCCEiLpxw3wvUKqWqlVLJwB3Ak5OueRK42/f124G/RaPeLoQQIjwzlmV8NfT7gO14h0L+VGtdr5R6ANintX4S+AnwC6VUE94e+x3RbDQRKO3E0GJpq7QzshZLO2HxtFXaOQtxm8QkhBAiehJu+QEhhBAS7kIIkZAWdLgrpW5QSh1TSjUppT4b4nyKUuox3/ndSqmqOLRxiVJqp1KqQSlVr5T6eIhrrlBKDSulDvr+3B/rdga1pU0p9YavHVOmCCuv7/ne08NKqZgvZqOUWhX0Xh1USo0opT4x6Zq4vKdKqZ8qpfqUUkeCjuUrpZ5VSp3w/Tdvmsfe7bvmhFLq7lDXxKCt31BKNfr+3/5RKZU7zWPP+jmJQTu/qJQ6FfT/96ZpHnvWjIhBOx8LamObUurgNI+N2fsZoLVekH/w3rxtBmqAZOAQsHbSNR8BfuD7+g7gsTi0swzY4vs6C+9SDZPbeQXw53i/p762tAGFZzl/E/AM3rkLFwK7F8DnoAdYthDeU+AyYAtwJOjY14HP+r7+LPB/QzwuH2jx/TfP93VeHNp6HWDyff1/Q7U1nM9JDNr5ReBTYXw2zpoR0W7npPPfAu6P9/vp/7OQe+6BZQ+01uOAf9mDYLcCP/d9/ThwtYrxNkVa626t9QHf1xagAe+M3cXqVuBR7bULyFVKhd4/MDauBpq11ifj2IYArfVLTJ3DEfw5/DlwW4iHXg88q7Ue1FoPAc8CN0StoYRuq9Z6h9ba5ft2F955K3E1zXsajnAyImLO1k5f7rwT+HW0Xn+2FnK4h1r2YHJoTlj2APAvexAXvrLQZmB3iNMXKaUOKaWeUUqti2nDJtLADqXUft9yEJOF877H0h1M/w9mobynJVrrbvD+sAdCLWq/0N5XgA/g/S0tlJk+J7Fwn6989NNpSl0L6T29FOjVWp+Y5nzM38+FHO4RW/YgFpRSmcDvgU9orUcmnT6At6ywEfgf4E+xbl+QS7TWW4AbgY8qpS6bdH4hvafJwC3A70KcXkjvaTgWzPsKoJT6POACfjXNJTN9TqLtIWA5sAnoxlvymGwhvad3cvZee8zfz4Uc7otm2QOlVBLeYP+V1voPk89rrUe01lbf108DSUqp+S1fOUda6y7ff/uAP+L91TZYOO97rNwIHNBaT9lUdiG9p0Cvv3Tl+29fiGsWzPvqu5n7FuDd2lcQniyMz0lUaa17tdZurbUH+NE0r78g3lNf9rwVeGy6a+Lxfi7kcF8Uyx74am0/ARq01t+e5ppS/70ApdRWvO/7QOxaGWhHhlIqy/813ptrRyZd9iTwXt+omQuBYX/JIQ6m7Q0tlPfUJ/hzeDfwRIhrtgPXKaXyfCWG63zHYkp5N975DHCL1np0mmvC+ZxE1aT7PLdP8/rhZEQsXAM0aq07Q52M2/sZy7u3s/2Dd+TGcbx3xD/vO/YA3g8mQCreX9mbgD1ATRza+Ca8vwoeBg76/twEfAj4kO+a+4B6vHfzdwEXx+n9rPG14ZCvPf73NLitCnjQ956/AdTFqa3peMM6J+hY3N9TvD9sugEn3p7jB/He53keOOH7b77v2jq8O5f5H/sB32e1CXh/nNrahLdO7f+s+keblQNPn+1zEuN2/sL3+TuMN7DLJrfT9/2UjIhlO33HH/F/LoOujdv76f8jyw8IIUQCWshlGSGEEHMk4S6EEAlIwl0IIRKQhLsQQiQgCXchhEhAEu5CCJGAJNyFECIB/X8CQlOHXIN2RgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11c9610b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(tlist,sim.expect[0])\n",
    "plt.ylim([0,1])"
   ]
  },
  {
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}