aare/python/tests/ClusterFinderCUDA.ipynb

{
    "cells": [
        {
            "cell_type": "markdown",
            "id": "title",
            "metadata": {},
            "source": [
                "# ClusterFinder — CPU vs CUDA\n",
                "\n",
                "Two finders over the same frames and the same pedestal:\n",
                "\n",
                "| | |\n",
                "|---|---|\n",
                "| **`ClusterFinderMT`** + `ClusterCollector` | multithreaded CPU baseline |\n",
                "| **`ClusterFinderCUDA.find_clusters_batched`** | batched across CUDA streams, over a pinned input buffer, so H2D / kernel / D2H overlap |\n",
                "\n",
                "**Part 1** times them. **Part 2** compares what they found."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 1,
            "id": "imports",
            "metadata": {},
            "outputs": [],
            "source": [
                "import sys; sys.path.append('/home/ferjao_k/aare/build')\n",
                "\n",
                "import time\n",
                "from pathlib import Path\n",
                "\n",
                "import boost_histogram as bh\n",
                "import matplotlib.pyplot as plt\n",
                "import numpy as np\n",
                "from tqdm import tqdm\n",
                "\n",
                "from aare import (File, ClusterFinderMT, ClusterCollector, ClusterFinderCUDA)"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "md-config",
            "metadata": {},
            "source": [
                "## Configuration"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 2,
            "id": "config",
            "metadata": {},
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "Image size:      (400, 400)   cluster (3, 3)\n",
                        "Pedestal frames: 600\n",
                        "Data frames:     100000  of 100000 in file\n"
                    ]
                }
            ],
            "source": [
                "base = Path('/mnt/sls_det_storage/moench_data/2603_MaxIVBeamtime/2026032408/process/xrf/')\n",
                "f = File(base / 'Cu_factor_10_data_master_0.json')\n",
                "pd = File(base / 'Cu_factor_10_pedestal_master_0.json')\n",
                "\n",
                "n_frames_pd  = 600                  # pedestal frames, taken from the head of the file\n",
                "N            = 100000                 # data frames, read after those\n",
                "cluster_size = (3, 3)\n",
                "rows, cols   = f.rows, f.cols\n",
                "image_size   = (rows, cols)\n",
                "\n",
                "N_SIGMA    = 5\n",
                "N_STREAMS  = 4                       # H2D / kernel / D2H overlap\n",
                "BATCH_SIZE = 2000                    # frames per find_clusters_batched call\n",
                "\n",
                "# CPU workers. ClusterFinderMT gives each thread its own ClusterFinder with its\n",
                "# own pedestal, and data frames are handed out round-robin -- so at high thread\n",
                "# counts each pedestal tracks on only 1/n of the frames while the CUDA finder's\n",
                "# single pedestal sees all of them, and the cluster counts drift apart. Keep it\n",
                "# low when comparing counts; raise it to measure CPU throughput.\n",
                "N_THREADS = 24\n",
                "\n",
                "# Fixed per-frame output slot the kernel writes into, copied WHOLE on every D2H.\n",
                "# Too low does not error -- it TRUNCATES. The CUDA cell prints the peak actually\n",
                "# reached, which is how you size this for a new dataset.\n",
                "CAP = 3000\n",
                "\n",
                "capacity = 3_000                    # initial ClusterVector capacity, CPU side\n",
                "N_BINS   = 200\n",
                "E_MIN, E_MAX = -2, 4000              # spectrum range, ADU\n",
                "\n",
                "print(f'Image size:      {image_size}   cluster {cluster_size}')\n",
                "print(f'Pedestal frames: {n_frames_pd}')\n",
                "print(f'Data frames:     {N}  of {f.total_frames} in file')"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "md-build",
            "metadata": {},
            "source": [
                "## Build the two finders"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 3,
            "id": "build",
            "metadata": {},
            "outputs": [],
            "source": [
                "cf_cpu = ClusterFinderMT(image_size, #only positional argument here i.e. need to use the keyword for next ones\n",
                "                         cluster_size=cluster_size, \n",
                "                         n_sigma=N_SIGMA,\n",
                "                         capacity=capacity, \n",
                "                         min_pedestal_samples=n_frames_pd, \n",
                "                         n_threads=N_THREADS)\n",
                "\n",
                "sink = ClusterCollector(cf_cpu)\n",
                "\n",
                "cf_cuda = ClusterFinderCUDA(image_size, \n",
                "                            cluster_size=cluster_size, \n",
                "                            n_sigma=N_SIGMA,\n",
                "                            max_clusters_per_frame=CAP, \n",
                "                            n_streams=N_STREAMS, \n",
                "                            min_pedestal_samples=n_frames_pd, \n",
                "                            algorithm=\"fixed_window\", \n",
                "                            time_kernels=False)"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "md-ped",
            "metadata": {},
            "source": [
                "## Pedestal\n",
                "\n",
                "Identical frames into both, so the implementation is the only thing that differs."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 4,
            "id": "train",
            "metadata": {},
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "Pedestal (600 frames): 1.48s\n"
                    ]
                }
            ],
            "source": [
                "pd.seek(0)\n",
                "t0 = time.perf_counter()\n",
                "for _ in range(n_frames_pd):\n",
                "    pd_img = pd.read_frame()\n",
                "    cf_cpu.push_pedestal_frame(pd_img.copy())\n",
                "    cf_cuda.push_pedestal_frame(pd_img.copy())\n",
                "print(f'Pedestal ({n_frames_pd} frames): {time.perf_counter() - t0:.2f}s')"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "md-io",
            "metadata": {},
            "source": [
                "## Read the data frames\n",
                "\n",
                "Kept out of both timing loops: the two finders run over the same in-memory array."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 5,
            "id": "io",
            "metadata": {},
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "Reading 100000 frames: 50.61s  (1,976 FPS, 32.00 GB at 0.63 GB/s)\n"
                    ]
                }
            ],
            "source": [
                "f.seek(0)\n",
                "t0 = time.perf_counter()\n",
                "data = f.read_n(N)\n",
                "t_io = time.perf_counter() - t0\n",
                "gb = f.bytes_per_frame * N / 1e9\n",
                "print(f'Reading {N} frames: {t_io:.2f}s  ({N / t_io:,.0f} FPS, '\n",
                "      f'{gb:.2f} GB at {gb / t_io:.2f} GB/s)')"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "md-cpu",
            "metadata": {},
            "source": [
                "## 1 — CPU, multithreaded"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 6,
            "id": "cpu",
            "metadata": {},
            "outputs": [
                {
                    "name": "stderr",
                    "output_type": "stream",
                    "text": [
                        "100%|█| 100000/100000 [00:2\n"
                    ]
                },
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "CPU MT (24 threads): 25.73s   3,887 FPS   233,104,694 clusters (2331.0/frame)\n"
                    ]
                }
            ],
            "source": [
                "t0 = time.perf_counter()\n",
                "# Pass the frame number: find_clusters(frame) defaults it to 0, and every\n",
                "# resulting cluster is then stamped frame 0. Part 2 groups by it.\n",
                "for i, frame in enumerate(tqdm(data)):\n",
                "    cf_cpu.find_clusters(frame, i)\n",
                "t_cpu = time.perf_counter() - t0\n",
                "\n",
                "cf_cpu.stop(); sink.stop()\n",
                "clusters_cpu = sink.steal_clusters()\n",
                "\n",
                "hist_cpu = bh.Histogram(bh.axis.Regular(N_BINS, E_MIN, E_MAX))\n",
                "n_clusters_cpu = 0\n",
                "for cv in clusters_cpu:\n",
                "    hist_cpu.fill(cv.sum())\n",
                "    n_clusters_cpu += cv.size\n",
                "\n",
                "print(f'CPU MT ({N_THREADS} threads): {t_cpu:.2f}s   {N / t_cpu:,.0f} FPS   '\n",
                "      f'{n_clusters_cpu:,} clusters ({n_clusters_cpu / N:.1f}/frame)')"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "md-cuda",
            "metadata": {},
            "source": [
                "## 2 — CUDA, batched + multi-streamed\n",
                "\n",
                "`register_input_buffer` page-locks the input once so H2D runs at full DMA\n",
                "bandwidth instead of through the driver's staging buffer. It is paid once,\n",
                "outside the loop, and released when the dataset is done."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 7,
            "id": "cuda",
            "metadata": {},
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "CUDA batched: 5.54s   18,057 FPS   233,095,737 clusters (2331.0/frame)\n",
                        "  input pinned in 0.82s, once, outside the loop\n",
                        "  vs CPU: 4.65x\n",
                        "  peak clusters in one frame: 2,545 of CAP 3,000\n"
                    ]
                }
            ],
            "source": [
                "t0 = time.perf_counter()\n",
                "cf_cuda.register_input_buffer(data)\n",
                "t_pin = time.perf_counter() - t0\n",
                "\n",
                "clusters_cuda = []\n",
                "t0 = time.perf_counter()\n",
                "for start in range(0, N, BATCH_SIZE):\n",
                "    stop = min(start + BATCH_SIZE, N)\n",
                "    clusters_cuda.extend(\n",
                "        cf_cuda.find_clusters_batched(data[start:stop], first_frame=start))\n",
                "t_cuda = time.perf_counter() - t0\n",
                "\n",
                "cf_cuda.unregister_input_buffer()\n",
                "\n",
                "hist_cuda = bh.Histogram(bh.axis.Regular(N_BINS, E_MIN, E_MAX))\n",
                "energies = [np.asarray(cv.sum()).ravel() for cv in clusters_cuda if cv.size > 0]\n",
                "if energies:\n",
                "    hist_cuda.fill(np.concatenate(energies))\n",
                "n_clusters_cuda = sum(cv.size for cv in clusters_cuda)\n",
                "\n",
                "print(f'CUDA batched: {t_cuda:.2f}s   {N / t_cuda:,.0f} FPS   '\n",
                "      f'{n_clusters_cuda:,} clusters ({n_clusters_cuda / N:.1f}/frame)')\n",
                "print(f'  input pinned in {t_pin:.2f}s, once, outside the loop')\n",
                "print(f'  vs CPU: {t_cpu / t_cuda:.2f}x')\n",
                "\n",
                "# One ClusterVector per frame, so this IS the per-frame maximum.\n",
                "peak = max((cv.size for cv in clusters_cuda), default=0)\n",
                "print(f'  peak clusters in one frame: {peak:,} of CAP {CAP:,}')\n",
                "if peak >= CAP:\n",
                "    print('  !! CAP REACHED -- clusters were dropped; raise CAP and re-run.')"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "md-summary",
            "metadata": {},
            "source": [
                "## Summary"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 8,
            "id": "summary",
            "metadata": {},
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "                             s        FPS   us/frame      clusters\n",
                        "CPU MT (24 thr)          25.73      3,887      257.3   233,104,694\n",
                        "CUDA batched              5.54     18,057       55.4   233,095,737\n",
                        "\n",
                        "speedup: 4.65x\n"
                    ]
                }
            ],
            "source": [
                "print(f'{\"\":<22} {\"s\":>7} {\"FPS\":>10} {\"us/frame\":>10} {\"clusters\":>13}')\n",
                "for name, t, n in (('CPU MT (%d thr)' % N_THREADS, t_cpu, n_clusters_cpu),\n",
                "                   ('CUDA batched', t_cuda, n_clusters_cuda)):\n",
                "    print(f'{name:<22} {t:7.2f} {N / t:10,.0f} {t * 1e6 / N:10.1f} {n:>13,}')\n",
                "print(f'\\nspeedup: {t_cpu / t_cuda:.2f}x')"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "md-part2",
            "metadata": {},
            "source": [
                "# Part 2 — Where the two disagree\n",
                "\n",
                "The two implementations are not bit-identical by construction: the CUDA finder\n",
                "updates its pedestal once per frame, the CPU finder per pixel, so once the\n",
                "pedestals drift apart a few pixels land on opposite sides of the threshold.\n",
                "\n",
                "Three views, cheapest first: the **spectra**, the **cluster sets**, and **where\n",
                "in the run** the difference accumulates."
            ]
        },
        {
            "cell_type": "markdown",
            "id": "md-spec",
            "metadata": {},
            "source": [
                "## Spectra\n",
                "\n",
                "Both curves should lie on top of each other; the ratio panel is the sensitive\n",
                "view."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 9,
            "id": "spectra",
            "metadata": {},
            "outputs": [
                {
                    "data": {
                        "image/png":
                            "iVBORw0KGgoAAAANSUhEUgAAA3kAAAKACAYAAAAhGrGwAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAlq5JREFUeJzt3Ql4E9XawPE3TVtAhMomgiziAgooLqCCVwUXEARxB68iKPiB4ILcq4JVSnHBFVFZFDfcRa+C3iuo4IaKKKAoAoooiCKyyybQNpnveU+ZkKTpnnYmk//veYZmmUxOTk7CvDnnvMdnWZYlAAAAAABPSHG6AAAAAACA+CHIAwAAAAAPIcgDAAAAAA8hyAMAAAAADyHIAwAAAAAPIcgDAAAAAA8hyAMAAAAADyHIAwAAAAAPIcgDAAAAAA8hyAPgGt99951cddVV0qxZM6latarsv//+cvzxx8v9998vmzdvDu3XsWNHs1WUiRMnypQpUyrs+EgctIVIt99+u3Tv3l0OPvhg8fl80q9fv0Lr7pdffpELL7xQDjjgAPNZPvvss+Xrr7+Oue+rr74qxx57rPncN2zYUIYOHSo7duwo8ftU0sfrbXqf7qP76mP0sbFoWc866yxTdn0N+lr0NZWEfl/17t1bDjzwQFNP559/folfCwDEA0EeAFd48skn5YQTTpD58+fLzTffLO+++65MmzZNLrnkEnn88celf//+lVYWTuxBW4jt4Ycflk2bNsl5550n6enphTaUDRs2yKmnnirLly+XZ555Rl577TXZvXu3+XHmxx9/jNj3pZdekssuu0zatWsnM2fOlKysLPMjiwZVJVGax+ttzz33nNlH99XH6GNffvnliP1++OEHU9acnBxTdn0N+lr0NelrK86dd95pvr+0vr744gvzQxUAVCoLABw2d+5cy+/3W+ecc461e/fuAvfv2bPHeuutt0LXTz/9dLNVlFatWsX9+MFg0Pr777+tRJFIZa1IpWkLOTk5Vm5uruVlgUAgdLl69epW3759Y+538803W2lpadaqVatCt23dutWqW7eudemll4Zuy8vLsxo0aGB17tw54vEvvfSSpacoM2bMKLI8pXn8O++8Y257+eWXI/Y9++yzrYYNG5pj2S655BJTVi2zTV+LvqZbbrnFKs5ZZ51lHXXUUcXup88Z6zsPAMqLnjwAjrvnnnvMkKbJkydLlSpVCtyvPQbac1CYjz/+2Dxe/4ZbtWqVuT186KUOt9JhVDpcS5+rfv36cuaZZ8qiRYvM/YcccogsWbJEPvnkE/NY3fQ227Zt2+Tf//63GVKq5dJhazr8a+fOnRHPrY+77rrrTC/kUUcdZZ5LexCKMnXqVGnfvr1Ur17dDBHr0qWLfPPNNxH76PA4vW/FihXSrVs3c7lx48byr3/9S/bs2ROxr/ZC3HXXXXLkkUea569Xr54ZDhvdE6GvT4fgvfnmm3LccceZYWzZ2dnmPq2Lzp07y3777WceP2TIEHnnnXci6lt7LVJTU+W3334r8JquvvpqqVOnjunFKUxx70l4GbV35JhjjjFlPPTQQ+XRRx8tcLySvkfBYFAee+wxM2SvWrVqZkjeySefLG+//XaxbcFucy+88IKpe30OLbu+L6NGjTL3RdN2qLdru4x+Xf/73/9M3Ws5tL3odfsxel3bxIknnigLFiwQJ6WklOy0Qd+nM844Q5o2bRq6rWbNmqYn7b///a/k5eWZ2+bNmydr16417TKc9uBr29bjFKU0j9fLepveF04f+8cff8iXX35prmvZtP4vuugiU2abvpZOnToVWSb7O2f27NmybNmyULvR9mLfp716+rnU9qlt5qOPPjKfD21H2hYzMjKkdu3a5rvgrbfeKvAc9nfLs88+Ky1atDBtpm3btqYuLMuSBx54wBxbX6u+B9omo2n59DOmr08/26eccop88MEHEfvo98T//d//me8X+/tD99PHAnC/VKcLACC5BQIB+fDDD81QTT2ZqGgaGOlz6olWkyZNZOPGjTJ37lz566+/zP16AnfxxRebEy0dtqnswPPvv/+W008/XX7//Xe57bbbTLChQcDIkSNl8eLF5uQn/OR++vTp8umnn5r7DzroIDM/p6hAV+c76Qmn/tUATU/WdHjYV199JS1btgztm5uba4JeHcKqJ4Zz5swxgZaWWZ/LDmB69uxpnv+WW26RDh06yK+//mqGqekwNA0W9OQwfP6RnpTqc+sJogYVevKsr1cvT5o0yZT/lVdeMSeY4QYOHCh33323PPHEE+bkNXxeks530v01KCvre2LToE+DNQ2itD51mN6NN95o6kqDutK+Rxowv/jii6YeR48ebQJCrQc7CCuqLdhGjBhhTsY1mNcAqKj3uDDffvutOU5mZqZ5Lg2wNRjS2/TE2/4R5NZbbzUB4cqVKyPeO7fZtWuX/Pzzz3LBBRcUuE/fD71fA/vmzZvL999/H7o9XFpamvlxwr6/MKV5vF7WgFl/kIguk32/fk607FrG6GPa+86aNcsEZbHadIMGDczwzMGDB8vWrVtNG1X6+bXnFesPE/raH3zwQRNkHXHEEeYHGr1f27H+YKBtWtuqtgMN5q688sqI59EgVH8Auvfee0Nt49xzz5W+ffuauh0/frx5/mHDhplgVT87drvXNq/H0+8H/eFJ60o/u/qj0nvvvWeCP9WnTx/zedDPtpZXP496XYfrAkgA5e4LBIBy+PPPP80Qqt69e5f4MdHDNT/66CNzDP0bbuXKleb2Z5991lzfuHGjuT5u3LgyDdEbM2aMlZKSYs2fPz/i9v/85z8Fhobp9YyMDGvz5s3Fvp7Vq1dbqamp1vXXXx9x+/bt262DDjooYnibDo/TY7/22msR+3br1s1q0aJF6Porr7xi9nvjjTci9tOy6+0TJ04M3da0aVMzXPbHH38sMOTO5/NZS5Ysibi9S5cuBepby3XggQeaobW2++67z9SXvg+FKel7omXUsixatKjAULuaNWtaO3fuLNV7NGfOHHM9MzOzTG3BbnOnnXZagfuysrLMfdG0Hert4fWhr6tatWrW77//HrpNX6Pup8MQ7delpk+fbm5/++23LTcobLjmmjVrTDn1vYimQyX1Ph2ire6++25zfe3atQX21SGYzZs3L7IMpXn8EUccYdputD/++MMc45577jHXP//8c3NdP0PRdB+9Tx9TFG0z2nZifR8ddthhZmhvccM4dehv//79reOOOy7iPj2Gfi/s2LGjQNs49thjzdBwm36u9PbvvvvOXNf2VLt2batHjx4FhuG2adPGOvHEE0O37b///tbQoUOLLCcA92K4JoCkoUOgDjvsMNNDNnbsWPNLuPZ4lZT+et66dWszpEqHdNmb/gIea7ioDpWqVatWscfVX8/1OPrrevhxtadAe6Wij6vP1aNHjwI9DNpTF15WHX6o+4UfU8uuvWDRx9TH66/14XSYor7e8F5EpYkqommP2vr16+X1118317VetfdPexfCh7uW5z1p1aqVtGnTJuK2f/7zn2Z4pp21saTvkSbdUDr8tDy0l6S8tKzae2PT3ialPa46lC769vD3OZbw113aLZ5iDVkt7L7C9i3qGGV5fDzKVJpyxaK98Np7Fk0/OzocUodZam+j7vP000+bHvZoOmxUe9ij20bXrl0jyhbdZrSHXHsMtccv/H3Xz9w555xjEl/Zw5p1eLAOF9beeR0KqiMIACQOgjwAjqpbt645kdUhaBVNT350+Jue8OvQQF2eQeeZ3HDDDbJ9+/ZiH79u3TqzzIOefIVvNWrUMHNhdJhh9NCtktDjKs30F31snacXfVytr+ihYjqMMHzemx5Th1fpEMToY/75558lKqsOy9L5cdFi3abzyXRo6YQJE0LBlg57jB7aWZ73RIPTaPZt9hCykr5HOt/I7/fHPGZplPQ9LooGuuHsrJWF3V7U/EYV/dpLs8WD/rCh72usYX32kEX7tel8TVXYvtF1EK00j9d941EmfW36A0o824zOh7300ktNsK/DKXXIpwZcOqc11vtd1jZjf9foMOTo9/6+++4znxG7PvS7R4PBp556ygxJ1mPrD1H6/QHA/ZiTB8BReqKtc0C0Z0XnUTVq1KjUx7ADnujEI9GBjJ08QX8dV5oSXdOj6xwvnQOj86qKC0h1LpSmUy/s/rL82m8/7j//+U9Eoory0GPqyaouRRGLBj3FlVUfb58UhivsJE8DM01qob1qOidIewZ1bbTilPQ9ifW89m32iXlJ3yMNJHUeoD6+PIFarHoLb4/hc/hitceKoMGBk7T+Dz/8cDMHMprepvdr0hx19NFHh24P7zHW3iVdxiBWr3G40jxe99U5pXpf+Lw8u5zaA6y0Z1nLWFj59bUVNce0LG1GAzudC6uBVfj90d9p5WW3f004pEmGYrF/xNF9x40bZ7bVq1ebhETDhw83PfaFfa8AcA+CPACO0wQTM2bMkGuuucZkk4tef0uHCelJRfQQRZs9HFB7cLRHyGZnSSyMBiGaaOSNN96IWKRZT8w18UI0TXqhSTA0oNATsnjRMutJpyZ8iMfwP7usmvREA5mTTjqpTMfQoaKaHGLp0qURJ9CFLR6tiTY0cYomg9GhnrpGWGmHtRX2nihNoKJJSsKHbOr6Zhqwag9gad4jHdY2ZswYM6RUk64UprC2UJTw9qi9szbNKlkZNNOi07QtaHCgGVfthEraM6s9Vjpc0Q6ytG1qkK3DAnv16hV6vP7goQuXF7dWXmker2XS9Ti1bYXvq8lHNLOr/TnRsul3jZZVe5ftH0Q00NFMmDfddJPEm35O9Hsv/POiP0DEyq5ZHjocVHsh9TNdXC97OP1c6/7a6/7555/HtUwAKgZBHgDH6VAgPdnWjHSaZfPaa6818680uNM5Wrq0gv7KXliQp0PuzjrrLHPSrkPFtGdIT0b0JC2cnnTriYr2NmlGOz2p0syeerv+Qh3+i78GMvqruvY46K/2eptmdtQTxNNOO82c6Ok8Np3Loid/77//vgluyhJQaVCggYZmV9TMeDo3Rl+H9qJpZk2de2MvaVBSuiSBZvbTzJU6X07n1+iQLO0t1RNVzawXK/thOH292iOmAZGWT3/h16BKe0hipdPXXlmd46aZ/rTMmr2yOCV9T5SeiGuAoL18emKvvR+a6VCHmdlz10r6HunQUs0eqPONtJ41ONSATtubHuv6668vsi0URetch7bZWTs1aNAgJNYSE4lGg3d7CQ79AUHnemlAZf8ooD2kSrNE6vISOidT60DrVjNB6rBBff/C24wGUvpeaJZW7Xn76aefTEZY7QXWz0I4DYLC56mW5vHajvU2/X7ReZzaI6c9e/oDkrYlPZZNP28aoGu70Hao5dYMrdq7pW0o3uwlTPQ7UIdSalvRjLnazvX1xIvO99NePB2GqcMy9bk0I6y+p/oDiv7V72LNzKnz/nTOq2Yp1UBXe4i1rkq6SD0Ahzmd+QUAwrMKara+Jk2aWOnp6SZ7n2aWGzlypLV+/foiF0PX7HoXX3yxyRynWS2vuOIKa8GCBRHZNdetW2f169fPOvLII82xNXvcMcccYz388MMRCyHroseama9GjRrm8ZoB0aYZ7W6//XaTyVLLqM919NFHWzfddJPJFGrTxw0ZMqRUb65myOvUqZPJFlmlShXzvPqaZs+eHdpH60fLXpKMjpqd78EHHzRZ86pWrWper772gQMHWj/99FNoP32ec889N2aZvv/+e7Owsz5e61az/T333HPmub799tsC+2vd6X2DBg0q0Wsu6Xtil1GzZGrWQq37Qw45xBo7dmyBY5b0PdKMgvo8rVu3Du3Xvn1767///W+xbcHOrvn666/HfF1fffWV1aFDB/OaDj74YPP+PPXUUzGza8aq+1jtx87O+MADD1hO0c+dliHWFp3ddsWKFdb5559v2vN+++1nnXnmmdbChQtjHlezbur7ru+DZo684YYbTHbZcHq9sEy8JXm8fQy9T/fRffUxsbJoKv3+0DJr2fU16GvR11TSeiosu2Zh79+9995r2rR+9nUh9SeffDLm57o0baOwdvrJJ5+YdqefaV3gXduoXrf30wXa9TOs9aOvXTPA6udJyxOe8RWAe/n0H6cDTQBA4tAFkrUHRBNTRA+t1V4CnZuna45pb2y8aG+n9ubai4Qj+eiQbu3x0h6n4npTASDZMVwTAFAoHWqnwyR1qKLOcdIgS7Pt6by58ABPhzlqhlTdX4eCxjPAA5QOM9ZhyAR4AFA8gjwAQKF0Hp+uYadz+TQroc6b0/XsdJ5fOJ3fp4kidK5bcVlKgbLQdggAKBmGawIAAACAh7AYOgAAAAB4CEEeAAAAAHgIQR4AAAAAeEjSJ17RRXL/+OMPs9CnLrIKAAAAAG6kq99t377dZL5OSSm8vy7pgzwN8Bo3blypbw4AAAAAlNVvv/0mjRo1KvT+pA3yJkyYYDZNCW5XVM2aNcUtEfrWrVslIyOD3kXq3JNo49R5MqCdU99eRxunvr3OcuE5+bZt20wHlY5CLErSL6GgFaVvnL6BBHnJy40fYi+jvqnzZEA7p769jjZOfXud5cLzw5LGLiReAQAAAAAPIcgDAAAAAA8hyAMAAAAAD0naxCsAAACIv0AgILm5uVRtBcwPy8nJkd27d7tmfpjXWQ7UeVpamvj9/nIfhyAPAAAAcTkh/vPPP+Wvv/6iNitwfedNmzZRvx6v8wMOOEAOOuigcgWWBHkAAAAoNzvAO/DAA2W//fajt6kCgmjtJdVeHnryvFnnlmXJ33//LevXrzfXGzRoUOZjEeQBAACgXPRE2A7w6tSpQ21WAIK85KjzatWqmb8a6OnnqaxDN0m8AgAAgHKx5+BpDx6A8rE/R+WZ25q0Qd6ECROkZcuW0q5dO6eLAgAA4AkMIwTc8TlK2uGaQ4YMMZu9ajwA79uZVV/SJM9cHha4Xt63TjKXb0l5Ua5MeVdyJVWqZ69zuJQAAADlk7Q9eQCSjwZ46b78IC+a3l7dt1tysuqENg0KAQDwstNOO01efvllp4uRNP73v//JcccdZ7J2ViSCPACe1Wrku9I8c6asyGppgjYN5HKsVEnP3iTj7xoly+/uarYBd74kO62q5j6b7vuedaKj5QcAVF5m0Ouvv14OPfRQqVKlijRu3Fh69OghH3zwQWifQw45xAyj003nTLVu3VqeeOKJ0P2jRo2SY489tsCxNSGNPubjjz8u9Pk7duxo9rn33nsL3NetWzdznx5/1apVkpKSEipHrE33K03Aoa+9d+/e5vrmzZtNPbRo0cK8xiZNmsgNN9wgW7dujXjceeedZ+6rWrWqyQDZp08f+eOPP0L365ID55xzjjRs2DBUn9ddd50ZQVeUyZMnm7qoWbOmeS2xluPYsmWLeT4diaebXi5s2Y5NmzZJo0aNCj1WtC+++ELOOOMMqV69ulnGoFOnTrJr167Q/V9//bWcffbZ5j5NMPR///d/smPHjlI9d/fu3c1tFR1YE+QB8KzcgCU5gX2/lGkQp0MyY9Fhmhr82Vvz3FdleGBwJZYWAOAEDZxOOOEE+fDDD+X++++XxYsXy7vvvmtO8HVqT7jRo0fL2rVr5bvvvpPzzz9fBg0aJFOnTo1LOTQQevbZZyNu08BJy2Wn0td99DYtg27/+te/pFWrVqHruv373/8u8XM++uijctVVV5nA0X4+3R588EFTD1OmTDF10b9//4jHad289tpr8uOPP8obb7whP//8s1x88cWh+/V4PXv2lLfffluWL19ujjN79mxTX0XR5QM0OLztttsK3eef//ynLFq0yJRLN72sgV4s/fv3l2OOOaZEdaEBnj53586d5auvvpL58+eb9z+8bs466yw5/PDD5csvvzTPvWTJEunXr1+pn1vr/LHHHpMKZSW5rVu3WloN+tctgsGgtWXLFvMX1LkXVVYbb3nbW2YriyNum2FdOvwBa8/I2mbbMfJAK5HxvUKdex1t3Nk637Vrl7V06VLzN9F07drVOvjgg60dO3YUuE9fn61p06bWww8/HHH/EUccYfXu3dtczsrKstq0aRPzGHqu+dFHHxVahtNPP9269tprrTp16lifffZZ6Pa7777b6tGjhznuyJEjrdzc3Ij/Owt7zpLYsGGD5fP5rO+//77I/V577TUrPT3dPHdh3nrrLXOsnJycQvd55JFHrEaNGpWobFpXWmfh9a+0jent8+bNC932xRdfmNt++OGHiH0nTpxo6vWDDz6IeaxoJ510knX77bdH3KZ1bdf5E088YR144IFWIBAI3f/NN9+YY//000+leu5Vq1aZ23/++eeYZSnq81TS2IWePACe9U3qVWYrizS/T1L9vtDQTTthCwDAO3R4ovbIaI+NDtGLpsPyiqLDFcuT5j5cenq6XH755RG9edoDdvXVV0tF+Oyzz8yQzKOOOqrI/XSopg6fTE1NLbQOX3rpJenQoYOkpaXF3Ed7wd588005/fTTy1Vm7W3TIZonnZSfOE2dfPLJ5ra5c+eGblu6dKnpdX3++edDPXFF0TXptHdO16XT11G/fn1TVq0j2549e8x7FH48e0278P1K8txNmzY1z/Xpp59KRUna7JoAvJ09U9lz8MpiyehzRES3f5n5fACAss2N1qHzlUl/pMv/Di/eihUrzILXRx55ZKmeIy8vT1588UUzpPHaa6+VeNEhfv/4xz/kkUcekYULF5oA69xzzy3VPLvSDFPVYKaoIEjnld15550ycODAAvfdeuutMn78eDPEUgMtnd8X7bLLLpO33nrLzGvTOY5PPfVUucqs8wc1OIqmt+l9djB22WWXyQMPPGDmDf7yyy/FHtfeR+tZh6rq3EoN0nR4pg4H1fahc/WGDRtmjnvjjTfKzp07Q8NKdZhsaZ/74IMPNu9BRaEnD4Bns2cWNQcPAAAN8EqzLpkGNvvvv7/pwdHev5tvvjlmAFRWOofriCOOkP/85z/yzDPPmLlmhfWOlZcGXtoTWRhNkqIBpq4rnZWVVeB+fe3ffPONvP/+++L3++XKK68M1aft4YcfNslKpk+fbubtaZBUXrHeK31e+/YRI0aY3skrrriixMe0M13qe6nz5TT7pZZdE9Bob6rSuY/PPfecPPTQQ6YH9KCDDjKJejRQ1tdf2ufWNqQBckXh7AeAZ/TMy89KNvPua0K3pcfp2Bo8aqbOsvxSDADJyu3fkxpQaXCwbNkyk0ilOBrYaKINPcnXZCjhAYcOaYzOQqnszIolXZdZh2dOmDDBDPvTBCAVpW7duiZTZSzbt283SUg0oJ02bVrMQFMfr1vz5s1NYKNJYebNmyft27cP7aOBkG7aE6bZKE899VS54447QolkSkuPtW5dwfVsN2zYYIItpYlqFi9ebAJlZQeeWtbMzEzJzs4u8Hi7PBrQhtPXtXr16oikL7ppGXR4r77/Y8eOlWbNmpX6uXWYa7169aSiEOQB8IyfpVGFHFd7A3da+ePu1bzUAeIT/Y/99wp5PgBA5ahdu7Z06dLFBFW6VED0vDwN0MLn5enJumZXjEUDmd9//90MG9RgxKZZGnVIZGGPi6ZBhGbIbNOmTYGgI560t0rLqoFerVq1InrwtE506QPNjllUb5/NDmZ0uGJ59imOBpAaSGvwe+KJ+csc6Vw6vU3n0inN9rkrbNkDrX8NnHX+22GHHRbzuLo8hi73oNlCw2lmUK2LaHZAqb2tWj+6rEJpnnv37t2mZ1Pfg4pCkAcAxdDlFfS//eV7r+dk7abOAMAjJk6caAIEDRo0YYYOmdQ5d7NmzZJJkyaZXr6S0NT72vOja87dfffdJmjQpRY0YNOlA2rUqFGi42jApXO8KmqYpk0DDO1J+vzzz83abXYPnr4OHUaocw414LPXttN9dViiBli66dxBLavOOxs5cqQJYuxevBkzZpjernbt2pneQO2VvOWWW+SUU04xAVVhNOjUTedKKu0V03rT+W0akGv9ag/jNddcE1qjUNeq0/Lr0EoVHcht3LjR/NXHFpZIR3vktJdWh6VqcK1z8nRo5g8//CCvvvpqaD+dg6htRV+Ttg99jK5taB+3pM+tPZ4aRIf3esZbws/J04hb3wh70/GtOu4XQPJ5wn+f2QAAKCkdaqfzxnTtN113Thc5154ZXQhdg7yS0uyTOj9N52lplkydwzV8+HAZMGCAGdJXGhoQxMr2GU8asGkvk2bGtGmyF+0Z0+BKex51GKO9/fbbb2YfPdfWTJlnnnmmCaz0GFpnn3zyiQlc7H2efPJJEwhqgDN06FATiEUnZ9Hgyp7zph5//HETfGoQp0477TRzXXsUbVreo48+2gSjumlQ/sILL5T69fuinlvLqHPqbrrpJhPo6fuv72d44KbBrbYNfX5duF0DTe0BLq1XXnnFtBEd9ltRfLqOgniErjivvw78+uuvJf5g6K8TOkbaTg/rBvqWaHm0XCWdCAzqPJHEs42HZ9S0s2nqYuYVyc62WdHPE098r1DnXkcbd7bOdQjeypUrTcBUkuF9KFt9BwIBE5zF6/xQe9s0GNXgTtP6VybNLKlzIrWXT/+68bmtCqhznT+oQ3sXLFgQmssXTYdzFvZ5KmnskvA9eeE0ytdfFSr6lw8A7syoWZnZNM1z3llPZNW+tXEAAEgkOrfs6aefjkguUll0fUIdalnZAZ7Tz63Bmw4RLizA88ycvDlz5pi1JPQXBB1/rBl8orMbaUXoPnq//towbtw4k50n2muvvWbStwJILuG9d+mV8XySKjoGwpcXlKue+lLmW9vJtgkASEg9e/Z05Hl1nqJTBjn43Dr3004a4+kgTxcS1HGvuibFRRddVOD+qVOnmjGyGujpZE0d+9q1a1fTvaqTMMO7LnXiaPjkSADeNy54qfl7SyU+50nWcxGL+96TMnHvsAh3pwoHAADJwfEgTwM23QqjE1X79+9vJq0q7cV77733zETYMWPGhPZ76623TIpTxoEDyeWp4HmVHuRFr/uUk1XyBVcBAAA8H+QVJScnxwzj1MxE4TSTzty5cwsM1dSxtcXRicHh63PYaWF1YqVbctDYZXFLeZIBdZ649e2T/GM4+XmxzJp5zpahOLRx6tzraOPO1rn9/cf5S+XVPbxb51YRn6eSlsPVQZ6uLaEZbewFB216XdfQsNmLIuoChMXR3r9YK93rMdzygdFyaKZQRXZN6tyL4tnGr9r/C/N369b8RVCdkFstf+h42tat4lZ8r1DnXkcbd7bOc3NzJRgMmvM23VAxtI7h/ToPBALmeXXdwujF4+0OqoQO8mzRJ4H6pRJ+m6YR1RSwJaHrXwwbNiyioho3bmyO4aYlFBRLKFDnXhXPNn7T7vw1jNIzRotTcnatzs+2+VhzkRP/T6TzneI2fK9Q515HG3e2zvVEdNOmTSbVvG6oONSv9+vc7/dLSkqKWQg+eipaSc+bXB3k1a1b17zI8F47tX79+gK9eyWlizTqNmHCBLPZvzZphbmp18wuj5vK5HXUeWLUd/i6eKrK3rXxnPys5Ipfz3bElxeQ5z9bKffPedeV2TZp49S519HGna9zzl0qTviIM84PvV3nviI+TyUth6vXyUtPT5cTTjhBZs2aFXG7Xu/QoXxDs4YMGWIydM6fP7+cpQRQmZZah4TWxavstfGKyrZ5dN6L0jrvRbk/eIWcYc2TTsF5jpYJAAAkL8d78nQs94oVKyIWCFy0aJHUrl3bLJGgQyv79Okjbdu2lfbt28vkyZPNgo1Orm8BwDmXB0aJBESW370vK29lrI1XtmyboxwpDwAg+Xz44YcyePBg04mhQ/1Q8dq1a2emgl144YXiNo63gAULFshxxx1nNqVBnV4eOXKkud6rVy+zbMLo0aPl2GOPNYunz5gxQ5o2bVqu59Whmi1btjRvDgAAAJKXTg26/vrr5dBDDzXTejRfQ48ePeSDDz6IGCY3ffr0Ao/V9Zw7duwYut6vX7/QMLu0tDQzxejss8+WZ555ptAkHpo5XqcozZtX9lEgt9xyi2RmZoYCvDfffNM8b7169UzeCe0s0WXIwuk+2pFywAEHSPXq1c259gsvvBCxjy5bdswxx5hj2MeZOXNmseX55JNPzIg8nVOm9fr4449H3K/JevT8/rDDDjP76LrZ7777bsQ+o0aNihgSrNtBBx1U5POG178vbGvVqlWpXvchhxxi6lLfQ/1rH0dHA9ruuOMOswqAGxPiOB7k6YciOgWvblOmTAnto79KrFq1ykzq1SUVTjvttHI/L8M1gcRUXzaZDQCAeNBzTA1GtCfs/vvvl8WLF5tgo1OnThEn9KVxzjnnyNq1a82xNSDSY914443SvXt3ycvbN+VA6Qi1L774Qq677jp5+umny/R8urTYTz/9JJdccknoNu0Y0SBPO0f0/FnLoIHrN998E9pHR85pYKjP/91338lVV11ltvBgsFGjRnLvvfeajhndzjjjDOnZs6csWbKk0PLoyLxu3brJqaeeap7vtttukxtuuCEiE/7tt98uTzzxhDz22GOm91FH6V1wwQUR5VManGld2pu+P0V55JFHIvb/7bffzOsMr5uSvG6d0vXHH3+Yx+tfe/pY+HHOPfdck6E/Onh2BSvJbd26VWdUmr9uEQwGrS1btpi/oM69qDxtfM/I2mZzMzeWke8V6tzraOPO1vmuXbuspUuXmr+JpmvXrtbBBx9s7dixo8B9+vpser44bdq0AvvceOON1umnnx663rdvX6tnz54F9vvggw/MMZ588smI20eNGmX17t3bWrZsmVWjRo2Y5VBaz7m5uTH/77z++uutiy++uNjX2rJlSys7O7vIfY477jjr9ttvL3KfWrVqWU899VSh999yyy3WkUceGXHbwIEDrZNPPjl0vUGDBtb48eMj9tF6u/zyy0PXs7KyrDZt2ljlMW3aNMvn81mrVq0q9esOr3N9nw877LAC9d+vXz+rT58+VjwV9XkqaezieE8eAHiRJodZkdVSmmfOlFYjI4efAADcYfPmzabXTnvsdNheNB3OFy/aA6ZDEnWooE1jx2effVauuOIKOfLII6V58+by2muvlfrY2munww+LYq+7pr1YsWhZdHjqjz/+WOioOc1K/+qrr8rOnTvNsM3CaA+ZDkEN16VLF9MTqMM0lY7Qi14eoFq1avLZZ59F3KY9lA0bNpRmzZpJ79695ZdffpHSePrpp+Wss84qdKpXSV53Tk6OvPjii3L11VcXyG554oknyqeffipu43jiFadEL6EAAPFisn3uzbqcE9Bx+vyeBiCJ3Vmv8PsGzBZp0Cb/8quXi/z0fuz9TrtZ5PRb8i9//YLIO/vWPI5w0NEi13xY4qJp8j89ydcAqzLo8+jwQNvs2bPl77//NgGQ0mBPgxIdOlgaOixUA6GiPPTQQyY4u/TSSyNu1+GGBx98sAm6dF7gxIkTzTDPcDpEUoO63bt3y/777y/Tpk0zuS2KmuMYvdyZXtehqhs3bpQGDRqY1zx27FgTWOm8PA203nrrrYhz85NOOkmef/55E/zqmth33XWXybCvQ0Xr1KlTbL2sXbvWDJd9+eWXC9xXktdt07mYf/31l5nvF02PoUNuNYh2U8Ib95SkkjEnD0BFqZ69TtKzN8nh2Uuluj8gaZL/qyUAwJ3roFXWGmj6fOHPpQGdJhlMTc3vd7nsssvkyy+/NL1KpbFr164CvWLhXnnlFZPEZOrUqXLggQdG3KcLbmtme52Ddvfdd5skiB9//HHEPi1atDD7aGKYa6+9Vvr27Wvm0RUluk6j61rnzh1xxBEm8NVl03ROoga34QuPd+3aVS666CI5+uijTW/cO++8Y25/7rnnSlQvU6ZMMb2x559/foH7SvK6bZo0R8sSK5DW3kcN8DRYdJOk7ckDkBh0qGNuYN9ipItTRZxb9rz0vkm1f40lWQyAJHXHhpLt1/ulku13fJ/8LQ40yNCgY9myZTEDgeigQHt/omkPT0ZGRomeT59Hhx3aQ0W1h0iHL2oGS5v2ZGlQcd9995X4ddStW1e2bNkS8z4N7Pr37y+vv/66CZSiae/T4Ycfbi5rlkkt45gxYyIyhmoQZu+jw0I1MNIgTROnxKIZMLU3L9z69etNMGv3wGnWT3392ju4adMmE0Bppkq7fmLRIbUa8OkQzuJYlmXqUZdi0/KX5XWrX3/91fS4hg+zDafv43777WeCPTdJ2p48AIlBA7z8IY8itWSbVAlbCB0AgPLQ+Wk6bFCn8OhQxlgBnE17nDS4iQ4kNHOl9nQVR7N36rBH7ZlSL730kslc+e2335oeJXvTpcO0pyo6C2dRdPmxWD1r2oOnQwx1uKJmgiwJfU3F9UoVt48O7bSzUdref/99EyDqkgThtAdShzzq69Xsm5q5szD6nBqM6XDPkizhsGLFChPglkRhr0nfC+39LKz+vv/+ezn++OPFbZK2J485eUBieN2fKeIXaTP6a5GdG0XGpku63+nlzwEAXqFzsXSelybQ0HXbdE04DTg0SNEeNg0q1L///W8zTFGDPU0qokMkJ0+eLD///HOBpRY0WNCeLO2V07lkmtxFe4l0CYUrr7wyNFTz4osvltatW0c8VhOE3HrrrWZoYlEBTzgNVKOHMGqAp8+lPW4nn3xyqGdNe5zsnkctkwZeOidOk4vocgs6By68Z1GXP9Chirp2oCZu0cQrOqwxek27cLocwvjx480QyGuuucYkYtHXq2Wy6bDUNWvWmF40/avDSXXYo673Z9M612UfmjRpYnoCdU7etm3bzPtQnKefftrM6Yuu35K+bqXl0XrVerSH1EbTpCvRSWZcwUpyLKEARept99a3G5cjKA23lJ82Tp17HW3c2TpP5CUU1B9//GENGTLEatq0qZWenm6WVDjvvPOsjz76KGK/V1991Wrbtq1Vs2ZN68ADD7S6dOliLViwIGIfXUJBT7F1S01NterVq2edddZZ1jPPPGMFAgGzjz5G7//qq69ilqdHjx5mK+kSCps3b7aqVatm/fDDD6HbdFkHuxzhm5bPlpmZaR1++OFW1apVzbII7du3N68x3NVXXx2qF30tZ555pvX+++9H7KNLHeg+4T7++GOzLIE+7pBDDrEmTZpU4P6jjjrKqlKlilWnTh2zDMGaNWsi9unVq5dZaiEtLc1q2LChdeGFF1pLliwp9rn/+usvUx+TJ0+OWb8led3q3XffNXUWXq/hfv/9d1O23377zXLbEgo+/UeSmP4aoL9m6BjrmjVrihvoW6Ll0XJV1kTgZEedu7e+c7Lyx+5rIpNEpOXX5RRyrFTpmXev/CyNJM3vkyWjz6nUctDGKx91Tn0nUxvXnitdAFvnUxWVAATlq2/tGdTEJLH+79QeMH0/CpsnV5HsrJOa6MRLz20VU+c333yzqXPt0Y0nnadY2OeppLFL0g7XBIDKXk5BsaQCAKAiZGZmhpYHC89QWRl0/puu1ecEJ5/7wAMPNENK3YggDwAqeDkF20wR+WiknbWrK/UOAIgb7d3R+XNO0F4npzj53DfffLO4VdIGeSReAeCEU3yLqXgAAFChknYJBRZDBwAAAOBFSduTByAxZAXy17cZ43RBAAAAEkTS9uQBSAxvWJ3M5iWabVPurCfy2TiniwIAcaXrigFw/nNETx4AV2k18l3JDVgR2SjT/d75PSpHUkUXrvHlBWXcez/KUzNnOrKkAgDEU3p6uqSkpMgff/wh9erVM9dZBiq+ikvnj8Svc8uyzOLsGzZsMJ8n/RyVFUEeAFf50tdX0lLzxP4q9acGZKJ1iWeyUZ5kPRcRxF4q74nf0ldLkAcgcekJqa7ptXbtWhPooeJ6eLSu4e0632+//aRJkyblet6kDfLIrgm4U7rk5Q9n9Nu/XqXI9f4Z4hXRPXY5WVfsvfSwI+UBgHjRXgc9Mc3LyzO9H4h/L8/27dulRo0a9OR5uM79fr+kpqaW+/lSkzm7pm72qvEA3CPHSpX0OzY4XQwAQCnpiWlaWprZEP+AY8+ePVK1alWCvEpiJXCdJ22QB8CdPraOM387O10QAACABEWQB8BVrgv8y/xd7nRBAAAAEhQzNwEAAADAQ+jJA+AqR8lKSTaaaKZ55kxzmeUUAABAeRHkAXCV11Pv2HtpsCSDXEmVb4OHm8u1AxskjQEWAACgnJI2yGMJBQBuUD17nbTbOwcxJ6vO3lvtZRUAAABKL2nn5OnyCUuXLpX58+c7XRQAAAAAiJukDfIAAAAAwIsI8gAAAADAQ5J2Th4A92g18l3JDVjm8uJUEZ/TBQIAAEhg9OQBcJwGeDmBoLmc7AGeLqcgOzc6XQwAAJDA6MkD4LgPUm8w30aNs38SkU2SrHIkVSxL5LR7PpAtUlNe92eKz+eTY7IXOl00AACQQAjyADiuvmx2ugiucJL1XGjYqjrS92vS92wCAIDSI8gDAJdYMvqciOs5WY4VBQAAJDDm5AEAAACAhxDkAQAAAICHJG2QN2HCBGnZsqW0a9fO6aIAAAAAQNwkbZA3ZMgQWbp0qcyfP9/pogAAAABA3JB4BYDjrgvcZP5OdrogLpMV6C95liVvZc6UwSlvyOCUaZIrqVI9e53TRQMAAC5GkAfAcR9bJzhdBFd6O+WMiCUVzELp+64CAADERJAHAAmxpEJXycma5mBpAABAoiDIA+C4zJQpey91dbgkAAAAiY8gD0ClOzrrXalXJSh//O0TS3yyOHW2+HgfAAAA4oIgD0Cl+8J3lezyN5SM1NWiYZ7ONcux+DoCAACIh6RdQgGAc/4X/Iek+wKSnppiNvGnS3qVqrwlxfjYOs5sAAAAReGncwCVblTwapmcZ8mHt3cT8TFQs6SuC/zL/F2+9luRp87Kv9GfLnLbmop6qwAAQAKiJw8AEkhOICg9x38uOXlBkUCO5OzZ7XSRAACAyxDkAah0HXzfSRvfT9R8KaX5fZLuT5Fl0kxa570oe5jHCAAAvDpcc+XKlXL11VfLunXrxO/3y7x586R69epOFwtAIR73PyjbUpqIyI3UUZnXzRPJyaL6AACAR4O8fv36yV133SWnnnqqbN68WapUqeJ0kQAAAADAEQkf5C1ZskTS0tJMgKdq167tdJEAAAAAIHnn5M2ZM0d69OghDRs2FJ/PJ9OnTy+wz8SJE6VZs2ZStWpVOeGEE+TTTz8N3ffTTz/J/vvvL+edd54cf/zxcs8991TyKwAAZ+i8PN0AAABcFeTt3LlT2rRpI+PHj495/9SpU2Xo0KGSmZkp33zzjemx69q1q6xevdrcn5uba4K+CRMmyBdffCGzZs0yGwAAAAAkI8eHa2rAplthxo4dK/3795cBAwaY6+PGjZP33ntPJk2aJGPGjJFGjRpJu3btpHHjxub+bt26yaJFi+Tss8+Oebw9e/aYzbZt2zbz17Iss7mBXRa3lCcZUOeVXN/iy99o4+Witfi+/0bZk7XFXM+VVKk+am3sOud7pdJR59S319HGqW+vs1z4f2dJy+J4kFeUnJwcWbhwoQwfPjzi9s6dO8vcuXPNZQ3wNKvmli1bJCMjwwz/HDhwYKHH1MAwOzu7wO1bt251zRuo5dixY4e5rENYQZ17zaZqR8quqgdK2tattPFyaFJDZLccLNukhqT5ApJr+SVv69aY+/K9Uvmoc+rb62jj1LfXWS48J7c7qBI6yNu4caMEAgGpX79+xO16/c8//zSXU1NTzTy80047zbwRGgB279690GOOGDFChg0bFlFR2guoAWLNmjXFDexgU8vklgblddR5xTo6613JDez7ESU3cIscUsMns2nj5fL+rToKIn8kRM6oelJNRNIzMmLuSxuvfNQ59e11tHHq2+ssF56Tl7Qcrg7yCnsxWuHhtxU35DOcLq8Qa4kFPZ5b3rzw8ripTF5HnVeceb5+kpaaJ9qacyRVTpBnxe933+cu0Ydumr9F1CdtvPJR59S319HGqW+v87nsnNwTQV7dunXN4uZ2r51t/fr1BXr3SksTteimPYUAKla65Em6L0/Eny7p/lT5YVRXM0QaAAAAHsyuWZT09HSzZEJ0tky93qFDh3Ide8iQIbJ06VKZP39+OUsJoCRyrFSROzaI3LaGCqsA1wVuMhsAAIDjPXk6mXHFihWh6ytXrjTZMXVR8yZNmpj5c3369JG2bdtK+/btZfLkyWb5hEGDBjlabgBwk4+tEyQnEJTn7ugll6XM3pdtM3ud00UDAADJFuQtWLBAOnXqFLpuJ0Xp27evTJkyRXr16iWbNm2S0aNHy9q1a6V169YyY8YMadq0abmel+GaALwkze+LGJxhhse6I2EwAACoZD7LLesGOESza2rGHJ0f5KbsmloeN2Xy8TrqvGJdcvuj5u/rd91AfVeSnKw65m969ibq3CF8r1DfXkcbp769znLhOXlJYxfHe/IAeN+31hFOFwEAACBpEOQBAIC42ZlVX9Ikz1z+n9VBhgcGm8vtfEvlWf895jLzRQGgYiVtkMecPKDyjPXnD9e0F+5GxXsj2FGCliXZmTPlFN+38rj/QdlVrZlk3PoV1Y8KpQGezgk1GXVj0Pt0WRW5s17+Df50su4CQJwlbZCnSyjoZo9rBVBxOvsILCrbXb5rJDdoRZxYp4p71wVtNfJdyQ3sK+/Xqf3yA4Ew862j5KpAprncRP6Ut1KHkz3URe/fTdYLcmXKu6EAT+eDXihitn0/8vxLduzt6fPlBc2+W/KqSC1HSw8A3pO0QR4AeNmS0eeEXesqe7IeFDfTAE+XgEj3p8iVKTNkf99uc3thvUGDU/4j1X27QwlmGP7n/PsX8Fmh90zfj/RC9j3Jei4U0Nvv+XK9cl8zkZwd9OwBQBwQ5AGo8J6Zxaki7shJBbfO3Zqfmi6PpVwst9/5mMi8VSKz0s3Jfvpta0L7nyKSHwyYx/YLBYAsF+Hse6fOsK6Th6SPDMx+Of89KeEPEM0zZ5pAT/9+mZortXw5kqM9fBVaegDwvqQN8piTB1QsDfCu970qA/3/K3J+DipPmi8g60YdKqflTTTr6kX29jk7d6ua5MhNKa+JyGMiJ1+bvxUhfJF37c3T4xyX+YbslP3yj+2C15dM8+5SU3ySluIr1/qOJ+U9KYtTr+AHIQCIg6Q962JOHlDx/D6fpKfqCVy6pGtyBThGh8/lWn5zWXtOwhdOd6J391n/3XJKyr65Wyq9HK8tfOH36f5bJX85I4K8ihT+3o0v4zGiA/GcrPyeWQ3crwrcJvOtlgTsAFAGSRvkAah444K9ZPAdU6hqF6g+aq3kbd0qdTMyJP32d01ik5ysfcPtKmJOW/SQvjutk+WR4IWyNqVhieZulZRd7m/2Xs/JuiIUKNiYs1f+9/LEvCclV9LM9eVpFdM7Hx2wO/WDBAAkOoI8ABWiu++zvZdYNsFt8ofI7VNRc9qih/R1T5knXdIX7A3K8ttFRfTvRgcKzNmLz3s5xPe6jLP+aa6vsA6WBrIp7u9f+A8NL4nImyN77L3G9wgAlAZBHoAKca//8b2X7qaGXSZ/iFzBOW2X3P6ofGsdUa45bfaQzOryt3yzt7fHHtKnKmPQbnSP5LDb85ddGFsJz+1lA/wzZHD2C5UadHX3zc0P0u019RTr6gFAsZI2yCPxCgDE7vl6MOUROS/lC8nJSi12qGP0kMzX5WDpFhgj1f3xG45ZXv+z/mGG/R1+Rz8ZkPJfhm4mkBxJFcsSs6aeMj3DZN8EgGIlbZBH4hUAyGcHcK8XsjyBLkquvX1PBXuYeZb2cFztra2uvSxh69lpwpPq6X75ZvRFInKR4wFedAZHhm6WXHgA71SG3PA19dS/fC+YhE4DKr0kAJBYkjbIAwDEFt5rF91TF83urbMfc5QOB3VZxe4bdtpVcrL+63BpEsfFeXfJkJTXpXPKQsd6ZaOHDDfXkbeWEOQBQDEI8gDEPT2+YgF07wV8g/du++Zk5c+3dENvHeLvZ2kkw4I3yfI7Ky5JDgCgYhDkAYiLL3195eTAeMnx72+uV2EBdLiUDj3838iuMixwg7nOwumJo7PvS7ku5T9mmQzF0hgAEBtBHoC40Hlb/0p7U66685X8GyYdLelbVlK7cHWSmVbBHyXVp/P1WDg9emiuroX3UfA4Vy1fcJ9/UsQ8wYpY+gMAvCBpgzyyawLxd3nK+/uuXGuvkwe4b/hp972b3SMkkt+rl8yi1zXUvyf6lolbhw+HL3YPAIiUtEEe2TUBAMoENmEBQzIPAXRiXUMAQPwlbZAHAED08E2nlgpwg3HBS83fWyQxrJb6Zg29czNnSprkytepV8uuas0k49avnC4aADguOf8nAwAgavifap450/xdnmRz8O4OXinP5XaXdH9KwgR5Pa2xERl9NUBPlYCjZQIAtyDIA1DuJROayJ8yOy15e0DgLTmBoKzLOkRqyXZPD9+MnoOnAV7+wvGJIXoNvT1ZfP8AgI1vRABlMlnukmCKJf2Dt5vrTi2WDMRTfpCj2TYlKTI42nPwskXMBgDwBoI8AGXSTrPu+SS0ULLIVQR4SHj7eodWmX81IYudmGVg4Gb5KuXYAj1IcI80X0DkrvoiA+eIHHik08UBAMckbZDHEgoAgNIkZskLWpKrmT4SeA7eOXkPyTqpE1oHz0tDrHMkVXIsv+TkBeX8Rz6Vn2UlC90DSFre+XYvJZZQAEo/B+8y6x25JeVlSfYshEge4XPxvsycKVm+J8PW1kus+Xo6B2+I73UZaQ0y1+cHj5SWvlWe6YFvbz0r9QJB+SPPJ5b4ZFLKfXsH3tLzCiD5cIYGoEQ0yUpAuzT2TldiDh6Scb5eirUvMYn+0JEu+UM5XwmeJXcH+7mu5yg8QdLiVJEL/Z9J7+y39t5rD7X2hsXZ58jWrVslIyNDfD5fRDAOAMmGIA9Aib0k50pW9vjQda/0AAAlkR+8nRNzCQI7K+c/rK8lJ+tyc117k3KrNZGcXatF+5ac6PX70tdX0lLzdPosve8AkEQI8gAUKvwkdn5qujwWvNhzv/4DZRUesPUVkftHviupwdhLEGiv3zqrtlSvhM/sCutguSRwt7l+u+9UuTJ1lohff5JJl3TzN3mY7Kh31su/oq/9tjVOFwkAKgVBHoASraNVTXLkppTXROQxagwosqcvf1kRy7IkbetWSc/IkOa3v2t6+tIzZ8r3qQWHEZa1ly98OKZanJonbVJ+FntN8AnWJXKR9alUvyMx5g3GOxGL5snx5QUlICkyMfc8+ZfThQKASkKQB6BE62ip5OoDACpu/b1w9ty+npkTZZk0M7eN9z8kHX3fxNx/YvACmRi8yFzuKR/KXanPmOGY9rH0M7v87vAe915J+VaeZD0XCoBNgO1PIcgDkDQI8gBEDPXqmXevrJaDPJliHXBKZDKW/B9NCpvbVxqpPl/+fLtUO4BMviGZJanz5pkzHS0LAFQ2zt4AhOiJ5uCU/8jw4HXm+ufBo+VY30/04AEVyB6maee8zFf43Nehe7d9+91fgaXzhstS3t97iTnFAJIDQR6QxArO5xHpnjJPLrzznb23cEIEIPFlpjy/99LDDpcEACoHQR6QxMLTqysWOAfgVeb7LatOQi5kDwCllbRB3oQJE8wWCOxNQQYkYe9demC8/CvtTbkqffbee5nPA8B7NKATKzLg+zZ4uLRzslAAUIGSNsgbMmSI2bZt2yYZGRlOFweo9N67X+Rg6Sn3yYO+vnLVHa/wDgDwrOgeOzsRy3KHygMAFS1pgzwgWWmqdv0V+6j6NWX5tcy5AwAA8BqCPCAJhmd2CX4q9/ofj5h3l37tZ04XDQAcUT9qGQsA8BqCPMBjwtfceirYQ3bmXiJ5KZaIP/9+DfB0fgoraQFIVh+k3mh+8JI7b9p3o64veNsaJ4sFAHFDkAd4jAZ44Vky0/0p8oH/NEnPvie0DwEegGSWI6liWSK+vKC5br4z84J8NwLwDII8wCNDMu+TR6Sz76t9wzGzN8lgEbMBAPY5yXouao3QK0JLyQCAFxDkAQk+JPPyQJbszDtMgmn6szTDMQGgOEtGnxNx/aTMiebvl1QdAI8gyAMSpKdOlz7QzJiqus4l2Tu/zh6SOcI3VM7Lzj9xYTgmAJTcFqlJdQHwFII8wMU9dVukhpyaN0lyAkGRKvvut5On6NpPrztZUAAAALgOQR4qpNdpvNwrp/gWx7x/XPBSeSp4nrl8Wcr7ckvKywUWqk2WerLnhJzi+1ae8D8Qus/uqZO9U0a0p07nkIQPMaK3DgDi43V/prRJ+VnkzrBvVrJtAkhgBHmI6/ywnnn3ys5AQwmmiZkfVpzrUt6Q6r7dkpNVJ3Sb3UPl5XpSS1Ly5KXgmZJtXVNgX7unrn72KlleyeUEgGRzmG+N7LFSybYJwDM8EeSlpqZK69atzeW2bdvKU0895XSRkqYHSi1OLZiy/wbfcFmyd35YtFv2bmpn1p7Q48xjfXlm3lnzzJnmeprfV2CCvJvr4uvUfqF5c9cFbpKPrRPM5cyUKdLXv9tctl+v/j3f/7lcnv2miHQVkeERx6anDgAqB9k2AXiNJ4K8Aw44QBYtWuR0MZKGCWoCe+T7Kv3N9fCU/fmhWclF99hpb9dS6xBzuXZggwz1vSE5WZfHfOypeRNDk+V1qM1Rvl9j7pcV6C9vWJ3yA8YzvxOZ84BZDymcJT6ZV+1suXF7H3P5KFkpr6feUWi5W+e9aP7qXLk56TfKQb4t5rpZXDcskAsXPo/ORiAHAM6L/jFxxO35/7+Ncag8AFBengjyUHlDDVdLfckJPCDV/SmSnpqy9950Sde5C3GgAVA7ETNEMXp4Y3loMJYTEBk3e7kMTgnaU93KvS6S9lr6fL6IutB5HOm3rZHJEXtqTx1BHQAkAv1hUBHkAUhUjgd5c+bMkQceeEAWLlwoa9eulWnTpsn5558fsc/EiRPNPnp/q1atZNy4cXLqqaeG7t+2bZuccMIJUq1aNbn77rvl9NNPd+CVeJsGW6aXysoPbMTvF7ljQ4U+Z3Hz8r6MEUTFov9Jv713WOXE4EVmix4KalmWnLx1q/yQkWGCtnyFLyMeOU/up2JfCwAAAJA0Qd7OnTulTZs2ctVVV8lFF+WffIebOnWqDB061AR6p5xyijzxxBPStWtXWbp0qTRp0sTss2rVKmnYsKF8//33cu6558rixYulZk3WvCnvXLNOwXky1v9YxJDMw7OXJmQiEDfP6wMAuMvglDfk+pT/iNypWcTItAkg8Tge5GnAplthxo4dK/3795cBAwaY69qL995778mkSZNkzJj8gRQa4ClNvtKyZUtZvny5ScASy549e8wW3gto9+To5gZ2WZwsT14gKHlaBn9+r9YeK83MJ0tzSR15sc6TCfVNnScD2nni1nf/lP+Z//Py8qz8Hzn1L/8/VGido3jUd+WzXNjGS1oWx4O8ouTk5JhhnMOHR2Yd7Ny5s8ydO9dc3rJli+y3335SpUoV+f33300P36GHHlroMTUwzM7OLnD71q1bXfMGajl27NhhLu8bOli5Gu5nyY9youwaclOBevIiN9R5MqG+qfNkQDtP3Pq+pNozkrc3c/JU/x1m/naaR///Kw/aOPXtdZYLzw/tDqqEDvI2btwogUBA6tevH3G7Xv/zzz/N5WXLlsnAgQMlJSU/AcYjjzwitWvXLvSYI0aMkGHDhkVUVOPGjSUjI8M1QzztYFPLVFkNaueoBhFJTj705cnTge6SkdFNkoETdZ7MqG/qPBnQzhO3vt+/dd8Io5xRfc3f9IyMcpbQe2jj1LfXWS48PyxpOVwd5BX2YrTC7ds6dOhg5uCVlPb46RbrOdzy5oWXp7LKlC65EWvd5Vp+uSxllqvqxGt1nuyob+o8GdDOE7++P7GONX87839DpdU5Ckd9Vz6fy9q4J4K8unXrit/vD/Xa2davX1+gd6+0JkyYYDbtKUQ+e607G2u4AQCS3XWBf5m/iZh0DEDyshf3cqX09HSzNMKsWbMibtfr2oNXHkOGDDHz9+bPn1/OUnrD/6wOZgMAAJGOCP4iOVl1Qpuu4woAbuZ4T55OZlyxYkXo+sqVK2XRokVmXp0ukaDz5/r06WOyZbZv314mT54sq1evlkGDBjlabq8sk6Brx6mcwCCz/t2FThcKAAAX0TVVU03qlXz2mrEA4GaOB3kLFiyQTp06ha7bSVH69u0rU6ZMkV69esmmTZtk9OjRZjF0XSZhxowZ0rRp03I9L8M1Rb709ZW01DzzX9dVvttkUUqrcr+fAAB4yb51Vgebf7UnDwDczvEgr2PHjsUuXTB48GCzxZMO19RNs2tqxpxklC55+b9I+tPlpatOEjnkH04XCQAAAECiB3lwQbKVOzbwNgAAAAAeQZAHAABQQq3zXtyXbXPcMSLb15oRMXLbGuoQgGskbZCXjHPywhOtqMWpEjaVHAAAlEROICjNM2fKB6m7pLEvR3Lygiw7BMBVXL2EQkVKxiUUNNHKvNQBoet/SzXZKdUcLRMAAImWbVOzUasz8x6VPVbS/l4OwMX4ZkqyRCv7+3bL8ru77r3ld4dLBABAombbzJeT5VhRAKBQSduTl8yJVgAAAAB4V9IGeTofr2XLltKuXTuniwIAABLYj1YTp4sAABGSNshLxjl5AAAgvnIkVR7Ou9AkYnnujl7yQlZvqhiA4xi7BwAAUEYnWc9Jri8/c3XvlNlkrQbgCgR5HsaSCQAAVF4iFpKwAHCLpA3ykmGdPF0yIS01T96xOsjwwGC5OnCbpKb45DmnCwYAAACgwqQm85w83bZt2yYZGRni1SUT0n15cuGxjeTCC3TZBHvpBAAAUCH/9/ryJCerjrwR7ChZwQFmXb3oZRcAoKIlbZCXTEsmpF8wyeliAADgebl6WpU/Pc/ICQSTOccdAAcR5AEAAMRB9ex1ocuXicj7t9+79xojaQBULoI8AACACvCE/4G9wzcfiOjtCw8GAaAiEOQBAABUwvBNDfjCr7vBzqz6kiZ5oev3B/8pzwe7mctXpsyQm1JeIygFkiXI+/rrryUtLU2OPvpoc/2tt96SZ599Vlq2bCmjRo2S9PR0cbtkyK75fDB/ovcApwsCAEASiu6xW5d1iPlbX5xfXulZ/93SzrdMqmvguXcOf7TrU/4j1X27TSIZGz2RgIeDvIEDB8rw4cNNkPfLL79I79695YILLpDXX39d/v77bxk3bpy4nVeza4avjZcT+Kek+1MI8gAAcIFT8yaZZCzpmTPl69R+ZuH0ih66Gd1TVzf3fvnVOkjEL6HgLjxwu33vlv/YYETwpz2Rmrk7ZNPPIs9fIXLTFxX6GgBUUpC3fPlyOfbYY81lDexOO+00efnll+Xzzz83AV8iBHlepQGe+Q/En2I2Td0MAACcl/9/cn62zf21hyxG71m8A7vonrq0lBSp7vfLKaM/C+1f2Pir6AB0R1Z9+SZ4uPTPnGmuP+CfIB32y6mQ1wCgfMr07WJZlgSDmhZYZPbs2dK9e3dzuXHjxrJx48ZyFgnlob8MpqfmSXpqikidw0U6DqdCAQBwgfD18nKyUkNr6pV3KKQd1P0sB0vPvPskPTBehqa+IVf6ZxXoqZtdjvKfZD0nucF9kwq7+ObJdmlSjiMCcFWQ17ZtW7nrrrvkrLPOkk8++UQmTcpfh23lypVSv76TI81hL4BufpfbtEJk+mCRlj2pGAAAXJyU5YtgSxmXd5EsyZwpY/2PSmffVxHBWfh0jO6+z+Re/+Ohx9q9dfbxcvz7y8Mp/WRA9qvmerwyJUQv6r4nSyTNF5CcUfXEJxbz9YBED/IefvhhueKKK2T69OmSmZkphx9+uLn9P//5j3To0CHeZURZFkC/YwP1BgCAS0X32A3WIC7FKjD/ze7pW5KSJ48Fz5fHrN4FjmX31h2V/Z0sl8qjz5lr+aWaSzOHAsnMZ+nYyzjZvXu3pKammi1R2IlXtm7dKjVr1hQ30LdEy6Pl8vlKN6fO/s8gPXtTBZXOm8pT56C+EwFtnDr3Oq+18eiEKW7MbBle58ff/qbkBoKmF/Gt1FvlMFnjuvImOq+18URgubDOSxq7lCkaO/TQQ2X+/PlSp86+ceR2kHf88cebjJtwxuq9iZnz+1YBAEAiihUcuXmBqhx/dckN68qjZw9wVpmCvFWrVsVcX27Pnj3y+++/SyLw6jp53fIeMn8rc7gGAABIbpHz9bqakUUm0LuzXv5NPSeKHHOJU8UDkk6pgry33347dPm9996LWF9Og6UPPvhAmjVrJonAq+vkAQAAOC1HUkUnBPny8rOxD39tkbz3nxoFkrcAcEGQd/7555u/Oia1b9++EfelpaXJIYccIg89lN+TBGekSS5VDwAAHGWWW9ibDVTlr+FLZhbAlUGevTae9tbpnLy6detWVLlQRt+kXrX3EolXAACAM6J77Cbe0U8G+6eZ9QEVSVkAF87J0/Xw4A7h6+aoxaki7sj9AwAAkK9vykyz1IMiKQtQ8cq81oHOv9Nt/fr1oR4+2zPPPBOPsqEEvvT1lbTUvFBgp1+c9pcoAACA27KF2klZ/jeyqwwL3GBuS/P7mK8HxFFKWR6UnZ0tnTt3NkHexo0bZcuWLREbKs8leXfJB1ZbSU9NMZv40yW9SlXeAgAA4Eo6VDP8B+lWwR+lZZC84EA8lanL5/HHH5cpU6ZInz594loYlN7P0kiGBoZKt9FdqT4AAJAwvXrd9245WVfsvSe/Vw+x7RzVQBZUO02GbO8nlvjkMPld3kodXmh1Me8xuZUpyMvJyZEOHTrEvzQAAABIOma6SVad0HUClHw7s+pLmuSZy9V9AWnr+7FE9flUoJuMy7tIJHOmzEj9lzSRddRpkilTkDdgwAB5+eWX5Y477oh/iVAqT/jv23uJnjwAAJB4NKCTsNUVkjm/QHhCvQEpb8st/t3mstaHbrusdPnhrq5mObN818Q8zoSR74qELVmhdZouBNLJpEyfoN27d8vkyZNl9uzZcswxx5g18sKNHTs2XuVDMU7xLaaOAACAJ5KyqOaZM/PX1cucKZ+mXiu1ZHvS9EJdbL0nAbHkNelirmtgZ792y7Kk2tatZVjComtEj6D6xWog44MXy/i4vwIkdJD33XffybHHHmsuf//99xH37ftlwd0mTJhgtkAg4HRRAAAAsJdm2gzPDejlJReig69sf57kpKTKXdnj9o7SmiTpcXie6AD5vMyZ+Rfev13kyyfyL/vTRW5bE4dnQ8IGeR999JEkuiFDhpht27ZtkpGR4XRxAAAAUKAXalVoyQV7zl4i9+ppUHdO3kOyTvJfy/K0fcMx7b/6+uIR2BVHe0uf+nSlXJkSzK/fvGClPC8qR3IOeAYAAEDCzdn7zjpUHsm7UL7MnCnZKU/JRSkfJ1TQp712Q3yvy0hrkLk+P3iktPStiih/eiX2lt4fvMJsi1OvCK25jCQO8jp16lTksMwPP/ywPGUCAAAAjPAAqK8mJvFZCTWUM3xIppb3Qv9n0jv7LUcT50X2luoyFvlle3NkDxkeGMzi9Mka5Nnz8Wy5ubmyaNEiMz+vb9++8SobAAAAUGhCkeihnK8Ez5L7fVcVCGKcypBZS7bJl2FDMitzOGZZe0t1GGf4nEgkUZD38MMPx7x91KhRsmPHjvKWCSX84lBXySXiT/HJLdQaAABI8uUXgpYlO3MD8n+33yXj/Q9XylDO6OQpS1LyZJEcJpcG7i6QIVO5LcBTdtkuFJE3bn9o760sz5XI4jon74orrpATTzxRHnzwwXgeFmE0wDNphf35v7A8Kz0lzUeQBwAAkk90APfgyHclfe/6cLHWhlsntaVx9k/7HnDPwfsySurlQI7IgNkiDdqU6Ef3i3wfyZjUgslTjvCtkeV320FSL1cGdoV51n/P3t7Re0K3JdK8R1RAkPfFF19I1apV43lIRPk6tZ+kp+ZJeurebvRzxoi0G0A9AQCApLdvmKauDfdYRA+bhn661pyuw6euT3ldrvfvCAWBZn6fiPQc/7kskz9kvP8h6ej7JmadXh04Xx4LXGjOpKN76syxPLY4vdvnPSJOQd6FF2pn7j76gVm7dq0sWLBA7rjjjrIcEiWkv0jlfwkl8tcHAABAxYrueSow5SXlnVDvW2jOnKTKMmlW7LF1qkx1v1/GjL5fRO731FlZdL1tyWpk/oZe433N8ns8WVPPe0Fe9LpyKSkp0qJFCxk9erR07tw5XmVDIfRLKP2ODdQPAABACRVMxlJw+KEGMstD1wqfkzZ075YMTsp7Mn+q0N4e0OVpm/PPRZ0uGOIf5D377LNleRgAAACABGKvqWfbY6Wypp7X5+QtXLhQli1bZtbMa9mypRx33HHxKxli+tY63PxtR/0AAADAgTX14NEgb/369dK7d2/5+OOP5YADDjBz8rZu3WoWSX/11VelXr16Utn+/vtvOeqoo+SSSy7xdHbPywOjzN99QwkAAAAAYJ8yrXR4/fXXy7Zt22TJkiWyefNm2bJli1kIXW+74YYbxAl33323nHTSSY48NwAAAJAsNAng5bc/ZDKVakIbeKQn791335XZs2ebnjObDtecMGGCI4lXfvrpJ/nhhx+kR48eJtj0svqyyekiAAAAIEmFL7GgCVnK2GeEClamdyUYDEpaWlqB2/U2va805syZY4Kzhg0bmrl906dPL7DPxIkTpVmzZmYNvhNOOEE+/fTTiPv//e9/y5gxYyQZfJB6o9kAAAAAJ5ZYSM/eJC/d9S8Zm/643OufyJvglSDvjDPOkBtvvFH++OOP0G1r1qyRm266Sc4888xSHWvnzp3Spk0bGT9+fMz7p06dKkOHDpXMzEz55ptv5NRTT5WuXbvK6tWrzf1vvfWWNG/e3GwAAAAAKkd331y5MGWOWVDe3nZm1af6E3W4pgZkPXv2lEMOOUQaN25seuA06Dr66KPlxRdfLNWxNGDTrTBjx46V/v37y4ABA8z1cePGyXvvvSeTJk0yvXfz5s0zyV5ef/112bFjh+Tm5krNmjVl5MiRMY+3Z88es9l0HqHS5DG6uYFdlljlsfYmrXVLWb2iqDoH9e0FtHHq3Oto49S517mxjedImlhW/rmpPVdPr7upjF6r85KWpUxBngZ2X3/9tcyaNcvMhdMn0zl5Z511lsRTTk6OWaZh+PDhEbfrvL+5c+eayxro2UM1p0yZYubkFRbg2ftnZ2cXuF2zg7rlDdRyaMCqNIAOl1utifmbtnWrI2XzqqLqHNS3F9DGqXOvo41T517nyjZ+0w+SF3b1lUm3m7/ne+Q81XJhndsdVHEN8j788EO57rrrTO+Z9padffbZZrODpFatWsnjjz9uhlTGw8aNGyUQCEj9+pHdvnr9zz//LNMxR4wYIcOGDYuoKA1aMzIyzGtyAzvY1DIdM+o9yQ3sCz6/TV1t+vLSMzIcLKH3hNe5Wz7EXkZ9U+fJgHZOfXsdbZz6jjZmRzfzt69HzlMtF54flrQcpQrydKjkNddcEzMY0hc/cOBAM7wyXkFeYS9GKzzWC+zXr1+xx6pSpYrZYj2HW9688PLkBDRzkSXp/vzpkyl70xm5qaxeYdc5dUt9exVtnDr3Oto4de51bm/jOq3oDGue5I7qE5GNU5O1JCqfy+q8pOUoVeKVb7/9Vs45J3LV++hhlDq8Ml7q1q0rfr+/QK+dLsYe3btXWrrcgw4xbdeunbjZ16n95K30O2T53V3Nln5iP0lPJVUtAAAA3CXN75PUlMg5emkRAzpRWUoVLaxbty7m0gm21NRU2bBhg8RLenq6WTJB5/6F0+sdOnQo17GHDBkiS5culfnz54ubpUuetEn5ed8NHUeI+NOdLBIAAABQwJLR58j4u0aZJRZ0y7HKlP4DcVCqmj/44INl8eLFcvjhh8e8/7vvvpMGDRqUqgA6mXHFihWh6ytXrpRFixZJ7dq1pUmTJmb+XJ8+faRt27bSvn17mTx5ssnkOWjQIEkW+gEJhXXV64rctsbZAgEAAADwRpDXrVs3k7lSlzzQhcnD7dq1S7KysqR79+6lKsCCBQukU6dOoet2UpS+ffuabJm9evWSTZs2yejRo2Xt2rXSunVrmTFjhjRt2lTKO1xTN03sAgAAACC+Vkt90dwl52bOlDTJla9Sr0no+XmJxGeVYt0AHa55/PHHm3lymmWzRYsWZvLfsmXLQgGTLq1Q3vlylUmza2rSGM0O6qbsmloeLVfuqLrmNu3yRuXUuVsm1noZ9U2dJwPaOfXtdbRx6rs4rUa+G8oSP9T3sgzwz0ioc1rLheeHJY1dStWTp8Gbrk937bXXmqUI7PhQX3SXLl1k4sSJCRXgAQAAAKi4OXq2nKwrqOZKVOrZkDpMUodLbtmyxcyl00DviCOOkFq1akkiYbgmAAAAAC8qc8obDercvvxAcdk1dbO7PN0qK9Df/B3jdEEAAAAAJATymrrcG1Z+UhqCPAAAACQyXTcvJ6uO9My7V36WRmZdvfAhnYgfgjwAAAAAFSpXw46wdI85gWBpl+xGKSRtkJcoc/IGp7yx91JXh0sCAAAAlE340gkzReSjkR33XuMctyIkbZCXKHPyBqdM23vpKYdLAgAAAMTHKb7FVGUFStogDwAAAIBz7Dl64UM6WSw9PgjyXOjorHelXpWg/PG3T75LFXHH0osAAABAxczRC0iKTAheKLdQwXFBkOdCuQFLcoPa6n0EeAAAAPCc6B675pk6U08I8uIkaVPaaNKVli1bunatv7QUn/xw9MumGzs9NWnfJgAAAACllLQ9eYmSeEX86fkbAAAA4FGXpbwv12lW+Tt37btRz4FvW+NksRJW0gZ5bnaUrJLakibS6wURHzPyAAAA4G03p7wsaZInOXlhSVnygkJXR9kQ5LnQ1NSRss3fREQGOF0UAAAAoMKdZD1n8lLYFqdeQW6KciDIAwAAAOCoJaPPibg+f+Th5q87s2e4H0EeAAAAAFe5PDDK/F3udEESVGoyZ9fULRAIOF0UAAAAAFFqBzaEFktnofTSSdrc/JpZc+nSpTJ//nyniwIAAAAgTJrfJ2n+lFASFk3KgpJL2p48AAAAAG6fo3dFqDcPJZe0PXkAAAAA4EX05LnQMXnPS8OAJR86XRAAAAAACYcgDwAAAICrvRroJHdlzjSXX/dnyuG+NVI9e53TxXIthmsCAAAAcC3NrDnBuiR0vU3KzyRiKQY9eS40K3Wo7PI3FJFuThcFAAAAcJT22H0Zdj0nixCmOElbQ25eJ6++bJZtsr/TxQAAAABcSZdVsLNuZgX6y9spZ4Rl5ETSDtdknTwAAAAgMYdv5lj7+qryLEtyA5ajZXKbpO3JAwAAAJB4ohOuvLU3IQv2IcgDAAAAkLAGp7yx91JXh0viHgR5AAAAABLW4JRpe+foTQsN56ye5MsrEOQBAAAASFga1MneKXka7AnT8wjy3Oj6wFCpGawqDzhdEAAAAMDlwnvtNOOmJSLNw+bppfl9SZd5k548F/rEOl4aWvwEAQAAAJTqPFqOE/s0+ihZKa+l3pHf0yfJNXyTIA8AAACAJ5ydPdv8Xa7/rP1W5Ik8qSJ5InfW27eTP13ktjXiZQR5LnRbyvPiTzlARLo5XRQAAAAgMTVoIzusqpImeeLLC5qbTIKWvKCki7clbZA3YcIEswUCAXGb3imzZZuvidPFAAAAABLaSdZzEQul5wSCku5Pye/p87CkDfKGDBlitm3btklGRobTxQEAAAAQZ0uiEq5oQhYN9H7LOkLqy2bPLrmQ4nQBAAAAAKAypPl9pifPpsM3q/t258/Z0+3Hd/ftPHO4yDOJmZUzaXvyAAAAACRrz95P5t8dWfUj5uxd98IC+Xhves4f056Q3KqHSiIiyAMAAACQlE6KmrMXbo+VuKFS4pYcAAAAAMphSRGLpO/JStyqJchzoTeDHcVn1ZJeThcEAAAAQMIhyHOhUcGrpWHQIsgDAAAAUGpk1wQAAAAAD6Enz4U6+L6TDF9Vp4sBAAAAJK03E3gKFUGeCz3uf1C2pTQRkRudLgoAAACQlEYl8BQqhmsCAAAAgIckfJC3fft2adeunRx77LFy9NFHy5NPPul0kQAAAAB4YApVG1/+oumJJuGHa+63337yySefmL9///23tG7dWi688EKpU6eO00UDAAAAkKAeT+ApVAnfk+f3+02Ap3bv3i2BQEAsK/aq9QAAAADgdY4HeXPmzJEePXpIw4YNxefzyfTp0wvsM3HiRGnWrJlUrVpVTjjhBPn0008j7v/rr7+kTZs20qhRI7nlllukbt26lfgKAAAAAMA9HA/ydu7caQK08ePHx7x/6tSpMnToUMnMzJRvvvlGTj31VOnatausXr06tM8BBxwg3377raxcuVJefvllWbdunSSyLVJDtkl+7yQAAAAAJFSQpwHbXXfdZebRxTJ27Fjp37+/DBgwQI466igZN26cNG7cWCZNmlRg3/r168sxxxxjegcT2el5E2RAYITTxQAAAACQgFydeCUnJ0cWLlwow4cPj7i9c+fOMnfuXHNZe+2qVasmNWvWlG3btpkA79prry30mHv27DGbTR+jdB6fW+by+cQym1vKkwzs9586p769ijZOnXsdbZw69zrauAN1LnvPyl10Tl7Ssrg6yNu4caNJpKI9dOH0+p9//mku//7776anz2741113nenNK8yYMWMkOzu7wO1bt251zRvYcD+RelXzy6TzFFHx9L3fsWOHuUydU99eRBunzr2ONk6dex1tvPJtqnak7Kp6oKS56Jzc7qBK6CDPFl2p2sjt2zQRy6JFi0p8rBEjRsiwYcMiKkqHf2ZkZJjeQDeYnjNAtvu1TPNc06C8zg7wtR1Q59S3F9HGqXOvo41T515HG698J28fLg0Dlsx20flhScvh6iBPs2TqEgl2r51t/fr1BXr3SqpKlSpmmzBhgtm0p9CuMLe8eVUkV/ZInqvKlAzs+qbOqW+voo1T515HG6fOvY427sRwTXHV+WFJy+F44pWipKenm566WbNmRdyu1zt06FCuYw8ZMkSWLl0q8+fPL2cpAQAAAMA9HO/J03lQK1asCF3XZRB0+GXt2rWlSZMmZmhlnz59pG3bttK+fXuZPHmyWT5h0KBBjpYbAAAAgHctSL3aTKES6SaJxvEgb8GCBdKpU6fQdXu+XN++fWXKlCnSq1cv2bRpk4wePVrWrl0rrVu3lhkzZkjTpk0dLDUAAAAAL0uXPEmT/KldicbxIK9jx47FZrUcPHiw2eIpek4eAAAAAHiBq+fkVSTm5AEAAADwIsd78lBQn8BIqRVMl8epHAAAAACllLQ9eTpUs2XLltKuXTtxm2+tw2W5pZM8AQAAAKB0kjbIY7gmAAAAAC9iuKYLPegfL+Kvk5DpWgEAAAAv6JPAU6gI8lyos+8r2SZNnC4GAAAAkLS+tQ6XhsWsAuBWSTtc081z8gAAAACgrJI2yGNOHgAAAICiplAN878qiYjhmgAAAADgoSlUSduTBwAAAABeRJAHAAAAAB6StEGemxOvPB3sLtOt05wuBgAAAIAElLRBnpsTrzwSvFReDp7tdDEAAAAAJKCkDfIAAAAAwIvIrulC5/o+l/19+4lIN6eLAgAAACSlp4PdJc06QAZI4iHIc6Ex/idkW4qma73N6aIAAAAASemR4KXSMGglZJDHcE0AAAAA8JCkDfLcnF0TAAAAgPNTqE71LUrItyFpgzw3Z9cEAAAA4PwUqutT3kjItyFpgzwAAAAA8CKCPAAAAADwELJrutDPcrDskQOlrtMFAQAAAJBwCPJc6IK8MdIwYMmHThcEAAAAQMJhuCYAAAAAeEjSBnluXkKhuuySarLb6WIAAAAAST2F6nc5UBJRajIvoaDbtm3bJCMjQ9zki9SBss3fREQudLooAAAAQFK6IIGnUCVtTx4AAAAAeFHS9uQBAAAAQNFTqBITQR4AAAAAeGgKFcM1AQAAAMBDCPIAAAAAwEMI8gAAAADAQ5iT50Ld8x6QugG/vOx0QQAAAAAkHII8F1ot9SVPLKeLAQAAACABMVwTAAAAADwkaXvyJkyYYLZAICBu87R/jAT8B4pIN6eLAgAAACSl7gk8hSppe/KGDBkiS5culfnz54vbtPMtk5ay0uliAAAAAEk9hepPqSOJKGmDPAAAAADwoqQdrgkAAAAAXpxCRZAHAAAAADGmUG2TnZKIGK4JAAAAAB5CkAcAAAAAHsJwTRd6IHiZVLVqyL+cLggAAACAhENPngu9EOwq7wRPcboYAAAAABIQQR4AAAAAeAjDNV2oT8pMqZpSIyHTtQIAAABe8EACT6GiJ8+Fbk55Rfr63nW6GAAAAEDSeiGBp1AlfJD322+/SceOHaVly5ZyzDHHyOuvv+50kQAAAADAMQk/XDM1NVXGjRsnxx57rKxfv16OP/546datm1SvXt3pogEAAABIUH0SeApVwgd5DRo0MJs68MADpXbt2rJ582aCPAAAAADlmkK1zddERO6WROP4cM05c+ZIjx49pGHDhuLz+WT69OkF9pk4caI0a9ZMqlatKieccIJ8+umnMY+1YMECCQaD0rhx40ooOQAAAAC4j+NB3s6dO6VNmzYyfvz4mPdPnTpVhg4dKpmZmfLNN9/IqaeeKl27dpXVq1dH7Ldp0ya58sorZfLkyZVUcgAAAABwH8eHa2rAplthxo4dK/3795cBAwaY6zr/7r333pNJkybJmDFjzG179uyRCy64QEaMGCEdOnQo8vl0X91s27ZtM38tyzKbG3xltZSAHCinuaQ8ycB+/93SBryO+qbOkwHtnPr2Oto49e11lvjyNxedH5a0LI4HeUXJycmRhQsXyvDhwyNu79y5s8ydOzf0Qvv16ydnnHGG9OnTp9hjamCYnZ1d4PatW7e65g3MqnKr1EuzpM3WrWYIKyqevvc7duwwl6lz6tuLaOPUudfRxqlzr6ONV76cak1kZ9WDJM1F5+R2B1VCB3kbN26UQCAg9evXj7hdr//555/m8ueff26GdOryCfZ8vhdeeEGOPvromMfU3r5hw4ZFVJTO4cvIyJCaNWuKG/zxd/5fLZNbGpTX2QE+dU59exVtnDr3Oto4de51tPHKt2fXatedH5a0HK4O8gp7MdrI7dv+8Y9/mGQrJVWlShWzxXoOt7x5jWW91BW/q8qUDOz6ps6pb6+ijVPnXkcbp869jjZeuRZYR+ZPoXLR+aEngry6deuK3+8P9drZdD286N690powYYLZtKfQbf6XerNs82u61uKHnwIAAACIv/6BEdIwYMmHCVi5jmfXLEp6erpZMmHWrFkRt+v14hKsFGfIkCGydOlSmT9/fjlLCQAAAADu4XhPnia7WLFiRej6ypUrZdGiRWZR8yZNmpj5c5pQpW3bttK+fXuzRIIunzBo0CBHyw0AAADAu5rIOjOFKhE5HuTpAuadOnUKXbeTovTt21emTJkivXr1MmvgjR49WtauXSutW7eWGTNmSNOmTT07XBMAAACAs/6XwFOoHA/yOnbsWOzSBYMHDzZbPOlwTd00u6ZmzAEAAAAAL3D1nDwAAAAAQOkQ5AEAAACAhzg+XNMpbp6T1z7vCTkoIPKO0wUBAAAAkHCStifPzUso7JRqskuqOl0MAAAAAAkoaYM8AAAAAPCipB2u6WbTUkfIHn8DEenmdFEAAACApNQ+gadQJW1Pns7Ha9mypbRr107c5jBZI41kvdPFAAAAAJLWzgSeQpW0QZ6b5+QBAAAAQFkxXBMAAAAAPDSFiiAPAAAAAGJModomfklESTtcEwAAAAC8KGmDPDcnXgEAAACAskraIM/NiVdGBAbKY8GLnC4GAAAAgATEnDwXesc6RRpaltPFAAAAAJCAkrYnDwAAAAC8iJ48F7ox5TVJSzkgIdO1AgAAAF4wIjBQ9g/uJ9mSeAjyXKh/yv9km6+J08UAAAAAktY7CTyFKmmHa5JdEwAAAIAXJW2Q5+bsmgAAAACcn0L1z5RZCfk2MFwTAAAAADw0hSppe/IAAAAAwIsI8gAAAADAQxiu6ULvWyeKSB3p7nRBAAAAACQcgjwX+nfgOmkYsAjyAAAAAJRa0g7XZAkFAAAAAF6UtEGem5dQaONbIc19vzldDAAAACCpp1DNk1aSiBiu6UIv+EfLthRN1zrQ6aIAAAAASenfCTyFKml78gAAAADAi+jJAwAAAIAYU6hq+dIlERHkAQAAAICHplAxXBMAAAAAPIQgDwAAAAA8hCDPhXIkVXLF73QxAAAAACQg5uS5UNu8Z0y61g+dLggAAACAhENPHgAAAAB4SNIGeRMmTJCWLVtKu3btnC4KAAAAAJfJSeApVEk7XHPIkCFm27Ztm2RkZIibfJI6RHb6DxaRbk4XBQAAAEhKbRN4ClXSBnluVku2i1/+droYAAAAABJQ0g7XBAAAAAAvoicPAAAAADw0hYogDwAAAAA8NIWK4ZoAAAAA4CEEeQAAAADgIQzXdKFBgX9LRrCqPOx0QQAAAAAkHII8F5prHSMNLcvpYgAAAABIQAzXBAAAAAAPoSfPhUalPCO+lFoJma4VAAAA8IJBCTyFyhM9eRdccIHUqlVLLr74YvGCC1M+ljN9C50uBgAAAJDUU6i+tY6QROSJIO+GG26Q559/3uliAAAAAIDjPBHkderUSWrUqOF0MQAAAAB4aArVwJTpkogcD/LmzJkjPXr0kIYNG4rP55Pp0wtW5MSJE6VZs2ZStWpVOeGEE+TTTz91pKwAAAAAksOFCTyFyvEgb+fOndKmTRsZP358zPunTp0qQ4cOlczMTPnmm2/k1FNPla5du8rq1asrvawAAAAA4HaOZ9fUgE23wowdO1b69+8vAwYMMNfHjRsn7733nkyaNEnGjBlT6ufbs2eP2Wxbt24N/bVcsjbdnj2WbPcFJW3rVtO7iYqn773dFqhz6tuLaOPUudfRxqlzr6ONV749Ljwn37Ztm/lbXNzieJBXlJycHFm4cKEMHz484vbOnTvL3Llzy3RMDQyzs7ML3N6kSRNxl80i2bqMAgAAAABnbHblOfn27dslIyMjMYO8jRs3SiAQkPr160fcrtf//PPP0PUuXbrI119/bYZ+NmrUSKZNmybt2rWLecwRI0bIsGHDQteDwaBs3rxZ6tSp46oIvXHjxvLbb79JzZo1nS5OUqDOqW+vo41T515HG6fOvY42Tp3bPXga4Gk+k6K4OsizRQdf+uLCb9PhmyVVpUoVs4U74IADxI00wCPIo869jDZOnScD2jn17XW0cerb62q67Jy8qB481yReKUrdunXF7/dH9Nqp9evXF+jdAwAAAAC4PMhLT083SybMmjUr4na93qFDB8fKBQAAAABu5fhwzR07dsiKFStC11euXCmLFi2S2rVrm2QoOn+uT58+0rZtW2nfvr1MnjzZLJ8waNAg8SodTpqVlVVgWCmoc6+gjVPnyYB2Tn17HW2c+va6Kgl8Tu6zHF434OOPP5ZOnToVuL1v374yZcqU0GLo999/v6xdu1Zat24tDz/8sJx22mkOlBYAAAAA3M3xIA8AAAAAkCRz8gAAAAAApUOQBwAAAAAeQpAHAAAAAB5CkAcAAAAAHkKQBwAAAAAeQpAHAAAAAB5CkAcAAAAAHkKQBwAAAAAeQpAHAAAAAB5CkAcAAAAAHkKQBwAAAAAeQpAHAAAAAB5CkAcAAAAAHkKQBwAAAAAeQpAHAAAAAB5CkAcAAAAAHkKQBwAAAAAeQpAHAAAAAB5CkAcAAAAAHuKaIG/OnDnSo0cPadiwofh8Ppk+fXqR+69du1b++c9/SosWLSQlJUWGDh1aaWUFAAAAALdyTZC3c+dOadOmjYwfP75E++/Zs0fq1asnmZmZ5nEAAAAAAJFUt1RC165dzVZShxxyiDzyyCPm8jPPPFOBJQMAAACAxOGaIK+yaA+gbrZgMCibN2+WOnXqmGGiAAAAAOBGlmXJ9u3bzRQ3nbJWmKQL8saMGSPZ2dlOFwMAAAAAyuS3336TRo0aFXp/0gV5I0aMkGHDhoWub926VZo0aSK//vqr1KxZU9xAexc3btwodevWLTJCB3WeqGjj1HkyoJ1T315HG6e+vS7ownPybdu2SdOmTaVGjRpF7pd0QV6VKlXMFu2AAw5wVZCXk5NjyuSWBuV11Dn17XW0cerc62jj1LnX0capc2XHBsVNMyOCAAAAAAAPcU1P3o4dO2TFihWh6ytXrpRFixZJ7dq1zXBKHWa5Zs0aef7550P76P32Yzds2GCup6enS8uWLR15DQAAAADgNNcEeQsWLJBOnTqFrtvz5vr27StTpkwxi5+vXr064jHHHXdc6PLChQvl5ZdfNmNUV61aVYklBwAAAAD3cE2Q17FjR5MStDAa6EUran8AAAAASEbMyQMAAAAADyHIAwAAAAAPIcgDAAAAAA8hyAMAAAAADyHIAwAAAAAPIcgDAAAAAA8hyAMAAAAADyHIAwAAAAAPIcgDAAAAAA8hyAMAAAAADyHIAwAAAAAPIcgDAAAAAA8hyAMAAAAADyHIAwAAAAAPIcgDAAAAAA8hyAMAAAAADyHIAwAAAAAPIcgDAAAAAA8hyAMAAAAADyHIAwAAAAAPSS3rA6+++uqYt2dkZEiLFi3kiiuukP333788ZQMAAAAAVFZP3pYtW2JuixYtkpEjR5pA75dffinr4QEAAAAAldmTN23atELv27Vrl1x55ZUyfPhwee2118r6FAAAAAAAN8zJq1atmtx6660yb968ijg8AAAAAKCyE6/Url1b/vrrr4o6PAAAAACgMoO8uXPnymGHHVZRhwcAAAAAxHNO3nfffRfz9q1bt8r8+fPlnnvukbvuuqushwcAAAAAVGaQd+yxx4rP5xPLsgrcV69ePTMnb9CgQWU9PAAAAACgMoO8lStXFrpO3gEHHFDWwwIAAAAAnAjymjZtav7u2bNH8vLypHr16uUpBwAAAADAycQrGzdulHPPPVf2339/qVmzpnTo0IHFzwEAAAAgUYO8ESNGyMKFCyU7O1seeOABE/QNHDgwvqUDAAAAAFTOcM333ntPnnnmGenWrZu5rn9bt24tubm5kpaWVtbDAgAAAACc6Mn7448/5LjjjgtdP/LIIyU9Pd3cDgAAAABIsJ48XTohNTXy4Xo9GAzGo1wAAAAAUD4v9xbZErUqQK1mIv981dM1W64g78wzz4wI9P7++2/p0aOH6dGzff311+UvJQAAAACU1paVIpt/Eal9aP51vZwEyhzkZWVlFbitZ8+e5S0PAAAAAMRP7UNFhnyZf3nCSUlRs3EN8gAAAAAACZp4Zffu3fL222/L9u3bC9y3bds2c58ulA4AAAAASIAg74knnpBHHnlEatSoUeA+XRz90UcflSeffLK85QMAAAAAVEaQ99JLL8nQoUMLvV/ve/7558t6eAAAAABAZQZ5P/30k7Rp06bQ+4855hizDwAAAAAgARKv5OXlyYYNG6RJkyYx79f7dB8AAAAAKFKSrmfnup68Vq1ayezZswu9f9asWWYfAAAAACjRenY2vRwd9KHie/KuvvpqGTZsmAnkunfvHnHff//7X7nrrrtk7NixZT08ACTvL5elwa+cAACvrmengZ69rh3/31VOkPd///d/MmfOHDnvvPPkyCOPlBYtWojP55Nly5bJ8uXL5dJLLzX7AEBSKy6I2/BD/t96R5b+2OG/eAIoO4aJAe6jQV1l/H/3ctT/0x4JJssc5KkXX3zRBHkvv/yyCewsyzLBXnZ2tgnyAECSPZArLojT28v6Hwq/cgLx+cxGf075AQWIz2dLFfZ/XPS++rnTnjxb+GPs3ryKHCZa+9D87wLdQs/nE6nbQeSSByWpgjylwRwBHYCkUtwJYryCOLf8ygkkoqJ+nQ8/qYv1OS3qhJJeP6Bw0Z+tov5vit5X/4b/v+bEMNGXowPPlSJVD5NEVO4gDwAk2YO6igzkihL9K2f43AUPDTkBSuTVy0V2pYjs/llErMjPaawTzfC5P7FEf55s9PrBi+I5ZDF6Xl1J93WDf0a95gknS6IiyAOA4hT3q78bRP/6Sc8ekrF3Lnd/kTQp+Dkt7VCvonoTStPr5zUFejETdygbivh/jv8/PIEgDwCKO5mx/+Nz06+Nxf76mEQnnkjOXvTouTObV4k0PEvkqvdFUlKK7pmLnvsTraw/4Hh9KGeBoXiJO5QNMdj/z5Xm/4/i5tXBMQR5AFDcyYyT8wTKI2K4WRL/4u7RzGmS7EOjo/et3Uxk//qxjxv9+Y33Z9r+rEWXt0ASBw+0v4iheHEayub14NiNShOcFfX+ODWvjjZTcUFew4YNpWfPnia75plnninp6ellPRQAuJPbe+9KNXwziX5xLypQiB6GxIlC4g6Njr4tGBRZvz72cSsyWAj/rEWXN9aJNMqXrKM4/KhTMsUFZ+E/EhY3D9WJ/yvj2WY8qsxBni6boIue33DDDbJu3Trp0qWLCfjOPfdcqV27dnxLCQAVzWtDTmJNHt++TmRSh/ykFF76lbw0vT/RCWpIouHOz10i/bhS1OcomYZRF7cmqCrqe6c0yToqam6ZlwLEkvyAVdhnLfpHwnjNQ433/7MsnF4xQV7Hjh3N9tBDD8mSJUvk7bfflgkTJsiAAQOkffv2oV6+ww5Lkl+OASQ2N6Vyrgj6WkzmwR3e+8WzNL0/JTl5KSyrYqKf9LmR1z93sYS3Ly+1p+j3Mlr4905pT/ZLG3wVNresuMDHS8lHytPTFe82Wdhw5nh+3llSqGLm5LVq1cpsI0aMML16GvDplpmZKYceeqjcd999poevKHPmzJEHHnhAFi5cKGvXrpVp06bJ+eefX+RjPvnkExk2bJgJMnX46C233CKDBg2Kx0tCZWGoFNwk0XoQSqP3S/lD2Q48UGRS+6IDmUQ4+SxrYpziXlNRJxxenF/ldA9PIvbclUd4+3KqPVXk/7tFvZfhr7O0wX28gq+SBD5lST7iVvHqHY3F/j+kuAC9qOHM8VRZC6cnc+KV+vXryzXXXGO2v//+W9577z2pUqVKsY/buXOntGnTRq666iq56KKLit1/5cqV0q1bN/M8L774onz++ecyePBgqVevXokeD5eI/sKN/k+PEyhUFK8NzyyN4n45deMv2EUNyYznr8FFnXzEKkP0SXoEX/48SHvdtmT9PiuqhycZeu7Chb//JWlPFdFm4j2XyR4GXtrv0NKsUxj+Y0BpT+ALO44qTxbJaE58vp36oTz8M1vcZzgZv/OSIbvmfvvtJxdccEGJ9u3atavZSurxxx+XJk2ayLhx48z1o446ShYsWCAPPvggQV4iz7+I/nW+NMcq6sstXl+E9Dx6RzIOE7MV1+6jhyy6IThxw1qF0c9VkjlItpIkfAnnhjqPp2TqrYtXe6rIH1uKmstUmrYYPgy8JN+hZen9KS4pSFFljGdW1aJ+rHAqoVNpA/bSLCNSFKe+m+JV/iSRsEsofPHFF9K5c+eI2zT5y9NPPy25ubmSlmavhhppz549ZrNt27bN/A0Gg2ZzAy2HZVmuKU+pvXp57C83HS62ZVX+Wkaa6lrVPiz/Pvu19n5532PML4MrC0/RvOHH/L/1Wuy7vmH5vv3t51TRz6vH1V/Z9XlfvVyCW1aJVfUwCYb/4m4/Nvz1RD9n+HHiWU8e54427stvf9fOjbw5UT938azzWvofpy/258qxdurC9yv8+6qwOt+wQYL16ok88Q/7xtjfSeFi1Xm0hPqu8FXKe+WO75U4tieTJCms3qL/zyhzG4h6P8I/78W1RbvN7n1s8NIX9rVxe23Cwuo//Hmi/+8v6WdL948ub1FlLOo4seoi4nqM+2J9B8V6r4o65yjHuUHBNh5VpohyFPM+F/ceuE2pyl/U+1qS6/sExSeW+Fz1vVLSsvgsbS0u4/P5ip2T17x5c+nXr5/cdtttodvmzp0rp5xyivzxxx/SoEGDmI8bNWqUZGdnF7h94sSJUq1aNXEDfUvy8vIkNTXV1EXCWbdEJG+PSOreYbr25fqt8u9Terk4m1bkP7Yoetw6hxfcP/w57TKFP294GXN3mQ9wXpUDJDW4W3y5u/L3SdvbHqKvhz9n9HFLUmZb9HHDyxzrOOHPm+Bc0cZL0xaTuc6j22J0u62othn9vNGf6USr8/VLC/9ejFbc90is98BN3xUOvXeu+F6Jp2L+3yrz57C8/w+HPY8r6rw8/1/GqmP7elH3leY4xT021rHCvyvC3ms9Yc9LrSGptRvn13dpyp9MSvt+rCu8rqx1SyQvpaqk1j3UNd8ru3btMlPUtm7dKjVr1vRmkKfz9zTZi03n5f3jH/8wiVsOOuigEvfkNW7cWLZs2VJkRVV2hL5hwwYzvzDF/mUskdi/JIX/sqS/YukvWvbfWL+ExbsM9nOp6OeN+qUsWOtQ2dDpofw6f61PyXvYinqtJRF+3PBjlbL8RZbRhRxr4+H1Vllt0Wt1Ht32onu3y9oWK+q4bqnz0nyvFCfW59/mhnYd63uwEt67hP+/M1pR/w9Et4Hi/o8o7P+aZK/zWP+H29eLuq80xynuscUdK+y9DG7+VTY0PFPq9Z2SX9+lKX8yKe37ManwugpOOkU2VD1sX527gMYutWrVKjbIq5Dhmvqhf+edd8zQyenTp1fEU5gg7s8//4y4bf369ebXpDp16hT6OE0CEysRjL5xbnnz7EDXbWUq9Tw7u+y1DskfAqlb7UPyr1f06wp/ThX9vP98JXL/YFB869fn13n0fUWy8l+vZiuMmGM4rwyFDjtW9HHMHIaf9z1PzLW9rPzXV961irzcxrf8sq9uK6stukhc6jz68xFqb1ZY8ohl+9pqUQtDh98f/t6oes1d0U7jVuel+l4pRlHHiv6ucCJjY7m+B5Po/87y/D8W3QbMPKW9/weo8M9T+P8Pxt7jxamOErvOo+si/HpR95XmOFHXi52vF/XY8Pd6wslm8OC++o7xPNHnEQn5vpRXKd6PmNcjjxVZ584raTniGuT99NNP8swzz8hzzz1nesZ0jlxF0bX4dDH2cO+//760bdu20Pl4KIei1qgpLoGFEydplfWc8ZzUXVS2qtIsTFqatYqSEQkgKi95RHS2wOgfJ8LvT7ZU+hUl+ruisj7v4d87yZTEyM3/j3lpKYBEUNKkIOXNbmpnM7UDuvDnKU3WS3heajzGhb722mum127evHkSCATk4Ycflquvvlr233//Eh9nx44dsmLFioglEhYtWiS1a9c2WTR1WOaaNWvk+eefN/frenjjx4836+TpMgqaiEXL8Morcfy1FLG/lKJP3JL55CyewWRRxyrJ80RnLStqraKSZk4s7QK0QFEp4qN/nAi/nxOS+Ij+fFbUCX5xmZLhLmQldNcPvkVlNy0qQAzPZhrrefj/GfEI8r766it56qmnZOrUqWZ+3BVXXCGvv/66NGrUSM4666xSBXhKlz/o1KlT6LoGb6pv374yZcoUM89u9erVofubNWsmM2bMkJtuukkmTJhgFkN/9NFHWT6hItlfStH/uXNy5rzS/HpX3GK84cJ7XhKxBzCZ18Jzg+JOODghqRxlTY9flGRegsStCgsU4jniBPH7IdZWmvdH51WuXy9y4IFJOgwTlRLkdejQQa6//noT7LVoETYxvow6duxosjQVRgO9aKeffrp8/fXX5X5ulBInZon9nhTV0xItvOcl1lpKbu/Z40QUya6ok/ny/nBDz517FBUouPk7OlkU9UMs70/F2MyaemUO8s444wwzPFKTnfTp08fMv3NLalHEAT0gyaE0/7k4NdenvDgRRTIr6jNemqGc/J/gbgQKFRMMJPri4cmK3uvyBXma5OS3336TZ599Vq699lozN69Xr17mPoI9D6AHBE7N9QHgjqGc4aKT5jDkD14PBggUEhdBdfkTr+j6ciNHjjTbrFmzTGZNXcKgZ8+ecvHFF5vt+OOPL89TwEn0gCDR0NsAlFxp5mVFJ80BvKC8Cc8AF4vbEgpnn3222XTphBdffNEEfPfdd5/JtgkAlYIeaKDkOIkFAM+K+2LougK7JmTRjaQoCYQeEJSFGxOx0AMNAACSPDFLuYK8YDBosl6++eabsmrVKjMXT5c20GGamoyFoZoJhB4QlFaiJmIBAADJq1ZyLCtS5iBPlzs477zzzFp1bdq0kaOPPtrctmzZMunXr58J/KZPnx7f0qLieu9YyBaJmoglVjsGAACIxekRR24P8rQHb86cOfLBBx9ELGKuPvzwQzn//PPl+eeflyuvvDIe5URF99559FcMJAHaMQAAQHyCvFdeeUVuu+22AgGevYbe8OHD5aWXXiLIczvmL8ELaMcAAAAhKVJG3333nZxzzjmF3t+1a1f59ttvy3p4AIk8kdnedCglAAAAEqMnb/PmzVK/fv1C79f7dDkFAEmishKxkAkWAACgYoI8Xf9OFz4vjN/vl7y8vLIeHkCiqaxELGSCBQAAqLjsmppFs0qVKjHv37NnT1kPDQBFYw4eAABA/IO8vn37FrsPmTVdhmFuAAAAgOeVOch79tln41sSVDyGuSER8eMEAABA5QR5SFAMc4MT2TbtxCxlWYCUHycAAAAqJ8jT9fF8Pl+B2zMyMqRFixYyZMgQady4cVkPD8BL2TbLm2mTHycAAAAqPsg79thjY97+119/yYwZM2T8+PHy2WefFbofAI8L77WrqEybAAAAiF+Q9/DDDxd5v/bk3XbbbSbgA4ASYw4eAACAO+fkDRw4ULp06VJRh0dJcLKMRMQcPAAAAHcGedWqVZPdu3dX1OFRksBuww/5f+sduW9eU/g8KcCpJCziE6nbQeSSBwu2W92POXgAAADuC/Lef/99ad68eUUdHiXpBdHgrqwZDYF4iv5xYfNKkaqHxW63/BgBAADgTJD39ttvx7x969atMn/+fHn66adlypQp5SkbyopeELhN9A8NE04W2b5OZFIHEbHovQMAAHBDkHf++efHvL1GjRpy5JFHmgDvkksuKU/ZAHi5Z29XisjuHfnX6b0DAABwPsgLBoPxKwWA5NL7JZH160UOPFAkJcXp0gAAAHgKZ1cAAAAA4CFlDvI+/PBDadmypWzbti3mvLxWrVrJp59+Wt7yAQAAAAAqI8gbN26cXHPNNVKzZs0C92VkZJh18saOHVvWwwMAAAAAKnNO3rfffiv33Xdfofd37txZHnxw7xpYqDgseA4AAAAgHj1569atk7S0tELvT01NlQ0bNpT18Cgpe30xG1kKAQAAgKRW5p68gw8+WBYvXiyHH354zPu/++47adCgQXnKhpJiXTwAAAAA5e3J69atm4wcOVJ2795d4L5du3ZJVlaWdO/evayHBwAAAABUZk/e7bffLm+++aY0b95crrvuOmnRooX4fD5ZtmyZTJgwQQKBgGRmZpb18AAAAACAygzy6tevL3PnzpVrr71WRowYIZZlmds10OvSpYtMnDjR7AMAAAAASIAgTzVt2lRmzJghW7ZskRUrVphA74gjjpBatWrFr4QAAAAAgMoJ8mwa1LVr1y4eh0JxWDIBAAAAQEUkXoFDWDIBAAAAQEX35KGSsWQCAAAAgELQkwcAAAAAHkKQBwAAAAAeQpAHAAAAAB5CkAcAAAAAHkLiFbdjyQQAAAAApUBPntuxZAIAAACAUqAnLxGwZAIAAACAEqInDwAAAAA8hCAPAAAAADyE4Zpu9OrlIrtSRHb/LLL5l/zhmgAAAABQAgR5bk22kru/SNre+Xi1mjldIgAAAAAJgiDPrWrUF7nqfZEURtQCAAAAKDkiCAAAAADwEII8AAAAAPAQgjwAAAAA8BCCPAAAAADwENcFeRMnTpRmzZpJ1apV5YQTTpBPP/20yP0nTJggRx11lFSrVk1atGghzz//fKWVFQAAAADcxlXZNadOnSpDhw41gd4pp5wiTzzxhHTt2lWWLl0qTZo0KbD/pEmTZMSIEfLkk09Ku3bt5KuvvpJrrrlGatWqJT169HDkNQAAAACAk1zVkzd27Fjp37+/DBgwwPTOjRs3Tho3bmyCuVheeOEFGThwoPTq1UsOPfRQ6d27t3n8fffdV+llBwAAAAA3cE2Ql5OTIwsXLpTOnTtH3K7X586dG/Mxe/bsMcM6w+mwTe3Ry83NrdDyAgAAAIAbuWa45saNGyUQCEj9+vUjbtfrf/75Z8zHdOnSRZ566ik5//zz5fjjjzdB4jPPPGMCPD1egwYNYgaGutm2bt1q/v71118SDAbFDYK7ArLNypP0v/6SFBZDr5w6DwZl27Ztkp6eTp1T355EG6fOvY42Tp17HW3cgTrf5b5zcj1fVZZlJUaQZ/P5fBHX9QVE32a74447TAB48sknm/00IOzXr5/cf//94vf7Yz5mzJgxkp2dXeD2pk2birssELn+NacLAQAAACSxBa48J9++fbtkZGQUer/PKi4MrMThmvvtt5+8/vrrcsEFF4Ruv/HGG2XRokXyySefFPpY7blbt26d6bmbPHmy3HrrraZnLlbEHd2Tp7+KbN68WerUqVNoMOlEhK5zEX/77TepWbOm08VJCtQ59e11tHHq3Oto49S519HGqXOloZsGeA0bNiyyd9E1PXk6TE6XTJg1a1ZEkKfXe/bsWeRj09LSpFGjRubyq6++Kt27dy/0RVepUsVs4Q444ABxIw3wCPKocy+jjVPnyYB2Tn17HW2c+va6mi47Jy+qB891QZ4aNmyY9OnTR9q2bSvt27c3vXKrV6+WQYMGmft1uYQ1a9aE1sJbvny5SbJy0kknyZYtW0x2zu+//16ee+45h18JAAAAADjDVUGeLoWwadMmGT16tKxdu1Zat24tM2bMCM2X09s06LNpopaHHnpIfvzxR9Ob16lTJ5OJ85BDDnHwVQAAAACAc1wV5KnBgwebLZYpU6ZEXNe19L755hvxGh1OmpWVVWBYKahzr6CNU+fJgHZOfXsdbZz69roqCXxO7prEKwAAAACA8nPHgg8AAAAAgLggyAMAAAAADyHIAwAAAAAPIchzoYkTJ0qzZs2katWqZu3ATz/91OkiJZxRo0aZxe3Dt4MOOih0v05F1X10Iclq1apJx44dZcmSJRHH2LNnj1x//fVSt25dqV69upx33nny+++/O/Bq3GnOnDnSo0cPU4dav9OnT4+4P151rMuj6NIquiaMbnr5r7/+kmRUXJ3369evQLs/+eSTI/ahzktuzJgx0q5dO6lRo4YceOCBcv7555tszuFo55Vb37Tx+Jo0aZIcc8wxoTXAdPmqmTNnhu6nfVd+ndPGK/57xufzydChQ73fzjXxCtzj1VdftdLS0qwnn3zSWrp0qXXjjTda1atXt3799Veni5ZQsrKyrFatWllr164NbevXrw/df++991o1atSw3njjDWvx4sVWr169rAYNGljbtm0L7TNo0CDr4IMPtmbNmmV9/fXXVqdOnaw2bdpYeXl5Dr0qd5kxY4aVmZlp6lC/SqZNmxZxf7zq+JxzzrFat25tzZ0712x6uXv37lYyKq7O+/bta+orvN1v2rQpYh/qvOS6dOliPfvss9b3339vLVq0yDr33HOtJk2aWDt27AjtQzuv3PqmjcfX22+/bb3zzjvWjz/+aLbbbrvNnIPoe6Bo35Vf57TxivPVV19ZhxxyiHXMMceY82ubV9s5QZ7LnHjiiaYhhTvyyCOt4cOHO1amRA3y9MMXSzAYtA466CDzobbt3r3bysjIsB5//HFz/a+//jJfuhp029asWWOlpKRY7777biW8gsQSHXDEq471hw499rx580L7fPHFF+a2H374wUpmhQV5PXv2LPQx1Hn56A9FWu+ffPKJuU47r9z6VrTxilerVi3rqaeeon07UOeKNl4xtm/fbh1xxBEmSDv99NNDQZ6Xv8cZrukiOTk5snDhQuncuXPE7XpdF3lH6fz000+m612Hvvbu3Vt++eUXc/vKlSvlzz//jKhnXf/k9NNPD9Wzvg+5ubkR++ixWrduzXtRAvGq4y+++MIMeTjppJNC++jwQ72Nz0RsH3/8sRnq1rx5c7nmmmtk/fr1ofuo8/LZunWr+Vu7dm3auQP1baONV4xAICCvvvqq7Ny50wwh5Hu88uvcRhuPvyFDhsi5554rZ511VsTtXm7nrlsMPZlt3LjRfODr168fcbte1waIktMP2fPPP29OdNetWyd33XWXdOjQwYyxtusyVj3/+uuv5rLuk56eLrVq1eK9KIN41bH+1YAlmt7GZ6Kgrl27yiWXXCJNmzY1/3HdcccdcsYZZ5j/oPQ/Leq87LTzdNiwYfKPf/zD/MdOO6/8+la08fhbvHixCTB2794t+++/v0ybNk1atmwZOjHle7zy6lzRxuPv1Vdfla+//lrmz5+fVOcrBHkupBNCo/+zi74NRdMvSdvRRx9tvkwPO+wwee6550KJKMpSz7wXpROPOo61P+9DbL169Qpd1hPjtm3bmoDvnXfekQsvvJA6L4frrrtOvvvuO/nss88K3Ec7r7z6po3HX4sWLWTRokUmQcQbb7whffv2lU8++SR0P+278upcAz3aeHz99ttvcuONN8r7779vEhoWxovtnOGaLqIZe/x+f4GIX4dbRf/CgNLRTEga7OkQTjvLZlH1rPvo8FnNlMR7UXrxqmPdR3tio23YsIHPRAk0aNDABHna7qnzstOMam+//bZ89NFH0qhRo9DttPPKre9YaOPlpz0Uhx9+uPlRSDMPtmnTRh555BHatwN1HgttvHwWLlxozis0W31qaqrZNKB+9NFHzWX7fMOL5ysEeS770GsjnDVrVsTtel2HGqLsNPXtsmXLzJelztHTD2N4PeuHVz/0dj3r+5CWlhaxz9q1a+X777/nvSiBeNWx9sDqvJyvvvoqtM+XX35pbuMzUbxNmzaZXzG13VPnpae/wGqP0ptvvikffvihade0c+fqOxbaeMW8D/p/Jt/jlV/nsdDGy+fMM880w2O159TeNLi+/PLLzeVDDz3Uu+crjqR7QbFLKDz99NMmU8/QoUPNEgqrVq2i1krhX//6l/Xxxx9bv/zyi8l0pClsNT2uXY+aRUkzJ7355psmXe5ll10WM11uo0aNrNmzZ5t0uWeccQZLKERlqvrmm2/Mpl8lY8eONZft5T7iVceakljTHWuWKt2OPvropF1Coag61/u03Wva5pUrV1offfSR1b59e5PymTovm2uvvda0Yf0uCV+W4u+//w7tQzuvvPqmjcffiBEjrDlz5pjvjO+++86k89eMge+//765n/ZduXVOG68cp4dl1/RyOyfIc6EJEyZYTZs2tdLT063jjz8+In00SsZe40QD5oYNG1oXXnihtWTJktD9mjJXl1nQtLlVqlSxTjvtNPPBDrdr1y7ruuuus2rXrm1Vq1bNfFBXr17NW7CXBhEaaERvmv45nnWs67xdfvnlJkjXTS9v2bIlKd+HoupcT4Q7d+5s1atXz7R7XV9Mb4+uT+q85GLVtW66lpuNdl559U0bj7+rr746dL6h3x1nnnlmKMBTtO/KrXPauDNBXtCj5ys+/ceZPkQAAAAAQLwxJw8AAAAAPIQgDwAAAAA8hCAPAAAAADyEIA8AAAAAPIQgDwAAAAA8hCAPAAAAADyEIA8AAAAAPIQgDwAAAAA8hCAPAIBK1rFjR/H5fGZbtGiR4/Xfr1+/UHmmT5/udHEAAOVEkAcAcLXwACR8O+eccySRXXPNNbJ27Vpp3bp1gfs6d+4sfr9f5s2bV2R9pKWlSf369eXss8+WZ555RoLBYMS+hQVtQ4cONYGm7ZFHHjFlAQB4A0EeAMD1NKDTICR8e+WVVyr0OXNycir0+Pvtt58cdNBBkpqaGnH76tWr5YsvvpDrrrtOnn766SLrY9WqVTJz5kzp1KmT3HjjjdK9e3fJy8srdVkyMjJMWQAA3kCQBwBwvSpVqpggJHyrVatWRI/VU089JRdccIEJno444gh5++23I46xdOlS6datm+y///6m96tPnz6ycePG0P3as6WB1bBhw6Ru3bqmd0zpcfR41apVM8HUc889Z57vr7/+kp07d0rNmjXlP//5T8Rz/fe//5Xq1avL9u3bS/1an332WROsXXvttTJ16lTzHIXVx8EHHyzHH3+83HbbbfLWW2+ZgG/KlCmlfk4AgLcQ5AEAPCE7O1suvfRS+e6770wwd/nll8vmzZvNfdrrdfrpp8uxxx4rCxYskHfffVfWrVtn9g+nAZz2rH3++efyxBNPmJ6yiy++WM4//3wzd27gwIGSmZkZ2l8Dud69e5vALJxe18fVqFGjVK/Bsizz2CuuuEKOPPJIad68ubz22msleuwZZ5whbdq0kTfffLNUzwkA8B6CPACA6/3vf/8zPXDh25133llgrtpll10mhx9+uNxzzz2mB+yrr74y902aNMn0eOntGjwdd9xxZg7bRx99JMuXLw8dQx97//33S4sWLcx+jz/+uLn8wAMPmL8a0OnzhBswYIC899578scff5jr2juo5b366qtL/Tpnz54tf//9t3Tp0sVc12CvsCGbsWiZNTAFACQ3gjwAgOvpMEntSQvfhgwZErHPMcccE9HDpr1o69evN9cXLlxoArrwIFEDIvXzzz+HHte2bduIY/7444/Srl27iNtOPPHEAtdbtWolzz//vLn+wgsvSJMmTeS0004r9evUgK5Xr16heXoatH755ZemHCXtCdShpACA5BY52xsAABfSoE172YqimSbDabBjZ5vUvz169JD77ruvwOMaNGgQ8TzFBU16WzTtzRs/frwMHz7cDLe86qqrSh1s6dBSzYSZm5treh5tgUDA9DrGKnu0ZcuWSbNmzULXNdDdunVrgf10PqEmWwEAeBM9eQAAz9OhmkuWLJFDDjnEBIvhW3RgF057++bPnx9xm87pi6bDKjUr5qOPPmqep2/fvqUu40svvSSNGjWSb7/9NqLHcty4cWauYHFZMz/88ENZvHixXHTRRUWWX4NU7dnU4acAAG8iyAMAuN6ePXvkzz//jNjCM2MWR4d2ak+ZDn/UeXq//PKLvP/++2benPaUFUYTrfzwww9y6623mrl7mgTFzl4Z3lOnmT4vvPBCufnmm80adxqslWWopiZr0XXzwjcto/a8vfPOOwXqY82aNfL111+buYY9e/Y0WTmvvPLK0H7//ve/zXG1l1HLrwGkZhDVIarRw10BAN5BkAcAcD3NhqnDKsO3f/zjHyV+fMOGDU3GTA3oNKmJBk+6rpwOWUxJKfy/Qh36qMsjaMZKnfOnwyjt7Jq6jEG4/v37m7X1ypJwRXvWNAAL74ULH3KpgWN4Aha7PrRnUtfM0/mG2ouoyyjoIuo2zR6qQan2BOrcQj2OBniffvqpNG3atNTlBAAkBp8Va3IBAACI6e677zZZN3/77bcCwy01cNQsm+np6UXWnq7Jp8s56FBMN9HeyWnTppklIwAAiYuePAAAijBx4kQzr02HeGrmTF1OIXzOnS55oPPwxowZY4Z3FhfghR9Xs3zqPDqnDRo0yJQFAOAN9OQBAFCEm266SaZOnWrm9OnSCH369JERI0aEljkYNWqU6d3TJRN0uGRJgiWdS7dr1y5zWY9Z0sCwouhSE9u2bTOXdRhoUcloAADuR5AHAAAAAB7CcE0AAAAA8BCCPAAAAADwEII8AAAAAPAQgjwAAAAA8BCCPAAAAADwEII8AAAAAPAQgjwAAAAA8BCCPAAAAADwEII8AAAAABDv+H88qGiB20n1CAAAAABJRU5ErkJggg==",
                        "text/plain": [
                            "<Figure size 900x650 with 2 Axes>"
                        ]
                    },
                    "metadata": {},
                    "output_type": "display_data"
                }
            ],
            "source": [
                "fig, (ax_spec, ax_ratio) = plt.subplots(\n",
                "    2, 1, figsize=(9, 6.5), sharex=True, gridspec_kw={'height_ratios': [3, 1]})\n",
                "\n",
                "edges = hist_cpu.axes[0].edges\n",
                "cpu_vals, cuda_vals = hist_cpu.values(), hist_cuda.values()\n",
                "\n",
                "ax_spec.stairs(cpu_vals, edges, label=f'CPU MT  ({n_clusters_cpu:,})',\n",
                "               color='C0', linewidth=1.4)\n",
                "ax_spec.stairs(cuda_vals, edges, label=f'CUDA  ({n_clusters_cuda:,})',\n",
                "               color='C1', linestyle='--', linewidth=1.4)\n",
                "ax_spec.set_ylabel('Counts')\n",
                "ax_spec.set_yscale('log')\n",
                "ax_spec.set_ylim(bottom=10)\n",
                "ax_spec.set_title(f'Cluster energy spectrum — {N:,} frames')\n",
                "ax_spec.legend()\n",
                "ax_spec.grid(alpha=0.2)\n",
                "\n",
                "with np.errstate(divide='ignore', invalid='ignore'):\n",
                "    ratio = np.where(cpu_vals > 0, cuda_vals / cpu_vals, np.nan)\n",
                "ax_ratio.stairs(ratio, edges, color='C1', linewidth=1.2)\n",
                "ax_ratio.axhline(1.0, color='gray', linewidth=0.5)\n",
                "ax_ratio.set_ylabel('CUDA / CPU')\n",
                "ax_ratio.set_xlabel('Energy [ADU]')\n",
                "ax_ratio.set_ylim(0.9, 1.1)\n",
                "ax_ratio.grid(alpha=0.3)\n",
                "\n",
                "plt.tight_layout()\n",
                "plt.show()"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "md-sets",
            "metadata": {},
            "source": [
                "## Cluster sets\n",
                "\n",
                "Every cluster is keyed by `(frame, x, y)`, so the two results can be differenced\n",
                "directly — no per-frame Python loop.\n",
                "\n",
                "A cluster whose centre moved by one pixel was *relabelled*, not gained or lost,\n",
                "so those are counted separately from genuine adds and drops."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "id": "sets",
            "metadata": {},
            "outputs": [],
            "source": [
                "def cluster_keys(vectors):\n",
                "    \"\"\"One int64 key per cluster: frame * rows*cols + y * cols + x.\n",
                "\n",
                "    Centres are unique within a frame, so keys are globally unique and the sets\n",
                "    below need only a sort, never a np.unique.\n",
                "    \"\"\"\n",
                "    n = len(vectors)\n",
                "    sizes = np.fromiter((cv.size for cv in vectors), np.int64, n)\n",
                "    frames = np.repeat(\n",
                "        np.fromiter((cv.frame_number for cv in vectors), np.int64, n), sizes)\n",
                "    xy = [np.asarray(cv) for cv in vectors if cv.size]\n",
                "    if not xy:\n",
                "        return np.empty(0, np.int64), frames\n",
                "    x = np.concatenate([v['x'] for v in xy]).astype(np.int64)\n",
                "    y = np.concatenate([v['y'] for v in xy]).astype(np.int64)\n",
                "    return frames * (rows * cols) + y * cols + x, frames\n",
                "\n",
                "\n",
                "t0 = time.perf_counter()\n",
                "k_cpu, fr_cpu = cluster_keys(clusters_cpu)\n",
                "k_cuda, fr_cuda = cluster_keys(clusters_cuda)\n",
                "\n",
                "# Sort ONCE. np.setdiff1d and np.isin each re-sort their argument on every call,\n",
                "# and the shift check below calls one eight times per direction -- about twenty\n",
                "# sorts of a 27-million-element array. Two sorts plus binary search instead.\n",
                "s_cpu = np.sort(k_cpu)\n",
                "s_cuda = np.sort(k_cuda)\n",
                "\n",
                "\n",
                "def present_in(vals, sorted_keys):\n",
                "    \"\"\"Which of `vals` appear in an already-sorted key array.\"\"\"\n",
                "    if sorted_keys.size == 0:\n",
                "        return np.zeros(vals.size, dtype=bool)\n",
                "    i = np.searchsorted(sorted_keys, vals)\n",
                "    np.clip(i, 0, sorted_keys.size - 1, out=i)\n",
                "    return sorted_keys[i] == vals\n",
                "\n",
                "\n",
                "cpu_only = s_cpu[~present_in(s_cpu, s_cuda)]\n",
                "cuda_only = s_cuda[~present_in(s_cuda, s_cpu)]\n",
                "\n",
                "\n",
                "def drop_shifted(only, sorted_other):\n",
                "    \"\"\"Drop mismatches that have a counterpart within 1 px in the other set.\"\"\"\n",
                "    keep = np.ones(only.size, dtype=bool)\n",
                "    for dy in (-1, 0, 1):\n",
                "        for dx in (-1, 0, 1):\n",
                "            if dx or dy:\n",
                "                keep &= ~present_in(only + dy * cols + dx, sorted_other)\n",
                "    return only[keep]\n",
                "\n",
                "\n",
                "cpu_real = drop_shifted(cpu_only, s_cuda)\n",
                "cuda_real = drop_shifted(cuda_only, s_cpu)\n",
                "touched = np.unique(np.concatenate([cpu_real, cuda_real]) // (rows * cols)).size\n",
                "\n",
                "print(f'CPU  {k_cpu.size:>12,} clusters')\n",
                "print(f'CUDA {k_cuda.size:>12,} clusters      net {k_cuda.size - k_cpu.size:+,}'\n",
                "      f'  ({abs(k_cuda.size - k_cpu.size) / k_cpu.size:.4%})')\n",
                "print()\n",
                "print(f'exact-position mismatches   CPU-only {cpu_only.size:>8,}   '\n",
                "      f'CUDA-only {cuda_only.size:>8,}')\n",
                "print(f'  of which shifted by 1 px  {cpu_only.size - cpu_real.size:>17,}   '\n",
                "      f'{cuda_only.size - cuda_real.size:>17,}')\n",
                "print(f'  genuine add / drop        {cpu_real.size:>17,}   {cuda_real.size:>17,}')\n",
                "print(f'\\nframes touched by a genuine difference: {touched:,} of {N:,}')\n",
                "print(f'[{time.perf_counter() - t0:.1f}s]')"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "md-drift",
            "metadata": {},
            "source": [
                "## Where in the run it accumulates\n",
                "\n",
                "Cumulative `CUDA − CPU` against frame index. A curve that rises and then flattens\n",
                "means the two pedestals settled into agreement and the difference is a start-up\n",
                "transient; a steadily sloping line would mean a persistent per-frame bias."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "id": "drift",
            "metadata": {},
            "outputs": [],
            "source": [
                "cnt_cpu = np.bincount(fr_cpu, minlength=N)\n",
                "cnt_cuda = np.bincount(fr_cuda, minlength=N)\n",
                "cum = np.cumsum(cnt_cuda.astype(np.int64) - cnt_cpu.astype(np.int64))\n",
                "\n",
                "fig, ax = plt.subplots(figsize=(9, 3.6))\n",
                "ax.plot(cum, color='C1', linewidth=1.2)\n",
                "ax.axhline(0, color='gray', linewidth=0.5)\n",
                "ax.set_xlabel('Frame')\n",
                "ax.set_ylabel('cumulative CUDA - CPU')\n",
                "ax.set_title(f'Cluster-count difference accumulating over {N:,} frames '\n",
                "             f'(final {cum[-1]:+,})')\n",
                "ax.grid(alpha=0.25)\n",
                "plt.tight_layout()\n",
                "plt.show()\n",
                "\n",
                "half = N // 2\n",
                "print(f'first half: {cum[half - 1]:+,}      second half: {cum[-1] - cum[half - 1]:+,}')"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "id": "e7b0da01-1787-4260-890d-533accda0e17",
            "metadata": {},
            "outputs": [],
            "source": []
        }
    ],
    "metadata": {
        "kernelspec": {
            "display_name": "Python 3 (ipykernel)",
            "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.11.15"
        }
    },
    "nbformat": 4,
    "nbformat_minor": 5
}