{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculation of control fields for QFT gate on two qubits using the CRAB algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alexander Pitchford (agp1@aber.ac.uk)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Example to demonstrate using the CRAB [1][2] algorithm in the control library \n",
    "to determine control pulses using the ctrlpulseoptim.create_pulse_optimizer function to \n",
    "generate an Optimizer object, through which the configuration can be\n",
    "manipulated before running the optmisation algorithm. In this case it is\n",
    "demonstrated by modifying the CRAB pulse parameters to show how pulse constraints\n",
    "for controls can be applied.\n",
    "\n",
    "The system in this example is two qubits in constant fields in x, y and z\n",
    "with a variable independant controls fields in x and y acting on each qubit\n",
    "The target evolution is the QFT gate. The user can experiment with the\n",
    "different:\n",
    "    phase options - phase_option = SU or PSU\n",
    "    propagtor computer type prop_type = DIAG or FRECHET\n",
    "    fidelity measures - fid_type = UNIT or TRACEDIFF\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 guess and ramping pulse parameters can be tried.\n",
    "The initial and final pulses are displayed in a plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import datetime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from qutip import Qobj, identity, sigmax, sigmay, sigmaz, tensor\n",
    "from qutip.qip.algorithms import qft\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",
    "import qutip.control.pulsegen as pulsegen\n",
    "\n",
    "example_name = 'QFT'\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining the physics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Sx = sigmax()\n",
    "Sy = sigmay()\n",
    "Sz = sigmaz()\n",
    "Si = 0.5*identity(2)\n",
    "\n",
    "# Drift Hamiltonian\n",
    "H_d = 0.5*(tensor(Sx, Sx) + tensor(Sy, Sy) + tensor(Sz, Sz))\n",
    "# The (four) control Hamiltonians\n",
    "H_c = [tensor(Sx, Si), tensor(Sy, Si), tensor(Si, Sx), tensor(Si, Sy)]\n",
    "n_ctrls = len(H_c)\n",
    "# start point for the gate evolution\n",
    "U_0 = identity(4)\n",
    "# Target for the gate evolution - Quantum Fourier Transform gate\n",
    "U_targ = qft.qft(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining the time evolution parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Number of time slots\n",
    "n_ts = 200\n",
    "# Time allowed for the evolution\n",
    "evo_time = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set the conditions which will cause the pulse optimisation to terminate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Fidelity error target\n",
    "fid_err_targ = 1e-3\n",
    "# Maximum iterations for the optisation algorithm\n",
    "max_iter = 20000\n",
    "# Maximum (elapsed) time allowed in seconds\n",
    "max_wall_time = 300"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Give an extension for output files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Set to None to suppress output files\n",
    "f_ext = \"{}_n_ts{}.txt\".format(example_name, n_ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create the optimiser objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "optim = cpo.create_pulse_optimizer(H_d, H_c, U_0, U_targ, n_ts, evo_time, \n",
    "                fid_err_targ=fid_err_targ, \n",
    "                max_iter=max_iter, max_wall_time=max_wall_time,\n",
    "                alg='CRAB', \n",
    "                dyn_type='UNIT', \n",
    "                prop_type='DIAG', \n",
    "                fid_type='UNIT', fid_params={'phase_option':'PSU'}, \n",
    "                log_level=log_level, gen_stats=True)\n",
    "                "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Configure the pulses for each of the controls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:qutip.control.pulsegen:The number of CRAB coefficients per basis function has been estimated as 3, which means a total of 6 optimisation variables for this pulse. Based on the dimension (4) of the system\n",
      "INFO:qutip.control.pulsegen:The number of CRAB coefficients per basis function has been estimated as 3, which means a total of 6 optimisation variables for this pulse. Based on the dimension (4) of the system\n",
      "INFO:qutip.control.pulsegen:The number of CRAB coefficients per basis function has been estimated as 3, which means a total of 6 optimisation variables for this pulse. Based on the dimension (4) of the system\n",
      "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"
     ]
    }
   ],
   "source": [
    "dyn = optim.dynamics\n",
    "\n",
    "# Control 1\n",
    "crab_pgen = optim.pulse_generator[0]\n",
    "# Start from a ramped pulse\n",
    "guess_pgen = pulsegen.create_pulse_gen('LIN', dyn=dyn, \n",
    "                                           pulse_params={'scaling':3.0})\n",
    "crab_pgen.guess_pulse = guess_pgen.gen_pulse()\n",
    "crab_pgen.scaling = 0.0\n",
    "# Add some higher frequency components\n",
    "crab_pgen.num_coeffs = 5\n",
    "\n",
    "# Control 2\n",
    "crab_pgen = optim.pulse_generator[1]\n",
    "# Apply a ramping pulse that will force the start and end to zero\n",
    "ramp_pgen = pulsegen.create_pulse_gen('GAUSSIAN_EDGE', dyn=dyn, \n",
    "                                    pulse_params={'decay_time':evo_time/50.0})\n",
    "crab_pgen.ramping_pulse = ramp_pgen.gen_pulse()\n",
    "\n",
    "# Control 3\n",
    "crab_pgen = optim.pulse_generator[2]\n",
    "# Add bounds\n",
    "crab_pgen.scaling = 0.5\n",
    "crab_pgen.lbound = -2.0\n",
    "crab_pgen.ubound = 2.0\n",
    "\n",
    "\n",
    "# Control 4\n",
    "crab_pgen = optim.pulse_generator[3]\n",
    "# Start from a triangular pulse with small signal\n",
    "guess_pgen = pulsegen.PulseGenTriangle(dyn=dyn)\n",
    "guess_pgen.num_waves = 1\n",
    "guess_pgen.scaling = 2.0\n",
    "guess_pgen.offset = 2.0\n",
    "crab_pgen.guess_pulse = guess_pgen.gen_pulse()\n",
    "crab_pgen.scaling = 0.1\n",
    "\n",
    "init_amps = np.zeros([n_ts, n_ctrls])\n",
    "for j in range(dyn.num_ctrls):\n",
    "    pgen = optim.pulse_generator[j]\n",
    "    pgen.init_pulse()\n",
    "    init_amps[:, j] = pgen.gen_pulse()\n",
    "\n",
    "dyn.initialize_controls(init_amps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the pulse optimisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial amplitudes output to file: ctrl_amps_initial_QFT_n_ts200.txt\n",
      "***********************************\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:qutip.control.optimizer:Optimising pulse(s) using CRAB with 'fmin' (Nelder-Mead) method\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting pulse optimisation\n",
      "Final amplitudes output to file: ctrl_amps_final_QFT_n_ts200.txt\n"
     ]
    }
   ],
   "source": [
    "# Save initial amplitudes to a text file\n",
    "if f_ext is not None:\n",
    "    pulsefile = \"ctrl_amps_initial_\" + f_ext\n",
    "    dyn.save_amps(pulsefile)\n",
    "    print(\"Initial amplitudes output to file: \" + pulsefile)\n",
    "\n",
    "print(\"***********************************\")\n",
    "print(\"Starting pulse optimisation\")\n",
    "result = optim.run_optimization()\n",
    "\n",
    "# Save final amplitudes to a text file\n",
    "if f_ext is not None:\n",
    "    pulsefile = \"ctrl_amps_final_\" + f_ext\n",
    "    dyn.save_amps(pulsefile)\n",
    "    print(\"Final amplitudes output to file: \" + pulsefile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Report the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "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:05:00.003775\n",
      "Wall time computing Hamiltonians: 0:00:22.547817 (7.52%)\n",
      "Wall time computing propagators: 0:04:24.789280 (88.26%)\n",
      "Wall time computing forward propagation: 0:00:03.857988 (1.29%)\n",
      "Wall time computing onward propagation: 0:00:03.634444 (1.21%)\n",
      "Wall time computing gradient: 0:00:00 (0.00%)\n",
      "\n",
      "**** Iterations and function calls ****\n",
      "Number of iterations: 8610\n",
      "Number of fidelity function calls: 10459\n",
      "Number of times fidelity is computed: 10459\n",
      "Number of gradient function calls: 0\n",
      "Number of times gradients are computed: 0\n",
      "Number of times timeslot evolution is recomputed: 10459\n",
      "\n",
      "**** Control amplitudes ****\n",
      "Number of control amplitude updates: 10458\n",
      "Mean number of updates per iteration: 1.2146341463414634\n",
      "Number of timeslot values changed: 2091599\n",
      "Mean number of timeslot changes per update: 199.99990437942245\n",
      "Number of amplitude values changed: 8350196\n",
      "Mean number of amplitude changes per update: 798.4505641614076\n",
      "------------------------------------\n",
      "Final evolution\n",
      "Quantum object: dims = [[4], [4]], shape = (4, 4), type = oper, isherm = False\n",
      "Qobj data =\n",
      "[[-0.47890815-0.25124383j -0.47862720-0.14293885j -0.47089104-0.19654053j\n",
      "  -0.39814103-0.19780087j]\n",
      " [-0.48524763-0.18754841j  0.19902638-0.45262542j  0.45666538+0.17490099j\n",
      "  -0.22150294+0.44348832j]\n",
      " [-0.44119050-0.18412439j  0.49542987+0.17619236j -0.42828736-0.17952947j\n",
      "   0.50360656+0.16023168j]\n",
      " [-0.41755379-0.18434167j -0.18861574+0.44037804j  0.50645054+0.16836502j\n",
      "   0.23840544-0.4695553j ]]\n",
      "\n",
      "********* Summary *****************\n",
      "Initial fidelity error 0.9483036893449602\n",
      "Final fidelity error 0.0032564547889728512\n",
      "Terminated due to Max wall time exceeded\n",
      "Number of iterations 8611\n",
      "Completed in 0:05:00.003775 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(\"Initial fidelity error {}\".format(result.initial_fid_err))\n",
    "print(\"Final fidelity error {}\".format(result.fid_err))\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": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztvXnYHFWZ9//5Zt9DSDCBkIU9BMOigaAyCBNwIAQYkWUE\nGcg7DuqrII6DIoxsMwgor4LgbwCRiCAICTBgDIsgi6JAAjKEEAKyBAJZSAhZSchy//6o6qSeTi/V\n3dVd1d3357qe6+mu5Zy7qrvPt+773OccmRmO4ziOkzU6pW2A4ziO4xTCBcpxHMfJJC5QjuM4TiZx\ngXIcx3EyiQuU4ziOk0lcoBzHcZxM4gLltAWS7pd0Won910n6fsyyHpP05eSsaxySDpE0P207HCcO\nLlBO0yLpTUmHxTnWzI40s5vD806X9Ke8/V81s/9MyK7dJU2RtETSckkvSPo3SZ1rLNfFxWkrXKAc\nJ0Ek7QI8DbwNjDGz/sAJwCeBvg2ov0u963CcRuEC5bQEOa9I0pWSlkl6Q9KRkf2PSfqypD2B64BP\nSVol6YNw/y8l/Vf4eoCkaZLeC8uaJmnHmKZcDPzZzP7NzBYAmNlcMzvFzHJ1HSNptqQPQrv2jNj5\npqR/D72u5ZLukNRDUm/gfmCH0O5VknaQdJGkqZJulbQCOF1Sd0lXSXo3/LtKUveY9/FqSW9LWiHp\nWUl/F9l3UegZ3ipppaRZobf4PUmLw/M+l3fPL5P0TFjevZK2Dff1CMtZGt6HGZIGx7zHTpvgAuW0\nEuOAucAg4IfALyQpeoCZzQG+CvzFzPqY2TYFyukETAZGAMOBD4FrY9pwGDC12E5JuwO3A2cD2wHT\ngd9K6hY57ETgCGAnYG/gdDNbDRwJvBva3cfM3g2PPzascxvg18D5wIHAvsA+wAHAf8S0f0Z43rbA\nbcAUST0i+48GbgEGAH8FHiS4X0OBS4Dr88r7Z+D/ANsDG4CfhttPA/oDw4CBBJ/JhzFtdNoEFyin\nlZhnZj83s43AzQSNYsVP5Wa21MzuMrM1ZrYSuBT4bMzTBwILSuw/Cfidmf3ezNYDVwI9gU9Hjvmp\nmb1rZu8DvyUQjFL8xcz+x8w2mdmHwCnAJWa22MzeI/DqTo1jvJndGl7/BjP7f0B3YI/IIX80swfN\nbAMwhUBkLw+v5TfASElR0b/FzF4MBfb7wIlhX9x6gnu1q5ltNLNnzWxFHBud9sEFymklFuZemNma\n8GWfSguR1EvS9ZLmhWGzJ4BtYiY5LCUQxmLsAMyL2LmJoL9qaOSYhZHXayh/DW+XqiN8vUOZMgAI\nw4tzwvDiBwRezqDIIYsirz8EloQPBLn35NkbtW0e0DUs7xYC7+s3YRjyh5K6xrHRaR9coJx2pNwU\n/t8m8BrGmVk/4OBwu4qfspmHgS+U2P8uQegwKDAIQQ4D3olRdjG787d3qIMgTPkuZQj7m75DEGIc\nEIY/lxPvuosxLM+O9QSitt7MLjaz0QTe40SCcKDjbMYFymlHFgE75vX7ROlL4A18EHbqX1hB2RcC\nn5b0I0lDACTtGiYEbAPcCRwlaXzoMXwbWAf8OabdAyX1L3Pc7cB/SNpO0iDgAuDWGOX3Jegneg/o\nIukCoF+M80rxJUmjJfUi6KOaamYbJR0qaUzola4gEK5NNdbltBguUE478gdgNrBQ0pIC+68i6Bda\nAjwFPBC3YDN7DfgUMBKYLWk5cBcwE1hpZnOBLwHXhOUfDRxtZh/FKPtlAvF5Pcx8Kxa2+6+wvheA\nWcBz4bZyPEhwra8QhOPWsnX4sFJuAX5JELbsAZwVbh9CkNixApgDPB4e6zibkS9Y6DhOPZD0GHCr\nmd2Yti1Oc+IelOM4jpNJYgmUpJ6S9ih/pOM4juMkQ9kQn6SjCcZqdDOznSTtSzDG4phGGOg4juO0\nJ3E8qIsIRqJ/AGBmzxOMcHccx3GcuhFnYsn1ZrY8f8aYOtmTCIMGDbKRI0embYbjOI5TgGeffXaJ\nmW1X7rg4AjVb0slAZ0m7EaSJxhmzkRojR45k5syZaZvhOI7jFEDSvPJHxQvxnQnsRTCY8HaCcQtn\nV2/aFsKZm2dJel7SVoqigJ9K+ls4u/MnkqjXcRzHyT5lPahwTrPzw796cKiZFRosCcHszbuFf+OA\n/w7/Oxln2R13smLaNAD6TZzIgJNOTNkix3GajaICJem3lOhralAW37HAryxINXxK0jaSts+ts+Nk\ni6gorZkxY/P2NTNmuFg5jlMxpTyoK8P/xxFMS5Kby+uLdJzRuBYMeFjSRuB6M7shb/9QOk61Mj/c\ntpVASToDOANg+PDhCZnnxCEnTDlR6rX//vTaf3/6TZwI0EG0XKwcx4lLUYEys8cBJP0/Mxsb2fXb\nQv1FVXKQmb0j6WPA7yW9bGZPVFNQKG43AIwdOzbTWYatQCFvKSdK+aKTe59/Tk6sXKgcxylEnCy+\n3pJ2NrPXASTtBPROonIzeyf8v1jSPQTjraIC9Q4dp+vfkXjLEjh1opS3VE5kBpx04lZi5V6V4zjF\niCNQ3wIek/Q6wbowI4Cv1FqxpN5AJzNbGb7+HMF0/FHuA74h6TcEyRHLvf8pPZbdcScLLwxWnogr\nSsXIiZV7VY7jFCPWbOaSugOjwrcvm9m6miuWdgbuCd92AW4zs0slfRXAzK4LF3O7FjiCYGXRSWZW\nNrw4duxY83FQyZHvNQ25+OK6iEch78yFynFaD0nP5nUdFT4uxlx8BVe5NLNfVWlb3XGBSoa0BKNR\ngug4TjokKVDXRN72AMYDz5nZ8bWZWD9coGonyXBeM9vgOE7yxBWoOAN1z8wreBvgNzXY5mSYLHkv\nuXrzkylcqBynPYiTJJHPanw285Yjq/0/+ckUOaHK7XMcp3WJE+KLzijRCRgNTDGz79bZtqrxEF9l\nNFMorZlsdRynMImF+NgyowTABmCemc2v2jInM2QpnBeXQmG/6HbHcVqHOB7UFfneUqFtWcI9qPK0\ngifSCtfgOO1Ikh7U4UC+GB1ZYJvTBDSj11SMYkkU4DNSOE4rUNSDkvQ14P8COwOvRXb1BZ40sy/V\n37zqcA9qa7KaBJEUlcwN6DhOutQ8DkpSf2AAcBlwbmTXSjN7PxEr64QLVEfaLRTW6mLsOM1OEiE+\nM7M3JX29QOHbZl2knNYK51WCp6Y7TmtQyoOaZmYTJb1BkGauyG4zs50bYWA1tLsH5R5ER9rNg3Sc\nrFOzB2VmE8P/Pii3ifDGeGs8mcJxaue2p9/i3ueD1Y5G79CPC4/eq+51llry/ROlTjSz55I3x6mU\naHIA0HbhvLj48h6OU56oCOXz9BtBr864nbZtmD2lQnyPljjPzOzv62NS7bRLiC/fW8rhjW08CoVC\nwe+f0/oUE6JyInTsvkM5edzwmutPbDbzeiFpGPArYDBBH9cNZnZ13jGHAPcCb4Sb7jaz/EUNt6LV\nBapdkx/qRbEUdXCxcpqfQmJUSoiSEqFSJLncRg+C8VAHEQjJH4HrzGxtjQZuD2xvZs9J6gs8C/yj\nmb0UOeYQ4N9z/WFxaUWB8nE+jaGUWIELlpNNqgnNNUKIipGkQN0JrARuDTedDGxjZifUbGXHeu4F\nrjWz30e2HUIbC5Q/2adLsf69qGCBfxZO40g7NJcUSQrUS2Y2uty2WpA0EngC+LiZrYhsPwS4G5gP\nvEMgVrOLlHEGcAbA8OHDPzlv3rykzGsY5RpEbwjTJf/zgeKiVQj//JxKyGJoLimSFKhbCTybp8L3\n44Cvm1nBpeCrMLQP8DhwqZndnbevH7DJzFZJmgBcbWa7lSuzGTyouI2dN2rZptDnWIhKhCyKf/6t\nTbOF5pIiSYGaA+wBvBVuGg7MJVh6w8xs7xqM7ApMAx40sx/HOP5NYKyZLSl1XNYEqpInb2+QWpO4\nQhYljqj596U5aJXQXFIkKVAjSu03s6piaZIE3Ay8b2ZnFzlmCLDIzEzSAcBUYISVMbqRAhWn4XEx\ncqqh3HfLv1fZpJVDc0mRaJq5pAHAMCIDe2sdqCvpIIKMwFnApnDzeQQeGmZ2naRvAF8j8NY+BP7N\nzP5cruykBKoW8cnHGw0nadwzT4+6huZmToZZUys3aszxMHZS5eelQJIe1H8CpxMsuZE7uGUH6i78\nwQ9YN+dlwMXHaT7iiJZ/X+NT99BcITGa96fg/4iD4hta6JwMC1aSAjUXGGNmHyVlXL1JSqDAf8xO\n81NubBf49xxSCM3NnAzTwt6NfDGqVFzyhS4nWBOvqlmkCj30dN9zFEPOO6/qMpMUqLuAr5nZ4qqt\naTBZS5JwnKzQ7h5WqllzdRSRgnVFxa9CwSv3UJMlgRpLMN3Qi8C63HYzO6Zq6+pM5gSqXEw5w664\n0/q04uwZmcuaK+Yt1fO3n2t3ouG/EvU1cmKAJAVqNnA9HZMZMLPHazWyXmRCoKKiVCqmnL/PxcpJ\nkWacHT+TWXON9Jbi2hKxYdlrvVP1pJMUqBlmVtnowpRJVaAKPbVAceEpJ2QuWE6KFJoxPw2PqqkG\ntKbhLRVgq3DuyoWw9G/Qoz9r3g6mUk2rLzJJgfoxQWjvPjqG+DK7HlRqAlVj3LfoU5d7V06KNCoE\nWMtaRJkYS1TAU2nE77XYcJiCCTELZ8FHq6Fbb/odfjADvnP1Vuc1giQFqtC6UC2bZl4Tk48KvpxJ\nfTELeVdphAjalCmvTGH669OrOnfCzhM4YfdE51POBEnOF5kvSE0hQsWo9eG0DKXGZJYaDrPV51Fh\nv1S9yPx6UPWk4QKV+9AXzoIhY2DS7+pTRwbCBq1OVJRmLgq+Q2MHl/0ddSD/vFYVKyjfsX7/yAMr\n8ooyLUKFSNhrqsgbilCxF9uINqsESc8kcRSwF9Ajty3OwoFp0VCBqvOT01Z1ZaXjtYUoJUrViEux\n8lpVqHKe0CdfeIwxc58GYOT8uQC8MHBnAPr17ArArD3G8ezeh2w+t+kEKUcNnkhi3lASJB31iUmS\nIb7rgF7AocCNwPHAM2b2L0kYWg8aKlApfcCAe1U1khOSeno8+XVc8KkLmlqk4mbMRcVqYJ/uDO7b\nvWVS2Ct5KK1miZaG3pPItUwZsS/T+/SCvkPKnjZq21F894DvVl1tkgL1gpntHfnfB7jfzP6uauvq\nTEMEKmUXeSs7cnhiRSymvDKFS/4SBAEa4d00ur5aqXQcURxPqGnXO4sZuWimuRE79K+uXAir32Om\nghy4OCHtLAnU02Y2TtJTwHHAUmC2me1atXV1pq4ClZFOxqIUS1vPko0pkqZHkzVvqtqsuaTCck0x\nBVORKMWy93ZlxaxlHQ7N7DVEKBY1AGDhLCYsfZcTPntp3duKJAXq+8A1wHjgZwQTxv7czC5IwtB6\nUFeBmnzUFq8p641+pWOyWphCP8y0vJg0vKlKB7NCY/uHMud9RH47y/7WixWrxnQIfTWDGOVT9nvX\nwP70umTxSeoO9DCz5bUYFynvCOBqoDNwo5ldnrdf4f4JwBrg9Djjr+omULkPcMRB6YX0qiHOrBYt\nLFpZDK/VSzDrEZpLizREa9kPv8mK3z8Ba8MmLgODWmulIs+9QRGizKeZS+oMvAIcDswHZgBfNLOX\nIsdMAM4kEKhxBEu+jytXduICldIAvGopOX4njDd3IPwxTui+AyfQZ8v2JhetrIXUClGNjWmH5tKk\n0qSDUiJSaKaFNS+/HZQ1rAf03m6z19QsYhSlpoegOrd5zSBQnwIuMrN/CN9/D8DMLosccz3wmJnd\nHr6fCxxiZgtKlV2TQN1/bhDCi9IEfTk1jd9ZuZCZa+YH51j3YNva5UxYvZoTtt2v/PkZuydZCufF\nJefljR08lslHTN68PeuhuSxQ7dihzftHDQse2sIHtX5HH53aDAtJkVjUoFim8JAxcOTlhc+JQVyB\n6lLugDoyFHg78n4+gZdU7pihwFYCJekM4AyA4cMT/nE2gTBFG+NqvpD5XtfMRTOZ2bMH020RABOs\nd0fvKse8PwV/1awAmiPBe5vFcF4c1i8bR69Nu/PcwtmMm/wF+m88gAEbDy4oRuN22rbtRKgUA046\nsaB3U3JF7JUL6TWsB/1G92HAgCAdPsu/87gkHjXI3Ytaft81UNSDkvSJUicmsOT78cARZvbl8P2p\nwDgz+0bkmGnA5Wb2p/D9I8B3zayke5SJ2czrSCMGgsYevFrt8tQ5yq0eGrPBaIZwXrnQXNdtnmbb\nwbNZ0+kVAHpt2p3+Gw/g9DFfdDFKihZNHGq2qEESHtT/K7HPgFrn4nsHGBZ5v2O4rdJj2oakvKU4\nnLD7CZvLzBermYtmMv316UG9YyfV9sMuJXDFvLNIY5LVH2alobnAK/oyJ48b3uGa1nR6ha4DdgZc\noKqmDYZeTH99OnPfn5uZ739SpNkH1YUgSWI8gejMAE42s9mRY44CvsGWJImfmtkB5cpuRQ8qK6Gr\nRsy+sJlC4jXvT0zp25vpA3cIDokMLkzjntQza65Yv5QTkxb1lvJpxu9JkuOgugJfAw4ONz0GXG9m\n6xMwcgJwFUGa+U1mdqmkrwKY2XVhmvm1wBEEaeaTyoX3oLUEKquhq7Tmm5vy0Le4ZMHDQX0fBum/\nE3aeyAmf+0nd6kwza27SA5OY+/5c9th2j5Z6Mq4bbeAt5chq2xCHJAXqRqArcHO46VRgY67vKIu0\nikBlxWsqR8nR6dTmXRVK3oDwx7hiVeJzEWYtay6rIcxMUMTDBlrWW4LW+E4kKVD/a2b7lNuWJZpd\noJr1yajQ+KtKRCvu+SWTNGKO22i2Aa3N+p2omXJ9lOCDzpvwe5CkQD0HnGBmr4XvdwammlnJLL80\naVaBaoUno3wqEa1EPLDIuI1F245lyaqgj+rJnofySK8Jmw9r1gGtTd1A5YnNFFYxXatLnxOZ0aEg\n4WDaproPNTDpgUnMXDSz6R9QkhSo8cBk4HVAwAiCvqBCK+1mgmYUqKZueCqk1EwX1V531CMav2Y6\nn/nwUVau3QDAgZ3mAHBD/7M6iFSWhagUmXuQiTnUYMr7f2V6796bxWZzgktucHgxIjM6FKy+DRaI\nzH3muf7IZkmGKEYiAiWpE3Ag8CywR7h5rpmtS8TKOtEMAlWyb6XFflxJUlVorvMjjVtUsoGkEvaL\n0+8TtTHiJRVaziEJMUl6wcms0YoPr0l6UH81sxjz3WSHrAtU/hcuRyt88ZKgLllzWV8mpQYKfZ8S\n+y6VW28sx5jjmdKvT9lwbiPW3So15Vcz/Mba4eE1SYG6EvgLcLelNWiqQrImUO3whauWhmfNtejY\nmMQa5jiCVGCgdLF6K6o7YYr97rLsXbXLw2uSArUS6A1sANYS9EOZmfVLwtB6kAWBaoUnuaSoxiOq\ne/9QC4+XKdcw59j8HSyzHMuUobszfdMHBetqhkY/RxZDge368Jr52czrSS0CdcUzV/Dy+y/XbENW\nfgCNpNL+oRypJys02XIqlVIyk9K6b50pl5eUUG52/Gb8bqfp+fnDa7Ie1CNmNr7ctiyRBYGC1v6i\nVRqaS12E4tDAFUVTIxTj/Iy6cplyrfxdrmQoRBK048NrPjULlKQeQC/gUeAQgtAeQD/gATMblYyp\nyZOFEF8rkMnQXL1pRW+qhcOZ9aLkop8J0I6iFCUJgfomcDawA8FkrjmBWgH83MyuTcjWxHGBqoym\nDc3Vk2ILtTVTo96iCSFO81PzchtmdjVwtaQzzeyaRK1zUqOS0FxbL4xXaKG23PIf0f1Zw70lp4WI\nlSQh6dPASCKCZma/qp9ZtdHuHlRbhuYaQVb7qEpl4WXFRseJkNiS75JuAXYBngc2hpsNyKxAtQvV\nhOba2iuqlahXlb+YYiOFoNRYpSwJp+PUSJwsvjnA6CQH6Ur6EXA08BHwGsHcflsNtJD0JrCSQBg3\nxFFcaE0PqiWz5pqZMmOHEheJRtfnOHUkyTTzKcBZZrYgQeM+B/zBzDZIugLAzL5b4Lg3gbFmtqSS\n8ptVoDw016RUMB1QLBFp03WOnPYhsRAfMAh4SdIzwOZJYs3smGqNM7OHIm+fAo6vtqxmxENzLcbY\nSR3FopjAREOCpSgkcB66c9qQOB7UZwttN7PHEzFA+i1wh5ndWmDfG8ByghDf9WZ2Q4lyzgDOABg+\nfPgn582bl4R5NeOhOQeIvSTFZlyMnBYm0amOJA0G9g/fPmNmi2Oc8zBQaGj6+WZ2b3jM+cBY4LhC\nfVyShprZO5I+BvweONPMnihXd6NDfB6acxzHiU+SWXwnAj8CHiMYrHuNpHPMrOTjoJkdVqbc04GJ\nwPhiCRhm9k74f7Gke4ADgLICVS88NOc4jtM44vRBnQ/sn/OaJG0HPAxUEK/oiKQjgO8AnzWzNUWO\n6Q10MrOV4evPAZdUW2el+IBWx3GcdIkjUJ3yQnpLgU411nst0B34vSSAp8zsq5J2AG40swnAYOCe\ncH8X4DYze6DGesty8W9n89K7KwqKkQuR4zhO44gjUA9IehC4PXx/EnB/LZWa2a5Ftr8LTAhfvw7s\nU0s9teBi5DiOky5lBcrMzpF0HJDLeb3BzO6pr1npceHRe6VtguM4jkMJgZK0KzDYzJ40s7uBu8Pt\nB0naxcxea5SRjuM4TvtRarmNacD3zGxW3vYxwA/M7OgG2FcVkt4DahkINQioaPaKFsXvg9+DHH4f\nAvw+BNR6H0aY2XblDioV4hucL04AZjZL0sgaDKs7cS68FJJmxp33r5Xx++D3IIffhwC/DwGNug+l\nsvG2KbGvZ9KGOI7jOE6UUgI1U9K/5m+U9GXg2fqZ5DiO4zilQ3xnE4xDOoUtgjQW6AZ8vt6GpUzR\nOf/aDL8Pfg9y+H0I8PsQ0JD7EGey2EOBj4dvZ5vZH+puleM4jtP2xJos1nEcx3EaTa1TFjmO4zhO\nXXCBiiDpCElzJf1N0rlp25MGkoZJelTSS5JmS/pm2jaliaTOkv4ajgtsSyRtI2mqpJclzZH0qbRt\najSSvhX+Hl6UdLukHmnb1Agk3SRpsaQXI9u2lfR7Sa+G/wfUq34XqBBJnYGfAUcCo4EvShqdrlWp\nsAH4tpmNBg4Evt6m9yHHN4E5aRuRMlcDD5jZKIL5MdvqfkgaCpwFjDWzjwOdgX9K16qG8UvgiLxt\n5wKPmNluwCPh+7rgArWFA4C/mdnrZvYR8Bvg2JRtajhmtsDMngtfryRojIama1U6SNoROAq4MW1b\n0kJSf+Bg4BcAZvaRmX2QrlWp0AXoKakL0At4N2V7GkK4QOz7eZuPBW4OX98M/GO96neB2sJQ4O3I\n+/m0acOcI5wxZD/g6XQtSY2rCNYt25S2ISmyE/AeMDkMdd4Yrs/WNoQLp14JvAUsAJab2UPpWpUq\ng81sQfh6IcHSSHXBBcopiKQ+wF3A2Wa2Im17Go2kicBiM2v3QeldgE8A/21m+wGrqWNIJ4uEfSzH\nEoj1DkBvSV9K16psEK6GXrdUcBeoLbwDDIu83zHc1nZI6kogTr8OZ7JvRz4DHCPpTYJw799LujVd\nk1JhPjDfzHJe9FQCwWonDgPeMLP3zGw9wcoOn07ZpjRZJGl7gPD/4jLHV40L1BZmALtJ2klSN4JO\n0PtStqnhKFjC+BfAHDP7cdr2pIWZfc/MdjSzkQTfhT+YWds9NZvZQuBtSXuEm8YDL6VoUhq8BRwo\nqVf4+xhPmyWK5HEfcFr4+jTg3npVFGdF3bbAzDZI+gbwIEGWzk1mNjtls9LgM8CpwCxJz4fbzjOz\n6Sna5KTLmcCvwwe314FJKdvTUMzsaUlTgecIslz/SptMeSTpduAQYJCk+cCFwOXAnZL+hWBZoxPr\nVr/PJOE4juNkEQ/xOY7jOJnEBcpxHMfJJC5QjuM4TiZxgXIcx3EyiQuU4ziOk0k8zdxxGoykgQST\nbAIMATYSTCcEsMbM2nkQqONsxtPMHSdFJF0ErDKzK9O2xXGyhof4HCdDSFoV/j9E0uOS7pX0uqTL\nJZ0i6RlJsyTtEh63naS7JM0I/z6T7hU4TnK4QDlOdtkH+CqwJ8HsHrub2QEEy3+cGR5zNfATM9sf\n+AJtvDSI03p4H5TjZJcZuWUNJL0G5JZ4mAUcGr4+DBgdTBEHQD9JfcxsVUMtdZw64ALlONllXeT1\npsj7TWz57XYCDjSztY00zHEagYf4HKe5eYgt4T4k7ZuiLY6TKC5QjtPcnAWMlfSCpJcI+qwcpyXw\nNHPHcRwnk7gH5TiO42QSFyjHcRwnk7hAOY7jOJnEBcpxHMfJJC5QjuM4TiZxgXIcx3EyiQuU4ziO\nk0lcoBzHcZxM4gLlOI7jZBIXKMdxHCeTuEA5juM4mcQFynEcx8kkLlBO5pA0XNIqSZ2rPH+VpJ0T\ntumXkv4ryTIbhaSRkkxS6uu/hUvZz4+8ny3pkATLf1PSYUmV56SLC5RTM5JOlzRL0hpJCyX9t6Rt\nKji/Q6NiZm+ZWR8z21iNPeG5r1dzbrVI2l7SLyQtkLRS0suSLpbUu8ZyGyYukh6TtExS93rXlcPM\n9jKzx8L6L5J0a6PqdrKPC5RTE5K+DVwBnAP0Bw4ERgC/l9QtTdsahaRtgb8APYFPmVlf4HCC+7FL\nA+qvWbwkjQT+DjDgmFrLc5wkcIFyqkZSP+Bi4Ewze8DM1pvZm8CJwEjgS+FxF0maKumO0Lt4TtI+\n4b5bgOHAb8PQ3HfyvYbwyf6/JP05POa3kgZK+rWkFZJmhA1szi6TtGv4eoKkl8J635H075HjJkp6\nXtIHYdl7R/btF9q5UtIdQI8St+LfgJXAl8Lrx8zeNrOzzeyFsLxPh3YuD/9/OlLXY5L+U9KTYX0P\nSRoU7n4i/P9BeO2fCj3WJyX9RNJS4CJJnST9h6R5khZL+pWk/hV8nP8MPAX8EjgtuiMMb/5/ku4P\nbXhS0hBJV4Ue18uS9osc/6ak74X3fZmkyZIK3r+c9yzpCOA84KSwjv+N7o8c38HLknRqeM1LJZ2f\nV3YnSedKei3cf2f4MIGkHpJuDbd/EH4mgyu4X04DcIFyauHTBA333dGNZrYKmE7gReQ4FpgCbAvc\nBvyPpK5mdirwFnB0GJr7YZG6/gk4FRhK4JX8BZgcljcHuLDIeb8AvhJ6NR8H/gCBAAE3AV8BBgLX\nA/dJ6h5trs1xAAAgAElEQVR6fv8D3BKWPwX4Qon7cBhwt5ltKrQzbBR/B/w0rOvHwO8kDYwcdjIw\nCfgY0A3ICenB4f9twvvzl/D9OOB1YDBwKXB6+HcosDPQB7i2hM35/DPw6/DvHwo01icC/wEMAtYR\n3P/nwvdTw2uKcgrwDwSf1e7huUUxsweAHwB3hNe5TzmDJY0G/pvge7EDwb3dMXLImcA/Ap8N9y8D\nfhbuO43Awx0WnvdV4MNydTqNxQXKqYVBwBIz21Bg34Jwf45nzWyqma0naMx6EIQD4zLZzF4zs+XA\n/cBrZvZwWPcUYL8i560HRkvqZ2bLzOy5cPsZwPVm9rSZbTSzmwka3gPDv67AVaFXOBWYUcK2geH1\nFuMo4FUzu8XMNpjZ7cDLwNF51/eKmX0I3AnsW6I8gHfN7JqwvA8JBOHHZvZ6+IDwPeCf4oT/JB1E\nEJa908yeBV4jEMwo95jZs2a2FrgHWGtmvwr7Ce9g6/t/behFvk8goF8sZ0cVHA9MM7MnzGwd8H0g\n+pDwVeB8M5sf7r8IOD68J+sJPrddw8//WTNbUQcbnRpwgXJqYQkwqEgjuH24P8fbuRehpzGf4Kk2\nLosirz8s8L5PkfO+AEwA5kl6XNKnwu0jgG+H4Z0PJH1A8DS9Q/j3jplZpJx5JWxbSnC9xdihwPnz\nCLzBHAsjr9eUuJ4cb+e9z69jHtCFwMMqx2nAQ2aW+7xuIy/MR+X3P2rfPCr7rOOyAx2/V6sJPosc\nI4B7Ip/vHGAjwT25BXgQ+I2kdyX9UFLXOtjo1IALlFMLfyHwOo6LbpTUBzgSeCSyeVhkfyeCUMy7\n4aaoECSKmc0ws2MJQmf/Q+CdQNCwXWpm20T+eoXezQJgqCRFihpeopqHgc+H11WIdwkayyjDgXfi\nXELM7fl1DAc20FFItkJST4Lw3WcVZGAuBL4F7KOwn7BKhkVeD2fLZ12KQte6GugVeT8k8noBHb9X\nvQi8ohxvA0fmfcY9zOyd0DO+2MxGE4SqJxKEOZ0M4QLlVE0YbrsYuEbSEZK6hskKdxJ4SLdEDv+k\npONCb+tsAmF7Kty3iKDfJFEkdZN0iqT+YWhxBVtCQD8HvippnAJ6SzpKUl8C4d0AnBVe03HAASWq\n+jHQD7hZ0oiw7qGSfhwmXkwHdpd0sqQukk4CRgPTYlzGe6HN5e7P7cC3JO0UPiDk+nMKhV+j/COB\nVzGaIKy4L7An8Edqa7C/LmnHsP/tfIIwYDkWASPzhP55glBlV0ljCcJ6OaYCEyUdFPYbXkLHNu06\n4NLIZ7KdpGPD14dKGqNgrN0KgpBfwT5EJz1coJyaCJMazgOuJPihP03w5Do+jPvnuBc4iaCj+lTg\nuFA0AC4D/iMMxfw7yXIq8KakFQR9EqeEds8E/pUgkWAZ8DeCJAPM7CMCr/B04P3Q7rspQtjP8mmC\nRu5pSSsJvMflwN/MbCnBE/q3CUJQ3wEmRkJqRTGzNQR9OE+G96dYv91NBA8ETwBvAGsJkgTKcRpB\n/9dbZrYw90dwX06J04dVhNuAhwgSOV4D4gxynhL+Xyop11f4fYJEi2UED0O35Q42s9nA18NtC8Jj\nNg8CBq4G7gMeCj+TpwiSSyDwxKYSfGfnAI/T8YHKyQDqGGZ3nOSRdBFBZ/SX0rbFqT+S3gS+bGYP\np22L09y4B+U4juNkkkQFSlJPSXskWabjOI7TniQW4pN0NEE/RDcz20nSvsAlZlbztClhyGAlQWfu\nBjMbW2uZjuM4TrZJcgLKiwgynR4DMLPnJe2UYPmHxulUdhzHcVqDJAVqvZkt7zh0pH7jW0oxaNAg\nGzlyZBpVO47jOGV49tlnl5jZduWOS1KgZks6GegsaTfgLODPCZVtwMOSNhJMT3ND/gGSziCYvobh\nw4czc+bMhKp2HMdxkkRSqZlZNpNkksSZwF4EAzBvJxhfcHZCZR9kZvsSzE7wdUkH5x9gZjeY2Vgz\nG7vddmWF2XEcx8k4iXlQ4YDC88O/RDGzd8L/iyXdQ9DX9UTpszLEzMkwa2rx/WOOh7GTGmeP4zhO\nE5DEQme/pURfU61ZfApWJO1kZivD158jmNKkOZg5GaaFjuSIg7beP+9Pwd+sqS5UjuM4EZLwoK4M\n/x9HMH1IbjGxL1JmosqYDCaYkRgCe28L145pDnKe08SrCotPzrtaOCt47wLlOI4DJDsOamb++KRC\n2xrB2LFjLfUkiajwDBkDk35X+vjJRwWeVDEhcxzHaREkPRtHG5JMkugtafOMy+EYqN4Jlt9cRMVp\nzPHlj88dM+3sQKxmTq6vfY7jOBknyTTzbwGPSXodEMHaNF9JsPzmoFLPKUfOa/Jwn+M4DpBsFt8D\n4finUeGml/OWW2gPKvWcooydFPxNPiooY/JRnjjhOE7bkphAScpf3GwfSZjZr5Kqo2moxHMqRE7Y\n6uBJLbvjTlZMK71OXr+JExlw0omJ1ek4jlMNSYb49o+87gGMB54D2kegZk4OEh0KpZNXQtSTmven\noNwKRKqUCK2ZMQOAXvvvX3T/mhkzip7v4uU4TqNIMsTXYfVOSdsAv0mq/KYgl1JeaWivGGOO3zJG\nqohAFRKjUiLUa//9S4pMOXErJF4uWo7j1IO6ragrqSvwopk1fH2o1NLMJx8V/K8lvFeozGif1thJ\nHUSkmBjVQzTiiKGLleM45YibZp5kH1R0RolOwGhgSlLlty2hN7bsiTmsuO2nMOTRDqJQziNKkgEn\nnbhVPfli6R6W4zhJkWQf1JWR1xuAeWY2P8Hys0t+anlCBI3/o8BA1szoDqylV++FDRWlckRFq5iH\nlRPULNjrOJmm1LydbZjRm+RMEleY2XfLbWsEDQ/xFQjDVUN+A98hfLZyIf36zGLA+E8kG0KsM8vu\nuJOFF14IlO//aiain1WrXJOTAvmCNO9Pwf/8RKv87U0uVnFDfEkK1HNm9om8bS+Y2d6JVFABqQgU\nVCQccZMbOjR+CQlho8lda/T6stqox0nDhy2fVY5iWZHgAuYUodhE0oV+21Ehi4pVE7UDURomUJK+\nBvxfYGfgtciuvsCTZvalmiqogqwKVM3JDdXOUpERCgkVpNeAV5oBmU+/iRMBSgpaI5NYnCYh9zvO\nCU2l82/mn5+gUN329Fvc+/w7ZY8bvUM/Ljx6r6rraaRA9QcGAJcB50Z2rTSz92sqvEoaKlC5p6AR\nB8Gk31U0BqnqRqoe2YINpJRQ16vhrkSMkrTBMx+dzSQtLBU+sMYRn6ffCJrscTttW/K4ZhKofma2\nQlLBK0pDpGoRqIU/+AHr5rxcwQmzYO1yGLgr9B1S9gk8kcaoSUN9hSjnVVZzv9ISo7gUu2YXqhYm\nP5xX5e+2kMhcsPQc9vpoFjf0P4tHek0oem5c8Tl236GcPG54xbZVQiMFapqZTZT0BkGauSK7zcx2\nLnJq3Wi4QEGH7L26NzRNHuorRskkkQrIkhiVI2thTydhCoTzbts4PlYYrRCFRGb8mumcsfynAMzu\nNoYnex5aVKgaIT5xaHiSRJZoaIgvzXBbi68hFTdhoRDN1sC7V9XcFAufXbD0HEauf503u+68WTji\nejLFKCgydXxoLfQ77L7nKIacd17VZTbSg/pEqf1m9lxNFVRBQwQqC15MXv+X0xq4V5Ud4iYN5IvO\n+DXT+cyHj24Wp0sG/qjD8XXzZOow5KVQRKKZBOrRErvNzP6+pgqqoCEClZV+oKzY4SSOe1X1Jcmk\nAYiITkL9TVVRw4NzI5OXGjbVkZkdWmsZTUsW+n/quDSHky6FZunwmTkqp5gQxRGfcTttW7m3kxuv\nlEboPeZKCOUSibLyIJTkQN0eBOOhDiJIlvgjcJ2ZrU2kggpomAcF6QtUDvek2oL8mTlyZKExyQKF\nxKiUECUaastC2D9qS174P61JpgvR8MliCdZ9WglcE74/GbgFOCHBOtKnTvPu1Yx7Um1BrvHIz3Zs\nR68qrhhV5QVVQy2raSfMstd6s+LJneHRV+G2/aD3dqx5+W0gWx5SOZL0oF4ys9HltjWCunpQWfdU\nWjyzz9maQl5VMzQ+lVCJZ5RKKnUKCUuxJgUY1gM+Wg3desOQMZn5XqThQT0n6UAzeyo0YByQwqJM\nDSBt970UMRY5dCpjyitTmP769KL7J+w8gRN2Ty9QkO9VRZc9yUqDVAmZ84xKkT/OKUHPqdwwi9gL\nk25+qF4Ku6xOzL5GkKQHNQfYA3gr3DQcmEuw9IY1ctLYuntQkF2Bgux7eU1AVJRmLgq+S2MHb/3A\nV2xfmqLVTGnqmfeMylHDb60WAcoR6zPNUt9YSBqzmY8otd/M5iVSUQzqJlDNMu4og1/IZqCUKBUT\nnELeVdxz6009ppFKgqgoNZUYRYn5G6tkbs5CJPoZZSj8n8pMEpIGAMOIhA5baqBuhj7gWLgnFYuc\nyCQlLMWELgteVY56e1elxhjli1LmxagQeb+tZa/1LihEDZmbMy4ZesBOw4P6T+B0giU3coW21kDd\nZgjvRal1Wv8Wp5AwJS0i+XVc8KkLUu2vyhEn5ThHuUa00tRuaFJRyjFzMsuuOo8Vi4ZszuQtdQ8z\nFVrNyENrGgI1FxhjZh8lUmANJC5QzR4yy9CTU1aY8soULvnLJUBjvJtG11cJccNQi1auY+mqdVsd\ns+LD9QD069l187b1fECnLqvp2llbHV+MgT0HsssJk7LTmIdsdX8WzmLN28HwziyETCsiI21ZGgJ1\nF/A1M1ucSIE1kLhAZeSpoyayfg35S1+XIgH7Jz0wiZmLZjbUo8mqN5WjkCf0yRceY8zcp4GOQrTN\nxqX03/TB5uO6du7E8s7GUjYCsFKbAOhrnaBzV+jcrWTdK9evZK8wveqtXfoysOdAtuu5XdHj6y0G\nBT3MUcNg9Xub07b7ffl72RekYqTcHqQhUGOBe4EXgc2PWWZ2TAJlHwFcDXQGbjSzy0sdXxeBgub2\nPuq4CmdN9uSI2lWK/OMqvIacSMx9fy57bLsHk4+YXIHRyZCGN5XEvHPj10zn2M5/ZnDfHh0+hyms\nYrpWM1PBz36sdQdggvXmhHnPbz4OKPp5TXllCu/cOpk9n13CyvUrAejbtbBQ1atfp2jYc+VC+o34\nkAEDXthyLVl8yKuElMP/aQjUbOB6YBawKbfdzB6vsdzOwCvA4cB8YAbwRTN7qdg5iQlURtzhREnz\nmqKiVEiQ4vzoC5UR48fViP6mSkjSniQnPS25lEO+KPXpBX2HlL6GYp95ic+6nKeZ1KrVsRNHsh59\nqIWUwv9pCNQMM6tsZbl45X4KuMjM/iF8/z0AM7us2Dk1CdT9526ZLigr3kY9aFRGYilRqvWeVjBr\n9KQHJm32mrLU/1OuMa67+JQj7x5PGbo70zd9UF3GY4VefDWeZiXJH7FS79uh/zYqwNCQti6NmST+\nKOky4D46hvhqTTMfCrwdeT8fGJd/kKQzgDMAhg9PKDuoFYUpR27GiWlnB41G0tdZqDFK+n7mysnV\nk5tBI6+OKa9MYeaimYwdPDaVkF4p1i8bx5p5O7J951Es6HorVz75G6Y+uuPm/XWbcbscBUJAU/r1\nqS00mZtpO1p2rvwC34lc2bmQbHRbMQrNAF+MkvPR1XGGiMwRvbbo76gUQ8bAkSV7WhIhSQ+q0LpQ\nNaeZSzoeOMLMvhy+PxUYZ2bfKHZOQ1fUbWYKiQhULyJVhHMSo0RMPY2EiHziLPnwZrcrWau36WHD\n6L/xAAZsPBhocEp2oXBenteU2H3M94Ch6Hel4R5wK4f1ShE3WalGgWqZJd8bHuJrR6oJwxX6Iicd\nwquGSKM3ZcS+TO/Ti7nrlzckIaKSwalRcgKUiX6yvIa5Zq+pHDEfahqWAdmK/c4ZJK2ZJI4C9gJ6\n5LaZ2SU1ltmFIEliPPAOQZLEyWY2u9g5zS5Q5SYnLUYijUe5RIYcxfZl4GlzykPfYvq7f9ySVZZg\nw1qtCEF8Tyi1dPRIf8uUz5zeeBti9FHl+qXqEq5NcyXcNiONJInrgF7AocCNwPHAM2b2LwmUPQG4\niiDN/CYzu7TU8c0oUHEnJy1GXeZ/K+fuZ/QHvDkc9NFHTFi+jBP6jarY1mpXYU0yHNewdPSIMEzp\n25vpw/dh5pr59a83hj3AVok8iYf7fMaVhpOGQL1gZntH/vcB7jezv0ukggpoJoFKah64LM7/lgYd\nnrAHHVyy4ak1JNcIkp4nsGAdkz/L9PWLoVvvunidVVPEo0ksFJq1sYF1ptLITD0//zQE6mkzGyfp\nKeA4YCkw28x2TaSCCsi6QNVbTBrRqGWVggkRYUM3u9sYLhn4o83HNtN8ceU87KoeaP73Jlj9XgdR\nqqasulJCRGoKhbZBOC9fkCqJzNS77UhDoL5PsNz7eOBnBBPG/tzMLkikggrIqkA1WjjayauKzhAx\noMtI+iw7s8P+b73zLUZrHm9334Unex7KI70mANkSobjU0vBsdc66DdCtNxN2+DtO+NxPkjU0SUpk\nnOYSOWL1S7V4OC+pB5ly5YzadhTfPeC7VduZahafpO5ADzNbnnjhMciSQGVFJDKRIZYwtz39Fr+c\ndTvLOz/Dmk6vANBr0+68v2gv1n8wroNntHmantWvtlx2VlVJNSsXMuGt/+WEbfdrrntRJIlnkhYx\nl4/Yg25bxLZcpmkLeU31fPgt9P1qaoFKmzQFqtzTbdqikPUJS4tRqL/ouWX302P7e4BAmGKNHWrX\n8S35NNvaZoWICFCh+QAnLH2XE1auzmSmaZJkeab8YrhAVSlQVzxzBS+//zJQWkyKPbUm0T/QCPK/\n1JAdO8utL7Ss8xMdvKaK+x5aOMRTlhYf55M/xOCC7Q/LduiyCoo9BDfLwya4QNUsUOXi+qX2Z6Wh\nL0fSHe+VUE0G3bH7DqXrgKeTeVpshznWCtEmHmQzehWlSPO3Wg8aJlCSPlFqf7Mu+R4nrt9sX4pS\n1Cs0Wet4oro+LbZJYw20vOdUiGbPZi0lSs1yDcVopEAVmoMvR2st+d5GxB04vHjFOpasXteh/ydK\nufFEXQc8XfJBoK5Pi+3UaLeTGOcR57ucdoNf6IG41UQpiof4XKBqIur55Pp8CrHiw/V06f0GECQq\nFGJQ7+58rF/3gvvipEjX/YfZ6o13u4YzCxBHCNIgq8JZLxq+3IakrsDXgNxj9GPA9Wa2Pqk6nMqI\ns5ZQMaKez4CNBxf0jgDoDMOHvMASe6qqejLRP5BbbiC3DlirCVQu1bqVl4yIyQm7n7DVd63auS+T\nJBO/gwyS5EDdG4GuwM3hplOBjbllMhpJq3hQtQgMxF/IrhjNOIi1Jloh9TpKO4UwnaYijQUL9zez\nfSLv/yDpfxMsv2WIKzy1CkxdFrJrZXKLOM6a2hoCFRUn956cJiRJgdooaRczew1A0s7AxgTLbwqS\nXKLbBabBjJ20pVGffFTz9ke55+S0CEmG+MYDk4HXAQEjgElmVirLry40IsRXbfp0jnoLz/r165k/\nfz5r166tWx0tybpVsH4NbPwIOneDPh9L26LKWbV4i/1de0H3Pmlb5IT06NGDHXfcka5du6ZtSqo0\nNMQnqRPwIbAbsEe4ea6ZrUui/LSoZjBpVrye+fPn07dvX0aOHImkVG1pSpa8Ch+tgv7bQe9BaVsT\nn9VLYPla6NYHBu2WtjVOBDNj6dKlzJ8/n5122iltc5qCRATKzDZJ+pmZ7Qe8kESZaXHxb2fz0rsr\ngNLeUFaEqBhr1651caqFngMCgfpwWXMJ1IfLgv89B6Rrh7MVkhg4cCDvvfde2qY0DUn2QT0i6QvA\n3dYig6uyLkLlcHGqgd6DgsZ+/YeBN9VzQLaFavWSLfZ265NtW9sY/01WRpIC9RXg34ANktYS9EOZ\nmfVLsI66c+HRe6VtgpMVcl7I+g+D/1lu9HPi1LWne09Oy9ApqYLMrK+ZdTKzbmbWL3zfVOLkNIb3\n33+fww8/nN12243DDz+cZcuWpW1SYXoPCvpxuvbc4kmtXlJxMVOmTGGvvfaiU6dO1CV5Z/WSIBzZ\ntWdgbwaE9JxzzmHUqFHsvffefP7zn+eDDz5I26S68/3vf5+9996bfffdl8997nO8++67aZvU9CQm\nUJIeibPNcS6//HLGjx/Pq6++yvjx47n88svTNqk0PQdsEakPKxfTj3/849x9990cfHCR2ThqJYP9\nTocffjgvvvgiL7zwArvvvjuXXXZZ2ibVnXPOOYcXXniB559/nokTJ3LJJZekbVLTU3OIT1IPoBcw\nSNIAgtAeQD9gaK3lO7UTTfxIitE79CsbDn3zzTeZOHEiL774IgBXXnklq1at4t577+Wxxx4D4LTT\nTuOQQw7hiiuuSNQ+7j93y9RFSbHNMBj3lcBjKeClFLveiy66KFk7ckT6na549Q5e/nBRosXHWTU1\nzjUfeOCBTJ06tUgJ1bPwBz9g3ZyXEy2z+56jGHLeeSWPiXPNq1ev9v6mBEiiD+orwNnADsCzbBGo\nFcC1CZTvtBiLFi1i++23B2DIkCEsWpRsw1o3uoQT3qad2ZcTpo9WBe+79YEuPdKzpww33XQTJ510\nUtpmNITzzz+fX/3qV/Tv359HH234ENCWo2aBMrOrgaslnWlm1yRgk5MwWU78kFSfJ80j6xQ2XPJq\n+pl90Wy90IbvDrqo8XbE4NJLL6VLly6ccsopiZddztNJg0svvZRLL72Uyy67jGuvvZaLL744bZOa\nmiSTJK6R9GlJJ0v659xfUuU7zUeXLl3YtGnT5ve5WS0GDx7MggULAFiwYAEf+1gTzdaQ64/6aBUs\nf7tD4kSx602UjCVElLrmX/7yl0ybNo1f//rXLRXuivM5n3LKKdx1112NNKslSTJJ4hbgSuAgYP/w\nL70FVpzUGTx4MIsXL2bp0qWsW7eOadOmAXDMMcdw883BpPc333wzxx57bJpmVkYus6//sMCDiSRO\nFLveRMlYQkSxa37ggQf44Q9/yH333UevXr1StjJZil3zq6++uvmYe++9l1GjRqVlYsuQ5DioscDo\nVhmk69RO165dueCCCzjggAMYOnTo5h/sueeey4knnsgvfvELRowYwZ133pmypVXQe1Dwl5sSafUS\nuvYeVPB677nnHs4880zee+89jjrqKPbdd18efPDByurL6EDcYp/xN77xDdatW8fhhx8OBIkS1113\nXZqmJkap7/XcuXPp1KkTI0aMaJnrTZMkJ4udApxlZgsSKbAGWmU9qFqYM2cOe+65Z9pmtD6rlwSh\nPujQJ5Ro+fkJEVmf1cIpif8201kPahDwkqRngM2TxJrZMQnW4TjZIicUOe8muq1W6i1+jpNxkhSo\nixIsy3Gah2i4L4kMv3yvqf8wFyanLUlMoMzscUmDCZIjAJ4xs8W1lCnpIuBfgdz0v+eZ2fRaymwn\nzKylsqcyTy5x4aNVW8SlEmHxcF7L4130lZGYQEk6EfgR8BjBYN1rJJ1jZrUOIf+JmV1Zq33tRo8e\nPVi6dCkDBw50kWoUOU8qF5pb/nbHqZFKiY2H81qe3HpQPXpkd1B11kgyxHc+sH/Oa5K0HfAwkPwc\nJ05ZdtxxR+bPn+9rz6TFug3ByryEU0xtCMfKFJvxIbe/57bQfQNB0MA/u1Yjt6KuE48kBapTXkhv\nKcmMszozHPA7E/i2mRWcrVPSGcAZAMOHN+f6TUnStWtXX7UzS8ycDLPKPKuNOR72/YfG2OM4TUCS\naeY/AvYGbg83nQTMMrPvlDnvYWBIgV3nA08BSwAD/hPY3sz+TzlbPM3ccRwnuzQ8zdzMzpF0HMFM\nEgA3mNk9Mc47LE75kn4O1GFovuM4jpNFklhuY1dgsJk9aWZ3A3eH2w+StIuZvVZD2dtHBv5+Hnix\nVnsdx3Gc5qDmEJ+kacD3zGxW3vYxwA/M7Ogayr4F2JcgxPcm8JU4M1VIeg+YV229BIOOK186tfXw\n++D3IIffhwC/DwG13ocRZrZduYOSEKgZZrZ/kX2zzGxMTRWkgKSZceKjrY7fB78HOfw+BPh9CGjU\nfUgiy26bEvt6JlC+4ziO04YkIVAzJf1r/kZJXyZYYddxHMdxKiaJLL6zgXskncIWQRoLdCNIbGhG\nbkjbgIzg98HvQQ6/DwF+HwIach+SHAd1KPDx8O1sM/tDIgU7juM4bUliAuU4juM4SZLYku+O4ziO\nkyQuUBEkHSFprqS/STo3bXvSQNIwSY9KeknSbEnfTNumNJHUWdJfw/F+bYmkbSRNlfSypDmSPpW2\nTY1G0rfC38OLkm6X1BZTkku6SdJiSS9Gtm0r6feSXg3/D6hX/S5QIZI6Az8DjgRGA1+UNDpdq1Jh\nA8GkvKOBA4Gvt+l9yPFNYE7aRqTM1cADZjYK2Ic2ux+ShgJnAWPN7ONAZ+Cf0rWqYfwSOCJv27nA\nI2a2G/BI+L4uuEBt4QDgb2b2upl9BPwGODZlmxqOmS0ws+fC1ysJGqOh6VqVDpJ2BI4CbkzblrSQ\n1B84GPgFgJl9ZGYfpGtVKnQBekrqAvQC3k3ZnoZgZk8A7+dtPha4OXx9M/CP9arfBWoLQ4G3I+/n\n06YNcw5JI4H9gKfTtSQ1rgK+A2xK25AU2YlgYarJYajzRkm90zaqkZjZO8CVwFvAAmC5mT2UrlWp\nMjgy5dxCYHC9KnKBcgoiqQ9wF3C2ma1I255GI2kisNjM2n2weRfgE8B/m9l+wGrqGNLJImEfy7EE\nYr0D0FvSl9K1KhtYkAZet1RwF6gtvAMMi7zfMdzWdkjqSiBOvw5nqG9HPgMcI+lNgnDv30u6NV2T\nUmE+MN/Mcl70VALBaicOA94ws/fMbD3Big2fTtmmNFkkaXsIVpwAFpc5vmpcoLYwA9hN0k6SuhF0\ngt6Xsk0NR5II+hvmmNmP07YnLczse2a2o5mNJPgu/MHM2u6p2cwWAm9L2iPcNB54KUWT0uAt4EBJ\nvcLfx3jaLFEkj/uA08LXpwH31quiJJd8b2rMbIOkbwAPEmTp3GRms1M2Kw0+A5wKzJL0fLjtPDOb\nnrPBg1UAAAGWSURBVKJNTrqcCfw6fHB7HZiUsj0NxcyeljQVeI4gy/WvtMmUR5JuBw4BBkmaD1wI\nXA7cKelfCJY1OrFu9ftMEo7jOE4W8RCf4ziOk0lcoBzHcZxM4gLlOI7jZBIXKMdxHCeTuEA5juM4\nmcTTzB2nwUgaSDDJJsAQYCPBdEIAa8ysnQeBOs5mPM3ccVJE0kXAKjO7Mm1bHCdreIjPcTKEpFXh\n/0MkPS7pXkmvS7pc0imSnpE0S9Iu4XHbSbpL0ozw7zPpXoHjJIcLlONkl32ArwJ7EszusbuZHUCw\n/MeZ4TFXAz8xs/2BL9DGS4M4rYf3QTlOdpmRW9ZA0mtAbomHWcCh4evDgNHBFHEA9JPUx8xWNdRS\nx6kDLlCOk13WRV5virzfxJbfbifgQDNb20jDHKcReIjPcZqbh9gS7kPSvina4jiJ4gLlOM3NWcBY\nSS9Ieomgz8pxWgJPM3ccx3EyiXtQjuM4TiZxgXIcx3EyiQuU4ziOk0lcoBzHcZxM4gLlOI7jZBIX\nKMdxHCeTuEA5juM4meT/Bz+rHEKUiISwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc8e09cc128>"
      ]
     },
     "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_xlabel(\"Time\")\n",
    "ax1.set_ylabel(\"Control amplitude\")\n",
    "for j in range(n_ctrls):\n",
    "    ax1.step(result.time, \n",
    "             np.hstack((result.initial_amps[:, j], result.initial_amps[-1, j])), \n",
    "             where='post')\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",
    "for j in range(n_ctrls):\n",
    "    ax2.step(result.time, \n",
    "             np.hstack((result.final_amps[:, j], result.final_amps[-1, j])), \n",
    "             where='post', label='u{}'.format(j))\n",
    "ax2.legend(loc=8, ncol=n_ctrls)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Versions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table><tr><th>Software</th><th>Version</th></tr><tr><td>QuTiP</td><td>4.1.0</td></tr><tr><td>Numpy</td><td>1.11.3</td></tr><tr><td>SciPy</td><td>0.18.1</td></tr><tr><td>matplotlib</td><td>2.0.0</td></tr><tr><td>Cython</td><td>0.25.2</td></tr><tr><td>Number of CPUs</td><td>4</td></tr><tr><td>BLAS Info</td><td>INTEL MKL</td></tr><tr><td>IPython</td><td>5.1.0</td></tr><tr><td>Python</td><td>3.6.0 |Anaconda 4.3.1 (64-bit)| (default, Dec 23 2016, 12:22:00) \n",
       "[GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]</td></tr><tr><td>OS</td><td>posix [linux]</td></tr><tr><td colspan='2'>Fri Jul 14 16:41:51 2017 BST</td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from qutip.ipynbtools import version_table\n",
    "\n",
    "version_table()"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "References:\n",
    "\n",
    "3.  Doria, P., Calarco, T. & Montangero, S. \n",
    "    Optimal Control Technique for Many-Body Quantum Dynamics. \n",
    "    Phys. Rev. Lett. 106, 1–4 (2011).\n",
    "\n",
    "4.  Caneva, T., Calarco, T. & Montangero, S. \n",
    "    Chopped random-basis quantum optimization. \n",
    "    Phys. Rev. A - At. Mol. Opt. Phys. 84, (2011)."
   ]
  }
 ],
 "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}