{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Peak Location Histogram\n",
    "\n",
    "Code to take scans from the OpenFlexure Microscope and plot the location of the sharpest image in each z stack. Allows comparison between the success of different scanning / focusing methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Standard imports\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import json\n",
    "import pandas as pd\n",
    "from collections import Counter\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following commented-out cells are used with the images to generate stack data. Given our data's size and nature, we have not made it available here, and included the results instead. If running this notebook on your own files, adapt these cells as necessary and remove the lines importing our results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Link to a directory containing scan folders using the z stack procedure\n",
    "\n",
    "# link = r\"G:\\smart_stacks\"\n",
    "# position_counter = []\n",
    "\n",
    "# # Go through each z stack, and note the position of the sharpest image in the stack\n",
    "\n",
    "# scan_list = [os.path.join(link, f) for f in os.listdir(link) if os.path.isdir(os.path.join(link, f))]\n",
    "# for scan in scan_list:\n",
    "#     coords_list = [os.path.join(scan, f) for f in os.listdir(scan) if os.path.isdir(os.path.join(scan, f)) and 'use' not in f]\n",
    "#     for coords in coords_list:\n",
    "#         try:\n",
    "#             with open(os.path.join(coords, 'data.json'), \"r\") as read_file:\n",
    "#                 data = json.load(read_file)\n",
    "#             sharpnesses = data['stack_{0}'.format(len(data)-2)][\"sharpness\"]\n",
    "#             position_counter.append(np.argmax(sharpnesses[-9:])-4)\n",
    "#         except:\n",
    "#             pass\n",
    "\n",
    "# smart_stacks = Counter(position_counter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Link to a directory containing scan folders using only the standard autofocus and taking a z stack of images\n",
    "\n",
    "# link = r\"G:\\autofocus_only\"\n",
    "# autofocused_stacks = []\n",
    "\n",
    "# scan_list = [os.path.join(link, f) for f in os.listdir(link) if os.path.isdir(os.path.join(link, f))]\n",
    "\n",
    "# # Sort a folder of unsorted images in a single directory into x-y coordinate subfolders containing a z stack\n",
    "\n",
    "# for scan in scan_list:\n",
    "#     img_list = [os.path.join(scan, f) for f in os.listdir(scan) if os.path.isfile(os.path.join(scan, f))]\n",
    "#     for img in img_list:\n",
    "#         print(img)\n",
    "#         name = os.path.split(img)[1]\n",
    "#         print(name)\n",
    "#         coord_link = os.path.join(scan, name.split('_')[2] + '_' + name.split('_')[3])\n",
    "#         if not os.path.exists(coord_link):\n",
    "#             os.mkdir(coord_link)\n",
    "#         os.replace(img, os.path.join(coord_link, name))\n",
    "\n",
    "# # Go through each z stack, and note the position of the sharpest image in the stack\n",
    "\n",
    "# scan_list = [os.path.join(link, f) for f in os.listdir(link) if os.path.isdir(os.path.join(link, f))]\n",
    "# for scan in scan_list:\n",
    "#     coords_list = [os.path.join(scan, f) for f in os.listdir(scan) if os.path.isdir(os.path.join(scan, f)) and 'use' not in f]\n",
    "#     for coords in coords_list:\n",
    "#         img_list = [os.path.join(coords, f) for f in os.listdir(coords) if os.path.isfile(os.path.join(coords, f))]\n",
    "#         img_list.sort(key= lambda i: int(re.findall(r'-?\\d+',str(i))[-1]))\n",
    "#         sharpnesses = [os.path.getsize(i) for i in img_list]\n",
    "#         autofocused_stacks.append(np.argmax(sharpnesses)-4)\n",
    "\n",
    "# autofocused_stacks = Counter(autofocused_stacks)\n",
    "\n",
    "# Print the total number of z stacks from each method, and the Counter - a dictionary of the locations of the sharpest image, with values equal to the frequency.\n",
    "# For example, '-1':7 would mean that in 7 of the stacks, the sharpest image is the image immediately before the central image.\n",
    "\n",
    "# a_total = sum(autofocused_stacks.values())\n",
    "# print('Number of old stacks is {0}'.format(a_total))\n",
    "# for i in autofocused_stacks.keys():\n",
    "#     autofocused_stacks[i] = autofocused_stacks[i] / a_total * 100\n",
    "# print(autofocused_stacks)\n",
    "\n",
    "# s_total = sum(smart_stacks.values())\n",
    "# print('Number of new stacks is {0}'.format(total))\n",
    "# for i in smart_stacks.keys():\n",
    "#     smart_stacks[i] = smart_stacks[i] / total * 100\n",
    "# print(smart_stacks)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of autofocused stacks is 4100.0\n",
      "Counter({-4: 57.414634146341456, -3: 31.024390243902438, -2: 6.902439024390245, -1: 2.219512195121951, 0: 1.1951219512195121, 1: 0.6829268292682927, 2: 0.1951219512195122, 4: 0.1951219512195122, 3: 0.17073170731707318})\n",
      "Number of smart stacks is 1227.0\n",
      "Counter({1: 65.8516707416463, 0: 30.888345558272206, -1: 1.9559902200488997, 2: 1.0594947025264874, -2: 0.16299918500407498, -4: 0.08149959250203749})\n",
      "Number of IHI stacks is 59427\n",
      "Counter({-4: 15.810995002271694, -2: 13.064768539552729, -1: 12.593602234674476, 0: 12.305854241338112, -3: 12.183014454709138, 4: 11.447658471738436, 1: 9.953388190553115, 3: 6.495364060107359, 2: 6.145354805054941})\n"
     ]
    }
   ],
   "source": [
    "autofocused_stacks = Counter({-4: 2354.0, -3: 1272.0, -2: 283.0, -1: 91.0, 0: 49.0, 1: 28.0, 2: 8.0, 4: 8.0, 3: 7.0})\n",
    "a_total = sum(autofocused_stacks.values())\n",
    "print('Number of autofocused stacks is {0}'.format(a_total))\n",
    "for i in autofocused_stacks.keys():\n",
    "    autofocused_stacks[i] = autofocused_stacks[i] / a_total * 100\n",
    "print(autofocused_stacks)\n",
    "\n",
    "smart_stacks = Counter({1: 808.0, 0: 379.0, -1: 24.0, 2: 13.0, -2: 2.0, -4: 1.0})\n",
    "S_total = sum(smart_stacks.values())\n",
    "print('Number of smart stacks is {0}'.format(S_total))\n",
    "for i in smart_stacks.keys():\n",
    "    smart_stacks[i] = smart_stacks[i] / S_total * 100\n",
    "print(smart_stacks)\n",
    "\n",
    "# Data for historical health clinic stacks obtained separately, but included here to recreate plot from publication.\n",
    "\n",
    "IHI = Counter({-4: 9396, -2: 7764, -1: 7484, 0: 7313, -3: 7240, 4: 6803, 1: 5915, 3: 3860, 2: 3652})\n",
    "I_total = sum(IHI.values())\n",
    "print('Number of IHI stacks is {0}'.format(I_total))\n",
    "for i in IHI.keys():\n",
    "    IHI[i] = IHI[i] / I_total * 100\n",
    "print(IHI)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    IHI Autofocus Stacks  Bath Autofocus Stacks  Bath Smart Stacks\n",
      "-4             15.810995              57.414634           0.081500\n",
      "-3             12.183014              31.024390                NaN\n",
      "-2             13.064769               6.902439           0.162999\n",
      "-1             12.593602               2.219512           1.955990\n",
      " 0             12.305854               1.195122          30.888346\n",
      " 1              9.953388               0.682927          65.851671\n",
      " 2              6.145355               0.195122           1.059495\n",
      " 3              6.495364               0.170732                NaN\n",
      " 4             11.447658               0.195122                NaN\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS0AAAIvCAYAAAAhw414AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA/FUlEQVR4nO3deZgcVbnH8e+PAGGNEhIWUZjLJqsEiKKyKqCIKELQy2URUAEFFRRErwKGTa4LopdNIyjIckWuwGVRwSirAhI2MRJiAgkEwUwEAlkgBN77xzmTdCo9M92T7umumd/nefqZqe3UW9Xdb586VXVKEYGZWVks1+oAzMzq4aRlZqXipGVmpeKkZWal4qRlZqXipGVmpeKkZQ0l6WxJ/5QUkg5vdTztSNJuef+8tdWxlJGTVj+RdKmk8RXDYyVN6WbeaZJOrhi+XdLFNa7nq5Jel/TdPsb51vyF2q0Py+4AfA04ClgXuLovMfQHSctJOlHSXyXNlfSipEcknVkxz8mSprUwTKti+VYHYI0jScCRwLeAoyV9IyIW9GMImwBvRMT/9eM6++pU4DjgC8A9wErAVsC7WxmU9c41rYFld2A14DRgFrBf5cTuDkskLaw4lHs6/70tzzutYr7DJP1N0gJJMySdKWn5PO1S4HJgubxc5PHKNZon8nJTJR1fWP/ykr6Zp70q6RlJ51VMD0mHFJYZn9fZNbyvpIckzcu1pj9L2raHffUx4JKIuCIipkbExIi4OiK+lMs7HDgD2KBreySNzdMOknSfpNmSZkm6WdKmhfjWkvSzfKj8iqTHJX2qWiC51ne+pKclbdHH7Rk0XNMaWI4GroyIhZIuy8P1HqJtBzwIjAH+BLwOIOnDwE+Bk4FfAdsCPwICOIVUa3kIOAeoTIrHkL78xwG3kRLrDyS9HBGX5HkuAT4EnJDXORJ4T60BS1oHuCbHdg2p1rQtsLCHxZ4FdpW0XkQ8U2X61cBmwMHAO/O4OfnvUOBM4G/AMNKPxM2StoyIBZJWBu4A5uflnwA2BoZXiX0l4Mq8rvdGxNN93J7BIyL86ocXcCkwvmJ4LPAG6YtQfL0BnFwx7+3Axb2UvxawANg6D69H+pBvUjHPbqQk89bCsguBw/P/b83z7FaY5y7gl4Vxx5G+mCvm4cOBhYV5nga+Uxh3LvBE/n/jvL4Deti2AA4pjBsPXJr/3zbP01HH+7EZ8Ne8rx8HLiMlmOUr5jkZmFZDWcPz+nfMw58GXinu5yrvw9bAncDdwBoV0+vensH08uFhaz0NjKry+kcfyjoCeDQiHgWIVHv4PalRvBG2JH3BKt1BqgVsVG0BScNISbDach2SViHV7ABuXYbY/gLcAvxV0nWSjpP0tp4WiIhJpKSxPXA+sCJwMXBvril1S9KovJ4nJb0MPJUnbZD/bg/8LSJm9BL3zfnvnhHxwrJsz2DipNVar0XElOKLOg8DKhrgt83tUwslLQT2BA6TtGKe9Y2uRSqWHUL7fw6CipizFRZNjHiddHj5fuB+0qHtZEn79Fho8lBEnBcR/0HaX9sDn+humZxob80xHQG8i3T4GKTEV48b8/qWOBTu6/YMFu3+YbXa7A50ADuyZI1tW2BlFjfIz8x/31Kx7CiWTAhdZxuHFNYxEdilMG5X0uHh1GpBRcRLwIxulnsyIuaR2s8APlCtjIq4F8UsaSiwRWFdERF/johvRcQupNrcET2UWc1j+e9a+e8Clt4Pm5Pa3L4REbdHxGPAGiy5Dx8Atiie8KjibOCbwE2Sltj+Bm3PgOSG+IHhaOCOiLinOEHSjSxukJ8CTAfGSvoSMIJ0eURlp2qzSO1qH5A0EXg1H7qcDdwo6WvAtaRkNxY4J3q+rOJs4BxJfye1zb0f+BxwLEBETJF0JXBhbpS+h9RG9N6I+GEuYzzwWUl3Ai8D36CiViPpvaTEfSupgX0T4B2kBv6qJP2K1Oj/J9Lh+HqkNqzXWHzY9iSwjqT3AH8H5uX99yrwBUnnkH4s/quwD/8HOAm4QdJJpKS+ITAiIpY4MRIR35P0GvB/ksZExK/7sj2DSqsb1QbLi+oN8VO6mXcaNTbEs7gB/uhupu9LOizcJA/vQKoJzAceAXamoiE+z/NJ0hd2IRUN0cBhpNrIAuAZ4CyWbLg+nKUb4gV8JZf3GulM2vGFeVYgnWGclsueAfygYvo6pEOpl0jtgJ9jyYb4LYFfA8+REsp04LvkEwTd7Jcjgd+RksKreXuuB95TiOsq4HlSUhqbxx9ASmKvkM6Y7lplH64D/Jz0I/AKMInFJzt2o3BCJG/TK/n9qnt7BtNLeYeZmZWC27TMrFSctMysVJy0zKxUnLTMrFSctMysVAbUdVojRoyIjo6OVodhNig98MADsyJiZLPXM6CSVkdHBxMmTGh1GGaDkqTp/bEeHx6aWak4aZlZqThpmVmpOGmZWak4aZlZqThpmVmpOGmZWak4aZlZqThpmVmpOGmZWak4aZlZqThpmVmpOGmZWam0PGlJOlDSY5LmSpoqaec8fndJkyTNk3SbpA16K8vMBr6WJi1JewLfJj2EcnXSQz2fkDSC9Gy9U0jPwJtAem6fmQ1yre5P6zTg9Ii4Nw8/AyDpKGBiRFyTh8cCsyRtFhGTWhKpmbWFltW0JA0BRgMjJU2RNEPS+ZJWJj2s8pGueSNiLukpvVu2JlozaxetPDxcm/QE3wNITzkeBWxLejT5asDswvyzSYeQS5B0lKQJkiZ0dnY2NWAza71WJq35+e95EfFsRMwCvg/sDcwBhhXmHwa8XCwkIsZFxOiIGD1yZNO7pzazFmtZ0oqIF4AZQFSOzn8nAtt0jZS0KrBRHm9mg1irL3n4GfAFSWtJWgP4EnATcB2wlaQxklYCTgX+4kZ4M2t10joDuB+YDDwGPAScFRGdwBjgLOAFYAfgwFYFaWbto6WXPETEa8Ax+VWcNh7YrN+DMrO21uqalplZXZy0zKxUnLTMrFSctMysVJy0zKxUnLTMrFSctMysVJy0zKxUnLTMrFSctMysVJy0zKxUnLTMrFSctMysVJy0zKxUnLTMrFSctMysVJy0zKxUnLTMrFSctMysVJy0zKxUnLTMrFSctMysVJy0zKxUnLTMrFSctMysVJy0zKxUnLTMrFSctMysVJy0zKxUnLTMrFSctMysVJy0zKxU2iJpSdpE0iuSrqgYd5Ck6ZLmSrpe0vBWxmhm7aEtkhZwAXB/14CkLYEfA4cCawPzgAtbE5qZtZPlWx2ApAOBF4E/ARvn0QcDN0bEnXmeU4DHJK0eES+3JFAzawstrWlJGgacDny5MGlL4JGugYiYCiwANq1SxlGSJkia0NnZ2cxwzawNtPrw8AzgkoiYURi/GjC7MG42sHqxgIgYFxGjI2L0yJEjmxSmmbWLlh0eShoF7AFsW2XyHGBYYdwwwIeGZoNcK9u0dgM6gKckQapdDZG0BfBbYJuuGSVtCAwFJvd7lGbWVlqZtMYBv6gYPpGUxD4HrAXcI2ln4EFSu9e1boQ3s5YlrYiYR7qUAQBJc4BXIqIT6JT0WeBKYE1gPHBESwI1s7bS8kseukTE2MLwVcBVrYnGzNpVq88empnVxUnLzErFScvMSsVJy8xKxUnLzErFScvMSsVJy8xKxUnLzErFScvMSsVJy8xKxUnLzErFScvMSsVJy8xKxUnLzErFScvMSsVJy8xKxUnLzErFScvMSsVJy8xKxUnLzErFScvMSsVJy8xKxUnLzErFScvMSsVJy8xKxUnLzErFScvMSsVJy8xKxUnLzErFScvMSsVJy8xKpWVJS9JQSZdImi7pZUkPS/pQxfTdJU2SNE/SbZI2aFWsZtY+WlnTWh54GtgVeBNwMvBLSR2SRgDXAqcAw4EJwNWtCtTM2sfyrVpxRMwFxlaMuknSk8D2wJrAxIi4BkDSWGCWpM0iYlJ/x2pm7aNt2rQkrQ1sCkwEtgQe6ZqWE9zUPN7MBrG2SFqSVgCuBC7LNanVgNmF2WYDq1dZ9ihJEyRN6OzsbH6wZtZSLU9akpYDLgcWAJ/Po+cAwwqzDgNeLi4fEeMiYnREjB45cmRTYzWz1mtp0pIk4BJgbWBMRLyWJ00EtqmYb1VgozzezAaxVte0LgI2Bz4SEfMrxl8HbCVpjKSVgFOBv7gR3sxaeZ3WBsDRwCjgOUlz8uvgiOgExgBnAS8AOwAHtipWM2sfrbzkYTqgHqaPBzbrv4jMrAxafXhoZlYXJy0zKxUnLTMrlZa1afWrc7ptOmuME6K55ZvZIq5pmVmpOGmZWak4aZlZqThpmVmpDI6G+DLzSQSzJThpDXbNTIpOiNYETlqN0OzaUFm5lmhN4DYtMysV17TMBpoBXsN1TcvMSsVJy8xKxYeHVl4+8zkouaZlZqXipGVmpeKkZWal4qRlZqXipGVmpeKkZWal4ksezFrB96v2mZOWWTVOKm3Lh4dmVipOWmZWKk5aZlYqTlpmVipOWmZWKk5aZlYqTlpmVipOWmZWKm2dtCQNl3SdpLmSpks6qNUxmVlrtfsV8RcAC4C1gVHAzZIeiYiJLY3KzFqmbWtaklYFxgCnRMSciLgbuAE4tLWRmVkrtXNNa1NgYURMrhj3CLBr5UySjgKOyoNzJD3egHWPAGY1oJz+LrvZ5Tv21pTfXrGf2O19mRs0IpjetHPSWg14qTBuNrB65YiIGAeMa+SKJU2IiNGNLLM/ym52+Y69NeWXOfZmaNvDQ2AOMKwwbhjwcgtiMbM20c5JazKwvKRNKsZtA7gR3mwQa9ukFRFzgWuB0yWtKmlHYF/g8n5YfUMPN/ux7GaX79hbU36ZY284RbTvQyklDQd+CuwJ/Av4WkRc1dqozKyV2jppmZkVte3hoZlZNU5aZlYqTlpmVipOWmZWKk5aZlYqTlpmVipOWiUi6VJJ41sdR19J+rikqZJel3Rpq+NpV5JC0iGtjqNdOWk1UE4qUfGaLekeSXvXWc5OefmOBse3Q04Y9y9DGVMkje3DckNIFwr/ElgfOK6vMfQHSR+VdLek53MnlFMkXSlpWJ7elPfIeuek1Xh3Aevm17uBB4HrJW3U0qiSo4GLgI0kjernda9L6rnj1xHxTETM7uf110zS+0m3kP0O2Al4B3AsqdeRoS0MzQAiwq8GvYBLgfGFcasDAexXMe444GFSTxbPAb8A1s3TOvL8la/bK8sn9R82nfQlugFYu4bY3gTMBbYmJa6LqswTwCGFceOBS/P/t1eJrSNPezdwJzAfeAG4ClgrTzu8ynK75Wl7Aw8ArwIzgQuBVQsx/Hue5xXS7Vy/AdaoiOniwvwnA9MqhrcEbgFezPvgMeDQHvbVD4AJPUzv6T3aLsc3M7+/9wN7FZZfHvgmMDVv9zPAed29D8AhpN5NPt6X7RloL9e0mkjSisCRpA/mg4XJJ5ISyH6kw6Vf5PFPk24MB3gXqYayf8Vy7wTeB3wY+GAu43s1hHMIMCkiHiUlv4Nz77D12B+YBpzD4trk05LWAW4FZuSYPwJsBfxvXu7qPJ68besCf5L0DlLSvZPUg8dhwD7Aj7pWKOkI4ArgelJCeB/wW2BIHXH/DynZvZe0v75MSqzdeRbYWNK7upne03s0jLS978vx3gLcIGnTiuUvIdXcxgJbkHrofaLaiiSdBJwPfDQirunj9gwsrc6aA+lFSgYLSb+wc4A38t/9e1luW9Kv63p5eCcqajGF8mcCQyvGfRV4tobYHga+UDE8CfhMYZ4ea1p5eAowtjDPGaSEtWLFuG1yebvk4Y48vFPFPJcDfy6UtW/ebxvk4aeA83vYrtvpvaY1Gzi8jvdxFVIyDVIC+z9S7XjNinmqvkfdlPcI8I38/8Z5uQN6mD+ATwI/BP4BbFOYXtf2DLSXa1qNdx/pIRyjgNGkh3P8XNKiniEl7SbpFklPS3oZuDtPqqW72kkR8WrF8D9ID/7olqQdgM1Jh2xdLiO1cTXClsC9EbGga0REPEL6cm3Zy3J3FsbdAQjYQtJawNtItbhl8T3gYkm3SxorabueZo6IeRHxUeDfgP8kHb79J/C4pM17WlbSSEkXSpok6UVJc0jb2fXedq27t206EzgIeG/el33enoHGSavx5kfElPx6MCK+SqqFHA8gaX3g16TDrANJie2jedkVayh/QWE4SF/ynhydy/6npIWSFpK+FKMLDfLVylqhhpha6Q16iTkiziA9c+CXpMPWeyWd2VvBETEtIi6NiGNIST+Ak3pZ7FJg5zzfzqQfr4ep7b2tNJ5U4zuwSlx92p6Bwkmrf7wOrJz/f2f+//iI+GNEPM7SNaWuxFRPu01Vkt5Easg+lsU1wFGkw7c7WbK2NRN4S8WyQ0ltLsXYinFNBN6d2/C6lt2G1Pj/1x7CmwjsUhi3Kyk5TIyImaSE/4Eeylgi5mypmkdEPBERF0bEAcCpwOd6KHMpEfEC6aTJWnlUd+/RLsCFEXFDpPbDZ4ENK6Z3tW32tE0AfyCdpDhZ0ilV4lmm7Smzdn6wRVmtmBumIZ05PJD0xT87j/s76Ut5gqQrScnj1EIZ00k1iL0lXQ28Gn2/ROCQXNbPImJ+5YS8/u9JOjFST7Hjgc9KupN0tuobLF1DeBLYMdcY5wHPkxqKjwMulfQt4M2ks4B3RcRdPcT2XeBBSecCPya1e50HXBkRT+V5TgMukvRPUsP+cqRG7l9ExKwc80WSPg48BBxAquG8mLdxNeDbwK9y7G8G9gL+1l1Q+Tq01YCbSTXi1UgnCbYitTNB9+/R46STHHeTEtrpVCS2iJiS9/uFklYC7gGGkw4Du8rumvcOSR8EfiNphYg4tS/bM+C0ulFtIL1IhwaVp8FfJh0aFBu8jyWdgZpPas/ai4rLAPI8J5HaUl6ncMlDoaxD0tvYbUwPA//TzbQRwGtd8QHrADeSLqV4mvTrXWyIH02qLcyn+0seXqTikoc8vYNCQ3weX3nJQyfpcoziJQ8HkxqzXyWdNbsZeHOetgLpEoWZeb0XkBLFtDx9pRzLk6RLJmaSzu69rYd99r48z7S8zCzgj8DBhfmqvUdbA3/K+2EacEyVfbgC6eTFNFKNbQbwg4rpxUsedsjb9l992Z6B9nLPpWZWKm7TMrNScdIys1Jx0jKzUnHSMrNSGVCXPIwYMSI6OjpaHYbZoPTAAw/MioiRzV7PgEpaHR0dTJgwodVhmA1Kkqb3x3p8eGhmpeKkZWal4qRlZqXipGVmpeKkZWal4qRlZqXipGVmpeKkZWalUtfFpZJWAd5O6r0xSP0fPR4R85oQm5nZUnpNWpLWID237uPA9lWWWSjpAVJ/1ZdF6pbWzKwpuk1auW/xU0g9L65E6kb2StIDJv9FepjAcNIjkd4NfB/4lqQLgDOjjZ8gbGbl1VNNq+vpt2cDV0TEkz0VJGlD4FDS04+PIHXla2bWUD0lrdOBH8eSz9jrVkQ8AZwm6b+AzzYiODOzom6TVkT8d18KzEnuh73OaGbWBy2/5EHSgZIekzRX0lRJO+fxu+en9M6TdJukWp6+bGYDXJ+TlqQDJP1G0l8l3SrpkD6UsSfpGW5HkJ4RuAvwhKQRwLWkEwHDgQmkxySZ2SDXp6Ql6Vjg56Tnwd1AarC/VFJvjwwvOg04PSLujYg3IuKZiHgG2J/0hOFrIuIVYCywjaTN+hKvmQ0cPSYtSd1N/xzwqYg4NCK+HhEfIT0Z+JhaVyxpCOnBnyMlTZE0Q9L5klYGtiQ9nBOASE8/nprHF8s5StIESRM6OztrXb2ZlVRvNa0HJG1XZfyqpCcQV3omj6/V2qQn7XY9xnwUsC1wMukx5MXrvGaTDiGXEBHjImJ0RIweObLp3VObWYv1lrTuAv4k6Xu5BtTlRtLh4BGS9pT0ZeAbeXyt5ue/50XEsxExi3SB6t7AHGBYYf5hpMfMm9kg1mPSiogvArsCewITJX0gTzoJGA9cANwCfAu4Hjiu1hXn231mkO5hXDQ6/50IbNM1UtKqwEZ5vJkNYr02xEfEfaR7Di8Brpd0ObBqRHyOdDi4NrBKRBwREfXWhH4GfEHSWvkexy8BNwHXAVtJGiNpJeBU4C8RManO8s1sgKnp7GFELIyIs0jtTm8DJkk6NJLOiHijj+s/A7gfmAw8BjwEnBURncAY4CzgBWAH4MA+rsPMBhBFRO9zFReSPgN8h5Rwjo6IaQ2Oq09Gjx4d1Z57qNPU1PXGN+vfh2YDjaQHImJ0s9fTa01L0n6SLpd0raQTJK0YERcDW5DO6P1V0ok9XB5hZtYwvV2n9UXgV6RDQkiHc9cBRMRzEfEJ4D+ALwITJG3bxFjNzHqtaR0HXBARu0XE/qSOAPfK3dAAEBE3kmpdfwTuaVqkZmb0nrRGsORlBhNZ3PnfIhExJyK+QLo8wsysaXrrbvke4EuSJpLO4p1Kasd6rNrM+fIIM7Om6S1pHUO6aPSOPPw8cHi+F9DMrN/1mLQi4glJ2wCbACuTnrzzSr9EZmZWRa9P44l0IdfkfojFzKxX3TbE59tq+mRZljUz60lPZw+nSTpV0pq1FiZppKQzgB6f3GNm1lc9Ja2vAccCz0i6TtKRkraRtFrXDJJWl7SdpGMk3UTqU+tI4KvNDdvMBquensZzkaQrSYnrKGBfctcxkhYWlhfwBKlPrR/1obcHM7Oa9Hb28CXg7Pwsw3eRLh7dAhhJSmCdwF+B2yPigSbHambW+9lDWHQG8b78MjNrGffMYGal4qRlZqXipGVmpeKkZWal4qRlZqXipGVmpVJz0pK0UzMDMTOrRT01rTsl/S0/3MLPnzezlqgnaXXdT/hdYIak/5W0l6TmPp/LzKxCzUkrIr4bEVsAOwNXAh8EbgamSzpNUkdzQjQzW6zuhviI+GNEfApYFzia1LPDKcAUSbdK+oSkFRocp5kZsAxnD/MTeC4G9geuyGXtAfyCdPj4FUlDGhOmmVlS0w3TRflp0vsAnwY+lMu5GxgHvAp8HvgvYIP8v5lZQ9SVtCRtQkpUnwTWJj2d5zzgJxExqWLWayRdSHr6tJOWmTVMzUlL0l3Ae0kd/t0BnAD8KiIWdLPIXcBnlzlCM7MK9dS03g58HxgXEX+vYf7xwPv6FJWZWTfqSVrrRcRrtc4cEZ0sfsirmVlD1HP28K2SPtLdREkf8bVaZtZs9SSts4CTeph+AnB6X4KQtImkVyRdUTHuIEnTJc2VdL2k4X0p28wGlnqS1k7ALT1MvxXYpY9xXADc3zUgaUvgx8ChpLOU84AL+1i2mQ0g9bRprQU818P0maQEUxdJBwIvAn8CNs6jDwZujIg78zynAI9JWt2PJzMb3Oqpab0IbNTD9I2BuhKKpGGkQ8ovFyZtCTzSNRARU4EFwKZVyjhK0gRJEzo7O+tZvZmVUD1J6y7gSEnrFCfkcZ8hXRVfjzOASyJiRmH8asDswrjZwOrFAiJiXESMjojRI0e6xxyzga6ew8OzgI8AD0k6B3g4jx9FaoRfDfhWrYVJGkW6V3HbKpPnAMMK44ZRZ03OzAaempNWRDws6QDgZ8B3SE+YhnSF/Czg4xExoY517wZ0AE/lLrlWA4ZI2gL4LbBN14ySNgSGApPrKN/MBqC67j2MiJskrU/qS2uTPHoycGtEzK9z3eNIPUJ0OZGUxD5HavS/R9LOwIOkdq9r3QhvZnX38pCT0/XLuuKImEe6lAEASXOAV/KV9J2SPkvqbHBN0i1BRyzrOs2s/PrUNU0zRMTYwvBVwFWticbM2lVdnQBKOlDSHyXNlPR6ldfCZgVqZgb1dU3zFVLHfv8C7s1/zcz6VT2Hh8cC9wG796HR3cysIeo5PFwHuMIJy8xaqZ6kNQV4c5PiMDOrST1J6xzg05JWa1YwZma9qadN63VSTw6TJP0UeDKPW0JE/LxBsZmZLaWepHVpxf8ndzNPAE5aZtY09SQtP6TCzFqunhum/ZAKM2u5uq6I7yJpqKT1JK3Y6IDMzHpS720820n6A6lfq6dI/cYjaS1Jv5e0RxNiNDNbpJ7beEaRei+dRWpsX9TrQkTMlLQycBipR4ZBRaepaWXHN6P3mcwGkXpqWqcD/yD13/41Uud/lX4PvKtBcZmZVVVP0toZ+ElEzGFxr6WVngLe0pCozMy6UU/SWomlHzZRqdinu5lZw9WTtKYC2/cw/f3A35YtHDOzntWTtK4CDi2cIQwASScAewGXNzA2M7Ol1HNF/PeAPYFbgEmkhHWupJGkbmt+hx9db2ZNVnNNKyIWkJLWicB84BXSE59nAScB+0TEG80I0sysS72PEFsInJtfZmb9rk+38ZiZtUo9V8R/spb53J+WmTVTvf1pBUtfCV+80NRJy8yaZln701oe2Ag4hvS06G80Iigzs+40oj+t30u6DPgzsB1wWyMCMzOrpiEN8RHxKnAFqcZlZtY0jTx7+CqwXgPLMzNbSkOSlqR1gc+SntBjZtY09Vzy8IduJg0HNgNWJHUCaGbWNPWcPdyQpS9vCOB54Frg/Ij4U6MCMzOrpp6zhx2NXLGkoaQbrPcg1damAv8ZEb/J03cHLgDWB+4DDo+I6Y2MwczKp5W38SwPPA3sCryJ9ADYX0rqkDSCVHs7hZTQJgBXtypQM2sfdd0w3UgRMRcYWzHqJklPkjoaXBOYGBHXAEgaC8yStFlETOrvWM2sfdRc05L0hqTX63wtrKP8tUld3UwkPTzjka5pOcFNzePNbBCrp6b1c9IV71sBjwOP5fFbkJLNo8CDfQlC0grAlcBlETFJ0mpAZ2G22cDqVZY9CjgKYP311+/L6s2sROpJWlcCY4CPRcQNlRMkfYzU1fIJEVHXcw8lLZeXXQB8Po+ew9IPyhhGekjsEiJiHDAOYPTo0X5IoNkAV09D/BnAj4sJCyAiricljjPrWbkkAZcAawNjIuK1PGkisE3FfKuSbsyeWE/5Zjbw1JO03kFqV+rOFGDrOtd/EbA58JGImF8x/jpgK0ljJK0EnAr8xY3wZlZP0noB+EAP0/ei5+ciLkHSBsDRwCjgOUlz8uvgiOgkHYqelde7A3BgHbGa2QBVT5vWVcAJki4hPZlnch6/KfAVYB/g+7UWli8ULXYoWDl9POn2IDOzRepJWicDGwNHAIcDXU/eWY6UfG7M85iZNU09t/G8Cuwn6QPAx4B/y5OeAP4vIm5tfHhmZkuq+4r4nJycoMysJfp076GkjSXtKOlNjQ7IzKwndSUtSftImkq6Iv5O0n2CSFpL0hRJBzQhRjOzReq593A30vVTzwOnUXHmLyJmkq7h8mUJZtZU9dS0TiXdxLwDqZ+rontI9yaamTVNPUnrncCVEfFGN9NnAOsse0hmZt2rJ2ktR3riTndGkG56NjNrmnqS1mPAzj1M34eKPrDMzJqhnqR1CXCApE9XLBeSVpH038B7yF3EmJk1Sz1XxF8kaUfgJ8A5pCfx/A+pa+QhwM8i4sqmRGlmltV1RXxEHCLpV8AhpJuZRXpSzs8j4ldNiM/MbAk1JS1JKwMfBx6PiOtI12uZmfW7Wtu0XiUdFm7bxFjMzHpVU9LK12Y9zdL9tpuZ9at6zh5eBhyanwxtZtYS9TTE/wnYH3hY0oXA34F5xZki4s4GxWZmtpR6ktbvKv7/IemSh0rK44Ysa1BmZt2pJ2kd0bQozMxq1GPSkvQuYEpEPB8Rl/VTTGZm3eqtpnUPcCjpSTzkx9WPA86MiL81OTYDdFq3DyxqiPimH8pt5dLb2cPiN2YoqaM/d0FjZi3Rpz7izcxaxUnLzErFScvMSqWWSx72ltTVhrUK6Vqsj0saVWXeiIhzGxWcmVlRLUnroPyqdHQ38wbgpGVmTdNb0npfv0RhZlajHpNWRNzRX4GYmdXCDfFmVip1dbdsA08zr7j31fbWDG1d05I0XNJ1kuZKmi6peELAzAaZdq9pXUB6AOzawCjgZkmPRMTElkZlNfF9k9YMbZu0JK0KjAG2iog5wN2SbiDdwP21lgZnbcGHtoNT2yYtYFNgYURMrhj3CLBr5UySjgKOyoNzJD3egHWPAGY1oJz+LrvZ5Q+a2DW27oTYNrG3sPwNGlBGr9o5aa0GvFQYNxtYvXJERIyjwU+2ljQhIkY3ssz+KLvZ5Tv21pRf5tiboZ0b4uew9NN/hgEvtyAWM2sT7Zy0JgPLS9qkYtw2gBvhzQaxtk1aETEXuBY4XdKqknYE9gUu74fVN/Rwsx/Lbnb5jr015Zc59oZTRPueJZE0HPgpsCfwL+BrEXFVa6Mys1Zq66RlZlbUtoeHZmbVOGmZWak4aZlZqThpmVmpOGmZWak4aZlZqThpDQKSLpU0vtVxDCbe583jpNUG8gc8Kl6zJd0jae86y9kpL9/RoLjWlPTfkp6U9KqkTkl3SfqPRpTfx5hq3kZJG0i6TNLTOf7nJI2XtGfFPFMkjW1mzNZY7dzLw2BzF/CJ/P8awOeB6yVtHhFTWxTTr4A3kx4Z9zipC5MdgDVbEYykFeuYdwVgPPA06RF4T5E6k9yNFsVvDRIRfrX4BVwKjC+MW530HMn9KsYdBzxM6gHjOeAXwLp5Wkeev/J1e2X5pH7HppO6/LkBWLuHmN6cy9inl9hvBy4BzgRmAi8CZ5Fq8acC/wQ6gbMKyx0E3EfqbmgWcDOwacX0ru05GPg1MBe4urttrBLXqDx9q15iL5bXAQj4CTAVmA88AXwLGFpYfg/Sj828vB13ABtVe09JfU09lt+zFYEVgO8DM4BXgWeBX7T6s1iGlw8P21CuURxJ+jA/WJh8IrA1sB+wPulLAKlGsW/+/13AusD+Fcu9k/Qcyw8DH8xlfK+HMOaQugHaN/ci25MDSF/CnYAvA18nJaHVgJ1zzF+X9KGKZYaSEt12pHtLXyd1p12sTX0buBLYKpfb0zZWmgm8ARzQQw1tf2AacE4ua13SflRe/iBgc+B44Ii8fgAk7QHcAjwAvIdUA/153g9LkLQNcE+e/z8iYgHwBVLN+hBgE+CjwL3dxGmVWp01/Vr0q7yQlCjmkL5sc4D9e1luW1LtYL08vFMe7qhS/kwqagrAV4Fneyl/P1ItaAEwAfgh8P7CPLcDDxfGTQQeLYx7BPheD+sanmPfMQ935OFTCvNV3cZuyvxs3o/zgT+SEuA7C/NMAcbWUNaXgL9XDN8F3NTLezoeeD+p9vnVwvQfAn8g3//rV+0v17Tax32kQ5pRwGjSQz1+LmlRj5KSdpN0S25Yfhm4O0+qpZvbSRHxasXwP0htPN2KiOuA9YC9SO1bWwC/l3RBYdZHCsPPAX+pMm6tim0ZlZ+09GTelqe62ZY/9xRjL/H/CFiH9KyB35G66r5P0ld7W1bSkZLuk/RPSXOAswuxbQ/c2ksxWwO/Ab4REd8uTPtZnj5F0o8kjamnzW4wc9JqH/MjYkp+PRgRXyW1dxwPIGl9UtvONOBAUmL7aF62lg/7gsJwkA6DehQRr0bEHyLi7IjYEzgFOKZw9u61KmVXG7ccgKRVSF/4IB12vYt0+BpVtmVubzH2Ev+ciPh1RIyNiHeTujo6vacEIenjpB+Nq4G9STXa06ly6NeLp0htkIdIelMhroeBfyMdOi8g1bwellTsrdcKnLTa2+vAyvn/d+b/j4+IP0bE4yxdU+pKTEOaGNNj+e/IZShj87z8NyLi9oh4jHTGtJanSSzrNj5GSoxdSWRBlbJ2AR6KiO9HxAMR8XfS4WqlB4AP9LKu2aT2ujeA8ZLWqJyYE+p1EfFF0o/Q5hQe3GJLc9JqHytKWie/NpF0Culw7Lo8/e+kmsgJkv5N0sdIZ+cqTSd9QfaWtFbx170e+Rqt2yUdlg/lOiTtQzpMepJUg+ir6aSTDF+QtJGk3Uk1jVo6d6tpGyVtK+lGSZ+QtJWkDSX9O3AS8MeI6MyzPgnsKGl9SSMkLUe6vGNrSfvm+I5j6Qb/M4APSfqBpHdIerukwyW9vXKmiHiJdOJjHunQes0c31ckHSxpS0n/BnyK9CM1GetZqxvV/FrUaFt52v1lUlL4TGG+Y0lnt+aT2rP2yvPvVjHPScAzpC/A7RXlFy+pOCS9/d3GNJR0mv/PwPMsPvX/I+BtFfPdDlxcWHY8cGlh3G+BKyqGDyAl4leAh0g1jIXA4Xl6R962narEttQ2VplnBHBu3o+zSYeZk4HvAMMr5htNOkM7n8WXPKwA/Dhv90vAVaTr5qKwjg+SzgrOz+u4Ddiw2j4HVgF+T2rrW4t07dsDufw5wP3Avq3+LJbh5Z5LzaxUfHhoZqXipGVmpeKkZWal4qRlZqUyoHp5GDFiRHR0dLQ6DLNB6YEHHpgVEcty/V5NBlTS6ujoYMKECa0Ow2xQkjS9P9bjw0MzKxUnLTMrFSctMysVJy0zKxUnLTMrlQF19tAGF9XSkU0f+Zbc9uWalpmVipOWmZWKk5aZlYqTlpmVipOWmZWKk5aZlYqTlpmVipOWmZVKy5OWpAMlPSZprqSpknbO43eXNEnSPEm3SarlKcpmNsC1NGlJ2hP4Nukpw6uTHpL5hKQRwLWkpxkPByaQnvZrZoNcq2/jOQ04PSLuzcPPAEg6CpgYEdfk4bHALEmbRcSklkRqZm2hZTUtSUNID8ocKWmKpBmSzpe0MrAl8EjXvBExF5iaxxfLOUrSBEkTOjs7i5PNbIBp5eHh2qQn+R4A7AyMArYFTgZWIz2xt9Js0iHkEiJiXESMjojRI0c2vXtqM2uxViat+fnveRHxbETMAr4P7E16TPiwwvzDSI+LN7NBrGVJKyJeAGYAlZ2AdP0/Edima6SkVYGN8ngzG8RafcnDz4AvSFpL0hrAl4CbgOuArSSNkbQScCrwFzfCm1mrk9YZwP3AZOAx4CHgrIjoBMYAZwEvADsAB7YqSDNrHy295CEiXgOOya/itPHAZv0elJm1tVbXtMzM6uKkZWal4qRlZqXipGVmpeKkZWal4qRlZqXipGVmpeKkZWal4qRlZqXipGVmpeKkZWal4qRlZqXipGVmpeKkZWal4qRlZqXipGVmpeKkZWal4qRlZqXipGVmpeKkZWal4qRlZqXipGVmpeKkZWal4qRlZqXipGVmpeKkZWal4qRlZqXipGVmpeKkZWal4qRlZqXipGVmpdIWSUvSJpJekXRFxbiDJE2XNFfS9ZKGtzJGM2sPbZG0gAuA+7sGJG0J/Bg4FFgbmAdc2JrQzKydLN/qACQdCLwI/AnYOI8+GLgxIu7M85wCPCZp9Yh4uSWBmllbaGlNS9Iw4HTgy4VJWwKPdA1ExFRgAbBp/0VnZu2o1YeHZwCXRMSMwvjVgNmFcbOB1YsFSDpK0gRJEzo7O5sUppm1i5YlLUmjgD2Ac6tMngMMK4wbBix1aBgR4yJidESMHjlyZMPjNLP20so2rd2ADuApSZBqV0MkbQH8Ftima0ZJGwJDgcn9HqWZtZVWJq1xwC8qhk8kJbHPAWsB90jaGXiQ1O51rRvhzaxlSSsi5pEuZQBA0hzglYjoBDolfRa4ElgTGA8c0ZJAzayttPyShy4RMbYwfBVwVWuiMbN21eqzh2ZmdXHSMrNScdIys1Jx0jKzUnHSMrNScdIys1Jx0jKzUnHSMrNScdIys1Jx0jKzUnHSMrNScdIys1Jx0jKzUmmbXh5s4El9O5o1lmtaZlYqTlpmVipOWmZWKk5aZlYqTlpmVipOWmZWKk5aZlYqTlpmVipOWmZWKk5aZlYqTlpmVipOWmZWKk5aZlYqTlpmVipOWmZWKk5aZlYqTlpmVipOWmZWKi1LWpKGSrpE0nRJL0t6WNKHKqbvLmmSpHmSbpO0QatiNbP20cqa1vLA08CuwJuAk4FfSuqQNAK4FjgFGA5MAK5uVaBm1j5a9mCLiJgLjK0YdZOkJ4HtgTWBiRFxDYCkscAsSZtFxKT+jtXM2kfbtGlJWhvYFJgIbAk80jUtJ7ipeXxxuaMkTZA0obOzs7/CNbMWaYukJWkF4ErgslyTWg2YXZhtNrB6cdmIGBcRoyNi9MiRI5sfrJm1VMuTlqTlgMuBBcDn8+g5wLDCrMOAl/sxNDNrQy1NWpIEXAKsDYyJiNfypInANhXzrQpslMeb2SDW6prWRcDmwEciYn7F+OuArSSNkbQScCrwFzfCm1krr9PaADgaGAU8J2lOfh0cEZ3AGOAs4AVgB+DAVsVqZu2jlZc8TAfUw/TxwGb9F5GZlUGrDw/NzOripGVmpeKkZWal4qRlZqXipGVmpeKkZWal4qRlZqXipGVmpeKkZWal4qRlZqXipGVmpeKkZWal4qRlZqXipGVmpeKkZWal4qRlZqXipGVmpeKkZWal4qRlZqXipGVmpeKkZWal4qRlZqXipGVmpeKkZWal4qRlZqXipGVmpeKkZWal4qRlZqXipGVmpeKkZWal4qRlZqXS1klL0nBJ10maK2m6pINaHdNAIzXvZdYMy7c6gF5cACwA1gZGATdLeiQiJrY0KjNrmbZNWpJWBcYAW0XEHOBuSTcAhwJfa2lwNuA1u6YY0dzyB7K2TVrApsDCiJhcMe4RYNcWxdMSPswyW1I7J63VgJcK42YDq1eOkHQUcFQenCPp8QasewQwqwHl9HfZzS7fsTeo/Dp/jNoq9h5s0IAyetXOSWsOMKwwbhjwcuWIiBgHjGvkiiVNiIjRjSyzP8pudvmOvTXllzn2Zmjns4eTgeUlbVIxbhvAjfBmg1jbJq2ImAtcC5wuaVVJOwL7Ape3NjIza6W2TVrZMcDKwEzgf4DP9dPlDg093OzHsptdvmNvTflljr3hFD73amYl0u41LTOzJThpmVmpOGn1QNLvJYWkhlwaIulASY9Lmi1ppqTLJBUv61iW8g+T9ICklyTNkPSdBsa+laRbJM2S1NA2hWbeYyrp85ImSHpV0qWNKjeXPVTSJTnmlyU9LOlDDSz/CknP5vdzsqTPNKrswno2kfSKpCuaUX6jOWl1Q9LBwAoNLvaPwI4R8SZgQ9J1cmc2sPxVgONJFwvuAOwOnNigsl8Dfgl8ukHlVaq8x/Rg4CJJWzao7H+Q9vFPG1RepeWBp0l3abwJOBn4paSOBpV/NtAREcOAjwJnStq+QWVXugC4vwnlNkU7X1zaMpLeBHwT+CRwT6PKjYinC6NeBzZuYPkXVQw+I+lK4H0NKvtx4HFJDYsXmn+PaURcm9czGnjrspZXKHsuMLZi1E2SngS2B6Y1oPzKM+WRXxsBDyxr2V0kHQi8CPyJBn4Wm8k1req+BVwEPNfogiXtJGk26cr+McAPGr2OCrvQ/hfjdnePaaNqWv1G0tqk7WnYPpd0oaR5wCTgWeDXDSx7GHA68OVGldkfnLQK8i/yjsB5zSg/Iu7Oh4dvBb5LA36Rq5H0KWA08L1mlN9ANd1j2u4krQBcCVwWEZMaVW5EHEPaFzuTLrZ+tVFlA2cAl0TEjAaW2XSDPmlJOljSnPz6DXAhcFxELGxC2YtExDPAb4FfNLp8SR8jtYd8KCL6dCNsT7E3WE33mLYzScuR7tRYAHy+0eVHxOsRcTfph+5zjShT0ihgD+DcRpTXnwZ9m1ZEXEn6hUTSm4HngauVbsMfkmebIenjEXFXX8vuxvKkNoo+qVa+pL2AnwAfjohHG1l2kyy6xzQi/p7HleYeU6UPyiWkkwh7R8RrTVzdMn1eCnYDOoCn8md9NWCIpC0iYrsGraMpBn3SKpgNvKVi+G3An0kNq53LWng+I3lXRDwlaQPgLOD3y1puRfnvJyWa/SLiz40qN5ctYCiwYh5eCYiIWKbDlYiYK6nrHtPPkHqo3Rd477JFnORLPpYn/QANyXEvbERNOrsI2BzYIyLmN6hMJK0FvB+4CZhPqhX9R341wjiWrOWfSEpiDanJNVVE+NXNi/QmBrB8g8o7C5gBzM1/xwFrNjDe24CFpEOurtdvGrwvKl/TGlT2cOD6vF+eAg5q4D4ZWyXusQ0qe4Nc3iuFfX5wA8oeCdxBOrP3EvAocGSj9ks3++mKZpXfyJfvPTSzUhn0DfFmVi5OWmZWKk5aZlYqTlpmVipOWmZWKk5aZlYqbZG0JH1M0p25j6n5uX+i6/PV3V3zHJ77tirFnej1kNQhaaykDWuc/3ZJtzc5rEFB0m75c7VbxbjjJe1fZd6xje5LrBaSpjW6L7Be1ld1+xu8jml97b+r5UlL0heB64C/k/pq+jCL+5h6f6vi6mcdpK5wakpapAd+HNO0aAaXB4H35L9djgeqfWkvzvP2t/1INzf3l+Opvv1toR1u4zkRuD4iKjuX+wPwk3wjar+RNDSW8baU/hARf2t1DANFRLwE3FvjvDNIdzL0q4h4qL/X2dZafUk+6baHi2qY73DSLRPvJt1f9xKpV8r/BlYqzHsa6ZfzJdLjvv8AvLswz265vP1JNxh3Ai9W3NIQwNakW2PmkfoyOh1YrsrtFj8CniF1GzIJOKowzzrAZTneV3NZNwFrVcRRfO3Ww764Hbi9yrZ8DPgx6abvF0l9dQ0B3gncTbpNZiLwwUJ57wT+l/SFnA88TupTbOXCfENIteBn8z75A7AZVW6NId30fAPwQi7zj8DONd5OUuu+fzuplv5iXse9wF6FeTbN88wk3W7zFHAN+dasin23Wx6eVuW9uLQytkL5w4DzK97bx4EvkZ90VVjHR/O8s/LrCuDNNeyTaV0x1Ptd6Ka844DH8j57AZhAul+1t+3fmNSbxZN52SdI916uUWUduwK/I93PO5fUR9qnC9t0RcXwENJtbS+R7uPsNv52qGn9GThM0hPA/8WSncFVcznpGYj7k6rqY0k7/psV86xH6nJjBrAqcAhwp6TtY+meD84DfkPqKXOlwrTrSd30ng18EDgFeCOvs6sTtbtJz2YcS3ozP0jqLnhoRHT1yXU56T61r5C6512b1BXyKqTkeiypy9svsrjb277Upn5A6nPp30kdAJ5M+jDsQeq765k87lpJG8TibmvWBx4GLiV1CbMlcCrpcPXAivJPA76eyxpPupH8hmIQkrYD7gIeAo4kJZ7PAuMlvTciaul583p63vdvIe37l0ndwcwm7cebJe0TEV3d6dxM+nx8jpQo1gP2pvumkf1IHe09wuJeSaveLJ+PBG4GtiPtr0dJzRvfJ/2Yfb2wyA9JP1YHkRLud0i91x7W047oQS3fhWLMBwPnkH4E7iJ9dt9Buv8Tet7+t5A+v8fn9WxI2sZfU3HYLGlf4FekH6qjSft9S9J3oFpMK+fteA/px+PBavMt0uyaVA1Zf1PgLyzO6rPyBnygMN/hefpphfE3AZN7KH8I6TD4ceCHVX79ruvh1/5rhfE/IX1J3pyHTyH9em9SZb5ZLP41nwN8sYcYu2Lp8RemYv7bqV7T+mlhvgfz+J0qxr0jjzusm7KV99chpCSxZh6/Rt6OCwvzf5lCTYvUc8VjwIqF9+ExUlNALTWt3vb990g3h29cWMfjwIN5eEQu66M17PvdKsZNo8rNwxRqWsA+ednDC/NdTKp1jSis47LCfOfnz4+6i68inkuX9btQsc4Ha1hfrzdP58/JTjmWbSs+P9NItbfleltH/lzdDUwFNqrl89/yhvhINattSdXJs0i/+PsBt0g6ucoiNxeGHyXVFBaRtIek2yT9i/TBfo2UHN9epbzregjvl4XhX5D6HdoqD+8F3Ac8KWn5rhdwC7AmsEWe737gK5KOk7R17ualGYqd9U0C5kbqQK5yHKRud4BUY5T0bUlTSV+210i/4gI2ybNtTaq1XlNYx/9WDuRfzV3zfG9U7BORame71Lgtve37XYB7I2JK1wwR8TrpB29UrgX/i3QI81+SjpS0CY21CymxX1UYfwWpC59io321z+5QUs27L3r9LlRxP2n/nJe/J6vUujJJK0r6uqRJkuaTPiddfcy9veLvBsDFEfFGL0V21ZZXAd4bEVNriaPlSQsW9cx4Z0ScHBF7kKqdjwLflLRGYfbnC8Ovkt54YNGhya9JtYJPk47730mq7hYP/yC1l3Tnn90Mr5f/rkX64L5WeHV9sdfMf/+ddBh1EqlW+YykU5twouGFwvACUnvPIhGxIP9buS9+Rjp8+29gT9L+OrYw37r578zCOor7aDipxnMKS++XzwNr1Ljdve374VR/754jJcg1Iv2k70n61T8bmCzpCUmN6jNqOPB8xT6tjKFreqVqn12o/rmsRY/fhW78nHSovAPpx/V5SdeqticInU3uwoZ0GPwuFp9l7NqGrs98LScs3kH6Yb86Iorvd7faoU1rKRHxD0kXk9oANiG1e9VqDKl2tX9U9CKZk9+L1VbXQ1lrk36pK4chtQ1B+iWfSWrYrOZxgIiYSUoCx0p6O6kN4zRSW8FF3SzbL3KnePuSDu9+WDF+68KsXQliLZbsVbRYS3iRVPu4gPQFWUoNv8Bd5fa0758nneAoWof0nr6Q1/UE8Mlcu92GlDgvlDQtFrd79dXzwHBJKxYS1zoV09tKTuQ/Bn6cvxMfILVxXU1KZD05EPh5RCx67J2k1QrzdLWTrkfvfkuqTHxb0iuVn7+etLymJWndbiZtlv/W+0ScVUiNm4uSUe7Rs7dqczWfKAwfSKrBdTXm/5YU51MRMaHKa6l+ziPi8Yj4OulL1XWo0/WLu3IfYlxWQ0k1o2I3wYcXhh8lnQX6eGH8EsORHqt1FylBPFhtv9QYV2/7/g7g3ZU1BElDSLXahyJdylAZV0TEwyx+8sxWdO9Vansv7iB9h4r75GBSLbdhj59rhoh4ISKuJh2KV+6P7rZ/FZb+nBxRGJ5Maq/6TC3NIBHxXdJlTz+Q9KVa4m6HmtZfJY0nHdI9STqFvDfpcOWXEfFUneX9lnR241JJPyO1ZZ3C4l/oehyZD2XuJ53B+gypRjI7Tz+X9CW5S9K5pJrVqqREtnNE7Kv0DMXxpFPTk0hv+r6kBshbczmTSbXDT0l6nnzqvFrSa7SImC3pXuAESc+Sfik/ReGXMiJekPQD4OuSXs7btB2LH95aWXv6MnAnqV3yElItbUSef0hE1PI8w1r2/eHA7yR9k3Sq/BjS+/1hAEnvINXWrwamkJLz4aR9/Yce1v03YGdJ+5B+NGdFxLQq8/2G1CbzI0kjSTXQvXOsZ0cfHyrSTJLGkU5o3EM6StiUdOb81orZutv+35LO9D9K2p/7U+gWOyJC0vGks9h/kPQj0hHF5sBaEbHUmc2I+L6k14FzJS0XEef0uBG1tNY380VKTjcA00lnUuaSTpWfxJJnnw4n1Z42Liw/lqWvnfkCi68luZ90yv92qp9xW+qMHYvPYG1FulZofn7zzmDpa4XWIH2BniT9us4k1TSOz9OHkqrjE0k1hZdyTAcVyjmadDi0kL5fp7VHYb5LgRlVlg/gzIrhDtIX8OUc//mkL37xrNoQ0smS5/I+uZ30oQ3SE4wq17E5qfF8JikJz8jv8969fB7q2fdvJ10aMTt/dpa4Tot0KHsZ6UdhHulw7Q4qrlOj+tnDzfJ7OI/ar9N6Nr//k+n+Oq3i+3N4Ht/Ryz6ZRvWzh71+F6qUdVh+37relydJn99hNWz/iPyevpBfV5LaP6udQX1/fv+6uqB+BDiisE1XFJY5lvTjd1JP2+DulquQNJZ0rcsK0bgHIAxIkg4gnXjYJep8WlE35Y3F+9560A6Hh1YSknYg1cDuI9Vstic9uv5e0mGSWdM5aVk95pAu8TiWdFg0k9SI+5/hKrv1Ex8emlmptPySBzOzejhpmVmpOGmZWak4aZlZqThpmVmpOGmZWan8PxWObZ4uRobsAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 288x648 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Convert the 3 Counters into a pandas dataframe for easier manipulation and plotting\n",
    "\n",
    "frame = pd.DataFrame([IHI, autofocused_stacks, smart_stacks]).T\n",
    "frame = frame.rename(columns={0: \"IHI Autofocus Stacks\", 1: \"Bath Autofocus Stacks\", 2: \"Bath Smart Stacks\"})\n",
    "frame = frame.sort_index()\n",
    "print(frame)\n",
    "\n",
    "# Plot 3 stacked bar charts with the frequency of each position, separated by scan method\n",
    "\n",
    "plt.rcParams.update({'font.size': 12})\n",
    "axes = frame.plot(kind='bar', stacked=False, width=1.0, figsize = (4,9), color = ['darkorange', 'green', 'blue'], subplots=True, sharex=True, sharey=True, legend=False, rot=0)\n",
    "axes[2].set_xlabel('Sharpest image position in stack', size = 16)\n",
    "axes[1].set_ylabel('Frequency (%)', size = 18)\n",
    "\n",
    "plt.plot()"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "aa28e4f14a7ff8b2d1fa149f2498635dcaf698094ea98853a01965ff2586af7e"
  },
  "kernelspec": {
   "display_name": "Python 3.8.2 64-bit ('.venv': poetry)",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  },
  "metadata": {
   "interpreter": {
    "hash": "aa28e4f14a7ff8b2d1fa149f2498635dcaf698094ea98853a01965ff2586af7e"
   }
  },
  "orig_nbformat": 2
 },
 "nbformat": 4,
 "nbformat_minor": 2
}