pixelav2fpga.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Goal is to translate data from Pixelav (which comes in the form of 8x32 array) into something the board can understand (string of bits 00, 01, 10, 11 that corresponds to different charge thresholds). Then want to map it into array2.\n",
    "\\\n",
    "pixelav: 0:3 y, 0:15 x (8x32 --> 256 pixels). then each pixel goes to 2 bits which classify its charge.  \n",
    "\\\n",
    "We realized that the current pixelav output is actually **13x21**, not **8x32**, so we are truncating values from the top and bottom rows of the original shape and padding with zeroes where we need in order to reshape the array.  \n",
    "\\\n",
    "I'm currently using example parquet files form Eliza in order to structure this script. Thse files have 28127 clusters, each with 20 13x21 arrays (one every 200ps until 4ns)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import random\n",
    "from matplotlib.colors import ListedColormap, BoundaryNorm\n",
    "import matplotlib.colors as mcolors\n",
    "import matplotlib.patches as mpatches"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_size = (8,32)\n",
    "thresh0 = 400\n",
    "thresh1 = 1600\n",
    "thresh2 = 2400"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_bit_string(cluster_input, cluster_id=None, time_id=None, image_cluster=True):\n",
    "    \"\"\"\n",
    "    Inputs: Reshaped 8x32 array of Pixelav data (or directly a 13x21 array if selecting reshape=True)\n",
    "    If reshaping: \n",
    "        Inputs: 13x21 array from Pixelav, labelled with the time slice\n",
    "        Outputs: 8x32 array with (hopefully) the same data which can be fed into Array2. \n",
    "        Truncates original rows on top and bottom and adds zeros to the extra columns on the right in order to fit 8x32 \n",
    "        shape. To avoid accidentally truncating data, we add the values of each deleted pixel to the nearest row. For \n",
    "        example, since we cut off row 12, if there were any charge in that row it would be added to the same column \n",
    "        in row 9, which is the last column we save. A message will appear if this happens.\n",
    "    Outputs: String of bits in (hopefully) correct order (string should start with LSB, going RT-RB-LT-LB and end with \n",
    "    MSB), ready to go into Array 2\n",
    "    \"\"\"\n",
    "    def classify_pixel_value(charge):\n",
    "        if charge <= thresh0: \n",
    "            return \"00\", 1\n",
    "        elif thresh0 < charge < thresh1:\n",
    "            return \"01\", 2\n",
    "        elif thresh1 < charge < thresh2:\n",
    "            return \"10\", 3\n",
    "        else:\n",
    "            return \"11\", 4\n",
    "\n",
    "    if cluster_input.shape != (8, 32):\n",
    "        pixelav_row = np.shape(cluster_input)[0]\n",
    "        pixelav_col = np.shape(cluster_input)[1]\n",
    "\n",
    "        if image_cluster:\n",
    "            original_heatmap = np.zeros((pixelav_row, pixelav_col))\n",
    "            for row in range(pixelav_row):\n",
    "                for col in range(pixelav_col):\n",
    "                    charge = cluster_input[row, col]\n",
    "                    _, colorval = classify_pixel_value(charge)\n",
    "                    original_heatmap[row, col] = colorval\n",
    "\n",
    "            fig, ax = plt.subplots(figsize=(3, 3))\n",
    "            cmap = plt.get_cmap('viridis', 4)\n",
    "            norm = mcolors.BoundaryNorm(boundaries=np.arange(-0.5, 4.5, 1), ncolors=cmap.N)\n",
    "\n",
    "            cax = ax.imshow(original_heatmap, cmap=cmap, interpolation='none', aspect='auto', vmin=1, vmax=4)\n",
    "            legend_labels = ['00', '01', '10', '11']\n",
    "            colors = [cmap(norm(i)) for i in range(4)]\n",
    "            patches = [mpatches.Patch(color=colors[i], label=legend_labels[i]) for i in range(4)]\n",
    "            ax.legend(handles=patches, bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)\n",
    "            \n",
    "            ax.set_xticks(np.arange(-0.5, 21, 1), minor=True)\n",
    "            ax.set_yticks(np.arange(-0.5, 13, 1), minor=True)\n",
    "            ax.grid(which='minor', color='w', linestyle='-', linewidth=0.5)\n",
    "            ax.tick_params(which='minor', size=0)\n",
    "            ax.set_title('Original 13x21 Cluster, Before Reshaping')\n",
    "            plt.show()\n",
    "            \n",
    "        chop_start = 1\n",
    "        chop_stop = 9\n",
    "        overflow_top = np.sum(cluster_input[:chop_start, :], axis=0)\n",
    "        overflow_bottom = np.sum(cluster_input[chop_stop+1:, :], axis=0)\n",
    "        print(np.shape(overflow_top))\n",
    "\n",
    "        final_input = np.zeros((8, 32))\n",
    "        final_input[0, :21] = cluster_input[chop_start, :] + overflow_top\n",
    "        final_input[-1, :21] = cluster_input[chop_stop, :] + overflow_bottom\n",
    "        final_input[1:7, :21] = cluster_input[chop_start+1:chop_stop-1, :]   \n",
    "\n",
    "        if np.any(overflow_top) or np.any(overflow_bottom):\n",
    "            print(f\"Reshaping: Overflow in cluster {cluster_id} at {time_id}ps:\")\n",
    "            if np.any(overflow_top):\n",
    "                print(f\"  Top overflow: {overflow_top}\")\n",
    "            if np.any(overflow_bottom):\n",
    "                print(f\"  Bottom overflow: {overflow_bottom}\")\n",
    "\n",
    "    big_row_start= 0\n",
    "    pixel_data = []\n",
    "    heatmap = np.zeros((8,32))\n",
    "    while big_row_start < 8:\n",
    "        for column in reversed(range(32)):\n",
    "            for mini_row in range(2):\n",
    "                row = mini_row + big_row_start\n",
    "                charge = final_input[row, column]\n",
    "                label, colorval = classify_pixel_value(charge)\n",
    "                pixel_data.append(label)\n",
    "                heatmap[row, column] = colorval\n",
    "        big_row_start += 2\n",
    "    bit_string = ''.join(pixel_data)\n",
    "\n",
    "    if image_cluster:\n",
    "        fig, ax = plt.subplots(figsize=(3, 3))\n",
    "        cmap = plt.get_cmap('viridis', 4)\n",
    "        norm = mcolors.BoundaryNorm(boundaries=np.arange(-0.5, 4.5, 1), ncolors=cmap.N)\n",
    "\n",
    "        cax = ax.imshow(heatmap, cmap=cmap, interpolation='none', aspect='auto', vmin=1, vmax=4)\n",
    "        legend_labels = ['00', '01', '10', '11']\n",
    "        colors = [cmap(norm(i)) for i in range(4)]\n",
    "        patches = [mpatches.Patch(color=colors[i], label=legend_labels[i]) for i in range(4)]\n",
    "        ax.legend(handles=patches, bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)\n",
    "        \n",
    "        ax.set_xticks(np.arange(-0.5, 32, 1), minor=True)\n",
    "        ax.set_yticks(np.arange(-0.5, 8, 1), minor=True)\n",
    "        ax.grid(which='minor', color='w', linestyle='-', linewidth=0.5)\n",
    "        ax.tick_params(which='minor', size=0)\n",
    "        ax.set_title('Final Cluster Plot')\n",
    "        plt.show()\n",
    "        \n",
    "    return bit_string"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_array2(bit_string, imageArray=True, zoomImage=False):\n",
    "    ar2row = 256\n",
    "    ar2col = 16\n",
    "    array2 = np.zeros((ar2row,ar2col))\n",
    "    len_w2_f2 = (104*16)+8\n",
    "    len_bits = (8*32)*2\n",
    "    labels = []\n",
    "\n",
    "    # fake numbers for shift register and dnn parts of the array\n",
    "    w2_f2 = [0]*len_w2_f2\n",
    "    unused = [0]*((256*16)-len_w2_f2-len_bits)\n",
    "\n",
    "    i=0\n",
    "    for row in range(ar2row):\n",
    "        for col in reversed(range(ar2col)):\n",
    "            if i < len_w2_f2: \n",
    "                array2[row,col] = w2_f2[i]\n",
    "            elif len_w2_f2 <= i < len_w2_f2 + len_bits: #filling array with individual bits\n",
    "                array2[row, col] = int(bit_string[i-len_w2_f2])\n",
    "            elif i >= len_w2_f2+len_bits:\n",
    "                array2[row, col] = unused[i - (len_w2_f2+len(bit_string))]\n",
    "            i += 1\n",
    "\n",
    "    if imageArray==True:\n",
    "        heatmap = np.zeros((ar2row, ar2col//2))\n",
    "        bit_i = 0\n",
    "        ii = 0\n",
    "        for row in range(ar2row):\n",
    "            for col in reversed(range(ar2col//2)):\n",
    "                if ii < len_w2_f2//2:\n",
    "                    heatmap[row, col] = 0\n",
    "                elif len_w2_f2//2 <= ii <len_w2_f2//2 + len_bits//2: \n",
    "                    if bit_i < len_bits: \n",
    "                        label = bit_string[bit_i:bit_i+2]\n",
    "                        labels.append(label)\n",
    "                        if label == \"00\":\n",
    "                            colorval = 1\n",
    "                        elif label == \"01\":\n",
    "                            colorval = 2\n",
    "                        elif label == \"10\":\n",
    "                            colorval = 3\n",
    "                        elif label == \"11\":\n",
    "                            colorval = 4\n",
    "                        heatmap[row, col] = colorval\n",
    "                        bit_i += 2\n",
    "                elif ii>=len_w2_f2//2 + len_bits//2:\n",
    "                    heatmap[row, col]=0 \n",
    "                ii += 1\n",
    "\n",
    "        if zoomImage==True:\n",
    "            heatmap = heatmap[98:150, :]\n",
    "        fig, ax = plt.subplots(figsize=(5, 5))\n",
    "        cmap = plt.get_cmap('viridis', 5) \n",
    "        norm = mcolors.BoundaryNorm(boundaries=np.arange(-0.5, 5.5, 1), ncolors=cmap.N)\n",
    "        legend_labels = ['N/A', '00', '01', '10', '11']\n",
    "        colors = [cmap(norm(i)) for i in range(5)]\n",
    "        cax = ax.imshow(heatmap, cmap=cmap, interpolation='none', aspect='auto', vmin=0, vmax=4)\n",
    "        patches = [mpatches.Patch(color=colors[i], label=legend_labels[i]) for i in range(5)]\n",
    "        ax.legend(handles=patches, bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)\n",
    "                \n",
    "        ax.set_xticks(np.arange(-0.5,heatmap.shape[1], 1), minor=True)\n",
    "        ax.set_yticks(np.arange(-0.5,heatmap.shape[0], 1), minor=True)\n",
    "        ax.grid(which='minor', color='w', linestyle='-', linewidth=0.5)\n",
    "        ax.tick_params(which='minor', size=0)\n",
    "\n",
    "        ax.set_title('Array2 (each box is 2 bits)')\n",
    "        plt.show()\n",
    "\n",
    "        return array2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example using one of Eliza's BIB parquet files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(28127, 13) (28127, 273) (28127, 5460)\n"
     ]
    }
   ],
   "source": [
    "# using eliza's BIB parquet file as a sample pixelav output\n",
    "truth = pd.read_parquet('labelsBIBtracks.parquet')\n",
    "recon2d = pd.read_parquet('recon2DBIBtracks.parquet')\n",
    "recon3d = pd.read_parquet('recon3DBIBtracks.parquet')\n",
    "print(truth.shape, recon2d.shape, recon3d.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEpCAYAAABx406aAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA1iUlEQVR4nO3deVhUZfsH8O+wzLBvKgIqioriFhiG4YrrSIqauZuBmVqhRqSmFSKmkktGGYkrGi6Zvi5pr5iRytubWmKUphIaKqmAYoCALMLz+8Mf8zIOy7DO4vdzXXNd8px7zrnPHOfcc56zPBIhhAARET3VDDSdABERaR6LARERsRgQERGLARERgcWAiIjAYkBERGAxICIisBgQERFYDIiICBouBkuWLIFEIqnVe7dt2waJRILr16/Xb1LlXL9+HRKJBNu2bWuwZegrfnZ198svv6BXr14wNzeHRCJBYmKiplPSmICAAFhYWGhk2Y2xr9EGtSoGf/zxB15++WW0aNECMpkMTk5OmDJlCv7444/6zk9nLF++HCNHjkTz5s0hkUiwZMmSCuMOHDgAuVwOJycnyGQytGzZEmPHjsXFixdrtdzU1FSEhYXBy8sLtra2aNq0KXx8fPD999+rxN65cwcLFy7EgAEDYGlpCYlEgpMnT9Z4mSdPnsSYMWPg4OAAqVQKe3t7+Pn5Yf/+/bVah9q4dOkSlixZohVfUIlEovQyNzdH586dsWzZMuTn59dqnsXFxRg3bhzu37+PTz75BDExMWjdunU9Z15zZTvGspeRkRFatGiBgIAA3Lp1S9PpUV2IGvrXv/4lpFKpcHBwEO+//77YvHmz+OCDD4Sjo6OQSqVi//79as+ruLhYPHz4sKYpCCGEePTokXj48KEoLS2t1fvVkZKSIgCI6OjoamMBCAcHByGXywUAERoaWmFcWFiYmDBhgvjoo4/E5s2bxbJly0Tbtm2FqampSExMrHGO69atE6ampmLSpEni888/FxEREeLZZ58VAMTWrVuVYk+cOCEACFdXV+Ht7S0AiBMnTtRoeYsXL1bMY/HixWLLli1i1apVwsfHRwAQO3fuFELU7LOrjb1799Yq/4YAQAwZMkTExMSImJgYsX79ejF58mQBQIwdO7ZW87x8+bIAIDZt2lTP2dZNdHS0ACCWLl0qYmJixKZNm8T06dOFoaGhaNeuXa2/z9Xx9/cX5ubmDTLv6jTGvkYb1KgYXL16VZiZmQk3NzeRkZGhNO3u3bvCzc1NmJubi2vXrlU5n9zc3JpnqgE12aGlpKQIIR5/DlUVg4qkpaUJIyMjMWvWrBrnePHiRXH37l2ltoKCAuHm5iZatmyp1J6TkyMyMzOFELXbmZa9Z+zYsaKoqEhlemxsrDh8+LAQQneLQW3+bwIQgYGBKu1jx44VBgYGtdpBnjp1SgAQe/furfF7K1Mf37uyYvDLL78otb/77rsCgNizZ0+dl1ERTRaDp0WNuolWr16N/Px8bNy4Ec2aNVOa1rRpU2zYsAF5eXlYtWqVor3svMClS5cwefJk2Nraok+fPkrTynv48CHmzp2Lpk2bwtLSEiNHjsStW7dUul4q6sdr06YNRowYgR9//BFeXl4wMTFB27Zt8eWXXyot4/79+5g3bx66desGCwsLWFlZwdfXF7/99ltNPg4lbdq0qfV77e3tYWZmhqysLEVbaGgoDAwMEBcXpxQ7c+ZMSKVSRa5dunRB06ZNlWJkMhleeOEF/P3333jw4IGi3dLSEnZ2drXOMyQkBHZ2dti6dSuMjY1VpsvlcowYMaLS9/v4+MDHx0elPSAgQOXz++qrr+Dp6QlLS0tYWVmhW7du+PTTTwE83vbjxo0DAAwYMEDRZVG+y+vo0aPo27cvzM3NYWlpieHDh6t0Y5b1Q1+7dg0vvPACLC0tMWXKFDU/jeo5ODgoulLKO3v2LIYNGwZra2uYmZmhf//++O9//6uUV//+/QEA48aNg0QiUfrcfvjhB8W62djYYNSoUbh8+bLSMqr63gHAjh074OnpCVNTU9jZ2WHixIlITU2t9br27dsXAHDt2jWl9itXrmDs2LGws7ODiYkJevTogW+++UYppri4GGFhYXB1dYWJiQmaNGmCPn364Pjx4yrLuXXrFkaPHg0LCws0a9YM8+bNQ0lJiVLMmjVr0KtXLzRp0gSmpqbw9PTEvn37VOYlkUgwe/Zs7Ny5Ex07doSJiQk8PT0RHx+vFFeXfQ0A/P777+jfvz9MTU3RsmVLLFu2DNHR0Vp3HqJGxeDw4cNo06aNYsM/qV+/fmjTpg2+/fZblWnjxo1Dfn4+VqxYgRkzZlS6jICAAKxbtw4vvPACVq5cCVNTUwwfPlztHK9evYqxY8diyJAh+Pjjj2Fra4uAgAClHcFff/2FgwcPYsSIEVi7di3mz5+PCxcuoH///rh9+7bay6qLrKws3L17FxcuXMBrr72GnJwcDBo0SDH9gw8+gIeHB6ZPn67YoR87dgybNm3C4sWL4e7uXuX809LSYGZmBjMzs3rJNzk5GVeuXMHo0aNhaWlZL/OszPHjxzFp0iTY2tpi5cqV+Oijj+Dj46PYYfbr1w9z584FALz33nuIiYlBTEwMOnXqBACIiYnB8OHDYWFhgZUrVyIkJASXLl1Cnz59VL58jx49glwuh729PdasWYOXXnqpVjkXFBTg3r17uHfvHm7cuIFdu3Zh+/btmDx5slIx+OGHH9CvXz/k5OQgNDQUK1asQFZWFgYOHIiff/4ZADBr1iy89957AIC5c+ciJiYG77//PgDg+++/h1wuR0ZGBpYsWYLg4GD89NNP6N27d4U7loq+d8uXL8crr7wCV1dXrF27FkFBQYiLi0O/fv2UfpDURNmybW1tFW1//PEHnn/+eVy+fBkLFy7Exx9/DHNzc4wePRoHDhxQxC1ZsgRhYWEYMGAAPv/8c7z//vtwdnbG+fPnlZZRUlICuVyOJk2aYM2aNejfvz8+/vhjbNy4USnu008/Rffu3bF06VKsWLECRkZGGDduXIX7pVOnTiEoKAgvv/wyli5diszMTAwbNkytc3jq7Gtu3bqFAQMG4I8//sCiRYvw9ttvY+fOnYofNlpF3UOIrKwsAUCMGjWqyriRI0cKACInJ0cIIURoaKgAICZNmqQSWzatTEJCggAggoKClOICAgJUul7KDlfLumeEEKJ169YCgIiPj1e0ZWRkCJlMJt555x1FW0FBgSgpKVFaRkpKipDJZGLp0qVKbahhV4e63UQdO3YUAAQAYWFhIT744AOVnC5cuCCkUql47bXXxD///CNatGghevToIYqLi6ucd3JysjAxMRFTp06tNKam3SyHDh0SAMQnn3yiVnxFn13//v1F//79VWL9/f1F69atFX+/9dZbwsrKSjx69KjG+T948EDY2NiIGTNmKLWnpaUJa2trpXZ/f38BQCxcuFCtdapM2XZ88jV69GhRUFCgiCstLRWurq5CLpcr9T/n5+cLFxcXMWTIEEVb2fmdJ7uJPDw8hL29vaK7TwghfvvtN2FgYCBeeeUVRVtl37vr168LQ0NDsXz5cqX2CxcuCCMjI5X2J5V9777//ntx9+5dkZqaKvbt2yeaNWsmZDKZSE1NVcQOGjRIdOvWTeUz6NWrl3B1dVW0ubu7i+HDh1e53LJtVf77KYQQ3bt3F56enkpt+fn5Sn8XFRWJrl27ioEDByq1l22nc+fOKdpu3LghTExMxIsvvqiyzrXZ18yZM0dIJBLx66+/KtoyMzOFnZ2dyjw1Te0jg7Jfp9X9KiybnpOTo9T++uuvV7uM2NhYAMCbb76p1D5nzhx100Tnzp2VjlyaNWuGjh074q+//lK0yWQyGBg8XvWSkhJkZmbCwsICHTt2VPk10lCio6MRGxuLL774Ap06dcLDhw9VDne7du2KsLAwbN68GXK5HPfu3cP27dtVuh3Ky8/Px7hx42BqaoqPPvqo3vIt254NfVQAADY2NsjLy6uwm6A6x48fR1ZWFiZNmqT4pX7v3j0YGhqiZ8+eOHHihMp73njjjTrnPGrUKBw/fhzHjx/HoUOHsGjRIsTGxmLy5MkQ/z9+VGJiIpKTkzF58mRkZmYqcsvLy8OgQYMQHx+P0tLSSpdx584dJCYmIiAgQKm775lnnsGQIUPw73//W+U9T37v9u/fj9LSUowfP17p83FwcICrq2uFn09FBg8ejGbNmqFVq1YYO3YszM3N8c0336Bly5YAHnfF/vDDDxg/fjwePHigWE5mZibkcjmSk5MVVx/Z2Njgjz/+QHJycrXLfXJ9+vbtq/TdBgBTU1PFv//55x9kZ2ejb9++FX63vb294enpqfjb2dkZo0aNwrFjx1S+j09SZ18TGxsLb29veHh4KNrs7OzqtTuyvlS+V3lC2U6gfB90RSorGi4uLtUu48aNGzAwMFCJbd++vbppwtnZWaXN1tYW//zzj+Lv0tJSfPrpp/jiiy+QkpKitNGbNGmi9rLqwtvbW/HviRMnKro41qxZoxQ3f/58fPXVV/j555+xYsUKdO7cudJ5lpSUYOLEibh06RKOHj0KJyenesvXysoKQPXbvz68+eab+Prrr+Hr64sWLVpg6NChGD9+PIYNG1bte8t2KAMHDqxwetl6lDEyMlLswOqiZcuWGDx4sOLvkSNHokmTJpg3bx6OHDkCPz8/RW7+/v6Vzic7O1upq6W8GzduAAA6duyoMq1Tp044duwY8vLyYG5urmh/8ruUnJwMIQRcXV0rXEZF54IqEhkZiQ4dOiA7Oxtbt25FfHw8ZDKZYvrVq1chhEBISAhCQkIqnEdGRgZatGiBpUuXYtSoUejQoQO6du2KYcOGYerUqXjmmWeU4k1MTFTOVT753QaAI0eOYNmyZUhMTERhYaGivaJ7mir6HDp06ID8/HzcvXsXDg4OlX4G6uxrbty4ofRdL1OTfVpjUbsYWFtbw9HREb///nuVcb///jtatGih8qUrX60bkqGhYYXtotzonitWrEBISAheffVVfPjhh7Czs4OBgQGCgoKq/GXWUGxtbTFw4EDs3LlTpRj89ddfip3IhQsXqpzPjBkzcOTIEezcubPSnWFtubm5qZVDVSQSidJ2KPPkLzB7e3skJibi2LFjOHr0KI4ePYro6Gi88sor2L59e5XLKNt+MTExFX6RnzyqKn+UWN/KzgHFx8fDz89Pkdvq1auVfimWV983Vj35vSstLYVEIsHRo0cr/K6ou3wvLy/06NEDADB69Gj06dMHkydPRlJSEiwsLBTrOm/ePMjl8grnUbZD7NevH65du4ZDhw7hu+++w+bNm/HJJ58gKioKr732miK+su92ef/5z38wcuRI9OvXD1988QUcHR1hbGyM6Oho7Nq1S611U5c6+xpdonYxAIARI0Zg06ZN+PHHH5WuTCjzn//8B9evX8esWbNqlUzr1q1RWlqKlJQUpYp99erVWs2vMvv27cOAAQOwZcsWpfasrCyVK3May8OHD5Gdna3UVlpaioCAAFhZWSEoKAgrVqzA2LFjMWbMGJX3z58/H9HR0YiIiMCkSZPqPb8OHTqgY8eOOHToED799NNa7bRsbW1VDumB//3iLU8qlcLPz0+xE33zzTexYcMGhISEoH379pXeud6uXTsAjwtK+V/qmvDo0SMAQG5uLoD/5WZlZVWr3MpuOktKSlKZduXKFTRt2lTpqKAi7dq1gxACLi4u6NChQ41zqIihoSHCw8MVJ4AXLlyItm3bAnh8pKHOutrZ2WHatGmYNm0acnNz0a9fPyxZskSpGKjjX//6F0xMTHDs2DGlI5Xo6OgK4yvqmvrzzz9hZmamchRSG61bt65w/1Xf+7T6UKOfRPPnz4epqSlmzZqFzMxMpWn379/H66+/DjMzM8yfP79WyZT9gvjiiy+U2tetW1er+VXG0NBQpXrv3bu3Ue6gzMjIUGm7fv064uLiFL+0yqxduxY//fQTNm7ciA8//BC9evXCG2+8gXv37inFrV69GmvWrMF7772Ht956q8FyDwsLQ2ZmJl577TXFjq687777DkeOHKn0/e3atcOVK1dw9+5dRdtvv/2mdFklAJX/WwYGBooug7LD/rKd3pNXv8jlclhZWWHFihUoLi5WyaH8shva4cOHAUBx5ZenpyfatWuHNWvWKApETXJzdHSEh4cHtm/frrTeFy9exHfffYcXXnih2pzGjBkDQ0NDhIWFqXwHhBAqn726fHx84OXlhYiICBQUFMDe3h4+Pj7YsGED7ty5oxJffl2fXKaFhQXat2+v1MWjLkNDQ0gkEqWjzevXr+PgwYMVxp8+fVrpXEJqaioOHTqEoUOHqnUkUh25XI7Tp08rPUrk/v372LlzZ53nXd9qdGTg6uqK7du3Y8qUKejWrRumT58OFxcXXL9+HVu2bMG9e/ewe/duxS+gmvL09MRLL72EiIgIZGZm4vnnn8epU6fw559/Aqi4z682RowYgaVLl2LatGno1asXLly4gJ07dyp+zdRGTEwMbty4oXj8QHx8PJYtWwYAmDp1quJXXbdu3TBo0CB4eHjA1tYWycnJ2LJlC4qLi5VO+F6+fBkhISEICAiAn58fgMfXO3t4eCj61IHHj7dYsGABXF1d0alTJ+zYsUMpryFDhqB58+aKv8tyKrv8LSYmBj/++COAx5ezVmXChAm4cOECli9fjl9//RWTJk1C69atkZmZidjYWMTFxVV5KP7qq69i7dq1kMvlmD59OjIyMhAVFYUuXbooXXDw2muv4f79+xg4cCBatmyJGzduYN26dfDw8FCcW/Hw8IChoSFWrlyJ7OxsyGQyDBw4EPb29li/fj2mTp2KZ599FhMnTkSzZs1w8+ZNfPvtt+jduzc+//zzKtcTePzIjQEDBiA0NLTSR4uU9+effyo++/z8fJw5cwbbt29H+/btMXXqVACPi9rmzZvh6+uLLl26YNq0aWjRogVu3bqFEydOwMrKSlFAKrN69Wr4+vrC29sb06dPx8OHD7Fu3TpYW1urlWe7du2wbNkyLFq0CNevX1dcKpySkoIDBw5g5syZmDdvXrXzqcj8+fMxbtw4bNu2Da+//joiIyPRp08fdOvWDTNmzEDbtm2Rnp6O06dP4++//1bcK9O5c2f4+PjA09MTdnZ2OHfuHPbt24fZs2fXOIfhw4dj7dq1GDZsGCZPnoyMjAxERkaiffv2FXZxd+3aFXK5HHPnzoVMJlP8EA0LC6vVZ/CkBQsWYMeOHRgyZAjmzJkDc3NzbN68Gc7Ozrh//3697dPqRW0uQfr999/FpEmThKOjozA2NhYODg5i0qRJ4sKFCyqxZZe4PXmXbPlp5eXl5YnAwEBhZ2cnLCwsxOjRo0VSUpIAID766CNFXGWXe1V0idqTlzQWFBSId955Rzg6OgpTU1PRu3dvcfr0aZW4mlxa2r9//0ovMSx/+WNoaKjo0aOHsLW1FUZGRsLJyUlMnDhR/P7774qYR48eieeee060bNlSZGVlKS3n008/VbrTs+wzVGfZQlR+GWRN/ivExcWJUaNGCXt7e2FkZCSaNWsm/Pz8xKFDh6r97Hbs2CHatm0rpFKp8PDwEMeOHVO5tHTfvn1i6NChwt7eXkilUuHs7CxmzZol7ty5ozSvTZs2ibZt2wpDQ0OVdT1x4oSQy+XC2tpamJiYiHbt2omAgAClywiruqv18OHDAoCIioqq9vN48nM0NDQULVu2FDNnzhTp6ekq8b/++qsYM2aMaNKkiZDJZKJ169Zi/PjxIi4uTil/VHIH8vfffy969+4tTE1NhZWVlfDz8xOXLl1SiqnqeyfE48fK9OnTR5ibmwtzc3Ph5uYmAgMDRVJSUpXrWtkdyEIIUVJSItq1ayfatWunuCz42rVr4pVXXhEODg7C2NhYtGjRQowYMULs27dP8b5ly5YJLy8vYWNjI0xNTYWbm5tYvny50l3ulW2rivYhW7ZsEa6urkImkwk3NzcRHR1dYRz+/87xHTt2KOK7d++u8p2py75GiMfbu2/fvkImk4mWLVuK8PBw8dlnnwkAIi0tTWUemiIRQvvPdiQmJqJ79+7YsWOHVl6SRfpnwYIF2L17N65evarU90z6QyKRIDAwUK0jxfoWFBSEDRs2IDc3t166o+qD1o1n8PDhQ5W2iIgIGBgYoF+/fhrIiJ5GJ06cQEhICAsB1dmT+7TMzEzExMSgT58+WlMIgBqeM2gMq1atQkJCAgYMGAAjIyPFpYUzZ85Eq1atNJ0ePSV++eUXTadAesLb2xs+Pj7o1KkT0tPTsWXLFuTk5FR6/4WmaF0x6NWrF44fP44PP/wQubm5cHZ2xpIlSxTPZiEi0iUvvPAC9u3bh40bN0IikeDZZ5/Fli1btK6nQyfOGRARUcPSunMGRETU+FgMiIhI+84ZaEJpaSlu376tGBeYiLSHEAIPHjyAk5NTjZ4jVVpaiqKiogbMTPtJpVK1PzMWAwC3b9/mlUpEWi41NVXtJ8wWFRUhJSVFIw+e1CZlT4GWSqXVxrIY4H+P2+6DF2CE6h/h+/7uICyfFKHWvBnb8LHvGn9TfSCAZR4Tsar9igbJQd3Yhpy3vsY+QjF+xL/VHktDCIE7d+7A0NAQrVq1arCn0mq7sh6PO3fuwNnZudpeDxYD/O+ZR0YwhpGk+mJgbmahVhxjGyfWUKrejWHmlprPtyHnrbex/3+9o7pduI8ePUJ+fj6cnJzqbdhXXdWsWTPcvn0bjx49qnasiqezZBKR3ip7Yqk6XSP6ruwzqG7UNoDFgIj0FC8GqdlnoDfFIDIyEm3atIGJiQl69uyJn3/+WdMpERHpDL04Z7Bnzx4EBwcjKioKPXv2REREBORyOZKSkmBvb6/p9IhIC2TcvIvsew0/hncZ66aWsHeu+2hpjUUvisHatWsxY8YMTJs2DQAQFRWFb7/9Flu3bsXChQs1nB0RaVrGzbsIcHsLxQWqo981FGMTY2y78mmNC0JkZCRWr16NtLQ0uLu7Y926dfDy8gIAFBQU4J133sFXX32FwsJCyOVyfPHFF0oDWNWWzncTFRUVISEhQWmcVQMDAwwePBinT5+u8D2FhYXIyclRehGR/sq+96BRCwEAFBcU1/hIpKyXIzQ0FOfPn4e7uzvkcrliuNy3334bhw8fxt69e3Hq1Cncvn27wjHRa0Pni8G9e/dQUlKiUhmbN2+OtLS0Ct8THh4Oa2trxYs3nBGRNijfy9G5c2dERUXBzMwMW7duRXZ2NrZs2YK1a9di4MCB8PT0RHR0NH766SecOXOmzsvW+WJQG4sWLUJ2drbilZqaqumUiOgpV10vR0JCAoqLi5Wmu7m5wdnZudJekJrQ+XMGTZs2haGhIdLT05Xa09PT4eDgUOF7ZDIZR7AiIq1SVS/HlStXkJaWBqlUChsbG5XplfWC1ITOHxlIpVJ4enoiLi5O0VZaWoq4uDh4e3trMDMiIt2h80cGABAcHAx/f3/06NEDXl5eiIiIQF5enuLqIiIibVddL4eDgwOKioqQlZWldHRQVS9ITej8kQEATJgwAWvWrMHixYvh4eGBxMRExMbG1svlVkREjaG6Xg5PT08YGxsrTU9KSsLNmzfrpRdEL44MAGD27NmYPXu2ptMgIqq1qno5rK2tMX36dAQHB8POzg5WVlaYM2cOvL298fzzz9d52XpTDIiIdN2ECRNw9+5dLF68GGlpafDw8FDq5fjkk09gYGCAl156Semms/rAYkBEes+6qSWMTYwb/Q5k66bqjcFQXlW9HCYmJoiMjERkZGRd01MhEUKIep+rjsnJyYG1tTWO7z4BczOLauM7erVH0s9X1Zo3YxmrjXnoUmxefi6GTBqA7OxsWFlZVRtfUFCAlJQUuLi4wMTERNH+ND6bqLLPoiI8Mihn+aQItQbbWHrwXSwevVKteTKWsdqYhy7FPhL182ve3rmZxnfO2kwvriYiIqK6YTEgIiIWAyIiYjEgIiKwGBAREVgMiIgILAZERATeZ0BET4k793OQlfuw0ZZnY2EKR7vqb5LTFiwGRKT37tzPwYuh21D0qKTRlik1MsSBsIAaF4TIyEisXr0aaWlpcHd3x7p16+Dl5QUA2LhxI3bt2oXz58/jwYMH+Oeff1QGu6ktdhMRkd7Lyn3YqIUAAIoeldT4SGTPnj0IDg5GaGgozp8/D3d3d8jlcmRkZAAA8vPzMWzYMLz33nv1ni+LARGRlli7di1mzJiBadOmoXPnzoiKioKZmRm2bt0KAAgKCsLChQvr5ZHVT2IxICLSAkVFRUhISFAa8N7AwACDBw+ulwHvq8NiQESkBe7du4eSkhKVERrra8D76rAYEBERiwERkTZo2rQpDA0NkZ6ertReXwPeV4fFgIhIC0ilUnh6eioNeF9aWoq4uLh6GfC+OhzpDBzpTBtjL6RnqhXbpY0Dkk9dbpAcONKZdsTWx0hnl2+mY0r4LrXyqk87F01GJ+fm1Qf+vz179sDf3x8bNmyAl5cXIiIi8PXXX+PKlSuKcwdpaWk4d+4cZsyYgfj4eFhaWsLZ2Rl2dnYq8+NIZ7XEkc60Jzbo3/FqxUa8MQrLtCBfjnTWcLH1NdKZLpgwYQLu3r2LxYsXIy0tDR4eHoiNjVWcVI6KikJYWJgivl+/fgCA6OhoBAQE1GnZLAZEpPdsLEwhNTJs9DuQbSxMa/y+2bNnY/bs2RVOW7JkCZYsWVLHzCrGYkBEes/RzgoHwgL4bKIq6HwxCA8Px/79+3HlyhWYmpqiV69eWLlyJTp27Kjp1IhIizjaWenUzrmx6fzVRKdOnUJgYCDOnDmD48ePo7i4GEOHDkVeXp6mUyMi0hk6f2QQGxur9Pe2bdtgb2+PhIQExckVIiKqms4fGTwpOzsbACq8zIqIiCqm80cG5ZWWliIoKAi9e/dG165dK40rLCxEYWGh4u+cnJzGSI+ISGvp1ZFBYGAgLl68iK+++qrKuPDwcFhbWyterVq1aqQMiYi0k94Ug9mzZ+PIkSM4ceIEWrZsWWXsokWLkJ2drXilpqY2UpZERNpJ57uJhBCYM2cODhw4gJMnT8LFxaXa98hkMshkskbIjohIN+h8MQgMDMSuXbtw6NAhWFpaKp77bW1tDVPTmt/9R0T6Ke1hFrKKGu+ScxupORxMbRpteXWl88Vg/fr1AAAfHx+l9vp4VgcR6Ye0h1kY+5+1KCp91GjLlBoYYV/f4BoVhPj4eKxevRoJCQm4c+cODhw4gNGjRyumCyEQGhqKTZs2ISsrC71798b69evh6upa53x1/pyBEKLCFwsBEZXJKspr1EIAAEWlj2p8JJKXlwd3d3dERkZWOH3VqlX47LPPEBUVhbNnz8Lc3BxyuRwFBQV1zlfnjwyIiPSFr68vfH19K5wmhEBERAQ++OADjBo1CgDw5Zdfonnz5jh48CAmTpxYp2Xr/JEBEdHTICUlBWlpaRg8eLCizdraGj179sTp06frPH8WAyIiHVB2cUzZ2AZlyga9qSuOdAaOdKaNsRzpjLFl6mOksyvZt/DK6Yr74RvSl96BcLNuUav3SiQSpRPIP/30E3r37o3bt2/D0dFRETd+/HhIJBLs2bNHZR4c6ayWONKZ9sRypDPGlnmaRjqrioODAwAgPT1dqRikp6fDw8OjzvNnNxERkQ5wcXGBg4MD4uLiFG05OTk4e/YsvL296zx/HhkQEWmJ3NxcXL36v66zlJQUJCYmws7ODs7OzggKCsKyZcvg6uoKFxcXhISEwMnJSelehNpiMSAi0hLnzp3DgAEDFH8HBwcDAPz9/bFt2zYsWLAAeXl5mDlzJrKystCnTx/ExsZWez5AHSwGRKT3bKTmkBoYNfodyDZS8xq9x8fHB1Vd0yORSLB06VIsXbq0rumpYDEgIr3nYGqDfX2D+WyiKrAYENFTwcHURqd2zo2NVxMRERGLARERsRgQERFYDIiICCwGREQEFgMiIgKLARERgfcZENFTQpTcBkr/abwFGthCYujUeMurIxYDItJ7ouQ2xN2hAIoacalSoNl3NSoI8fHxWL16NRISEnDnzh2l8QwAYP/+/YiKikJCQgLu37+PX3/9tV4eXw1wcBsAHNxGG2M5uA1jy9TH4Dai+A+IzBfVyqs+SZocgMS4i9rxR48exX//+194enpizJgxKsUgJiYGKSkpcHJywowZM6otBhzcppY4uI32xHJwG8aWeZoGt/H19YWvr2+l06dOnQoAuH79er0vmyeQiYiIxYCIiFgMiIgIelgMPvroI0gkEgQFBWk6FSIinaFXxeCXX37Bhg0b8Mwzz2g6FSIinaI3xSA3NxdTpkzBpk2bYGtrq+l0iIhqLDc3F4mJiUhMTAQApKSkIDExETdv3gQA3L9/H4mJibh06RIAICkpCYmJiUhLS6vzsvWmGAQGBmL48OEYPHhwtbGFhYXIyclRehERadq5c+fQvXt3dO/eHQAQHByM7t27Y/HixQCAb775Bt27d8fw4cMBABMnTkT37t0RFRVV52XrxX0GX331Fc6fP49ffvlFrfjw8HCEhYU1cFZEpDUMbAFI0eh3IBvUrJfCx8cHVd0HHBAQgICAgDrmVTGdLwapqal46623cPz48WrvsCuzaNEiBAcHK/7OyclBq1atGipFItIwiaET0Ow7PpuoCjpfDBISEpCRkYFnn31W0VZSUoL4+Hh8/vnnKCwshKGhodJ7ZDIZZDJZY6dKRBokMXQCdGjn3Nh0vhgMGjQIFy5cUGqbNm0a3Nzc8O6776oUAiIiUqXzxcDS0hJdu3ZVajM3N0eTJk1U2omIqGJ6czUREVF5fCBzzT4DnT8yqMjJkyc1nQIRaUhZ13BRURFMTU01nI1mFRU9vnpKne5yvSwGRPT0MjIygpmZGe7evQtjY2MYGDydHSClpaW4e/cuzMzMYGRU/a6exYCI9IpEIoGjoyNSUlJw48YNTaejUQYGBnB2doZEIqk2liOdgSOdaWMsRzpjbJmajnRWprS0VNFN8rSSSqVqHxnxyKAcjnSmPbEc6YyxZWo70pmBgYHaN6ISryYiIiKwGBAREVgMiIgILAZERAQWAyIiAosBERGBxYCIiMBiQEREYDEgIiKwGBAREVgMiIgILAZERAQWAyIiAosBERGBxYCIiMBiQERE4EhnADjSWWPF6uvoZRzprGFjazvSGdUMRzorhyOdcfSyho7Vljx0Kba2I51RzbCbiIiIWAyIiEhPisGtW7fw8ssvo0mTJjA1NUW3bt1w7tw5TadFRKQzdP6cwT///IPevXtjwIABOHr0KJo1a4bk5GTY2tpqOjUiIp2h88Vg5cqVaNWqFaKjoxVtLi4uGsyIiEj36Hw30TfffIMePXpg3LhxsLe3R/fu3bFp0yZNp0VEpFN0vhj89ddfWL9+PVxdXXHs2DG88cYbmDt3LrZv317pewoLC5GTk6P0IiJ6mul8N1FpaSl69OiBFStWAAC6d++OixcvIioqCv7+/hW+Jzw8HGFhYY2ZJhGRVtP5IwNHR0d07txZqa1Tp064efNmpe9ZtGgRsrOzFa/U1NSGTpOISKvp/JFB7969kZSUpNT2559/onXr1pW+RyaTQSaTNXRqREQ6Q+ePDN5++22cOXMGK1aswNWrV7Fr1y5s3LgRgYGBmk6NiEhn6HwxeO6553DgwAHs3r0bXbt2xYcffoiIiAhMmTJF06kREekMne8mAoARI0ZgxIgRmk6DiEhn6fyRARER1R2LARER6Uc3EekG49F31YqTOBY2cCZE9CSOdAaOdNZYsRcNbqsV28W6Ff48llR9YC1y0HSstuShS7Ec6axx8MigHI501rCx86T/Uit2zbNT8aEW5MuRzrQjliOdNQ6eMyAiIhYDIiJiMSAiIrAYEBERWAyIiAgsBkREBBYDIiICiwEREYHFgIiIwGJARERgMSAiIrAYEBERWAyIiAgsBkREBBYDIiICB7cBwMFtGiv2QnqmWrFd2jgg+dTlBslB07HakocuxXJwm8bBwW3K4eA2DRsb9O94tWIj3hiFZVqQLwe30Y5YDm7TONhNRERELAZERMRiQERE0INiUFJSgpCQELi4uMDU1BTt2rXDhx9+CJ4XJyJSn86fQF65ciXWr1+P7du3o0uXLjh37hymTZsGa2trzJ07V9PpERHpBJ0vBj/99BNGjRqF4cOHAwDatGmD3bt34+eff9ZwZkREukPnu4l69eqFuLg4/PnnnwCA3377DT/++CN8fX0rfU9hYSFycnKUXkRETzOdPzJYuHAhcnJy4ObmBkNDQ5SUlGD58uWYMmVKpe8JDw9HWFhYI2ZJRKTddP7I4Ouvv8bOnTuxa9cunD9/Htu3b8eaNWuwffv2St+zaNEiZGdnK16pqamNmDERkfbR+SOD+fPnY+HChZg4cSIAoFu3brhx4wbCw8Ph7+9f4XtkMhlkMlljpklEpNV0/sggPz8fBgbKq2FoaIjS0lINZUREpHt0/sjAz88Py5cvh7OzM7p06YJff/0Va9euxauvvqrp1IiIdIbOF4N169YhJCQEb775JjIyMuDk5IRZs2Zh8eLFmk6NiEhn6HwxsLS0REREBCIiIjSdChGRztL5cwZERFR3LAZERMSRzgCOdNZYsRzpTHvy0KVYjnTWOHT+nEF94khnHOmsoWO1JQ9diuVIZ42D3URERMRiQERELAZERAQWAyIiAosBERGBxYCIiMBiQEREYDEgIiKwGBAREVgMiIgILAZERAQWAyIiAosBERGBxYCIiMBiQEREYDEgIiJwpDMAHOmssWI50pn25KFLsRzprHFwpLNyONIZRzpr6FhtyUOXYjnSWeNgNxEREbEYEBGRDhSD+Ph4+Pn5wcnJCRKJBAcPHlSaLoTA4sWL4ejoCFNTUwwePBjJycmaSZaISEdpfTHIy8uDu7s7IiMjK5y+atUqfPbZZ4iKisLZs2dhbm4OuVyOgoKCRs6UiEh3af0JZF9fX/j6+lY4TQiBiIgIfPDBBxg1ahQA4Msvv0Tz5s1x8OBBTJw4sTFTJSLSWVp/ZFCVlJQUpKWlYfDgwYo2a2tr9OzZE6dPn9ZgZkREukXrjwyqkpaWBgBo3ry5Unvz5s0V0ypSWFiIwsJCxd85OTkNkyARkY7Q6SOD2goPD4e1tbXi1apVK02nRESkUTpdDBwcHAAA6enpSu3p6emKaRVZtGgRsrOzFa/U1NQGzZOISNvpdDFwcXGBg4MD4uLiFG05OTk4e/YsvL29K32fTCaDlZWV0ouI6Gmm9ecMcnNzcfXq/55fkpKSgsTERNjZ2cHZ2RlBQUFYtmwZXF1d4eLigpCQEDg5OWH06NGaS5qISMdofTE4d+4cBgwYoPg7ODgYAODv749t27ZhwYIFyMvLw8yZM5GVlYU+ffogNjYWJiYmmkqZiEjnaH0x8PHxQVUPVpVIJFi6dCmWLl3aiFkREekXnT5nQERE9YPFgIiItL+biLTbvZmVX7VVXlFrWzTdqN5d4dIX+tUlJSKqBY50Bo50VpdYjl7Gkc4aOpYjnTUOHhmUw5HOOHpZQ8dqSx66FMuRzhoHzxkQERGLARERsRgQERFYDIiICCwGREQEFgMiIgKLARERgcWAiIjAYkBERGAxICIisBgQERFYDIiICCwGREQEFgMiIgKLARERgcWAiIjAkc4AcKSzusRypDOOdNbQsRzprHFwpLNyONIZRzpr6FhtyUOXYjnSWeNgNxEREbEYEBGRDhSD+Ph4+Pn5wcnJCRKJBAcPHlRMKy4uxrvvvotu3brB3NwcTk5OeOWVV3D79m3NJUxEpIO0vhjk5eXB3d0dkZGRKtPy8/Nx/vx5hISE4Pz589i/fz+SkpIwcuRIDWRKRKS7tP4Esq+vL3x9fSucZm1tjePHjyu1ff755/Dy8sLNmzfh7OzcGCkSEek8rS8GNZWdnQ2JRAIbG5tKYwoLC1FYWKj4OycnpxEyIyLSXlrfTVQTBQUFePfddzFp0qQqr0cODw+HtbW14tWqVatGzJKISPvoTTEoLi7G+PHjIYTA+vXrq4xdtGgRsrOzFa/U1NRGypKISDvpRTdRWSG4ceMGfvjhh2rvUpTJZJDJZI2UHRGR9tP5YlBWCJKTk3HixAk0adJE0ykREekcrS8Gubm5uHr1f88vSUlJQWJiIuzs7ODo6IixY8fi/PnzOHLkCEpKSpCWlgYAsLOzg1Qq1VTaREQ6ReuLwblz5zBgwADF38HBwQAAf39/LFmyBN988w0AwMPDQ+l9J06cgI+PT2OlSUSk07S+GPj4+KCqB6vyoatERHWnN1cTERFR7Wn9kUFjKDu6eIRiQI0Djbz8XLUfq6vvsSVFBerF5j7Qinw1HasteehS7CM8jmEvQMPi4DYA/v77b954RqTlUlNT0bJlS02nobdYDACUlpbi9u3bsLS0hEQiAfD4ERWtWrVCamqq3o2upM/rBuj3+unzugEVr58QAg8ePICTkxMMDNiz3VDYTQTAwMCg0l8cVlZWevmlA/R73QD9Xj99XjdAdf2sra01mM3TgWWWiIhYDIiIiMWgUjKZDKGhoXr5DCN9XjdAv9dPn9cN0P/102Y8gUxERDwyICIiFgMiIgKLARERgcWAiIjAYlCpyMhItGnTBiYmJujZsyd+/vlnTadUZ0uWLIFEIlF6ubm5aTqtWouPj4efnx+cnJwgkUhw8OBBpelCCCxevBiOjo4wNTXF4MGDkZycrJlka6i6dQsICFDZlsOGDdNMsjUUHh6O5557DpaWlrC3t8fo0aORlJSkFFNQUIDAwEA0adIEFhYWeOmll5Cenq6hjJ8OLAYV2LNnD4KDgxEaGorz58/D3d0dcrkcGRkZmk6tzrp06YI7d+4oXj/++KOmU6q1vLw8uLu7IzIyssLpq1atwmeffYaoqCicPXsW5ubmkMvlKChQ7+F6mlTdugHAsGHDlLbl7t27GzHD2jt16hQCAwNx5swZHD9+HMXFxRg6dCjy8vIUMW+//TYOHz6MvXv34tSpU7h9+zbGjBmjwayfAoJUeHl5icDAQMXfJSUlwsnJSYSHh2swq7oLDQ0V7u7umk6jQQAQBw4cUPxdWloqHBwcxOrVqxVtWVlZQiaTid27d2sgw9p7ct2EEMLf31+MGjVKI/nUt4yMDAFAnDp1SgjxeDsZGxuLvXv3KmIuX74sAIjTp09rKk29xyODJxQVFSEhIQGDBw9WtBkYGGDw4ME4ffq0BjOrH8nJyXByckLbtm0xZcoU3Lx5U9MpNYiUlBSkpaUpbUdra2v07NlTL7YjAJw8eRL29vbo2LEj3njjDWRmZmo6pVrJzs4G8HioWgBISEhAcXGx0rZzc3ODs7Oz3mw7bcRi8IR79+6hpKQEzZs3V2pv3ry5YnxlXdWzZ09s27YNsbGxWL9+PVJSUtC3b188ePBA06nVu7JtpY/bEXjcRfTll18iLi4OK1euxKlTp+Dr64uSkhJNp1YjpaWlCAoKQu/evdG1a1cAj7edVCqFjY2NUqy+bDttxaeWPkV8fX0V/37mmWfQs2dPtG7dGl9//TWmT5+uwcyopiZOnKj4d7du3fDMM8+gXbt2OHnyJAYNGqTBzGomMDAQFy9e1OlzV/qCRwZPaNq0KQwNDVWuXEhPT4eDg4OGsmoYNjY26NChA65evarpVOpd2bZ6GrYjALRt2xZNmzbVqW05e/ZsHDlyBCdOnFB6hLyDgwOKioqQlZWlFK+v205bsBg8QSqVwtPTE3FxcYq20tJSxMXFwdvbW4OZ1b/c3Fxcu3YNjo6Omk6l3rm4uMDBwUFpO+bk5ODs2bN6tx2Bx6P1ZWZm6sS2FEJg9uzZOHDgAH744Qe4uLgoTff09ISxsbHStktKSsLNmzf1cttpC3YTVSA4OBj+/v7o0aMHvLy8EBERgby8PEybNk3TqdXJvHnz4Ofnh9atW+P27dsIDQ2FoaEhJk2apOnUaiU3N1fpl3BKSgoSExNhZ2cHZ2dnBAUFYdmyZXB1dYWLiwtCQkLg5OSE0aNHay5pNVW1bnZ2dggLC8NLL70EBwcHXLt2DQsWLED79u0hl8s1mLV6AgMDsWvXLhw6dAiWlpaK8wDW1tYwNTWFtbU1pk+fjuDgYNjZ2cHKygpz5syBt7c3nn/+eQ1nr8c0fTmTtlq3bp1wdnYWUqlUeHl5iTNnzmg6pTqbMGGCcHR0FFKpVLRo0UJMmDBBXL16VdNp1dqJEycEAJWXv7+/EOLx5aUhISGiefPmQiaTiUGDBomkpCTNJq2mqtYtPz9fDB06VDRr1kwYGxuL1q1bixkzZoi0tDRNp62WitYLgIiOjlbEPHz4ULz55pvC1tZWmJmZiRdffFHcuXNHc0k/BfgIayIi4jkDIiJiMSAiIrAYEBERWAyIiAgsBkREBBYDIiICiwEREYHFgIiIwGJARERgMSAiIrAYEBERWAyIiAjA/wG/qSQUYjrmIgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 300x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(21,)\n",
      "Reshaping: Overflow in cluster 0 at 2000ps:\n",
      "  Top overflow: [   0.    0.    0.    0.    0.    0.    0.    0.    0.    9. 2278.    0.\n",
      "    0.    0.    0.    0.    0.    0.    0.    0.    0.]\n",
      "  Bottom overflow: [   0.    0.    0.    0.    0.    0.    0. 3114. 1416.    0.    0.    0.\n",
      "    0.    0.    0.    0.    0.    0.    0.    0.    0.]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV4AAAEpCAYAAADbDEpjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAmF0lEQVR4nO3de1xUZR4/8M+AMMhtQAWBnyDkDS9JikKsqXnPzEuWbWYFqbka1KLVGtvmLRVNU7qgpCl2U7xsWFlqaYq1eUXd1BR180IioikzOMqAM8/vD3/Mr4nhcuZyBmY+79frvF475znneb7z7O7Hw5lzUQghBIiISDZuji6AiMjVMHiJiGTG4CUikhmDl4hIZgxeIiKZMXiJiGTG4CUikhmDl4hIZgxeIiKZMXht7Pz581AoFFizZo1dx4mMjERSUpJN+5w1axYUCoVN+2wskpKSEBkZ6egyyEUweCVas2YNFAqF2eW1115zdHlmlZeXY+nSpYiPj4dKpYKXlxfat2+PlJQUnD59WrY61q5di4yMDNnGA4AHH3zQ5L+jZs2aoWfPnli9ejUMBoNNxpg/fz42b95sk77INTRxdAGN1Zw5cxAVFWWyrkuXLmjdujVu374NDw8PB1Vm6tq1a3jooYeQn5+PRx55BE899RR8fX1RUFCAnJwcrFixAhUVFbLUsnbtWhw/fhypqamyjFelVatWSE9PBwBcvXoVH3/8MSZMmIDTp09jwYIFVvc/f/58PP744xg1apTVfZFrYPBaaOjQoejRo4fZNi8vL5mrqVlSUhKOHDmCTZs24bHHHjNpe/PNN/H66687qDLbMBgMqKioqHXOVSoVnn76aePnv/3tb+jQoQPef/99vPnmmw3mH0lyHTzVYGPmzvEmJSXB19cXly5dwqhRo+Dr64ugoCC88sor0Ov1JvsvXrwYf/nLX9C8eXM0bdoUsbGx2LRpk0W17N+/H19//TUmTJhQLXQBQKlUYvHixZK+SxWFQoFZs2YZP5eVlSE1NRWRkZFQKpUIDg7GoEGDcPjwYQB3/+T/+uuvceHCBeOf/X88p6rT6TBz5ky0bdsWSqUS4eHh+Mc//gGdTldt3JSUFHz22Wfo3LkzlEoltm3bJmlevL29cf/990Or1eLq1as1bqfVavHyyy8jPDwcSqUSHTp0wOLFi/HHB/opFApotVp89NFHxu9l63Pv5Hx4xGshtVqNa9eumaxr0aJFjdvr9XoMGTIE8fHxWLx4MXbs2IG3334bbdq0wZQpU4zbvfPOOxgxYgTGjRuHiooK5OTkYMyYMdiyZQuGDRsmqcYvv/wSAPDMM89I2s8SkydPxqZNm5CSkoJOnTrh999/x48//oiTJ0+ie/fueP3116FWq/Hbb79h6dKlAABfX18Ad49aR4wYgR9//BGTJk1Cx44dcezYMSxduhSnT5+udv70+++/x4YNG5CSkoIWLVpY9KPYr7/+Cnd3dwQEBJhtF0JgxIgR2LVrFyZMmID77rsP27dvx6uvvopLly4Zv8Mnn3yCiRMnIi4uDpMmTQIAtGnTRnI95GIESZKdnS0AmF2EEOLcuXMCgMjOzjbuk5iYKACIOXPmmPTVrVs3ERsba7Lu1q1bJp8rKipEly5dRP/+/U3Wt27dWiQmJtZa66OPPioAiBs3btTru82cOVP88X8S5r5LFQBi5syZxs8qlUokJyfX2v+wYcNE69atq63/5JNPhJubm/jhhx9M1mdlZQkA4j//+Y/JuG5ubuLEiRP1+k59+/YV0dHR4urVq+Lq1avi5MmT4qWXXhIAxPDhw43bJSYmmtS2efNmAUDMnTvXpL/HH39cKBQKcfbsWeM6Hx+fOv+7IPojHvFaKDMzE+3bt5e0z+TJk00+9+7dG5988onJuqZNmxr/840bN6DX69G7d2+sW7dOco0ajQYA4OfnJ3lfqQICArB//34UFRUhLCxM0r4bN25Ex44dER0dbfJXRP/+/QEAu3btwl/+8hfj+r59+6JTp0717v/UqVMICgoyflYoFBg2bBhWr15d4z7ffPMN3N3d8dJLL5msf/nll7Fp0yZs3boVKSkp9a6B6I8YvBaKi4ur8cc1c7y8vEz+zw8AgYGBuHHjhsm6LVu2YO7cuTh69KjJ+U1Lrq/19/cHcPf8a01/UtvKW2+9hcTERISHhyM2NhYPP/wwnn32Wdxzzz117nvmzBmcPHmy2vxUKSkpMfn856tJ6hIZGYmVK1dCoVDAy8sL7dq1Q3BwcK37XLhwAWFhYdX+0erYsaOxnchSDF6ZuLu717nNDz/8gBEjRqBPnz5YtmwZQkND4eHhgezsbKxdu1bymNHR0QCAY8eOoXfv3pL3ryns//yDIAA88cQT6N27N3Jzc/Htt99i0aJFWLhwIT7//HMMHTq01nEMBgPuvfdeLFmyxGx7eHi4yec//lVQHz4+Phg4cKCkfYjsicHbgPz73/+Gl5cXtm/fDqVSaVyfnZ1tUX/Dhw9Heno6Pv30U4uCNzAwEABQWlpqsr6mo73Q0FC88MILeOGFF1BSUoLu3btj3rx5xuCtKcjbtGmD//73vxgwYECDuXOudevW2LFjB8rKykyOek+dOmVsr9JQaqbGg5eTNSDu7u5QKBQmR5Tnz5+3+K6ohIQEPPTQQ/jwww/N9lFRUYFXXnmlxv39/f3RokUL7Nmzx2T9smXLTD7r9Xqo1WqTdcHBwQgLCzM5XeLj41NtO+Du0fKlS5ewcuXKam23b9+GVqutsUZ7efjhh6HX6/H++++brF+6dCkUCoXJUbyPj0+1f5yIasMj3gZk2LBhWLJkCR566CE89dRTKCkpQWZmJtq2bYuff/7Zoj4//vhjDB48GKNHj8bw4cMxYMAA+Pj44MyZM8jJycHly5drvZZ34sSJWLBgASZOnIgePXpgz5491W4zLisrQ6tWrfD4448jJiYGvr6+2LFjBw4ePIi3337buF1sbCzWr1+PadOmoWfPnvD19cXw4cPxzDPPYMOGDZg8eTJ27dqFXr16Qa/X49SpU9iwYQO2b98u6Xy6LQwfPhz9+vXD66+/jvPnzyMmJgbffvstvvjiC6SmpppcMhYbG4sdO3ZgyZIlCAsLQ1RUFOLj42WtlxoZR19W0dhUXU528OBBs+01XU7m4+NTbds/X74lhBCrVq0S7dq1E0qlUkRHR4vs7Gyz29XncrIqt27dEosXLxY9e/YUvr6+wtPTU7Rr1068+OKLJpdFmRvn1q1bYsKECUKlUgk/Pz/xxBNPiJKSEpPLyXQ6nXj11VdFTEyM8PPzEz4+PiImJkYsW7bMpK+bN2+Kp556SgQEBAgAJpdvVVRUiIULF4rOnTsLpVIpAgMDRWxsrJg9e7ZQq9XG7QDUednaH/Xt21d07ty5zu3+fDmZEEKUlZWJqVOnirCwMOHh4SHatWsnFi1aJAwGg8l2p06dEn369BFNmzYVAHhpGdVJIcQfbsMhIiK74zleIiKZMXiJiGTG4CUikhmDl4hIZgxeIiKZMXiJiGQm+w0UBoMBRUVF8PPz462WRA2MEAJlZWUICwuDm5u04zK9Xo/Kyko7VdbweXp61nvOZA/eoqKiag89IaKGpbCwEK1atarXtkIIFBcXu/xt025uboiKioKnp2ed28oevFUPHHkAD6MJqr/r6vV1qZg3NsPsvpa22atfe4453ePLGvtdWDmiUX0XVx/TXv3aY8w7qMSP+EbSM5yrQjc4OBje3t4u+Zds1V/yly9fRkRERJ1zIHvwVhXUBB5ooqgevD7evmbXW9Nmr37tOaa7p9JsGwD4VDSu7+LqY9qrX7uM+f/uY61veOr1emPoNm/evF77OKugoCAUFRXhzp07db5AlT+uEZHFqs7pent7O7gSx6s6xWDuedV/xuAlIqu54umFP5MyBwxeIiKZ8Xm8RGRzJRevQn2tTLbxVC38EBxh/p19DRGDl4hsquTiVSRF/x2V5fJd0+vh5YE1p96RHL6ZmZlYtGgRiouLERMTg/feew9xcXEAgPLycrz88svIycmBTqfDkCFDsGzZMrRs2dLqenmqgYhsSn2tTNbQBYDK8krJR9hVb0OZOXMmDh8+jJiYGAwZMsT4VuupU6fiq6++wsaNG5GXl4eioiKMHj3aJvUyeInIJS1ZsgTPP/88nnvuOXTq1AlZWVnw9vbG6tWroVarsWrVKixZsgT9+/dHbGwssrOz8dNPP2Hfvn1Wj21R8GZmZiIyMhJeXl6Ij4/HgQMHrC6EiEguFRUVyM/Px8CBA43r3NzcMHDgQOzduxf5+fmorKw0aY+OjkZERAT27t1r9fiSg7euw3Mioobu2rVr0Ov11c7XtmzZEsXFxSguLoanpycCAgLMtltLcvDWdnhORER1kxS8dR2eExE1Bi1atIC7uzuuXLlisv7KlSsICQlBSEgIKioqqj34p6rdWpKCt67Dc3N0Oh00Go3JQkTkSJ6enoiNjcXOnTuN6wwGA3bu3ImEhATExsbCw8PDpL2goAAXL15EQkKC1ePb/Tre9PR0zJ49297DEBFJMm3aNCQmJqJHjx6Ii4tDRkYGtFotnnvuOahUKkyYMAHTpk1Ds2bN4O/vjxdffBEJCQm4//77rR5bUvDWdXhuTlpaGqZNm2b8rNFo+DxeInK4v/71r7h69SpmzJiB4uJi3Hfffdi2bZvxL/qlS5fCzc0Njz32mMkNFLYg6VRDXYfn5iiVSvj7+5ssROS8VC384OFV+2MRbc3DywOqFvV/hnCVlJQUXLhwATqdDvv370d8fLyxzcvLC5mZmbh+/Tq0Wi0+//xzm5zfBSw41VDb4TkRUXBEENaceofPaqiF5OCt6/CciCg4IqhRBaHcLPpxLSUlBSkpKbauhYjIJfBZDUREMmPwEhHJjMFLRCQzBi8RkcwYvEREMmPwEhHJjO9cIyKbu3xdg9Kbt2UbL8C3KUKbNZ67Yhm8RGRTl69r8OjMNai4o5dtTM8m7sidnSQ5fGt72eWKFSuwdu1aHD58GGVlZbhx40a1B6NbSiGEEDbpqZ40Gg1UKhW+W7cLPt6+1do7xLVFwYGzZve1tM1e/Vo75rErv5tt6xwZgjN5Jy3u11Xmr7GMaa9+7TGm9tZNDBrbD2q1ul7PVSkvL8e5c+cQFRUFLy8vAMDJi1cwLn1tnfva2mdpT6FjRP3voF2/fj2effZZZGVlIT4+HhkZGdi4cSMKCgoQHByMjIwMlJeXA7j7sK+6gtfcXNTEYUe888ZmoImi+oM05myejhmjFprdx9I2e/Vr7Zip3+wx25YxZSTmNrLvwjGd47vcEfK+HdiR/vg2HQDIysrC119/jdWrV+O1115DamoqAGD37t02H5s/rhGRy3H023QYvETkcix5m44tMXiJiGTG4CUil2PJ23RsicFLRC7Hkrfp2BKv4yUil1TX23SKi4tRXFyMs2fvXnZ37Ngx+Pn5ISIiAs2aNbNqbAYvEbmkut6mk5WVZfKG9D59+gAAsrOzkZSUZNXYDF4isqkA36bwbOIu+51rAb5NJe9X29t0Zs2ahVmzZllZmXmSg3fPnj1YtGgR8vPzcfnyZeTm5mLUqFF2KI2IGqPQZv7InZ3EZzXUQnLwarVaxMTEYPz48Rg9erQ9aiKiRi60mX+jCkK5SQ7eoUOHYujQofaohYjIJfByMiIimdn9xzWdTgedTmf8rNFo7D0kEVGDZvcj3vT0dKhUKuMSHh5u7yGJiBo0uwdvWloa1Gq1cSksLLT3kEREDZrdTzUolUoolUp7D0NE1GhIDt6bN28ab6EDgHPnzuHo0aNo1qwZIiIibFocEZEzkhy8hw4dQr9+/Yyfp02bBgBITEzEmjVrbFYYETVexbdLUVqhlW28AE8fhDQNkG08a0kO3gcffBAyv6aNiBqR4tulePyHJagw3JFtTE+3JtjUe5qk8K3rLlwhBGbOnImVK1eitLQUvXr1wvLly9GuXTur6+V1vERkU6UVWllDFwAqDHckH2FX3YWbmZlptv2tt97Cu+++i6ysLOzfvx8+Pj4YMmSI8QWY1uBDcojIJdV2F64QAhkZGfjXv/6FkSNHAgA+/vhjtGzZEps3b8aTTz5p1dg84iUi+pNz586huLjY5GWYKpUK8fHxNnkZJoOXiOhPql54aa+XYTJ4iYhkxuAlIvqTqhde2utlmAxeIqI/iYqKQkhIiMnLMDUaDfbv32+Tl2HyqgYickl13YWbmpqKuXPnol27doiKisIbb7yBsLAwm7xxRyFkvhtCo9FApVLhu3W74OPtW629Q1xbFBw4a2ZPy9vs1a+1Yx678rvZts6RITiTd9Lifl1l/hrLmPbq1x5jam/dxKCx/aBWq+HvX/cbJMrLy3Hu3DlERUXBy8sLAHBKfQnP7jV/baw9fZyQjGjV/6n39rt37za5C7dK1V24VTdQrFixAqWlpXjggQewbNkytG/f3mx/5uaiJg474p03NgNNFB7V1s/ZPB0zRi00u4+lbfbq19oxU7/ZY7YtY8pIzG1k34VjOsd3uSMqa+yvvgI8feDp1kT2O9cCPH0k7VPXXbgKhQJz5szBnDlzrC2vGp5qICKbCmkagE29p/FZDbVg8BKRzYU0DWhUQSg3XtVARCQzBi8RkcwYvEREMmPwEhHJjMFLRCQzBi8RkcwkBW96ejp69uwJPz8/BAcHY9SoUSgoKLBXbURETknSdbx5eXlITk5Gz549cefOHfzzn//E4MGD8csvv8DHR9pdI0TkvIS+CDDckG9At0Ao3MPkG89KkoJ327ZtJp/XrFmD4OBg5Ofno0+fPjYtjIgaJ6Evgrg6GECFjKN6AkHfSgrful52+fnnnyMrKwv5+fm4fv06jhw5gvvuu88m1Vp1jletVgMAmjVrVuM2Op0OGo3GZCEiJ2a4AXlDF3fHk3iEXdfLLrVaLR544AEsXFjz8zAsZfEtwwaDAampqejVqxe6dOlS43bp6emYPXu2pcMQEdlFbS+7BIBnnnkGAHD+/Hmbj23xEW9ycjKOHz+OnJycWrdLS0uDWq02LoWFhZYOSUTkFCw64k1JScGWLVuwZ88etGrVqtZtlUollEqlRcURETkjScErhMCLL76I3Nxc7N69G1FRUfaqi4jIaUkK3uTkZKxduxZffPEF/Pz8jK85VqlUaNq0qV0KJCJyNpLO8S5fvhxqtRoPPvggQkNDjcv69evtVR8RkdORfKqBiMgZ1PWyy+vXr+PixYsoKioCAONduiEhIVa/4p3PaiAil3To0CF069YN3bp1AwBMmzYN3bp1w4wZMwAAX375Jbp164Zhw4YBAJ588kl069YNWVlZVo/NV/8QkW25BQLwhOx3rrkFStqjrpddJiUlISkpycq6zGPwEpFNKdzDgKBv+ayGWjB4icjmFO5hQCMKQrnxHC8RkcwYvEREMmPwEpHVeKmptDlg8BKRxTw8PAAAt27dcnAljldRcfcqDnd39zq3VQiZ/6nSaDRQqVT4bt0u+Hj7VmvvENcWBQfOmtnT8jZ79WvtmMeu/G62rXNkCM7knbS4X1eZv8Yypr36tceY2ls3MWhsP6jVavj7+9fY9x9dvnwZpaWlCA4Ohre3NxQKRb32cyYGgwFFRUXw8PBAREREnXPgsKsa5o3NQBOFR7X1czZPx4xR5h88bGmbvfq1dszUb/aYbcuYMhJzG9l34ZjO8V3uiMoa+6tJ1V1cJSUlkvd1Jm5ubvUKXYCXkxGRlRQKBUJDQxEcHIzKSunB7Sw8PT3h5la/s7cMXiKyCXd393qd3yT+uEZEJDsGLxGRzBi8REQyY/ASEcmMwUtEJDMGLxGRzCS/c61r167w9/eHv78/EhISsHXrVnvVRkTklCQFb6tWrbBgwQLk5+fj0KFD6N+/P0aOHIkTJ07Yqz4iIqcj6QaK4cOHm3yeN28eli9fjn379qFz5842LYyIyFlZfOeaXq/Hxo0bodVqkZCQUON2Op0OOp3O+Fmj0Vg6JBGRU5D849qxY8fg6+sLpVKJyZMnIzc3F506dapx+/T0dKhUKuMSHh5uVcFERI2d5ODt0KEDjh49iv3792PKlClITEzEL7/8UuP2aWlpUKvVxqWwsNCqgomIGjvJpxo8PT3Rtm1bAEBsbCwOHjyId955Bx988IHZ7ZVKJZRKpXVVEhE5Eauv4zUYDCbncImIqHaSjnjT0tIwdOhQREREoKysDGvXrsXu3buxfft2e9VHROR0JAVvSUkJnn32WVy+fBkqlQpdu3bF9u3bMWjQIHvVR0TkdCQF76pVq+xVBxGRy+CzGoiIZMbgJSKSGYOXiEhmDF4iIpkxeImIZMbgJSKSGYOXiEhmDF4iIpkphBBCzgE1Gg1UKhW+W7cLPt6+1do7xLVFwYGzZve1tM1e/Vo75rErv5tt6xwZgjN5Jy3u11Xmr7GMaa9+7TGm9tZNDBrbD2q1Gv7+/jX2Tdax+EHo1po3NgNNFB7V1s/ZPB0zRi00u4+lbfbq19oxU7/ZY7YtY8pIzG1k34VjOsd3uSMqa+yPbIenGoiIZMbgJSKSGYOXiEhmDF4iIpkxeImIZMbgJSKSGYOXiEhmDF4iIplZFbwLFiyAQqFAamqqjcohInJ+FgfvwYMH8cEHH6Br1662rIeIyOlZFLw3b97EuHHjsHLlSgQGBtq6JiIip2ZR8CYnJ2PYsGEYOHBgndvqdDpoNBqThYjIlUl+SE5OTg4OHz6MgwcP1mv79PR0zJ49W3JhRETOStIRb2FhIf7+97/js88+g5eXV732SUtLg1qtNi6FhYUWFUpE5CwkHfHm5+ejpKQE3bt3N67T6/XYs2cP3n//feh0Ori7u5vso1QqoVQqbVMtEZETkBS8AwYMwLFjx0zWPffcc4iOjsb06dOrhS4REVUnKXj9/PzQpUsXk3U+Pj5o3rx5tfVERGQe71wjIpKZ1a/+2b17tw3KICJyHTziJSKSGYOXiEhmDF4iIpkxeImIZMbgJSKSGYOXiEhmDF4iIplZfR0v1ayidSCuTUqosa3Fir1m2zwf7mPPsojIwRRCCCHngBqNBiqVCt+t2wUfb99q7R3i2qLgwFmz+1raZq9+6xqzXd+OOHG+2Gxb58gQnMk7aVG/jvguHLPh9WuPMbW3bmLQ2H5Qq9Xw9/evsW+yjsOOeOeNzUAThUe19XM2T8eMUQvN7mNpm736rWvMfx1ZgNTlX5hty5gyEnMb0XfhmA2vX3uMeUdU1tgf2Q7P8RIRyYzBS0QkMwYvEZHMGLxERDJj8BIRyYzBS0QkMwYvEZHMGLxERDKTFLyzZs2CQqEwWaKjo+1VGxGRU5J851rnzp2xY8eO/99BEz7ugYhICsmp2aRJE4SEhNijFiIilyD5HO+ZM2cQFhaGe+65B+PGjcPFixdr3V6n00Gj0ZgsRESuTFLwxsfHY82aNdi2bRuWL1+Oc+fOoXfv3igrK6txn/T0dKhUKuMSHh5uddFERI2ZpOAdOnQoxowZg65du2LIkCH45ptvUFpaig0bNtS4T1paGtRqtXEpLCy0umgiosbMql/GAgIC0L59e5w9W/MzQZVKJZRKpTXDEBE5Fauu47158yb+97//ITQ01Fb1EBE5PUnB+8orryAvLw/nz5/HTz/9hEcffRTu7u4YO3asveojInI6kk41/Pbbbxg7dix+//13BAUF4YEHHsC+ffsQFBRkr/qIiJyOpODNycmxVx1ERC6Dz2ogIpIZg5eISGYMXiIimTF4iYhkxuAlIpIZg5eISGYMXiIimfEp5nakCNXBY9TVGtuIyDUphBBCzgE1Gg1UKhW+W7cLPt6+1do7xLVFwQHzD92xtM1e/dY1ZvshHXBCbf5pbJ1V4Ti9vcCifh3xXThmw+vXHmNqb93EoLH9oFar4e/vX2PfZB2HHfHOG5uBJgqPauvnbJ6OGaMWmt3H0jZ79VvXmG8Uz8Yrhz8x27a4+zN4sxF9F47Z8Pq1x5h3RGWN/ZHt8BwvEZHMGLxERDJj8BIRyYzBS0QkMwYvEZHMGLxERDJj8BIRyYzBS0QkM8nBe+nSJTz99NNo3rw5mjZtinvvvReHDh2yR21ERE5J0p1rN27cQK9evdCvXz9s3boVQUFBOHPmDAIDA+1VHxGR05EUvAsXLkR4eDiys7ON66KiomxeFBGRM5N0quHLL79Ejx49MGbMGAQHB6Nbt25YuXJlrfvodDpoNBqThYjIlUkK3l9//RXLly9Hu3btsH37dkyZMgUvvfQSPvrooxr3SU9Ph0qlMi7h4eFWF01E1JhJCl6DwYDu3btj/vz56NatGyZNmoTnn38eWVlZNe6TlpYGtVptXAoLzT8mkYjIVUgK3tDQUHTq1MlkXceOHXHx4sUa91EqlfD39zdZiIhcmaTg7dWrFwoKTB/effr0abRu3dqmRREROTNJwTt16lTs27cP8+fPx9mzZ7F27VqsWLECycnJ9qqPiMjpSArenj17Ijc3F+vWrUOXLl3w5ptvIiMjA+PGjbNXfURETkfyq38eeeQRPPLII/aohYjIJfBZDUREMmPwEhHJjMFLRCQzBi8RkcwYvEREMmPwEhHJjMFLRCQzhRBCyDmgRqOBSqXCd+t2wcfbt1p7h7i2KDhw1uy+lrbZq98OcW0RGJBf45jXFXE4oTb/UKDOqnCc3l5gts1R34VjNq5+7TGm9tZNDBrbD2q1ms9VsSPJN1DYyryxGWii8Ki2fs7m6ZgxaqHZfSxts1e/czZPR3zc5BrHPOH2BV45/InZtsXdn8GbDey7cMzG1a89xrwjKmvsj2yHpxqIiGTG4CUikhmDl4hIZgxeIiKZMXiJiGTG4CUikhmDl4hIZgxeIiKZSQreyMhIKBSKagvfuUZEVH+S7lw7ePAg9Hq98fPx48cxaNAgjBkzxuaFERE5K0nBGxQUZPJ5wYIFaNOmDfr27WvTooiInJnF53grKirw6aefYvz48VAoFLasiYjIqVn8kJzNmzejtLQUSUlJtW6n0+mg0+mMnzUajaVDEhE5BYuPeFetWoWhQ4ciLCys1u3S09OhUqmMS3h4uKVDEhE5BYuC98KFC9ixYwcmTpxY57ZpaWlQq9XGpbDQ/PNpiYhchUWnGrKzsxEcHIxhw4bVua1SqYRSqbRkGCIipyT5iNdgMCA7OxuJiYlo0sRhz1EnImq0JAfvjh07cPHiRYwfP94e9RAROT3Jh6yDBw+GzK9pIyJyKnxWAxGRzBi8REQyY/ASEcmMwUtEJDMGLxGRzBi8REQyY/ASEclM9lvPqq4BvoNKwMzlwNpbN3FHVJrd19I2e/WrvXUTmjJDzWO63YReqzPfVtbwvgvHbFz92mPMO7i7jtfq25dCyDzDv/32G59QRtTAFRYWolWrVo4uw2nJHrwGgwFFRUXw8/NDWVkZwsPDUVhYCH9/fznLaDQ0Gg3nqA6co9pJmR8hBMrKyhAWFgY3N56JtBfZTzW4ubkZ/yWtenOFv78//w9TB85R3ThHtavv/KhUKhmqcW38J42ISGYMXiIimTk0eJVKJWbOnMkHpdeCc1Q3zlHtOD8Nj+w/rhERuTqeaiAikhmDl4hIZgxeIiKZMXiJiGTm0ODNzMxEZGQkvLy8EB8fjwMHDjiyHIfas2cPhg8fjrCwMCgUCmzevNmkXQiBGTNmIDQ0FE2bNsXAgQNx5swZxxTrAOnp6ejZsyf8/PwQHByMUaNGoaCgwGSb8vJyJCcno3nz5vD19cVjjz2GK1euOKhi+S1fvhxdu3Y13iiRkJCArVu3GttdfX4aEocF7/r16zFt2jTMnDkThw8fRkxMDIYMGYKSkhJHleRQWq0WMTExyMzMNNv+1ltv4d1330VWVhb2798PHx8fDBkyBOXl5TJX6hh5eXlITk7Gvn378N1336GyshKDBw+GVqs1bjN16lR89dVX2LhxI/Ly8lBUVITRo0c7sGp5tWrVCgsWLEB+fj4OHTqE/v37Y+TIkThx4gQAzk+DIhwkLi5OJCcnGz/r9XoRFhYm0tPTHVVSgwFA5ObmGj8bDAYREhIiFi1aZFxXWloqlEqlWLdunQMqdLySkhIBQOTl5Qkh7s6Hh4eH2Lhxo3GbkydPCgBi7969jirT4QIDA8WHH37I+WlgHHLEW1FRgfz8fAwcONC4zs3NDQMHDsTevXsdUVKDdu7cORQXF5vMl0qlQnx8vMvOl1qtBgA0a9YMAJCfn4/KykqTOYqOjkZERIRLzpFer0dOTg60Wi0SEhI4Pw2M7A/JAYBr165Br9ejZcuWJutbtmyJU6dOOaKkBq24uBgAzM5XVZsrMRgMSE1NRa9evdClSxcAd+fI09MTAQEBJtu62hwdO3YMCQkJKC8vh6+vL3Jzc9GpUyccPXqU89OAOCR4iayRnJyM48eP48cff3R0KQ1Ohw4dcPToUajVamzatAmJiYnIy8tzdFn0Jw451dCiRQu4u7tX+0X1ypUrCAkJcURJDVrVnHC+gJSUFGzZsgW7du0yeVB3SEgIKioqUFpaarK9q82Rp6cn2rZti9jYWKSnpyMmJgbvvPMO56eBcUjwenp6IjY2Fjt37jSuMxgM2LlzJxISEhxRUoMWFRWFkJAQk/nSaDTYv3+/y8yXEAIpKSnIzc3F999/j6ioKJP22NhYeHh4mMxRQUEBLl686DJzZI7BYIBOp+P8NDSO+lUvJydHKJVKsWbNGvHLL7+ISZMmiYCAAFFcXOyokhyqrKxMHDlyRBw5ckQAEEuWLBFHjhwRFy5cEEIIsWDBAhEQECC++OIL8fPPP4uRI0eKqKgocfv2bQdXLo8pU6YIlUoldu/eLS5fvmxcbt26Zdxm8uTJIiIiQnz//ffi0KFDIiEhQSQkJDiwanm99tprIi8vT5w7d078/PPP4rXXXhMKhUJ8++23QgjOT0PisOAVQoj33ntPRERECE9PTxEXFyf27dvnyHIcateuXQJ3X/9psiQmJgoh7l5S9sYbb4iWLVsKpVIpBgwYIAoKChxbtIzMzQ0AkZ2dbdzm9u3b4oUXXhCBgYHC29tbPProo+Ly5cuOK1pm48ePF61btxaenp4iKChIDBgwwBi6QnB+GhI+FpKISGZ8VgMRkcwYvEREMmPwEhHJjMFLRCQzBi8RkcwYvEREMmPwEhHJjMFLRCQzBi8RkcwYvEREMmPwEhHJjMFLRCSz/wso8yRBK5PBswAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 300x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhAAAAHDCAYAAAB1QiOVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABBPElEQVR4nO3de1hVddr/8Q/IySPgCWQU0zRP5WEwPJsHlMyxPFTq2KTW1OQjjkpPYzSNp6bQLKMMdSrCetKwk1ramEalUx7DLMs0LUtNgTRFxASE9fujn3vagcpCNmt/4/26rn1dsfbae39uZIabtdb9XT6WZVkCAACwwdfpAAAAwDw0EAAAwDYaCAAAYBsNBAAAsI0GAgAA2EYDAQAAbKOBAAAAttFAAAAA22ggAACAbTQQqBJuuOEG3XXXXU7H0BVXXKE//OEP5Xrtt99+Kx8fHz322GMVnKp8+vTpoz59+jgdQ0uWLJGPj48+/vjjS+5bkZkLCwvVpEkTLVy4sELeDzANDUQVtXDhQvn4+KhLly5OR7mo9PR03XHHHbrqqqtUo0YNNW/eXH/+85919OjRMr/HRx99pHXr1mnatGkeTIryOnTokGbNmqXo6GiFhoaqfv366tOnj959912Pf/aRI0c0c+ZM7dy50/Zr/f39FR8fr4cfflhnz56t+HCAl6OBqKKWLl2qK664Qtu2bdP+/fudjnNB06ZN0wcffKBhw4bpqaee0qhRo/TKK6+oU6dOyszMLNN7zJs3T/3791eLFi08nLZqWbdundatW3fZ77Nq1SrNnTtXLVq00D//+U/94x//UG5urgYMGKDU1NQKSPpfv8585MgRzZo1q1wNhCSNHz9ex44d07JlyyooIWAOGogq6MCBA9q0aZPmz5+vBg0aaOnSpWV63blz51RQUFDqc3l5eRUZ0WX+/Pnav3+/5s6dqz//+c965JFHtHr1amVlZenpp5++5Ouzs7O1Zs0a3XrrrR7JV5UFBAQoICDgst+nb9++OnjwoJYtW6aJEydq8uTJ2rRpk1q3bq3p06dXQNL/qqjM54WEhGjgwIFasmRJhb0nYAoaiCpo6dKlCg0N1eDBg3XzzTeX2kD88nx7UlKSrrzySgUGBmr37t2aOXOmfHx8tHv3bv3xj39UaGioevbsKUn67LPPNG7cODVv3lxBQUEKDw/XHXfcoePHj7ve+/3335ePj49WrFhR4nOXLVsmHx8fbd68WZLUu3dv+fq6/5j27t1bdevW1ZdffnnJWtesWaNz584pJiamxHMnT57UlClT1KRJEwUGBqpFixaaO3euiouL3fZ77LHH1L17d9WrV0/Vq1dXVFSUXnvttVI/76WXXlJ0dLRq1Kih0NBQ9e7du9S/0j/88ENFR0crKChIzZs314svvnjJWn7piSeeUNOmTVW9enVdd911+vzzz0vs895776lXr16qWbOmQkJCdNNNN7l9z1JTU+Xj46Pnn3/e7XWPPPKIfHx89Pbbb180Q2nXEyxYsEDt2rVz1d+5c+dL/nXerl071a9f321bYGCgbrjhBh0+fFi5ubkXff15Z86c0V/+8hfVq1dPderU0e23364TJ05cMPMHH3yga6+9VtLPRxJ8fHzk4+Pjagb27dunESNGKDw8XEFBQWrcuLFGjRqlnJwct/ccMGCAPvzwQ/34449lygn8Vvg5HQCVb+nSpRo+fLgCAgI0evRoLVq0SNu3b3f9n+kvpaam6uzZs7r77rsVGBiounXrup675ZZb1LJlSz3yyCM6f1f49evX65tvvtH48eMVHh6uL774Qs8884y++OILbdmyRT4+PurTp4+aNGmipUuXatiwYSWyXXnllerWrdsF858+fVqnT58u8UunNJs2bVK9evXUtGlTt+1nzpzRddddp++//15/+ctfFBkZqU2bNikhIUFHjx5VUlKSa98nn3xSN954o8aMGaOCggKlpaXplltu0erVqzV48GDXfrNmzdLMmTPVvXt3zZ49WwEBAdq6davee+89DRw40LXf/v37dfPNN+vOO+/U2LFj9fzzz2vcuHGKiopSu3btLlnTiy++qNzcXE2cOFFnz57Vk08+qX79+mnXrl0KCwuTJL377rsaNGiQmjdvrpkzZ+qnn37SggUL1KNHD+3YsUNXXHGFxo8frzfeeEPx8fEaMGCAmjRpol27dmnWrFm68847dcMNN1wyyy89++yz+utf/6qbb75ZkydP1tmzZ/XZZ59p69at+uMf/2jrvSQpMzNTNWrUUI0aNcq0f1xcnEJCQjRz5kzt3btXixYt0nfffacPPvhAPj4+JfZv06aNZs+erenTp+vuu+9Wr169JEndu3dXQUGBYmNjlZ+fr0mTJik8PFzff/+9Vq9erZMnTyo4ONj1PlFRUbIsS5s2bSr3BbKAkSxUKR9//LElyVq/fr1lWZZVXFxsNW7c2Jo8ebLbfgcOHLAkWXXq1LGys7PdnpsxY4YlyRo9enSJ9z9z5kyJbS+//LIlydq4caNrW0JCghUYGGidPHnStS07O9vy8/OzZsyYcdEaHnroIUuSlZ6efqlyrZ49e1pRUVGlvkfNmjWtr776ym37/fffb1WrVs06ePDgBWsqKCiwrr76aqtfv36ubfv27bN8fX2tYcOGWUVFRW77FxcXu/67adOmJb4X2dnZVmBgoHXvvfdetJbz/ybVq1e3Dh8+7Nq+detWS5I1depU17aOHTtaDRs2tI4fP+7a9umnn1q+vr7W7bff7tp29OhRq27dutaAAQOs/Px8q1OnTlZkZKSVk5Nz0SyWZVnXXXeddd1117m+vummm6x27dpd8nVlsW/fPisoKMj605/+dMl9U1NTLUlWVFSUVVBQ4Nr+6KOPWpKsVatWXTDz9u3bLUlWamqq23t+8sknliTr1VdfveTnHzlyxJJkzZ0799KFAb8hnMKoYpYuXaqwsDD17dtXkuTj46ORI0cqLS1NRUVFJfYfMWKEGjRoUOp73XPPPSW2Va9e3fXfZ8+e1bFjx9S1a1dJ0o4dO1zP3X777crPz3c7FbB8+XKdO3dOt9122wXzb9y4UbNmzdKtt96qfv36XaJa6fjx4woNDS2x/dVXX1WvXr0UGhqqY8eOuR4xMTEqKirSxo0bS63pxIkTysnJUa9evdzqWblypYqLizV9+vQSp1x+/ddv27ZtXX/tSlKDBg3UqlUrffPNN5esR5KGDh2q3/3ud66vo6Oj1aVLF9cph6NHj2rnzp0aN26c2xGj9u3ba8CAAW6nJsLDw5WcnKz169erV69e2rlzp55//nnVqVOnTFl+KSQkRIcPH9b27dttv/aXzpw5o1tuuUXVq1fXnDlzyvy6u+++W/7+/q6vJ0yYID8/v0ueiinN+SMM77zzjs6cOXPRfc//fB07dsz25wAmo4GoQoqKipSWlqa+ffvqwIED2r9/v/bv368uXbooKytL6enpJV7TrFmzC75fac/9+OOPmjx5ssLCwlS9enU1aNDAtd8vzx23bt1a1157rdv1F0uXLlXXrl0vOC2xZ88eDRs2TFdffbWee+65Mtdt/f/TK7+0b98+rV27Vg0aNHB7nL9WIjs727Xv6tWr1bVrVwUFBalu3bpq0KCBFi1a5FbP119/LV9fX7Vt2/aSeSIjI0tsCw0NLXG+/kJatmxZYttVV12lb7/9VpL03XffSZJatWpVYr82bdro2LFjbhe9jho1SoMHD9a2bdt01113qX///mXK8WvTpk1TrVq1FB0drZYtW2rixIn66KOPbL1HUVGRRo0apd27d+u1115TREREmV/76+9LrVq11KhRI9f3xY5mzZopPj5ezz33nOrXr6/Y2FglJyeXuP5B+u/PV2mnSYDfMq6BqELee+89HT16VGlpaUpLSyvx/NKlS93O1Uvuf33/WmnP3Xrrrdq0aZPuu+8+dezYUbVq1VJxcbGuv/76Ehcn3n777Zo8ebIOHz6s/Px8bdmy5YKTFYcOHdLAgQMVHByst99+W7Vr1y5LyapXr16pv5iLi4s1YMAA/e1vfyv1dVdddZUk6T//+Y9uvPFG9e7dWwsXLlSjRo3k7++v1NTUco/uVatWrdTtpTU6leH48eOuRZh2796t4uLiEkdRyqJNmzbau3evVq9erbVr1+r111/XwoULNX36dM2aNatM73HXXXdp9erVWrp0aZmOMHnS448/rnHjxmnVqlVat26d/vrXvyoxMVFbtmxR48aNXfud//kqyzU5wG8JDUQVsnTpUjVs2FDJycklnnvjjTe0YsUKLV68+KJNw8WcOHFC6enpmjVrltv43b59+0rdf9SoUYqPj9fLL7+sn376Sf7+/ho5cmSJ/Y4fP66BAwcqPz9f6enpatSoUZkztW7dWq+//nqJ7VdeeaVOnz5d6nTGL73++usKCgrSO++8o8DAQNf2X69PcOWVV6q4uFi7d+9Wx44dy5yvPEr7fn711Ve64oorJMl1wejevXtL7Ldnzx7Vr19fNWvWdG2bOHGicnNzlZiYqISEBCUlJSk+Pr5c2WrWrKmRI0dq5MiRKigo0PDhw/Xwww8rISFBQUFBF33tfffdp9TUVCUlJWn06NG2P3vfvn2uU3PSzxfbHj169KIXg17qqME111yja665Rg8++KA2bdqkHj16aPHixfrnP//p2ufAgQOSfm6ggKqEUxhVxE8//aQ33nhDf/jDH3TzzTeXeMTFxSk3N1dvvvlmuT/j/F/Wv/5L+pcTDb9Uv359DRo0SC+99JKWLl2q66+/vsRfcXl5ebrhhhv0/fff6+233y718P3FdOvWTSdOnChxfcGtt96qzZs365133inxmpMnT+rcuXOumnx8fNyuD/n222+1cuVKt9cMHTpUvr6+mj17dokjLRV9ZGHlypX6/vvvXV9v27ZNW7du1aBBgyRJjRo1UseOHfXCCy/o5MmTrv0+//xzrVu3zu0X6muvvably5drzpw5uv/++zVq1Cg9+OCD+uqrr2zn+uWorvTzmgtt27aVZVkqLCy86GvnzZunxx57TA888IAmT55s+7Ml6ZlnnnH7nEWLFuncuXOu70tpzjdSv/w+SdKpU6dcPwPnXXPNNfL19VV+fr7b9oyMDPn4+Fx0cgj4LeIIRBXx5ptvKjc3VzfeeGOpz3ft2tW1qFRpRwHKok6dOurdu7ceffRRFRYW6ne/+53WrVvn+gutNLfffrtuvvlmSdJDDz1U4vkxY8Zo27ZtuuOOO/Tll1+6rWNQq1YtDR069KKZBg8eLD8/P7377ru6++67Xdvvu+8+vfnmm/rDH/7gGqHMy8vTrl279Nprr+nbb79V/fr1NXjwYM2fP1/XX3+9/vjHPyo7O1vJyclq0aKFPvvsM9f7tWjRQn//+9/10EMPqVevXho+fLgCAwO1fft2RUREKDExsazfxktq0aKFevbsqQkTJig/P19JSUmqV6+e2+mYefPmadCgQerWrZvuvPNO1xhncHCwZs6cKenn6zwmTJigvn37Ki4uTpL09NNP6/3339e4ceP04Ycf2jqVMXDgQIWHh6tHjx4KCwvTl19+qaefflqDBw++6CmnFStW6G9/+5tatmypNm3a6KWXXnJ7fsCAAa7x1IspKChQ//79deutt2rv3r1auHChevbsecGfeennI0chISFavHixateurZo1a6pLly769NNPFRcXp1tuuUVXXXWVzp07p//7v/9TtWrVNGLECLf3WL9+vXr06KF69epdMiPwm+LkCAgqz5AhQ6ygoCArLy/vgvuMGzfO8vf3t44dO+YaGZw3b16J/c6Pcf7www8lnjt8+LA1bNgwKyQkxAoODrZuueUW15hbaeOZ+fn5VmhoqBUcHGz99NNPJZ4/P/ZY2qNp06Zlqv3GG2+0+vfvX2J7bm6ulZCQYLVo0cIKCAiw6tevb3Xv3t167LHH3MYBU1JSrJYtW1qBgYFW69atrdTUVNf34Neef/55q1OnTlZgYKAVGhpqXXfdda6R2fP1DB48uMTrfj1eWJpf/ps8/vjjVpMmTazAwECrV69e1qefflpi/3fffdfq0aOHVb16datOnTrWkCFDrN27d7ueHz58uFW7dm3r22+/dXvdqlWryjSW+OvM//rXv6zevXtb9erVswIDA60rr7zSuu+++y45Enr+e3mhx/vvv3/R158f49ywYYN19913W6GhoVatWrWsMWPGuI2xlpb5fL1t27a1/Pz8XCOd33zzjXXHHXdYV155pRUUFGTVrVvX6tu3r/Xuu++6vfbkyZNWQECA9dxzz100I/Bb5GNZDl25Bejn5bEjIiI0ZMgQpaSkeOQz/vOf/6hPnz7as2eP7VMgwMUkJSXp0Ucf1ddff13ua4cAU3ENBBy1cuVK/fDDD7r99ts99hm9evXSwIED9eijj3rsM1D1FBYWav78+XrwwQdpHlAlcQQCjti6das+++wzPfTQQ6pfv77bokwAAO/HEQg4YtGiRZowYYIaNmxo+0ZSAADncQQCAADYxhEIAABgGw0EAACwzesWkiouLtaRI0dUu3Ztbk4DAFWYZVnKzc1VREREue7PciFFRUWXXB21qvL397/g/Xp+zWMNRHJysubNm6fMzEx16NBBCxYsUHR09CVfd+TIETVp0sRTsQAAhjl06JDbDczKy7IsZWZmlli6HO5CQkIUHh5+yT/iPdJALF++XPHx8Vq8eLG6dOmipKQkxcbGau/evWrYsOFFX3t+ydueukF+8vdEvErx95en6OHRSU7HKDfT80vm12B6fsn8GsjvrHMq1Icq+913L+V889CwYUPVqFGDo9y/YlmWzpw5o+zsbEm65I0LPdJAzJ8/X3fddZfGjx8vSVq8eLHWrFmj559/Xvfff/9FX3v+H9RP/vLzMbeBqFmjFvkdZnoNpueXzK+B/A77/zOCFfGLvqioyNU8cN+SCzu/KFp2drYaNmx40dMZFX4RZUFBgTIyMtxuk+zr66uYmBht3ry5xP75+fk6deqU2wMAgIp0/pqHGjVqOJzE+53/Hl3qOpEKbyCOHTumoqKiEnfPCwsLU2ZmZon9ExMTFRwc7Hpw/QMAwFM4bXFpZf0eOT7GmZCQoJycHNfj0KFDTkcCAACXUOHXQNSvX1/VqlVTVlaW2/asrCyFh4eX2D8wMFCBgYEVHQMAgDLJPviDco7lVtrnBdevrYaRDSrt8zylwhuIgIAARUVFKT09XUOHDpX089oO6enpiouLq+iPAwCg3LIP/qBxrSer8GzlrQvhH+SvJXueLHMTMW7cOL3wwgtKTEx0G0RYuXKlhg0bpl/ekWLDhg267bbb3I7mb968WT179tT111+vNWvWVFgdHpnCiI+P19ixY9W5c2dFR0crKSlJeXl5rqmMsvj7y1NUs0YtT8SrFK2iW2j2ymlOxyg30/NL5tdgen7J/BrI76y8M6f1wehVHv2MnGO5ldo8SFLh2ULlHMu1dRQiKChIc+fO1V/+8heFhoZecL9Vq1ZpyJAhbttSUlI0adIkpaSk6MiRI4qIiCh39l/ySAMxcuRI/fDDD5o+fboyMzPVsWNHrV27tsSFlRfz8Ogko8ePZq+cpulD5zodo9xMzy+ZX4Pp+SXzayC/s85ZrBZ5XkxMjPbv36/ExEQ9+uijF9zvzTff1NNPP+36+vTp01q+fLk+/vhjZWZmasmSJXrggQcqJJPHLqKMi4vTd999p/z8fG3dulVdunTx1EcBAPCbVq1aNT3yyCNasGCBDh8+XOo+X3zxhbKzs9WvXz/XtldeeUWtW7dWq1atdNttt+n5559XRd2E2/EpDAAAcGnDhg1Tx44dNWPGjFKfX7VqlWJjYxUQEODalpKSottuu02SdP311ysnJ0cbNmyokDw0EAAAGGLu3Ll64YUX9OWXX5Z4btWqVbrxxhtdX+/du1fbtm3T6NGjJUl+fn4aOXKkUlJSKiQLDQQAAIbo3bu3YmNjlZCQ4Lb96NGj+uSTTzR48GDXtpSUFJ07d04RERHy8/OTn5+fFi1apNdff105OTmXncXrbucNAAAubM6cOerYsaNatWrl2vbWW2+pe/fuqlu3riTp3LlzevHFF/X4449r4MCBbq8fOnSoXn75Zd1zzz2XlcNrGwjGOJ1len7J/BpMzy+ZXwP5nVUZY5wmuuaaazRmzBg99dRTrm1vvvmm2+mL1atX68SJE7rzzjsVHBzs9voRI0YoJSXlt9tAMMbpLNPzS+bXYHp+yfwayO8sxjgvbPbs2Vq+fLkkKS8vT+np6UpKSnI9n5KSopiYmBLNg/RzA/Hoo4/qs88+U/v27cudwWsbCAAAPC24fm35B/lX+kqUwfVrl3n/JUuWlNh2xRVXKD8/X5L0xhtvqFmzZmrRooXr+bfeeuuC7xcdHV0ho5w0EACAKqthZAMt2fOk0ffCqFWrlubOrfwjTTQQAIAqrWFkA6NvbvXriyQrC2OcAADANhoIAABgm49VUYtiV5BTp04pODhY619+3/gxzr3b9jsdo9xMzy+ZX4Pp+SXzayC/s/LOnNaA0X2Vk5OjOnXqXNZ7nT17VgcOHFCzZs0UFBRUQQl/m8r6vfLaayAY43SW6fkl82swPb9kfg3kdxZjnN6NUxgAAMA2GggAAGAbDQQAALDNa6+BAACgMmT+cEo5uT9V2ucF166u8AaXd1GoN6CBAABUWZk/nNIfJ6WooLCo0j4zwL+ali2403YTkZycrHnz5ikzM1MdOnTQggULFB0dLennyYl7771XaWlpys/PV2xsrBYuXKiwsDBPlCCJUxgAgCosJ/enSm0eJKmgsMj2EY/ly5crPj5eM2bM0I4dO9ShQwfFxsYqOztbkjR16lS99dZbevXVV7VhwwYdOXJEw4cP90R8F9aB8BDT569Nzy+ZX4Pp+SXzayC/sypjHYi932Tpzvv+73Kj2pYy709q1bzsRwe6dOmia6+9Vk8//bQkqbi4WE2aNNGkSZM0YcIENWjQQMuWLdPNN98sSdqzZ4/atGmjzZs3q2vXrraysQ6Ew0yfvzY9v2R+Dabnl8yvgfzOYh2InxUUFCgjI0MJCQmubb6+voqJidHmzZsVHR2twsJCxcTEuJ5v3bq1IiMjy9VAlBWnMAAA8GLHjh1TUVFRiesZwsLClJmZqczMTAUEBCgkJKTU5z2FBgIAANhGAwEAgBerX7++qlWrpqysLLftWVlZCg8PV3h4uAoKCnTy5MlSn/cUGggAALxYQECAoqKilJ6e7tpWXFys9PR0devWTVFRUfL393d7fu/evTp48KC6devmsVxeexElAAD4WXx8vMaOHavOnTsrOjpaSUlJysvL0/jx4xUcHKw777xT8fHxqlu3rurUqaNJkyapW7duHruAUvLiBuLvL08xfoxz9sppTscoN9PzS+bXYHp+yfwayO+svDOn9cHoVU7H8AojR47UDz/8oOnTpyszM1MdO3bU2rVrXRdWPvHEE/L19dWIESPcFpLyJK9dB6KPbmKM00Gm55fMr8H0/JL5NZDfWeesQn2gVR5dB8KklSgri/HrQAAA4GnhDepo2YI7uRdGOdBAAACqtPAGdX4Tv9ArG1MYAADANhoIAABgGw0EAACwzWuvgTgzpLP8/C989ae3K2oUojPDujgdo9xMzy/9NmoAAG/ltWOca9/boZq1ajsdp9zatGikL/cfdTpGuZmeXzK/hrZN6ht9K2bJ/NtJk99ZlXE7b5Rk/Bjn9CfWGH0EYu79QzVtzkqnY5Sb6fkl82t4bGwfo2f4JfPXISC/s7idt3fjGggAAGCb1x6BAACgMnyfe0onfqq8haRCq1fX72qbv+4EDQQAoMr6PveUYl58XvlFlbeUdWC1anr39jtsNxHJycmaN2+eMjMz1aFDBy1YsEDR0dGSpGeeeUbLli3Tjh07lJubqxMnTigkJMQD6f+LUxgAgCrrxE8/VWrzIEn5RUW2j3gsX75c8fHxmjFjhnbs2KEOHTooNjZW2dnZkqQzZ87o+uuv1wMPPOCJyKWigQAAwMvNnz9fd911l8aPH6+2bdtq8eLFqlGjhp5//nlJ0pQpU3T//fd79Pbdv+a1pzBmTx1s/Bjn3PuHOh2j3EzPL5lfw1UtGmnG1kSnY1wW02sgv7PyTufqg37czrugoEAZGRlKSEhwbfP19VVMTIw2b97sWC6vbSAY43SW6fkl82swPb9kfg3kd9a5wrNOR/AKx44dU1FRkcLCwty2h4WFac+ePQ6l4hQGAAAoBxoIAAC8WP369VWtWjVlZWW5bc/KylJ4eLhDqWggAADwagEBAYqKilJ6erprW3FxsdLT09WtWzfHcnntNRAAAOBn8fHxGjt2rDp37qzo6GglJSUpLy9P48ePlyRlZmYqMzNT+/f/fO+TXbt2qXbt2oqMjFTdunU9kokGAgAALzdy5Ej98MMPmj59ujIzM9WxY0etXbvWdWHl4sWLNWvWLNf+vXv3liSlpqZq3LhxHslEAwEAqLJCq1dXYLVqlb4SZWj16rZfFxcXp7i4uFKfmzlzpmbOnHmZyezhdt4eYvqtpE3PL5lfg+n5JfNrIL+z8k7n6vp+v/f47by5F4Y7buftMNPnr03PL5lfg+n5JfNrIL+zKmsdiN/VruPVv9C9FVMYAADANhoIAABgGw0EAACwjQYCAADYRgMBAABsY4zTQ0wfnzI9v2R+Dabnl8yvgfzOqqwxTrhjjNNhpo9PmZ5fMr8G0/NL5tdAfmdxO2/v5rUNBAAAleF4/nHlnjtdaZ9X26+W6gXWq7TP8xQaCABAlXU8/7imffZ3FVqFlfaZ/j7+mtv+4TI3ERs3btS8efOUkZGho0ePasWKFRo6dKjrecuyNGPGDD377LM6efKkevTooUWLFqlly5YequBnXEQJAKiycs+drtTmQZIKrUJbRzzy8vLUoUMHJScnl/r8o48+qqeeekqLFy/W1q1bVbNmTcXGxursWc+eAuIIBAAAXmzQoEEaNGhQqc9ZlqWkpCQ9+OCDuummmyRJL774osLCwrRy5UqNGjXKY7k4AgEAgKEOHDigzMxMxcTEuLYFBwerS5cu2rx5s0c/2/YRiMo6FzN76mDjxzjn3j/U6RjlZnp+yfwaTM8vmV8D+Z2VdzpX17/1D6djeLXMzExJUlhYmNv2sLAw13OeYruBOH8u5o477tDw4cNLPH/+XMwLL7ygZs2a6R//+IdiY2O1e/duW7O3jHE6y/T8kvk1mJ5fMr8G8juLMU7vZruB8NZzMQAAVDXh4eGSpKysLDVq1Mi1PSsrSx07dvToZ1foNRDlOReTn5+vU6dOuT0AAMClNWvWTOHh4UpPT3dtO3XqlLZu3apu3bp59LMrdAqjPOdiEhMTNWvWrIqMAQDAb8bp06e1f/9+19cHDhzQzp07VbduXUVGRmrKlCn65z//qZYtW7ouHYiIiHC7PtETHB/jTEhIUHx8vOvrU6dOqUmTJg4mAgDAe3z88cfq27ev6+vzvzPHjh2rJUuW6G9/+5vy8vJ099136+TJk+rZs6fWrl3r8Xt+VGgDUZ5zMYGBgQoMDKzIGAAAlEltv1ry9/Gv9JUoa/vVKvP+ffr00cXue+nj46PZs2dr9uzZFRGvzCq0gfjluZjzDcP5czETJkyoyI/yai2m7VbNFv3VYtpup6NchqFOBwAAj6sXWE9z2z/MvTDKwfbtvH95LqZTp06aP3+++vbt6zoXM3fuXM2ZM8dtjPOzzz4r8xjnb+F23jVb5Kp5zeb6Ju8bp6OUW6SuMvo2wJL5tzI2Pb9kfg3kdxa383aGx27nXVnnYkxeB6LFtN2a0nKSkvYtcDpKud3l84DR8+OS+TPwpueXzK+B/M5iHQjvZruB8NZzMQAAoPJwLwwAAGAbDQQAALCNBgIAANhGAwEAAGyzPcbpab+FMU7J/PEp0/NL5tdgen7J/BpatW2kT7M8e0tkT+oUGm70958xTmd4bIyzspg8ximZPz5len7J/BpMzy+ZX8OsecN011srnY5Rbv/X/Sajv/+VNcZpFR2Rik9UymdJknxD5VMtovI+z0O8toEAAMDTrKIjsn4YKKmgEj81QGqwrsxNxMaNGzVv3jxlZGTo6NGjWrFihduNst544w0tXrxYGRkZ+vHHH/XJJ594/FbeEtdAAACqsuITqtzmQT9/no0jHnl5eerQoYOSk5Mv+HzPnj01d+7cigpYJhyBAADAiw0aNEiDBg264PN/+tOfJEnffvttJSX6GUcgAACAbTQQAADANk5hAMAFBB23FLHRqybd7enudAD8lnltAzF76mDj14GYe/9Qp2OUm+n5JfNrMD2/ZH4N5HdW3ulcXf/WP5yOgQvw2gaCdSCcZXp+yfwaTM8vmV8D+Z3F7by9m9c2EAAAQDp9+rT279/v+vrAgQPauXOn6tatq8jISP344486ePCgjhw5Iknau3evJCk8PFzh4eEey8VFlAAAeLGPP/5YnTp1UqdOnSRJ8fHx6tSpk6ZPny5JevPNN9WpUycNHjxYkjRq1Ch16tRJixcv9mgujkAAAKou31BJAar0lSh9Q8u8d58+fXSx21aNGzdO48aNq4Bc9tBAAACqLJ9qEVKDddwLoxxoIAAAVZpPtQjpN/ALvbJxO28PMf02xqbnl8yvwfT8kvk1kN9Z3M7bGdzO22Gmj0+Znl8yvwbT80vm10B+ZzHG6d2YwgAAVBledtDdK5X1e0QDAQD4zfP395cknTlzxuEk3u/89+j89+xCvPYUBgAAFaVatWoKCQlRdna2JKlGjRry8fFxOJV3sSxLZ86cUXZ2tkJCQlStWrWL7k8DAQCoEs6vyni+iUDpQkJCyrSCJQ0EAKBK8PHxUaNGjdSwYUMVFhY6Hccr+fv7X/LIw3mMcXqI6eNTpueXzK/B9PyS+TWQ31kVOcaJiue1RyAY43SW6fkl82swPb9kfg3kdxZjnN6NKQwAAGAbDQQAALCNBgIAANhGAwEAAGyjgQAAALbRQAAAANtYB8JDTJ+/Nj2/ZH4NpueXzK+B/M5iHQjvxjoQHmL6/LXp+SXzazA9v2R+DeR3FutAeDdOYQAAANtoIAAAgG00EAAAwDYaCAAAYBsNBAAAsI0xTg8xfXzK9PyS+TWYnl8yvwbyO4sxTu/GGKeHmD4+ZXp+yfwaTM8vmV8D+Z3FGKd34xQGAACwjQYCAADYRgMBAABso4EAAAC20UAAAADbaCAAAIBtrAPhIabPX5ueXzK/BtPzS+bXQH5nsQ6Ed2MdCA8xff7a9PyS+TWYnl8yvwbyO4t1ILwbpzAAAIBtNBAAAMA2GggAAGAbDQQAALCNBgIAANjGGKeHmD4+ZXp+yfwaTM8vmV8D+Z3FGKd3Y4zTQ0wfnzI9v2R+Dabnl8yvgfzOYozTu3EKAwAA2EYDAQAAbKOBAAAAttFAAAAA22ggAACAbbbGOBMTE/XGG29oz549ql69urp37665c+eqVatWrn3Onj2re++9V2lpacrPz1dsbKwWLlyosLCwMn0GY5zewfT8kvk1mJ5fMr+GVm0b6dOsTKdjlFuHsHCj8+fl5mpox06McXopW2OcGzZs0MSJE3Xttdfq3LlzeuCBBzRw4EDt3r1bNWvWlCRNnTpVa9as0auvvqrg4GDFxcVp+PDh+uijj2wFY4zTWabnl8yvwfT8kvk1zJo3THe9tdLpGOX27JChRucvPssYpzez1UCsXbvW7eslS5aoYcOGysjIUO/evZWTk6OUlBQtW7ZM/fr1kySlpqaqTZs22rJli7p27VpxyQEAgGMu6xqInJwcSVLdunUlSRkZGSosLFRMTIxrn9atWysyMlKbN28u9T3y8/N16tQptwcAAPBu5W4giouLNWXKFPXo0UNXX321JCkzM1MBAQEKCQlx2zcsLEyZmaWfh0tMTFRwcLDr0aRJk/JGAgAAlaTcDcTEiRP1+eefKy0t7bICJCQkKCcnx/U4dOjQZb0fAADwvHLdCyMuLk6rV6/Wxo0b1bhxY9f28PBwFRQU6OTJk25HIbKyshQeHl7qewUGBiowMLA8MQAAgENsNRCWZWnSpElasWKFPvjgAzVr1szt+aioKPn7+ys9PV0jRoyQJO3du1cHDx5Ut27dKi41AFSCoOOWIjZ61Q2LbQnqbnb+c4WWDjodAhdkax2I//mf/9GyZcu0atUqt7UfgoODVb16dUnShAkT9Pbbb2vJkiWqU6eOJk2aJEnatGlTmT6DdSC8g+n5JfNrMD2/ZH4N5HcWt/P2brYaCB8fn1K3p6amaty4cZL+u5DUyy+/7LaQ1IVOYfza+QYieshDrAPhINPzS+bXYHp+yfwayO+sc4Vnte2tf9BAeCnbpzAuJSgoSMnJyUpOTi53KAAA4N24FwYAALCNBgIAANhGAwEAAGyjgQAAALbZmsKoDIxxegfT80vm12B6fsn8GsjvLMY4vVu5VqKsDNzO21mm55fMr8H0/JL5NZDfWecKuZ23N+MUBgAAsI0GAgAA2EYDAQAAbKOBAAAAttFAAAAA22ggAACAbawD4SGmz1+bnl8yvwbT80vm10B+Z7EOhHdjHQgPMX3+2vT8kvk1mJ5fMr8G8juLdSC8G6cwAACAbTQQAADANhoIAABgGw0EAACwjQYCAADYxhinh5g+PmV6fsn8GkzPL5lfA/mdxRind2OM00NMH58yPb9kfg2m55fMr4H8zmKM07txCgMAANhGAwEAAGyjgQAAALbRQAAAANtoIAAAgG2McXqI6eNTpueXzK/B9PyS+TWQ31mMcXo3xjg9xPTxKdPzS+bXYHp+yfwayO8sxji9G6cwAACAbTQQAADANhoIAABgGw0EAACwjQYCAADYRgMBAABsYx0IDzF9/tr0/JL5NZieXzK/BvI7i3UgvBvrQHiI6fPXpueXzK/B9PyS+TWQ31msA+HdOIUBAABso4EAAAC20UAAAADbaCAAAIBtNBAAAMA2xjg9xPTxKdPzS+bXYHp+yfwayO8sxji9G2OcHmL6+JTp+SXzazA9v2R+DeR3FmOc3o1TGAAAwDYaCAAAYBsNBAAAsI0GAgAA2EYDAQAAbPPaKQygqjtbz0dHevs4HaPcIjZ61YR4lWT6z1DxWR/pLadT4EJYB8JDTJ+/Nj2/ZH4Nrdo20qdZmU7HKLeg45bx/wam5zf9ZygvN1dDO3ZiHQgv5bVHIFgHwlmm55fMr2HWvGG6662VTscot4iNlvH/BqbnN/1nqPgs60B4M66BAAAAttFAAAAA22ggAACAbTQQAADANhoIAABgG2OcHmL6+Jfp+SXzazA9v2R+DeR3Frfz9m6McXqI6eNfpueXzK/B9PyS+TWQ31ncztu7cQoDAADYRgMBAABso4EAAAC20UAAAADbaCAAAIBttsY4Fy1apEWLFunbb7+VJLVr107Tp0/XoEGDJElnz57Vvffeq7S0NOXn5ys2NlYLFy5UWFhYmQMxxukdTM8vmV+D6fkl82sgv7MY4/RutsY4GzdurDlz5qhly5ayLEsvvPCCbrrpJn3yySdq166dpk6dqjVr1ujVV19VcHCw4uLiNHz4cH300Ue2gzHG6SzT80vm12B6fsn8GsjvLMY4vZutBmLIkCFuXz/88MNatGiRtmzZosaNGyslJUXLli1Tv379JEmpqalq06aNtmzZoq5du1ZcagAA4KhyXwNRVFSktLQ05eXlqVu3bsrIyFBhYaFiYmJc+7Ru3VqRkZHavHnzBd8nPz9fp06dcnsAAADvZruB2LVrl2rVqqXAwEDdc889WrFihdq2bavMzEwFBAQoJCTEbf+wsDBlZmZe8P0SExMVHBzsejRp0sR2EQAAoHLZbiBatWqlnTt3auvWrZowYYLGjh2r3bt3lztAQkKCcnJyXI9Dhw6V+70AAEDlsH0vjICAALVo0UKSFBUVpe3bt+vJJ5/UyJEjVVBQoJMnT7odhcjKylJ4ePgF3y8wMFCBgYH2kwMAAMdc9joQxcXFys/PV1RUlPz9/ZWenu56bu/evTp48KC6det2uR8DAAC8iK11IBISEjRo0CBFRkYqNzdXy5Yt09y5c/XOO+9owIABmjBhgt5++20tWbJEderU0aRJkyRJmzZtKnMg1oHwDqbnl8yvwfT8kvk1kN9ZrAPh3WydwsjOztbtt9+uo0ePKjg4WO3bt3c1D5L0xBNPyNfXVyNGjHBbSKo8WAfCWabnl8yvwfT8kvk1kN9ZrAPh3Ww1ECkpKRd9PigoSMnJyUpOTr6sUAAAwLtxLwwAAGAbDQQAALCNBgIAANhGAwEAAGyzNcZZGRjj9A6m55fMr8H0/JL5NZDfWYxxejfbK1FWFsY4nWV6fsn8GkzPL5lfA/mdxRind+MUBgAAsI0GAgAA2EYDAQAAbKOBAAAAttFAAAAA22ggAACAbawD4SGmz1+bnl8yvwbT80vm10B+Z7EOhHdjHQgPMX3+2vT8kvk1mJ5fMr8G8juLdSC8G6cwAACAbTQQAADANhoIAABgGw0EAACwjQYCAADYxhinh5g+PmV6fsn8GkzPL5lfA/mdxRind2OM00NMH58yPb9kfg2m55fMr4H8zmKM07txCgMAANhGAwEAAGyjgQAAALbRQAAAANtoIAAAgG2McXqI6eNTpueXzK/B9PyS+TWQ31mMcXo3xjg9xPTxKdPzS+bXYHp+yfwayO8sxji9G6cwAACAbTQQAADANhoIAABgGw0EAACwjQYCAADYRgMBAABsYx0IDzF9/rrH74Olws+cjlFuH+1qavy/wTUdGuibvG+cjnFZmtdsbnQNkbrK6J8h0/83wDoQ3o11IDzE9PnrD5d2knXyHqdjlNu0OX8x/t/g6ef6K2nfAqdjXJYpLScZXcNdPg8Y/TNk+v8GWAfCu3EKAwAA2EYDAQAAbKOBAAAAttFAAAAA22ggAACAbYxxeojp41Om55fMr8H0/JL5NZDfWYxxejfGOD3E9PEp0/NL5tdgen7J/BrI7yzGOL0bpzAAAIBtNBAAAMA2GggAAGAbDQQAALCNBgIAANjmtVMYNd76WH4+/k7HKLdqY/uoxoqtTscoN9PzS7+BGu4f6nQCALggr10HYv3L76tmjVpOxym3VtEttHfbfqdjlJvp+SXza7iqVxujZ/gl89chIL+zWAfCu3ntEYiHRycZfQRi9sppmj50rtMxys30/JL5NczYmmj0DL9k/joE5HcW60B4N66BAAAAttFAAAAA22ggAACAbTQQAADANhoIAABgG2OcHmL6CKHp+SXzazA9v2R+DeR3Vt6Z0xowui9jnF6KMU4PMX2E0PT8kvk1mJ5fMr8G8jvrnFXodARcBKcwAACAbTQQAADANhoIAABgGw0EAACwjQYCAADYdlljnHPmzFFCQoImT56spKQkSdLZs2d17733Ki0tTfn5+YqNjdXChQsVFhZWpvdkjNM7mJ5fMr8G0/NL5tdAfmcxxundyj3GuX37dv3rX/9S+/bt3bZPnTpVa9as0auvvqrg4GDFxcVp+PDh+uijj2y9P2OczjI9v2R+Dabnl8yvgfzOYozTu5XrFMbp06c1ZswYPfvsswoNDXVtz8nJUUpKiubPn69+/fopKipKqamp2rRpk7Zs2VJhoQEAgLPK1UBMnDhRgwcPVkxMjNv2jIwMFRYWum1v3bq1IiMjtXnz5lLfKz8/X6dOnXJ7AAAA72b7FEZaWpp27Nih7du3l3guMzNTAQEBCgkJcdseFhamzMzMUt8vMTFRs2bNshsDAAA4yNYRiEOHDmny5MlaunSpgoKCKiRAQkKCcnJyXI9Dhw5VyPsCAADPsdVAZGRkKDs7W7///e/l5+cnPz8/bdiwQU899ZT8/PwUFhamgoICnTx50u11WVlZCg8PL/U9AwMDVadOHbcHAADwbrZOYfTv31+7du1y2zZ+/Hi1bt1a06ZNU5MmTeTv76/09HSNGDFCkrR3714dPHhQ3bp1q7jUAADAUZd9O+8+ffqoY8eOrnUgJkyYoLfffltLlixRnTp1NGnSJEnSpk2byvR+rAPhHUzPL5lfg+n5JfNrIL+zWAfCu1X47byfeOIJ+fr6asSIEW4LSdnFOhDOMj2/ZH4NpueXzK+B/M5iHQjvdtkNxAcffOD2dVBQkJKTk5WcnHy5bw0AALwU98IAAAC20UAAAADbaCAAAIBtNBAAAMC2yx7jrGiMcXoH0/NL5tdgen7J/BrI7yzGOL1bhY9xVhTGOJ1len7J/BpMzy+ZXwP5ncUYp3fjFAYAALCNBgIAANhGAwEAAGyjgQAAALbRQAAAANtoIAAAgG2sA+Ehps9fm55fMr8G0/NL5tdAfmexDoR3Yx0IDzF9/tr0/JL5NZieXzK/BvI7i3UgvBunMAAAgG00EAAAwDYaCAAAYBsNBAAAsI0GAgAA2MYYp4eYPj5len7J/BpMzy+ZXwP5ncUYp3djjNNDTB+fMj2/ZH4NpueXzK+B/M5ijNO7cQoDAADYRgMBAABso4EAAAC20UAAAADbaCAAAIBtjHF6iOnjU6bnl8yvwfT8kvk1kN9ZjHF6N8Y4PcT08SnT80vm12B6fsn8GsjvLMY4vRunMAAAgG00EAAAwDYaCAAAYBsNBAAAsI0GAgAA2EYDAQAAbGMdCA8xff7a9PyS+TWYnl8yvwbyO4t1ILwb60B4iOnz16bnl8yvwfT8kvk1kN9ZrAPh3TiFAQAAbKOBAAAAttFAAAAA22ggAACAbTQQAADANsY4PcT08SnT80vm12B6fsn8GsjvLMY4vRtjnB5i+viU6fkl82swPb9kfg3kdxZjnN6NUxgAAMA2GggAAGAbDQQAALCNBgIAANhGAwEAAGyjgQAAALaxDoSHmD5/bXp+yfwaTM8vmV8D+Z3FOhDejXUgPMT0+WvT80vm12B6fsn8GsjvLNaB8G6cwgAAALbRQAAAANtoIAAAgG00EAAAwDYaCAAAYBtjnB5i+viU6fkl82swPb9kfg3kdxZjnN6NMU4PMX18yvT8kvk1mJ5fMr8G8juLMU7vxikMAABgGw0EAACwjQYCAADYRgMBAABso4EAAAC22RrjnDlzpmbNmuW2rVWrVtqzZ48k6ezZs7r33nuVlpam/Px8xcbGauHChQoLCytzIMY4vYPp+SXzazA9v2R+DeR3FmOc3s32GGe7du307rvv/vcN/P77FlOnTtWaNWv06quvKjg4WHFxcRo+fLg++ugj28EY43SW6fkl82swPb9kfg3kdxZjnN7NdgPh5+en8PDwEttzcnKUkpKiZcuWqV+/fpKk1NRUtWnTRlu2bFHXrl0vPy0AAPAKtq+B2LdvnyIiItS8eXONGTNGBw8elCRlZGSosLBQMTExrn1bt26tyMhIbd68+YLvl5+fr1OnTrk9AACAd7PVQHTp0kVLlizR2rVrtWjRIh04cEC9evVSbm6uMjMzFRAQoJCQELfXhIWFKTMz84LvmZiYqODgYNejSZMm5SoEAABUHlunMAYNGuT67/bt26tLly5q2rSpXnnlFVWvXr1cARISEhQfH+/6+tSpUzQRAAB4ucsa4wwJCdFVV12l/fv3Kzw8XAUFBTp58qTbPllZWaVeM3FeYGCg6tSp4/YAAADe7bIaiNOnT+vrr79Wo0aNFBUVJX9/f6Wnp7ue37t3rw4ePKhu3bpddlAAAOA9bK0D8b//+78aMmSImjZtqiNHjmjGjBnauXOndu/erQYNGmjChAl6++23tWTJEtWpU0eTJk2SJG3atKnMgVgHwjuYnl8yvwbT80vm10B+Z7EOhHezdQ3E4cOHNXr0aB0/flwNGjRQz549tWXLFjVo0ECS9MQTT8jX11cjRoxwW0iqPFgHwlmm55fMr8H0/JL5NZDfWawD4d1sNRBpaWkXfT4oKEjJyclKTk6+rFAAAMC7cS8MAABgGw0EAACwjQYCAADYZvteGJ52fijknAqlMs+HeJ+8M6eNvgDI9PyS+TWYnl8yvwbyO+ucfs5uY1gQlcjWGGdlOHz4MCtRAgBcDh06pMaNGzsdA7/idQ1EcXGxjhw5otq1a8vHx8cjn3F+uexDhw4ZOVtMfueZXoPp+SXzazA9v+T5GizLUm5uriIiIuTryxl3b+N1pzB8fX0rrdM0fels8jvP9BpMzy+ZX4Pp+SXP1hAcHOyR98Xlo6UDAAC20UAAAADbqmQDERgYqBkzZigwMNDpKOVCfueZXoPp+SXzazA9v/TbqAHl53UXUQIAAO9XJY9AAACAy0MDAQAAbKOBAAAAttFAAAAA26pcA5GcnKwrrrhCQUFB6tKli7Zt2+Z0pDLbuHGjhgwZooiICPn4+GjlypVOR7IlMTFR1157rWrXrq2GDRtq6NCh2rt3r9OxbFm0aJHat2/vWjinW7du+ve//+10rHKbM2eOfHx8NGXKFKejlMnMmTPl4+Pj9mjdurXTsWz7/vvvddttt6levXqqXr26rrnmGn388cdOxyqTK664osS/gY+PjyZOnOh0NFSyKtVALF++XPHx8ZoxY4Z27NihDh06KDY2VtnZ2U5HK5O8vDx16NBBycnJTkcplw0bNmjixInasmWL1q9fr8LCQg0cOFB5eXlORyuzxo0ba86cOcrIyNDHH3+sfv366aabbtIXX3zhdDTbtm/frn/9619q376901FsadeunY4ePep6fPjhh05HsuXEiRPq0aOH/P399e9//1u7d+/W448/rtDQUKejlcn27dvdvv/r16+XJN1yyy0OJ0Ols6qQ6Ohoa+LEia6vi4qKrIiICCsxMdHBVOUjyVqxYoXTMS5Ldna2JcnasGGD01EuS2hoqPXcc885HcOW3Nxcq2XLltb69eut6667zpo8ebLTkcpkxowZVocOHZyOcVmmTZtm9ezZ0+kYFWby5MnWlVdeaRUXFzsdBZWsyhyBKCgoUEZGhmJiYlzbfH19FRMTo82bNzuYrOrKycmRJNWtW9fhJOVTVFSktLQ05eXlqVu3bk7HsWXixIkaPHiw2/8eTLFv3z5FRESoefPmGjNmjA4ePOh0JFvefPNNde7cWbfccosaNmyoTp066dlnn3U6VrkUFBTopZde0h133OGxmx/Ce1WZBuLYsWMqKipSWFiY2/awsDBlZmY6lKrqKi4u1pQpU9SjRw9dffXVTsexZdeuXapVq5YCAwN1zz33aMWKFWrbtq3TscosLS1NO3bsUGJiotNRbOvSpYuWLFmitWvXatGiRTpw4IB69eql3Nxcp6OV2TfffKNFixapZcuWeueddzRhwgT99a9/1QsvvOB0NNtWrlypkydPaty4cU5HgQO87m6cqBomTpyozz//3Ljz15LUqlUr7dy5Uzk5OXrttdc0duxYbdiwwYgm4tChQ5o8ebLWr1+voKAgp+PYNmjQINd/t2/fXl26dFHTpk31yiuv6M4773QwWdkVFxerc+fOeuSRRyRJnTp10ueff67Fixdr7NixDqezJyUlRYMGDVJERITTUeCAKnMEon79+qpWrZqysrLctmdlZSk8PNyhVFVTXFycVq9erffff7/Sbt1ekQICAtSiRQtFRUUpMTFRHTp00JNPPul0rDLJyMhQdna2fv/738vPz09+fn7asGGDnnrqKfn5+amoqMjpiLaEhIToqquu0v79+52OUmaNGjUq0Wy2adPGuFMx3333nd599139+c9/djoKHFJlGoiAgABFRUUpPT3dta24uFjp6enGnb82lWVZiouL04oVK/Tee++pWbNmTkeqEMXFxcrPz3c6Rpn0799fu3bt0s6dO12Pzp07a8yYMdq5c6eqVavmdERbTp8+ra+//lqNGjVyOkqZ9ejRo8T48ldffaWmTZs6lKh8UlNT1bBhQw0ePNjpKHBIlTqFER8fr7Fjx6pz586Kjo5WUlKS8vLyNH78eKejlcnp06fd/tI6cOCAdu7cqbp16yoyMtLBZGUzceJELVu2TKtWrVLt2rVd154EBwerevXqDqcrm4SEBA0aNEiRkZHKzc3VsmXL9MEHH+idd95xOlqZ1K5du8Q1JzVr1lS9evWMuBblf//3fzVkyBA1bdpUR44c0YwZM1StWjWNHj3a6WhlNnXqVHXv3l2PPPKIbr31Vm3btk3PPPOMnnnmGaejlVlxcbFSU1M1duxY+flVqV8j+CWnx0Aq24IFC6zIyEgrICDAio6OtrZs2eJ0pDJ7//33LUklHmPHjnU6WpmUll2SlZqa6nS0Mrvjjjuspk2bWgEBAVaDBg2s/v37W+vWrXM61mUxaYxz5MiRVqNGjayAgADrd7/7nTVy5Ehr//79Tsey7a233rKuvvpqKzAw0GrdurX1zDPPOB3JlnfeeceSZO3du9fpKHAQt/MGAAC2VZlrIAAAQMWhgQAAALbRQAAAANtoIAAAgG00EAAAwDYaCAAAYBsNBAAAsI0GAgAA2EYDAQAAbKOBAAAAttFAAAAA22ggAACAbf8PPf4DMDRaV/8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 500x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "clusterSlice_ex = recon3d.iloc[0].to_numpy().reshape(20,13,21)[9]\n",
    "time_label_ex = 200*9+200\n",
    "#final_in_ex = reshape_pixelav(clusterSlice_ex, 0, time_label_ex)\n",
    "bitstring_ex = make_bit_string(clusterSlice_ex, 0, time_label_ex)\n",
    "array2_ex = make_array2(bitstring_ex, imageArray=True, zoomImage=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(256, 16)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(array2_ex)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Scan chain to Shift Register conversion:**  \n",
    "The Shift Register data is a series of 2-bit representations per pixel that starts in the top right corner of the array (figure 5) and goes RT-RB-LT-LB. Meanwhile, the Scan Chain takes 3 bits per pixel (will add this transition soon), starts with the right side of the array (pixel 63), and goes in a LB-RB-RT-LT order. This code will translate between those two formats so that we can tell which spot in the array is coming from which pixel. \n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_string = [0]*512\n",
    "test_string[0], test_string[-1] = 1, 1\n",
    "test_string = ''.join(map(str, test_string))\n",
    "\n",
    "def shiftreg_to_scanchain(bit_string):\n",
    "    big_row_len = 16*2*2*2\n",
    "    row_len = int(big_row_len/2)\n",
    "    big_row_start = 0\n",
    "    scanchain_map = []\n",
    "    while big_row_start < big_row_len*4:\n",
    "        section_row0 = bit_string[big_row_start:(big_row_start + row_len)][::-1]\n",
    "        section_row1 = bit_string[big_row_start+row_len:big_row_start+row_len*2][::-1]\n",
    "\n",
    "        move_right = 0\n",
    "        while move_right < 64:\n",
    "            idx = 0 + move_right\n",
    "            lb = section_row1[idx] + section_row1[idx+1]\n",
    "            rb = section_row1[idx+2] + section_row1[idx+3]\n",
    "            rt = section_row0[idx+2] + section_row0[idx+3]\n",
    "            lt = section_row0[idx] + section_row0[idx+1]\n",
    "\n",
    "            square = lb + rb + rt + lt\n",
    "            scanchain_map.append(square)\n",
    "            move_right += 4\n",
    "        big_row_start += big_row_len\n",
    "    scanchain_map = ''.join(scanchain_map)\n",
    "    \n",
    "    # assuming thresholds are the same as in the 2-bit version \n",
    "    scanchain_map_final = []\n",
    "    for i in range(0, len(scanchain_map), 2):\n",
    "        pixel = scanchain_map[i]+scanchain_map[i+1]\n",
    "        if pixel == '00':\n",
    "            newpixel = '000'\n",
    "        elif pixel == '01':\n",
    "            newpixel = '001'\n",
    "        elif pixel == '10':\n",
    "            newpixel = '011'\n",
    "        else:\n",
    "            newpixel = '111'\n",
    "        scanchain_map_final.append(newpixel)\n",
    "    return(''.join(scanchain_map_final))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}