examples/Simulated_calibration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulated calibration\n",
    "\n",
    "Calibration of control pulses is the process of fine-tuning parameters in a feedback-loop with the experiment. We will simulate this process here by constructing a black-box simulation and interacting with it exactly like an experiment.\n",
    "\n",
    "We have manange imports and creation of the black-box the same way as in the previous example in a helper `single_qubit_blackbox_exp.py`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2022-05-18 12:39:33.694052: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcudart.so.11.0'; dlerror: libcudart.so.11.0: cannot open shared object file: No such file or directory\n",
      "2022-05-18 12:39:33.694082: I tensorflow/stream_executor/cuda/cudart_stub.cc:29] Ignore above cudart dlerror if you do not have a GPU set up on your machine.\n",
      "2022-05-18 12:39:35.129198: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcuda.so.1'; dlerror: libcuda.so.1: cannot open shared object file: No such file or directory\n",
      "2022-05-18 12:39:35.129237: W tensorflow/stream_executor/cuda/cuda_driver.cc:269] failed call to cuInit: UNKNOWN ERROR (303)\n",
      "2022-05-18 12:39:35.129256: I tensorflow/stream_executor/cuda/cuda_diagnostics.cc:156] kernel driver does not appear to be running on this host (knecht): /proc/driver/nvidia/version does not exist\n",
      "2022-05-18 12:39:35.129516: I tensorflow/core/platform/cpu_feature_guard.cc:151] This TensorFlow binary is optimized with oneAPI Deep Neural Network Library (oneDNN) to use the following CPU instructions in performance-critical operations:  AVX2 FMA\n",
      "To enable them in other operations, rebuild TensorFlow with the appropriate compiler flags.\n"
     ]
    }
   ],
   "source": [
    "from single_qubit_blackbox_exp import create_experiment\n",
    "\n",
    "blackbox = create_experiment()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This blackbox is constructed the same way as in the C1 example. The difference will be in how we interact with it. First, we decide on what experiment we want to perform and need to specify it as a python function. A general, minimal example would be\n",
    "\n",
    "`def exp_communication(params):\n",
    "    # Send parameters to experiment controller\n",
    "    # and recieve a measurement result.\n",
    "    return measurement_result\n",
    " `\n",
    "\n",
    "Again, `params` is a linear vector of bare numbers. The measurement result can be a single number or a set of results. It can also include additional information about statistics, like averaging, standard deviation, etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ORBIT - Single-length randomized benchmarking\n",
    "The following defines an [ORBIT](https://arxiv.org/abs/1403.0035) procedure. In short, we define sequences of gates that result in an identity gate if our individual gates are perfect. Any deviation from identity gives us a measure of the imperfections in our gates. Our helper `qt_utils` provides these sequences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from c3.utils import qt_utils"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[['rx90p[0]',\n",
       "  'ry90p[0]',\n",
       "  'ry90p[0]',\n",
       "  'ry90p[0]',\n",
       "  'rx90p[0]',\n",
       "  'ry90p[0]',\n",
       "  'rx90p[0]',\n",
       "  'rx90p[0]',\n",
       "  'rx90m[0]',\n",
       "  'ry90m[0]',\n",
       "  'rx90p[0]',\n",
       "  'rx90p[0]']]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qt_utils.single_length_RB(\n",
    "            RB_number=1, RB_length=5, target=0\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The desired number of 5 gates is selected from a specific set (the Clifford group) and has to be decomposed into the available gate-set. Here, this means 4 gates per Clifford, hence a sequence of 20 gates."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Communication with the experiment\n",
    "Some of the following code is specific to the fact that this a *simulated* calibration. The interface of $C^2$ to the experiment is simple: parameters in $\\rightarrow$ results out. Thus, we have to wrap the blackbox by defining the target states and the `opt_map`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import tensorflow as tf\n",
    "\n",
    "def ORBIT_wrapper(p):\n",
    "    def ORBIT(params, exp, opt_map, qubit_labels, logdir):    \n",
    "        ### ORBIT meta-parameters ###\n",
    "        RB_length = 60 # How long each sequence is\n",
    "        RB_number = 40  # How many sequences\n",
    "        shots = 1000    # How many averages per readout\n",
    "\n",
    "        ################################\n",
    "        ### Simulation specific part ###\n",
    "        ################################\n",
    "\n",
    "        do_noise = False  # Whether to add artificial noise to the results\n",
    "\n",
    "        qubit_label = list(qubit_labels.keys())[0]\n",
    "        state_labels = qubit_labels[qubit_label]\n",
    "        state_label = [tuple(l) for l in state_labels]\n",
    "\n",
    "        # Creating the RB sequences #\n",
    "        seqs = qt_utils.single_length_RB(\n",
    "                RB_number=RB_number, RB_length=RB_length, target=0\n",
    "        )\n",
    "\n",
    "        # Transmitting the parameters to the experiment #\n",
    "        exp.pmap.set_parameters(params, opt_map)\n",
    "        exp.set_opt_gates_seq(seqs)\n",
    "\n",
    "        # Simulating the gates #\n",
    "        U_dict = exp.compute_propagators()\n",
    "\n",
    "        # Running the RB sequences and read-out the results #\n",
    "        pops = exp.evaluate(seqs)\n",
    "        pop1s, _ = exp.process(pops, labels=state_label)\n",
    "\n",
    "        results = []\n",
    "        results_std = []\n",
    "        shots_nums = []\n",
    "\n",
    "        # Collecting results and statistics, add noise #\n",
    "        if do_noise:\n",
    "            for p1 in pop1s:\n",
    "                draws = tf.keras.backend.random_binomial(\n",
    "                    [shots],\n",
    "                    p=p1[0],\n",
    "                    dtype=tf.float64,\n",
    "                )\n",
    "                results.append([np.mean(draws)])\n",
    "                results_std.append([np.std(draws)/np.sqrt(shots)])\n",
    "                shots_nums.append([shots])\n",
    "        else:\n",
    "            for p1 in pop1s:\n",
    "                results.append(p1.numpy())\n",
    "                results_std.append([0])\n",
    "                shots_nums.append([shots])\n",
    "\n",
    "        #######################################\n",
    "        ### End of Simulation specific part ###\n",
    "        #######################################\n",
    "\n",
    "        goal = np.mean(results)\n",
    "        return goal, results, results_std, seqs, shots_nums\n",
    "    return ORBIT(\n",
    "                p, blackbox, gateset_opt_map, state_labels, \"/tmp/c3logs/blackbox\"\n",
    "            )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimization\n",
    "We first import algorithms and the correct optimizer object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import copy\n",
    "\n",
    "from c3.experiment import Experiment as Exp\n",
    "from c3.c3objs import Quantity as Qty\n",
    "from c3.parametermap import ParameterMap as PMap\n",
    "from c3.libraries import algorithms, envelopes\n",
    "from c3.signal import gates, pulse\n",
    "from c3.optimizers.calibration import Calibration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Representation of the experiment within $C^3$\n",
    "At this point we have to make sure that the gates (\"RX90p\", etc.) and drive line (\"d1\") are compatible to the experiment controller operating the blackbox. We mirror the blackbox by creating an experiment in the $C^3$ context:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_final = 7e-9   # Time for single qubit gates\n",
    "sideband = 50e6\n",
    "lo_freq = 5e9 + sideband\n",
    "\n",
    " # ### MAKE GATESET\n",
    "gauss_params_single = {\n",
    "    'amp': Qty(\n",
    "        value=0.45,\n",
    "        min_val=0.4,\n",
    "        max_val=0.6,\n",
    "        unit=\"V\"\n",
    "    ),\n",
    "    't_final': Qty(\n",
    "        value=t_final,\n",
    "        min_val=0.5 * t_final,\n",
    "        max_val=1.5 * t_final,\n",
    "        unit=\"s\"\n",
    "    ),\n",
    "    'sigma': Qty(\n",
    "        value=t_final / 4,\n",
    "        min_val=t_final / 8,\n",
    "        max_val=t_final / 2,\n",
    "        unit=\"s\"\n",
    "    ),\n",
    "    'xy_angle': Qty(\n",
    "        value=0.0,\n",
    "        min_val=-0.5 * np.pi,\n",
    "        max_val=2.5 * np.pi,\n",
    "        unit='rad'\n",
    "    ),\n",
    "    'freq_offset': Qty(\n",
    "        value=-sideband - 0.5e6,\n",
    "        min_val=-53 * 1e6,\n",
    "        max_val=-47 * 1e6,\n",
    "        unit='Hz 2pi'\n",
    "    ),\n",
    "    'delta': Qty(\n",
    "        value=-1,\n",
    "        min_val=-5,\n",
    "        max_val=3,\n",
    "        unit=\"\"\n",
    "    )\n",
    "}\n",
    "\n",
    "gauss_env_single = pulse.EnvelopeDrag(\n",
    "    name=\"gauss\",\n",
    "    desc=\"Gaussian comp for single-qubit gates\",\n",
    "    params=gauss_params_single,\n",
    "    shape=envelopes.gaussian_nonorm\n",
    ")\n",
    "nodrive_env = pulse.Envelope(\n",
    "    name=\"no_drive\",\n",
    "    params={\n",
    "        't_final': Qty(\n",
    "            value=t_final,\n",
    "            min_val=0.5 * t_final,\n",
    "            max_val=1.5 * t_final,\n",
    "            unit=\"s\"\n",
    "        )\n",
    "    },\n",
    "    shape=envelopes.no_drive\n",
    ")\n",
    "carrier_parameters = {\n",
    "    'freq': Qty(\n",
    "        value=lo_freq,\n",
    "        min_val=4.5e9,\n",
    "        max_val=6e9,\n",
    "        unit='Hz 2pi'\n",
    "    ),\n",
    "    'framechange': Qty(\n",
    "        value=0.0,\n",
    "        min_val= -np.pi,\n",
    "        max_val= 3 * np.pi,\n",
    "        unit='rad'\n",
    "    )\n",
    "}\n",
    "carr = pulse.Carrier(\n",
    "    name=\"carrier\",\n",
    "    desc=\"Frequency of the local oscillator\",\n",
    "    params=carrier_parameters\n",
    ")\n",
    "\n",
    "rx90p = gates.Instruction(\n",
    "    name=\"rx90p\",\n",
    "    t_start=0.0,\n",
    "    t_end=t_final,\n",
    "    channels=[\"d1\"],\n",
    "    targets=[0]\n",
    ")\n",
    "QId = gates.Instruction(\n",
    "    name=\"id\",\n",
    "    t_start=0.0,\n",
    "    t_end=t_final,\n",
    "    channels=[\"d1\"],\n",
    "    targets=[0]\n",
    ")\n",
    "\n",
    "rx90p.add_component(gauss_env_single, \"d1\")\n",
    "rx90p.add_component(carr, \"d1\")\n",
    "QId.add_component(nodrive_env, \"d1\")\n",
    "QId.add_component(copy.deepcopy(carr), \"d1\")\n",
    "QId.comps['d1']['carrier'].params['framechange'].set_value(\n",
    "    (-sideband * t_final * 2 * np.pi) % (2*np.pi)\n",
    ")\n",
    "ry90p = copy.deepcopy(rx90p)\n",
    "ry90p.name = \"ry90p\"\n",
    "rx90m = copy.deepcopy(rx90p)\n",
    "rx90m.name = \"rx90m\"\n",
    "ry90m = copy.deepcopy(rx90p)\n",
    "ry90m.name = \"ry90m\"\n",
    "ry90p.comps['d1']['gauss'].params['xy_angle'].set_value(0.5 * np.pi)\n",
    "rx90m.comps['d1']['gauss'].params['xy_angle'].set_value(np.pi)\n",
    "ry90m.comps['d1']['gauss'].params['xy_angle'].set_value(1.5 * np.pi)\n",
    "\n",
    "parameter_map = PMap(instructions=[QId, rx90p, ry90p, rx90m, ry90m])\n",
    "\n",
    "# ### MAKE EXPERIMENT\n",
    "exp = Exp(pmap=parameter_map)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we define the parameters we whish to calibrate. See how these gate instructions are defined in the experiment setup example or in `single_qubit_blackbox_exp.py`. Our gate-set is made up of 4 gates, rotations of 90 degrees around the $x$ and $y$-axis in positive and negative direction. While it is possible to optimize each parameters of each gate individually, in this example all four gates share parameters. They only differ in the phase $\\phi_{xy}$ that is set in the definitions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "gateset_opt_map =   [\n",
    "    [\n",
    "      (\"rx90p[0]\", \"d1\", \"gauss\", \"amp\"),\n",
    "      (\"ry90p[0]\", \"d1\", \"gauss\", \"amp\"),\n",
    "      (\"rx90m[0]\", \"d1\", \"gauss\", \"amp\"),\n",
    "      (\"ry90m[0]\", \"d1\", \"gauss\", \"amp\")\n",
    "    ],\n",
    "    [\n",
    "      (\"rx90p[0]\", \"d1\", \"gauss\", \"delta\"),\n",
    "      (\"ry90p[0]\", \"d1\", \"gauss\", \"delta\"),\n",
    "      (\"rx90m[0]\", \"d1\", \"gauss\", \"delta\"),\n",
    "      (\"ry90m[0]\", \"d1\", \"gauss\", \"delta\")\n",
    "    ],\n",
    "    [\n",
    "      (\"rx90p[0]\", \"d1\", \"gauss\", \"freq_offset\"),\n",
    "      (\"ry90p[0]\", \"d1\", \"gauss\", \"freq_offset\"),\n",
    "      (\"rx90m[0]\", \"d1\", \"gauss\", \"freq_offset\"),\n",
    "      (\"ry90m[0]\", \"d1\", \"gauss\", \"freq_offset\")\n",
    "    ],\n",
    "    [\n",
    "      (\"id[0]\", \"d1\", \"carrier\", \"framechange\")\n",
    "    ]\n",
    "  ]\n",
    "\n",
    "parameter_map.set_opt_map(gateset_opt_map)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As defined above, we have 16 parameters where 4 share their numerical value. This leaves 4 values to optimize."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rx90p[0]-d1-gauss-amp                 : 450.000 mV \n",
      "ry90p[0]-d1-gauss-amp\n",
      "rx90m[0]-d1-gauss-amp\n",
      "ry90m[0]-d1-gauss-amp\n",
      "\n",
      "rx90p[0]-d1-gauss-delta               : -1.000  \n",
      "ry90p[0]-d1-gauss-delta\n",
      "rx90m[0]-d1-gauss-delta\n",
      "ry90m[0]-d1-gauss-delta\n",
      "\n",
      "rx90p[0]-d1-gauss-freq_offset         : -50.500 MHz 2pi \n",
      "ry90p[0]-d1-gauss-freq_offset\n",
      "rx90m[0]-d1-gauss-freq_offset\n",
      "ry90m[0]-d1-gauss-freq_offset\n",
      "\n",
      "id[0]-d1-carrier-framechange          : 4.084 rad \n",
      "\n"
     ]
    }
   ],
   "source": [
    "parameter_map.print_parameters()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is important to note that in this example, we are transmitting only these four parameters to the experiment. We don't know how the blackbox will implement the pulse shapes and care has to be taken that the parameters are understood on the other end. Optionally, we could specifiy a virtual AWG within $C^3$ and transmit pixilated pulse shapes directly to the physiscal AWG."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Algorithms\n",
    "As an optimization algoritm, we choose [CMA-Es](https://en.wikipedia.org/wiki/CMA-ES) and set up some options specific to this algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Maximum Evaluations\n",
    "We set the maximum number of evaluations to 30 for quick demonstration. Ideally around 250 function evaluations are useful for a good optimization. See [docs](https://c3-toolset.readthedocs.io/en/latest/Simulated_calibration.html#analysis) for a longer run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "alg_options = {\n",
    "    \"popsize\" : 10,\n",
    "    \"maxfevals\" : 300,\n",
    "    \"init_point\" : \"True\",\n",
    "    \"tolfun\" : 0.01,\n",
    "    \"spread\" : 0.25\n",
    "  }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define the subspace as both excited states $\\{|1>,|2>\\}$, assuming read-out can distinguish between 0, 1 and 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "state_labels = {\n",
    "      \"excited\" : [(1,), (2,)]\n",
    "  }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the real world, this setup needs to be handled in the experiment controller side.\n",
    "We construct the optimizer object with the options we setup:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import tempfile\n",
    "\n",
    "# Create a temporary directory to store logfiles, modify as needed\n",
    "log_dir = os.path.join(tempfile.TemporaryDirectory().name, \"c3logs\")\n",
    "\n",
    "opt = Calibration(\n",
    "    dir_path=log_dir,\n",
    "    run_name=\"ORBIT_cal\",\n",
    "    eval_func=ORBIT_wrapper,\n",
    "    pmap=parameter_map,\n",
    "    exp_right=exp,\n",
    "    algorithm=algorithms.cmaes,\n",
    "    options=alg_options\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And run the calibration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = parameter_map.get_parameters_scaled()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C3:STATUS:Saving as: /tmp/tmpgwxh9a4c/c3logs/ORBIT_cal/2022_05_18_T_12_39_37/calibration.log\n",
      "(5_w,10)-aCMA-ES (mu_w=3.2,w_1=45%) in dimension 4 (seed=987495, Wed May 18 12:39:37 2022)\n",
      "C3:STATUS:Adding initial point to CMA sample.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/ubuntu/.local/lib/python3.8/site-packages/cma/utilities/utils.py:343: UserWarning: input x0 should be a list or 1-D array, trying to flatten (4, 1)-array ()\n",
      "  warnings.warn(msg + ' (' +\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iterat #Fevals   function value  axis ratio  sigma  min&max std  t[m:s]\n",
      "    1     11 5.018774219790417e-01 1.0e+00 2.12e-01  2e-01  2e-01 0:20.5\n",
      "    2     21 5.257349992937843e-01 1.3e+00 2.09e-01  2e-01  2e-01 0:39.1\n",
      "    3     31 5.147175357658132e-01 1.4e+00 1.92e-01  1e-01  2e-01 0:57.7\n",
      "    4     41 5.211210083125237e-01 1.6e+00 1.61e-01  1e-01  2e-01 1:16.2\n",
      "    5     51 4.042151811973714e-01 1.7e+00 1.77e-01  1e-01  2e-01 1:34.8\n",
      "    6     61 2.671535011626892e-01 1.9e+00 1.99e-01  1e-01  2e-01 1:53.4\n",
      "    7     71 3.449590715074073e-01 1.8e+00 2.22e-01  1e-01  2e-01 2:12.0\n",
      "    8     81 2.830766798739788e-01 2.4e+00 2.01e-01  1e-01  2e-01 2:30.6\n",
      "    9     91 3.698922041864815e-01 2.5e+00 2.18e-01  1e-01  3e-01 2:49.3\n",
      "   10    101 2.901095186143745e-01 3.0e+00 2.09e-01  1e-01  2e-01 3:07.9\n",
      "   11    111 3.390430736823523e-01 3.1e+00 2.18e-01  1e-01  3e-01 3:26.5\n",
      "   12    121 3.183914858462594e-01 3.3e+00 1.81e-01  7e-02  2e-01 3:45.1\n",
      "   13    131 2.737803095670904e-01 3.0e+00 1.54e-01  6e-02  2e-01 4:03.8\n",
      "   14    141 2.941330001634861e-01 3.1e+00 1.36e-01  4e-02  2e-01 4:22.3\n",
      "   15    151 2.698400240491868e-01 3.6e+00 1.14e-01  4e-02  1e-01 4:40.9\n",
      "   16    161 2.616957775258201e-01 3.9e+00 1.17e-01  4e-02  1e-01 4:59.5\n",
      "   17    171 2.339378829266300e-01 4.3e+00 1.09e-01  4e-02  1e-01 5:18.1\n",
      "   18    181 2.370179348635957e-01 5.2e+00 9.11e-02  3e-02  1e-01 5:36.8\n",
      "   19    191 2.289690729072825e-01 5.0e+00 8.29e-02  2e-02  1e-01 5:55.4\n",
      "   21    211 2.273259675209152e-01 5.3e+00 6.82e-02  2e-02  8e-02 6:33.6\n",
      "   23    231 2.313063798702339e-01 6.1e+00 7.75e-02  2e-02  1e-01 7:10.9\n",
      "   25    251 1.929736338163719e-01 6.0e+00 6.08e-02  2e-02  8e-02 7:48.1\n",
      "   27    271 2.123614169199560e-01 5.8e+00 4.40e-02  1e-02  5e-02 8:25.3\n",
      "   29    291 2.281454056009735e-01 6.4e+00 3.92e-02  7e-03  4e-02 9:02.6\n",
      "   30    301 2.039788087280569e-01 6.9e+00 4.22e-02  7e-03  5e-02 9:21.4\n",
      "termination on maxfevals=300\n",
      "final/bestever f-value = 2.039788e-01 1.929736e-01\n",
      "incumbent solution: [-1.000649165824243, 0.026755245328456176, 0.1536340057227788, 0.3963135893003173]\n",
      "std deviation: [0.0073345213313633926, 0.010271003725077475, 0.054213069570127025, 0.025259204881618324]\n"
     ]
    }
   ],
   "source": [
    "opt.optimize_controls()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis\n",
    "The following code uses matplotlib to create an ORBIT plot from the logfile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install -q -U pip\n",
    "!pip install -q matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fd3b7accf40>]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAERCAYAAABRpiGMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABVaklEQVR4nO2dfXxU1Z3/318SEgJJSAgmmBAJaEABDQ2xEKRoUbelu30U69Zua5+0Wqlt1634sNtWq1a0/vqEFWu71ba2a2Xttt0ta1UENhrSEiTIc5TnICCQAAkhySTn98edm8Y4mTk3OTf3Tjjv12tekDtn7nzmnJn7ved8H44opbBYLBaLpb8MC1qAxWKxWJIba0gsFovFMiCsIbFYLBbLgLCGxGKxWCwDwhoSi8VisQwIa0gsFovFMiBSgxYQBGPHjlUlJSXdf7e3t5OWlhacIE2sTrNYnWaxOs0SRp21tbVHlFJn9T5+RhqSkpIS1q1b1/13Y2Mjubm5ASrSw+o0i9VpFqvTLGHUKSJ7Yh23S1vAyZMng5aghdVpFqvTLFanWZJFJ1hDAsDOnTuDlqCF1WkWq9MsVqdZkkUnWENisVgslgFiDQmOzyQZsDrNYnWaxeo0S7LoBGtIABgzZkzQErSwOs1idZrF6jRLsugEa0gAWL9+fdAStLA6zWJ1msXqNEuy6IQADYmI5IvIUyKyPfpYLiLjNV73GRE5KCIbej02i4gSkcsHQ78pGppauePZjTQ0tQYtxWKxWPpFIIZERNKA54E0YBowFWgBXhKRTI1TLFNKzej5AO4FDgCrvOoJKla7oamVT/xkLUea2/nET9YmNCZhiynvC6vTLFanWaxO80gQG1uJyPXAT4BzlVI7o8fGAQ3A7Uqph+K8dhKQqZTa2Ov4i8ArSql/S/T+FRUVqmdCYlDc8exGjjS38/inK7j+F+sYm5nGdz52UdCyLBaLJSYiUquUquh9PKilrauAva4RAVBKHQS2RJ/rE6XUzhhGZBJwKfDT/ohZvXp1f142YBbNL2X7wZNc/4t1bD94kkXzS+O2D0qnV6xOs1idZrE6zROUIbkI2BXj+C7gwn6c7wvAn5VSMdP3ExHUdsNFORn85obZjM1M4zc3zKYoJyNu+2TZFtnqNIvVaRar0zxB1doaC9TGOH4CGCkiGUopLe+ziKQA1wE3J2h3A3ADQGFhIatWrQJg0qRJdHZ2dv+dl5fHtGnTWLNmDQCpqanMnTuX9evXc+LECQAqKio4dOgQ+/btA6C0tJT09HQ2bdoEQH5+PpMnT6aqqgqA9PR0KisrWbduHc3NzQDMmjWL/fv309DQwPvGwLDWJg63n2TLli0AjBs3jokTJ1JdXQ1ARkYGIkJNTQ2trU7XVFZWsmvXLg4ePAjA1KlT6ezsZPv27QAUFRUxfvx4ampqAMjMzKSiooLq6mra2toAmDt3Ljt27ODw4cMATJ8+nba2Nurr6wEoLi6moKCguzZZdnY25eXlVFVVEYlEAJg3bx6bN2/m6NGjAHR2drJ3797uzNySkhLGjBnTHYWSm5tLWVkZq1evRimFiHDppZdSV1dHY2MjAOXl5Rw7dozdu3d3j1NWVhZ1dXXGxun06dNEIhHtcQKYMmUKKSkpccdp1qxZRsepubmZSCRifJzKyso4efKksXFqbm5m7969xsfJ6+8p0Tg1NzdTU1NjfJzA7O+pubmZVatWGR8n6P/vqS+C8pG0A88ppT7Y6/ivgE8CIz0Ykg8BjwITlFIRndeExUfihYamVpaurGfR/NKEM5dkYSh+JotlKBM2H8kRICvG8WzglK4RifIF4Oe6RiQWrlUOK250186Gt7Siu4JGpz+9Rqz5QdjH3cXqNIvVaZ6gDMlGoCTG8YnAa7onEZGzgfcBjw9EjDv9CytLV9YzZVwWN12YwpRxWSxdWR+0pLjo9Kf7mR7/dEVgnyns4+5idZrF6jRPUIbkWWCCiJS4B0SkALgA+M+eDUWkQET60vlZYGV/nexeCSp50I3u+sH601rRXUHS0NTKzze1JewjrxFrFoslvARlSJ7AmXksEZHUqKF4ACdq61G3kYhcgpNk+EjvE4iIAJ/DyUcZEOXl5QnbBLkU40Z3TSo8Syu6KyjcPho2cnTCPvIaseYHOuMeBqxOs1id5gnEkCil2oErgU6c3JGtOP6R+Uqp5h5Nm4HjwJsxTvNeYCTwx4HqOXbsWMI2QS/FFOVk8OXZeaE1IvC3Pvrm/HFafVSUk8F3PnZRYJ9JZ9zDgNVpFqvTPIHV2lJKHVJKXauUmqyUmqKUukopta9Xmzql1Bil1D0xXr9SKVU4ECe7ixsSF48wLMXo6AwSt4++9uy2pFiuCnt/ulidZrE6zWOr/2oShqWYsOP2UXaa2D6yWM4ggkpIDBWTJk3SaucuxQSFrs4gKcrJ4NsfuiApjEgy9CdYnaaxOs1jZyRAVlaslJbwYXWaxeo0i9VplmTRCdaQAMmT+GN1msXqNIvVaZZk0QnWkAxJ7GZZFotlMLGGBKdgWTKgozMMpUeGUn+GAavTLFaneawhAaZNmxa0BC10dAad7wLJ0Z8NTa08tUMlxawtGfoTrE7TJItOsIYEoLt0ctjR0RmGfJew96c7a9u+50BSFMEMe3+6WJ1mSRadYA3JkMPmuyTGnbV9pXxEUhTBtFjCjjUkOJu4mMQvZ7euzqBLj5juT9O4s7YfvdqWFBn4Ye9PF6vTLMmiEwLa2Cpo/NzYyl02mTIui+0HT9pZQUixm2pZLN4J28ZWocLdstIEfjq7Ter0k2TQWZSTwdUlkaQwIsnQn2B1miZZdII1JADdexKbwE9nt0mdfmJ1msXqNIvVaR5rSAxjnd0Wi+VMw/pIgObmZjIzMwNUpIfVaRar0yxWp1nCqNP6SOJw6NChwN7bS4RXkDq9YHWaxeo0i9VpHmtIgH379iVu5IHaPY1c/vAqavc0xm3npZxJQ1Mr96x4PfTJc2C+P/3C6jSL1WmWZNEJ1pAYp3ZPIx9f9goNTa18fNkrcY2JboRXQ1MrVy97hU1HOrk6em6LxWIJC9aQAKWl5iKrvvb0qyjgPaVnoaJ/94VuhNeSFdtobOlg6tlZNLZ0sGTFNmN6/cBkf/qJ1WkWq9MsyaITrCEBID09Xaudjj9jckE2SsGGfU0o5fzdF/oRXk5AxLCUlLf9HUYamlr57uqGpJg16Y570FidZrE6zWMNCbBp06aEbXT9GXd/eBr52em0tEXIz07n7g8PvILn4gUXkDtyOJv2N5E7cjiLF1ww4HP6gdtHuw8eS4piiDrjHgasTrNYneaxhkQTXX9GUU4Gz37pEj48o5Bnv3RJ3DwSXeNUlJPBMzfNYfrYFJ65aU5oc1NsMUSL5czEGhIgPz8/YRsvGeu6RROXrqxnQt5IxmamMSFvZNwLb1FOBovnF4fWiMDf+mjZZpUUxRB1xj0MWJ1msTrNYxMSgUgkolVp03Shv9o9jVzzWDW5o9JobGnn6S9WMnNCbp/tdXUGSUNTKz98YQe3XDE51EYPkqM/weo0jdXZf2xCYhyqqqq02pkuz768dh+zz83jigvymX1uHstr48eN6+oMkqKcDBaMbQy9EYHk6E+wOk1jdZonMEMiIvki8pSIbI8+lovIeA+vLxOR34vIehHZFj3Hg35qNs2i+aXsPXqKI83t7D16KvRLQRaLxRKLQAyJiKQBzwNpwDRgKtACvCQiCYvLiMgcYAXwoFKqXCl1PvBD4OP90WM6zE637InXAo/JEg6YDDobmlr5xdZI6CPLIDn6E6xO0ySLTgjIRyIi1wM/Ac5VSu2MHhsHNAC3K6UeivNaAbYA/96znYgMB65QSq1I9P793dhKx0fi18ZWdiMmc9jNxyyW/hE2H8lVwF7XiAAopQ7iGIirErx2LnA+8N89DyqlOnSMSCx0jIpbpmTNjiNxy5T4sbGVe+F7ff/hwPIzvBSX9Gv3SVO4Y/TFqSRFmHLY+9PF6jRLsuiE4AzJRcCuGMd3ARcmeO2c6L+joz6SzSKyUUTuFZF+3VY2NzcnbNNdpqQwO26ZEj82tlq6sp6C7HT2H2+nIDt90C98XopLgl5/en1/XSOmgztG9605khRhyqb70y+sTrMki06AoGLLxgK1MY6fAEaKSIZSqq+rRnH0398A1yil/iIiZcD/AO8G/i7Wi0TkBuAGgMLCQlatWgXApEmTiEQi3X/n5eUxbdo01qxZA0Bqaipz587lWONROjs7OXLkCApFc/PJ7teUlpaSnp7enYl6z2V5rNir+MqFivoNNexNT6eyspJ169Z1fzlmzZrF/v372fj6Pv7wRgdfvfJ8inIz2LJlCwDjxo1j4sSJVFdXM7a9g9/sbkeAN1samZ3TwqpVq6isrGTXrl0cPHgQgKlTp9LZ2cn27dsBKCoqYvz48dTU1ACQmZlJRUUF1dXVtLW1ATB37lx27NjB4cOHAZg+fTptbW3U1zvGqri4mEfXHScvtY1PniP8uKmTpSvr+fv8E0QiEQDmzZvH5s2bOXr0KOCELe7du5edO50JZ0lJCWPGjOneOjQ3N5eysjJWr16NUgoR4dJLL6Wuro7GRqfIZXl5OceOHaN2604e+MtppozL4pplL/O1smHkZQzrc5zWr1/fvbNcRUUFhw4d6q6i6o5T/aZNfOVCxe92dPGDT1dQv6GGepw16b7GqaGhAYApU6aQkpISc5wAMjIymDVrFjU1NbS2Ol/hgY5Tc3MzkUgk4TgVFBR038VmZ2dTXl5OVVVVn+NUVlbGyZMnjYzT7t27aW5uZu/evWRlZVFXVxf396QzTu7vKT8/n8mTJ3dHMQ10nJqbm6mpqTE+Trq/J91xam5uZtWqVcbHCZzrXn/GqS+C8pG0A88ppT7Y6/ivgE8CI/syJCLyU+DzwL8ppe7tcfzLOA73y5RSq+O9f28fSWtrKxkZ8SczDU2tfPSRl2lui5CZnsrvbo6fta6D7lr9pQ++xL7GU4zLTufgiTaKc0ey+rb3Dui9/dDpotOfutzx7EaONLfz+KcruP4X6xibmcZ3PnaRkXOb1OknVqdZrM7+EzYfyREgK8bxbOBUnNkIwMnovxt6HXfL7F7sVcz+/fu12qWmCLkj00hNkbjtdJdidP0ppQWZKAWn2yMo5fw9mHiNLtPtTx38WCp0ManTT6xOs1id5gnKkGwESmIcnwi8luC1rnOit/bOPo4nxJ0Ox2PpynqmFY7m5dvnM61wdNy9Q3T9CboXybs/PJ2zstJpae/krKx07v7wdL0PFhA6/amLVyPmBZM6/cTqNIvVaZ6gDMmzwAQRKXEPiEgBcAHwnz0bikiBiPTU+Scco9F7fcO9uv7VuFr0L/peorZ0L5JFORn87uZLmH22mSU1r3h1tns5r26+jcmKAhaLxSxBGZIncGYeS0QkNWooHsCJ2nrUbSQilwAHgEfcY0qpfTi+kJtFpDTargi4DXheKfWSVzFTpkxJ2Eb3ou/XUkxRTgYPXKV3MTUd5eQ1pFmnP/0yTl7Q0RkGrE6zWJ3mCcSQKKXagStxZhZbgK04/pH5SqmeMW/NwHHgzV6n+DrwI+BPIrINWAP8Dvhwf/SkdG8YFR+dO2MvSzFeL6Y6Ov24QHs1jjo6/ci38YruuAeN1WkWq9M8gdXaUkodUkpdq5SarJSaopS6Kjrb6NmmTik1Ril1T6/jnUqp+5VSpUqp85VS5yqlbkvgpO8TN0QwEaaXYrxeTHV0Ll1ZzznR0vTnJChNr4tXP4WOTj+d6LrojnvQWJ1msTrNY6v/ahKGO30dFs4sZu0bR3lh62HWvnGUhTOLE79IA9N+Cj+d6BaLZXCxhgQnWSkRfizFeL2Y6uhcXruPymhp+kqN0vR+oKMTgnei6+oMGqvTLFaneezGVkBbW1vCSpthKPQ3lHSGAavTLFanWcKoM2wJiaHCLW8Rj6CXYhqaWrnh8VVa/pnvXTODnW81871rZgRyt6/Tn37hJWItSJ1esDrNYnWaxxoSDwS1FOPOMk60q4T+mYamVr729AYmnZXJ157ekLCtyTDhoAlDSLHFciZiDQmErp5Nb1z/zO2X5Cb0z+j6cvy86AbVn179WGEfdxer0yxWp3msj8QDfmwuVbunkduW1/HgwjJmTsjt8311/R66bf0shhgUYfAPWSxDGesjiYNbFjoeftzB1+5p5JrHqjlxOsI1j1VTu6cxZjvX77F571sJ/R5BZ+A3NLXy+WUvBLKs5NWPpTPuYcDqNIvVaR5rSKB7P4J4+BH+e9vyOnJHpfHXu64gd1Qaty2vi9nO9XsUjJSEfg/Qz8D34pTX8ae4xvbYqUhgPgovfiydcQ8DVqdZrE7zWEOiiR938A8uLKOxpZ2L73uBxpZ2HlxYFrOda8S+Uj7C6Pa9XpzyOrMxP3R6ZagFEFgsyYD1kaAfr+2Hj+R/Nx3k68vreGhhGe+fHjsByb2Ql+aPov5wi5G1fy8+kjue3ci+xlaKczO6/43V1g+dXmhoauXqR19BRFBK8cxNc+K+fxjj9GNhdZrF6uw/1kcSh127Ym0f/05Mh/82NLVy/5+2MntSHvf/aWufd9Hu2n9a12ljF2cvMyzdsit+6PTCkhVbaTzVwdTCbBpPdbBkxda47XXHPWisTrNYneaxhgS692gebLzuXXLNpK5Aal15KbtiWqc3eu9cGX8ny6DG3StWp1msTvNYQxIgfkZOmaxSvGh+KXuOnuJIczt7jp4KpFKvDosXnE/uqOFsOXCC3FHDWbzg/EF9f+ufsZypWEMCTJ06NZD39RquqqPTjzBlP3T6QVFOBs/cOId5k8fyzI3x/SNgVqefCZ5B9adXrE6zJItOsIYEgM7OzsSN0L/j9HJn6sXvoqPTrw2jTOv0Cz906oynnxt1BdmfXrA6zZIsOsEaEgC2b9+esI3uHWdDUytXL3uFNTuOcPWyV4zemeroDMOGUTo6w4DJcfez34dSf4YBq9M81pBoonvHuWTFNhpbopFDLR0sWbFtUHWGofrvUGLpynomRHecnBBnx8mgq0NbLEFiDQlQVFSUsI3+HWfvvBxzeTo6Or0kGvqFjs4woKNz4cxiqqOhz9UGd5z0wlDqzzBgdZrHGhJg/PjxCdvo3nEuXnABWRmpvPz6EbIyUlm84AIjGlvW1qBu+CIta+PX3/Fjz3av6PSnF/yKhtLRubx2H7Ojoc+z44Q+++lsN92ffmF1miVZdII1JIB+cTRdR26qCLkj00iV+HkMurSsrWHfjTfS+eab7LvxxrjGxK89272g259e6nf5cYHW0blofil7o6HPe+OEPvvpbE+W4n1Wp1mSRSdYQ2KcpSvreX9HA0/8+T7e39GQ8ILSsraG+vmX92kcXCOiTp8GQJ0+HdeYhGHPdh281u/y4wKtQ9DVlC2WZMAaEiAzM9PYuW4c3cSHnn6YyIEDfOjph7lxdFOfbV0jETlwIKZx6G1EXOIZkzAkD+r0p66B8PMCrTvuutWU/XK2m/x++onVaZZk0Qm2aKNRYl34ZcQIipctY9TsWZ7b1s+/nMiBA32+X2phIaUrX3zHcT+KS5rG62ZdYf88FsuZgC3aGIfq6mqtdvHW9L3MHnTbFt5/PzJiREwtMmIEhfffH/O5oPaWd9HpTy9hyn58noamVj73aDAbcHlF9/sZNFanWZJFJ2gYEhFZGX08YPKNRSRfRJ4Ske3Rx3IR0QpTEJHdIrIhxuOK/mhpa2tL2CbRmv6BO+98h2FwUadPc+DOOz23HTV7FsXLlr3DmPQ1ywkLuv0ZVJiyO5aNrcFtwOUFnf4MA1anWZJFJ+jNSIqBu4GnTL2piKQBzwNpwDRgKtACvCQiWguDSqkZMR4vmNLYm6Ur67ni9D5ue2IxV5ze9441fS+zBy9texuTsBsRXfxyonspZxLkBlwWy1AioY9ERF5VSr3L6JuKXA/8BDhXKbUzemwc0ADcrpR6KMHrdyulSvr7/r19JJFIhNTU1Liv2fPiGhpv+TLpne20paSR+8MfMeHyeW9rY9pH0rv9gTvuoPA73wm9EdHpTy8+El10N7Zy33tyQSY7DjWHPhNdpz/DgNVpljDqHIiPRNsbLyL/qdn0KmCva0QAlFIHgS3R5waVHTt2xH2+ZW0Nrbd+hfTOdgDSO9tpvfUr74ia8jJ78DrTGDV7Fh2PLDVqRGr3NHL5w6uo3dNo7JyQuD/Bnygn3Y2t3PdOjbSG3oiAXn+GAavTLMmiE/QMSZ6IfEpEPp3oAczUfN+LgFjbf+0CLtQ5gYg8KCKviMgOEfmziHxI873fweHDh/t8zmsIrmsgUgsLEy5BeWmbSKdXavc0cs1j1Zw4HeGax6qNGhNdnead6PobWxXlZPCP5wa1AZc3TI67n1idZkkWnQA686Zi4AkSbTfnoDt7GQvUxjh+AhgpIhlKqXge0MPAeuB2IAW4Afi9iHxZKbVUU4MWOo7x3iG4o2bPihmWGwsvbU1y2/I6ckel8de7ruDi+17gtuV1vHjrZYOuwySLF5zPuj3HtDa2amhq5eeb2iid0ZoUxsRiCTM6PpJ64As65wIeV0olzBgTkXbgOaXUB3sd/xXwSWBkAkMS65z/A7wHyFdKvePKLyI34BgcCgsLZz71lBM7MGnSJDo7O9mzZw8AeXl5TJs2jTVr1gAw4vXXyVn6SExjotLSaLz5S5T8wz+Qnp7Opk2bAMjPz2fy5MlUVVUBkJ6eTmVlJevWraO5uRmAWbNmsX//fhoaGgCYMmUKKSkpbNmyBYBx48YxceLE7hDAjIwMzj33XN544w1aW52uqaysZNeuXd1bck6dOpXOzs7u8tNFRUWMHz++u9RCZmYmFRUVVFdXs/lgC/fXtJGWOoz2zi7ufHc65+WmMn36dNra2qivdxzQxcXFFBQU4PqUsrOzKS8vp6qqikgkAsC8efPYvHkzR48eBWDChAmkpKSwc6ezcllSUsKYMWNYv349ALm5uZSVlbF69WqUUogIl156KXV1dTQ2OjOj8vJyjh07xu7du7vHKSsri7q6upjjlJqayty5c/nNn9fyw5pjfH56Gv945WwOHTrEvn1Odn9paSnp6ems/utG7q85zbBhwvDhadw6Q8jLGGZsnGbNmkVNTY2RcWprayMSiXDZZZexY8eO7rtUE+NUVlbGyZMnjY1TJBJh8uTJWuO0fv16Tpw4AUBFRUXMcfLr9xSJRMjKyjI+TgBz5841Nk6HDh0iNTXV+Dh5+T31HqesrKyYPhIdQ7JeKVUet9Hf2i5VSi3SaHcA2KGUuqzX8T8AlyulRum8X6/XfgMnuqxCKRVrttNNb2d7Q0ND3EqbXh3jfpFIp6dzNbXy0R+/TPPpCJkjUvndly4xdmeuq1M30dBLOx0H/i2/eZXntxyiYvwo1u1v4cqpBfzwE0bjSYxictz9xOo0Sxh1DkpCoo4RibIRKIlxfCLwWrwXikhGHyHC7nZiKZoaunHvFvoiLCG4iXR6YenKesrG57DlnvdTNj7HaAisjk4vG4XpFm3UDyl2bp6aW1re9ndYMTnufmJ1miVZdIK+s/3TInK5wfd9FpggIiXuAREpAC4A3hb5JSIFItJT5zXAwzHOORNow4n8Mo5Xx3jYCbrIoO5FX3djKdD/TIsXXEDuyOHsPdFF7sjhxkr9WyxnKjqGZCXwXpxIK1M8gTPzWCIiqVFD8QBO1NajbiMRuQQ4ADzS6/WfEJGLe7S7BvgI8KBSqtmrmOJivVLrrmM8KCOiq1MHv4oMNjS18szu1ITZ4roXfS8bS+l+pqKcDJ65aQ6zS7L7zDUJEybH3U+sTrMki07QMyRPK6U+q5T6XqKGIvJ+nTdVSrUDV+IsR20BtgLZwPxehqAZOA682ePYCuAh4MfRsii7caK3blRKfUPn/XtTUFDQn5cNOmHX6S5DtXalJlyG0r3o624s1fO8OiHFRTkZLFmYHFsRh33cXaxOsySLTtAzJLErAw6wrVLqkFLqWqXUZKXUFKXUVUqpfb3a1Cmlxiil7un1um8rpS6OlkUpUUq9Syn1Ew8634YflYD9wKROPzaMcperPj3ptFbpEZ2Lvu7GUv1Btz/92qFRlzPx++knVqd5dAzJNBHZqfPAqZllSQK81rrSuZi6y1U/WH/amN/FS5VgP/Bzh0aLZaigY0jagT2aj3Z/ZPpLdnZ20BK0MKlz0fxSNh84ziUPrGTzgeNxL/q6F1OvfhfdrXb9qBLc0NTKUzu6Ep4v6B0a4cz8fvqJ1Wkeo0Ub/Sjw6Ad+bWyVTOgWOAS449mNHGlu5/FPV3D9L9YxNjON73xsYLEXuu/v13t72VTLdHFJiyVZGUgeyWc8vI+XtqHBzZgNO7o6dUupTysazcu3z2da0WgjYbVedOoWWPQjTNmdZVyn4cvxcwtdXYba9zNorE7zJDQkSqk63ZN5aRsm3NIEYUdHp+4ylJcLtNeLqV5/6hVY9ONC7n72h//SrGWcgt5xcih9P8OA1Wkercx2ERkmIheIyJQex84XkWdEZLOIvCgig17+3fJOdNf0vV6gTV9MFy84n9xRw7UKLJp+b/ezZ6eJXaqyWAyg4yM5D/gjMDl66A/A54BNQAHwFk4132HAR5RSf/RNrSF6+0i6uroYNiz829fr6AzDmr5uf+rW0PKLoaYzaKxOs4RR50B8JA/gJAv+CFgGzMPJPl8N5CqlzgYygaXALcYUDyKbN28OWoIWOjrDsKav259BLxnp6AxD+O9Q+n6GAavTPDqG5N3AJUqpryqlbgYWAAuBm5VSJwGUUm3AvwDn+abUR9yy2mFHV6cfF2gvSXlB9qdpnV7Cf/1KXBxq38+gsTrNo2NITiuldrt/KKX+AuxXSr1tSz2lVAfOxlSWIUYY7sp18KLT3djKVE0wr30UdLa8xWISHUMSqwjisT7advZxPNSUlZUFLUGLoHR6TcrT1Wn6Yqqr073oS0a2sZpgXmcuXoyO/X6axeo0j44hSRORYhE5x33EOuYe91mvL5w8eTJoCVoEpdNrLoeOTj9mOYvml7K5IZqt39B3tr570b/nyiJj2epe+sirYbbfT7NYnebRMSRTgd04Jd7dR6xju3D2E0k63C0sw05QOr068HV0+lZ6RHr9GwP3ov+V5VsSXvT9KA/j1TDb76dZrE7z6BiSQ8A9vR53xzj27WhbyxDEtAPfr4z1aYXRbP3CvrP1veSReDF4XkrYBx1ZZ7GYRCeP5K9KqYvjNupH2yDpnUeye/duSkpKghOkyVDTWbunkduW1/HgwjJmTsgd8Pt6zaHR0RmGvJyhNu5BY3X2n37nkXg0DD/wpCokjBkzJmgJWgwlnV6q+uo65b3e6evoDMPsYSiNexiwOs1jOm3ynw2fb1BYv3590BK0GEo6vUZY6TrlvSzB6fZn0ImTQ2ncw4DVaR7dWlvpIjJPRD4iIkUxnq8UkT8AyROvZgkUXR9JGPYDsVgs8UloSESkBNgAvAT8J1AvIv8Qfe5yEVkNVAEX4zjdk47c3IGvzw8GQ0mn7pKRH055LzrDgK7OoJMch1p/Bk2y6AQ9Z/vTwPnAz4HhwPVAF/AIjk9kffTf/4hmt4ceu7GVd4IsXBh00cQg0f3sYQgKsAx9BlK08WLgcqXU95VSDwEfwKkE/EngCqVUhVLql8liRGKxevXqoCVooavT9J2p19Ij1z3ynNG7Yr98FGEfd7fft+8+kLDfw7AEGPb+dLE6zaNba+uI+4dS6nXgKPAPSqmVvikbRBLNysKCjk4/Msa9OsZPtKnAanJ5MaJhH3e3328pH5HQOPi5BKhL2PvTxeo0j44haYtxrKGncXERkfsHLmnwEYmTBh0idHT6cWfq1TH+1ZkZgdwVezWiYR93t99/uP608V0s/SDs/elidZpHx0eyXilVnuhYvONhYyj7SPxaK9dZqw96nf6OZzdypLmdxz9dwfW/WMfYzDS+87GLBnxe6x+yWBwG4iOZISKdPR+xjkWPJ2X4b11dcmw1r6PTrztTHT+F+97S1hzIXbHX5R2d/gy6hH5RTgb/eK5KCiMylH5HYSBZdIKeIWkEftHr8WSMY7+Mth1UROReEVEi8pn+nqOxcdBl9wtdnTnbNvK5pV8jZ9tGnxW9k6KcDK4tJZALn1cjqtOfYXBiD7XvZ9BYneZJ1WizVyn1WZ2Ticirum8sIvnA9wB3mvQa8FWl1H4P5xhPkmbT+0XL2hr23Xgj6vRp9t14I8XLljFq9qygZQ0a7szJFIvml/KJn6ztnuV864bZxs5tsQwVdGYkf9f7gIiMFZE8nbaxEJE04Hmc/Uum4ZSlbwFeEpFMnXNEuR8YcORYeXno3TpAYp09jQjQbUxa1tYMhjzAWQr648GspNj5T2fcw+DEHirfz7BgdZpHp2jjWwAiMkNEnhGR4zjl4g+LSJOI/FpEpvVsq8F1wEXAYqVURCnVCSwGJgE36ZxARGYClwA/0nzPPjl2rK8NH8NFPJ29jYjLYBqThqZWrl72ClWvH+XqZa+E3pjojnvQtbZ0dQad2T4UfkdhIll0gn6trU8BNcAVwDrg6eijFlgArBORqz2871U4S2bdO7copQ4CW6LP6fAwcBexw5M9sXv37oGeYlCIp/PAnXe+w4i4qNOnOXDnnT6p+htLVmyjsaWDwpFdNLZ0sGTFNt/fcyDojnvQF2gdnV6TRv34PEPhdxQmkkUn6NXamoFz1/8lIF8pdblS6tro43IgH7gFWCYi0zXf9yKcHRV7swu4UEPTR4AMHGNmAQrvvx8ZMSLmczJiBIX3D0aKT+9Q8uRJqOqLoKO2Gppa+fmmtoTv6yVp9OpHX2HNjiNc/Wj4Z42W5EAnj+TXwB+UUv+RoN0/Ah9USn0y4ZuKtAPPKaU+2Ov4r3BKr4xUSsX8hovIcGAT8HmlVJWIXIZTUPKzSqkn4rznDcANAIWFhTOfeuopACZNmkRLSwuHDjmbO+bl5TFt2jTWrFkDQGpqKnPnzmX9+vWcOHECgIqKCg4dOsS+ffsAKC0tJT09nU2bNgGQn5/P5MmTqaqqAiA9PZ3KykrWrVtHc3MzALNmzWL//v00NDQAMGXKFFJSUtiyZQsA48aNY+LEiVRXVwOQkZHB2WefzZtvvklrq9M1lZWV7Nq1i4MHDzo6Tp/mxOLboe1vkzQZMYKjN91Ix5QpZGZmUlFRQXV1NW3RNnPnzmXHjh0cPnwYgOnTp9PW1kZ9vXMhKi4upqCggOf+7y/84Y0Ori0bw/veczFVVVVEIhEA5s2bx+bNm9mx/y3uW9tKamoKyDAWz0wlL2MYJSUljBkzprssdm5uLmVlZaxevRqlFCLCpZdeSl1dXXekSnl5OceOHeu+K5s0aRJZWVndIZEmxmnEiBFUVFTEHaf/OTyanYeayJR2jrR2cd7Zudz1d5PijtOsWbOoqanpc5ymTp1KZ2cn27dvB6CoqIjx48dTU+MsP7rj9McXq7in6gSFo+BY5wi+/d6xcOpYzHFKyx3HbSsaOCutg/0nu7h3/lkxx+kzj61i7b5Wpo1NYWuj4pKSLD51rvP8QMepvb2d888/3/g4mf49tbe3M3r0aGPj1J/fk5vHlp2dTXl5eczf05tvvklaWhplZWWcPHmye+vdIH9PWVlZMfNIdAzJa0opnVmCAK8ppRLOSgZoSG4BLlNKfSz692VoGJKe9E5IbGxsTIpKmzo6e/pKZMQII1FbXhING5paeXjFZm5dMC30uQ86/Vm7p5FrHqsmd1QajS3tPP3FSiO7OergJlg++MFzue2PbyRMsNRJXrzlN6/y/JZDzC0dS1X9Ea6cWsAPP/EuI3qH0u8oDIRR50ASErWKMSrHIrVr6jkCZMU4ng2cimNEcoA7cBzzxkiWxB8dnaNmz6J42TJSCwuNhf563bf8Q2c3h96IgF5/Lq/dx+xz87jignxmn5vH8tp9g6DMwU2w/NzPa4zVz1q84HxyRw1ny4ET5I4azuIF5xtQ6jCUfkdhIFl0gp4hSdM5UXRGotUW2AiUxDg+ESefpC9mAxHgGRHZICIbgJ9Gn7sneuwbmhqGLKNmz6J05YvG8ke8ZIzrruknC4vml7L36CmONLez9+ipQS2G6IYeZ6dJwtBjXV9OUU4Gz9w4h3mTx/LMjXOSwuBbwo+OIdkgIgs12i0kvhHoybPAhOimWQCISAFwAc7mWfQ8LiLDAJRS/6uUKlZKzXAfwBeiTb8RPeZ5c628vFgpMeEjKJ26uRTuxew0wwOr/usFnf4MOo+kKCeDf7m0MOH7Ll1Zzzl5IxmbmcY5eSMTzhr9CGe2vyOzJItO0PORlOH4IL4KPBXN+ej5fArwT8D/A96rlEpYlyOakLgO2IrjE+kCfgbMBd6llGqOtrsEWAP8RCkVM7/EhI+kq6uLYcNMb19vnrDrdNf0H/uncr74q/XGiib6Rdj700VHZ5C+HJeh1J9hIIw6++0jUUrVATcDP8FJQvyziPwq+vgzcBhYBtysY0Si52wHrgQ6cXJHtuL4R+a7RiRKM3AceDPGB8rvY2nrHR8yEW6kQtgJu053Ceyq7/85sD0xvBD2/nTR0bm8dh+VUV9O5SD7clyGUn+GgWTRCXq1tlBK/UZEtgD/hlMGxS1j0gI8B3w7anC0UUodAq5N0KYOGNPHc4eBGV7e0+Iv7jLQnb9aw9LP2a1eBxO3JtiI4SnsOXqKB64K70zQMvTQnjcppeqUUguB0UBB9DFaKbXQqxEJG6mpWvY0cJJBZ1FOBtfPGKVlRFrW1lA///KE5Vv8ysROhv4EPZ1B+3JgaPVnGEgWnaDhIxmKDOWNrfzC9AZLuvkuQW+WZbFY/sZA8kiGPG6GaNgJSqfXMiGJdHqpUuznfiBn6rj7NcM7U/vTL5JFJ1hDAtBdAiDsBKXT68U8nk6vVYq97nrohTNx3N0KzWt2HDFeoflM7E8/SRadYA3JGY3unanJhESvVYrDsPY/lHArNE8tzE6KCs2W5MAaEpxiZMmASZ1elqu8JiSmjBzd5zn7U6XYrwS6M3HcQaFQbDlwAoXCZIXmM7M//SNZdII1JADdlX/DjkmdXperdC7m7jnvuvSsPs/p1gLrbUxMFZj0gulx98v3YFLndXMm0hHp4mhLGx2RLq6bM9HYuc/E35GfJItOsIYEoLt8ddgxqdMP34N7zn/5fX3cc/Y2JkEYETDbn37u82FS5/LafVSeN5aPvquIyvPGJkxc9GIcz8TfkZ8ki06whuSMpSgng+9dM4OdbzXzvWtmGFk28lJk0I8qxV4wXVxyyYqtNJ6K+h5OdbBkxVYj5/WCzkXfSxHKoDf1siQP1pDgbKSTDJjU2dDUytee3sCkszL52tMbEl4kdJMHi3IyuP+j07UMk+kqxbq4F8iu4SMNXiCFLuX4HrqUAsTAOR10xt1L9V/d4AUvy58NTa08u29EUhibM/H37jfWkODsuJYMmNTp5SLhhuxGDhzoM9/DL51+4H727111gbHclOvmlBDpUhxtbiPSpbhuTsnAhUbR6U+ve8boBC/oLn+6Rux4W1dSzFzC/v10SRadYA0JQPeWnmHHpE7di4SX5EE/dHpFZ+bkh39oee0+ys/JpSg3g/Jzco0WTdTpTz8+k+7sxTVinz2v3XjSqB+cib93v7GG5AxF5yLhNXkwaHRnTl58ObosnFnMq3sbOXE6wqt7G1k4s3jA5/SCX/k2OrMX14j9YP3ppKj6bDGPNSRAfn5+0BK0MK0z0UXCa/KgSxD96XXmVJSTweL5xcYuuH5uyavbn37l2+i8729umE3B6JGhTxptaGrlP94YFvrlN0ie6xJYQwLA5MmTg5agxWDr7E/yIAy+zv7OnEzq9Lolr5ewWl2dfuSx6J6zKCeDpdeFe+te15cTSc1ICl9OslyXwBoSAKqqqoKWoMVg6+xv8uBg6+zvzMmkTi9LSw1NrXz0kZf5/YYDfPSRlxNe0HR0+hGq6/WcYf8dub6cf5rQkhS+nLD3Z0+sIbHEJSzJg/Hoz8zJdB4J6C8tffP3m3jrZBuj0lN562Qb3/z9wJ2qflRJ9rPychBYX45/WENC8oTZBaXTa/LgYOv0OnNy77RbIhLIEkf9oWZEYEZxDiLO3/HQ6c9F80vZfOA4lzywks0HjhutVKAbCRb235E7a8zNSA29LwfC3589sRtbWYYMuptl3fHsRo40t/P4pyu4/hfrGJuZxnc+Nnhb09buaeTjj1UzfJjQ0aX47RcrmTkhd0DndEu0iAhKKZ65yYy/wvSGZpbkxm5sFYdkMSpWZ3x0Z06L5peiXl3HKxVzUK+uG/QljpkTcvn2h6fTqRTf/vD0hEZEpz+XrqxnUn4m8yaPZVJ+5qAvQzU0tfKFx1aG3oEN9nfkB9aQAM3N8ZcWwoLVmRidsis52zZy60uPkdvcyK0vPUbOto2DqNCZkXzj95vIGZnGN36/ido9jX22bWhq5Udrjya8QC+cWczaN47ywtbDrH3jqJE8Fl1nu7tZ1qtvnja+WZZpavc08uUVh+P2eVhIlt87WENiSQJMhrW6y1+4UV4BJFfetryO3FFp/PWuK8gdlcZty+titnMv5CfaVUJfzvLafVRG81gqDeWx6Drb3c2yzskeFurNsmr3NHLNY9Wc6lBc81h1UhiTZMEaEmDWrPBEIMXjTNRpMqw1LJn6Dy4so7GlnYvve4HGlnYeXFgWs517If/1TZcmjJpaNL+UPdE8lj0aeSw6LJpfyuaGqAO/IZ4D3/GzjsnNfdvffeHXvi2JcA342jvmxzXgYSFZfu9gDQkA+/fvD1qCFmeiTpMhqP3NNzHNzAm5PP3FSrJHpPJ0HEe7GzX1+Z/XJIya8m1LYun1bwwWL7iA3JHD2bS/idyRw1m84II+23q5MTBtcFwDfsmDq+MacD/x8pmS5fcO1pAA0NDQELQELc5EnSaLEfY3U98PZk7I5cVbL4vraHeNw/CuNi3jYLpEytKV9UwrHM3Lt89nWuHoPo14UU4Gz9w0h/NzSRgtpntj4EeCpWvA04d1xTXgXtE1Dl4/U7L83iFAQyIi+SLylIhsjz6Wi8h4jddlisjtIlIlIrUiskVENonIP4uIuU0gLKHA5J12fzP1g1qKAefzf3Z6eiCht37uopnonF5norpjNHNCLt95z0gtI6JzTi/GYaglePYkEEMiImnA80AaMA2YCrQAL4lIZoKXlwD3Ag8opWYqpaYCdwLfBe7qj54pU6b052WDzpmq0+SdttdMfb92CfRinEz3p5f6WTpG3O0jlTYqYR/pntOLEXOjxtbsOJIwaqyhqZXf7U+8AZfuuHsxDl4Nc7L83iG4Gcl1wEXAYqVURCnVCSwGJgE3JXhtC7BMKfXf7gGl1B+AjcDH+iMmJSWlPy8bdKxOM7jGZFhBQcJMfT/uIr0aJ5P96YdhdPvoux+ZrNVHOjcGXmaibtTY1MLsuFFj7mdvbO00NnvwYhy8zq7D/jvqSVCG5Cpgr1Jqp3tAKXUQ2BJ9rk+UUruUUotiPJUNvNUfMVu2bOnPywYdq9Mco2bP4s27v5Ww3IuXC0XtnkYuf3hVwrBSr8bJZH963T5Xx+i4ffTFX60PqIZV7yix2FFj7mf/XGmHsdmDV+PgZXadDL8jl6AMyUXArhjHdwEXejmRiIwSkXuBUcAdBrRZQkbQPgqdC0XtnkbuvfeXfOvXd3Lvvb+Ma0z88D3ooh/Sq290inIy+N41MzjY0sX3rplhrDSL7szJjRrbcuBE3KgxL0UbvRgIL8ZB92Yj2Qik1paItAPPKaU+2Ov4r4BPAiOVUgmvGiJSBbwbeAP4glLq5ThtbwBuACgsLJz51FNPATBp0iSOHj3K8ePHAcjLy2PatGmsWbMGgNTUVObOncv69es5ceIEABUVFRw6dIh9+5ykr9LSUtLT07u3xszPz2fy5MndZaDT09OprKxk3bp13dmqs2bNYv/+/d2RGVOmTCElJaX7LmTcuHFMnDiR6upqADIyMhg9ejTHjx+ntdXpmsrKSnbt2sXBgwcBmDp1Kp2dnWzfvh2AoqIixo8fT02Nkx+RmZlJRUUF1dXVtLW1ATB37lx27NjB4cOHAZg+fTptbW3U1zsXjeLiYgoKCrrLNWRnZ1NeXk5VVRWRSASAefPmsXnzZo4ePQrA6NGjycvLY+dOZ8JZUlLCmDFjWL9+PQC5ubmUlZWxevVqlFKICJdeeil1dXU0Njo/sPLyco4dO0bt1p38/o+buOmvT/PTOZ/kA+8/n7yMYUbGKSUlhcrKSiPjdOM/P8ZXVv2U9M4O2lKG84P3foHrFl7U5zg1NLby/ee38aFzh3PRecV9jtOBplZ+t+M093/yPTQf3jfgcXrpL3Xc9r9vAjA8bTg/+vBEWt7aH3OcImlZ3FN1grGpbew72cUdszL42Pvf+45x2rrnIIuefZ3CkYq3OtL46Scv5PCubQP6PX13dQO7Dx7jK+UjWLZZMaFgDAvGNvY5TsXnl/Hwis1cPq6dvIxhfY7Tf698md/tOM0nLsrh7987h5qamkH7Pb3eGOE7f2kjd+Rwjp3q4M53p1M+YUyf4/Tzmjf56OQRzJ/9Lk6ePNnn7ymSlsWLB4dz8cijjBkhcX9Pu3fvBpzrXlZWFnV1dZ7GKSsrK2atraQ2JNHXDAc+AfwE+LpS6keJXtO7aGNbW1tSVNo8E3V+/6FfM//JJQyPtNORmsbK6xbz1a9fa+TcpnS2rK1hzxdvRNr+lqOi0kcw4bGBldt378pL80dRf7jFSH6I14KVtXsauW15HQ8uLOsz0sk959JrLmTR068lPKdOIUj3s08Zl8X2gyeN5sYE9Tu6/OFVnDgd4a93XcHF971A9ohUXrz1sne0cwtwusQLqfazn2IRtqKNR4CsGMezgVO6RgRAKdWhlPoF8CTwXREZ51WMe9cfds40nS1ra3jfUw8xPNIOwPBIO+976iFjGegmdLrZ8j2NCIC0DTxb3l1a+lTJKWOOfq/RUF97egOTzsrka09viOsj2XzgOHO+82LCEva6S1buctnOt5qNLZe5BPU70q1osGTFVhpPdVCQ3kHjqQ6WrNja5znDElIclCHZiBPG25uJwGvxXigiaSKSGuOpOpxwYk8+FotZWtbWMPbOuwZ8sX9HTSyXAGpjxcPPbHnXn3HrqlMJ/Rm6eFn7X7qynnPyRjI2M41z8kbGv0ipXv/GOaduQqKOEQsDXnJYdCoavLOMQN/pcUH623oSlCF5FpggIiXuAREpAC4A/rNnQxEpEJGeOu8Evh7jnO65jnoVk5GRHPsshF2ne/FPOXZswBf7wShnYqI//c6WjyhFc4ciEsAStG5F4aUr65lWNJofL8hjWlHfGfDgLSFR24ih78RuaGrll9siRvOBvIRT61Q0WLzgfHJHDWffSUXuqOEsXnB+n219K43jkaAMyRM4M48lIpIaNRQP4ERtPeo2EpFLgAPAI71ef7OITO3Rbi5O/ska4FWvYpKlOFqYdfYuiDjQQoiDUc7ERH/2N1tehyUrtnKyNcK8KQWcbI3EXeLQxcuFT7eisGscfro9xVg01MKZxVS/cZTfrW+gOkFZfLeq74nTkbhVfd3PPmxkjvEcGpNLS0U5GfzoE+WMznT+HezSOP0hEEOilGoHrgQ6cXJHtuL4R+YrpXoW4W8GjgNv9jj2JPAU8GsRqRORzcCPge8Af6/6ET3gRmGEnbDq9KOqrp8XaBfd/mxZW0P9/Mv7/Bz92ddebznEWdJwI2/iVk7UxGsmtk5FYdc4dJ1qMpZL8eQruxk+TMjLTGf4MOHJV3b32Va3LL/72b8wpdOoz0k3nFoXd1kvN7Uj9Mt6LoHV2lJKHVJKXauUmqyUmqKUukopta9Xmzql1Bil1D09ju1SSi1WSs1QSpUppaYppS5SSt3fywhp44b/hZ2w6vRrGao/F2gv6PSnayQjBw7ENYpe9rV3o3LW7DjC1Y/2XdLDXeLY1diRcInDPW8i4+RXJnZRTgafOj/VYC6F6lV5uO/7Q10ntvvZH3i5UcufoOv7iChF46l2Y8uPrsFbVDbcaJ0xP7HVfy0Dxs9lKC8XaNN4Xa7T2Z0RnCWrot1bePjZb1C0e0ufS1buEkdaCgmXOLxEQ/mVia2D7jLUdXMm0hHp4mhLGx2RLq6bM7HPc+o6sd3Pnp0mCT+7bn+6y4+XnDfW2PKjl8RJv2rBecUaEpxEpGTAtM5ESza6+L0MpXuB9kq8/vRzE6zCXVu4s+pxRp84yp1Vj1O4K3YpjIamVr78m/VIajpf/s16Y5VlvRgHL3e7Ot9P3WWo5bX7qDxvLB99VxGV541NuOPjuNEjePfEMYwbHfuGxqUoJ4OfXH9Zws+u7+wXOpXi5deP0KkUJpYf3dDn413pCUOfvQYl+IU1JMCuXbGqtYQPkzp1l2x08XsZyg/i9adfy3Uta2v48G8fZkRnBwAjOjv48G8fjtn/bjHCCRpb2PoRBur1blfn+/ngwjKOtbRz/r+u4FiCZaidh5tZs+MIOw83G70r19GpG7H2gQvPpj3SxemOTtojXXzgwrMTnjsRro+kIIOEPpKFM4upfv0Iv3u1gerXj8QNSvATa0iguyRC2DGl03SElYtrTDrHjAm9EYH4/enHcp333BjF9MP1fPnndzH9cD3x/AR+hIF6jUjS+X6OGz2CvMw0hg0T8jLT4s8gNHZn9Kqzdk8jn31mZ8Iw4eW1+3jXOblkj0jlXefk9jkjeui5bYwZlcY1FxczZlQaDz0Xf796nRmeO8tIU+0JZxlPvrKL4anDyBuVzvDUYTz5SjA3xdaQnGH4vW/5qNmzOHL/faE3IonwY7nO6yzn1vxm7qp6nNzmRu6qepxb8+PHkpj2Z3jNgv/5praEs4GlK+spG5/DlnveT9n4nD4vkktX1nPuWZnMmzyWc8/KTBhdtvlANHIqTma965851aHi+mfAudN/dW8jJ05HeHVvY593+g8uLONEawcvbD3MidaOuNv36s6c3NnQhsORuLMhB0EQphZmIwgmltb6gzUkOMXZkgETOgcj0W+o9Kfp5Tovs5yWtTWcuvUrpEeXwNI7Ozh161cGNaPf88ZW6ZkJl5Z0jZObR/LC1sMJ80gArcx61z/z/KKL4/pnwJmRzI7m0MyOk0Ojn62uP3Ny83feO+WsuPk78LfIvi0HTmhF9vmFNSRAZ2dn0BK0MKFzMBL9hlJ/mowa053lhKk8jM4sZ+nKeibkjSQ3I4UJCZZidI2T7oXcff9pRdG95eNk1rthwguWvRo3TBiiPpq3oj6at+L7aHSy1d1z6hhRN3/nWEv8/B34W2TfiOHDtJIX/cIaEuguEx12TOgcjES/odafJqPGdGY5gzFrNIk7e1i146jW7EHHOC2aX8reaDLk3gQXU90L9MwJuSy9tpxT7RGWXlue8MIf6Yrmh3TFzw8xvXWx207aE1d8DktNMmtIzkCSMcJqKJFoljMYs0aTLK/dx4xzchg5XJhxTk7CUF2dC6/XZEjdJbj7/7SV88ekcP+ftsZ9/yUrtvXKD4m/fa/pPI6inAw+Oz1dK0zZy46XfiUuWkOCs2FNMmBSp5+Jfmdif3ol3ixnMGaNJlk4s5gNe5to6xrGhr1NcWckXi68XoIHdJfgpozL4sEPnqsRheZt+16TWxe76Hw/dUu0NDS1cvWyaDWFZX1XU+gv1pAA48ePD1qCFqZ1+pXod6b2p0mSadboOocvv6AgoXM4yAS67hIpVccSRqF53b5XJ7LN9SWNzUxL6EsC/e9nZ7RES2ecEi1LVmzjaLOzr8/R5va4OUn9wRoSwlsMsTdWp1nCrjNZ8nJc5/Dr+w8ldA7rJvr5gbsEFmlJXFyyKCeDH10bdWJf27cT28sGXAtnFvNKtKLxKxqfXef7uWTFNk5El+BOxFmCO9UeoT3SRXtnF+2RLk61RxKe2wvWkFgsISYZ8nK81LDSLU0P/qzp6/oedJ3YXpzdT76yi9RoRePUYWIoeVBvCW5kWippKUJLW4S0FGFkWqy9AfuPNSRAZmZm0BK0sDrNcqbq9FJjzUtE0pdn5yW8QOuWpm9oauWbdz/Jhx/4Et+8+8mE7+9lY6tfbe/SSpzU8X14y/4XhomTPDhMEicP6oz74gUXkDUilZdfP0LWiNQ+l+Cum1NCR5eiraOTji7FdXNKEp7bC9aQABUV79jLPpRYnWYZajp1DETL2hr2Rmus7U2Qk+LVOayjUzfC6pnHf8ctzz/K6BNHueX5R3nm8d/1ec7aPY18/LFqGhpb+bjGxlaSkW0scdKLj8RL8mDtnkYWr25OaBgBUlOE3JFppKb0bZiefGU3aSnDGDc6g7SUYXH3d+kP1pAA1dXVQUvQwuo0y1DSqVOE0zUi3YmOp0/HNSZea23p9ufB46f5y65jHDweO1emZW0N73vqIYZHHOfw8Eg773vqoT51/vPTG7jwcD1P/Pk+Ljxczz8/vSHu5/lcaYdWhWQd34fXMGWd5EHXMO472hLXMLqfaVK0lMykuKVkVK/kf7NbN1tDArS1tQUtQQur0yxDRadOEc7+ZMsvml+KenUda989F/Xqurh32y1ra8i69V8SLpfV7mnk3nt/ybd+fSf33vvLd1wk+6Nz/um9fLP6Z4w+cZRvVv+M+af3xnxvd/bwYPVxrdphur4P3TBl3XP+89MbUEoxbWwKSqk+DSM4DvyWtTV84D6nb/py4F83ZyIdnV0cPN5KR2f8/V36gzUkFksSo1uEsz/Z8jnbNnLrS48x+sRRbn3pMXK2bYyrIeXYsYQlXB7/4W/55is/ZWxLI9985ac8/sPfvu15rzpb1tbw8We//7ay/B9/9vsxNXgJCvBjL3bdc5YWZNKlYNORTrqU83df/N9vV3B39c8Y29LI3dU/4/9+uyJmO7u0NQjMnTs3aAlaWJ1mGQo6dS+8XrPl3zEz6GNG4GVLgpa1NdzywrK3FaK85YVlb2vrtbCl19lLUU4GP7vxcq2gAC+Vj3UCEnSrFN902XkME+hSMEycv2PhbQlQIVFHv4hgl7Z8YMeOHUFL0MLqNMtQ0Kl74fWSLa87y/GyJYHbVtre3lba3t7Wi87+1iTTGXddH4nnEikaVYqX1+7jkvPGcuV5WVzSx+6QXo2oboJlf7GGBDh8+HDQErSwOs0yFHR6ufDqZsvrXqC9XMi9tNXV2d+aZDrjruvP8LIEpluleNH8UtI3vcqnH72d9E2vxpy5eDWiRTkZPHPTHOZNHsszN80xXiXYGhKLJcnxUk5Fp8aa7gXay4Xc60VfR6efNcl0DYSXJTDdtq5vKre5sU/fVH+MqOmNz972nipOfZahSkVFhVq3bl3330eOHGHs2LEBKtLD6jTLUNPZsraGA3feSeH99w84Ez7WspXuMlhfF3IvbfurVed8Ov3pLllNGZfF9oMnE1YVXrqynkXzSxNepPe8uIbDd91F/n33MeHyeXE/i4up/vSisy9EpFYp9Y6EITsjYeiEgYYFq9MsujoHe98UL+28tu2PVp1K1i1razhy1VUJw5S95ofo3Om3rK2h9davkNl0hNYYu1163QbbS3/6Ve7exRoSoL5+8CqQDgSr0yxWZ3x0L9Beikv6tX2BjhHtvlAfOqy106TJpSCd6Lb+BA/o9qcf4cw9sYbEYrH0ie4sx0txSb+2L4iHlzBlv9/bpbeG/gYP6PSnl+TS/hCYIRGRfBF5SkS2Rx/LRSRhAX4ROVtE7haRjSKySUS2icizInJhf7UUFw9eKeuBYHWaxeo0S1h1el0yMo3uTMPP4AHd5NL+EoghEZE04HkgDZgGTAVagJdEJFHJy28CnwD+Xik1HZgBdAI1/TUmBQUF/XnZoGN1msXqNEtYdfY338QUXmYafviRdJNLB0JQM5LrgIuAxUqpiFKqE1gMTAJu0nj9g0qpfQBKqdPA7UAGcEN/xPSM4AozVqdZrE6zhFVnf5eMdElUddnrTMPkhmaDNRsLypBcBexVSu10DyilDgJbos/FYxHw772OHYj+m2tMocViGRL0d8lItyx/oqrLsTQkem9TG5oN1mwsKENyERBre7BdQNzlqegMpqvX4cnRf1f1R0x2dnZ/XjboWJ1msTrNEmadXi/kumX5vTjwvUasmehPv2dj3ecKIiFRRNqB55RSH+x1/FfAJ4GRSintQGcReQS4FJiplIoZdC8iNxBd+iosLJz51FNPATBp0iSysrKoq6sDIC8vj2nTprFmzRoAUlNTmTt3LuvXr+fEiROAs4HPoUOH2LfPqYFTWlpKeno6mzZtAiA/P5/JkydTVVUFQHp6OpWVlaxbt47m5mYAZs2axf79+2loaABgypQppKSksGXLFgDGjRvHxIkTu/d4yMjIYNasWdTU1NDa6nRNZWUlu3bt4uDBgwBMnTqVzs5Otm/fDkBRURHjx4/v3vs5MzOTiooKqquru3MT5s6dy44dO7rLRkyfPp22trbukNPi4mIKCgq6ly2ys7MpLy+nqqqKSMTZ93nevHls3ryZo0ePAlBWVsbJkyfZudOZcJaUlDBmzBjWr18PQG5uLmVlZaxevRqlnGJyl156KXV1dTQ2OmXFy8vLOXbsGLt377bjZMfJ2Dite+IJRj/5C07fcD0Vn/lMzHE6uno1uY/8GGl3iiECqLQ0Ur75Dc75wAeoqalh+Pbt72jT3TY9ncYv3UTHlCmhGKfh27eT9+iyt89M0tPJffi7bEtN9TROWVlZMRMSk96QiMjlwNPAPKXUFp3X9M5sr6qqSopKsFanWaxOswwFnX35FODts5j6+ZcTOXAgxhkcUgsLKV35om86veI1+78vwpbZfgTIinE8GzjlwYiUAU8CH9I1IrFw7wTCjtVpFqvTLENBp19l+U3r9IpfiaAuQRmSjUBJjOMTgdd0TiAiFwH/BfyjUuoVY8osFssZix9l+cOCn4mgQS1t3QA8BkxUSu2OHisAGoA7lFIP9WhbALzV08EeNSK/Bz6llKqKHjsb+JZS6ouJ3r/30lZXVxfDhoU/yd/qNIvVaZahorO/hRNNG5Ew9mfYlraewJl5LBGRVBEZBjyAE7X1qNtIRC7BCe19pMexC4EXgeeAEhH5JxH5J+AaYEp/xGzevLmfH2NwsTrNYnWaZajoNF2W3y+dYSIQQ6KUageuxMlI3wJsxfGPzFdKNfdo2gwcB97scexuYCzwReCXPR7f668eNzoi7FidZrE6zTKUdHoxEH4tGSVLfwKkBvXGSqlDwLUJ2tQBY3od+5ifuiwWiwX+ZiAsiQnXAlxAlJWVBS1BC6vTLFanWaxOsySLTrCGBIDHHnssaAlaWJ1msTrNYnWaJVl0gt1qF3CyYN0s4zBjdZrF6jSL1WmWMOoMW9SWxWKxWIYIZ+SMRETeAvb0ODQWJ9s+7FidZrE6zWJ1miWMOicopc7qffCMNCQWi8ViMYdd2rJYLBbLgLCGxGKxWCwDwhqSJENE/k9ElIiUBK0lGRCRs0Xkf0Uk1Gu4yaJzKCIiv4r+pi4LWkuycsYaEhHJF5GnRGR79LFcRMYHrSseInIVEMoNH0SkQkRWiMhWEXlNRP4iIlcHrOljQDVwbpw2Z4vI3SKyUUQ2icg2EXk2WtMtNDp7tC0Tkd+LyPqo1u0i8uAgaJwhIo/3GN8tIvJDETmrV7tMEVka1bVFRP4sItP81udVZ4/2FSSosOEHHvpzsog8Ex3r10Rkg4jcONh6E6KUOuMeQBpQBzyDUyYmBWdfk3ogM2h9cTTXA/8DKKAkaE09tJXg1ET7JZAaPXZjVOcHA9RVA5TiFAlVfbRZBuwAiqN/j4h+L04BF4ZFZ7TdHJwippf0OHYzsHsQNG4D/hMYFf27KHpsB5DRo90KoApnczqAbwNvAUWD1JdaOnu0Xw38d/S7etkgfjcT6gRGA3txitS6/bkA6AIWDZZWrc8TtIBAPjRcH/3iTOpxbBxOEcmvB62vD81fA34NfCuEhuRLUU3v6nX8OPCbAHW5Ri2RIflCr2PnRj/Pj0KkU3CKm3691/HhwIJB0LgNOK/Xsc9H++mq6N9XRv+e36NNGnAMeGSQ+jKhzh7HPwq8DHwmIEOSqD8/EP37o73a1QHVg6VV53GmLm1dBexVSu10DyilDuJUIr4qMFV9ICJjgK8DdwStpQ/crdy6i4CKiOAsnaYEoghQSulsMbcI+Pdex9w9VHPNKoqNps65wPk4d889X9uhlFrhi7C3c5FS6vVex3r301VAB86MxNXXjnOxHqzflY5ORGQ4sAS4dZB09UZH5zt+Vz3+Dux3FYsz1ZBchLP3SW92AYO2Nu6BbwC/UkrtSdgyGP4D5w7rX6Nr5MOAO4F0nDv+0KKUiqgem6ZFmRz9d9Ugy4nHnOi/o6M+ks1Rv869IpLh95tHDUJvJuPcMa+J/n0RcCBG211AgYjk+ygR0NYJzpLgBqXUWr81xUJT58ro/291fSci8ingAmDpYOjUJbAy8gEzFqiNcfwEMFJEMpTmvvF+IyKlwMdxvjyhRCl1QkQuB36Ok4nr7iNzpVJqdaDi+scNwGYcn09YKI7++xvgGqXUX0SkDMdn9m7g7wZTjIik4CzF/EwptSN6eCxwMkbzE9F/84DDgyCvm1g6RSQXWAxcMpha4hFLp1IqIiL/gLOb7AERORZt/nGl1PKApMbkTJ2RJBNLgAeUUseDFtIXIjIF+AtO2ZkxQD5wF/CsiCwIUptXogbxGpwfa1vQenrgbg7+M6XUX6B7v54lwJUicukg6/k3nGWsrw7y+3olls5/w/Hd7Yz5imB4h87oLGQtkAnkK6UKcCLMlonIZwLQ2CdnqiE5AmTFOJ4NnArRbOQ9wHR6bD8cUr4N5ABfUUqdUkp1KaX+A2da/qSIJMXMN3qH/yTwIaXUlqD19MK909/Q6/ir0X8vHiwhIvJZnFnyAqVUS4+n4v2uAAZ1y79YOkXkXJyL8bcHU0s84vTn14GpwM1KqUYApdSLODPlZSJSMOhi+yApfuA+sBHHcdmbiTh7yYeFK3Gcan91fNeAE10G8CcRaQfuVEr9KQhxPbgQ2B/DAO8APoLTr/WDLcoLInIR8F/APyqlXglYTiy2Rf/tffPX2cdxX4iu0d+KE5nVe5lqI1AhImm9fAATgUMx2gehcz5wGnipx2/K3YX1pyLSDCxVSv00YJ3g/K7alFL7eh3fgeN/vAh43n+VGgQdNhbEA2cN/G0htEABTpREKMN/e+j8Vm/tQT9wYvGPEw1j7XH8P3Bi3scGrO8J4udnuMEXc3scOxt4LCw6cXwkEeBfex13Q9nfOwj6/gnYBIzrcewfgBui//87eoXRMsjhvzo6Y7T/TG/dYdCJMztWOMtaPV/3QPR4+WDqjfc4U2ckT+CEfS4RkU/iXOwewLmYhH0ZKYz8CCeJ7x4RuUsppUTkvcDHgKeVUmErhd1NNIP9RZzksJIepWfGAlOC0tUbpdQ+EfkhcLOIPK2UqheRIuA24Hml1Et+vn/0d/I4zlr+FT3u5t8DvBnV+GcReQ74toi8Tyl1CsdX1gnc76c+LzrDgKbOR3GW4R4SkeuVUu3R7+sNwCv8bVkzeIK2ZEE9cGYgv8aZJm7HuZAUB60rjt4P4KyPH8S5G9mCE74YuLaovvcBL+EkzW3CWeb4OpAeoKaHon12LNpnG6KPtB5tno0+F+uxKiw6o+1ScMKq63GWut4AHiRGxrYPGo/F6adv9WiXCTwS/V1txVl6mTaIY66lM9p2RrSf90affz36d2FYdOJE5P13dLxfw4kmvB8YPVh9qvOw+5FYLBaLZUCcqVFbFovFYjGENSQWi8ViGRDWkFgsFotlQFhDYrFYLJYBYQ2JxWKxWAaENSQWi8ViGRDWkFgsFotlQFhDYrFYLJYBYQ2J5YxDRPJFZIOIHBMRFf3/F6LPfVVEPhKwvo+IyFdjHJ8V1fzuAGRZLH1iDYnljEMpdVgpNQP4Q/TvGepv1V6/ilOxOEg+Qux9Plpw9nxpifGcxRIYZ2rRRosl6VBKbQLeFbQOi6U3dkZisQAiUiwiG4BC4EPR5a4NInJFjzY3iMgWEdkuIm+IyP0iMrzH8+5y2W4ReZ+IrBKRhujyWY6IzBCRp3uce72IXNdLx3PAh4DCHu1uj55vQ/Rc3+r1muki8t/R990lIn8WkfIez98Y1a1E5GYR+YmI1EXbL+p1ruzo86+JyKsiUisi3xaRkSb72zLECLpqpH3YR1APYuz/AewGnojR9jagjeieJTj7ldQD/x7jnCeA7wMCjMSp2JwD3A78gui+LcBknF0FPxbjHLv70Ny7Oux5OHvBfBe6i7B+C2gGLuzRriT62o3AxOixG3C2UJjSo91PgRU9NM4EWgnR/jf2Eb6HnZFYLAkQkdHAN4FnlFJVAEqpN4GHgc+IyMReL8kC7lcOp4BKHOPyBHCLUioSPccO4AWczan6y7ei//6bUsot5X0fjh/lvhjtVyqldkX//yyOseu53/tsYF8PjbXAv0b1WywxsT4SiyUxlTgzi5d7Hd/E3y7Eu3ocP6p6bJvqXrhF5ATwLyLy99HzdQLnAG8NQNsVwGbVY5tjpVSHiLyKs2GS9DAw4OwT4nIs+m/Pvb/XADeKSCbwc+AlpdTDA9BnOQOwhsRiSczY6L+3icgXexxPBQ7hzEB60tzHef4duAxnS9dtACLyRPTYQLTVxjh+DMjAMVg9o7xOuf9RSnVFd+ZL6fH8LTgbvd0IfAJ4S0T+H/CgUqprADotQxhrSCyWxLhbBX9TKfWL/pxARDKAhcCPXSNiiCPAmBjHx+D4Nk7FeK5PoktaPwB+EM1XuRP4TvR9fhrvtZYzF+sjsVjeTgfOchUiMkFE5uDsj90ClPVuLCL/LiLTNM47HOfOv/eWpOMSaBglIh+Kc94XgGk9o6pEJBVnG9kXei1rJUREfuaeSyn1F+AqoAm4yMt5LGcW1pBYLG9nFzA++v8bgC8opU7gONu/ICKzAMTh6zh5HQlnGNFz/B/wcREZHz3HHODyPjSMFZF0YA5OBFhf3I1jnL4t0XUqnFlEFnBXIl0xuBzouXx3UfRcL/XjXJYzhaDDxuzDPgb7AeQDG3D8CCr6/y9En6sEtgKvAWuByT1e99no8R3R1/wMyO/x/EvRc7ZHn7+51/sWA/+FEw68Jvr6/+nRvqSHvpei77MJJ6/kfdE2Kvr6P/Q47/ToefbghC8/D8zs8fw1wJboa/fiGJipvc73i2jbz0W1vRZ9fgPw+aDHzD7C/XDjzi0Wi8Vi6Rd2actisVgsA8IaEovFYrEMCGtILBaLxTIgrCGxWCwWy4CwhsRisVgsA8IaEovFYrEMCGtILBaLxTIgrCGxWCwWy4CwhsRisVgsA8IaEovFYrEMiP8P7hJWPmeBz+MAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import json\n",
    "from matplotlib.ticker import MaxNLocator\n",
    "from  matplotlib import rcParams\n",
    "from matplotlib import cycler\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "rcParams['xtick.direction'] = 'in'\n",
    "rcParams['axes.grid'] = True\n",
    "rcParams['grid.linestyle'] = '--'\n",
    "rcParams['markers.fillstyle'] = 'none'\n",
    "rcParams['axes.prop_cycle'] = cycler(\n",
    "    'linestyle', [\"-\", \"--\"]\n",
    ")\n",
    "\n",
    "# enable usetex by setting it to True if LaTeX is installed\n",
    "rcParams['text.usetex'] = False\n",
    "rcParams['font.size'] = 16\n",
    "rcParams['font.family'] = 'serif'\n",
    "\n",
    "logfilename = opt.logdir + \"calibration.log\"\n",
    "with open(logfilename, \"r\") as filename:\n",
    "    log = filename.readlines()\n",
    "    \n",
    "\n",
    "options = json.loads(log[7])\n",
    "\n",
    "goal_function = []\n",
    "batch = 0\n",
    "batch_size = options[\"popsize\"]\n",
    "\n",
    "\n",
    "eval = 0\n",
    "for line in log[9:]:\n",
    "    if line[0] == \"{\":\n",
    "        if not eval % batch_size:\n",
    "            batch = eval // batch_size\n",
    "            goal_function.append([])\n",
    "        eval += 1\n",
    "        point = json.loads(line)\n",
    "        if 'goal' in point.keys():\n",
    "            goal_function[batch].append(point['goal'])\n",
    "\n",
    "# Clean unfinished batch\n",
    "if len(goal_function[-1])<batch_size:\n",
    "    goal_function.pop(-1)\n",
    "\n",
    "fig, ax = plt.subplots(1)\n",
    "means = []\n",
    "bests = []\n",
    "for ii in range(len(goal_function)):\n",
    "    means.append(np.mean(np.array(goal_function[ii])))\n",
    "    bests.append(np.min(np.array(goal_function[ii])))\n",
    "    for pt in goal_function[ii]:\n",
    "        ax.plot(ii+1, pt, color='tab:blue', marker=\"D\", markersize=2.5, linewidth=0)\n",
    "\n",
    "ax.xaxis.set_major_locator(MaxNLocator(integer=True))\n",
    "ax.set_ylabel('ORBIT')\n",
    "ax.set_xlabel('Iterations')\n",
    "ax.plot(\n",
    "    range(1, len(goal_function)+1), bests, color=\"tab:red\", marker=\"D\",\n",
    "    markersize=5.5, linewidth=0, fillstyle='full'\n",
    ")\n"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "8fc56ae400e717d872a76f4d6b257151d16696a9d0a72e6998d355f9b43887c7"
  },
  "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}