{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run an EasyVVUQ campaign to analyze the sensitivity for the Ishigami function\n",
    "\n",
    "This is done with SC."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-06-07T14:15:21.463066Z",
     "start_time": "2021-06-07T14:15:16.010020Z"
    },
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Run an EasyVVUQ campaign to analyze the sensitivity for the Ishigami function\n",
    "# This is done with SC.\n",
    "import os\n",
    "import easyvvuq as uq\n",
    "import chaospy as cp\n",
    "import pickle\n",
    "import numpy as np\n",
    "import matplotlib.pylab as plt\n",
    "import time\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-06-07T14:15:21.476480Z",
     "start_time": "2021-06-07T14:15:21.468639Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'1.20.3'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-06-07T14:15:21.500137Z",
     "start_time": "2021-06-07T14:15:21.486267Z"
    },
    "code_folding": [
     0,
     1
    ]
   },
   "outputs": [],
   "source": [
    "# Define the Ishigami function\n",
    "def ishigamiSA(a,b):\n",
    "    '''Exact sensitivity indices of the Ishigami function for given a and b.\n",
    "    From https://openturns.github.io/openturns/master/examples/meta_modeling/chaos_ishigami.html\n",
    "    '''\n",
    "    var = 1.0/2 + a**2/8 + b*np.pi**4/5 + b**2*np.pi**8/18\n",
    "    S1 = (1.0/2 + b*np.pi**4/5+b**2*np.pi**8/50)/var\n",
    "    S2 = (a**2/8)/var\n",
    "    S3 = 0\n",
    "    S13 = b**2*np.pi**8/2*(1.0/9-1.0/25)/var\n",
    "    exact = {\n",
    "            'expectation' : a/2,\n",
    "            'variance' : var,\n",
    "            'S1' : (1.0/2 + b*np.pi**4/5+b**2*np.pi**8.0/50)/var,\n",
    "            'S2' : (a**2/8)/var,\n",
    "            'S3' : 0,\n",
    "            'S12' : 0,\n",
    "            'S23' : 0,\n",
    "            'S13' : S13,\n",
    "            'S123' : 0,\n",
    "            'ST1' : S1 + S13,\n",
    "            'ST2' : S2,\n",
    "            'ST3' : S3 + S13\n",
    "            }\n",
    "    return exact\n",
    "\n",
    "Ishigami_a = 7.0\n",
    "Ishigami_b = 0.1\n",
    "exact = ishigamiSA(Ishigami_a, Ishigami_b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-06-07T14:15:21.510035Z",
     "start_time": "2021-06-07T14:15:21.505751Z"
    },
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# define a model to run the Ishigami code directly from python, expecting a dictionary and returning a dictionary\n",
    "def run_ishigami_model(input):\n",
    "    import Ishigami\n",
    "    qois = [\"Ishigami\"]\n",
    "    del input['out_file']\n",
    "    return {qois[0]: Ishigami.evaluate(**input)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-06-07T14:15:21.520594Z",
     "start_time": "2021-06-07T14:15:21.515444Z"
    },
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Define parameter space\n",
    "def define_params():\n",
    "    return {\n",
    "        \"x1\":       {\"type\": \"float\",   \"min\": -np.pi,     \"max\": np.pi,      \"default\": 0.0},\n",
    "        \"x2\":       {\"type\": \"float\",   \"min\": -np.pi,     \"max\": np.pi,      \"default\": 0.0},\n",
    "        \"x3\":       {\"type\": \"float\",   \"min\": -np.pi,     \"max\": np.pi,      \"default\": 0.0},\n",
    "        \"a\":        {\"type\": \"float\",   \"min\": Ishigami_a, \"max\": Ishigami_a, \"default\": Ishigami_a},\n",
    "        \"b\":        {\"type\": \"float\",   \"min\": Ishigami_b, \"max\": Ishigami_b, \"default\": Ishigami_b},\n",
    "        \"out_file\": {\"type\": \"string\",  \"default\": \"output.csv\"}\n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-06-07T14:15:21.530671Z",
     "start_time": "2021-06-07T14:15:21.525407Z"
    },
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Define parameter space\n",
    "def define_vary():\n",
    "    return {\n",
    "        \"x1\":   cp.Uniform(-np.pi, np.pi),\n",
    "        \"x2\":   cp.Uniform(-np.pi, np.pi),\n",
    "        \"x3\":   cp.Uniform(-np.pi, np.pi)\n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-06-07T14:15:21.547702Z",
     "start_time": "2021-06-07T14:15:21.534015Z"
    },
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Set up and run a campaign\n",
    "def run_campaign(sc_order=2, use_files=False):\n",
    "\n",
    "    times = np.zeros(7)\n",
    "\n",
    "    time_start = time.time()\n",
    "    time_start_whole = time_start\n",
    "\n",
    "    # Set up a fresh campaign called \"Ishigami_sc.\"\n",
    "    my_campaign = uq.Campaign(name='Ishigami_sc.')\n",
    "\n",
    "\n",
    "    # Create an encoder and decoder for SC test app\n",
    "    if use_files:\n",
    "        encoder = uq.encoders.GenericEncoder(template_fname='Ishigami.template',\n",
    "                                             delimiter='$',\n",
    "                                             target_filename='Ishigami_in.json')\n",
    "\n",
    "        decoder = uq.decoders.SimpleCSV(target_filename=\"output.csv\",\n",
    "                                        output_columns=[\"Ishigami\"])\n",
    "\n",
    "        execute = uq.actions.ExecuteLocal('python3 %s/Ishigami.py Ishigami_in.json' % (os.getcwd()))\n",
    "\n",
    "        actions = uq.actions.Actions(uq.actions.CreateRunDirectory('/tmp'), \n",
    "                          uq.actions.Encode(encoder), execute, uq.actions.Decode(decoder))\n",
    "    else:\n",
    "        actions = uq.actions.Actions(uq.actions.ExecutePython(run_ishigami_model))\n",
    "\n",
    "    # Add the app (automatically set as current app)\n",
    "    my_campaign.add_app(name=\"Ishigami\", params=define_params(), actions=actions)\n",
    "\n",
    "    # Create the sampler\n",
    "    time_end = time.time()\n",
    "    times[1] = time_end-time_start\n",
    "    print('Time for phase 1 = %.3f' % (times[1]))\n",
    "\n",
    "    time_start = time.time()\n",
    "    # Associate a sampler with the campaign\n",
    "    my_campaign.set_sampler(uq.sampling.SCSampler(vary=define_vary(), polynomial_order=sc_order))\n",
    "\n",
    "    # Will draw all (of the finite set of samples)\n",
    "    my_campaign.draw_samples()\n",
    "    print('Number of samples = %s' % my_campaign.get_active_sampler().count)\n",
    "\n",
    "    time_end = time.time()\n",
    "    times[2] = time_end-time_start\n",
    "    print('Time for phase 2 = %.3f' % (times[2]))\n",
    "\n",
    "    time_start = time.time()\n",
    "    # Run the cases\n",
    "    my_campaign.execute(sequential=True).collate(progress_bar=True)\n",
    "\n",
    "    time_end = time.time()\n",
    "    times[3] = time_end-time_start\n",
    "    print('Time for phase 3 = %.3f' % (times[3]))\n",
    "\n",
    "    time_start = time.time()\n",
    "    # Get the results\n",
    "    results_df = my_campaign.get_collation_result()\n",
    "\n",
    "    time_end = time.time()\n",
    "    times[4] = time_end-time_start\n",
    "    print('Time for phase 4 = %.3f' % (times[4]))\n",
    "\n",
    "    time_start = time.time()\n",
    "    # Post-processing analysis\n",
    "    results = my_campaign.analyse(qoi_cols=[\"Ishigami\"])\n",
    "    \n",
    "    time_end = time.time()\n",
    "    times[5] = time_end-time_start\n",
    "    print('Time for phase 5 = %.3f' % (times[5]))\n",
    "\n",
    "    time_start = time.time()\n",
    "    # Save the results\n",
    "    pickle.dump(results, open('Ishigami_results.pickle','bw'))\n",
    "    time_end = time.time()\n",
    "    times[6] = time_end-time_start\n",
    "    print('Time for phase 6 = %.3f' % (times[6]))\n",
    "\n",
    "    times[0] = time_end - time_start_whole\n",
    "\n",
    "    return results_df, results, times, sc_order, my_campaign.get_active_sampler().count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-06-07T14:17:40.849526Z",
     "start_time": "2021-06-07T14:15:21.555355Z"
    },
    "code_folding": [
     0
    ],
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 8/8 [00:00<00:00, 940.64it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time for phase 1 = 0.044\n",
      "Number of samples = 8\n",
      "Time for phase 2 = 0.045\n",
      "Time for phase 3 = 0.027\n",
      "Time for phase 4 = 0.016\n",
      "Time for phase 5 = 0.023\n",
      "Time for phase 6 = 0.002\n",
      "Time for phase 1 = 0.019\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "100%|██████████| 27/27 [00:00<00:00, 2785.54it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of samples = 27\n",
      "Time for phase 2 = 0.076\n",
      "Time for phase 3 = 0.019\n",
      "Time for phase 4 = 0.004\n",
      "Time for phase 5 = 0.040\n",
      "Time for phase 6 = 0.003\n",
      "Time for phase 1 = 0.018\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 64/64 [00:00<00:00, 2152.86it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of samples = 64\n",
      "Time for phase 2 = 0.145\n",
      "Time for phase 3 = 0.044\n",
      "Time for phase 4 = 0.009\n",
      "Time for phase 5 = 0.105\n",
      "Time for phase 6 = 0.003\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time for phase 1 = 0.044\n",
      "Number of samples = 125\n",
      "Time for phase 2 = 0.210\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 125/125 [00:00<00:00, 2963.33it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time for phase 3 = 0.060\n",
      "Time for phase 4 = 0.009\n",
      "Time for phase 5 = 0.131\n",
      "Time for phase 6 = 0.002\n",
      "Time for phase 1 = 0.017\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 216/216 [00:00<00:00, 3375.84it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of samples = 216\n",
      "Time for phase 2 = 0.238\n",
      "Time for phase 3 = 0.078\n",
      "Time for phase 4 = 0.008\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time for phase 5 = 0.293\n",
      "Time for phase 6 = 0.003\n",
      "Time for phase 1 = 0.024\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 343/343 [00:00<00:00, 4233.40it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of samples = 343\n",
      "Time for phase 2 = 0.348\n",
      "Time for phase 3 = 0.103\n",
      "Time for phase 4 = 0.011\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time for phase 5 = 0.457\n",
      "Time for phase 6 = 0.005\n",
      "Time for phase 1 = 0.034\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 512/512 [00:00<00:00, 4486.04it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of samples = 512\n",
      "Time for phase 2 = 0.563\n",
      "Time for phase 3 = 0.137\n",
      "Time for phase 4 = 0.014\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time for phase 5 = 0.545\n",
      "Time for phase 6 = 0.008\n",
      "Time for phase 1 = 0.037\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 729/729 [00:00<00:00, 4553.03it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of samples = 729\n",
      "Time for phase 2 = 0.456\n",
      "Time for phase 3 = 0.194\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time for phase 4 = 0.018\n",
      "Time for phase 5 = 0.721\n",
      "Time for phase 6 = 0.005\n",
      "Time for phase 1 = 0.026\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 44%|████▍     | 442/1000 [00:00<00:00, 4418.19it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of samples = 1000\n",
      "Time for phase 2 = 0.612\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1000/1000 [00:00<00:00, 4381.37it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time for phase 3 = 0.266\n",
      "Time for phase 4 = 0.029\n",
      "Time for phase 5 = 1.033\n",
      "Time for phase 6 = 0.006\n",
      "Time for phase 1 = 0.025\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 30%|███       | 405/1331 [00:00<00:00, 4035.10it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of samples = 1331\n",
      "Time for phase 2 = 0.931\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1331/1331 [00:00<00:00, 3897.26it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time for phase 3 = 0.390\n",
      "Time for phase 4 = 0.031\n",
      "Time for phase 5 = 1.407\n",
      "Time for phase 6 = 0.012\n",
      "Time for phase 1 = 0.033\n",
      "Number of samples = 1728\n",
      "Time for phase 2 = 1.003\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1728/1728 [00:00<00:00, 3370.05it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time for phase 3 = 0.753\n",
      "Time for phase 4 = 0.053\n",
      "Time for phase 5 = 2.527\n",
      "Time for phase 6 = 0.011\n",
      "Time for phase 1 = 0.028\n",
      "Number of samples = 2197\n",
      "Time for phase 2 = 1.260\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 2197/2197 [00:00<00:00, 4021.69it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time for phase 3 = 0.773\n",
      "Time for phase 4 = 0.079\n",
      "Time for phase 5 = 2.388\n",
      "Time for phase 6 = 0.007\n",
      "Time for phase 1 = 0.019\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 19%|█▊        | 514/2744 [00:00<00:00, 5130.82it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of samples = 2744\n",
      "Time for phase 2 = 1.476\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 2744/2744 [00:00<00:00, 5031.89it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time for phase 3 = 0.630\n",
      "Time for phase 4 = 0.135\n",
      "Time for phase 5 = 3.027\n",
      "Time for phase 6 = 0.055\n",
      "Time for phase 1 = 0.041\n",
      "Number of samples = 3375\n",
      "Time for phase 2 = 1.584\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 3375/3375 [00:00<00:00, 3976.55it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time for phase 3 = 1.068\n",
      "Time for phase 4 = 0.079\n",
      "Time for phase 5 = 4.012\n",
      "Time for phase 6 = 0.018\n",
      "Time for phase 1 = 0.025\n",
      "Number of samples = 4096\n",
      "Time for phase 2 = 2.015\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 4096/4096 [00:00<00:00, 4329.78it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time for phase 3 = 1.173\n",
      "Time for phase 4 = 0.240\n",
      "Time for phase 5 = 6.876\n",
      "Time for phase 6 = 0.022\n",
      "Time for phase 1 = 0.029\n",
      "Number of samples = 4913\n",
      "Time for phase 2 = 3.027\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 4913/4913 [00:01<00:00, 2775.14it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time for phase 3 = 2.001\n",
      "Time for phase 4 = 0.284\n",
      "Time for phase 5 = 8.112\n",
      "Time for phase 6 = 0.025\n",
      "Time for phase 1 = 0.034\n",
      "Number of samples = 5832\n",
      "Time for phase 2 = 2.935\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 5832/5832 [00:01<00:00, 4304.04it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time for phase 3 = 1.595\n",
      "Time for phase 4 = 0.317\n",
      "Time for phase 5 = 11.989\n",
      "Time for phase 6 = 0.048\n",
      "Time for phase 1 = 0.064\n",
      "Number of samples = 6859\n",
      "Time for phase 2 = 4.160\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 6859/6859 [00:01<00:00, 3865.79it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time for phase 3 = 2.020\n",
      "Time for phase 4 = 0.296\n",
      "Time for phase 5 = 11.716\n",
      "Time for phase 6 = 0.030\n",
      "Time for phase 1 = 0.034\n",
      "Number of samples = 8000\n",
      "Time for phase 2 = 4.110\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 8000/8000 [00:02<00:00, 3496.31it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time for phase 3 = 2.624\n",
      "Time for phase 4 = 0.206\n",
      "Time for phase 5 = 15.051\n",
      "Time for phase 6 = 0.043\n",
      "Time for phase 1 = 0.034\n",
      "Number of samples = 9261\n",
      "Time for phase 2 = 4.660\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 9261/9261 [00:02<00:00, 3531.32it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time for phase 3 = 2.939\n",
      "Time for phase 4 = 0.264\n",
      "Time for phase 5 = 18.864\n",
      "Time for phase 6 = 0.139\n"
     ]
    }
   ],
   "source": [
    "# Calculate the stochastic collocation expansion for a range of orders\n",
    "\n",
    "R = {}\n",
    "for sc_order in range(1, 21):\n",
    "    R[sc_order] = {}\n",
    "    (R[sc_order]['results_df'], \n",
    "     R[sc_order]['results'], \n",
    "     R[sc_order]['times'], \n",
    "     R[sc_order]['order'], \n",
    "     R[sc_order]['number_of_samples']) = run_campaign(sc_order=sc_order, use_files=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-06-07T14:17:41.284212Z",
     "start_time": "2021-06-07T14:17:40.854979Z"
    },
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# save the results\n",
    "\n",
    "pickle.dump(R, open('collected_results.pickle','bw'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-06-07T14:17:41.753581Z",
     "start_time": "2021-06-07T14:17:41.285527Z"
    },
    "code_folding": [
     0
    ]
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Total</th>\n",
       "      <th>Phase 1</th>\n",
       "      <th>Phase 2</th>\n",
       "      <th>Phase 3</th>\n",
       "      <th>Phase 4</th>\n",
       "      <th>Phase 5</th>\n",
       "      <th>Phase 6</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.158249</td>\n",
       "      <td>0.044463</td>\n",
       "      <td>0.045157</td>\n",
       "      <td>0.027455</td>\n",
       "      <td>0.015654</td>\n",
       "      <td>0.023364</td>\n",
       "      <td>0.001825</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.159667</td>\n",
       "      <td>0.018620</td>\n",
       "      <td>0.075858</td>\n",
       "      <td>0.018592</td>\n",
       "      <td>0.003819</td>\n",
       "      <td>0.039893</td>\n",
       "      <td>0.002594</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.323880</td>\n",
       "      <td>0.017593</td>\n",
       "      <td>0.145160</td>\n",
       "      <td>0.044086</td>\n",
       "      <td>0.008541</td>\n",
       "      <td>0.104988</td>\n",
       "      <td>0.003109</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.454608</td>\n",
       "      <td>0.043664</td>\n",
       "      <td>0.209539</td>\n",
       "      <td>0.059554</td>\n",
       "      <td>0.008784</td>\n",
       "      <td>0.131073</td>\n",
       "      <td>0.001588</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.637004</td>\n",
       "      <td>0.016862</td>\n",
       "      <td>0.237678</td>\n",
       "      <td>0.078054</td>\n",
       "      <td>0.008347</td>\n",
       "      <td>0.293213</td>\n",
       "      <td>0.002509</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.949501</td>\n",
       "      <td>0.023694</td>\n",
       "      <td>0.348445</td>\n",
       "      <td>0.103378</td>\n",
       "      <td>0.011489</td>\n",
       "      <td>0.457311</td>\n",
       "      <td>0.004618</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>1.301590</td>\n",
       "      <td>0.034174</td>\n",
       "      <td>0.562767</td>\n",
       "      <td>0.137497</td>\n",
       "      <td>0.013694</td>\n",
       "      <td>0.545244</td>\n",
       "      <td>0.007882</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>1.432124</td>\n",
       "      <td>0.037351</td>\n",
       "      <td>0.456433</td>\n",
       "      <td>0.193846</td>\n",
       "      <td>0.018405</td>\n",
       "      <td>0.720531</td>\n",
       "      <td>0.004730</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>1.974014</td>\n",
       "      <td>0.025784</td>\n",
       "      <td>0.611962</td>\n",
       "      <td>0.265501</td>\n",
       "      <td>0.029102</td>\n",
       "      <td>1.033048</td>\n",
       "      <td>0.005515</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>2.797863</td>\n",
       "      <td>0.024832</td>\n",
       "      <td>0.930756</td>\n",
       "      <td>0.390384</td>\n",
       "      <td>0.031481</td>\n",
       "      <td>1.407211</td>\n",
       "      <td>0.011949</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>4.380687</td>\n",
       "      <td>0.032555</td>\n",
       "      <td>1.002752</td>\n",
       "      <td>0.752780</td>\n",
       "      <td>0.053494</td>\n",
       "      <td>2.526922</td>\n",
       "      <td>0.010974</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>4.536749</td>\n",
       "      <td>0.028437</td>\n",
       "      <td>1.259950</td>\n",
       "      <td>0.773154</td>\n",
       "      <td>0.078727</td>\n",
       "      <td>2.388321</td>\n",
       "      <td>0.007400</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>5.343889</td>\n",
       "      <td>0.018555</td>\n",
       "      <td>1.476173</td>\n",
       "      <td>0.630431</td>\n",
       "      <td>0.135218</td>\n",
       "      <td>3.027458</td>\n",
       "      <td>0.054535</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>6.803688</td>\n",
       "      <td>0.040932</td>\n",
       "      <td>1.584220</td>\n",
       "      <td>1.068493</td>\n",
       "      <td>0.079370</td>\n",
       "      <td>4.012036</td>\n",
       "      <td>0.017829</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>10.351328</td>\n",
       "      <td>0.024677</td>\n",
       "      <td>2.015312</td>\n",
       "      <td>1.172575</td>\n",
       "      <td>0.240019</td>\n",
       "      <td>6.875905</td>\n",
       "      <td>0.021886</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>13.479355</td>\n",
       "      <td>0.028786</td>\n",
       "      <td>3.027048</td>\n",
       "      <td>2.000618</td>\n",
       "      <td>0.284229</td>\n",
       "      <td>8.112038</td>\n",
       "      <td>0.024901</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>16.928440</td>\n",
       "      <td>0.033623</td>\n",
       "      <td>2.935383</td>\n",
       "      <td>1.595398</td>\n",
       "      <td>0.317334</td>\n",
       "      <td>11.989267</td>\n",
       "      <td>0.048347</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>18.287801</td>\n",
       "      <td>0.064279</td>\n",
       "      <td>4.160151</td>\n",
       "      <td>2.020112</td>\n",
       "      <td>0.295770</td>\n",
       "      <td>11.715884</td>\n",
       "      <td>0.030149</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>22.068392</td>\n",
       "      <td>0.033853</td>\n",
       "      <td>4.109619</td>\n",
       "      <td>2.624097</td>\n",
       "      <td>0.205618</td>\n",
       "      <td>15.050779</td>\n",
       "      <td>0.043160</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>26.902022</td>\n",
       "      <td>0.034284</td>\n",
       "      <td>4.659992</td>\n",
       "      <td>2.939169</td>\n",
       "      <td>0.264200</td>\n",
       "      <td>18.864318</td>\n",
       "      <td>0.139383</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        Total   Phase 1   Phase 2   Phase 3   Phase 4    Phase 5   Phase 6\n",
       "1    0.158249  0.044463  0.045157  0.027455  0.015654   0.023364  0.001825\n",
       "2    0.159667  0.018620  0.075858  0.018592  0.003819   0.039893  0.002594\n",
       "3    0.323880  0.017593  0.145160  0.044086  0.008541   0.104988  0.003109\n",
       "4    0.454608  0.043664  0.209539  0.059554  0.008784   0.131073  0.001588\n",
       "5    0.637004  0.016862  0.237678  0.078054  0.008347   0.293213  0.002509\n",
       "6    0.949501  0.023694  0.348445  0.103378  0.011489   0.457311  0.004618\n",
       "7    1.301590  0.034174  0.562767  0.137497  0.013694   0.545244  0.007882\n",
       "8    1.432124  0.037351  0.456433  0.193846  0.018405   0.720531  0.004730\n",
       "9    1.974014  0.025784  0.611962  0.265501  0.029102   1.033048  0.005515\n",
       "10   2.797863  0.024832  0.930756  0.390384  0.031481   1.407211  0.011949\n",
       "11   4.380687  0.032555  1.002752  0.752780  0.053494   2.526922  0.010974\n",
       "12   4.536749  0.028437  1.259950  0.773154  0.078727   2.388321  0.007400\n",
       "13   5.343889  0.018555  1.476173  0.630431  0.135218   3.027458  0.054535\n",
       "14   6.803688  0.040932  1.584220  1.068493  0.079370   4.012036  0.017829\n",
       "15  10.351328  0.024677  2.015312  1.172575  0.240019   6.875905  0.021886\n",
       "16  13.479355  0.028786  3.027048  2.000618  0.284229   8.112038  0.024901\n",
       "17  16.928440  0.033623  2.935383  1.595398  0.317334  11.989267  0.048347\n",
       "18  18.287801  0.064279  4.160151  2.020112  0.295770  11.715884  0.030149\n",
       "19  22.068392  0.033853  4.109619  2.624097  0.205618  15.050779  0.043160\n",
       "20  26.902022  0.034284  4.659992  2.939169  0.264200  18.864318  0.139383"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# produce a table of the time taken for various phases\n",
    "# the phases are:\n",
    "#   1: creation of campaign\n",
    "#   2: creation of samples\n",
    "#   3: running the cases\n",
    "#   4: calculation of statistics including Sobols\n",
    "#   5: returning of analysed results\n",
    "#   6: saving campaign and pickled results\n",
    "\n",
    "Timings = pd.DataFrame(np.array([R[r]['times'] for r in list(R.keys())]), \n",
    "             columns=['Total', 'Phase 1', 'Phase 2', 'Phase 3', 'Phase 4', 'Phase 5', 'Phase 6'], \n",
    "             index=[R[r]['order'] for r in list(R.keys())])\n",
    "Timings.to_csv(open('Timings.csv', 'w'))\n",
    "display(Timings)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-06-07T14:17:45.416273Z",
     "start_time": "2021-06-07T14:17:41.757959Z"
    },
    "code_folding": [
     0
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEGCAYAAACdJRn3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABDa0lEQVR4nO3dd3hUZfbA8e9JIwklEQglCUhRUKoBxIIdFVQUxIq9Ynddd1FxLazurgV35Ye9ACIWREVERbGDIiq9SW8SggTQhAAJaef3x51gSGaSSTIzd5Kcz/PMk5l37p17MoQ5c+/7vucVVcUYY4ypigi3AzDGGFP7WPIwxhhTZZY8jDHGVJklD2OMMVVmycMYY0yVRbkdQCg0b95c27Vr53YYxhhTqyxYsGCnqiZ5e65eJI927doxf/58t8MwxphaRUQ2+3rOLlsZY4ypMksexhhjqsyShzHGmCqrF30exhhTmYKCAtLT08nLy3M7lJCLjY0lNTWV6Ohov/ex5GGMMUB6ejqNGzemXbt2iIjb4YSMqrJr1y7S09Np37693/tZ8qjAtEVbGT1zNRlZuSQnxjFiQGeGpKW4HZYxJgjy8vLqXeIAEBGaNWvGjh07qrSfJQ8fpi3aysipy8gtKAJga1YuI6cuA7AEYkwdVd8SR4nq/N7WYe7D6JmrDySOErkFRYyeudqliIwxJnzYmYcPGVm5Xtu3ZuVy95TFHNmqCUe0bsyRrZvQvFEDr9vaZS9jTF1lycOH5MQ4eu/+gnuippAsO8nQ5jxZeDEzI05izrqdTF249cC2zRs14EhPIjmiVWOOaNWEXzKy+XH6i7zDZJIb7CRjX3PGfHApcKslEGPqgPr+5dCShw9juqyl24JXiZN8AFJlJ09Ev8qVvdtx9Hk38fvefFZt283K33JYuW03q37bzWs/bCK/sBiA8yK+5/HoV4kvtf8j+jJPfhLFkLR/uvZ7GWNqLlh9ops2bWLgwIGccMIJ/Pjjj/Ts2ZNrr72Whx9+mMzMTN588026du3KHXfcwbJlyygsLGTUqFEMHjyYTZs2ceWVV7J3714Ann32WY4//ni+/fZbRo0aRfPmzVm+fDm9e/fmjTfeqHH/jtTGZWhFpCHwPJAPfKuqb1a0fZ8+fbTKta2e7gbZW7wdHRq1hJh4iG7o/IxpCNHxFEfHk1MUw878KJLWvkMTKX/pK724OamPrK9aLMaYoFu5ciVHHnkkAP/8aAW/ZOz2ue2iX7PILyou1x4TGUFa20Sv+3RJbsLD53atMIZNmzZx2GGHsWjRIrp27crRRx9Nz549GTduHNOnT2fChAl06dKFLl26cMUVV5CVlUXfvn1ZtGgRIkJERASxsbGsXbuWYcOGMX/+fL799lsGDx7MihUrSE5Opl+/fowePZoTTjjB5+9fQkQWqGofb7GGzZmHiIwHBgGZqtqtVPtA4P+ASOBVVX0cGAq8p6oficg7QIXJo1qy0308odDpTMjfBwX7IH8P5O2GnN+IyN9LQv5eEgr2gZfEAZAsu1iyJYuebRIDHrIxJjS8JY6K2quiffv2dO/eHYCuXbvSv39/RITu3buzadMm0tPTmT59Ok899RTgDDH+9ddfSU5O5vbbb2fx4sVERkayZs2aA6/Zt29fUlNTATjqqKPYtGlTueRRVWGTPIDXgGeB10saRCQSeA44A0gH5onIdCAVWObZ7OAhUYGSkOr9zCOhDZz3TKW773viCOJzt5Vr305TBj83h9OPbMFdp3eiW0pCIKI1xgRQZWcI/R7/mq1eBtWkJMbxzk3H1ejYDRr8OQAnIiLiwOOIiAgKCwuJjIzk/fffp3PnzgftN2rUKFq2bMmSJUsoLi4mNjbW62tGRkZSWFhYoxghjIbqqups4PcyzX2Bdaq6QVXzgcnAYJxEkurZJji/Q/+HIDru4LboOKfdD/FnPUJhZOxBbQq0jCtmzNFZ/LzxdwY98z03T1rAqt98nx4bY8LPiAGdiYuOPKgtLjqSEQM6+9gjcAYMGMAzzzxDSZfDokWLAMjOzqZ169ZEREQwadIkioqC8726RNgkDx9SgNJf/9M9bVOBC0TkBeAjbzuKyHARmS8i86s6cxKAHhfDuWOdMw3E+XnuWKfdz/2jBj9z0P5y8n1ExDdlyLJbmddrJn87JZU563YycMx33PbWQtZl5lQ9zoosneL03YxKdH4unRLY1zemnhqSlsJjQ7uTkhiH4JxxPDa0e0hGWz344IMUFBTQo0cPunXrxoMPPgjArbfeysSJEzn22GNZs2YNDRs2DGocYdVhLiLtgI9L+jxE5CJggKre4Hl8JdBXVe+oyutWq8M8WPL3wVePwE8vwCHtyTlrLC9tbMmEORvZV1DE4J7J3Nn/cDokNarZUMClU+CjO6Gg1Kl1dFzVEuDSKU6s2enOZbz+D/m/rzG1jLcO4/qk1naY+5AOtCn1OBXIcCmWwIiJh7MehyMHwbRbafzWefz9uNu47u57eGluBq//sJnpSzLoc+ghLEnPZr9n6G+VhgLu3gYzRhycOMB5/OHtsPZzaNwKGid7fraGJq2hUSuI9lxqK5t8src4j8ESiDEm7JPHPOBwEWkPbAUuBS5zN6QAaXcC3PIDfPEgzH2Wpms/Z+SQF7nhhFN5cdZ6xn2/sdwuJeVRvCaPfb/DLx/C8vdh0/c4PSxeFO2HLT9Dzm/O/bLiDnGSye8boLBMaeqCXOdMxJKHMfVe2CQPEXkbOAVoLiLpwMOqOk5Ebgdm4gzVHa+qK1wMM7AaNIJBT8MRg2D6HTDudJJO+CsPDryX8d9v9Prxf1DZlP05sGoGLH8P1n8NxYXQ7DA4+V5Y+JqTIMpKaAN3LQVVyP3D2SYnw/Nzm+fnb5D5i/eYfQ5hNsbUJ2GTPFR1mI/2GcCMEIcTWof1d85CZt4P3/0XVn/GyU2up8medeXKoyxtdBKs/AiWvQdrZkJhLjRJhWNvhe4XQqseIALNOnrv8ygZLSYC8U2dW8su5WPyOUlS4YObodfV0PZY53WMMfVOWHWYB0tYdZhXZvVn8NGdFO/JpEgjiJY/h9sVagSFEkksBRDfHLoOgW4XQptjIMLLwLmadHh763CPagBtjoWtCyE/B5p3ht5XQ89hThIyphazDvO61WFe/3QeCG1+JOLprkQU7DvoqSgpJp9orisaySnHX8iV/TpWXJ+mx8XV758o2c9b8snfCys+gAWvOWdLX46CLoOds5F2J9jZiDH1gCWPcBTftPxIKY848inucBoPfbyaWev+4IkLe/gsCV9jvpJPTENIu8K5bV8BCybC0smw7F2nz6XXVdDzMtjwjQ31NaaGxowZw/Dhw4mPjy/33Guvvcb8+fN59tlnQx5XuE8SrL8SUr02S0IqE645mlHnduE7zwTDWWuqMQkyUFp2hbOfhL+thvNfgoZJ8MVD8FQnp28kewugfw71tYmKpq4I0STcMWPGsG/fvso3DDFLHuGqgvIoIsI1/doz/fZ+NG0YzdXjf+aRj34hryC45QgqFB0HPS+F6z6D23525rNomXhKhvoaU9uV9AkG+MvR3r17Oeecc+jZsyfdunXjn//8JxkZGZx66qmceuqpAEyYMIFOnTpx8sknM2fOnAD8MtVjl63CVUV9Dh5HtGrC9NtP4LEZKxk/ZyM/rN/JM8PSOLxlY5eC9kjq7PSLeJOd7gwTtn4RE84+vQ9+W+b7+fR55edJlUzCXTDR+z6tujsThCvw2WefkZyczCeffAI49aomTJjAN998Q/Pmzdm2bRsPP/wwCxYsICEhgVNPPZW0tLSq/GYBY2ce4azHxfDX5TAqy/nppb8gNjqSfw7uxvhr+rAjZz+DnvmeSXM34fooOh+X3UDh1dNh3VdOEjGmNvI2wbaidj91796dL7/8knvvvZfvvvuOhISDq27/9NNPnHLKKSQlJRETE8Mll1xSo+PVhJ151BGnHdGST+86kRHvLuXBD1cwa80OTumcxAvfbnBnmcz+D3mfZ9LtQtjwLbwx1Bn2e+pIaH+ynYmY8FLJGYLPeVAJbeDaT6p92E6dOrFgwQJmzJjByJEjOfPMM8ttU9MVAAPFzjzqkBaNY5lwzdE8NKgL36zK5IFpK9ialYvyZ22saYu2Vvo6AeGrKvHgZ+GOhXDO/yDrV3h9MLx2jqekijG1RA2XbPAlIyOD+Ph4rrjiCv7+97+zcOFCGjduTE6OU3H7mGOO4dtvv2XXrl0UFBTw7rvv1uh4NWFnHnVMRIRw3QnteXHWejJzDj6FrrA2VjD4GuobFQNHXw9HXQ4LX3dm1b92DrQ/CU65Hw6t2WI6xgSdH32S1bFs2TJGjBhBREQE0dHRvPDCC8ydO5ezzjqL1q1b88033zBq1CiOO+44WrduTa9evYK+bocvNsO8jmp/3ydea2MJsPHxc0IdTsUKcmH+BPj+f7B3B3Q8zUkif2y0eSImZGyGedVmmFd62UpEvhCRxFKPDxGRmTUN1ARXcmKc1/ZWCbFe210VHQfH3Qp/WQJnPArblsC40+GDm2yeiDFhyp8+j+aqmlXyQFX/AFoELSITEN6WyQRQ1YMr84aTmIbQ7074y1KITQAtPvh5mydiTNjwJ3kUi0jbkgcicig+F4sw4cLbMpm3nNyRvfuLGPLcHJZvzXY7RN8aNII8H+u6W0l4E0T14TK+N9X5vf3pMP8H8L2IzPI8PgkYXuUjmZAbkpZSrnN8cFoy102YxyUvzeXZy3txaucwPYlMSPU+FDIyBnatd0rOGxNAsbGx7Nq1i2bNmoXNcNhQUFV27dpFbGzVLmn71WEuIs2BY3H6W+eq6s5qRemS+thhXpHtu/O4dsI8Vm/P4dHB3bjsmLaV7xRq3krCR0YDkc5f4cn3wPF3etqMqbmCggLS09PJy8urfOM6JjY2ltTUVKKjD/7/VFGHuc/kISJHqOoqEenl7XlVXVjTgKtLRIYA5+D0vTynqp9XtL0lj/L27C/k9rcW8u3qHdxySkdGnNmZiIgw+7blbT2SdifCp/fAyunQoiuc+3/Q5mi3IzWmTqpu8nhZVYeLyDdenlZVPa2awYwHBgGZqtqtVPtA4P9wlpt9VVUrmeLpjPwCnlLV6yvazpKHd4VFxTw0fQVv/fQr5/ZM5qmLetAgqnwne1ha9Ql88ndn6dy+N8JpD0JsE7ejMqZOqVbyKLVzrKrmVdZWhWBOAvYAr5ckDxGJBNYAZwDpwDxgGE4ieazMS1ynqpme/f4LvFnZWZAlD99UlRdnbeCJz1bRt11TXr6qN4nxMW6H5Z+83fD1o/DzK9C4NZzzFBwRZnNYjKnFajTPA/jBzza/qOps4PcyzX2Bdaq6QVXzgcnAYFVdpqqDytwyxfEE8Kmbl8/qAhHhllM6MnZYGou3ZDH0hR/4dVf4rR3gVWwTOHs0XP8FxCXC5MvgnSth9za3IzOmzvOZPESklYj0BuJEJE1EenlupwDll7SqmRSg9NCadE+bL3cApwMXisjN3jYQkeEiMl9E5u/Y4eJiSbXEeT2TeeOGY9i1J5/zn5/D4i1ZbofkvzZHw02znT6RNTPhub4wbxwseScki/UYUx9V1OdxNXAN0Acofc1nNzBRVadW+6Ai7YCPS122uggYoKo3eB5fCfRV1Tuqe4zS7LKV/9bv2MM1E35mR85+LuvblpkrtrtTlbe6dq2Hj++CjbNBIg6eaBgd5xRntBInxvilWpetVHWiqp4KXKOqp5a6Da5J4vAhHWhT6nEqkBHgYxg/dExqxAe39iOpcQPGz9nkXlXe6mrWEa6aDnFNbYa6MUHkT5/HUyIyWkSCWTFsHnC4iLQXkRjgUmB6EI9nKtC8UQOKisqfkZZU5Q17IpD7h/fnbIa6MQHhT/LogTMSapyI/OjpS6j2mEgReRuYC3QWkXQRuV5VC4HbgZnASmCKqq6o7jFMzW3L9j6YLmzrYpXlayXD+GahjcOYOqrS5KGqOar6iqoeD9wDPAxsE5GJInJYVQ+oqsNUtbWqRqtqqqqO87TPUNVOqtpRVf9d5d/EBJSvqrzJiWFYldcbb4v1ILBvp7M+dWHNlgs1pr7zpyR7pIicJyIf4Ezi+y/QAfgImBHk+IxLfFXlPaJ1k9pRPM7bSoaDn4NjboafXoBX+8POtW5HaUyt5c8kwQ3AN8A4Vf2hzHNjVfXOIMYXEDbaqnqmLdrK6JmrycjKpXViLIc2i2fu+t+55ZSO3DOgc+0tHrf6U5h2KxTmOfNEjrrc1lA3xouazjBvpKp7ghJZiFjyCIziYuXBD5fz5k+/cvPJHbl3YC1OILszYOpw2PQddLsABj3trCFijDmgouThsyS7iDyDZ90Obx8QteGMwwRWRITw6OBuiMCLs9ajKPcNPKJ2JpAmyXDVh/D90/DNfyB9Hlww3oosGuOnitbzsK/qppySBALw0qwNoHDfWbU0gUREwkl/h/YnwXvXw/gBcNo/oN9fIcKfgYjG1F8+k4eqTgxlIKb2ECmVQGZvQIGRtTWBALTpCzd/58xM/+oR2DALzn8JmrR2OzJjwlalKwmKSBJwL9AFODBOs7ol2U3dUJJABOHl2RuAWp5A4hLhwgnQ8TT49F54sR/0vAx+mXbweiJW2sQYwL9laN8E3sFZfOlm4GrAKg0aRIRHBndFBF6evQFV5f6zj6y9CUQEel0FbY6BSUNh7jN/Ppe9xVnZECyBGIN/M8ybeSbyFajqLFW9DmdJWmMQEf55XleuPu5QXvluI//+ZGXtmAdSkaTO3tutNpYxB/hz5lHg+blNRM7BKVjoo/aDqY9EhFHndQXg1e83osAD59TiMxCA3T4KQFptLGMA/5LHv0QkAfgb8AzQBPhrUKMytU5JAhERxn2/EVV4cFAtTiAJqc6lqrLiEkMeijHhqNLkoaofe+5mA6cGNxxTm4kID5/bBYDxczayLjOH9Tv2kJGVV3vWAynR/yGnj6OgVCFIiXCq9X56L5z5L4iMdi8+Y1zm72irG4F2pbf39H0Yc5CSBLJhxx5mr915oL1kPRCgdiSQkk7xrx75c7TVqQ/Ab0vhx+dg+wq4aCI0tCq9pn7y57LVh8B3wJdAUXDDMXWBiLB+R/mKNiXrgdSK5AFOAik3supSaNUdPvoLvHIKXPo2tOrmRnTGuMqf5BGvqvcGPRJTp2Rk1fL1QCpy1DBo3gneuRzGnQFDnoeu57sdlTEh5c9Q3Y9F5OygR2LqFN/rgXhvr3VSe8Pwb6FlN3j3Gvj6X1BcXNlextQZ/iSPv+AkkFwR2S0iOSKyO9iBmdrN13ogg3rWoZIfjVvBNR9D2pUwezRMvgzy7L+GqR/8WUmwsapGqGqcqjbxPK72MrSBIiINRWSBiAxyOxZT3pC0FB4b2p2UxDgEaJ0QS3JCLJPmbmbxliy3wwucqAZw3jNw9lOw9nN49XTYuc7tqIwJukrX8wAQkUOAwzm4ttXsah1QZDwwCMhU1W6l2gfirFQYCbyqqo9X8jqPAHuBFaWGE3tl63mEh8zdeVzw4g/s3V/EuzcfR8ekRm6HFFgbv4N3r4aiQrhwPBx+utsRGVMjNV0M6gacS1epwGKc0iRzq1sYUUROAvYAr5ckDxGJBNYAZwDpwDxgGE4ieazMS1wH9ACa4ySznZY8ao+NO/dy4Qs/EBsdydRbj6dlk1qyJrq//tgMky+H7cudTvT0eVZY0dRaFSUPf/s8jgY2q+qpQBo1KIzoOWP5vUxzX2Cdqm5Q1XxgMjBYVZep6qAyt0ycyYrHApcBN4pIud9DRIaLyHwRmb9jh9VxDBftmzfktWv7krUvn6vH/0x2bkHlO9UmhxwK18+ElN6wYqpnlrr+WVhx6RS3IzQmIPxJHnmqmgcgIg1UdRXgo3JctaUApWtBpHvavFLVf6jqXcBbwCuqWm6Yi6q+rKp9VLVPUlJSgMM1NdE9NYGXruzD+h17uHHifPIK6tj0oZiGsGd7+XYrrGjqEH+SR7qIJALTgC9E5EOc4oiB5K0AUqWdMar6WmWXrEx4OuHw5vzv4qOYt/l37nx7EUXFtbwSb1m+CihaYUVTR/gz2up8Vc1S1VHAg8A4YEiA40gH2pR6nErgE5QJM+f2TObhQV34/JftPDBtee0v5V5ago/C043sLNjUDVVaqNmznsd0T79EIM0DDheR9iISA1wKTA/wMUwYuqZfe247tSNv//wrT3+51u1wAqf/QxBddkKkwN4/YPWnroRkTCBVKXkEgoi8DcwFOotIuohcr6qFwO3ATGAlMEVVV4Q6NuOOv5/ZmYv7pDL2q7VMmrvJ7XACo8fFcO5YSGgDiPPz7KegdTdnMuH88W5HaEyN+DXPo7azobrhr7ComJvfWMhXq7bz3GW9OLt7HZqJXlr+Xnj3Wlg7E078G5z2oLP8rTFhqKZDdRGRliIyyHNrEdjwjIGoyAievSyNPocewl2TF/PDup2V71QbxTSES9+CXlfDd/+FD26GwkBfBTYm+PyZJHgxMBr4FmdU1InACFV9L+jRBYidedQe2fsKuOilH8jIyuOmkzsw+ectZGTl1r7FpCqjCrOfgm/+BR1OgYsnQazrVX+MOUhNZ5gvAc7wTM4rWRzqS1XtGfBIg8SSR+3yW3YeA8fMJqvMBMK46EgeG9q97iQQgEVvOpMHk46Ay9+FJsluR2TMATW9bBVRkjg8dvm5nzHV0iohlpio8n9iJYtJ1Slpl8NlU+CPTfDqGZC50u2IjPGLP0ngMxGZKSLXiMg1wCeAjTU0QbUjZ7/X9jqxmFRZh/WHa2dAcQGMHwCbvnc7ImMq5c8kwRHASzjFCHsCL6vqPcEOzNRvdX4xqbJa94QbvoRGrWDS+bD8fbcjMqZClSYPEXlCVaeq6t2q+ldV/UBEnghFcKb+8raYVFSEMGJAoMuqhZHEtnDdZ5DSB967Dt67Hp7uBqMSnZ9WVNGEEX8uW53hpe2sQAdiTGllF5OKi46ksFi99oXUKfFN4coPILkXLH/PqvKasBXl6wkRuQW4FeggIktLPdUYmBPswIwZkpZyYGTV/sIiLnvlJ+6espg2h8TTPTXB5eiCKDoW9maWby+pymtrgpgwUNHXuLeAc3FqTJ1b6tZbVa8IQWzGHNAgKpKXruxNs4YNuOH1efyWned2SMGVvdVHu1XlNeHBZ/JQ1WxV3aSqw1R1c6lb2YWcjAmJ5o0aMO6aPuzJK+TG1+eTm1/H1gEpzVdVXpsHYsJEHb+AbOqaI1o1YeywNJZnZPO3dxdTXNfWASnhtSovENkAcrNCHo4xZVnyMLVO/yNbcv9ZRzJj2W+M+XKN2+EEh7eqvMfe4nScTxwEe+to7S9Ta/jsMDcmnN1wYnvWZuYw9ut1dGzRiMFH1aGSJSV6XFy+c/yw02HyFTDhLLjqQ7uMZVzj88xDRHJEZLevWyiDNKYsEeFfQ7rTt31TRry3lIW//uF2SKFx2OlwxfuwexuMH+iUNTHGBRV1mDdW1SbAGOA+IAVnedh7gX+FJDpjKhATFcGLV/SmVZNYhr++gK11sXSJN+36wdUfwv7dTgLZUcfqfZlawZ8+jwGq+ryq5qjqblV9Abgg2IFVREQiROTfIvKMiFztZizGXU0bxjDu6j7sLyjihonz2bu/0O2QQiOlN1wzA4qLnEtY25a4HZGpZ/xJHkUicrmIRHo+tC8Hqj1GUkTGi0imiCwv0z5QRFaLyDoRua+SlxmMcyZUANjA93ru8JaNefbyXqz+bTd/mVyHR2CV1bKLU84kKg5eOxe2/Ox2RKYe8Sd5XAZcDGz33C7ytFXXa8DA0g0iEgk8h1P2pAswTES6iEh3Efm4zK0F0BmYq6p3A7fUIBZTR5zcKYmHBnXhy5XbeWLmKrfDCZ1mHZ0E0rAZvD4ENsxyOyJTT1Q62kpVN+F80w8IVZ0tIu3KNPcF1qnqBgARmQwMVtXHgEFlX0NE0oGStTu9ngWJyHBgOEDbtm0DE7wJa1cf3461mXt4adYGDktqxEV92rgdUmgktoFrP4NJQ+DNi+Di16HzwEp3M6Ym/Kmq20lEviq5zCQiPUTkgQDHkQJsKfU43dPmy1RggIg8A8z2toGqvqyqfVS1T1JSUuAiNWFLRBh1Xlf6HdaM+z9Yxv++WE2/x7+m/X2f0O/xr5m2yEfJj7qgcUu45hPnUtY7l1tJdxN0/ly2egUYidO/gKouBS4NcBzipc3nhWtV3aeq16vqHar6XIBjMbVYdGQEz1/Wm8S4aMZ+tY6tWbkosDUrl5FTl9XtBBLfFK6aDql94f0bYPqdVtLdBI0/ySNeVcv2xAV6SEs6UPoaQyqQEeBjmHoiIT4akfLfR+rkMrZlxTZx5oEkHQELJ1pJdxM0/iSPnSLSEc+ZgIhcCGwLcBzzgMNFpL2IxOCc2UwP8DFMPVKvlrEtKybemQNSVklJd2MCwJ/kcRvOMrRHiMhW4C7g5uoeUETeBuYCnUUkXUSuV9VC4HZgJrASmKKqK6p7DGPq3TK2ZVlJdxNkFY628gyhvUVVTxeRhkCEqubU5ICqOsxH+wxgRk1e25gSIwZ0ZuTUZeQW/DkYLyYyom4vY1taQqrnkpWXdmMCoMIzD1UtAnp77u+taeIwJlTKLmMbFSFERkDvQw9xO7TQ8FXSvc2xoY/F1En+VNVdJCLTgXeBvSWNqjo1aFEZEwCll7HdtHMv5z77Pbe8uYD3bj6e2OhIl6MLspJqvF894lyqSkiBhi1h+buQ0guOu9Xd+Eyt50/yaArsAk4r1aY4cy2MqRXaNW/ImEuO4vqJ83lw2nKevLCH1xFZdUrZku5FBfDedTBzJEREwjE3uRebqfX8mWF+bSgCMSbY+h/ZkjtPO4yxX6/jqLaJXH7MoW6HFFqR0XDheHj3Gvj0HpAI6Huj21GZWqrS5CEiscD1QFcgtqRdVa8LYlzGBMVfTu/E4vRsRk1fQZfWTUhrW0/6QEpERsOFE+C9a2HG3502SyCmGvwZqjsJaAUMAGbhTOCzjnNTK0VGCGMvPYqWTWK55Y2F7NzjfT5InRYV4ySQzmc7CWTeOLcjMrWQP8njMFV9ENirqhOBc4DuwQ3LmOBJjI/hxSt688e+fO54axGFRcVuhxR6UTFw0UToNBA+uRvmT3A7IlPL+JM8Cjw/s0SkG5AAtAtaRMaEQLeUBP59fnfmbthV90uW+BIV41TgPfxM+PguWDDR7YhMLeJP8nhZRA4BHsQpGfIL8ERQozImBC7sncoVx7blpdkbmLEs0BV3aomoBnDxJDjsDPjoL7BwktsRmVrCn9FWr3ruzgI6BDccY0LroUFdWZGxmxHvLqFTy0Yc1qKx2yGFXnQsXPIGTL4Mpt/hjMJKu9ztqEyY82c9j2aetcIXisgCERkjIs1CEZwxwRYTFcHzl/ciLiaS4ZMWkJNXUPlOdVF0LFz6FnQ4BT68DRa/7XZEJsz5c9lqMpAJXABcCOwE3glmUMaEUuuEOJ4Z1ovNu/Yx4t2lqNaTNdDLio6FYW9Dh5Nh2i3w8d22HojxyZ/k0VRVH1XVjZ7bv4DEIMdlTEgd17EZ9w08gs9W/MZLsze4HY57ouPg0reheSeYP87WAzE++ZM8vhGRS0UkwnO7GPgk2IEZE2o3nNiec3q05snPVjFn3U63w3FPTDzk7y3fbuuBmFL8SR43AW8B+Z7bZOBuEckRES8rzhhTO4kIT17Qg45Jjbjj7UVsrQ8LR/my29YDMRXzZ7RVPRx+Yuqrhg2iePHK3gx+dg6XvjSXomJlW3YeyYlxjBjQ+UCV3jrP1gMxlfDnzAMR6SEi54nI0JJbsAOrJJ62IjJdRMaLyH1uxmLqno5Jjbi4Typb/sglIzsPBbZm5TJy6jKmLfLxjbyu8bUeSNqVoY/FhCV/huqOB8bjjLY613MbVN0Dej7wM0VkeZn2gSKyWkTW+ZEQOgGfeIozdqluLMb4MnPF9nJtuQVF9Wc2eo+L4dyxkNAGEGicDHHN4KcXYPsvbkdnwoA/63kcq6qB/IB+DXgWeL2kwbPc7XPAGUA6MM+zAFUk8FiZ/a8DFgH/EJFLcAo3GhNQGT76O3y110ll1wP5fSNMOAsmDYHrPoOmNme4PvPnstVcEQlY8lDV2cDvZZr7AutUdYOqlnTKD1bVZao6qMwtE7gWeFhVT8Mp1FiOiAwXkfkiMn/Hjh2BCt/UE8mJXi7ZVNBeLzRtD1dOcxaVen0wZNeTS3jGK3+Sx0ScBLJaRJaKyDIRWRrgOFKA0r1z6Z42Xz4D7hSRF4FN3jZQ1ZdVtY+q9klKSgpYoKZ+GDGgM3FllqqNjBBGDOjsUkRhosURcOVU2PeHcwaytx4Paa7n/LlsNR64ElgGBKt2tbf1QH1O81XV5Tiz3Y0JipJRVaNnriYjK5eGDaLYs7+Qur5yrV+S0+Cyd+CNoTDpfLjmY4hNcDsqE2L+JI9fVXV6kONIB9qUepwKZAT5mMZUaEhayoEkUlhUzKUv/8j9U5fRPSWBDkmNXI7OZe36OdV4Jw+Dty6BK6Y6kwtNveHPZatVIvKWiAwL4lDdecDhItJeRGKAS3HKvxsTFqIiIxg7LI3oqAhuf2sReQVFbofkvk5nwtBXYMtPMOVKKMx3OyITQv4kjzhgP3AmgRmq+zYwF+gsIukicr2qFgK3AzOBlcAUVV1R3WMYEwzJiXH896Ke/LJtN/+ZsdLtcMJDt6Fw7v/Bui9h6g1QVOh2RCZE/Jlhfm0gD6iqw3y0zwBmBPJYxgRa/yNbcuOJ7Xnlu40c16EZZ3Vv7XZI7ut1FezPgZn3Q8xf4LxnIMKv+cemFvNnkmCqiHzgmdi3XUTeFxGrUWDqrREDjqBnm0TueX8pW37f53Y44eG42+Dke2HxG04Sqa9l7esRf74eTMDpf0jGGT77kafNmHopJiqCZ4elAXD724vILwzWIMRa5pSRcMwtziz0bx93OxoTZP4kjyRVnaCqhZ7ba4BNnDD1Wpum8Tx5QQ+WbMli9MxVbocTHkRgwH/gqCtg1uPw/g22mFQd5k/y2CkiV4hIpOd2BbAr2IEZE+7O6t6aq447lFe+28hXK8vXwqqXIiLgvLGQ3AuWvWuLSdVh/iSP64CLgd+AbTiT864LZlDG1Bb3n30kXVo34W/vLmFbdj2qe1WRiEjYm1m+3RaTqlMqTR6q+quqnqeqSaraQlWHqOrmUARnTLiLjY7k2cvSKCgs5s63F1FYZP0fgO+6V7aYVJ3hz2iriSKSWOrxIZ4y7cYYoENSI/4ztDvzNv3B01+ucTuc8OBr0ShbTKrO8OeyVQ9VzSp5oKp/AGlBi8iYWmjwUSlc0qcNz3+7ntlrrIqzz8Wkjr019LGYoPAneUSIyCElD0SkKf7VxDKmXhl1XlcOS2rE3VMWk7k7z+1w3FV2MalGrSAqDuaPh7023qYu8Cd5/Bf4QUQeFZFHgB+AJ4MbljG1T1xMJM9d3os9+wu5653FFBXX84lyPS6Gvy6HUVnw99VOKfesX+GtiyF/r9vRmRoS9WMmqGcxqNNwSqd/paq1ah3KPn366Pz5890Ow9QTU+Zt4Z73l3JWt1YsTc8mIyuX5MQ4RgzofKBKb7218iOYchUcdgZc+hZE2kWMcCYiC1S1j7fn/CpAo6q/qOqzqvpMbUscxoTaRX1S6d02kU+X/8bWrFwU2JqVy8ipy5i2qJ6vvnfkuXD2U7B2Jnz8FytjUotZ9TJjAkxEyMgu3+eRW1DE6JmrXYgozBx9PZx0Dyx6A775t9vRmGqyc0ZjguA3L8kDICPLJhICcOr9kLMNZo+GRi2h741uR2SqqMIzD085ki9DFYwxdUVyopdhqhW01zsiMGgMdBoIM0bAL7b2W21TYfJQ1SJgn4jYAsXGVMGIAZ2Ji448qK1BVAQjBnR2KaIwFBkFF06A1D5OEcVNc9yOyFSBP30eecAyERknImNLbsEOrISIdPAc+71SbQ09M99fEZHLQxWLMf4akpbCY0O7k5IYhwARAkmNYjjbFo86WEw8XDYFEtvC28Ngu43HqS38SR6fAA8Cs4EFpW6VEpHxnkWklpdpHygiq0VknYjcV9FrqOoGVb2+TPNQ4D1VvRE4z59YjAm1IWkpzLnvNDY+fg4vXNGb9Kw8/veFlS8pJ76pMwckOg7euMDqX9US/hRGnAi8zZ9J4y1Pmz9eAwaWbhCRSOA54CygCzBMRLqISHcR+bjMrYWP100FtnjuF/kZizGuGdC1FcP6tuGl2euZu95mWJeT2BaueB/y98CkobDvd7cjMpXwpzDiKcBanA/854E1InKSPy+uqrOBsn8FfYF1njOKfGAyMFhVl6nqoDI3L3WdAUjHSSA+fwcRGS4i80Vk/o4dVmvIuO/BQV1o16whd09ZTPa+ArfDCT+tujkTB//Y6FzCKrCRaeHM3/IkZ6rqyap6EjAAeLoGx0zhz7MGcBKBz2m3ItJMRF4E0kRkpKd5KnCBiLyAsyxuOar6sqr2UdU+SUm28KFxX3xMFGMuOYodOfu5f9oy/KnuUO+0PxGGvgxbfoJxA+HprrYSYZjyZ55HtKoemNmkqmtEJLoGxxQvbT7/F6nqLuDmMm17gWtrEIMxrujZJpG/ntGJ0TNXc1rnFlzQ20qUl9P1fFj9GSyd/GdbyUqE4NTMMq7z58xjgWe00yme2yv42WHuQzrQptTjVCCjBq9nTK1y88kd6duuKQ9PX8GW3/e5HU542uxl2K6tRBhW/EkeNwMrgDuBvwC/UOZMoIrmAYeLSHsRiQEuBWyGkKk3IiOE/13SEwHuemexrT7oja8RVzYSK2xUNsM8Aligqv9T1aGqer6qPq2q+/15cRF5G5gLdBaRdBG5XlULgduBmcBKYIqqrqjh72FMrZJ6SDyPDunGgs1/8Nw3690OJ/zYSoRhr7IZ5sXAEhFpW50XV9VhqtpaVaNVNVVVx3naZ6hqJ1XtqKpWGc3US0PSUhh8VDJjv17Lwl//cDuc8OJrJcK+w0Mfi/HKn8tWrYEVIvKViEwvuQU7MGPqg0cGd6NVk1j++s5i9uwvdDuc8FF2JcLGrSG6EcwfB3t3uh2dwY/FoETkZG/tqjorKBEFgS0GZcLZTxt2cekrP3JR71SevLCn2+GEry3zYOIgaN0TrpoO0bFuR1TnVXsxKE+fx3OqOqvsLSiRGlMPHdOhGbee0pEp89P5dNk2t8MJX22OhvNfcuaATLsZim2ggZuC2udhjPHPXad3okdqAvdNXeZzLRADdB0CZzwCKz6Arx91O5p6zfo8jAkD0ZERjLnkKPILi/nbu4spLrbZ5z4dfyf0vga+/x8s8LfMngk0f2aY/zPoURhj6JDUiAcHdeH+D5Yx7vuN3HhSB7dDCk8icPZ/IWsLfPxXSGwDHU9zO6p6x5+qurOATThlSmbhTPJbGOS4jKmXhvVtwxldWjJ65mp+ydjtdjjhKzIKLnoNko6AKVfbOiAu8Ge01Y3AcKCpqnYUkcOBF1W1fygCDAQbbWVqk9/35jNgzGwicGajb8vOIzkxjhEDOjMkzWcN0fopOx1e6Q+R0XDDl9C4ldsR1SnVHm3lcRvQD9gNoKprAV/rbBhjaqhpwxiG9kphe85+MrLzUGBrVi4jpy5j2qKtbocXXhJS4bJ3nPU/3roE8ve6HVG94U/y2O9ZdwMAEYmigiq4xpia+3hJ+SG7uQVFjJ652svW9VzyUXDhePhtKbx/IxTb+nCh4E/ymCUi9wNxInIG8C4+1tAwxgRGRpb3hZB8tdd7nQfCwMdh9Sfw+QNuR1Mv+JM87gN2AMuAm4AZgP3rGBNEyYle6jpV0G6AY26CY26BH5+Hn152O5o6r9Khup6Jgq+IyESgK7BVbQk0Y4JqxIDOjJy6jNyCPy/BxERGMGJAZxejqgUG/BuyfoXP7oWszfDLh06nekKqU2zRFpIKGJ9nHiLyooh09dxPABYDrwOLRGRYaMIzpn4akpbCY0O7k5IYhwBREUJ8TCSnd2npdmjhLSISLngFmqTC3GedFQjRP1citKVsA6aiy1Ynllpn41pgjap2B3oD9wQ9MmPquSFpKcy57zQ2Pn4Ok4cfS3ZeAU9+tsrtsMJfTEMo9lKh2FYiDKiKkkd+qftnANMAVPW3YAZkjCmvT7umXHN8O16fu5kfN+xyO5zwl+OjwKStRBgwFSWPLBEZJCJpOPM8PoMDQ3VD2msnIh0866i/V6ptiIi8IiIfisiZoYzHGDeMGNCZtk3juff9peTm23DUCtlKhEFXUfK4CWe52AnAXaXOOPoDn/h7ABEZLyKZIrK8TPtAEVktIutE5L6KXkNVN6jq9WXapqnqjcA1wCX+xmNMbRUfE8UTF/Rg8659PPW5zfeokLeVCCOinHYTED5HW6nqGmCgl/aZOOuP++s14FmcznYARCQSeA7nclg6MM9TqTcSeKzM/tepamYFr/+A57WMqfOO69iMK45ty/g5Gzm7eyt6H9rU7ZDCU8moqq8ecS5VRcdDwV5AXA2rLvFZ20pExla0o6re6fdBRNoBH6tqN8/j44BRqjrA83ik5zXLJo6yr/Oeql7ouS/A48AXqvqll22H49Tkom3btr03b97sb7jGhLU9+wsZ8PRsGkRHMOPOE4mNjnQ7pPBXmA+Tzof0eXDtDEj1Wq7JlFHd2lY3AycAGcB8YEGZW02kAFtKPU73tHklIs1E5EUgrSTRAHcApwMXisjNZfdR1ZdVtY+q9klKSqphuMaEj0YNonhsaHc27NjLmC/Xuh1O7RAVAxe/Dk1aw+TLrOM8ACqaJNgauAinP6EQeAd4X1X/CMBxvZ07+px4qKq7cJJZ6baxQIVnR8bUVSd1SuKSPm14efZ6zurWip5tEt0OKfw1bAbDJsOrZ8Dbw+C6z5xhvaZafJ55qOouVX1RVU/F6ZROxFlR8MoAHDcdaFPqcSrOGY4xxk/3n3MkSY0bcM97S9lfaKOv/NLiSKeI4vbl8IGtg14Tlda2EpFewF3AFcCn1PySFTgLSh0uIu1FJAa4FLClbY2pgoS4aP5zfndWb8/huW/Wux1O7dHpTDjjUVg5Hb6tsJvVVKCi8iT/FJEFwN3ALKCPql6vqlVasktE3gbmAp1FJF1ErlfVQpxhwDOBlcCUUrPZjTF+6n9kS85PS+H5b9axIiPb7XBqj+Nug7QrYPaTsOy9yrc35VQ02qoY2ACU1IAu2VAAVdUewQ8vMGwlQVOXZe3L5/T/zaZlkwZMu60f0ZH+FMs2FO6H1wdDxiK4Zgak9nY7orBT3dFW7XEmBA7y3M713EruG2PCQGJ8DP8a0o0VGbt5aZZdvvJbVAO45A1o1MIZgbXbul2roqIO883ebjid3SeELkRjTGUGdmvFOT1aM/ardazZnuN2OLVHw+Yw7B3I3+OMwMrf53ZEtUZFfR5NRGSkiDwrImeK4w6cS1lWFN+YMPPIeV1pFBvFiPeWUlhko4j81rILXDAOti2BD28FW67ILxVdtpoEdMZZQfAG4HPgQmCwqg4OQWzGmCpo1qgBo87rypItWYyfs9HtcGqXzgPhjH/Cig9g1hNuR1MrVDRJsINn/Q5E5FVgJ9BWVe2c2JgwdW6P1ny8JIP/fr6G049sSYekRm6HVHscfydkrnKG7yZ1hq7nux1RWKsoeRSU3FHVIhHZaInDmPAmIvxrSDfOeHo21782j/1FxWzLyiM5MY4RAzozJM1nFSAjAueOgd83wAe3wCHtIDnN7ajCVkXJo6eI7PbcFyDO87hkqG6ToEdnjKmyFk1iOad7K976+c/ycVuzchk5dRmAJZCKlIzAeuU0mDjYKV+Ss83WQPeiotFWkaraxHNrrKpRpe5b4jAmjM1as6NcW25BEaNn2joglWqUBH2uhf3ZkJOBrYHunc0mMqYOysjK89Ge67XdlDF/fPk2WwP9IJY8jKmDkhO9rxTdOiE2xJHUUr5Ktlsp9wMseRhTB40Y0Jk4L4tERYiwcedeFyKqZXytdd6weWjjCGOWPIypg4akpfDY0O6kJMYhQEpiHNf1a0fO/kLOGfsd7y+wb9AV8rYGOgJ7d8KCia6EFG58FkasS6wwojGOjKxc7npnMT9v/J3z01J4ZHBXGsdGux1WeFo65c810BNS4cS/O2Xc138Fx97qlHWPrGjAau1XUWFESx7G1DNFxcpz36xjzJdraNM0nrGXptlKhP4qKoTP/wE/vQiHne4sLBWb4HZUQVPdqrrGmDooMkK4s//hTLnpOAqLlAte+IGXZq2nuLjuf5GsscgoOOsJGDQGNnzrLGn7+wa3o3KFJQ9j6qk+7Zoy484TOaNLSx77dBVXT/iZzBzvQ3xNGX2uhSunwd5MZ0Lhxu/cjijkwj55iEgHERknIu+VaW8oIgtEZJBbsRlT2yXER/P85b14bGh35m36nbP/7zu+XZ3pdli1Q/sT4YavoGELmDQE5k9wO6KQCmqfh4iMx1k8KlNVu5VqHwj8HxAJvKqqj/vxWu+p6oWlHj8C7AVWqOrHFe1rfR7GVG7t9hzueHsRq37L4YYT2nNEq8Y8/eVaMrJyrTZWRfKy4b3rYd0XcMzNcOa/w6IjfdqirYyeubpG/34V9XkE+zd8DXgWeL1UMJHAc8AZOAtLzROR6TiJpOxq9NeparmvQSJyOvALYDOejAmQw1s2Ztpt/fjPjJW8+v1GRP5c2sJqY1UgNgEuewc+fxB+fA52rnU60uMSa/SyNfnwn7ZoKyOnLiO3oAgIzr9f0EdbiUg74OOSMw8ROQ4YpaoDPI9HAqhq2cRR9nUOnHmIyL+BhkAXnDXWz1fV4jLbDweGA7Rt27b35s2bA/lrGVOn9Xr0C37fm1+uPSUxjjn3neZCRLXEgonwyd1wSHuWtRxC0i+v0UJ3kClJbOk1gqPPu8mvl3E+/JeSW/Dnx1pMVARXH3co3VIS2LO/kD15hezZX0iO5+eBx/sLWbE1m0IvAyCq+u/n5pmHNynAllKP04FjfG0sIs2AfwNpIjJSVR9T1X94nrsG2Fk2cQCo6svAy+Bctgpc+MbUfX94SRzgfINdvCWLnqkJiEiIo6oFel8NzTpSMOlCuu0cjQgg0IodJCx4gJ9V6dD/Wnbu2c/OnHx27d3Pjpz97NyTz849+9m1x7m/IiObsp/9+YXFvPLdwYt8RQg0ahDl3GKdnwlx0V4TBwS2tpkbycPbX5zPD3dV3QXc7OO51wIUkzGmlOTEOLb6+KAZ8twckhNiGdCtFWd1a03vQw8hMsISSYm8lOPYXRRLCzn4/YuTfJIXjKbP3Dbl9omOFJo1bEDzxjE0b9SgXOIoIcCXfzuZxp5kERcd6TWJ93v8a6//fr5qnlWHG8kjHSj97qUCGS7EYYzxYcSAzgddMweIi47kwXOPJDYqkk+X/8abP/3KhDmbaN6oAQO6tuSsbq05pkNToiOdQZyB6LB1kz/x5+YX8cu23azIyGZZejbLtmazNnMPa6P/8Po1OUV2MuPwj9iXfCy0PY7EpBSSGjWgSVzUQUmgog//jn6sDjliQGe+/+B57mIyybKTDG3OGC7lhAG3Vv2N8MGNPo8oYA3QH9gKzAMuU9UVwYrBRlsZU3WVfXju2V/IN6sy+WzFb3yzKpN9+UUkxkdzxpEtSYyPZtKPm8krdc0+LjqSx4Z2rxUJpGyHM0BsdATDT+rAIfExLN+6m+Vbs1m3Yw9FntOEZg1j6JaSQLeUJlw+9xyS2VnudfcTTYOoKCj0JIbmnaFdPzi0H7Q7ARq3OnB8rx/+59/q3/u3dAqFH95BVNGf83YKI2OJGvxMlRa0cq08iYi8DZwCNAe2Aw+r6jgRORsYgzPCaryq/jtoQWDJw5hgyysoYvaaHXy6/De+XLmdnLxCr9u1bNKAWSNOJdZLxd+y3DpzUVWOfewrtu/e73Ob5o0a0D2lCd1TEjwJI4HWCbEHzh7mTX+JbgseIE7+7DvK1RiW9/4XR599LWxbDJu+h81z4NcfIX+Ps1HTjnDo8RAZTdHCN4ks/jOGwshYos563Ek0ezNh7w7Ys8O5vyfTKdpYcj/rV7z2BiS0gb8u9/u9sNpWljyMCZn8wmI6PfBphdskxkfTsnEsLZo0oGWTWFp6frZoHEurhFgWb/mDxz9dVaMzl8qSj6qyLTuPNdtzWLt9D2szc1izfQ/rMvewZ7/35Afw0/39admk8lkC86a/RJuFo2mhO8mU5r5HWxUVwm9LYPMPsGkO/PqDM3fEb+KUim/YwvnZqAUse9f3tqOy/H9lSx6WPIwJJV/X7BPjo7nhhPZs372f7bvz2J6zn8zdeWTm7D9w+acicdGRXHJ0G5rERdMk1hlZ1CQu2vkZG01CvNP+5S/buf+D5QdddoqJjOCsbq2IiYpgbWb5JNG8UQM6tWzE4S0aMW1xBtm5BeWOH5KhysVF8EgzfI4jGvoKNExybo1aQHwziChzJvd0N2fp3LICeObh/jRIY0yd46vDfdS5Xb2eORQVK7/vzWf77jwyc/K47jXvX/ZyC4qYujCdnP2FVPV7b35RMR8uyTiQJC7sncphLRrRqWVjDm/RiEMaxhzYNq3tIV7jHzGgc9UOWh0RkU4JeF8f/v70WfR/yFlzvaBUAo+Oc9oDxJKHMSbgShKEv30WkRFCUuMGJDVuACSQ4mOocMk3/+JiJWd/IbtzC8jOLWB3XgG7cwvYnVtIdm4B/56x0utxBJj/wOkBjz/gavrhX5JgSq9H0v+hKnWWV8YuWxljwo630U5V6fPwddmsVs2QL7sYVYA//P1hl62MMbVKTb/5+7psFpLLToHS4+KQJ4uqsORhjAlLQ9JSqn2ZyPXLTvWAJQ9jTJ1Uk+RjKhf2i0EZY4wJP5Y8jDHGVJklD2OMMVVmycMYY0yVWfIwxhhTZfVikqCI7ADCeR3a5uClfnP4sPhqxuKrGYuvZmoS36GqmuTtiXqRPMKdiMz3NYszHFh8NWPx1YzFVzPBis8uWxljjKkySx7GGGOqzJJHeHjZ7QAqYfHVjMVXMxZfzQQlPuvzMMYYU2V25mGMMabKLHkYY4ypMkseISAibUTkGxFZKSIrROQvXrY5RUSyRWSx5xa49SL9i3GTiCzzHLvcylniGCsi60RkqYj0CnF8nUu9N4tFZLeI3FVmm5C+hyIyXkQyRWR5qbamIvKFiKz1/DzEx74DRWS15/28L4TxjRaRVZ5/ww9EJNHHvhX+PQQxvlEisrXUv+HZPvZ16/17p1Rsm0RksY99Q/H+ef1cCdnfoKraLcg3oDXQy3O/MbAG6FJmm1OAj12McRPQvILnzwY+xVnJ81jgJxdjjQR+w5nA5Np7CJwE9AKWl2p7ErjPc/8+4Akf8a8HOgAxwJKyfw9BjO9MIMpz/wlv8fnz9xDE+EYBf/fj39+V96/M8/8FHnLx/fP6uRKqv0E78wgBVd2mqgs993OAlUBtW2hgMPC6On4EEkWktUux9AfWq6qrVQNUdTbwe5nmwcBEz/2JwBAvu/YF1qnqBlXNByZ79gt6fKr6uaoWeh7+CKQG+rj+8vH++cO196+EiAhwMfB2oI/rrwo+V0LyN2jJI8REpB2QBvzk5enjRGSJiHwqIl1DGxkKfC4iC0RkuJfnU4AtpR6n414CvBTf/2ndfA8BWqrqNnD+cwMtvGwTLu/ldThnk95U9vcQTLd7LquN93HJJRzevxOB7aq61sfzIX3/ynyuhORv0JJHCIlII+B94C5V3V3m6YU4l2F6As8A00IcXj9V7QWcBdwmIieVeV687BPycd4iEgOcB7zr5Wm330N/uf5eisg/gELgTR+bVPb3ECwvAB2Bo4BtOJeGynL9/QOGUfFZR8jev0o+V3zu5qWtSu+hJY8QEZFonH/gN1V1atnnVXW3qu7x3J8BRItI81DFp6oZnp+ZwAc4p7WlpQNtSj1OBTJCE91BzgIWqur2sk+4/R56bC+5nOf5mellG1ffSxG5GhgEXK6eC+Bl+fH3EBSqul1Vi1S1GHjFx3Hdfv+igKHAO762CdX75+NzJSR/g5Y8QsBzfXQcsFJV/+djm1ae7RCRvjj/NrtCFF9DEWlcch+nU3V5mc2mA1eJ41ggu+TUOMR8fuNz8z0sZTpwtef+1cCHXraZBxwuIu09Z1KXevYLOhEZCNwLnKeq+3xs48/fQ7DiK92Pdr6P47r2/nmcDqxS1XRvT4bq/avgcyU0f4PBHA1gtwMjG07AOSVcCiz23M4GbgZu9mxzO7ACZ9TDj8DxIYyvg+e4Szwx/MPTXjo+AZ7DGaGxDOjjwvsYj5MMEkq1ufYe4iSxbUABzje564FmwFfAWs/Ppp5tk4EZpfY9G2d0zPqS9ztE8a3DudZd8nf4Ytn4fP09hCi+SZ6/r6U4H2atw+n987S/VvI3V2pbN94/X58rIfkbtPIkxhhjqswuWxljjKkySx7GGGOqzJKHMcaYKrPkYYwxpsoseRhjjKkySx7GVIGI/MNTwXSpp2LqMZ72aBF53FPJdLmI/CwiZwXomKeIyMeBeC1jAiXK7QCMqS1E5Dicmdm9VHW/Z/Z6jOfpR3GqnHbzPNcSOLmax4lU1aIaxBmlfxY/NCYoLHkY47/WwE5V3Q+gqjsBRCQeuBFoX+q57cCUsi8gIv2Bp3D+780DbvEkm03AeJzZyM+KSBYwBtiJU7OrZP+GOHW7unteY5Sqfigi1wDnALFAQ+C0wP7qxhzMLlsZ47/PgTYiskZEnheRkjOLw4BftZKidCISizM7+RJVLfnwv6XUJnmqegJOQcdXgHNxqre2KrXNP4CvVfVo4FRgtCehABwHXK2qljhM0FnyMMZP6hRd7A0MB3YA73i+8furM7BRVdd4Hk/EWXCoREmhvSM8261VpwTEG6W2ORO4z7OC3bc4ZxptPc99oarVWR/DmCqzy1bGVIGnL+Jb4FsRWYZTeG4K0FZEGquzKI8v3spgl7a39KEqeI0LVHX1QY1Ox/1e77sYE3h25mGMn8RZR/3wUk1HAZvVqU47DhjrqVCKiLQWkSvKvMQqoJ2IHOZ5fCUwy8uhVgHtRaSj5/GwUs/NBO4oVT04rSa/kzHVZcnDGP81AiaKyC8ishRnvehRnucewLmU9YuILMfpt9hRemdVzQOuBd71nLUUAy+WPYhnu+HAJyLyPVB6ud1HgWhgqec4jwbstzOmCqyqrjHGmCqzMw9jjDFVZsnDGGNMlVnyMMYYU2WWPIwxxlSZJQ9jjDFVZsnDGGNMlVnyMMYYU2X/D8jRLteu6wr8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the convergence of the mean and standard deviation to that of the highest order\n",
    "\n",
    "mean_analytic = exact['expectation']\n",
    "std_analytic = np.sqrt(exact['variance'])\n",
    "\n",
    "O = [R[r]['order'] for r in list(R.keys())]\n",
    "plt.figure()\n",
    "plt.semilogy([o for o in O], \n",
    "             [np.abs(R[o]['results'].describe('Ishigami', 'mean') - mean_analytic) for o in O],\n",
    "             'o-', label='mean')\n",
    "plt.semilogy([o for o in O], \n",
    "             [np.abs(R[o]['results'].describe('Ishigami', 'std') - std_analytic) for o in O],\n",
    "             'o-', label='std')\n",
    "plt.xlabel('SC order')\n",
    "plt.ylabel('RMSerror compared to analytic')\n",
    "plt.legend(loc=0)\n",
    "plt.savefig('Convergence_mean_std.png')\n",
    "plt.savefig('Convergence_mean_std.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-06-07T14:17:45.822277Z",
     "start_time": "2021-06-07T14:17:45.419336Z"
    },
    "code_folding": [
     0
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEGCAYAAACdJRn3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABLaElEQVR4nO3dd3hUZfbA8e9JL/RQQwcVpUkXEAVFKWJBRLBQrKhrQf3ZdRV1V2xrwYaogIiyIqIgRRQVlAVpgkgTkSKhE3p6Ob8/7gRDMpNMkikJOZ/nmWfm3rnlzGQy79y3nFdUFWOMMaYoQoIdgDHGmLLHCg9jjDFFZoWHMcaYIrPCwxhjTJFZ4WGMMabIwoIdQCBUr15dGzVqFOwwjDGmTFm5cuUBVa3h7rlyUXg0atSIFStWBDsMY4wpU0Rku6fnCq22EpFYEQnJtRwiIjG+Cs4YY0zZ402bx3dA7sIiBpjvn3CMMcaUBd4UHlGqejxnwfXYrjyMMaYc86bNI0lE2qnqLwAi0h5I8W9YnolILPA2kA4sUNWPgxWLMcbkyMjIICEhgdTU1GCHUmRRUVHUq1eP8PBwr/fxpvC4F/hMRHa5lusAg4senmciMh64FNinqi1zre8DvA6EAu+r6vPAAGCaqn4lIp8CVngYY4IuISGBihUr0qhRI0Qk2OF4TVVJTEwkISGBxo0be71foYWHqi4XkTOBZoAAG1U1o/ihujUReBOYlLNCREKBt4CLgQRguYjMBOoBv7k2y/JxHCd8uWonL837nV2HU4ivEs2DvZvRv21df53OGFPGpaamlrmCA0BEiIuLY//+/UXaz2PhISIXqur3IjIgz1OniwiqOr04gbqjqj+KSKM8qzsBm1V1iyue/wJX4BQk9YDVFNBmIyIjgBEADRo0KFI8X67ayaPTfyMlwymbdh5O4dHpTnllBYgxxpOyVnDkKE7cBTWYd3fdX+bmdmmRz1R0dYEduZYTXOumA1eJyDvAV552VtVxqtpBVTvUqOF2jItHL837/UTBkSMlI4uX5v1epOMYY8ypyuOVh6o+5Xr4jKpuzf2ciHhfMVZ87opCVdUk4EZ/nnjX4RQuD1nEQ2FTiZcD7NLqvJg5iK8Od/PnaY0xxuf69OnDzz//TLdu3Zg1a5bPjutNg/nnQLs866YB7X0WhXsJQP1cy/WAXR629anhFZbxUMb7xEi6c2I5wPPh7xOhIRw4fhHVK0QGIgxjzCksUO2qDz74IMnJybz77rs+PW5BbQZnishVQGURGZDrdgMQ5dMo3FuO077SWEQigGuAmQE4Lw+Ff3qi4MgRI+ncK//lolcW8tmKHdgMjMaY4sppV915OAXl73bVL1ftLPYxly9fTuvWrUlNTSUpKYkWLVqwdu1aevbsScWKFX0XvEtBVx7NcNo2quC0c+Q4BtzqyyBEZArQA6guIgnAU6r6gYjcBczD6ao7XlXX+fK8nsSk7HG7vq4kclqNCjw4bQ1frt7Jc1e2omFcbCBCMsaUIU9/tY71u456fH7VX4dJz8o+aV1KRhYPTVvDlGV/ud2neXwlnrqshcdjduzYkcsvv5wnnniClJQUhgwZQsuWLT1uX1IFtXnMAGaISDdVXeS3CJxzXeth/Rxgjj/P7VblenBkR77VUrkeU2/rwifL/uL5uRvp/dqP3HvRGdzSrTFhoZbd3hjjnbwFR2HrvfXkk0/SsWNHoqKiGDNmTImOVRhv2jwmiMhqYAIwV8tDfU3PJ+GreyAj90B6ge6PEBIiDOnckIvOqsWTM9by/NyNzFy9ixeuak2repWDFnI+a6bCd8/AkQSnMOz5JLQeFLj9jSnHCrpCADj3+e/ZeTh/oo66VaL59LYuxT7vwYMHOX78OBkZGaSmphIb67+aEW9+Lp8BjAOGAptF5DkROcNvEZUGrQfBZWOgcn1AILYmoHA04cQmtStHMW5YB8YOacf+42lc8dYi/j17Pcnpmb6JYc1UeLUljKri3K+ZWrR9v7rHdfWkzv1X93h/jJLub4wp0IO9mxEdHnrSuujwUB7s3axExx0xYgTPPvss119/PQ8//HCJjlUYb0aYK/At8K2IXABMBv4hIr8Cj6jqEr9GGCytB538S/uzG2DRa9DmeqjydyewPi3r0KVpdZ6fu5H3ftrK3LV7eO7KVhxMSi9+T4qcL++cK5+cL++cuPLKzoK0Y84t/TjMezzPVRPO8pwH4OBWyEqDzDTISnfdZ5y8bstCZznv/t89Y1cfxvhAzneBL3tbTZo0ibCwMK677jqysrLo2rUr33//PU899RQbN27k+PHj1KtXjw8++IDevXuX+DVIYbVQIhIHDMG58tgLfIDT66kN8JmqBmLMR4l06NBBSzwZ1OG/4M2OcGY/GDje7SZLtyTy6Be/sWV/EqEiZOV6b6PDQxk9oJV3H45XW7ptcyEsGup3+rugyLllJBXttUgohEVCaIRzy3mcc797ted9H9oKMdWKdj5jyoENGzZw1llnBTuMYnMXv4isVNUO7rb3ps1jCfAR0F9VE3KtXyEiY4sdaVlTpQGcOxIWvgAdb4GGXfNtck6TOObccx4d/zWfY2knV1+lZGQxeu4GejSrQWxkGOEFNLDrkQT3IyQzU5DMVOfLu2pDiKwIkZVc965bRAVSZ/4fUekH8+2fHF2bmAfXQ0iom6Pn2u6FM4lJ2e3+yVdbQNuh0OVOJwZjTLnkTeHRzFMjuaq+4ON4Srdz74VVk2HuwzBigdsv4ajwUI6nuW/32Hs0jTbPfOvaLoQKkeFUjAqjQqTrFhVGxcgw/k/jqCsH8u2/m+r82Op9RJxcNAKEiCDy971kCUvSh/BPHXvSWJVkjeDZlEGcseQv0jOznVuW65aZfdK6sKSreFbG5dt/Qtgg7myeDSvGw/L3oUV/6HoPxLcpwZtqjCmLCkqM+BWgrsf5nlfVy/0XVikVEQMXPwOf3+wUIu2Hu90svkq0254UVWLCuefC0zmelsnxtEyOpTr3x1MzOJ6WyY6DyRxPy2RmVmduD51F7rc9WSN4PmMQM6f/lu+4+XUmKSTTlV4lkV0ax4uZg5iZ3Rm+Wn/ySwoLITI0hIgw5xYeGsJf6V3JCMnOv39aN27s15uYC5+Ape/Aiomw9nNo3B3OvQea9oQymhjOGFM0Hts8RKS72ydcVHWhXyLyA5+0eeRQhfF9IHEz3PMLROXvnps3Ky8Uoc3jSAKHX+1MUnY4INSRgye+vJdVuIgv7uxKtjo5+FWdcLJVUVz3qlz33lL2HUvLd+jalaKYO/K8E4VEeKi4/WHgqRshQNWYcIZ2acTwLg2JC0uFlRPh53fg2G6o1dK5Emk5ANZ9YV19TblibR4uZalwCCgR6Ps8jLsAFr4Ivf+db5Ni96TIyoDPbiQ2NItrs55hQ0atE09Fh4cyuu+Z1KkcXWiIj11yltvC65G+Z1I1NqLQ/R/s3czN/iHcen4T1u86xpjv/uDdhX9ydYd63NLtFhqdcwf89hksfgO+GOFU66Ufh2zXtC+F9RYzxpQ5hbZ5iMjpwGigOblyWqlqEz/GVbrFt4W2Q2DpWGh/A1Q/Pd8m/dvWLXq3u2+fgoRlhA8cz20ZnYvdja+k3QAL23/zvmO89+NWpi5P4OOlf9G3ZW1GnN+PNndcC5u/hU+H/l1w5LCuvsacUrzpqrsIeAp4FSfH1Y2u/Z4qcMdSxKfVVjmO74M32kODLnC9DwbPrZ8JU4dCx1uh38slP14A7DuayoTF25j883aOpWbSqXE1bju/CRd+2gzB3edKYNThQIdpTECUxmqr1atXc8cdd3D06FFCQ0N5/PHHGTzY/SziRa228maEebSqfodTYGxX1VHAhUV7CaegCjWh+0Pwxzz449uSHSvxT5hxJ8S3c1sNVlrVrBTFw33OZMmjPXmi31kkHEzm5g9XsIc4D3soTB7oDEIsB1lujClQSbJIeCkmJoZJkyaxbt06vv76a+69914OHz7sk2N7U3ikikgI8IeI3CUiVwI1fXL2sq7TbVCtKXz9qNNeURwZKfDZcJAQuHqiM1CvjKkQGcYt5zVh4UMX8NrgNryYOZhkPbltJUUj+E66OAMQJ10O47rDb9OK/74ZU5b5IQWQu5Ts6enpnH66U60eHx9PzZo1izxXuSfejPO4F4gB7gGexbnqcN9HtbwJi4A+o+GTQbBsnDNwrqjmPgx7foNrPy3zg+7CQ0Po37Yu9316LlnZmq+r71fZ3dj6WE9Y8yksedPp8jx/FHS+A9oNcwY5GnMqmPuI83/tScJy9ymAZtwFKz90v0/tVk5nHQ8KS8m+bNky0tPTadq0aVFeiUfe5LZa7np4HD9P/1omnd4LTrsIFrwArQZBhSLMl/7rf+GXD6HbfdCsj/9iDLD4KtHMPNyNmen5p+0d+fkGbj2vPy3/MdSp8lv8Bsx7zHn/OtwI59wGleKDELUxAZS34ChsvZc8pWTfvXs3Q4cO5cMPPyQkxDfTR3jT2+oM4EGgYe7tVdXaPcDputt7NLzTBX74F1z2unf77dsAs+6Dht3ggif8G2OAuevqGxkWQpcm1Zi/fi8zVu/i3NPiuPW89nS/YTay8xdY8gYsHgNL3oJWV0P1M2DFBzZOxJRNBVwhAJ7z11WuDzfOLvZp3aVkP3r0KP369eNf//oXnTt3Lvax8/Km2uozYCzwHpBVyLblU40znPaPn9+GDjdDndYFb592HKYOg4gKMPADCPXmz1B2FNTV90hKBv9d9hcT/reNGyYsp1mtitx6fhMuv3I8ET23OwMOV4w/uauvjRMxpxp3cwaFRzvrSyAnJfvWrVt5+OGHeeWVV7jyyisZNmwYV199dQmDPpk3XXVXqmp7n541wPzSVTevlMPwRjuocSbcMNtzmg5V+PwWWDcdhs2Axuf7N65SKj0zm69+3cV7P21h455j1KoUyQ1dG3PdOQ2o/M7ZcNTNXM6V68N9awMfrDFeKHJXXR9PuDZp0iS+/PJLpk+ffiIl+5133snNN99MixZ/T041ceJE2rRp41X8BXXV9abwGAXsA74ATlTIqWr+tK2lVEAKD4AVE2DWvU6vqRZXut9m+fsw+//gwifg/Af9H1Mpp6r89McBxv24hUWbDxAbEcrakGs8jBMBHt/j/EIzppQpjeM8isIf4zyG47R5LAZWum4B+CYug9oNg1qt4Jt/Qnpy/ud3/uJ06z3tYuj2f4GPrxQSEc4/owaTbzmH2fd0o1eL2uzM9jROBHjrHPj968AFaIxxq9DCQ1Ubu7mV39QkBQkJdRrKjuxwehHllnLIGc8RWxMGjAMf9Xg4lbSIr8yrg9vwbvj1+caJJGsEE0MGOlcdUwbDJ9c4syIaY4LCq28wEWkpIoNEZFjOzd+BlVmNukHz/rDoVacuE5x2ji//AUd3O1VaNhNfgSYnncMjGbeQkF2dbBUSsqvzSMYtjEoeQOatP0Kvf8G2n5yrkAXP559y1xjjd9501X0K6IGTGHEO0BdYBEzya2RlWa9nYcMsZ9rajBSIqgSpR6DP81C/Y7CjK/UKGidy6dtLeeaK6+l011XwzROwYDT8OgX6vghnlHxeZmOMd7y58hgI9AT2qOqNwNlA2cuhEUh//QwCZCQD6hQcEgoxBdTlmxMe7N2M6PCTZ2mMDg/hpnMbcSw1k0HvLuH+ufvY1/ttGDYTQiOdUf5TroVD24ITtDHljDeFR4qqZgOZIlIJp+eVtXkU5LtnIDvPVLSa5aw3herfti6jB7SibpVoBKhbJZrRA1rz5GUtmH9/d+664DRmrdlNz5cXMmF3AzJH/OjM8LhloVOVtfBFWPWx35POGVOeeTM6bYWIVMEZJLgSJ03JMn8GVebltHV4u97k42k+lOiIUB7o3YwB7ery1Mx1PP3Vej5dvoNn+w+l410DnaqsH/6Nc+nn6u5rgwxNObV9+3YGDBhAVlYWGRkZ3H333dx+++0+ObY3va3+oaqHVXUscDEw3FV9ZTypXK9o602RNalRgUk3dWLskHYcTcng6rFL+L95B9jfZyzE1oC840RyJqMyppSYvWU2vab1ovWHrek1rReztxQ/LYknderUYfHixaxevZqlS5fy/PPPs2vXLp8cu0j9RVV1m6qu8cmZS0BEzhKRsSIyTUTuCHY8+fR8Mv9ANh+kHjAnExH6tKzD/P/rzj96NGXmrzu58D8L0KQD7newKz9TSszeMptRi0exO2k3irI7aTejFo8qUQHiLiX7pk2biIx0mqjT0tLIzs721UvwqtrKp0RkPHApsE9VW+Za3wd4HQgF3ldVj5nFVHUDcLtrnpH3/Bxy0eVUjfgw9YDxLCYijIf6nMlV7esxauY6dm6Po15I/gIkMzSSsKREiLWOC8a/Xlj2AhsPbvT4/Jr9a0jPTj9pXWpWKk/+70mmbZrmdp8zq53Jw50e9nhMTynZd+zYQb9+/di8eTMvvfQS8fG+yVodjJFqE4GT8o+LSCjwFk434ObAtSLSXERaicisPLearn0ux+ky/F1gw/dS60FOHqZRh517Kzj8rqmrKuvtkOvyDTLM0FAn3fXbneH3uUGK0BhH3oKjsPXeevLJJ/n2229ZsWIFDz30EAD169dnzZo1bN68mQ8//JC9e/eW6Bw5vLryEJGzgfNciz+p6q/FPaGq/igijfKs7gRsVtUtrvP9F7hCVUfjXKW4O85MYKaIzAY+cRPzCGAEQIMGDYobriljRIQpqZ05HpKZbzKqTdqAr2t9DFOugTZDnIm8oioFO2RzCiroCgGg17Re7E7anW99ndg6TOgzodjndZeSPUd8fDwtWrTgp59+YuDAgcU+R45CrzxEZCTwMc7UszWBySJyd4nPfLK6QO7k9gmudZ5i6iEiY0TkXZyBi/mo6jhV7aCqHWrUKMIETabMi68SzczsbnRLH0OTtI/plj6Gmdnd+Cu8MceHfwPn/R/8+gm809Xp3mtMgI1sN5Ko0KiT1kWFRjGy3cgSHTcnJfv111/Pww8/TEJCAikpTgaGQ4cO8b///Y9mzZqV6Bw5vLnyuBk4R1WTAETkBWAJ8EaBexWNu/zlHtP9quoCYIEPz29OIe4mowoLEZLTs+j75s+8PPBOzjmjL3x5uzOf+jm3Q8+nICImiFGb8qRfk34AvP7L6+xJ2kPt2NqMbDfyxPrimDRpEmFhYVx33XUnUrKvW7eOBx98EBFBVXnggQdo1aqVT16DN4WHcPIkUFm4/7IviQSgfq7leoBv+pOZcsfTZFT1q0Vz/9Rfuea9n7n53MY8cPMCohY8C0vHwub5cOW7UM9t9mljfK5fk34lKizyGjZsGMOGOWkHQ0NDWbp0KQC9e/snbY83hccEYKmIfOFa7g+M93Ecy4HTRaQxsBO4BrjOx+cw5YinQYZzR57H6DkbeX/RVhZs2s8rgx6j9ZmXwJd3wgcXQ7f7ofvDEBbh5qjGmBzeDBJ8BbgROAgcAm5U1VeLe0IRmYJT7dVMRBJE5GZVzQTuAuYBG4CpqrquuOcwxpOYiDCe7d+SSTd14nhqJle+vZhX/4wn47ZFcPa18NPL8P6FsHedk9LEUpwY45Y3Mwl+pKpDC1tXmgVsJkFTphxJyeDpmeuYvmonrepW5pVBZ3P6oZ+cVCZJic78LLnnUg+PhsvGWLdr49aGDRs488wzEU9TUJdiqsrGjRt9PpNgi9wLrjEZZXpOc2MAKkeH88rgNowd0o5dh1Po98Yixu1rRtbtSyAs8uSCAyzFiSlQVFQUiYmJFPaDvLRRVRITE4mKiip841w8tnmIyKPAY0C0iBzNWQ2kA+OKG6gxpU2flnXo0Kgaj03/jefmbGT++mp8mpnqvleIpTgxHtSrV4+EhAT2798f7FCKLCoqinr1ipZ7z2Ph4RqgN1pERqvqoyUNzpjSrHqFSN4d2p7pv+xk1Mx17CKOupI/xUlydG2sQ69xJzw8nMaNGwc7jIDxpsHcCg5TLogIV7Wvx7z7zueV7GvypTgB2JcicGxPEKIzpnQJRm4rY0q1+CrRTM/omm8e9Y8ye1JTE2HsebD1x2CHaUxQBTyrrjFlgad51L+OuYyPo96CSVfABY8740JC7DeYKX88fupFpFpBt0AGaUyguZtHHeBQbFOODPkGWlwJ3z8LUwZD8sEgRGhMcBX0k2klsMJ1vx/YBPzherzS/6EZEzz551GPYkDbeP7Yd5x+41az5pz/wCUvw58/wLvnQ4L9S5jyxZtBgmOBmao6x7XcF7hIVf8vAPH5hA0SNL6y6q9D3PXJKvYfS+OflzVnSL0DyGc3wLHd0Ps56HQrlMFBYsa4U9JBgh1zCg4AVZ0LdPdVcMaUJW0bVGXW3d3oeloc//xyLSN/CiHpxu/htJ4w90GYdiOkHQt2mMb4nTeFxwEReUJEGolIQxF5HEj0d2DGlFZVYyMYP7wjD/Q6g1lrdnH5B+vYdOE4J637+hkwroflxjKnPG+qraoBTwHn48yx8SPwjKqWmVZCq7Yy/rL4zwPcM2UVSWlZPDegJVdW3QbTboKkgxAikJVrWlHLjWXKmIKqrQotPHIdpIKqHvdpZAFihYfxp71HU7l7yiqWbT3ItZ0a8NQFcUS92caZMz2vyvWdOe2NKQNK1OYhIl1FZD2w3rV8toi87eMYjSmzalWK4pNbzuH27k2ZsuwvrvroTzT3FUdulhvLnCK8afN4FeiNq51DVX/FqcIyxriEhYbwSN8zeX9YB3YcTGaXxrndLjm6doAjM8Y/vBoaq6o78qzKcruhMeXcRc1rMfue83g5a3C+3Fiq8FWazWZgTg3eFB47RKQroCISISIP4Mz2Z4xxo361GL7MPPek3Fi7s6uxTWsyOHsWLHzJKUmMKcO8yW11O/A6UBdIAL4B7vRnUMaUde5yY0WSzusxE+jzw79g/0a44k2nB5YxZVCBVx6uWQNfU9XrVbWWqtZU1SGqauM8jCmAu9xYaUSwusPzzniQtZ/DhEvg6O4gRWhMyRRYeKhqFlBDRPJPbGCM8ShvbqxalSKpXSmS9xZt48PQAejgybD/d3jvQti1KtjhGlNk3gwSfBdoB8wEknLWq+or/g3Nd2ychykNjqdlcu9/VzF/wz6uP6cBozplEz71ekg6AFe+42TqNaYUKWluq13ALNe2FXPdjDFFUCEyjHeHduC27k34eOlfDJudwuHrv4Y6reGzG2DBC9aQbsoMr0eYl2V25WFKm2krE3hs+m/UqRLFB0Nac9rPj8OvU5yrjyvehgibKd0EX0lHmNcQkZdEZI6IfJ9z832YxpQfA9vXY8qIc0hKy+TKd1ewsPkzcPEzsO5LmHgJHN0V7BCNKZA31VYfAxuBxsDTwDZguR9jMqZcaN+wGl/eeS51q0Rz48TljNfL0Ws+gQN/wLgLYOELlpXXlFreFB5xqvoBkKGqC1X1JqCzn+MyplyoVzWGz+/oykVn1eKZWet5bH09Mm742snG+8NzcGQHoM79V/dYAWJKDW8KjwzX/W4R6ScibYF6fozJmHIlNjKMsUPa848eTZmybAdDvjpOVmhU/g0zUuC7ZwIfoDFueDPC/F8iUhn4P+ANoBJwn1+jKoCIhADPuuJYoaofBisWY3wlJER4qM+ZnFGrIg99vgYJ8zB40LLymlKi0MJDVWe5Hh4BLijJyURkPHApsE9VW+Za3wcnBUoo8L6qPl/AYa7ASZVyECddijGnjP5t69IgLoY9H8QRLwfyPZ8cXRvrh2VKA296WzURka9E5ICI7BORGSLSpJjnmwj0yXP8UOAtoC/QHLhWRJqLSCsRmZXnVhNoBixR1fuBO4oZhzGlVrsGVXkn7Hq3WXn3pwDJZWYST3MK86bN4xNgKlAbiAc+A6YU52Sq+iPOFUNunYDNqrpFVdOB/wJXqOpvqnppnts+nKuNQ659PaaGF5ERIrJCRFbs37+/OOEaEzSTk845KStvQnZ1Jmb1prbuh/F94HDeWRKMCSxv2jxEVT/KtTxZRO7yYQx1gdz/CQnAOQVsPx14Q0TOw5lP3S1VHQeMA2eQoA/iNCZg3GXlBVge0423j70MH1wM10+D2i09HMEY//LmyuMHEXlERBqJSEMReQiYLSLVRKSaD2IQN+s8ftmrarKq3qyqd6vqWz44vzGljrusvABbYtqQPGQWIDChL2z9KfDBGYN3hcdg4DbgB2ABTjvDTcBKwBc5PxKA+rmW6+Hk0zKm3MqblbdulWgGd6jHH/uOc82Moxy8djZUiofJA5z07sYEmDe9rRr7OYblwOki0hjYCVwDXOfncxpT6vVvW5f+beuetK5Xi9rc+ckvXPnxdj66djoNvrkFpt0Ex/ZCl38EKVJTHnnT2ypURC4XkXtE5P6cW3FOJiJTgCVAMxFJEJGbVTUTuAuYhzO97VRVXVec4xtzqut5Vi2m3NqZoykZXDlhPWsumABnXQbzHoV5j0N2drBDNOWEN/N5zAFSgd+AE59MVX3av6H5jmXVNaeaLfuPM3zCMg4cS+ft687mgi3/geXvQaurnay8YTZ/mym5grLqetPbqp6qtvZxTMaYEmhSowKf39GVGycs55aPVvH8lfdwdc94+O5pOL4PBk+GqErBDtOcwrxpMJ8rIr38HokxpkhqVozi09u60LVpHA9+/htvZlyG9n8Htv/PmR992XuWldf4jTeFx8/AFyKSIiJHReSYiBz1d2DGmMJViAzjg+EdubJtXV7+ZhP/3NaKrGs+deZHn/OgZeU1fuNN4fEfoAsQo6qVVLWiqtr1sDGlRERYCK8MOpvbuzdl8s9/ccfPVdDoKuQbLmVZeY0PedPm8QewVsvDfLXGlFEiwiN9z6R2pUienrUejdzvdvStZeU1vuJN4bEbWCAic4G0nJWq+orfojLGFMsN5zamZqUodk+Lo65l5TV+5E211VbgOyACqJjrZowphS5pVYe3Q65zm5V3Vlq7IEVlTjXejDB/GkBEKjqLetzvURljSuSTlM4cC8nkobCpxEsie7QqKURwVdYcWP0JtLEkDqZkCi08RKQl8BFQzbV8ABhmo8CNKb3cZeWNJpVJMa/T8cs7ID0JOt0axAhNWedNtdU44H5VbaiqDXGmo33Pv2EZY0rCXVbeFKL4of0YaHYJzHkAFr0apOjMqcCbBvNYVf0hZ0FVF4hIrB9jMsaUUE5CxZfm/c6uwynUqhxFCPD+4t2cPfhFeofHwPxRkHYMLvwniNu+WcZ45E3hsUVE/olTdQUwBKcR3RhTiuXNynskJYMbJyzjH//9jZcHPsmVEbHw03+cKqzeoyHEm4oIYxzefFpuAmrgzOA3HagO3OjPoIwxvlc5OpyPbj6HcxpX4/7P1jK5xv3Q+U5YOha+uhuyPc7qbEw+3vS2OgTcE4BYjDF+FhsZxvgbOnLnx7/wxIx1JPe9iRHdK8DCF5wrkCvHWUZe4xVv5vP4VkSq5FquKiLz/BqVMcZvosJDGTu0PZe2rsNzc3/nlcyB6MXPwrovYOpQyEgNdoimDPCm2qq6qh7OWXBdidT0W0TGGL8LDw3h9WvaMqhDPcZ89wf/OnQR2u8V2DQPPrka0mw4lymYNw3m2SLSQFX/AhCRhuTLuGaMKWtCQ4TnB7QmJiKMDxZtJaljR57rP5aQGf+Aj/pDm+udBvUjCVC5HvR8EloPCnbYppTwpvB4HFgkIgtdy+cDI/wXkjEmUEJChKcua06FyDDe/GEzyenNeGXgRMI+Gw4JKzjxOzEnpTtYAWIA7xrMvxaRdkBnQID7VDV/xjVjTJkkIjzQuxmxkWG88PVGUjLiGRdTDUnO82+ek9LdCg+Dd1ceuAqLWX6OxRgTRHf0aEpsZChPzlgHUYnuN7KU7sbFq8LDGFM+DOvSiJiIMHbOiKNeiKV0N57ZkFJjzEkGtq/nNqV7tsI76f2CFJUpbTwWHiJSraBbIIM0xgTWlNTOPJJxCwnZ1clWYZ9WJp0wBmTMgqO7gx2eKQUKqrZaidPVwl3GNAWa+CUiY0zQuUvp3k42MTnyeZjYD26YDZXqBDFCE2werzxUtbGqNnHd571ZwWHMKcxdSvdVegY/njMOju91ChC7AinXvGrzEJHLReRl1+1SfwdljAmu/m3rMnpAK+pWiUaAuNgIQgRe31SNowM/tQLEIKoFDxYXkeeBjsDHrlXXAitU9VE/x+YzHTp00BUrVgQ7DGPKtB837efWSStoUqMCn/YVKk0bDBVqWRXWKUxEVqpqB3fPeXPlcQlwsaqOV9XxQB8gaF0uRKS5iEwVkXdEZGCw4jCmvDn/jBq8N6wDW/YfZ/Bc/fsK5MNL7QqkHPK2q26VXI8rF/dkIjJeRPaJyNo86/uIyO8isllEHinkMH2BN1T1DmBYcWMxxhRdTgHyZ04BctV/4dgeK0DKIW8Kj9HAKhGZKCIf4vTCeq6Y55uIc+VygoiEAm/hFArNgWtdVxetRGRWnltNnBkNrxGRl4C4YsZhjCmm88+owfs5BcjXWAFSThXa5gEgInVw2j0AlqnqnmKfUKQRMEtVW7qWuwCjVLW3a/lRAFUdXchxQoHpqnqFh+dH4Erg2KBBg/bbt28vbsjGGDd+3LSfWyatoGmNCnzaByp9fg1UrA3DZ1kbyCmipG0eAF2AHkB312NfqgvsyLWc4Frnlog0EpFxwCTgJU/bqeo4Ve2gqh1q1Kjhs2CNMQ67AinfvJlJ8G3gduA3YC1wm4i85cMYPA1CdEtVt6nqCFW9XlUX+TAOY0wRFViALB0Hr7aEUVWc+zVTgx2u8SFvrjy6A71VdYKqTsDpfdXDhzEkAPVzLdcDdvnw+MYYP8pdgFyTU4Ac3gFzH3LmAUH/ng/ECpBThjeFx+9Ag1zL9YE1PoxhOXC6iDQWkQjgGmCmD49vjPGznF5Ym10FSHZkJfJVIOTMB2JOCQUlRvxKRGbi9GjaICILRGQBsAEoViOCiEwBlgDNRCRBRG5W1UzgLmCe69hTVXVdcY5vjAme7rkKEPJOJJXD5gM5ZRSUGPFlX59MVa/1sH4OMMfX5zPGBFZOAbL74zjqis0HciorKDHiwpwbsBGo6LptcK0zxph8up9Rw+18IFkqvJZ+ZZCiMr7mTW+rQcAy4GpgELDU0oIYYwryScrJ84EkakVA6ZG+wGn7MGWeN9PQPg50VNV9ACJSA5gPTPNnYMaYssvdfCD9QxbxSsQ7MHUYDP4YwiIKOIIp7bzpbRWSU3C4JHq5nzGmnHI3H8hM7cYvrZ6EP76B6bdAVmaQojO+4E0h8LWIzBORG0TkBmA21rhtjClA3vlAqsSEk60w9vj5ZF38b1g/A2beBdnZwQ7VFJO3ua0GAN1wRoP/qKpf+DswX7L5PIwJvo+WbOOfM9Zx2dnxvF7nG0IWPAcdboJ+r4C4SzRhgq2g3FaFtnmISCwwQ1Wni0gznDEa4aqa4etAjTGnrqFdGpGcnsXouRuJDuvNC12TkMWvQ3gM9PqXFSBljDcN5j8C54lIVZyG8hXAYOB6fwZmjDn13Na9KUnpWYz57g9iugziqY7JyJI3IbIi9ChsKh9TmnhTeIiqJovIzTiTML0oIqv8HZgx5tR030Wnk5yWyfuLthLT/SYeapMMC0Y7VyDn3hPs8IyXvCo8XHNuXA/cXIT9jDEmHxHh8X5nkZKRxdsLtxJ78T3c2SIZvv0nhEdDp1uDHaLxgjeFwEjgUeALVV0nIk2AH/wbljHmVCYiPHtFS1LSs3jp2z+JueQxbsxIgTkPQEQstLku2CGaQhRaeKjqjzjtHjnLWwC7tjTGlEhIiPDiwNakZGTx9Jw/iLniXwzOSIEZdzpXIC0slUlpZtVPxpigCQsN4fVr2pLy0QoembmJmAH/4bLMu+DzWyBhJaz/0snEW7ke9HwSWg8KdsjGxUaKG2OCKiIshLFD2nNO42rc+8UfzG/3JlSsA0vesMmkSjFvEiOe6806Y4wprqjwUN4f3pHW9Spzx7Q/SE13M4zMJpMqVby58njDy3XGGFNsFSLDmHhjJ06vWZGI5L3uN7LJpEoNj20eru65XYEaInJ/rqcqAaHu9zLGmOKrHB3ORzd3Yu/LcdTBJpMqzQq68ogAKuAUMBVz3Y4CNp+HMcYv4ipE8k7okHyTSWUrvJl+WZCiMnl5vPJwzRa4UEQmqup2ABEJASqo6tFABWiMKX8+SurE4ZB0HgqbSrwkkkglKnOM3hnzIe2Yk87EBJU3bR6jRaSSK0HieuB3EXnQz3EZY8qx+CrRzMzuRrf0MTRJ+5iOae9wR8Z9tAjZBv+9HjLTgh1iuedN4dHcdaXRH2cejwbAUH8GZYwp39xNJvV9dntWtH4ati6E6SMgOytI0RnwrvAIF5FwnMJjhisVe+GTgBhjTDHlnUwqLtZp/3hlf0cyej7tDB6c8wB4MR+R8Q9vRpi/C2wDfgV+FJGGOI3mxhjjN/3b1qV/27onlmes3sm9n67mH9HdeLfLPYQsGQOxNeGCR4MYZfnlTW6rMcCYnGUR+Qu4wJ9BGWNMXle0qcuhpHRGfbWeR6Kv4oU2icjC5yEmDs4ZEezwyp0ipydRZ95aa/MwxgTcDec25p4LT2Pqyp28GPEPaHYJzH0IfpsW7NDKneImRnwamODLQIwxxhv3XXwGiUnpvPPjdmr0/ic3pRyGL26H6KpwWs9gh1duFDTCfI2np4Ba/gnHGGMKJiI8c0VLDiWn88y8rVTv/x8uTxsBnw6F4TOhXodgh1guFHTlUQvoDRzKs16AxX6LKO/JnMmnHgcqq+pA17r+QD+gJvCWqn4TqHiMMcEXGiK8OrgNR1NWcN/MbVQdOJbzfroePr4abvoaajQLdoinvILaPGbhjCbfnue2DVjgzcFFZLyI7BORtXnW9xGR30Vks4gUOOu9qm5R1ZvzrPtSVW8FbgAGexOLMebUEhkWytih7WkRX4lbpu9gzQUTISQMPrrSEigGgMfCQ1VvVtVFHp7zdo7IiUCf3CtEJBR4C+gLNAeuFZHmItJKRGbludUs5PhPuI5ljCmHKkSGMeGGjtStGs310/expc8kJ33JR1dC8sFgh3dKE/XzIBsRaQTMUtWWruUuwChV7e1afhRAVUcXcpxpuaqtBHge+FZV53vYfgQwAqBBgwbtt2/f7psXZIwpdRIOJTPwnSVkqzLrMqg54zqoGA/Z6XB0l81EWEwislJV3TYiBWMmwbrAjlzLCa51bolInIiMBdrmFDTA3cBFwEARud3dfqo6TlU7qGqHGjVq+Ch0Y0xpVK9qDJNu7kRaZjaD5oWSfPZwOLwVju7EZiL0j2AUHuJmncfLH1VNVNXbVbVpztWJqo5R1fau9WP9Fqkxpsw4o1ZFxt/Qkb1H0zi2+ov8G9hMhD5VYOEhIqEi4rZaqAQSgPq5lusBu3x8DmNMOdS+YVXeGdKOGln73W9gDek+U2DhoapZQLKIVPbhOZcDp4tIYxGJAK4BZvrw+MaYcqxHs5ociXA/FC05unaAozl1eVNtlQr8JiIfiMiYnJs3BxeRKcASoJmIJIjIzaqaCdwFzAM2AFNVdV1xX4AxxuT1ul6bbyZCVZid1jZIEZ16vElPMtt1KzJVvdbD+jk4c4MYY4zPfXi8EwdzzUS4W6uRSjiXZ30D2xZBo27BDrHM8yar7oeu6qUzXKt+d83pYYwxpVJ8lWhmHu7GzPS/C4nKHOfL6GdoPOU6uGku1GoRxAjLvkKrrUSkB/AHzmC8t4FNInK+f8MyxpjiczcT4REqMOfsNyEiBiYPtMbzEvKmzeM/QC9V7a6q5+Pku3rVv2EZY0zx5Z2JsE7lKOIrR/HWL2n8cfFESD8Ok6+ClLyp+4y3Ch1hLiJrVLV1YetKsw4dOuiKFSuCHYYxJoj2Hk1lwNuLScvMZs5lSs2Z10HdDjD0CwiPCnZ4pVJJR5ivcPW06uG6vQes9G2IxhjjX7UqRfHhTR3JyMrmmvnhHL/kTfhrMUy/FbKzgh1emeNN4XEHsA64BxgJrAfcpgQxxpjS7LSaFXl/eAcSDqUwbGk9Mi76F2yYCV8/4vTlNV7zWHiIyHeuh8+o6iuqOkBVr1TVV1U1LUDxGWOMT3VsVI3XB7dh1Y7D3LmlC9md74Jl4+B/rwU7tDKloCuPOiLSHbhcRNqKSLvct0AFaIwxvta3VR2evLQ536zfy9Opg9GWA2H+KPj1v8EOrcwoaJzHk8AjOLmnXsnznAIX+isoY4zxtxvPbczuI6mM+3EL8b0e4LbG+2DGnRBbw+ZC94LHwkNVpwHTROSfqvpsAGMyxpiAeKTPmew+ksrob7ZQd8BLXJp8M0wdBjfMhvg2wQ6vVCu0wdwKDmPMqSokRHj56tZ0aRLHfTO2sKzrOIiuBh8PhCVvwastYVQV597mAjlJMObzMMaYUiNnLvQm1Stw0/QE/ug9EdKTYN7jziRSNpmUW1Z4GGPKvcrR4Uy8qSMVo8K4/otDZIXHkm+OOptM6iSFTQYVIiJrAxWMMcYES53K0Uy8sRMpGVlI8gH3G1k+rBMKmwwqG/hVRBoEKB5jjAmaZrUrMm5oB3ZrnPsNKtcLbEClmDfVVnWAdSLynYjMzLn5OzBjjAmGLk3jWNzoznyTSaVoBMub3h2kqEofbyaDetrvURhjTCny2t42/JRxCw+FTaWuHECB5zKu5fv1p/O/y4MdXengTVfdhcBGoKLrtsG1zhhjTkm7DqcwM7sb3dLHcH76q2QTQpOQPew6nBLs0EoNbyaDGgQsA64GBgFLRWSgvwMLttlbZtNrWi9af9iaXtN6MXtLsWbiNcaUQfFVok883qG1+CyrO9eFfsfZlZOCGFXp4k2bx+NAR1UdrqrDgE7AP/0bVnDN3jKbUYtHsTtpN4qyO2k3oxaPsgLEmHIi70yEb2b2R1AeiZ0VxKhKF28KjxBV3ZdrOdHL/cqs1395ndSs1JPWpWal8vovrwcpImNMIOWdiVAr12deZC/aJc5i4dLlwQ6vVPCmwfxrEZkHTHEtDwbm+C+k4NuTtKdI640xp57+bevSv23dE8spB05D3vyWvbOfY1G19+h2evUgRhd8hQ0SFGAM8C7QGjgbGKeqDwcgtqCpHVvb7foaMTUCHIkxprSIrt6Q7HY3clXIAp79aBa//FW+5z8vbJCgAl+q6nRVvV9V71PVLwIUW9CMbDeSqFA3cxorHE0/GviAjDGlQuQFDxASGsG9EV9y44Tl/L7nWLBDChpv2i5+FpGOfo+kFOnXpB+juo6iTmwdBKFObB1uaXkLB9MOMvL7kaRl2USKxpRLFWsjnW6hT9ZCTg/dw9APlvJXYnKwowoK0ULm7RWR9cAZwHYgCZz2I1Vt7f/wfKNDhw66YsWKEh9nzpY5PPzTw/Rq2IuXur9EiJzS/QaMMe4c3w+vt+Zoo16c/+f1VIwKY9rtXalVyU1tRRknIitVtYO757xp87gdaIozc+BlwKWu+4AQkSYi8oGITMu1roeI/CQiY0WkR6BiuaTJJTzQ4QG+2f4NLy1/icIKXmPMKahCDeg0gkp/zGBK/yokHk9n2AfLOJycftJmJR0rVuKxZmum+nU+Em/aPF5V1e15b94cXETGi8i+vJl5RaSPiPwuIptF5JFCYtiiqjfnXQ0cB6KAgKa5HN5iOEObD2Xyhsl8uO5Dv53HBikaU4p1vQciYjnr97d5b1gHth5I4oYJy0lKywRKPlasxGPN1kx15h/x43wk3lRbvQVMVNUid24WkfNxvuQnqWpL17pQYBNwMc4X/3LgWiAUGJ3nEDfljDERkWmqOtD1OERVs0WkFvCKql5fUBy+qrbKka3ZPPzjw3y97WtGnzeaS5tc6rNjw98fnNxjTaJCoxjVdRT9mvTz6bmMKa1mb5nN67+8zp6kPdSOrc3IdiMD+vkv9PzfPQs/vQy3/495idW5Y/JKujStxouDGzNk7rUcSMmf1j02LJZLm15KSmYKqZmppGalkpqZSkpmyknrElMS0bzziQBVI6sy7fJp1IiugVMx5KLqpIvf8xvsXQs/vQKZblKpVK4P93k/y0ZB1Vbetnk0A7ZRjDYPEWkEzMpVeHQBRqlqb9fyozgHzFtw5D3OicIj17oI4JO8613PjQBGADRo0KD99u1eXSx5LT0rndvn386qfat4u+fbdInv4rNj95rWi91Ju/OtrxNbh28GfuOz85hTW7C/fEsi2D+g3J0/MjSSG1rcwOlVT+dAygEOHE3gwC8fsD+2OolV67Lj6F6OZxxGpODv1GpR1YgKjSIqLIrosGiiwlyPQ/9+PG3TtAKPERdRmbMi4zgrO4Tmx49w1v4txCcdwilOBFBmx8bwetUq7AkLpXZmFiMPHaZfUgqMOuz1+1DSwqOhu/VFqLpqxMmFx0Cgj6re4loeCpyjqnd52D8O+DfOlcr7qjpaRAYAvYEqwDuquqCgGHx95ZHjWPoxhn89nF3HdzGxz0TOrHamT47b+sPWbn91ADzS6RHOjT+XhpUanvzLw8fK8hcPlP34SyrYX74l5ekHVO3Y2nw78Fu/nVdV2Xp0K8PmDuNI2pECtw2TMKqFhFMj+TDV655D9apNSTgQxsL1aVSs8x0Z5O/G6+0PwF6fdGN3Rv7zx2Urtx4+yobwUDZERPBnRDhZru+ByqFRnFWpCWfVakfq6o+YHqGkhfzdMhGVnc2oZKHfnb658ih0hLmqbheRbsDpqjpBRGoAFbw+u5t43J2mgPMn4jTa5143HZheghh8omJERd7p+Q5D5g7hjvl3MPmSydStULfwHT04kHKAcWvGeSw4QiWU55c9D0B8bDxd4rvQNb4r59Q5h8qRlYt93rzyfvHk1LcCAfviKcmXf2mIP9gKSrFT2t+DbUe2uS04wMnyMHzucNrUbEPbmm1pU6MNVaKqFPtcaVlprE9cz6p9q1i1bxWr963mcNrhAvf5/PLPqR5dnSqRVQhJOwavtYYK6dBvFACvRG/i7eWRxMR/gcrfjejhEsnIdiO9imvkocOMilFS83z5P3joGP1a3wS1WkLt1qRWrsvmo1tZn7ie9Ynr2XBwA5M3fUpGlJD3qzY1JITXq1bGV3/9QgsPEXkK6IBTdTUBCAcmA+cW85wJQP1cy/WAXcU8VtDViq3F2IvGMmzuMG7/9nY+6vtRkT/MR9KOMGHtBD7Z+AnpWel0qt2JX/f/etJ4kpxfja1rtGbJriUs3rWYedvm8fkfnxMiIbSs3pKu8V3pGt+VVtVbMW/bvCJ9+WZmZ3I47TCJKYm8uPzFoH7xFOXLX1U5mn6UxNREElMSSUxN5Lmlz7mNf/TS0USHRVM9ujrVo6sTFx1HZGhkgXGU1asXT6l0dift5pMNn9Cjfg/iK8QHOKqCJWckM27NOD5c/yGCuP0RFRsWS2Z2JpPWTWL82vEANK7cmHY1250oUBpUbICIuP37dYnvwup9q1m9bzWr9q1iXeI6MrIzAGhUqRE96vegbc22vLnqTfan7M93/jqxdTij6hl/r4iqDF3vhu+fhYSVUK899110Oqv+upgluyCyxjwk/DCaUYXUg33JONKm8DdClX77d4CnaqeLRv19eqBl9Za0rN7yxLqMrAzaT27v9v3bk+G7Qc7eVFutBtoCv6hqW9e6NSVo8wjDaTDvCezEaTC/TlXXFfM1FMpf1Va5/bL3F2795lbOijuL93q9R3RYdKH7JGck8/GGj5mwdgLHM47Tt3Ff7mxzJw0qNfDqiysjO4O1B9ayeNdiFu9czNrEtWRrNpEhkWRoBtmafWLb8JBwLm1yKXUr1OVg6kESUxM5mHqQgykHOZh6kMNphz1e8eQ2vPlw2tdqT7ta7Xx6tZObpyqL2PBYejbo6cSe8vdryMzOLPa5KkVUOqkwqR5dnRrRNdhxbAczNs8gPfvvX46nQrVPqISSpVkANKvajB71e3BB/QtoHtfcr9WgBVHVE93f9ybv5fKml9Oqeiv+s+I/HqvdUjJTWHtg7YlCYPX+1RxLd6qJqkVVo3ZMbTYd2kSm/v3ZyF0ghYWE0SKuhXP1UrMNbWq0IS7676lni1Ttl3P1Ed8WhjoVIl2f/45dh0/+AQMQGxnKP3qcRo2KkdSsGOm6j6JabAShIQLpyTDjH7DOfSKP5Og6xDy8sdD31FftpiVt81imqp1E5BdVbSciscASbwoPEZkC9ACqA3uBp1T1AxG5BHgNp4fVeFX9t9evphgCUXgAzN8+n/sX3E/3+t15tcerhIW4v7BLz0rns02fMW7NOA6mHqRH/R7c1eYumlVrVqLzH0k7wtLdS3nif0+Q4q6nhUvF8IrERcdRLaoa1aKq5Xv875//TWJqYr79IkKcaTnTs9MRhDOqnkGH2h1oX6s97Wu1p1pUtRPbevurXVXZm7yXzYc38+fhP/nz8J98sdlzBpzasbWJi4o7EfNJj6PjiIuK4475d7A3eW++fWvG1GTMhWNITEl0Gjw93Ap678pKp4XZW2bz5OInSc/KX/g1j2vOwh0L+WHHD6zev5pszaZmdE161O9Bj/o96FSnE5GhkQG58tpyeAvPLXuOpbuX0qxqMx7v/Dhta7Y98Rq8PX+2ZrPl8BZW7XeqnmZvmX2ikMytYnhF3uz5Ji2qtyjwqrOo52fRazD/KbhpHjToTONHZnvxU+xvoSHCmTHHeE1fpGnWn8zJ7sKFsoKYXNVeyRrB6NA7uHPkY1SKDiM6PNRjgf/09x/x2fZXkZCME+s0O5yrG97HUxcO9TqukhYeDwCn4zRYjwZuwunh9IbXEQRZoAoPgCkbp/Dc0ucYeMZAnuz85El/3MzsTL768yve+fUddiftpmPtjtzT9h7a1Gzj0xg8NbgLwoohK4gIjXCz198K+tV1UcOL+G3/b6zYu4KVe1fy6/5fT3zZNq3clPa12hMiIXyx+Yt81W73tr+XhpUanigk/jz8J38e+ZOkjL8n2KkWVY3j6cdP+sWfw9sv7pI2FidlJNHlky4er8Se6vIUvRr1olJEpUKPFUzvrXmPMavGAM575+7L71DqIX7a+RMLdixg0c5FpGSmEB0WTZNKTdh0eNOJKh3w7ZVXUkYSY38dy+T1k4kOj+butndz9RlXe/zBVVQF/Q+sGb7GJ+c4SXoSvH421GwOw2dy7vPfs9PNrIN1q0Qz//7u7D+Wxr5jqa77NEJ3reTyjQ8QnpXCC7EPMvHAmVwesoiHwqYSL4ns0jhezBzEzOxuJ44VGiJUigqjYlQ4laLDqBQVTsUo537O2t2kRa44qdosbX9vaoV05X+PXOj1yypR4eE6wMVAL9fiN6rqv+4OfhDIwgOcxsr3f3ufCuEVSMpIonZsbXrU78GSXUvYdnQbLeNack+7e+hcp7Nfqgp8ccnq7a+ujKwM1iWuY+XelazYu4JV+1adVBh4EhcVR9MqTWlapSmnVTmNJpWb0LRKU6pGVfVJT6GS/mourNonIiSCCxpcwGVNLqNr3a6Eh4R7fexA+f3g7wz8aiCv9XiNng17Frp9WlYay/csZ8GOBXy26bOTqj1z1Iyuyfyr5xf7c6uqzN06l/+s+A/7UvZx5WlXMrLdyJOqjHwhKN3dl7wF8x6DG2bz5aHGPDr9N1Iy/r76iQ4PZfSAVieleQecgXsz7oKKteHa/0Kt5h4Ln6ox4TzY+0yOpmZwLDWDoymZzn1qJkdTMjiWmsnR1Ax2H8lfZQZOE/rW573/Pyhx4eE6SBxwPvCXqq70+uylQKALj1l/zuKJ/z2R77K5ZnRNHuv8GBfWv9Dv3WyD1U0zMzuTdh+18/irfULvCZxW5bRCOxUEu7Ha03v4VJenaFy5MTP/nMncrXM5lHaIalHVuKTxJVzW9DLOqnbWib9tsF/DxoMbufqrq70uPHIrqLt4zZiatK3Z9kR7QbOqzdxeMeR9/YObDWbRzkWs2LuCs6qdxeOdH+fsGmcX67UVJij/Axkp8HobiGvqFCCrd/HSvN/ZdTiF+CrRPNi72ckFR3a209C+6BVo2A0GTYJYpxD9ctVO7wsfNwq68vHVlYfHa0QRmQU8oqprRaQO8AuwAmgqIuNU9TWvIyhnxqwa47a+NTQklJ4NivZPXBw5/xzB+OIKCwmjdmxtj7/6OtR2+znMp1+TfkFtmC7sPWxRvQUPdHiARTsX8dWWr/j090+ZvGEyTSs35bKmlxEVFsVrK18rs92FPf0NK0VUon2t9qzet5p52+YBEB0WTevqrU/0djq7xtksTFiYr8fca7+8RnRoNP/s/E+uOv0qQkNC8x3fV4LyPxAeDef9H8x9ELYupH/bHp6/6NOOwfTb4PfZ0G44XPIyhP1dnZyzX4GFTwEe7N3MbeHzYO+Stavm5vHKQ0TWqWoL1+PHgDNVdZiIVAT+Vx6z6nor4PWtpUxZH6BWHEfSjjBv2zy++vMrVu9f7XG7QDa4n7jyuOC1Iv9o8eZvuCdpz0njI34/9DvZmo0ghIaEuu0FVyumFvOvnl+yF1aaZaTCG+2gUl24+RtwV8NwaDtMuRb2b4Q+o6HTCPfbldCXq3YWu/DJUdxBghm5HvcE3gNQ1WMikr8y1Jzg6VebpxkKTzXBvPIJlsqRlRnUbBCDmg3ir6N/0e8L9681kFMZlyTrszd/w9qxtenbuC99G/cFnEbwNfvXsHrfat7+9W23x92XvK/YMZUJ4VHO1cfs+2Hzd3D6RSc/v30xfDoEsjNhyDRo6n0VUlHlnUbX1woqPHaIyN04g/raAV8DiEg0zkBB48HIdiPd/mrzdnTpqSDY1U7B1KBSA+rE1inzPyCK+jeMDY+lS3wXusR34YvNX5T5119sbYc6XXd/+Dec1vPvq4pfJsGs+6FqQ7j2U6h+WlDDLKmCUrLfDLQAbgAGq+ph1/rOOCPNjQfuZiI8latsTH4j2408MS4mR7B+QIjbjED+5W4q53LzAyosAro/CLt+gZdOc+bTeK4uzLwbGp8Ht8wv8wUHFHDl4UqFfrub9T8AP/gzqFNBef7lbZy//19H/zpRfeNpnIU/eZMxwF/KY9XlSULCAYFkV1r29OMQEgatBkN01aCG5isF9bY6HXgMOAS8gtPmcT6wGbhZVQPXAm1MGdStbjfe/vVt3ur5FufXOz9ocQTjygPK+Q+oH/5Nvnyv2Znww7+gzTVBCcnXCqq2mgAswUlauBQYD8QBDwBv+T80Y4wpo454mODU0/oyqKDCo4KqjlPVl4EUVf1MVVNdo8sLTgpjjAm6YFZblXuV6xVtfRlUUOGRuztu3jy+1lXXmEIEK0ttXqUljnKl55POoMHcwqOd9aeIgrrqnikia3DSoTR1Pca13MTvkRljTFnVepBz/90zTlVV5XpOwZGz/hRQUOFxVsCiMOYUVpLBeiU7cXBOa1xaDzqlCou8Cuqq63aOchE5F7gOuNNfQRlzKghWL6e8Sksc5tTiVfJ8EWmDU2AMArZSCuYPN8YYEzwFjfM4A7gGuBZIBD7FSaR4QYBiM+aUEKxeT9bbyvhTQVceG4GfgMtUdTOAiNwXkKiMORWUktoi621l/KGgrrpXAXuAH0TkPRHpSan5dzCm7Ahag7kxfuTNHOaxQH+c6qsLgQ+BL1Q1MJMS+ICI7AfcdgAoBaoDB4IdRAEsvpIr7TFafCVzKsfXUFVruHvC62loAUSkGnA1TpZd/yWiL0dEZIWnyVZKA4uv5Ep7jBZfyZTX+AqqtspHVQ+q6rtWcBhjTPlWpMLDGGOMASs8SoNxwQ6gEBZfyZX2GC2+kimX8RWpzcMYY4wBu/IwxhhTDFZ4GGOMKTIrPAJAROqLyA8iskFE1onISDfb9BCRIyKy2nULaOJ/EdkmIr+5zp1vimFxjBGRzSKyRkTaBTC2Zrnel9UiclRE7s2zTUDfPxEZLyL7RGRtrnXVRORbEfnDde92smoR6SMiv7vey0cCHONLIrLR9Tf8QkSqeNi3wM+DH+MbJSI7c/0dL/Gwr9/fQw/xfZortm0istrDvn59/zx9pwT0M6iqdvPzDagDtHM9rghsAprn2aYHMCuIMW4Dqhfw/CXAXJwsA52BpUGKMxQn80HDYL5/wPlAO2BtrnUvAo+4Hj8CvOAh/j9x5sSJAH7N+1nwc4y9gDDX4xfcxejN58GP8Y0CHvDiM+D399BdfHme/w/wZDDeP0/fKYH8DNqVRwCo6m5V/cX1+BiwAagb3KiK7Apgkjp+BqqISJ0gxNET+FM9TBkQKKr6I3Awz+orcDIw4Lrv72bXTsBmVd2iqunAf137BSRGVf1GVTNdiz8DQZsX1cN76I2AvIcFxSdOwrBBwBRfn9cbBXynBOwzaIVHgIlII6AtsNTN011E5FcRmSsiLQIbGQp8IyIrRWSEm+frAjtyLScQnALwGjz/wwbz/QOopaq7wfnnBmq62aa0vI8AN+FcTbpT2OfBn+5yVauN91DtUhrew/OAvar6h4fnA/b+5flOCdhn0AqPABKRCsDnwL2qmnde+F9wqmLOBt4AvgxweOeqajugL3CniJyf53l3STED2s9bRCKAy4HP3Dwd7PfPW0F/HwFE5HEgE/jYwyaFfR785R2gKdAG2I1TNZRXaXgPr6Xgq46AvH+FfKd43M3NuiK/f1Z4BIiIhOP8kT9W1XyTaanqUVU97no8BwgXkeqBik9Vd7nu9wFf4Fza5pYA1M+1XA/YFZjoTugL/KKqe/M+Eez3z2VvTlWe636fm22C/j6KyHDgUuB6dVWC5+XF58EvVHWvqmapajbwnofzBvU9FJEwYADOHEduBeL98/CdErDPoBUeAeCqH/0A2KCqr3jYprZrO0SkE87fJjFA8cWKSMWcxziNqmvzbDYTGCaOzsCRnMvjAPL4ay+Y718uM4HhrsfDgRlutlkOnC4ijV1XUte49gsIEekDPAxcrqrJHrbx5vPgr/hyt6Nd6eG8QX0PgYuAjaqa4O7JQLx/BXynBO4z6K/eAHY7qXdDN5zLwjXAatftEuB24HbXNncB63B6PvwMdA1gfE1c5/3VFcPjrvW54xPgLZxeGr8BHQL8HsbgFAaVc60L2vuHU4jtBjJwfsndDMQB3wF/uO6rubaNB+bk2vcSnN4xf+a81wGMcTNOfXfO53Bs3hg9fR4CFN9Hrs/XGpwvtDrBeg/dxedaPzHnc5dr24C+fwV8pwTsM2jpSYwxxhSZVVsZY4wpMis8jDHGFJkVHsYYY4rMCg9jjDFFZoWHMcaYIrPCw5giEJHHXVlM17gypp7jWh8uIs+7spmuFZFlItLXR+fsISKzfHEsY3wlLNgBGFNWiEgXnJHZ7VQ1zTWCPcL19LM4mU5bup6rBXQv5nlCVTWrBHGG6d/JD43xCys8jPFeHeCAqqYBqOoBABGJAW4FGud6bi8wNe8BRKQn8DLO/95y4A5XYbMNGI8zGvlNETkMvAYcwMnblbN/LE7urlauY4xS1RkicgPQD4gCYoELffvSjTmZVVsZ471vgPoisklE3haRnCuL04C/tJDEdCIShTM6ebCq5nz535Frk1RV7YaT1PE94DKc7K21c23zOPC9qnYELgBechUoAF2A4apqBYfxOys8jPGSOokX2wMjgP3Ap65f/N5qBmxV1U2u5Q9xJhzKkZNo70zXdn+okwJicq5tegGPuGawW4BzpdHA9dy3qlqc+TGMKTKrtjKmCFxtEQuABSLyG07yualAAxGpqM7EPJ64S4WdW1LuUxVwjKtU9feTVjoN90nudzHG9+zKwxgviTOX+um5VrUBtquTnfYDYIwrSykiUkdEhuQ5xEagkYic5loeCix0c6qNQGMRaepavjbXc/OAu3NlEG5bktdkTHFZ4WGM9yoAH4rIehFZgzNn9CjXc0/gVGWtF5G1OO0W+3PvrKqpwI3AZ66rlmxgbN6TuLYbAcwWkUVA7il3nwXCgTWu8zzrs1dnTBFYVl1jjDFFZlcexhhjiswKD2OMMUVmhYcxxpgis8LDGGNMkVnhYYwxpsis8DDGGFNkVngYY4wpsv8HGNo1cryxu8IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the convergence of the first sobol to that of the highest order\n",
    "\n",
    "sobol_first_exact = {'x1': exact['S1'], 'x2': exact['S2'], 'x3': exact['S3']}\n",
    "\n",
    "O = [R[r]['order'] for r in list(R.keys())]\n",
    "plt.figure()\n",
    "for v in list(R[O[0]]['results'].sobols_first('Ishigami').keys()):\n",
    "    plt.semilogy([o for o in O],\n",
    "                 [np.abs(R[o]['results'].sobols_first('Ishigami')[v] - sobol_first_exact[v]) for o in O],\n",
    "                 'o-',\n",
    "                 label=v)\n",
    "plt.xlabel('SC order')\n",
    "plt.ylabel('ABSerror for 1st sobol compared to analytic')\n",
    "plt.legend(loc=0)\n",
    "plt.savefig('Convergence_sobol_first.png')\n",
    "plt.savefig('Convergence_sobol_first.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-06-07T14:17:47.032599Z",
     "start_time": "2021-06-07T14:17:45.823359Z"
    },
    "code_folding": [
     0
    ]
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1000/1000 [00:00<00:00, 3661.17it/s]\n"
     ]
    }
   ],
   "source": [
    "# prepare the test data\n",
    "test_campaign = uq.Campaign(name='Ishigami.')        \n",
    "test_campaign.add_app(name=\"Ishigami\", params=define_params(), \n",
    "                      actions=uq.actions.Actions(uq.actions.ExecutePython(run_ishigami_model)))\n",
    "test_campaign.set_sampler(uq.sampling.quasirandom.LHCSampler(vary=define_vary(), count=100))\n",
    "test_campaign.execute(nsamples=1000, sequential=True).collate(progress_bar=True)\n",
    "test_df = test_campaign.get_collation_result()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-06-07T14:18:01.960202Z",
     "start_time": "2021-06-07T14:17:47.037407Z"
    },
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# calculate the SC surrogates\n",
    "if __name__ == '__main__':\n",
    "\n",
    "    test_points = test_df[test_campaign.get_active_sampler().vary.get_keys()]\n",
    "    test_results = np.squeeze(test_df['Ishigami'].values)\n",
    "    test_predictions = {}\n",
    "    for i in list(R.keys()):\n",
    "        test_predictions[i] = np.squeeze(np.array(R[i]['results'].surrogate()(test_points)['Ishigami']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-06-07T14:18:02.613973Z",
     "start_time": "2021-06-07T14:18:01.963819Z"
    },
    "code_folding": [
     0
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEJCAYAAACKWmBmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAucUlEQVR4nO3deXxU5d3//9cnCyRAIOyQALIHUIoI0optRb0VtFYprVb0V621t0trbe+vpYJ2sbW9tcXa1taltm69bV2L1B2t+9JWQISwBRBBEvYlEEKAkHx+f8wEhzAzGUhmziR5Px+P88ica87yyWGYT67rnOu6zN0RERGJJSPoAEREJL0pUYiISFxKFCIiEpcShYiIxKVEISIicSlRiIhIXEoUIiISlxKFiIjElRV0AA0xs4HAjUAnd/9KIvt069bN+/fvn9S4RERamvnz52919+71ywNJFGZ2P3AOsNndj4sonwT8DsgE/uzut7r7auByM3sy0eP379+fefPmNXXYIiItmpmtjVYeVNPTg8CkyAIzywTuBM4CRgBTzWxE6kMTEZFIgSQKd38T2F6veBywyt1Xu/t+4FHgvJQHJyIih0inm9mFwLqI9VKg0My6mtk9wGgzmxFrZzO7wszmmdm8LVu2JDtWEZFWI51uZluUMnf3bcBVDe3s7vcC9wKMHTtWQ+KKSKtWXV1NaWkpe/fuPey9nJwc+vTpQ3Z2dkLHSqdEUQr0jVjvA6wPKBYRkWattLSUvLw8+vfvj9knf4e7O9u2baO0tJQBAwYkdKx0ShRzgSFmNgAoAy4ELkrVyWcvKGPmnBLWl1dRkJ/LtIlFTB5dmKrTi4g0qb179x6WJADMjK5du3IkTfRBPR77CDAB6GZmpcBP3P0+M7sGmEPo8dj73X1JKuKZvaCMGbOKqaquAaCsvIoZs4oBEk4WSjQikm7qJ4mGymMJJFG4+9QY5c8Dz6c4HGbOKTmYJOpUVdfw8+eWUtQrj87t2pDfLpuc7Myo+zdFohERSVfp1PQUmPXlVVHLt+7ez1m/e+vgem52Jp3bZdO5fZuDyaNzuzbMXlAWNdHMnFOiRCEizZ4SBVCQn0tZlGTRrUMbbj7vOHbsqWbHnv3sqNzPjj3VlO/Zz449+ykrr2LHnv1U7DsQ9bhl5VX88sXlDOuVx/DeHRnQrT3ZmdGfSFbTlYg0NXeP2szkfmQPhipRANMmFh3SdASh2sMPvzCCs0b2bnD/8be+wvrywx9By8ow/vzWaqprQv8obTIzGNyjA8N65zG8V0eG9c5jWK+OvLNqq5quRKRJ5eTksG3bNrp27Rr1qaecnJyEj6VEwSdfxkf7F/0PJg6LmmhumTKSs0f2ZvXW3SzfUMGyjbtYvqGCd1ZtZdb7ZQe3zTCorZfg1XQlIo3Rp08fSktLoz7dVNePIlF2pFWQ5mDs2LGe6kEBj7TpaHvlfpZv3MWyDRXc/OzSmNtdctIxHFfQiWMLOzK0Z56arkQkacxsvruPPaxciSJ4J9/6atR7JG0yM8jONCr31xxcH9Y7j2MLOjGysBPHFXakqFceLxRvjFmjUbIQkUTFShRqekoDse6R3DJlJOeOKmDNtkqKy3ayZP0uFpft5LlF63nkvY+B0H0QgAP12q7UdCUiTUWJIg00dI9kYPcODOzegfOOD627O+u2V7F4/U6Ky3Zy9+sfRj1u3VNXRT3zKOqVx8Du7WmbFbsviJquRCQaNT21ALGarurXNrIyjAHd2jO0Vx7DeuaFfvbK4/01O7hh9mI1XYm0cmp6asHiNV2dPbI3H22tpGRTBSUbd1GycTeLSst5btGGg9saUP/PBTVdiUgdJYoWoKGmq6JeoaYnRhUc3Kdy3wFWbKpgxaYKrv97cdTjxuqxLiKtixJFCzF5dOER/fXfvm0Wo/t1ZnS/ztzxyqqoTVdtszJYu62SY7q2b8pQRaSZSacZ7iQg0yYWkVtvwMOsDKPGnTNuf5NbX1jO7hjDlIhIy6dEIUweXcgtU0ZSmJ+LAYX5udx2/ijevv40vjiqgHve+JBTb3udJ+ato7Z+F3IRafH01JM06IN15dz09BI+WFfOp/p04idfPJYxx3QOOiwRaWKxnnpSjUIadHzffGZdPZ7bLxjFxp17+fLd7/K9RxewYadudou0BkoUkpCMDGPKCX147fsTuObUwTy/eCOn3fYGv39lJXvrzcUhIi2Lmp7kqHy8bQ//+/wyXlyykcL8XM4Y0YOXl25iffle9ewWaabU9CRNql/XdtzztTH87Zufpra2lgffXUtZ+V6cT+bTmL2grMHjiEj6Uz8KaZTxg7thGYfPoFVVXcP0WYtYvbWSop55DO3Zgf6a4U+kWVKikEbbEGV2P4C91bX84dWVBydlys40BnbrwNBeeQztEf7ZM48Fa3dwY8RYU5rhTyS9KFFIo8Wac7wwP5dXrjuFVZt3s3JzBSUbd7NyUwULPt7BMwvXxz2mxpoSSR9KFNJosQYlnDaxiJzsTI4r7MRxhZ0O2ady3wFWbt7Nik0V/ODJRVGPq7GmRNKDbmZLo0Xr2d3QEOXt22ZxfN98Lhjbl8L83Kjb9OqU+OTvIpI8CdUozCwX6OfuJUmOR5qpIx2UMFK0GgmEhj/fuHOvEoZIwBqsUZjZF4EPgBfD68eb2dNJjktakWg1kqtOGcjOqmqm3PUOKzZVBB2iSKvWYIc7M5sPnAa87u6jw2WL3P1TKYjvqKjDXcuwZP1OLntgLlXVNfzxa2MYP6hb0CGJtGiN6XB3wN13JiEmkbiOLejEU98+mV4dc7j0/vf4xwfqwCcShEQSxWIzuwjINLMhZvZ74N0kxyUChJqhnrxqPCf068x3H/2Au15fRUscdkYknSWSKL4DHAvsA/4G7AS+m8ygIpnZQDO7z8yeTNU5Jb10apfNXy4fxxdHFfCrF0v40T8Wc6CmNuiwRFqNRBLFF9z9Rnc/Mbz8EDg3kYOb2f1mttnMFtcrn2RmJWa2ysymxzuGu69298sTOZ+0XG2zMvndV4/nylMG8vC/P+aqh+ezZ79m3RNJhUQSxYwEy6J5EJgUWWBmmcCdwFnACGCqmY0ws5Fm9my9pUeC55FWICPDmHHWcH523rG8unwzU+/9N1t37ws6LJEWL2Y/CjM7CzgbKDSzOyLe6ggk9Kecu79pZv3rFY8DVrn76vB5HgXOc/dbgHOOIHZppS45qT+9OuZw7aMLmHLXuzx42YkM7N4h6LBEWqx4NYr1wDxgLzA/YnkamNiIcxYC6yLWS8NlUZlZVzO7BxhtZjFrMmZ2hZnNM7N5W7ZsaUR40hyceWwvHvnvz7B73wG+fPe7/OblEk6+9VUGTH+Ok299VUOcizShRPpRZLt79VGfIFSjeNbdjwuvnw9MdPdvhte/Boxz9+8c7TnqUz+K1mPN1kq+cvc7bK089COam53Z4DAiInKoxvSj6G9mT5rZUjNbXbc0IpZSoG/Eeh9CtReRI9a/W3uysjIPK68bfVZEGi+RsZ4eAH4C/AY4FbiM0DA8R2suMMTMBgBlwIXARY04nrRym3ZGnw+jrLyKa/72PsN7d2RYrzyG9e5IQacczA7/+GriJJHYEkkUue7+ipmZu68FbjKztwglj7jM7BFgAtDNzEqBn7j7fWZ2DTAHyATud/clR/8rSGsXaz6MnKwMPlhXzrOLNhwsy8vJCiWNXh0p6pXH8N55rNxUwU+fWaaJk0RiSCRR7DWzDGBl+Au+DEjosVV3nxqj/Hng+YSjFIkj1nwYdfcoKvZWs2JTBcs2VLB84y5KNlYwe0EZFftiP7yniZNEPpFIovge0A64FriZ0ACBlyYxJpEjUvdlHqvpKC8nmzHHdGHMMV0O7uPulJVXUbKxgssfiv7ggyZOEglpMFG4+9zwy92E7k+IpJ0jnQ/DzOjTuR19OrejMEbTVUGMCZVEWpsGE4WZPQPUf4Z2J6E+Fn909+h3EkWaiVgTJ006tmdAEYmkl0Qej11NqDbxp/CyC9gEDA2vizRr9SdO6t0ph75dcnn4Px8zf+2OoMMTCVwiHe7edPfPRyszsyXufmxSIzwK6nAnjbV19z6+fPe77Kqq5u9Xj9cQIdIqNKbDXXcz6xdxoH5A3VRj+5soPpG00q1DWx66bBxmxqUPvMeWCg0+KK1XIoniOuBtM3vNzF4H3gKmmVl74KFkBicSpP7d2nPfpWPZUrGPyx+aS2Wcx2lFWrIGE0W4z8MQQo/Jfg8ocvfn3L3S3X+b1OhEAja6X2fuvOgEFpft5Jq/va8Jk6RVajBRmFk2cCXwI+CHwDfDZSKtwunDe3Lz5ON4rWQLP5y9WFOxSquTSIe7u4Fs4K7w+tfCZd9MVlAi6ebiTx/DhvK9/OG1VfTulMt3/2tI0CGJpEwiieJEdx8Vsf6qmS1MVkAi6eq6M4eyfmcVv/nnCnrn53DB2L4N7yTSAiRyM7vGzAbVrZjZQKAmzvYiLZKZceuUT/G5Id2YMauY10s2Bx2SSEokkii+D7xmZq+b2RvAq4SehBJpddpkZXDXxSdQ1DOPb/31fRaX7Qw6JJGki5sozCwTGEXoqadrw0uRu7+WgthE0lJeTjYPXHYindu14esPzGXd9j1BhySSVHEThbvXAOe6+z53X+TuC91dPY+k1evZMYeHvnEi1TW1XPrAe+yoVN9TabkSGcLjF0An4DGgsq7c3d9PbmhHT0N4SKrMXbOdi//8Hwo65bD/QC0bdu7VDHnSbMUawiORp57Gh3/+LKLMCc1LIdKqndi/CxeN68uD7649WKYZ8qSliZsowvconnb336QoHpFm5+Wlhz/9FJohb7kShbQIcROFu9eY2bmAEoVIDLFmwisr38uX736Xwd07MKRnBwb16MCQHh0o6JRLRoYdsu3sBWUxZ+gTCVoiTU/vmtkfaEb3KERSqSDGDHnt22SSlWH8c9kmHpu37mB5bnYmg8NJY1CPDmzfvY+H//Mx+w6ExpFS05Wkm0RuZkd7FNbdPW3vUehmtqTS7AVlh82Ql5udyS1TRh78ot9euZ9Vm3ezavNuVm6uOPh6w87YE0QW5ufyzvS0/W8mLdBR38x291OTE5JIy1CXDOI1HXVp34ZxA7owbkCXQ/at2FvNp2566bC5hiF2k5ZIqiUyZ/aPo5W7+8+ilYu0RpNHFx5VM1FeTnbMpquC/NymCE2k0RIZwqMyYqkBzgL6JzEmkVZl2sQicrMzDyv/0uiCAKIROVwiTU+/jlw3s9uAp5MWkUgrU7/pqlenHA7U1vLYvFIuGd+fHnk5AUcorV0iTz3V1w4Y2NSBiLRm9ZuuSjZWcN6db3PtIwt4+PJPk5WZSOVfJDkSmeGu2MwWhZclQAnwu+SHJtJ6FfXK43+/NJJ/r97O7S+vCDocaeUSqVGcE/H6ALDJ3TXLvEiSTTmhD3PX7OCu1z9kzDGdOX14z6BDklYqkfpsFrDR3dcSGm78W2aWn9SoRASAn3xxBMcWdOR/HvtAw5lLYBJJFH8nNMvdYOA+YADwt6RGFcHMhpvZPWb2pJldnarziqSDnOxM7r54DA5866/vs++AJpeU1EskUdSGm5qmAL919/8BeidycDO738w2m9nieuWTzKzEzFaZ2fR4x3D3Ze5+FXABcFiPQZGWrl/Xdvz6/FEUl+3k5meXBh2OtEKJJIpqM5sKXAI8Gy7LTvD4DwKTIgvCI9LeSag/xghgqpmNMLORZvZsvaVHeJ9zgbeBVxI8r0iLcuaxvbjy8wN5+N8fM3tBWdDhSCuTSKK4DDgJ+IW7f2RmA4CHEzm4u78JbK9XPA5Y5e6r3X0/8ChwnrsXu/s59ZbN4eM87e7jgYtjncvMrjCzeWY2b8uWLYmEJ9KsfH9iEeP6d2HGrGJWbqoIOhxpRRpMFO6+1N2vdfdHwusfufutjThnIbAuYr00XBaVmU0wszvM7I/A83HivNfdx7r72O7duzciPJH0lJ2Zwe8vGk37tplc9fB8Kvfp4UNJjSB68ViUsphD2Lr76+FEdaW735nEuETSXs+OOdwxdTQfba1kxqxiGhr9WaQpBJEoSoG+Eet9gPUBxCHSLI0f1I3/d8ZQnl64nof/83HQ4UgrEDNRmFmOmR3WhmNmPcysMYPPzAWGmNkAM2sDXIjGjhI5It+aMJgJRd25+ZmlLFxXHnQ40sLFq1HcAXwuSvkZJDg1qpk9AvwLKDKzUjO7PPyo7TXAHGAZ8Li7LzmysEVat4wM4zcXHE/3vLZ866/vU75nf9AhSQsWc4Y7M1vq7iNivLfE3Y9NamSNoBnupLX4YF0559/zLp8b0p0/XzL2sLm4RY7E0cxwF+8Tp6EsRdLA8X3z+dE5I/jxP5Zw/M9eomLvgagz7Ik0Rrwv/M1mNq5+oZmdCKijgkiayGubRaYZu/YewIGy8ipmzCpWxzxpMvFqFNOAx83sQWB+uGwsoR7aFyY5LhFJ0G0vraCmXhNyVXUNP31mCUN75jGoR3vaZh0+g55IomImCnd/L1yj+Dbw9XDxYuDTdT2mRSR466PMtw2wY081Z9/xFlkZxsDu7Snq1ZFhvfIY1iuPol55FObnYhZqYZ69oOzgDHtqupL64s5HEU4IPwEIP8qatjewRVqrgvxcyqIkix55bfnhOSMo2biLko0VvL92B88s/KTLUl7bLIb2yqNtljF3zQ6qa0K1krqmK0DJQoA4icLM7gF+7+5LzKwTocdca4AuZvb9uiE9RCRY0yYWMWNWMVXVnwxBnpudyQ1nD+fcUQUwquBgecXealZsqmD5xgqWb6igZGMF//pw+2FDI1RV1zBzTokShQDxaxSfCw/vDaGBAVe4+2Qz6wW8AChRiKSBui/zRJqO8nKyGXNMF8Yc0+Vg2YDpz0U9bqwmLWl94iWKyB48ZwBPALj7xrp2TRFJD5NHFx71X/+xmq4K8nMbG5a0EPEejy03s3PMbDRwMvAigJllAfoEibQQ0yYWkZt9+FNRF47rG2VraY3iJYorCQ218QDwPXffGC4/HYheVxWRZmfy6EJumTIy9BQU0LtTDnltM3l24QZNvSpAnCE8mjMN4SHSOK8t38xlD87l6gmDuH7SsKDDkRSJNYSHhuIQkcOcOqwHU8f15Y9vfMj8tTuCDkcCpkQhIlHd+IURFOTn8v0nFrJnv2bTa82UKEQkqg5ts5j5lVF8tLWSX71YEnQ4EqAGE4WZ9TSz+8zshfD6CDO7PPmhiUjQThrUlctO7s+D767h3VVbgw5HApJIjeJBQpMM1XXvXAF8L0nxiEia+cHEYQzs1p5pTy6iYm910OFIABJJFN3c/XGgFiA8Q52emRNpJXLbZHLbBaPYsLOKnz+7LOhwJACJJIpKM+sKoeFgzOwzwM6kRiUiaeWEfp25esIgHpu3jleXbwo6HEmxRBLF/wOeBgaZ2TvAX4DvJDUqEUk7154+hGG98rj+78XsqNQc3a1Jg4nC3d8HTgHGE+qtfay7L0p2YCKSXtpmZXL7BcdTvmc/P356SdDhSAol+njsOGAUcAIw1cwuSV5IIpKuRhR05LunD+GZhet5dtH6hneQFiHuxEUAZvZ/wCDgAz65ie2EmqBEpJW56pRBvLx0Ez+avZhxA7rQIy8n6JAkyRKpUYwFTnb3b7n7d8LLtckOTETSU1ZmBr++4Hj27K/hhlmLaYnjxcmhEkkUi4FeyQ5ERJqPwT06MG1iEf9ctokn55cGHY4kWbypUJ8h1MSUByw1s/eAfXXvu/u5yQ9PRNLVN04ewMtLN/GzZ5YyfnA3CjXRUYsV7x7FbSmLQkSanYwM47bzRzHxt29y/ZOL+Ms3xpGRodkvW6KYicLd3wAws1+6+/WR75nZL4E3khybiKS5vl3a8cMvjOCGp4oZffPL7KqqjjtntzRPidyjOCNK2VlNHYiINE+52RlkGOysqsaBsvIqZswqZvaCsqBDkyYSM1GY2dVmVgwUmdmiiOUjQB3uRASA215aQW29B5+qqmv41YvLgwlImly8exR/A14AbgGmR5RXuPv2pEYVwcwmADcDS4BH3f31VJ1bRBq2vrwqevnOvUyY+RpDeuYxpEcHhvbMY3CPDgzu0YGc7MxDtp29oIyZc0pYX16lpqs0FO8exU5Cg/9NPdqDm9n9wDnAZnc/LqJ8EvA7IBP4s7vfGucwDuwGcgA9hyeSZgrycymLkizycrIYUdCRlZt289ryzRwIVzvMoF+XdgzpkceQnh3YtXc/T84rY9+BWuCTpitAySJNWDI7y5jZ5wl9yf+lLlGYWSahOS3OIPTFP5dQMsokVHuJ9A1gq7vXmllP4HZ3v7ih844dO9bnzZvXdL+IiMQ0e0EZM2YVU1X9yewDudmZ3DJl5MEv+v0Halm7rZIVm3azcnMFK8M/P9paSXVN9O+gwvxc3pl+Wkp+Bwkxs/nuPrZ+eYNDeDSGu79pZv3rFY8DVrn76nBgjwLnufsthGofsewA2sZ608yuAK4A6NevX2PCFpEjUJcM4jUdtcnKCDVB9cwDeh8sr66pZeiNLxAtVcRq0pLUi5sown/9z3H3/2rCcxYC6yLWS4FPx4lhCjARyAf+EGs7d78XuBdCNYqmCFREEjN5dOFRNRNlZ2bEbLoqUAe+tBH38Vh3rwH2mFmnJjxntB45Mb/Y3X2Wu1/p7l/VjWyRlmfaxCJy693cNuC7pw8JJiA5TCJNT3uBYjN7GaisK2zEwIClQN+I9T6AxisWaaXqN1117dCGrbv389G2ygb2lFRJJFE8F16aylxgiJkNAMqAC4GLmvD4ItLM1G+6mvbEQv705momH19IUa+8ACMTSGyGu4eAR4D54eVv4bIGmdkjwL8IddorNbPL3f0AcA0wB1gGPO7umi5LRA664ezhdMzN5oaniqmt35tPUi6RiYsmAA8Bawg1HfY1s0vd/c2G9nX3qH0w3P154PkjCVREWo/O7dtw49nDue6JhTw6dx0XfVpPMgYpkbGefg2c6e6nuPvnCT2B9JvkhiUird2UEwo5aWBXbn1hGZsr9gYdTquWSKLIdveSuhV3XwFkJy8kEREwM37xpePYW13Lz59dFnQ4rVoiiWKemd1nZhPCy58I3asQEUmqgd078O1TB/P0wvW8sWJL0OG0WokkiqsJDch3LfBdYClwVTKDEhGpc9WEgQzs3p4fzi6man9NwztIk4s3zPgr4Zc/c/fb3X2Ku3/J3X/j7vti7Sci0pTaZmXyv18aybrtVfz+1ZVBh9MqxatR9DazU4BzzWy0mZ0QuaQqQBGRzwzsyvlj+nDvm6sp2VgRdDitTrzHY39MaB6KPsDt9d5zQMM6ikjK3HD2cF5ZvpkbnirmiStP0vzcKRSzRuHuT7r7WcCv3P3UeouShIikVF3fivlrd/Do3HUN7yBNJpGe2TenIhARkYaob0UwEnnqSUQkLahvRTCUKESkWVHfitSLmyjMLMPMFqcqGBGRRKhvRWo1NHFRLbDQzDQil4ikDfWtSK1E5qPoDSwxs/c4dOKic5MWlYhIAyL7VpyneSuSKpFE8dOkRyEichTUtyI1Enk89g1gOZAXXpaFy0REAqW+FanRYKIwswuA94DzgQuA/5jZV5IdmIhIIqacUMjg7u258aliBkx/jpNvfZXZC8qCDqtFSaTp6UbgRHffDGBm3YF/Ak8mMzARkUT844P1rNtRRd2EqWXlVcyYVQxwyDzccvQSSRQZdUkibBvqfyEiaWLmnBL2Hag9pKyquobpsxaxdtseRhR0ZHjvPArzczGLfg9j9oIyZs4pYX15FQX5uUybWKQkEyGRRPGimc0BHgmvfxXNdy0iaWJ9eVXU8r3Vtfz2lRV4uKrRMSeL4b07MqKgIyN6d2R4744M6dmBF4o3MmNWMVXVof4YqpEcLm6isFD6vQM4EfgsYMC97v5UCmITEWlQQX4uZVGSRWF+Li/9z+dZvrGCZRt2sXTDLpau38Wj7607mBSywk9JHaj1Q/atqq5h5pwSJYqwuInC3d3MZrv7GGBWimISEUnYtIlFh9QIAHKzM5k2sYj2bbMYc0xnxhzT+eB7NbXOmm2VoeSxfhd3vf5h1OPGqqm0Ronca/i3mZ2Y9EhERI7C5NGF3DJlZOgeBKGaxC1TRsasDWRmGIO6d+CcTxXwg0nDKMzPjbpdQYzy1iiRexSnAlea2VpCPbONUGXjU0mNTEQkQZNHFx51M1G0GklWhjFtYlFThdfsJXKP4ipgbWrCERFJrboEU/fUU9usDNyd8YO7BhxZ+jB3j7+B2fzwPYpmY+zYsT5v3rygwxCRZuijrZWc+Zs3OO/4Qm47f1TQ4aRU+Pt+bP1y3aMQEYkwoFt7vvHZATw5v5QP1pUHHU5aSCRRnEooWXxoZovMrNjMFiU7MBGRoHzntCF0z2vLTU8vobY2fqtLa5DIzeyzkh5FHGb2OeBiQrGOcPfxQcYjIi1fh7ZZXD9pGN9/YiFPLSjjy2P6BB1SoBIZPXYt0Bc4Lfx6TyL7AZjZ/Wa2uf4seWY2ycxKzGyVmU1v4PxvuftVwLPAQ4mcV0SksaaMLmRU33xufXE5u/cdCDqcQCUyeuxPgOuBGeGibODhBI//IDCp3vEygTsJ1VRGAFPNbISZjTSzZ+stPSJ2vYhPhhEREUmqjAzjpi+OYEvFPv7w6qqgwwlUIjWDLwHnEp7dzt3XE5qXokHu/iawvV7xOGCVu6929/3Ao8B57l7s7ufUW+pGrO0H7HT3XYn9WiIijTe6X2e+fEIf7n/7I9ZsrWx4hxYqkUSx30PP0DqAmbVv5DkLgcgZRkrDZfFcDjwQbwMzu8LM5pnZvC1btjQyRBGRkOsnFZGdafz8uaVBhxKYRBLF42b2RyDfzP6b0FwUf2rEOaON8xv3sQJ3/4m7v9vANve6+1h3H9u9e/dGhCci8okeHXP4zulD+OeyzbyxonX+EZrIzezbCE1S9HegCPixu/++EecsJXRzvE4fYH0jjiciklSXndyf/l3b8bNnllBdU9vwDi1MQk8vufvL7j7N3b/v7i838pxzgSFmNsDM2gAXAk838pgiIknTNiuTH50zgg+3VPLQu2uCDiflkjpTnZk9AvwLKDKzUjO73N0PANcAc4BlwOPuviSZcYiINNZpw3pwytDu/O6fK9m6e1/Q4aRUUhOFu091997unu3ufdz9vnD58+4+1N0HufsvkhmDiEhTMDN+dM4IqqpruG1OSdDhpFTCicLMss1sdL2+DSIircbgHh34+vj+PDZvHcWlO4MOJ2ViJgozu8fMjg2/7gQsBP4CLDCzqSmKT0QkrVz7X0Po2r4NP31mCQ2Nvt1SxKtRfC7i3sFlwAp3HwmMAX6Q9MhERNJQx5xspk0sYt7aHTy9sHU8sBkvUeyPeH0GMBvA3TcmMyARkXR3/pi+jCzsxC3PL2fP/pY/DlS8RFFuZueY2WjgZOBFADPLAjSZrIi0WhkZxk3njmDjrr3c9dqHQYeTdPESxZWEHmN9APheRE3idOC5ZAcmIpLOxhzThfOOL+Det1bz8bY9QYeTVDEThbuvcPdJ7n68uz8YUT7H3a9LSXQiImls+lnDyDTjF8+37HGgYk5cZGZ3xNvR3a9t+nBERJqP3p1y+fapg7jtpRW8s2orJw/uFnRISRFvhrurgMXA44TGYoo2mJ+ISKv2zc8N5P53PuLS+9+jptYpyM9l2sQiJo9uaFDs5iNeougNnA98FTgAPAb83d13pCIwEZHm4MXFG9m99wAHwnNrl5VXMWNWMUCLSRbx7lFsc/d73P1U4OtAPrDEzL6WothERNLezDkl7K85tONdVXUNP39uKTv3VAcUVdOKV6MAwMxOAKYS6kvxAjA/2UGJiDQX68uropZv3b2fUT97iR55bRnaM48hPTswtGceQ3t2YEjPPDrmZB/cdvaCMmbOKWF9eVVaNl3Fu5n9U+AcQiO8PgrMCI/8KiIiYQX5uZRFSRZd27fhis8PZMWm3azcXMGj762jqrrm4Pu9OuYwpGcHMgze/XAb1TXp23RlscYqMbNaYDVQdwXqNjTA3f1TyQ/v6IwdO9bnzZsXdBgi0grMXlDGjFnFhySB3OxMbpky8pAv+tpap6y8ihWbKijZVMHKTbtZsamCJet3RT1uYX4u70w/LenxRzKz+e4+tn55vKanAUmMR0SkRahLBg01HWVkGH27tKNvl3acPrznwfIB05+LOhd0rCatIMRMFO6+Nlq5mWUSmpUu6vsiIq3N5NGFR91MFKvpqiA/fUZKijfMeEczm2FmfzCzMy3kO4Saoy5IXYgiIi3XtIlF5GZnHlJmwHVnDAkmoCjiNT39H7CD0FSm3wSmAW2A89z9g+SHJiLS8tVvuurcPpvtldWHPXIbpHg3s4vD80/UNTdtBfq5e0UK4zsqupktIs2Vu/OVe/7Fuu17eH3aBNq1abAXQ5OJdTM73uixB3uKuHsN8FFzSBIiIs2ZmXHD2cPYXLGPP7/1UdDhAPETxSgz2xVeKoBP1b02s+jPc4mISKONOaYLk47txR/f+JAtFfuCDifuEB6Z7t4xvOS5e1bE646pDFJEpLX5waQi9h2o5Y5XVgYdStwahYiIBGRg9w5c9Ol+/O29j/lwy+5AY1GiEBFJU9eePoTc7Ex+9eLyQONQohARSVPdOrTlqlMGMmfJJuau2R5YHEoUIiJp7PLPDqRnx7b87/PLiNWdIdmUKERE0lhum0yuO6OIBR+X88LijYHEoEQhIpLmvjymD0N7duBXLy5n/4HalJ9fiUJEJM1lZhgzzhrOmm17eOS9j1N+/rRPFGY2wsweN7O7zewrQccjIhKECUXdGT+oK797ZSW79qZ2itWkJgozu9/MNpvZ4nrlk8ysxMxWmdn0Bg5zFvB7d78auCRpwYqIpDGzUK1ie+V+/vjGhyk9d7JrFA8CkyILwgMM3kkoAYwApoZrDSPN7Nl6Sw9Co9heaGYzga5JjldEJG2N7NOJyccX8Oe3PmLDztRNbJTUROHubwL1H/4dB6xy99Xuvp/QfNznuXuxu59Tb9kcXr4NTCc0gq2ISKt13ZlFuMPtL61I2TmDuEdRCKyLWC8Nl0VlZv3N7F7gL8DMONtdYWbzzGzeli1bmixYEZF00rdLOy4dfwxPvl/Ksg2pGZ81iERhUcpi9iJx9zXufoW7X+zub8fZ7l53H+vuY7t3794kgYqIpKNrTh1Cx5xsbn0hNUN7BJEoSoG+Eet9gPUBxCEi0ix1apfNNacO5o0VW3h7ZfJb5INIFHOBIWY2wMzaABcCTwcQh4hIs3XJ+GPo0zmXW15YRm1tcof2SPbjsY8QmnO7yMxKzexydz8AXAPMAZYBj7v7kmTGISLS0rTNymTaxCKWrN/FPxaWJfVcMefMbs40Z7aItAa1tc65d77NjspqXrnuFHKyMxt1vKOZM1tERNJYRoZxw9nDKSuv4qF31yTvPEk7soiIJN34Qd0Y3iuPW19YzoDpz3Hyra8ye0HTNkVlNenRREQkpWYvKGP11sqDfQzKyquYMasYgMmjY3ZROyKqUYiINGMz55Swr97Q41XVNcycU9Jk51CiEBFpxtaXRx/zKVb50VCiEBFpxgryc4+o/GgoUYiINGPTJhaRW++x2NzsUB+LpqKb2SIizVjdDeuZc0pYX15FQX4u0yYWNdmNbFCiEBFp9iaPLmzSxFCfmp5ERCQuJQoREYlLiUJEROJSohARkbiUKEREJK4WOcy4mW0B1gYdRwzdgORPSXX0FF/jKL7GUXyN09j4jnH3w+aSbpGJIp2Z2bxo472nC8XXOIqvcRRf4yQrPjU9iYhIXEoUIiISlxJF6t0bdAANUHyNo/gaR/E1TlLi0z0KERGJSzUKERGJS4kiCcysr5m9ZmbLzGyJmX03yjYTzGynmX0QXn6c4hjXmFlx+NzzorxvZnaHma0ys0VmdkIKYyuKuC4fmNkuM/tevW1Sev3M7H4z22xmiyPKupjZy2a2Mvyzc4x9J5lZSfhaTk9hfDPNbHn43+8pM8uPsW/cz0IS47vJzMoi/g3PjrFvUNfvsYjY1pjZBzH2TcX1i/qdkrLPoLtraeIF6A2cEH6dB6wARtTbZgLwbIAxrgG6xXn/bOAFwIDPAP8JKM5MYCOh57sDu37A54ETgMURZb8CpodfTwd+GSP+D4GBQBtgYf3PQhLjOxPICr/+ZbT4EvksJDG+m4DvJ/DvH8j1q/f+r4EfB3j9on6npOozqBpFErj7Bnd/P/y6AlgGJG8M4OQ4D/iLh/wbyDez3gHEcTrwobsH2oHS3d8EttcrPg94KPz6IWBylF3HAavcfbW77wceDe+X9Pjc/SV3PxBe/TfQp6nPm6gY1y8RgV2/OmZmwAXAI0193kTF+U5JyWdQiSLJzKw/MBr4T5S3TzKzhWb2gpkdm9rIcOAlM5tvZldEeb8QWBexXkowye5CYv8HDfL6AfR09w0Q+o8M9IiyTbpcx28QqiFG09BnIZmuCTeN3R+j2SQdrt/ngE3uvjLG+ym9fvW+U1LyGVSiSCIz6wD8Hfieu++q9/b7hJpTRgG/B2anOLyT3f0E4Czg22b2+XrvW5R9UvqInJm1Ac4FnojydtDXL1HpcB1vBA4Af42xSUOfhWS5GxgEHA9sINS8U1/g1w+YSvzaRMquXwPfKTF3i1J2RNdQiSJJzCyb0D/oX919Vv333X2Xu+8Ov34eyDazbqmKz93Xh39uBp4iVD2NVAr0jVjvA6xPTXQHnQW87+6b6r8R9PUL21TXHBf+uTnKNoFeRzO7FDgHuNjDDdb1JfBZSAp33+TuNe5eC/wpxnmDvn5ZwBTgsVjbpOr6xfhOSclnUIkiCcJtmvcBy9z99hjb9Apvh5mNI/RvsS1F8bU3s7y614Ruei6ut9nTwCUW8hlgZ10VN4Vi/iUX5PWL8DRwafj1pcA/omwzFxhiZgPCNaQLw/slnZlNAq4HznX3PTG2SeSzkKz4Iu95fSnGeQO7fmH/BSx399Job6bq+sX5TknNZzCZd+pb6wJ8llDVbhHwQXg5G7gKuCq8zTXAEkJPIPwbGJ/C+AaGz7swHMON4fLI+Ay4k9DTEsXA2BRfw3aEvvg7RZQFdv0IJawNQDWhv9AuB7oCrwArwz+7hLctAJ6P2PdsQk+pfFh3rVMU3ypCbdN1n8F76scX67OQovj+L/zZWkToi6t3Ol2/cPmDdZ+5iG2DuH6xvlNS8hlUz2wREYlLTU8iIhKXEoWIiMSlRCEiInEpUYiISFxKFCIiEpcShUgcZnZjeLTOReHRQT8dLs82s1vDo3YuNrP3zOysJjrnBDN7timOJdIUsoIOQCRdmdlJhHo1n+Du+8I9v9uE376Z0Iiex4Xf6wmccpTnyXT3mkbEmeWfDP4n0uSUKERi6w1sdfd9AO6+FcDM2gH/DQyIeG8T8Hj9A5jZ6cBthP6vzQWuDieWNcD9hHry/sHMyoHfAlsJjWNVt397QmNZjQwf4yZ3/4eZfR34ApADtAdOa9pfXeQTanoSie0loK+ZrTCzu8ysrsYwGPjYGxiUzcxyCPXs/aq7133RXx2xyV53/yyhAQ3/BHyR0EilvSK2uRF41d1PBE4FZoaTB8BJwKXuriQhSaVEIRKDhwYdHANcAWwBHgv/JZ+oIuAjd18RXn+I0AQ5deoGmhsW3m6lh4ZKeDhimzOB6eHZ1V4nVIPoF37vZXc/mjkeRI6Imp5E4gjfO3gdeN3MigkNvPY40M/M8jw0iUws0YZ3jlQZeao4x/iyu5ccUhi6qV4ZfReRpqUahUgMFpq7e0hE0fHAWg+NxHofcEd4NE7MrLeZ/X/1DrEc6G9mg8PrXwPeiHKq5cAAMxsUXp8a8d4c4DsRI+WObszvJHI0lChEYusAPGRmS81sEaE5im8Kv/dDQs1RS81sMaH7DFsid3b3vcBlwBPh2kgtcE/9k4S3uwJ4zszeBiKnfb0ZyAYWhc9zc5P9diIJ0uixIiISl2oUIiISlxKFiIjEpUQhIiJxKVGIiEhcShQiIhKXEoWIiMSlRCEiInEpUYiISFz/P0qM3Xrp58PpAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the convergence of the surrogate\n",
    "if __name__ == '__main__':\n",
    "    _o = []\n",
    "    _RMS = []\n",
    "    for r in R.values():\n",
    "        _RMS.append((np.sqrt((((test_predictions[r['order']] - test_results))**2).mean())))\n",
    "        _o.append(r['order'])\n",
    "\n",
    "    plt.figure()\n",
    "    plt.semilogy(_o, _RMS, 'o-')\n",
    "    plt.xlabel('SC order')\n",
    "    plt.ylabel('RMS error for the SC surrogate')\n",
    "    plt.legend(loc=0)\n",
    "    plt.savefig('Convergence_SC_surrogate.png')\n",
    "    plt.savefig('Convergence_SC_surrogate.pdf')"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all",
   "executable": " /usr/bin/env python",
   "main_language": "python",
   "notebook_metadata_filter": "-all"
  },
  "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.9.4"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}