{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "179f6bb1-ef6f-4328-b857-d85740b69b9b",
   "metadata": {},
   "source": [
    "# Compare 2 videos with SSIM and PSNR metrics\n",
    "\n",
    "The classic metrics **PSNR** (Peak Signal to Noise Ratio) and SSIM (Structural SIMilarity) are implemented in C for improved performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "be964f70-cf5f-4fd5-871e-188f4955e2d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import pathlib\n",
    "import timeit\n",
    "\n",
    "import cv2\n",
    "import numpy as np\n",
    "import tqdm\n",
    "\n",
    "from cutcutcodec.core.io import read\n",
    "from cutcutcodec.core.opti.parallel import starmap\n",
    "from cutcutcodec.core.signal.metric import psnr as compute_psnr, ssim as compute_ssim  # the main functions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5c31cef-e0e1-4ced-8aa1-8e5e6c3c9b26",
   "metadata": {},
   "source": [
    "## Preparing videos to be compared\n",
    "Here we'll compare the original video with a noisy version."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9c2dab0a-8e62-4d65-b86f-bd15fed6426c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/data/Documents/cutcutcodec_git/media/video/intro.webm at 30000/1001 fps, of size 1280x720\n"
     ]
    }
   ],
   "source": [
    "from cutcutcodec.core.analysis.stream.rate_video import optimal_rate_video\n",
    "from cutcutcodec.core.analysis.stream.shape import optimal_shape_video\n",
    "from cutcutcodec.utils import get_project_root\n",
    "\n",
    "file = get_project_root().parent / \"media\" / \"video\" / \"intro.webm\"\n",
    "stream_ref = read(file).out_select(\"video\")[0]\n",
    "rate = optimal_rate_video(stream_ref)\n",
    "shape = optimal_shape_video(stream_ref)\n",
    "print(f\"{file} at {rate} fps, of size {shape[1]}x{shape[0]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "85eacf9c-3256-4437-9e27-84bc94b8cb91",
   "metadata": {},
   "outputs": [],
   "source": [
    "from cutcutcodec.core.filter.video.equation import FilterVideoEquation\n",
    "from cutcutcodec.core.generation.video.noise import GeneratorVideoNoise\n",
    "\n",
    "noise = GeneratorVideoNoise(0).out_streams[0]\n",
    "stream_comp = FilterVideoEquation([stream_ref, noise], \"(9*b0+b1)/10\", \"(9*g0+g1)/10\", \"(9*r0+r1)/10\").out_streams[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a45ac754-6dc3-42be-90a4-a76f7a554750",
   "metadata": {},
   "source": [
    "## Compare each frame"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a33c4311-3cab-4291-9fe9-949c47e6e5b4",
   "metadata": {},
   "source": [
    "### Convert all streams into YUV in the same colorspace\n",
    "Image comparisons are always made in yuv, not rgb."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8e7939ef-ab7e-4343-8b14-ec971c74169c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from cutcutcodec.core.filter.video.colorspace import FilterVideoColorspace\n",
    "from cutcutcodec.config import Config  # to get the default used colorspace\n",
    "\n",
    "dst_color = f\"y'pbpr_{Config().target_trc}_{Config().target_prim}\"\n",
    "ref_yuv = FilterVideoColorspace([stream_ref], dst=dst_color).out_streams[0]\n",
    "comp_yuv = FilterVideoColorspace([stream_comp], dst=dst_color).out_streams[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "88798d3c-8fcf-4c13-8ac3-2bdc4a58ac63",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compare(stream1, stream2, time, shape):\n",
    "    \"\"\"Compare 2 images.\"\"\"\n",
    "    frame1, frame2 = stream1.snapshot(time, shape), stream2.snapshot(time, shape)\n",
    "    frame1, frame2 = frame1.numpy(force=True), frame2.numpy(force=True)\n",
    "    psnr = compute_psnr(frame1, frame2, weights=(6, 1, 1))  # factor 6 on y\n",
    "    ssim = compute_ssim(frame1, frame2, weights=(6, 1, 1), data_range=1.0)\n",
    "    return psnr, ssim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d715ff1e-3c0d-4586-9e88-162310007efc",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 294/294 [00:46<00:00,  6.31it/s]\n"
     ]
    }
   ],
   "source": [
    "times = np.arange(0, ref_yuv.duration, 1/rate).tolist()\n",
    "all_psnr, all_ssim = zip(*starmap(compare, ((ref_yuv, comp_yuv, t, shape) for t in tqdm.tqdm(times))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "fbf03e40-a89e-44b9-87a1-dd562214a845",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAASi5JREFUeJzt3Xd4lfX9//HnGTkne5FFSELYeyMRUasVV611i0pFscO28K1KtUoVqFqltnUPUL9V+/OrgraiVhEHioqizCB7QwLZZO/knPv3x0kORAKEcMidc/J6XFeu5pzc55z3SZHz4v1ZFsMwDEREREQChNXsAkRERER8SeFGREREAorCjYiIiAQUhRsREREJKAo3IiIiElAUbkRERCSgKNyIiIhIQLGbXUBHc7vd5OTkEBERgcViMbscERERaQPDMKioqCA5ORmr9di9mS4XbnJyckhNTTW7DBEREWmH7OxsUlJSjnlNlws3ERERgOeXExkZaXI1IiIi0hbl5eWkpqZ6P8ePpcuFm+ahqMjISIUbERERP9OWKSWaUCwiIiIBReFGREREAorCjYiIiAQUhRsREREJKAo3IiIiElAUbkRERCSgKNyIiIhIQFG4ERERkYCicCMiIiIBReFGREREAorCjYiIiAQUhRsREREJKAo3IiIictLqG90UVNRS3+g2u5Sudyq4iIiI+MYrX+/h4835lNc2sD2/0htsxvSM4T+/PcO0uhRuRERE5ISt2HWQP/93c6s/c9rNHRhSuBEREZETUlPvYubb3wPw0+HduWxkD/olhJMaG0p5TQP1LnOHphRuRERE5IS8+u1e9h6sJikymLlXDiMiOMj7s5gwh4mVeWhCsYiIiLSZYRi8uXo/ALdN7Nci2HQWCjciIiLSZuv3l7GzoJLgICuXDO9udjmt0rCUiIiIHFeDy83KPcW8/l0WABcOSSKyE3ZtQOFGRERE2uCVr/fy0OIt3ttXjU4xsZpjU7gRERGR49qSVw5AbJiDCX3jmNA3zuSKjk7hRkRERI6rsKIOgD/9ZBBXj+m8XRvQhGIRERFpg4JyT7hJiHCaXMnxKdyIiIjIcRVWesJNvMKNiIiI+Lv6RjfFVfWAOjciIiISAIqaujZ2q4WYUPN3ID4ehRsRERE5poKKQ0NSVqvF5GqOT+FGREREjqmgvBbwjyEpULgRERGR4zjUuQk2uZK2UbgRERGRY2re4yYhUp0bERERCQDNnRsNS4mIiEhAKKxonnOjYSkREREJAOrciIiISEDxHr2gOTciIiLi79xuw7uJnz8cvQAKNyIiInIMb63JptFt4LBZiQv3j3BjN7sAERER6VwMw+CvH25lc2453+w6CMBvzulDkM0/eiIKNyIiItLC5txynv9yt/f25SOTuWNiPxMrOjEKNyIiItLC3qJqAPrEh3HfTwdzdr94LJbOf6ZUM4UbERERaSGr2BNuhqdEc+6ABJOrOXH+MXgmIiIiHSaruAqA1NhQkytpH4UbERERaWHfQU/npqfCjYiIiAQCb7jppnAjIiIifq6+0U1uWQ0AaQo37fPss8+Snp5OcHAwGRkZrFy58pjXP/HEEwwYMICQkBBSU1O54447qK2t7aBqRUREAtv+kmrcBoQ6bMT7yaZ9P2RquFm4cCEzZsxgzpw5rF27lhEjRnDhhRdSUFDQ6vWvv/4699xzD3PmzGHLli3885//ZOHChfzpT3/q4MpFREQCU/NKqbTYUL9a/n04U8PNY489xq9+9SumTp3K4MGDmT9/PqGhobz00kutXv/NN98wYcIEbrjhBtLT07ngggu4/vrrj9vtERERkbZpDjf+ulIKTAw39fX1rFmzhokTJx4qxmpl4sSJrFixotXHnHHGGaxZs8YbZnbv3s3ixYv5yU9+ctTXqauro7y8vMWXiIiItGQYBm+szOK/63MA/10pBSZu4ldUVITL5SIxMbHF/YmJiWzdurXVx9xwww0UFRVx5plnYhgGjY2N/OY3vznmsNTcuXO5//77fVq7iIhIoFm5p5iZb2/w3k6PCzOxmpNj+oTiE7Fs2TIefvhhnnvuOdauXcvbb7/NBx98wIMPPnjUx8ycOZOysjLvV3Z2dgdWLCIi4h/WZJUAnuXf149L5Wcjk02uqP1M69zExcVhs9nIz89vcX9+fj5JSUmtPmbWrFnceOON/PKXvwRg2LBhVFVV8etf/5p7770Xq/XIrOZ0OnE6/XO2t4iISEdZn10KwM8zevKrs3ubW8xJMq1z43A4GDNmDEuXLvXe53a7Wbp0KePHj2/1MdXV1UcEGJvNBnjGCkVERKR91meXATAiNdrcQnzA1IMzZ8yYwU033cTYsWMZN24cTzzxBFVVVUydOhWAKVOm0KNHD+bOnQvApZdeymOPPcaoUaPIyMhg586dzJo1i0svvdQbckREROTE5JXVkldei9UCQ3tEml3OSTM13EyaNInCwkJmz55NXl4eI0eOZMmSJd5JxllZWS06Nffddx8Wi4X77ruPAwcOEB8fz6WXXspDDz1k1lsQERHxe+v3lwLQPzGCUIep0cAnLEYXG88pLy8nKiqKsrIyIiP9P52KiIicrL9/tJVnP9/FpLGpPHL1cLPLadWJfH771WopERER8b3v93vm2wxPjTK5Et9QuBEREeni8so8ZzSmd/PfvW0Op3AjIiLSxZXWNAAQHRpkciW+oXAjIiLShRmGQWl1PQDRoQ6Tq/ENhRsREZEurLreRYPLs7YoRp0bERER8XclTV0bh81KSFBg7BmncCMiItKFlVYfmm9jsVhMrsY3FG5ERES6sOZwExMg821A4UZERKRLax6WigqQ+TagcCMiItKlNS8DD5TJxKBwIyIi0qWVVjUtAw/RsJSIiIgEAO8GfmHq3IiIiEgAaJ5zownFIiIiEhC8S8FD1LkRERGRABBoRy+Awo2IiEiXdvgmfoFC4UZERKQLO7QUXJ0bERER8XNu96ETwbXPjYiIiPi9itpG3J4DwbVDsYiIiPi/0hpP1ybUYcNpD4wTwUHhRkREpEsyDIP9JTVAYM23AbCbXYCIiIh0vD+8tZ631x4AICqA9rgBdW5ERES6pBW7Dnq/P3dgvImV+J46NyIiIl1QRW0jAEtuP4uBSZEmV+Nb6tyIiIh0MW63QWWdJ9x0C3OaXI3vKdyIiIh0MVX1jd7vI4IDbxBH4UZERKSLae7aBNksOO2BFwUC7x2JiIjIMTXPtwl32rFYLCZX43sKNyIiIl2MN9wE4JAUKNyIiIh0Oc3DUhHOwNrfppnCjYiISBdTUes5CVydGxEREQkIlbXNnRuFGxEREQkAzcNS6tyIiIhIQGieUByIe9yAwo2IiEiXc2gpuCYUi4iISACorPNMKFbnRkRERAKCdym4wo2IiIgEgsN3KA5ECjciIiJdjMKNiIiIBBQtBRcREZGA0ryJX2SwVkuJiIhIAPAev6BhKREREfF3LrdBVb0L0LCUiIiIBICq+kbv9+rciIiIiN9rnm/jsFkJDrKZXM2poXAjIiLShXiXgQfokBQo3IiIiHQpzUcvBOqQFCjciIiIdCmBfiI4KNyIiIh0KYG+OzEo3IiIiHQp5bWBfSI4KNyIiIh0KfnldQAkRAabXMmpo3AjIiLSheSV1QDQXeFGREREAkFuWS0A3aNDTK7k1FG4ERER6UK84SZKnRsREREJAHlN4SZJ4UZERET8XUVtA5V1nqXg6tyIiIiI32sekooKCSLUoaXgIiIi4ue6wnwbULgRERHpMpqXgQfyfBtQuBEREekyDnVuAncZOCjciIiIdBm5pRqWEhERkQCSWx74y8BB4UZERMSvLdmYy3UvrKCwou6o19Q2uHhjZRZbc8sBSNawlIiIiHRWv/m/tXy7u5h5y3Yd9ZoFK7OY+fYGCpoCUFpsaEeVZ4rAXeQuIiIS4LKLq73fWy1Hv25bfiUAo9OiuW5cGmndFG5ERESkE/pyR6H3+waX+6jX7S/xhKDrx6VxzdjUU16X2TQsJSIi4qe+2HYo3BRV1R/1uuYOT2qAD0c1U+dGRETEz+w7WMXHm/L5emeR977iytbDjdttcKDUs3mfwo2IiIh0Svcu2sjyw4INwMGq1ldL5VfU0uAysFstJEUG9hLwZhqWEhER8TPNnZgJfbvxh/P7A1B8lGGp7GLPtcnRIdiONes4gJgebp599lnS09MJDg4mIyODlStXHvP60tJSpk2bRvfu3XE6nfTv35/Fixd3ULUiIiLmK69pAGDWTwdz7WmeCcLFVfW43cYR1x6abxPYe9scztRhqYULFzJjxgzmz59PRkYGTzzxBBdeeCHbtm0jISHhiOvr6+s5//zzSUhI4N///jc9evRg3759REdHd3zxIiIiJjAMg7KmcBMVEkRMqAMAtwGlNQ3EhjlaXJ/dtFIqNaZrzLcBk8PNY489xq9+9SumTp0KwPz58/nggw946aWXuOeee464/qWXXqK4uJhvvvmGoKAgANLT0zuyZBEREVNV17tobOrQRIUE4bBbiQy2U17bSHFV3ZHhpmlYKiWm63RuTBuWqq+vZ82aNUycOPFQMVYrEydOZMWKFa0+5r333mP8+PFMmzaNxMREhg4dysMPP4zL5Trq69TV1VFeXt7iS0RExF81d22CbBZCgmwAdAt3AnCwlRVT3s5NF1kpBSaGm6KiIlwuF4mJiS3uT0xMJC8vr9XH7N69m3//+9+4XC4WL17MrFmzePTRR/nLX/5y1NeZO3cuUVFR3q/U1MDfvEhERAJXea0n3EQGB2GxeCYIN3drDh42qbiqrpE3V2Wzs8CzO3FKFxqWMn1C8Ylwu90kJCTwwgsvMGbMGCZNmsS9997L/Pnzj/qYmTNnUlZW5v3Kzs7uwIpFRER8q6z60HybZt1aCTcP/Hczf/zP995VVD0D/MiFw5k25yYuLg6bzUZ+fn6L+/Pz80lKSmr1Md27dycoKAibzea9b9CgQeTl5VFfX4/D4TjiMU6nE6fT6dviRURETNI8LBV5eLgJbwo3lZ69bvYUVfHvtfsB+MmwJMb3iSMuvOt8FprWuXE4HIwZM4alS5d673O73SxdupTx48e3+pgJEyawc+dO3O5D52ds376d7t27txpsREREAk15bSPww86NJ7g0d2me/HQ7LrfBjwcm8NzkMdx4es+OL9REpg5LzZgxgxdffJF//etfbNmyhd/+9rdUVVV5V09NmTKFmTNneq//7W9/S3FxMbfddhvbt2/ngw8+4OGHH2batGlmvQUREZEO1Vrn5vA5NzvyK3h3fQ4AM5o2+OtqTF0KPmnSJAoLC5k9ezZ5eXmMHDmSJUuWeCcZZ2VlYbUeyl+pqal89NFH3HHHHQwfPpwePXpw2223cffdd5v1FkRERDrUoT1uDn2EHz4s9fin2zEMuGhIEkN7RJlSo9lMP1tq+vTpTJ8+vdWfLVu27Ij7xo8fz7fffnuKqxIREemcymtam1DsGZbalFNORW0jFgvc0UW7NuBnq6VERES6utbCTWKkJ9xUNM3HuXR4MgOSIjq+uE7C9M6NiIiItJ13zk3woXDTNyGcOy/oz9a8CoKDbNx14QCzyusUFG5ERET8SFkrnRuLxcL0H/czq6ROR8NSIiIifqS1cCMtKdyIiIj4Ee/xCwo3R6VwIyIi4kfUuTk+zbkREZFOb2dBJX9+bxO1DS6Gp0Qz66eDvIdGnmrV9Y2s3FNMRW0jY3rGkBwd0iGv25q6Rhe1DZ5d+tW5OTqFGxER6fT+79t9LN9ZBMDqfSVcObpHh21Qd9uCTD7Z7DkHcUBiBB/dcbZPn98wDJ75bCd7D1aTFhvKTWf0JDq09SOFmrs2FgtEOPURfjT6zYiISKeXmV3a4vbarJIOCTebcsq8wQZgR0EFjS43dpvvZnW89l0Wj36y3Xu70e3mDxe0vpS7vMazj01kcBBWa8d0rvyR5tyIiEinVt/oZnNuOQBXjOoBwNp9JR3y2vOW7QLgkuHdCbJZcBuQX1Hns+ffW1TFQx9sAWBg06Z7W/Mqjnr9oXOl1Js4FoUbERHp1LbnV1Df6CYqJIjLRiYDsO4HnZxTIbu4msUbcgGYdk5fEiODAcgtrfFe0+hy89v/W8PYv3zCuIc+5eNNeW1+fsMwuO+djdQ0uBjfuxt/+skgAHYXVh71Ma3tTixHUrgREZFOrXlIanhKFKPSYgDYd7CaokrfdVBa882uItwGjO0Zw+DkSJKjPBOJc8pqvdc8t2wXH27Mo6iynoKKOl79dl+bn/+zrQUs31mEw2blr1cNo09COABZxdU0utytPmbF7oMA9OwW1t631SUo3IiISKf2/f5SAEakRBMVEkS/phCwLqv0lL7u+v1lAIzp6QlU3aNbdm425ZTx1NIdAPz67N4ArN5bQn1j68Gk2Ueb8rj82a+569/fA3DLmb3o2S2M7pHBBAdZaXAZZJfUHPG4RpebResOAHDZiOSTfXsBTeFGREQ6pbyyWuZ/sYvlOzyrpIaneCYQj27q3qzLOrXzbryhKjUagO5NnZvcps7N/32bRaPb4MIhidxz0UBiwxzUNLjYcKD0mM/79Gc7yMwupbiqnrhwJ9PO7QOA1WqhV5wnuLU2NPX1roMUVtQRExrEOQMSfPAOA5dmJImISKezPb+CG//5Hfnlh4aemkPG8NQoFq7O9k4yPhVqG1xszfVM7G0OVclNnZucps7NxgOezs7lI3tgtVrI6BXLhxvz+HZ3MWN6xgKeeTWLN+Sx92AVkcF2rhydwpam55175TAm9Ikj4rADMHvHh7Elt5xdhZWcNygRgD1FVfzx3+vZe7AagJ+NSMZhV2/iWBRuRESkUymsqGPS8ysoqW6gd1wYQ3pEMTI12juht3dTd2NvUdUpq2FzbjmNboO4cAc9mjbtO7xzU9/oZlvTqqbmJemn9+7WFG4OMu3cvgB8vq2Aaa+v9T7vjoJKXG6DxEgn152WesRGhH3iPHNpdhd63pvbbXDXW+tZ3bQ6zGKBa8amnqq3HTAUbkREpFN5/NPtlFQ3MDApggW/Pv2IDe3S40IB2F9SQ4PLTZAP95xp9r13EnO0N4B0j2qac1NWw46CCupdbiKD7aTEeELP6b27AZ55N811rdrrCSU2qwWX2+C177IAz9Baazss9473BLedBZWs2VfMZ1sLWL2vhFCHjUevGUFqbGiHbV7ozxRuRESkU9h3sIpVe0tYsNITAB64bGirO/UmRngm3tY2uDlQUkN6XOsrhxauyqJntzBv6GiL2gYXf/9oG59vLQA8k5ibNR+7UFRZ791nZ0hylDek9EsIJyY0iJLqBr7fX8aYnjHeoatJp6Xy+ndZuNwGAKPSDj3v4XrHe97L6n0lXDVvhff+P1wwgIuHdW/z++jqNGgnIiKmKyiv5bxHv+DOt9bjNuDioUmM6xXb6rVWq4WesZ4QsPdg60NTGw+Ucfd/NvDLf62mpt7V5jrezTzAP5fvYXfTkNdpvWK8P4sJDSI4yPOx+ckWT/gZ2iOyRV0ZvTxB6tvdBzEMgw1N4ea601JJiw31Xts8KfqHeseH07zxcJjDxrj0WH5/Xj9uPiO9ze9B1LkREZFOYH9pDY1uA6fdysRBidz300HHvD49LpRt+RWeeTetnFTQvDdOZV0jH2/O47KRPdpUx3e7iwGYOCiRq8f0YPxhXR+LxUJyVAi7i6r4cnshwBFDRKf3jmXJJs+8m5+NSKa0uoEgm4UBSRFcPCyJ57/YTZDNctShpXCnnQcvH0pRRf0xz5iSY1PnRkRETFfb4OmupMWG8uzk0d7Ju0fTPBTVvILohzbllHm//8/aA22u47s9nnBz8xnpXDS0+xHzYpr3umk2JPkH4aaPJwyt2Vfi3UW5f2IETruNq0an4LRb+VH/BIKDbEetYXJGT26b2E/B5iSocyMiIqara9r4zhnUtn9zp3c73rDUoWXiy3cU8uSnO+gRE8JZ/eIwDM/xBSGOlgFjf0k1B0prsFstjO4Z3erzDkmO4uudnl2C4yOc9PrBfJ/+CRHeeTdvNE0eHtbUpemfGMHyu39MRLA+ek81/YZFRMR0dQ2ecBNsP3pH43DecNPKcvDDl2mnxoaQXVzD459ub3FNUmQwi287i9iwQ92RlU1dm6E9ogh1tP7x+IcL+jOhbxzVdY0M7RGF7QcnczfPu1myKc97VMLhQ1DxEc42vT85ORqWEhER09U1eoal2tq5ae6YZJfUsCO/osVZTIcv0543eQxXjurBNWNSWkz+zSuv5YH/bmrxnM3zbTJ6tz6RGcBpt/Gj/vFcPKw7qYdNED7chH5xLW6fln7055NTQ50bERExXXPnxtnGzk1ChNO7HPz8x7/kmjEp/P2aEQBsyvEMSQ1JjmJojygemzTy0Os0utiSW8GVz33NO5k5/GRYdy4YkkRFbQOfbfOsgMo4yiqttpo0NhUMg3qXwcCkCAYkRZzU88mJU+dGxE/8vxV7ufDxL7nw8S95N7PtEyRF/EFtU+cmuI2dG6vVwiXDDh0e+cGGXIqr6vnda2t48lPPYZaHd2qaOe02RqZG88uzPAddznhzPVvzypn74VYKK+pIiw3ljD5xRzzuRDjsVm4cn84vzuzFhL4n91zSPgo3In7ixa92sy2/gm35Ffzrm71mlyPiUyfauQF49NoR7HzoYhIjnVTXu7h9YSaLN+RxoOnsp+bznVpz5wUDyOgVS2VdIxc/+RWvN03+feSq4cdcyST+QeFGxE/UNhyaU1Be22hiJSK+17wUvK2dm2Z2m5UfD/SckN2898yU8T2Z//PRXDA48aiPc9itPH/jGAZ3j8TwbBrMLRN6Mb5P23czls5Lc25E/ER946FwU1HbYGIlIr7nXQp+Ap2bZj8emMgbK7MBzyZ4f7xoIOHO43+8RYc6eP9/ziS/opYgm5VuYdpXJlAo3Ij4iQbX4eFGnRsJLCe6WupwE/p2w2G3Ut/o5vJRyW0KNs2sVstxNwwU/6NhKRE/cXjnprre1WLpq4i/q23HnJtmoQ47V41OITLYztQJvXxdmvghdW5E/IDbbdDYdJpws8q6Rm3PLgHD27mxt+/f3HOvHMZDlw/F+oNN9aRrUudGxA/Ut9Kl0dCUBJLmzs3JrFRSsJFm7e7c5OTksHz5cgoKCnC7W/7F+/vf//6kCxORQ+oOG5KKDg2itLqBck0qlgBysp0bkcO1K9y88sor3HrrrTgcDrp169bi1FSLxaJwI+Jjh08mjg11UFrdoM6NBBRfdG5EmrUr3MyaNYvZs2czc+ZMrFalbJFTrXkyscNmJSIkCNCwlAQWdW7El9r1p6i6uprrrrtOwUakg3jDjd1KRNMyV+11I4FEnRvxpXalk1/84he89dZbvq5FRI6ieUKxw24lIrg53KhzI4Hj0CZ++keznLx2DUvNnTuXn/70pyxZsoRhw4YRFBTU4uePPfaYT4oTEY8Ww1LB6txI4NGwlPhSu8PNRx99xIABAwCOmFAsIr7V3LkJsluICNacGwk8dRqWEh9qV7h59NFHeemll7j55pt9XI6ItKa1zo0Oz5RAcjLHL4j8ULv+FDmdTiZMmODrWkTkKA5NKLZ5OzeVdQo3Eji8E4rbcfyCyA+1K9zcdtttPP30076uRUSOosVqKc25kQCkzo34UruGpVauXMlnn33G+++/z5AhQ46YUPz222/7pDgR8Wiec+O0WYnUaikJMC63QYPLc3aaOjfiC+0KN9HR0Vx55ZW+rkVEjqK5c9NyQrE6NxIYmrs2oM6N+Ea7ws3LL7/s6zpE5Bi8+9zYtM+NBJ7m+TYATnVuxAcUkUX8QMs5N1oKLoGluXMTZLNg08ne4gNt7tyMHj2apUuXEhMTw6hRo465n83atWt9UpyIeLRcLeX5z7ayrhGX29CHgfi95j1u1LURX2lzuLnssstwOp0AXH755aeqHhFpRWvDUuAJOFEhQUd7mIhfqG3q3ARrvo34SJvDzZw5c1r9XkROvUOdGwtOuw2H3Up9o5uK2gaFG/F76tyIr7UrJmdnZ7N//37v7ZUrV3L77bfzwgsv+KwwETmk4bDODaDl4BJQahu0x434Vrv+JN1www18/vnnAOTl5TFx4kRWrlzJvffeywMPPODTAkWk5YRiQJOKJaAcOhFcnRvxjXaFm40bNzJu3DgA3nzzTYYNG8Y333zDa6+9xiuvvOLL+kSEQ3/5N4eb5s5NaXW9aTWJ+Epz50ZzbsRX2vUnqaGhwTu5+NNPP+VnP/sZAAMHDiQ3N9d31YkIcPiEYs+/bLtHhQBwoLTGtJpEDldUWcfSLfl88H1ui0352uJQ50bhRnyjXZv4DRkyhPnz53PJJZfwySef8OCDDwKQk5NDt27dfFqgiLTcoRggNdYTbrKLFW7kkJKqeg5W1dM3IbxDX3dLbjnXzF/hPcz19+f1Y8b5/dv8eA1Lia+1KyY/8sgjPP/885xzzjlcf/31jBgxAoD33nvPO1wlIr7zwwnFqbGhAGSXVJtWk3QuuworOe+xL7jg8S/4ZldRu5/n0Y+3ceurq3l77X4+2ZzP2qySY15fXd/I9NfXUlnXSGyYA4AFK7NodB3adTjrYDW3vrqar3e2XpeGpcTX2tW5ycjIYP/+/dTX1xMTE8O+fftYtGgRPXv25I477vB1jSJdXv0P2vapMU3hpljhxiyl1fVszi0nzGEnPNhOhNPzv1aLhX0Hq3n/+xwqahuJDg2iuKqeYT2iuHpMyjE3QD3cwco6Xvsui5AgGxcNTaLB5ebrnUWszSolNsxBuNNOVV0jew9WUdPgYmdBJcVVnjlY9/xnA2//7gzCHHZCHEfvhry3PocXvtzFnRcM4JwBCXy1o5CnP9sJwEeb8r3XPTFpJFtyy/nv+hzcBoxJj2FyRhpOu42HPtjMrsIqEiOdvDvtTC556isKKur4fFsh5w9OxO02mPFmJqv3lbDxQDmf33mOd+5YM3VuxNfaFW4uu+wyrrzySn7zm99QWlrKuHHjcDgcFBUV4Xa7+e1vf+vrOkW6tB+ulmoeltpfUoNhGG3+wBTfqG1wcc38FewoqDyhx23JrWDWTwcd9f+vgopa7liYSXW9i92FVZTVeA5HfWjxljY9f6+4MOoaXGQVVzP2L58SEmRj0bQzGJgU2eK6ukYXK3YdZMbCTBrdBr97bS3/75ZxPLx4KwCnpcdQ2+Cmur6RXYVV3L4ws8XjP/g+lw++PzS/Msxh4+nrR5MUFcxVY1J44cvdzFu2k8q6BrbkVrB6n6f7c6C0hrfWZDM5o2eL51PnRnytXeFm7dq1PP744wD8+9//JikpiXXr1vGf//yH2bNnK9yI+Jh3QnFTuElp6txU1jVSWt1ATNNwgHSMxz/Zzo6CSkIdNmJCHZTXNlBZ14hheH4eZLNwzoAEeseHUd4UUN5Ymc1LX+/BGWTl7osGtvq8c97dxNc7D3pvD0yKIMRhY11WKaEOGwOSIvhR/3gqaxupaXAREmQjrVsoEcF2ahvcXDgkiY0HyrjllVU0ug1qGlw8vHgr/++WQ9MFnlu2k8c+3k6j21NshNNORV0jV89f4bkdbOeFG8cSE+bA5Ta49dXVfLqlAIfdysNXDCM1JoQ3V+/n651F1DS4GNcrlj//bAg9oj2Be9Jpqbzw5W7WZpWyNqvU+7qj0qJZl1XK3z/axpKNeVw1OoXLR/UA1LkR32tXuKmuriYiIgKAjz/+mCuvvBKr1crpp5/Ovn37fFqgiBz6yz+oac5NcJCN+AgnhRV17C+pUbg5xVbsOshLX+8ht6yG4sp6cstrAXjyulGcPzgRAMMwqKpv2ozObvX+f9VsREo097y9gXnLdhETGsSQ5CjG9IwhOMjzgf7hhlw+3JiH3WrhL5cPpVu4k3MHxGO3WWl0ubHb2tbVOLt/POtmn8+B0houfXo5X24v5ONNeZw/OJEvthfytyXbvNeePziRh68Yxu0L1/HNLk+ouvuigd4/TzarhSevG8X/fbuP8X26MTwlGoCM3kdfONInPpxHrhrGVzuKKKmux4KFAUkR/OGC/pz36BfkltXy1Y4ivtpRxDe7ikiNCeX7/aWAOjfiO+0KN3379uWdd97hiiuu4KOPPvLOsykoKCAyMvI4jxaRE/XDCcUAKTEhFFbUkV1SzbCUKLNKC1ir9hbz8OItlFU3sLuo6oif3zS+pzfYAFgsFsKdR/8r9bpxaewrrmbesl3e4Z+02FCmn9uXXYWV/O/yPQD8+uzeXDcurcVj2xpsmkUEBzEwKYifn96Tl7/ey69fXUNksJ2apuGfn5+exgM/G4q16dDV1355Om63QYPbfUT3JMxp59Yf9Tmh1590WhqTTks74v43bx3P6n3FbMmt4IUvd/Pm6v0tfq7OjfhKu8LN7NmzueGGG7jjjjs477zzGD9+PODp4owaNcqnBYrIkXNuwDOpeF1WqSYVnwI5pTX8+v+tpqTaM6RktcDkjJ6cMyCebuFO4iOc3mGYE3HXBQNwGwZfbCuksKKOrOJq/vif770/v2R4d35/Xj+fvY/bz+vP/pIavthWSHnTbtZjesZw3yWDvcGmmdVqwWk9teEiNTaU1NhQrhgF49Jj+WhTHu9m5niHXdW5EV9pV7i5+uqrOfPMM8nNzfUuAwc477zzuOKKK3xWnIh4tBpumve60XLwk5ZdXM0D728mt6yGcKedfQerKaluYEhyJPdeMojUmFDv8vuTYbVamHnxIGZePIiK2gaeWrqDzGzP6qcrRvXgoqHdffBuDokKDeLFKWOpqfdMMrbbLPTqFnZEsDHDxMGJTBycSFGlZ2UVqHMjvtOucAOQlJREUlJSi/u0x43IqdH8L1vnDzo3AMt3FPHMZzv45Vm9vfM3fKnR5cZlGDhs1oBalVVd38i9izZSWFHH9/tLvZ2NZlEhQTw3eTQ9u4WdktePCA7i3ksGn5Ln/qGQpsnIndH5g5MOhRt1bsRH2h1uRKTj1P9gQjFA73jPLrR7D1bzj4+347Tb+NXZvX36ui8t38MjS7ZS1+hmZGo0C289Hafdhstt0NjK/Ax/8vRnO1m07oD39sjUaH53Th9qG92EOWyMSovxbkonp87EwQn8aZHn+9KmYUCRk6VwI+IHGlxHDkudlh7DrJ8O5rvdB/l4cz5vrs7ml2f1OmZ35UBpDRv2lxEZbGdMegxlNQ08/MEWsoqraXAZNLjc1LvchDnsnDMgnmc+3+ld3pyZXcpfP9zKppxyVu4pBmDmxQNPeLJpZ7CrsJL//Wo3ADPO70/fhHAmDko8YnM5OfUSIoK936f5YOhPBBRuRPyC91Twwzo3FouFX5zZi2vGpjDuoU/ZUVDJ+v1ljEyNbvU5Glxurp73DbllnmXMg7tHEmSzsH5/WavXbzjguX/S2FTG9YrlD2+t5+Wv97a45u8fbePMfnEMSfas1qpvdLNmXwnZxdUE2S1cNqJHp5jfcTjDMPjze5tocBn8eGCCTyfwSvssu/Mcvtl1kJ8M8+2cI+m6FG5E/EBrE4qbRQYHcdGQJN7JzOHRj7cxoW8c3+0+SO/4cO675NBuuJ9tLSC3rJZQhw271cLm3HIAokODePCyoYQH2wmyWgmyWVi9r4TnPt9J/6QI7r9sCE67lXcyD/DVjiJSYkL435vG8sQnO1iyKY9bXllF/8QIaupdbM+vaDF3Ja+sjt+e4/vOTnV9I/O/2E1hRS39EyO4dmwqYU47hmHw6rf7qK538euzercarD7cmMdXO4pw2K3MubRj5rzIsaXHhZEed2rmNknXpHAj0skZhnHEDsU/dO3YVN7JzPFujgbw+bZC0mJDuemMdADeatpT5MbTe/Lz03tyyyur2F9Sw7zJYxjfp+WmbBm9u3Hr2b2xWCzYmgLCMzeM5r31OVw0JIn4CCd/uWIoq/YWk19eR355nfexceFOesWFsmpvCY9/sp2JgxLol3hyk1kXrMxi3he7KKyoI71bGA0ud4ujD97NzOGVqafxwpe7eW7ZLsCzfPvXZ3uCVYPLTWZ2KeuzS/nfrzz7yfzm7N6nbLKwiJjLYhjNI+pdQ3l5OVFRUZSVlWnDQfELDS43/e79EIDM2ecTHXrkJFfDMHjlm71sOFBGbYOLkCA7/1m7H4fdyoQ+3YgJdfDu+hxcboNPZ5xN34QI3G6DyvpGIoOD2l1bdnE167JLaXS5CQmykRgVzIiUaKwWuOWVVXy+rZA+8WH8v19kEBFsJyTIdsTOvcdTXd9IxsNLqfjBaqa4cCfXjE3hjZVZrU5EtVktXDs2laLKOlbsOkhl3aHHp8SE8MkdPzrmoZIi0rmcyOe3OjcinVzzZGI4eufGYrEwdUIv723DMCipruezrQXeZbbgWRHUN8HTRbFaLScVbODQpmytmXvlcC5/9mt2FVbx438so67RzbAeUbz9uzOOGXBKq+uZt2wXOwoqSYhw0isujIraRtJiQ3np5tNYtbeYHfmV/Prs3iRFBfOzEclMfXkVeeW1OOxW7rloIJnZpby3Poc3VmZ5nzcmNIhxvWIZkBjB9RlpCjYiAaxThJtnn32Wv//97+Tl5TFixAiefvrpNu2Zs2DBAq6//nouu+wy3nnnnVNfqIgJmufbQMsJxcdisVh46vpRfPB9DhYs7CqqZHNOOf/z446bPJsUFczbvzuDm19eyfZ8zxDShgNlvLpiH1MnpFNZ14jLbbToRNXUu5j6yirWHXbgYrPJGWn0TQinb0J4i/sHdY9k+d3nUlLdQFRIEA67lbpGF2f2jWN/aQ2hDhsT+sQxJDmy001uFpFTw/RhqYULFzJlyhTmz59PRkYGTzzxBG+99Rbbtm0jISHhqI/bu3cvZ555Jr179yY2NrbN4UbDUuJvCsprGffwUqwW2D33ErPLOWHV9Y1s2F/GhgNl/OWDLYQ5bIQ47BRVeubpnN0/npvG96S+0c3zX+4mM7uUqJAgbjuvH08u3UFZTQMOu5VvZ56nfWdEurAT+fw2fVOHxx57jF/96ldMnTqVwYMHM3/+fEJDQ3nppZeO+hiXy8XkyZO5//776d3bt5uWiXQ2dcdYKeUPQh12Mnp3Y+qEXgzuHklVvcsbbAC+3F7IL/61mt++tpbM7FJCgmy8dPNYbjmzFy9PPY2UmBBuPbu3go2ItJmpw1L19fWsWbOGmTNneu+zWq1MnDiRFStWHPVxDzzwAAkJCfziF7/gq6++OuZr1NXVUVd36C/S8vLyky9cpAPVt3IiuD+yWS3M+/loFqzKZlyvWE7v1Y2CilqeWrqTTTlluNwG5wyI5+YJvbyHUo5Oi2H53T82uXIR8TemhpuioiJcLheJiYkt7k9MTGTr1q2tPmb58uX885//JDMzs02vMXfuXO6///6TLVXENK3tTuyvenYL4+6LBra4/ei1I47xCBGRE+dXf1tWVFRw44038uKLLxIXF9emx8ycOZOysjLvV3Z29imuUsS36lvZnVhERI7O1M5NXFwcNpuN/Pz8Fvfn5+cfceI4wK5du9i7dy+XXnqp9z632/MXv91uZ9u2bfTp03I3VKfTidPpPAXVi3SMY+1OLCIiRzL1b0uHw8GYMWNYunSp9z63283SpUsZP378EdcPHDiQDRs2kJmZ6f362c9+xrnnnktmZiapqakdWb5Ih1C4ERE5MabvczNjxgxuuukmxo4dy7hx43jiiSeoqqpi6tSpAEyZMoUePXowd+5cgoODGTp0aIvHR0dHAxxxv0igqAugOTciIh3B9HAzadIkCgsLmT17Nnl5eYwcOZIlS5Z4JxlnZWVhteovdem6Gpo6Nyd6bIGISFdl+iZ+HU2b+Im/ef/7HKa/vo6MXrEsvPXI4VoRka5AZ0uJ+KGy6gYWb8ylvtFNg8tNRW0jDruVlXuKAQ1LiYi0lcKNiAncboNPtuSTW1pDiMNGfIST+/+7mX0Hq4/6GKfCjYhImyjciJxChmGwYtdBokMdDEyKwAC+3FHI459s5/v9ZUdcnxwVzKi0GOw2C+FOO/WNbtZmlbCrsIqRqdEdXr+IiD9SuBHxAZfbYF1WCWU1DfSICWFgkmc8eN4Xu/jbkm0ABAdZsVksVNW7AIhw2jmrfxyVdS525lcwsHskj1w1nPiII/dlqqhtICI4qOPekIiIH1O4ETkBr67Yyxsrs7n/siGclh4LeObK/M+CdXy5vdB73WUjkxmYFMk/PjoUbGobPKueIoPtXDUmhd+d07fVINMaBRsRkbbTaimRoyiuquerHYVU17sICbKx92AVT3y6A4C4cCeLbzuTkqoGfv3qavYdrMZpt9IvMZxNOeUc/l/VtWNT+Mvlw8gtq6HBZZASE0JwkM2kdyUi4p+0WkqkjSrrGmlodBMT5gBgd2Elc97bxO7CKvLKa3G5j8z+EcF2iirruOiJr6iub6S2wU1KTAjP3ziGIclRbNhfxktf76GitpF+ieHcdl4/HHYrPbuFdfTbExHpkhRupMuqqXdx6dPLyS6u5nfn9GF4SjR/WrSBgoo67zWDu0eSHB1CbYOL2gYXEwcnMnFQIlfP/4biqnoAzujTjWduGE1sU0AalhLF45NGmvGWREQEhRsJcGuzSnjxy924DYOkyGBG94zhwiFJBAfZ+Ofy3ewpqgLgqc92eh8zMCmCv1w+lJSYUJKiglt93i//eC478iuwWCyMSInGZrV0yPsREZHjU7iRgOV2G9z97+/ZUVDpve9fK/aRFhvKz09PY96yXQD8/PQ01meXUd/opl9iOA9cNtTbhTmayOAgxvSMPaX1i4hI+yjcSMD6eHM+OwoqiXDa+eNFA9hTVM0HG3LIKq7m4cVbARiREsUDPxuKVZ0XEZGAoXAjAWffwSo+3VLAgpVZAEw5oyc3jk8HYMYF/Xnxy93sKKggzGHn1h/1VrAREQkwCjcSUF77bh8Pvr/Zu6dMcJCVWyb08v483GnnjvP7m1WeiIh0AIUbCRjrskq4d9FGAEanRdMnPpwLhiTRLbxtG+WJiEhgULiRgPHNroMA/HhgAv87ZayGm0REuigdMywBY11WCQAT+sYp2IiIdGEKNxIQDMNgXVYpAKPSok2tRUREzKVwIwEhu7iGg1X1OGxWhiTrzDARka5M4UYCwtqmIanByZE47TqUUkSkK9OEYvFrhmGwYtdB3lufA8DotBiTKxIREbMp3IjfMgyDP7+3iX+t2Oe9T/NtRERE4Ub81uOf7vAGmwl9u5ESHcr5gxNNrkpERMymcCN+qa7RxYtf7gbgwcuHcuPpPU2uSEREOgtNKBa/tC6rlJoGF3HhDn6ekWZ2OSIi0oko3IhfWr6jCPBs2GexaMM+ERE5ROFG/NLynYfCjYiIyOEUbsTvlNU08P3+UkDhRkREjqQJxeJX/rRoA+9l5uA2oHdcGD2iQ8wuSUREOhmFG/EbZdUNvP5dlvf2laN7mFiNiIh0Vgo34jc25pQB0CM6hHemTSA+wmlyRSIi0hlpzo34je/3e8LNyLRoBRsRETkqhRvxGxsOlAIwvEeUuYWIiEinpnAjfmN9tqdzMyxF4UZERI5O4Ub8wsHKOg6U1gAwTJ0bERE5BoUb8QsbDni6Nr3jw4gIDjK5GhER6cwUbsQvrNh1ENB8GxEROT6FG+n0KmobeH2lZ3+bi4d1N7kaERHp7BRupFOrqG3gX9/spaK2kT7xYZw/KNHskkREpJPTJn7Sab25Opu7//M9huG5feuP+mC16gRwERE5NnVupNP6aGOeN9gM6xHF5SN13IKIiByfOjfSaW3OLQfgjV+dzum9Y7FY1LUREZHjU+dGOqXiqnpyy2oBGNojUsFGRETaTOFGOqUtTV2bnt1Cta+NiIicEIUb6ZQ253jCzZDkSJMrERERf6NwI51S83ybwd0VbkRE5MQo3EinU9vgYlOO57iFwerciIjICdJqKelUnv9iF3M/3Oq9Pbi7jlsQEZETo86NdBout8E/l+/x3h6QGEFipNPEikRExB+pcyOdxne7D1JQUUdUSBBv/+4MekSHaAm4iIicMIUb6TTezcwB4CfDutMnPtzkakRExF9pWEo6hbpGF4s35gJw2chkk6sRERF/pnAjncKqPSVU1DaSEOFkXHqs2eWIiIgfU7iRTuHb3QcBOLNvnE7+FhGRk6JwI53Cd3s84Sajt7o2IiJychRuxHQ19S4ys0sBOL13N3OLERERv6dwI6Zbl1VCg8sgKTKYtNhQs8sRERE/p6XgYqrPtxXw+ndZAJzeO1b72oiIyElTuBHTrM0qYerLq7y3x/fRkJSIiJw8hRsxzaYDnsMxU2JCuGxkMpeN7GFyRSIiEggUbsQ0u4uqALh4aBJ3XTjQ5GpERCRQaEKxmGZ3oSfc9NZRCyIi4kMKN2KaPU2dm15xYSZXIiIigUThRkxR1+hif0k1AL3jFW5ERMR3FG7EFFkHq3EbEO60Ex/uNLscEREJIAo3Yordhw1JaW8bERHxJYUbMUXzfBsNSYmIiK8p3IgpdhdWAppMLCIivqdwI6bYUaBwIyIip4bCjXS4nNIa7yngo9NizC1GREQCjsKNdLhF6w5gGJDRK5ZUnQIuIiI+1inCzbPPPkt6ejrBwcFkZGSwcuXKo1774osvctZZZxETE0NMTAwTJ0485vXSuRiGwX/W7Afg6jEpJlcjIiKByPRws3DhQmbMmMGcOXNYu3YtI0aM4MILL6SgoKDV65ctW8b111/P559/zooVK0hNTeWCCy7gwIEDHVy5nKg/v7eJAbOWsLuoipAgGxcP6252SSIiEoAshmEYZhaQkZHBaaedxjPPPAOA2+0mNTWV//mf/+Gee+457uNdLhcxMTE888wzTJky5bjXl5eXExUVRVlZGZGRkSddv7RNo8vN4DkfUd/oBuBXZ/Xi3ksGm1yViIj4ixP5/Db1VPD6+nrWrFnDzJkzvfdZrVYmTpzIihUr2vQc1dXVNDQ0EBsb2+rP6+rqqKur894uLy8/uaKlXfYVV1Pf6CY4yMo395xHbJjD7JJERCRAmTosVVRUhMvlIjExscX9iYmJ5OXltek57r77bpKTk5k4cWKrP587dy5RUVHer9TU1JOuW07c9rwKAPonRijYiIjIKWX6nJuT8de//pUFCxawaNEigoODW71m5syZlJWVeb+ys7M7uEoB2JZ/KNyIiIicSqYOS8XFxWGz2cjPz29xf35+PklJScd87D/+8Q/++te/8umnnzJ8+PCjXud0OnE6dTCj2bY3hZuBSQo3IiJyapnauXE4HIwZM4alS5d673O73SxdupTx48cf9XF/+9vfePDBB1myZAljx47tiFLlJG3LU+dGREQ6hqmdG4AZM2Zw0003MXbsWMaNG8cTTzxBVVUVU6dOBWDKlCn06NGDuXPnAvDII48we/ZsXn/9ddLT071zc8LDwwkPDzftfcjR1Ta42HuwGoAB6tyIiMgpZnq4mTRpEoWFhcyePZu8vDxGjhzJkiVLvJOMs7KysFoPNZjmzZtHfX09V199dYvnmTNnDn/+8587snRpo12FlbjcBlEhQSREaIhQREROLdPDDcD06dOZPn16qz9btmxZi9t79+499QWJT63cUwzAgMQILBaLydWIiEig6xThRgJTYUUdmdmlzP1wKwDnDIw3uSIREekKFG7klPjg+1ymv7GW5v2vJw5K5Naz+5hblIiIdAkKN3JKLN6Qi2FATGgQZ/WL569XDcNm1ZCUiIicego34nOGYbB6n2eezfyfjyGjdzeTKxIRka7Er3cols5pf0kN+eV1BNksjEiNNrscERHpYhRuxOeauzZDkqMIDrKZXI2IiHQ1Cjfic6v3lgAwtmeMyZWIiEhXpDk34jP55bW8/l0Wn28tAGBsusKNiIh0PIUb8Zl/fLSNt9bsB8BigTE9Y02uSEREuiKFG/EJwzBYtr0QgAsGJzJxUCLxOmpBRERMoHAjPrElt4LCijpCgmw8fcMonHZNJBYREXNoQrH4xJc7PF2b03vHKtiIiIipFG7EJ75sGpI6u7/OjxIREXMp3MhJq65v9C7/VrgRERGzKdzISduwv4x6l5vESCe948LMLkdERLo4hRs5aev3lwIwKjUGi0WHY4qIiLkUbuSkrc8uA9A5UiIi0iko3MhJy8wuBWBEapS5hYiIiKBwIyepsKKOA6U1WCwwrIfCjYiImE/hRk7K903zbfrGhxMRHGRuMSIiImiHYmmnRpeb+97ZyHd7igHNtxERkc5D4UbaZcXugyxYle29fXrvbiZWIyIicojCjbTL9/s9K6ROS4/hF2f2ZuKgBJMrEhER8VC4kXZZ37RC6oLBSVw0NMncYkRERA6jCcXSLs2dm+EpWiElIiKdi8KNnLCC8lryymuxWmColn+LiEgno3AjJ2x9U9emb0I4YU6NbIqISOeicCMnrHlvm+Ep0abWISIi0hqFGzkhLrfB0i0FAIzQfBsREemEFG7khLzyzV4255YT4bRzoVZJiYhIJ6QJE9ImTy/dwUtf76G8thGAmT8ZREJEsMlViYiIHEnhRo7L5TZ44avdVDQFmzP7xnHdaakmVyUiItI6hRs5ri255VTUNhLhtLNo2gR6dgvFarWYXZaIiEirFG7kuL7dfRCA03rF0jch3ORqREREjk0TiuW4vt3tOfk7o1esyZWIiIgcn8KNHJPbbbBqryfc6ORvERHxBxqWkqN6dcVe3s3MoaymgXCnnSHJkWaXJCIiclwKN9KqXYWVzHlvE27Dc3tC327YbWr0iYhI56dwI616/JPtuA0Y37sb12ekcVbfOLNLEhERaROFGznCltxy3v8+F4BZPx3MYA1HiYiIH9E4gxzhuWW7ALhkeHcFGxER8TsKN9JCdnE1izd4uja/O6ePydWIiIicOA1LCQAVtQ18sb2QdzNzcLkNzuoXx5BknfotIiL+R+GmiyqvbaCkqp6ymgYWrTvAW6v3U1nX6P35b36kro2IiPgnhZsuKKe0houf/IqymoYW9/eOD2NAYgQjUqM5o4827BMREf+kcNMF/W3JVspqGrBbLTjsVs4ZEM+1Y1M5u1+8DsQUERG/p3DTxXy5vZB3MnOwWGDR7yYwLEXzakREJLAo3AQgwzCoqGukrLqBukYX//dtFl/vLKKmwcX+khoArhyVomAjIiIBSeHGR2qbgoPFAhbAYrHgchs0uNw0uNw0ug0Mw8BteA6jdBscum0YuA0Dw/s9TbcP/7nn+toGF1nF1ZTXNOIyDNxug/LaBvLKaml0G1TVNZJXVktVvavVOoNsFi4e2p1ZPx3Usb8gERGRDqJw4yObcsq5at43ZpfRgsNuBQPGpsdw0xnpxIQ66JcQTkyYw+zSREREThmFGx+xWy1EhQRhGAYGgAE2mwW71UqQzYLdZsFq8XxZLDR9T9PtQ99bLfzg9qHrbVbP86TEhBAb5sRmsWCzQqjDTveoYJxBVoLtNpKigkmKCibUof97RUSk69Gnn4+MSI1m/ZwLzC5DRESky9PxCyIiIhJQFG5EREQkoCjciIiISEBRuBEREZGAonAjIiIiAUXhRkRERAKKwo2IiIgEFIUbERERCSgKNyIiIhJQFG5EREQkoCjciIiISEBRuBEREZGAonAjIiIiAUXhRkRERAKK3ewCOpphGACUl5ebXImIiIi0VfPndvPn+LF0uXBTUVEBQGpqqsmViIiIyImqqKggKirqmNdYjLZEoADidrvJyckhIiICi8Xi0+cuLy8nNTWV7OxsIiMjffrcciT9vjuWft8dT7/zjqXfd8c60d+3YRhUVFSQnJyM1XrsWTVdrnNjtVpJSUk5pa8RGRmp/zA6kH7fHUu/746n33nH0u+7Y53I7/t4HZtmmlAsIiIiAUXhRkRERAKKwo0POZ1O5syZg9PpNLuULkG/746l33fH0++8Y+n33bFO5e+7y00oFhERkcCmzo2IiIgEFIUbERERCSgKNyIiIhJQFG5EREQkoCjc+Mizzz5Leno6wcHBZGRksHLlSrNLClhz587ltNNOIyIigoSEBC6//HK2bdtmdlldxl//+lcsFgu333672aUErAMHDvDzn/+cbt26ERISwrBhw1i9erXZZQUkl8vFrFmz6NWrFyEhIfTp04cHH3ywTecXSdt8+eWXXHrppSQnJ2OxWHjnnXda/NwwDGbPnk337t0JCQlh4sSJ7Nix46ReU+HGBxYuXMiMGTOYM2cOa9euZcSIEVx44YUUFBSYXVpA+uKLL5g2bRrffvstn3zyCQ0NDVxwwQVUVVWZXVrAW7VqFc8//zzDhw83u5SAVVJSwoQJEwgKCuLDDz9k8+bNPProo8TExJhdWkB65JFHmDdvHs888wxbtmzhkUce4W9/+xtPP/202aUFjKqqKkaMGMGzzz7b6s//9re/8dRTTzF//ny+++47wsLCuPDCC6mtrW3/ixpy0saNG2dMmzbNe9vlchnJycnG3LlzTayq6ygoKDAA44svvjC7lIBWUVFh9OvXz/jkk0+MH/3oR8Ztt91mdkkB6e677zbOPPNMs8voMi655BLjlltuaXHflVdeaUyePNmkigIbYCxatMh72+12G0lJScbf//53732lpaWG0+k03njjjXa/jjo3J6m+vp41a9YwceJE731Wq5WJEyeyYsUKEyvrOsrKygCIjY01uZLANm3aNC655JIWf9bF99577z3Gjh3LNddcQ0JCAqNGjeLFF180u6yAdcYZZ7B06VK2b98OwPr161m+fDkXX3yxyZV1DXv27CEvL6/F3ytRUVFkZGSc1Gdolzs409eKiopwuVwkJia2uD8xMZGtW7eaVFXX4Xa7uf3225kwYQJDhw41u5yAtWDBAtauXcuqVavMLiXg7d69m3nz5jFjxgz+9Kc/sWrVKn7/+9/jcDi46aabzC4v4Nxzzz2Ul5czcOBAbDYbLpeLhx56iMmTJ5tdWpeQl5cH0OpnaPPP2kPhRvzatGnT2LhxI8uXLze7lICVnZ3NbbfdxieffEJwcLDZ5QQ8t9vN2LFjefjhhwEYNWoUGzduZP78+Qo3p8Cbb77Ja6+9xuuvv86QIUPIzMzk9ttvJzk5Wb9vP6ZhqZMUFxeHzWYjPz+/xf35+fkkJSWZVFXXMH36dN5//30+//xzUlJSzC4nYK1Zs4aCggJGjx6N3W7HbrfzxRdf8NRTT2G323G5XGaXGFC6d+/O4MGDW9w3aNAgsrKyTKoosN11113cc889XHfddQwbNowbb7yRO+64g7lz55pdWpfQ/Dnp689QhZuT5HA4GDNmDEuXLvXe53a7Wbp0KePHjzexssBlGAbTp09n0aJFfPbZZ/Tq1cvskgLaeeedx4YNG8jMzPR+jR07lsmTJ5OZmYnNZjO7xIAyYcKEI7Y22L59Oz179jSposBWXV2N1dryo9Bms+F2u02qqGvp1asXSUlJLT5Dy8vL+e67707qM1TDUj4wY8YMbrrpJsaOHcu4ceN44oknqKqqYurUqWaXFpCmTZvG66+/zrvvvktERIR3XDYqKoqQkBCTqws8ERERR8xnCgsLo1u3bprndArccccdnHHGGTz88MNce+21rFy5khdeeIEXXnjB7NIC0qWXXspDDz1EWloaQ4YMYd26dTz22GPccsstZpcWMCorK9m5c6f39p49e8jMzCQ2Npa0tDRuv/12/vKXv9CvXz969erFrFmzSE5O5vLLL2//i57Eii45zNNPP22kpaUZDofDGDdunPHtt9+aXVLAAlr9evnll80urcvQUvBT67///a8xdOhQw+l0GgMHDjReeOEFs0sKWOXl5cZtt91mpKWlGcHBwUbv3r2Ne++916irqzO7tIDx+eeft/p39k033WQYhmc5+KxZs4zExETD6XQa5513nrFt27aTek2LYWgbRhEREQkcmnMjIiIiAUXhRkRERAKKwo2IiIgEFIUbERERCSgKNyIiIhJQFG5EREQkoCjciIiISEBRuBEREZGAonAjIqZZtmwZFouF0tJSU15/6dKlDBo0qE2Hfy5ZsoSRI0fqzCERP6BwIyId4pxzzuH2229vcd8ZZ5xBbm4uUVFRptT0xz/+kfvuu69Nh39edNFFBAUF8dprr3VAZSJyMhRuRMQ0DoeDpKQkLBZLh7/28uXL2bVrF1dddVWbH3PzzTfz1FNPncKqRMQXFG5E5JS7+eab+eKLL3jyySexWCxYLBb27t17xLDUK6+8QnR0NO+//z4DBgwgNDSUq6++murqav71r3+Rnp5OTEwMv//971sMJdXV1XHnnXfSo0cPwsLCyMjIYNmyZcesacGCBZx//vkEBwd771u/fj3nnnsuERERREZGMmbMGFavXu39+aWXXsrq1avZtWuXT38/IuJbdrMLEJHA9+STT7J9+3aGDh3KAw88AEB8fDx79+494trq6mqeeuopFixYQEVFBVdeeSVXXHEF0dHRLF68mN27d3PVVVcxYcIEJk2aBMD06dPZvHkzCxYsIDk5mUWLFnHRRRexYcMG+vXr12pNX331FTfccEOL+yZPnsyoUaOYN28eNpuNzMxMgoKCvD9PS0sjMTGRr776ij59+vjotyMivqZwIyKnXFRUFA6Hg9DQUJKSko55bUNDA/PmzfOGh6uvvppXX32V/Px8wsPDGTx4MOeeey6ff/45kyZNIisri5dffpmsrCySk5MBuPPOO1myZAkvv/wyDz/8cKuvs2/fPu/1zbKysrjrrrsYOHAgQKvBKDk5mX379p3w70BEOo7CjYh0KqGhoS26IomJiaSnpxMeHt7ivoKCAgA2bNiAy+Wif//+LZ6nrq6Obt26HfV1ampqWgxJAcyYMYNf/vKXvPrqq0ycOJFrrrnmiA5NSEgI1dXV7X5/InLqKdyISKdy+DAQgMViafW+5iXZlZWV2Gw21qxZc8Sqp8MD0Q/FxcVRUlLS4r4///nP3HDDDXzwwQd8+OGHzJkzhwULFnDFFVd4rykuLiY+Pr5d701EOobCjYh0CIfD0ab9ZE7UqFGjcLlcFBQUcNZZZ53Q4zZv3nzE/f3796d///7ccccdXH/99bz88svecFNbW8uuXbsYNWqUz+oXEd/TaikR6RDp6el899137N27l6KiIp9thte/f38mT57MlClTePvtt9mzZw8rV65k7ty5fPDBB0d93IUXXsjy5cu9t2tqapg+fTrLli1j3759fP3116xatYpBgwZ5r/n2229xOp2MHz/eJ7WLyKmhcCMiHeLOO+/EZrMxePBg4uPjycrK8tlzv/zyy0yZMoU//OEPDBgwgMsvv5xVq1aRlpZ21MdMnjyZTZs2sW3bNgBsNhsHDx5kypQp9O/fn2uvvZaLL76Y+++/3/uYN954g8mTJxMaGuqz2kXE9yyGYRhmFyEiYoa77rqL8vJynn/++eNeW1RUxIABA1i9ejW9evXqgOpEpL3UuRGRLuvee++lZ8+ebRoi27t3L88995yCjYgfUOdGREREAoo6NyIiIhJQFG5EREQkoCjciIiISEBRuBEREZGAonAjIiIiAUXhRkRERAKKwo2IiIgEFIUbERERCSgKNyIiIhJQ/j+HkBCqqLjswgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAYAAABB4NqyAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAVB5JREFUeJzt3Xt8k/Xd//FXkjbpudBzCy2UM8hRQETAoSCU+VNRnCdU2PT2VkGHeLiH83zq1M0pm8PpFNwU2dzGUJwwZRREATmIgGjlUCintrTQpsc0Ta7fH6WBSgstpE2bvJ+PRx4jV65c+SRm5M33aDIMw0BEREQkgJh9XYCIiIhIa1MAEhERkYCjACQiIiIBRwFIREREAo4CkIiIiAQcBSAREREJOApAIiIiEnCCfF1AW+R2uzl06BCRkZGYTCZflyMiIiJNYBgGpaWlpKSkYDafvo1HAagBhw4dIjU11ddliIiIyFnYv38/nTt3Pu05CkANiIyMBGo/wKioKB9XIyIiIk1ht9tJTU31/I6fjgJQA+q6vaKiohSARERE2pmmDF/RIGgREREJOApAIiIiEnAUgERERCTgKACJiIhIwFEAEhERkYCjACQiIiIBRwFIREREAo4CkIiIiAQcBSAREREJOApAIiIiEnAUgERERCTgKACJiIhIwFEAEhEREa8wDIOj5dXUuNy+LuWMtBu8iIiIeMWjS7bzzrpcTCbonRjJhH6JzLi0B7Ygi69LO4UCkIiIiHjFuj1HATAM+C6vlO/ySnEZBg9O7OPjyk6lLjARERHxiuIKJwB//tkF/PLHfQH402c5HC6p9GVZDVIAEhERkXNmGAYlldUA9EiI4PYx6VzQNQZHjZuJv13N+U9/wsrvCnxc5QkKQCIiInLOKp0unC4DgA5hwZhMJh6+vC9mE9irajhaXs3MhZv5Ls/u40praQyQiIiInLO67q9gi4nQ4NpBz4NTO/DhPaOxV9bwu//u5IvdRVw+dw0DOkVz/fBUbrwgzWf1qgVIREREzlldAIoOtWIymTzHz0uJZmT3WP4w9XyGpHXA5TbYsr+YvJIqX5UKqAVIREREvKCksjYAdQgLbvDxDmFWFt89igPHKli/5ygDO0e3ZnmnUAASERGRc1Y3ADo6tOEAVKdzxzA6Dw1rjZJOS11gIiIics7qusA6nCEAtRUKQCIiInLOio93gUU30gXW1igAiYiIyDmrGwN0pi6wtkIBSERERM7ZiS4wq48raRoFIBERETlndYOgG5sF1tYoAImIiMg587QAKQCJiIhIoKgbAxSlMUAiIiISKDQNXkRERALOiZWgNQhaREREAoDT5abMUQNoGryIiIgECPvx1h+AqJD2scuWApCIiIick7pVoCNDggiytI9o0T6qFBERkTarbgB0e+n+AgUgEREROUe7j5QB7WcNIID20VEnIiIibUphmYOnl+6gqKyaL3YXAjCgUwffFtUMPm0BmjdvHgMHDiQqKoqoqChGjhzJxx9/DMDRo0e555576N27N6GhoaSlpXHvvfdSUlJy2mtOnz4dk8lU75aRkdEab0dERCRgfLztMEu2HGLNrkLcBlwzpBOPXN7X12U1mU9bgDp37syvfvUrevbsiWEYvP3221x11VV89dVXGIbBoUOH+PWvf02/fv3Yt28fd955J4cOHeLvf//7aa+bkZHB/PnzPfdtNltLvxUREZGAUlhWu/fXqB6x3HNpT0akx2AymXxcVdP5NABdccUV9e4/++yzzJs3j3Xr1nHbbbfxj3/8w/NY9+7defbZZ7n55pupqakhKKjx0m02G0lJSU2uw+Fw4HA4PPftdnsz3oWIiEjgKa6oDUCDUztwYbdYH1fTfG1mELTL5WLRokWUl5czcuTIBs8pKSkhKirqtOEHICsri4SEBHr37s1dd91FUVHRac/PzMwkOjrac0tNTT3r9yEiIhIIjh2f+dWxnaz8/EM+D0Dbtm0jIiICm83GnXfeyeLFi+nXr98p5xUWFvL0009zxx13nPZ6GRkZ/PnPf2bFihU8//zzrFq1ikmTJuFyuRp9zpw5cygpKfHc9u/ff87vS0RExJ8Vt7OtL37I57PAevfuzZYtWygpKeHvf/8706ZNY9WqVfVCkN1u5/LLL6dfv3488cQTp73eDTfc4PnzgAEDGDhwIN27dycrK4tx48Y1+BybzaZxQiIiIs1Q1wXWXjY//SGftwBZrVZ69OjB0KFDyczMZNCgQbzyyiuex0tLS8nIyCAyMpLFixcTHNy8D7pbt27ExcWxa9cub5cuIiISsI4dD0AdwxWAvMLtdnsGJNvtdiZMmIDVauWDDz4gJCSk2dc7cOAARUVFJCcne7tUERGRgFVc3r67wHwagObMmcPq1avZu3cv27ZtY86cOWRlZTF16lRP+CkvL+fNN9/EbreTl5dHXl5evfE8ffr0YfHixQCUlZXx4IMPsm7dOvbu3cuKFSu46qqr6NGjBxMnTvTV2xQREfErTpeb0uO7v7fXLjCfjgEqKCjg1ltv5fDhw0RHRzNw4ECWL1/OZZddRlZWFuvXrwegR48e9Z6Xk5ND165dAcjOzvYsjmixWNi6dStvv/02xcXFpKSkMGHCBJ5++mmN8REREfGSkpN2f29P+3+dzKcB6M0332z0sbFjx2IYxhmvcfI5oaGhLF++3Cu1iYiISMPqBkBHtaPd33+ofVYtIiIiPlO3+3vH8PY5/gcUgERERKSZ6hZBbK/jf0ABSERERJqpbgp8e50BBgpAIiIi0kx1Y4A6hqkFSERERAJE3RggtQCJiIhIwPCMAVILkIiIiASKE11gagESERGRAFGsFiAREREJNP4wC8ynK0GLiIhI+5FbVMEzH+1gb1E50L5ngSkAiYiISJMs2pDLf3bkA2Axm0iKDvFxRWdPAUhERESaZP+xSgAm9U/i5gu7kBDZfgOQxgCJiIhIkxw4VgHAlYNSGNUjzsfVnBsFIBEREWmSA8dbgDp3DPNxJedOAUhERETOqMrp4kipA4DOHUN9XM25UwASERGRMzpYXNv6E261tOv1f+ooAImIiMgZ7T9aO/6nc8cwTCaTj6s5dwpAIiIickYnxv+0/+4vUAASERGRJqgLQKkx7X8ANCgAiYiISBPUTYFXC5CIiIgEDHWBiYiISMDxpzWAQAFIREREzsBR46KwrHYNoJQOagESERGRAGCvrPH8OTq0/a8BBApAIiIicgalVU4AIm1BWMztfw0gUAASERGRMyitqm0BigwJ8nEl3qMAJCIiIqdlr2sBCvGP7i9QABIREZEzUAuQiIiIBBzPGCAFIBEREQkUdS1AUX4yAwwUgEREROQM7JVqARIREZEAY/eMAVILkFfMmzePgQMHEhUVRVRUFCNHjuTjjz/2PF5VVcWMGTOIjY0lIiKCKVOmkJ+ff9prGobBY489RnJyMqGhoYwfP56dO3e29FsRERHxWxoE7WWdO3fmV7/6FZs2bWLjxo1ceumlXHXVVXzzzTcA3HfffXz44Ye8//77rFq1ikOHDnHNNdec9povvPACc+fO5bXXXmP9+vWEh4czceJEqqqqWuMtiYiI+J26QdBRftQCZDIMw/B1ESeLiYnhxRdf5NprryU+Pp6FCxdy7bXXAvDdd9/Rt29f1q5dy4UXXnjKcw3DICUlhfvvv58HHngAgJKSEhITE1mwYAE33HBDk2qw2+1ER0dTUlJCVFSU996ciIhIO3Tj6+tYu6eIV24YzFWDO/m6nEY15/e7zYwBcrlcLFq0iPLyckaOHMmmTZtwOp2MHz/ec06fPn1IS0tj7dq1DV4jJyeHvLy8es+Jjo5mxIgRjT4HwOFwYLfb691ERESklt0PW4B8HoC2bdtGREQENpuNO++8k8WLF9OvXz/y8vKwWq106NCh3vmJiYnk5eU1eK2644mJiU1+DkBmZibR0dGeW2pq6rm9KRERET+iMUAtoHfv3mzZsoX169dz1113MW3aNHbs2NGqNcyZM4eSkhLPbf/+/a36+iIiIm1ZqR9uheHzKGe1WunRowcAQ4cOZcOGDbzyyitcf/31VFdXU1xcXK8VKD8/n6SkpAavVXc8Pz+f5OTkes8ZPHhwozXYbDZsNtu5vxkRERE/YxjGSQsh+jw2eI3PW4B+yO1243A4GDp0KMHBwaxYscLzWHZ2Nrm5uYwcObLB56anp5OUlFTvOXa7nfXr1zf6HBEREWlcldNNjbt2vpRagLxkzpw5TJo0ibS0NEpLS1m4cCFZWVksX76c6OhobrvtNmbPnk1MTAxRUVHcc889jBw5st4MsD59+pCZmcnVV1+NyWRi1qxZPPPMM/Ts2ZP09HQeffRRUlJSmDx5su/eqIiISDtVNwDabIJwq8XH1XiPTwNQQUEBt956K4cPHyY6OpqBAweyfPlyLrvsMgB++9vfYjabmTJlCg6Hg4kTJ/KHP/yh3jWys7MpKSnx3H/ooYcoLy/njjvuoLi4mNGjR7Ns2TJCQkJa9b2JiIj4g7rxPxG2IEwmk4+r8Z42tw5QW6B1gERERGptzj3GNX/4gs4dQ1nzf5f6upzTapfrAImIiEjbU+qH+4CBApCIiIicxokp8P4zAwwUgEREROQ07JXHp8ArAImIiEig8MdFEEEBSERERBrxp8/28M/NBwH/awHyr3cjIiIiXlFQWsUzH33ruZ8Y7V/LySgAiYiIyClKKmq7vkKDLTyU0ZspQzv7uCLvUgASERGRU5Q5agc/x4Rb+emodB9X430aAyQiIiKnqKh2AbUrQPsjBSARERE5RV0LULjNf/b/OpkCkIiIiJyi3BOA1AIkIiIiAcITgKwKQCIiIhIgyo+PAVILkIiIiASMuhagCI0BEhERkUBRpjFAIiIiEmg0CFpEREQCTrnj+Bggq7rAREREJECUV6sFSERERALMiUHQCkAiIiISIMqOd4GFKQCJiIhIoNA0eBEREQk4FRoDJCIiIoGmTFthiIiISCCpcbmpcroBtQCJiIhIgKjbBwwgXGOAREREJBDUjf8JtpiwBSkAiYiISADw920wQAFIREREfqDMsw2GApCIiIgEiBMtQP7Z/QUKQCIiIvID6gITERGRgFO3Eaq/7gMGCkAiIiLyA559wKzqAhMREZEAoS6wFpaZmcnw4cOJjIwkISGByZMnk52d7Xl87969mEymBm/vv/9+o9edPn36KednZGS0xlsSERFp9yoc6gJrUatWrWLGjBmsW7eOTz75BKfTyYQJEygvLwcgNTWVw4cP17s9+eSTREREMGnSpNNeOyMjo97z3nvvvdZ4SyIiIu2eZxq8Hwcgn76zZcuW1bu/YMECEhIS2LRpExdffDEWi4WkpKR65yxevJjrrruOiIiI017bZrOd8tzGOBwOHA6H577dbm/iOxAREfE/ni4wjQFqHSUlJQDExMQ0+PimTZvYsmULt9122xmvlZWVRUJCAr179+auu+6iqKio0XMzMzOJjo723FJTU8/uDYiIiPgBe5UTgMiQYB9X0nJMhmEYvi4CwO12c+WVV1JcXMyaNWsaPOfuu+8mKyuLHTt2nPZaixYtIiwsjPT0dHbv3s3DDz9MREQEa9euxWI5Nc021AKUmppKSUkJUVFR5/bGRERE2pmfvPYFG/Ye49Wbzufygcm+LqfJ7HY70dHRTfr9bjOdezNmzGD79u2Nhp/KykoWLlzIo48+esZr3XDDDZ4/DxgwgIEDB9K9e3eysrIYN27cKefbbDZsNtvZFy8iIuJHjpTWNgrER/rvb2Ob6AKbOXMmS5cuZeXKlXTu3LnBc/7+979TUVHBrbfe2uzrd+vWjbi4OHbt2nWupYqIiPi9wrJqAOIirD6upOX4tAXIMAzuueceFi9eTFZWFunp6Y2e++abb3LllVcSHx/f7Nc5cOAARUVFJCe3n2Y8ERERX6isdlF2fBC0WoBayIwZM3jnnXdYuHAhkZGR5OXlkZeXR2VlZb3zdu3axerVq7n99tsbvE6fPn1YvHgxAGVlZTz44IOsW7eOvXv3smLFCq666ip69OjBxIkTW/w9iYiItGeFZbXdX7Ygs9YBainz5s2jpKSEsWPHkpyc7Ln99a9/rXfeW2+9RefOnZkwYUKD18nOzvbMILNYLGzdupUrr7ySXr16cdtttzF06FA+++wzjfMRERE5g4Lj43/iImyYTCYfV9NyfN4F1hTPPfcczz33XJOuExoayvLly8+5NhERkUBU1wLkz91f0EYGQYuIiEjbUBeA4iIUgERERCRABMIUeFAAEhERkZN4usD8eAo8KACJiIjISepagOLUAiQiIiLt2YFjFewqKG3SuXWLIMZrDJCIiIi0Vy63wfV/XMcVv/ucY+XVZzxfLUAiIiLS7u04ZOdgcSWVThe5RyvOeP6JMUAKQCIiItJOfbG70PPnutadxpQ7aqiodgH+3wLkv2tci4iICGv3FHn+XNe605APvz7Ei8uzAQgJNhNutbR4bb6kFiARERE/5XS5+TLnqOf+6VqAFnyx19NFdmG3WL/eBgPUAiQiIuKXKqprWLY9z9OlBadvAaoLP2/cOoxL+yS0eH2+pgAkIiLiZ1xug0mvfMa+otpQYzKBYcCRRgJQZbXL0zp0QdcYLGb/bv0BdYGJiIj4naPl1Z7wM6pHLP97cXcACksbngZ/4FjtuZEhQUSHBbdOkT6mACQiIuJnjh5f7ycm3Mq7t1/IJb3jgcZbgPYfD0BpMWGtU2AboAAkIiLiZ4qOB53Y8Nr9vOqmtBc2Mgg693hrUWpHBSARERFppwpPagGCEzu7lzpqqHK6MAyD1d8f8QSl3KOVAKTFKgCJiIhIO3X0eLCJO76ac6QtCGtQ7U/+kVIHn35bwK1vfcm9i74CTnSBpXYM9UG1vqEAJCIi4meKftACZDKZPFtbHClzsOLbfAA+31XEgWMV7D8+BT41gMYAndU0+NzcXPbt20dFRQXx8fGcd9552Gz+vWS2iIhIe1EXgGIjrJ5jcZE2DhZXUljq4POTtsdYsuWQAtDp7N27l3nz5rFo0SIOHDiAYRiex6xWK2PGjOGOO+5gypQpmM1qWBIREfGVHw6ChhObm27KPcb+42N+AN7+Yi/l1S5MJujUQV1g9dx7770MGjSInJwcnnnmGXbs2EFJSQnV1dXk5eXx73//m9GjR/PYY48xcOBANmzY0NJ1i4iISCOOelqATvTOxEfWhqElXx0CoE9SJFaLmYLjM8MSI0MICfbv/b9O1qQWoPDwcPbs2UNsbOwpjyUkJHDppZdy6aWX8vjjj7Ns2TL279/P8OHDvV6siIiInFlRWf0xQHCiBSjPXgVARv8kyh01vPFZDiYTXDU4pfUL9aEmBaDMzMwmXzAjI+OsixEREZFzV+iZBXYiAP1wfM+YnvEM7dKROZP6YjLh95uf/tBZ7wVWUFBAdnY2AL179yYhwf83ThMREWnrqmvc2KtqAIgJP9EFdsWgFCqqXRSVV5PaMZShXToCYA6Afb8a0uwAVFpayt13382iRYtwuWp3mLVYLFx//fW8+uqrREdHe71IERERaZpjFbXdX2YTdAg9sa9XSLCFaRd19VFVbU+zp2vdfvvtrF+/nqVLl1JcXExxcTFLly5l48aN/O///m9L1CgiIiJNdGL8jy1gW3eaotktQEuXLmX58uWMHj3ac2zixIm88cYbGv8jIiLiY0Xlp06Bl1M1uwUoNja2wW6u6OhoOnbs6JWiRERE5OwcbWARRDlVswPQI488wuzZs8nLy/Mcy8vL48EHH+TRRx/1anEiIiLSPIUNTIGXUzWpC2zIkCH1psft3LmTtLQ00tLSgNqtMWw2G0eOHNE4IBERkRZQXFFNdY2bmHArQZbG2y8Kf7ARqjSsSQFo8uTJLVyGiIiINGb7wRKu/P0a3AZ0DAtm6b1jGty2wlHj4sOva1d67hYf3tpltitNCkCPP/54S9chIiIijcjKLsB9fAvOYxVO/vttPreM7FrvnCqni3fW7ePAsUoSIm1cO7Rz6xfajvh019LMzEyGDx9OZGQkCQkJTJ482bO4Yp2xY8diMpnq3e68887TXtcwDB577DGSk5MJDQ1l/Pjx7Ny5syXfioiISIv59nApACHBtT/bX+UW13v8T5/toc+jy3jmo28BmDW+F2HWs17rOCA0KQB17NiRmJiYJt2aY9WqVcyYMYN169bxySef4HQ6mTBhAuXl5fXO+5//+R8OHz7sub3wwgunve4LL7zA3Llzee2111i/fj3h4eFMnDiRqqqqZtUnIiLSFnx72A7A1BFdANice8zzmNtt8OaaHM/9QakduG6YWn/OpEnx8OWXX/b8uaioiGeeeYaJEycycuRIANauXcvy5cubPQts2bJl9e4vWLCAhIQENm3axMUXX+w5HhYWRlJSUpOuaRgGL7/8Mo888ghXXXUVAH/+859JTEzkX//6FzfccEOzahQREfGliuoacopqGwZuvCCNN9fksLeogqIyB7ERNjbsPcrhkioiQ4JY9eAldAwLDrh9vc5Gk1qApk2b5rl9/vnnPPXUU7z33nvce++93Hvvvbz33ns89dRTrFq16pyKKSkpATilJendd98lLi6O/v37M2fOHCoqKhq9Rk5ODnl5eYwfP95zLDo6mhEjRrB27doGn+NwOLDb7fVuIiIibUF2XimGUbuxaY+ECHokRAAnusE+OD7oOeO8JGLCrQo/TdTsMUDLly9vcMXnjIwMPv3007MuxO12M2vWLEaNGkX//v09x2+66SbeeecdVq5cyZw5c/jLX/7CzTff3Oh16tYnSkxMrHc8MTGx3tpFJ8vMzCQ6OtpzS01NPev3ISIi4k1143/6JkcBMCS1AwCvrdrN//x5I0u21AagKwen+KS+9qrZI6RiY2NZsmQJ999/f73jS5YsITY29qwLmTFjBtu3b2fNmjX1jt9xxx2ePw8YMIDk5GTGjRvH7t276d69+1m/3snmzJnD7NmzPfftdrtCkIiItAl143/qAtD5XTry/qYDbNx3YhxQYpSNkd3O/jc4EDU7AD355JPcfvvtZGVlMWLECADWr1/PsmXLeOONN86qiJkzZ7J06VJWr15N586nH7hV95q7du1qMADVjRXKz88nOTnZczw/P5/Bgwc3eE2bzYbNpgWjRESkbXG5Db7aXxt0+iZHAnBZv0Te/mIv4bYgJvVPIthiZmT32NMujiinanYAmj59On379mXu3Ln885//BKBv376sWbPGE06ayjAM7rnnHhYvXkxWVhbp6elnfM6WLVsA6oWbk6Wnp5OUlMSKFSs8gcdut7N+/XruuuuuZtUnIiLSGpwuN8UVTqJCg7AFWdhXVM76PUf5ePthth+0YzGbGJpWOz42LsLGslkXn+GKciZntUjAiBEjePfdd8/5xWfMmMHChQtZsmQJkZGRnjE60dHRhIaGsnv3bhYuXMiPf/xjYmNj2bp1K/fddx8XX3wxAwcO9FynT58+ZGZmcvXVV2MymZg1axbPPPMMPXv2JD09nUcffZSUlBStaC0iIm3O0q2HmP3Xr6l2uencMZQlM0ZxzR++oOj4pqbWIDO/+ckg0mLDfFypf2lSACovLyc8vOlLajf1/Hnz5gG1ix2ebP78+UyfPh2r1cqnn37Kyy+/THl5OampqUyZMoVHHnmk3vnZ2dmeGWQADz30EOXl5dxxxx0UFxczevRoli1bRkhISJPfg4iISGtYvPkg1S43AAeOVfLTBRsoKq8mLsLKpX0SuOXCrgzoHO3jKv2PyTAM40wnJScn8/Of/5xp06Y12vVkGAaffvopL730EhdffDFz5szxerGtxW63Ex0dTUlJCVFRUb4uR0RE/JRhGAx/dgWFZQ6uOb8T/9x80PPYE1f0Y/qoMw8NkROa8/vdpBagrKwsHn74YZ544gkGDRrEsGHDSElJISQkhGPHjrFjxw7Wrl1LUFAQc+bM0Y7wIiIijdhXVE5USDAdw60cLqmisMxBkNnE4//vPFZ/X0hhmYPIkCB+MkyzkVtSkwJQ7969+cc//kFubi7vv/8+n332GV988QWVlZXExcUxZMgQ3njjDSZNmoTFYmnpmkVERNqlLfuLuXbeF3SPj+Djn4/h6/3FAPROiiQ6LJgZl3TnyQ938LNR6YTbtJdXS2rWp5uWlsb9999/yhpAIiIicnqGYfDkh99Q4zbIzi9l1c4jfH2gdvzqwM4dAJh+UVfG9k6gS4wGPLc0xUsREZGztGZnIXe/u4nYCBtdY8OICbdx+5h0z6KFAFVOF6+s2MmeI2X1dnH/8xd7qXLWDn4enFo7yNlkMpEe1/RJR3L2FIBERETO0t837cdeVYO9qoacwtoNS3cdKWPJjFGec57797f8ee0+z/2fDO3M+5sOkPX9EazHFy+sawGS1qMAJCIicpY2H2/R+b+MPkSHBvPEB9/w9f5ivso9xpC0jnyyI98Tfv734m6kdAjlphFpHClzkJV9BEeNm6iQIHoe3+BUWo8CkIiIyFkoKK0i92gFJhPcfGEakSHBbNp3jH9sPsD8z2u7t+597ysAbh+dzpwf9/U89/kpA/nw60PUuA1GdtM2Fr7QrE+8pqaGp556igMHDrRUPSIiIu3C5n3FAPROjCQyJBiAn47qCsAHXx/ixjfWUel0MaZnHA9m9K733MSoEG4f0407f9SdQcd3d5fW1awAFBQUxIsvvkhNTU1L1SMiItIubM6t3aT0/C4dPcf6d4rm8gEnFgy+rF8ib9w6DFuQlohpa5rdBXbppZeyatUqunbt2gLliIiInJ2DxZX897sCokODGds7nqjjrTIt4VBxJev2FAEwNK1jvcdenXo+zztqcDhdxEbYWqwGOTfNDkCTJk3iF7/4Bdu2bWPo0KGn7Pl15ZVXeq04ERGRprr/b1tYt+coAF1iw/j452MIs3p3qOuBYxX8cvF2Vn1/xHNsaJeOp5wXYQsiQgsZtmlN2gvsZGZz471mJpMJl8t1zkX5mvYCExFpX0oqnAx5+j+4DegYFsyxCid3/qg7v5jUp9nXKq1yUlLpJCkqpN7g5P1HK/h/v1tDSaUTi9lEx7BgRveI47fXD8ZkMnnz7chZ8vpeYCdzu91nXZiIiEhLWL3zCG4DeiZE8H8Zfbj9zxv502d7uOb8TvRKjGzwOVVOF/M/38uVg1Po1CEUgMeWbPdMWx/bO57504djMplwuw0e/PvXlFQ66Zccxe9uGkL3eE1db880705ERNq9rOzaLqlL+iQwvl8il/VLpMZt8MqKnY0+5801OTy/7DseWbwNgIrqGhZt2F/vmn/I2s3kVz9n2LOfsm7PUUKDLcy7+XyFHz9wVh2UK1asYMWKFRQUFJzSIvTWW295pTAREZHTcbrcFJY5qKh2ecbkjO0VD8B943vxyY58Pt52mJzC8ga3l6h7zppdhZRUOFmfU0R1jZu0mDCuHtKJV1bs5MXl2fWe88vL+9IlVltV+INmB6Ann3ySp556imHDhpGcnKx+TxERaXU7Dtm5/o9rKXWcWJYl3GphWNcYAPqlRHFJ73hWZh/hob9/zcjucYzvm8CATtGYTCbKHTV8dXwau9Nl8Om3+XyZUzuA+tI+Cdw1tjv/2HyAA8cqGZzagSevPI8wq4WejXSnSfvT7AD02muvsWDBAm655ZaWqEdEROSM/rZxP6WOGixmE2FWC1VOFzeP7II16MTIjrsv6cHK7CNs2HuMDXuPMXfFTsb0jONP04axPqcIp+vEHKB/bzvs2Zl9fN9EQoItLLz9QlbtPMKU8zt5fTaZ+F6z/4tWV1dz0UUXtUQtIiIiZ+R2GyzbngfA67cMZVzfRAzDOKVHYnjXGF65YTDfHi5l/7EKPtmRz2c7C3n+42wMasPPsC4d2bjvGCu+KwBqp69fkF7bipQWG8YtsV1a8Z1Ja2r2IOjbb7+dhQsXtkQtIiIiZ/T1gWLy7FWEWy2M6hEH0OhwjKsGd+IXk/rw6k3n89rN5wPw1uc5LPqydrDzz0anM6BTtOf8Hw9IqteKJP6r2S1AVVVVvP7663z66acMHDiQ4OD6K22+9NJLXitORJrnUHEl2fmlhAZbGN41BotZY/TEvxwqruRvG2vDy6XHu6qa6tI+ifxsVDpvfZ5DpdNFVEgQo3rE8aNe8ewsKCPYYqK3xvgEjGYHoK1btzJ48GAAtm/fXu8xDYgW8Z1yRw2XvbSK8uraxUifmdyfmy9U8734j6zsAqbP3+C5n3FeUrOv8cjlfZk0IAmzCdLjIogOrf1H/GBtSBpwmh2AVq5c2RJ1iMg5Kiqr9oQfgN1HynxYjYj3fXu4FKid7TWiWyzj+iY0+xpms4nhx2eKSWA7545Ou93Ov/71L7777jtv1CMiZ8n5gzW57JU1jZwp0j5V19R+x68a0om3pg9vVveXyA81OwBdd911/P73vwegsrKSYcOGcd111zFgwAD+8Y9/eL1AEWmaGlf9bf1KKp0+qkSkZVQf32vSatEgZTl3zf4WrV69mjFjxgCwePFiDMOguLiYuXPn8swzz3i9QBFpGqfrBy1AVQpA4l8cztrvuE2ztMQLmv0tKikpISamtv902bJlTJkyhbCwMC6//HJ27mx8zxURaVk17votQHa1AImfqXYpAIn3NPtblJqaytq1aykvL2fZsmVMmDABgGPHjhESEuL1AkWkaWp+2AKkACR+pm4MkNbpEW9o9iywWbNmMXXqVCIiIujSpQtjx44FarvGBgwY4O36RKSJ6pb1D7aYcLoMjQESv6MAJN7U7AB09913c8EFF7B//34uu+wyzObaL2K3bt00BkjEh2qOzwKLCbeSb3dQXu2ixuUmSANGxU84auq6wDT7S87dWe3uNmzYMIYNG4ZhGJ79Vy6//HJv1yYizVA3Cywm3Ea+3QGAvaqGmHCrL8sS8RqHWoDEi87qW/Tmm2/Sv39/QkJCCAkJoX///vzpT3/ydm0i0gx1s8BCgs2EW2v/haxxQOJP6gZBaxq8eEOzW4Aee+wxXnrpJe655x5GjhwJwNq1a7nvvvvIzc3lqaee8nqRInJmdbPAgs1mokODKa92aRyQ+BWH8/g6QGoBEi9odgCaN28eb7zxBjfeeKPn2JVXXsnAgQO55557FIBEfKSuBSjIYiIqNJhDJVVaC0j8iqbBizc1+1vkdDoZNmzYKceHDh1KTU3zlt7PzMxk+PDhREZGkpCQwOTJk8nOzvY8fvToUe655x569+5NaGgoaWlp3HvvvZSUlJz2utOnT8dkMtW7ZWRkNKs2kfambgxQkMVM1PENHtUCJP5Es8DEm5r9LbrllluYN2/eKcdff/11pk6d2qxrrVq1ihkzZrBu3To++eQTnE4nEyZMoLy8HIBDhw5x6NAhfv3rX7N9+3YWLFjAsmXLuO2228547YyMDA4fPuy5vffee82qTaS9qZsFFmw2ERVSG4C0H5j4Ew2CFm86q1lgb775Jv/5z3+48MILAVi/fj25ubnceuutzJ4923PeSy+9dNrrLFu2rN79BQsWkJCQwKZNm7j44ovp379/vf3FunfvzrPPPsvNN99MTU0NQUGNl2+z2UhKSmrS+3E4HDgcDs99u93epOeJtCVOTwuQiQibWoDE/1RrGrx4UbMD0Pbt2zn//PMB2L17NwBxcXHExcWxfft2z3kmk6nZxdR1bdVttdHYOVFRUacNPwBZWVkkJCTQsWNHLr30Up555hliY2MbPDczM5Mnn3yy2fWKtCU1njFAZqJCa///oTFA4k9OBCC1AMm5a3YAWrlyZUvUgdvtZtasWYwaNYr+/fs3eE5hYSFPP/00d9xxx2mvlZGRwTXXXEN6ejq7d+/m4YcfZtKkSaxduxaL5dR/OcyZM6dey5Xdbic1NfXc3pBIKzsxC8xEtMYAiR/yTINXABIvOKsusJYwY8YMtm/fzpo1axp83G63c/nll9OvXz+eeOKJ017rhhtu8Px5wIABDBw4kO7du5OVlcW4ceNOOd9ms2Gz2c6pfhFfc548CNozBkgBSPxH3TR4tQCJN7SJb9HMmTNZunQpK1eupHPnzqc8XlpaSkZGBpGRkSxevJjg4OBmXb9bt27ExcWxa9cub5Us0ubUTYMPtqgFSPyTWoDEm3z6LTIMg5kzZ7J48WL++9//kp6efso5drudCRMmYLVa+eCDD85qx/kDBw5QVFREcnKyN8oWaZM8Y4DMJ6bB26s0C0z8g9tteFo5tRK0eINPv0UzZszgnXfeYeHChURGRpKXl0deXh6VlZXAifBTXl7Om2++id1u95zjcrk81+nTpw+LFy8GoKysjAcffJB169axd+9eVqxYwVVXXUWPHj2YOHGiT96nSGtwuk/MAqtrAVIXmPiLutYfUAuQeIdPxwDVrSc0duzYesfnz5/P9OnT2bx5M+vXrwegR48e9c7Jycmha9euAGRnZ3tmkFksFrZu3crbb79NcXExKSkpTJgwgaefflrjfMSv1Xi6wE6aBaYAJH6ibg0g0DR48Q6fBiDDME77+NixY894zg+vExoayvLly8+5NpH2pq574IdjgAzDOKtlKUTakuqTAlCwRd9nOXdqRxTxE3UrQQcd3wy19phBqUPjgKT9O3kAtAK9eIMCkIifqDmpBSjMGkRMuBWA/UcrfFmWiFdoCrx4m75JIn7i5HWAALrEhgGwr0gBSNo/7QQv3qZvkoifONEFVts9kB4bDkBOYbnPahLxFs9O8JoCL16ib5KInzjRBVbXAlQbgPYVKQBJ+1c3C8wWrBlg4h0KQCJ+wunZDLW2BahrXG0X2F51gYkfUAuQeJu+SSJ+4sRmqGoBEv/jCUAaAyReom+SiJ84pQXo+CDofLuDimpNhZf2zVFTOwtMAUi8Rd8kET9R84NZYB3CrJ71gDQTTNo7zxggBSDxEn2TRPxE3SywYPOJReK6eqbCqxtM2jd1gYm36Zsk4id+uA4QQNe42nFA7325n5XZBT6pS8QbPCtBaxC0eIlP9wITEe/xrAN00j5JvRIjAVj1/RFWfX+EF68dyE+GpTZ6Dbfb4B+bD5CdV4rZbOJHveLJKSzn39sOc6zCSWiwmV6JkSRE2ugWH8Fl/RIJtwXhdLkpq6qhQ1iwtimQFuFwahq8eJcCkIif8KwDZD7xL+RbRnbBYjbxVe4xln+Tz8OLt7FlfzHVNW7ySx2M7hHL9IvSsQaZMQyDxz7Yzjvrcj3Pf331nlNeZ3NusefPocEWgiwmSqtqB1kP7BzNzEt6UFhWTUV1DYlRIUzqn1SvVUrkbKgFSLxNAUjET/xwFhhAVEgwd/6oO263wYyFm/l4ex7vrj8RcFZ/f4TXVu0h2GKistqFvaoGkwluHtEFR42L/+zIJzIkiGkju9IzMZKSSie7C8ooKnewZmdh7RpDzhM1bD1Qwh1/2VSvrhsvSOO5q/urZUjOicYAibcpAIn4Cc86QJZTg4bZbOK31w/mR70OcrikiiCzCVuwmddX76GwrNpzntVi5pmr+3Pd8W6yF07zeoZhsLOgDIvZREyYlWqXm+eXfcemfcfoGhtOZEgQH207zHtf5lJU5iAk2EKV04XFXLtZa7jNQvf4CCYP7kR0WPA5v/+tB4r5MucoldUuhqfH0LljKF/lFlPldJEQFcKYHnGYzQph7VXdNHjNAhNvUQAS8ROeafDmhn8gQoIt3HBBWr1jU0d04fv8UoItZkKCLcRH2JocRkwmk2eMUZ2Xrhtc7/7A1bt57t/f8Z8d+Y1eJ/PjbxnXN5GR3WLpEBbM6B5xdAizNnp+jcvNuj1H2bD3KIZh0LljGCWVTjI//pbjGbBBAzpFM7RLR6xBZnomRHBxr3gSo0Ka9F7F96o1DV68TAFIxE801AV2JuG2IIakdWypkvifMd1Iiwln95EybEG1IcswDMqrXZRWOVnxbQHf5ZXy0dbDfLT1MAAdw4L5ybBUDhVXYq+qwWoxMbRLDJEhQXyfX8q/tx2u12p1stE94ugYbiUru4ByRw0DOkXTIczKpn3H2HawhG0HSzznWoPMXDOkE11iw7mkTzx9kqJa7HOQc6cuMPE2BSARP3GiC6zt/ECYTCYy+ic1+vgDE3qz7WAJ/96Wx878UnYfKWNvUcUpg68//bb+FP4OYcFc2ieBMKuF7QftfHvYzr3jenL32O6YTCacLjc1LoNQa+2MoaIyB4s27KfcUUNFtYvNucfYeqCERRv2A/Del7msfugSL7978SaH9gITL1MAEvETzuM/EEHtaJyLyWRiYOcODOzcAahtxXrvy1y2HyyhR0IEseE2jlVUs2nfMVxug/hIG+P7JTK6R9xpg16wxczJs6VjI2zMuKSH575hGHyxu4hl2/P4y7p97D9WgdPlblPhUerzdIEF67+ReIcCkIifcNatBN2Of8SDLWZuHdn1lOO3j/Hu65hMJkb1iGNkt1je+zKXGrdBYZmD5OhQ776QeI1D0+DFy/RNEvETJ/YCaz8tQL5mNpuIj7QBUGB3+LgaOZ0TY4C0EKJ4hwKQiB8wDKNNjgFqDxLqAlCpAlBbps1QxdvUBSbiB2pOmv8d3Mg0eGlYfGQIUEJBaZWvS/FbhmGw8MtcfvOf76msdtE3OZK5Nw6hc8ewJl+j+vg6QJoFJt6iACTiB+q6v0BdYM2VGFXbApTvx11ghmFQ5XSz6vsCXl+9h5BgCzdf2IW0mDC6x0d4Zss19Vpf5hwlLtJG9/iIeseXbDnE1gMllDtqyM4vpaTSidViprTKyaGSEwFzc24xd72zmffvHIktyIzTZVDtchNkNhHSyF5fmgYv3qYAJOIH6gZAgwJQcyVE1i6GeKSdtgBVOV04XW7MJhMuw2D/0Qp2FZSxu6CM3UfK2X2kjNyjFVRUu+o974vdRQDEhlt57ZahDO8ac8q1DcNg6dbDfJlzlPhIG8nRIWRlH+GjbbVrNg3vWruwZFJUKKVVztMueGm1mHkoozcj0mO55a31bDtYQv/Hl9drvQy3WsicMpArB6Wc8nyHApB4mQKQiB84uQVIXWDNkxDVcoOg3W6DUkcNFdU17CuqYG9hOfuPVVDldFNc4aSwzMHQLh35fwOTCbMGYa9yEhpsITWm8a6h0ionn36bT0mFk9U7C8nKLjjtCtgnC7damD6qKy43/GdHHoWlDorKq7npjXVcPaQT56d1xFHj5nBJFcUV1eQUlrM+5+gp17GYTbgNgw17j9U7HmQ2cdOINGLCrfRIiCAxKgSH002o1ULX2DBiI2o/69/dOIS73tlMmaOm3vPLq13c/7ctHCl10C0unDE94zwb6WolaPE2BSARP1BzfIqw2YT2u2omTxfYWbYAudwGZlPt1Poal5uj5dVUOl1s2neMF5Zlk2c//XVXfX+Elz75vt6x20anc/uYdFxug04dQjGZTBiGQVb2EX65eFu97qQf6hAWTI/4CHokRNA9PoLuCeF0iQ0nMSqEcKvFsyntLyb1oaK6hgfe/5p/b8vjbxsP8LeNB065ntVi5vrhqThqXJ595O67rBfRocGs33OUIIuJ7LxSdhaUcffY7gxroCXph8b0jGfdw+M8XWTWIHNtC9E/tvLh14d4eukOACael8irN51PkMXs2Q1eAUi8RQFIxA843XVT4PXj0Fx1XWDNbQHaf7SCx5ZsZ2X2EQBMJjAaaYkxm6BzxzDS48LpEhtGhC2IcFsQEbYglm49xNcHSnC5DSJsQZRUOnlzTQ5vrskBoHPHUPolR5Fnr2LrgRLPsUGdO5AaE8ZPhnWmU4fa9YsMg2aN5wmzBvHqTeezbs9R/r3tMAeOVWALspAUHUJchJVgi5nx/RLrjfU5WZfY8Ca/1g9FHH//J/v1TwaS0iGEHYfsrN9zlOXf5DP5D5+TGBnCvqIKAKwWTYMX71AAEvEDdS1AwWr9aba6afCFZQ5cbgNLI5/hgWMV/Hp5NnsKy8krqaKwzFGv66ku/JhNEBpsIcwWxPSLuvLTUV0JDT7R8vJD0y7qWu/+8m/y+OXibRwtr8ZsMnHgWCUHjlUCta0ft1zYhdkTehFm9c5f3yaTiZHdYxnZPdYr1zsXtiALcyb1BWo/h7vf3cz2g3a2Yz9xjlaCFi9RABLxA06XWoDOVmyEDbMJ3EbtnmEJJ+0QX1Fdw3++yafGbfCb/2Rz+AddTxekx/DklecRH2nD7TYwm010DLM2GqKaYuJ5SVzWNxGAqhoX6/YUcai4yvNY3cKN/m7ieUn8+94xfH2gmF99/B1Hy2s3wNVK0OItCkAifqDGsw2GWoCay2I2ERdho6DUQUHpiQBUXFHNtPkb+Hp/sefcHgkRPDSxN8nRoSRG2zzdZ95WN44rzBrEpX0SW+Q12oPeSZH0TooE4KG/bwXUAiTeowAk4gc822BoBthZSYiqDUBPfPAN/TtFA/Dx9sPk2x2eQcUpHUJ5/Ip+nplM0nquPb8zm/Yeo6y6hqSolgmdEnh8+rdlZmYmw4cPJzIykoSEBCZPnkx2dna9c6qqqpgxYwaxsbFEREQwZcoU8vMbX2sCateueOyxx0hOTiY0NJTx48ezc+fOlnwrIj7lPD4GSGsAnZ26lpyN+46x4Iu9LPhiL/l2BynRIbz/vyP5+10XMffGIQo/PmI2m3j+2oG8etP5jY6lEmkunwagVatWMWPGDNatW8cnn3yC0+lkwoQJlJeXe8657777+PDDD3n//fdZtWoVhw4d4pprrjntdV944QXmzp3La6+9xvr16wkPD2fixIlUVbXPhc5EzkT7gJ2bAcdbffolR3Hnj7oz/aKuvHbzUP77wFh6Jkb6uDoRaQkmw2hs4mbrO3LkCAkJCaxatYqLL76YkpIS4uPjWbhwIddeey0A3333HX379mXt2rVceOGFp1zDMAxSUlK4//77eeCBBwAoKSkhMTGRBQsWcMMNN5yxDrvdTnR0NCUlJURFRXn3TYq0gC92F3LTG+vpmRDBJ7N/5Oty2h232yD3aAVdYsPUwiDSjjXn97tN/XOxpKR2jYuYmNqFtDZt2oTT6WT8+PGec/r06UNaWhpr165t8Bo5OTnk5eXVe050dDQjRoxo9DkOhwO73V7vJtKe1GgW2Dkxm010jQtX+BEJIG3mb0u3282sWbMYNWoU/fv3ByAvLw+r1UqHDh3qnZuYmEheXl6D16k7nphYf+bE6Z6TmZlJdHS055aamnqO70akdWkWmIhI87SZADRjxgy2b9/OokWLWv2158yZQ0lJiee2f//+Vq9B5Fx41gHSQogiIk3SJgLQzJkzWbp0KStXrqRz586e40lJSVRXV1NcXFzv/Pz8fJKSkhq8Vt3xH84UO91zbDYbUVFR9W4i7Ym6wEREmsenf1sahsHMmTNZvHgx//3vf0lPT6/3+NChQwkODmbFihWeY9nZ2eTm5jJy5MgGr5menk5SUlK959jtdtavX9/oc0TaO3WBiYg0j08D0IwZM3jnnXdYuHAhkZGR5OXlkZeXR2Vl7b430dHR3HbbbcyePZuVK1eyadMmfvrTnzJy5Mh6M8D69OnD4sWLgdp9bWbNmsUzzzzDBx98wLZt27j11ltJSUlh8uTJvnibIi3OqYUQRUSaxacrQc+bNw+AsWPH1js+f/58pk+fDsBvf/tbzGYzU6ZMweFwMHHiRP7whz/UOz87O9szgwzgoYceory8nDvuuIPi4mJGjx7NsmXLCAnRCqLin+oWQlQLkIhI07SpdYDaCq0DJO3NX9bu5dEl35BxXhKv3TLU1+WIiPhEu10HSETOzond4NUCJCLSFApAIn6gbhC0VbPARESaRH9bivgBtQCJiDSPApCIH9A6QCIizaO/LUX8gGcdIK0ELSLSJApAIn7AqRYgEZFm0d+WIn6g5vg6QBoDJCLSNApAIn6gxl3bAhSslaBFRJpEf1uK+AGnWoBERJpFAUjED9TNAgvWGCARkSbx6V5gItI8FdU17CuqINhi4khpNUfKHFTXuFm+Iw+AIM0CExFpEgUgkTaq3FHD57sK2X+skrySSnIKK1iz6whVTnejz9EsMBGRplEAEmkjXG6DPHsV3+eVsmTLQZZ/k0+l03XKedGhwQB0DAsmMSoEW7CFI6UOnC43P+oV39pli4i0SwpAIj52qLiS11fvYcmWgxyrcNZ7rEtsGAM6RZMcHUJSdCgXdI2hf6coTCZ1dYmInAsFIBEv21VQRkllNUNSO2I2mzAMg+0H7Ww9WMy+ogpKKpxYg8x0DAvGUePmnXX7KK+ubekJtphIjg7lkt7xTB7SicGpHRR2RERagAKQyDkqqXTy+a5CvjlUwtYDJXy2sxCAHgkRDO/akW8Pl7Jlf/Fpr3F+WgfuGdeTMT3iNI5HRKQVKACJNIFhGNS4DSocLj74+iCffluA2zDIt1exq6CM4+sQAmA2QWiwhV0FZewqKAPAGmRmZLdYusdHEBMeTHWNm2MVTmrcBkNSO3Dt0M6YNYNLRKTVKACJnMb2gyXcu+gr9hwpP+153ePDGdEtlq6xYUzol0RshJWPt+WRZ68iNNjC1ed3Ii7C1kpVi4jImSgAScDbtO8Yn+8qxOlyU+1yU+MyqDn+5w+2HPKMz6nTNTaMm0akERtuo0NYMP07RZMYFXLKda8bntpab0FERJpJAUgC0v6jFew4bOer3GL+uHo3htH4uaN6xPLCtYOIsAVhtZgJCTZrYLKISDunACQBxVHj4qVPvufNz3I8G4gCXNYvkZToEIItZoIsZoItJoItZhKjbFw9pDPWIA1MFhHxJwpAElBe/nQnf1y1B4A+SZF0DLNy04g0rhiU4uPKRESkNSkAScAod9Twzrp9AGReM4AbL0jzcUUiIuIrateXgPH+xv2UVtWQHhfO9cM0QFlEJJCpBUj8Wr69ij99toePt+dRUOoA4GejumrNHRGRAKcAJH6rsMzBj1/5jKLyas+xpKgQpgzt7MOqRESkLVAAEr/11w37KSqvpktsGI9c3o+UDiGkxYQRZtXXXkQk0OmXQPxSjcvNu8cHPN97aU8u65fo44pERKQtUQASv/PpjnzW7CrkUEkVMeFWLh+Y7OuSRESkjVEAEr+ydncRt/95o+f+dcNSCQm2+LAiERFpixSAxK+89XkOAAM6RTOsa0fuGtvdxxWJiEhbpAAkfiO3qIJPv80H4LfXD6ZHQoSPKxIRkbZKCyGK35j/RQ6GARf3ilf4ERGR0/JpAFq9ejVXXHEFKSkpmEwm/vWvf9V73GQyNXh78cUXG73mE088ccr5ffr0aeF3Ir629UAxf15bO+vrttHpPq5GRETaOp92gZWXlzNo0CB+9rOfcc0115zy+OHDh+vd//jjj7ntttuYMmXKaa973nnn8emnn3ruBwWpp89frd9TxAdfH2LNrkJcboPLBybzo17xvi5LRETaOJ8mg0mTJjFp0qRGH09KSqp3f8mSJVxyySV069bttNcNCgo65bmn43A4cDgcnvt2u73JzxXfcbsN7l30Ffn22v92cRFWnr6qv4+rEhGR9qDdjAHKz8/no48+4rbbbjvjuTt37iQlJYVu3boxdepUcnNzT3t+ZmYm0dHRnltqqjbKbA825x4j3+4gwhbEfeN78ZfbRhATbvV1WSIi0g60mwD09ttvExkZ2WBX2clGjBjBggULWLZsGfPmzSMnJ4cxY8ZQWlra6HPmzJlDSUmJ57Z//35vly8t4OPteQCM75vAz8f3pG9ylI8rEhGR9qLdDI556623mDp1KiEhIac97+QutYEDBzJixAi6dOnC3/72t0Zbj2w2Gzabzav1SssyDINlxwNQRn+t9CwiIs3TLgLQZ599RnZ2Nn/961+b/dwOHTrQq1cvdu3a1QKVSWtz1Lh47qNvWfFdAQeLKwmzWhjbW4OeRUSkedpFF9ibb77J0KFDGTRoULOfW1ZWxu7du0lOVitBe1fmqOHWN7/k7bX7OHCsEoD/NzBZW12IiEiz+bQFqKysrF7LTE5ODlu2bCEmJoa0tDSgdkbW+++/z29+85sGrzFu3DiuvvpqZs6cCcADDzzAFVdcQZcuXTh06BCPP/44FouFG2+8seXfkLSot9bksD7nKBG2IJ67ZgDpseH0SY70dVkiItIO+TQAbdy4kUsuucRzf/bs2QBMmzaNBQsWALBo0SIMw2g0wOzevZvCwkLP/QMHDnDjjTdSVFREfHw8o0ePZt26dcTHq5ukPTMMgw++PgTAY1f048pBKT6uSERE2jOTYRiGr4toa+x2O9HR0ZSUlBAVpZlFbUF2XikTX16N1WJmwyPjiQ4N9nVJIiLSxjTn97tdjAGSwOZ0ufnweOvPxb3iFH5EROSctYtZYBK4lm0/zD3vfYXTVdtQ+f8GqutLRETOnVqApE17+4t9nvCTFBXC+H6JPq5IRET8gVqApM0qKnOwPqcIgA9mjqJ3UiS2IE15FxGRc6cAJG3Wp9/m4zbgvJQoBnbu4OtyRETEj6gLTNqsur2+JvVP8nElIiLibxSApE0qsFfx+a7a9Z0yFIBERMTLFICkTXp15S6cLoPz0zrQI0GrPYuIiHdpDJC0Kd8cKmFnfhnvfbkfgPsn9PZxRSIi4o8UgKTN2H+0gsmvfu6Z9n5htxgu6h7r46pERMQfKQBJm/HB14dwugxiw610T4jgiSvPw2Qy+bosERHxQwpA0mbUbXfxUEZvrh+e5uNqRETEn2kQtLQJuwpK+S6vlCCziYnnadaXiIi0LLUAic8U2KvIyj5Cjdtg6dba1p8xPePoEGb1cWUiIuLvFIDEZ+5+dzMb9x2rd+zaoak+qkZERAKJApD4xPaDJWzcd4xgi4kxPePpFhfOxP5JDO8a4+vSREQkACgAiU+8uz4XgInnJfH7m873cTUiIhJoFICkVS3/Jo+lWw/z6Y58AG6+sIuPKxIRkUCkACStJt9exaxFW6h0ugDolRjBiHR1eYmISOtTAJJW8+vl2VQ6XZyXEsU153dmXJ8ELXQoIiI+oQAkreK7PDt/33wAgKcn9+f8tI4+rkhERAKZFkKUVvHe+lwMAzLOS1L4ERERn1MAkhbndLn5cOthAG64QOv8iIiI7ykASYtb/f0RjpZXExdhY3SPOF+XIyIiojFA0nIK7FW8sDybzcdXe75yUApBFmVuERHxPQUgaTHPL8vmH8cHPgNMGdrJh9WIiIicoAAkLeJgcSVLthwEYPZlvRjYOZrzUqJ9XJWIiEgtBSBpEW+s3kON2+Ci7rHcO66nr8sRERGpRwFIvGbTvqN8e7iUbw/bPXt93TW2u4+rEhEROZUCkHjFweJKbnh9HU6X4Tn2s1HpmvUlIiJtkgKQeMV763Nxugw6dQhlQKdopl6Yxpie8b4uS0REpEEKQHLOqmvcLNqwH4CHf9yXywcm+7giERGR0/PpoiyrV6/miiuuICUlBZPJxL/+9a96j0+fPh2TyVTvlpGRccbrvvrqq3Tt2pWQkBBGjBjBl19+2ULvILC53QbzP8/hjr9spLDMQXykjQnnJfq6LBERkTPyaQAqLy9n0KBBvPrqq42ek5GRweHDhz23995777TX/Otf/8rs2bN5/PHH2bx5M4MGDWLixIkUFBR4u/yA95tPsnnywx1kZR8B4KYL0gjWQociItIO+LQLbNKkSUyaNOm059hsNpKSkpp8zZdeeon/+Z//4ac//SkAr732Gh999BFvvfUWv/jFL86pXjlhyZaDvLpyNwB3XNyNngkRXDEoxcdViYiINE2b/+d6VlYWCQkJ9O7dm7vuuouioqJGz62urmbTpk2MHz/ec8xsNjN+/HjWrl3b6PMcDgd2u73eTRq3t7Cch/+5Daid5v7wj/vyk2GphARbfFyZiIhI07TpQdAZGRlcc801pKens3v3bh5++GEmTZrE2rVrsVhO/bEtLCzE5XKRmFh/HEpiYiLfffddo6+TmZnJk08+6fX627PCMgcvLssmO78Uq8XMFYNTcDhd7CooY8Peo5RXu7ggPYYHJvT2dakiIiLN1qYD0A033OD584ABAxg4cCDdu3cnKyuLcePGee115syZw+zZsz337XY7qampXrt+e+N2G9z31y18trPQc+zLvUfrnRMVEsTL1w/GYja1dnkiIiLnrE0HoB/q1q0bcXFx7Nq1q8EAFBcXh8ViIT8/v97x/Pz8044jstls2Gw2r9fbnlQ5XRTYHdirnCzZcpDPdhZiCzLzwrUDOVLqYOnWw3QMC2ZQagdsQRYu7hVHSodQX5ctIiJyVtpVADpw4ABFRUUkJze8zozVamXo0KGsWLGCyZMnA+B2u1mxYgUzZ85sxUobdrC4ktyiilOOGxgNnA1uN1RU12AymbAGmbEen2FV43bjdLkBCDKbCbaYMZvAZRi43XX/a+ByG9S4DZwuN9U1bqpdtc9zugxqXG5qjp9zuKSKpVsPUVpVU+/1H7m8L1cNrt3B/fYx3bz5UYiIiPiUTwNQWVkZu3bt8tzPyclhy5YtxMTEEBMTw5NPPsmUKVNISkpi9+7dPPTQQ/To0YOJEyd6njNu3DiuvvpqT8CZPXs206ZNY9iwYVxwwQW8/PLLlJeXe2aF+dKSLQd5YVm2r8tolC3ITLgtiL7JkWScl8TNF3bxdUkiIiItwqcBaOPGjVxyySWe+3XjcKZNm8a8efPYunUrb7/9NsXFxaSkpDBhwgSefvrpet1Vu3fvprDwxFiV66+/niNHjvDYY4+Rl5fH4MGDWbZs2SkDo32hY5iVngkRTT7fbDIRYrVgona1ZUeNC5PJRJC5tkUI8LTmuAwDi8mExWzCXPe/ZhMWE7WtR0EWrBYz1iATwRYzFnPtdSxmMyHBZsb1SeSi7rGYNaZHREQCgMkwjIb7XwKY3W4nOjqakpISoqKifF2OiIiINEFzfr/b/DpAIiIiIt6mACQiIiIBRwFIREREAo4CkIiIiAQcBSAREREJOApAIiIiEnAUgERERCTgKACJiIhIwFEAEhERkYCjACQiIiIBRwFIREREAo4CkIiIiAQcBSAREREJOApAIiIiEnCCfF1AW2QYBgB2u93HlYiIiEhT1f1u1/2On44CUANKS0sBSE1N9XElIiIi0lylpaVER0ef9hyT0ZSYFGDcbjeHDh0iMjISk8nk1Wvb7XZSU1PZv38/UVFRXr22nEqfd+vS59269Hm3Pn3mrau5n7dhGJSWlpKSkoLZfPpRPmoBaoDZbKZz584t+hpRUVH6P08r0ufduvR5ty593q1Pn3nras7nfaaWnzoaBC0iIiIBRwFIREREAo4CUCuz2Ww8/vjj2Gw2X5cSEPR5ty593q1Ln3fr02feulry89YgaBEREQk4agESERGRgKMAJCIiIgFHAUhEREQCjgKQiIiIBBwFoFb06quv0rVrV0JCQhgxYgRffvmlr0vyS5mZmQwfPpzIyEgSEhKYPHky2dnZvi4rYPzqV7/CZDIxa9YsX5fi1w4ePMjNN99MbGwsoaGhDBgwgI0bN/q6LL/kcrl49NFHSU9PJzQ0lO7du/P00083ab8paZrVq1dzxRVXkJKSgslk4l//+le9xw3D4LHHHiM5OZnQ0FDGjx/Pzp07z+k1FYBayV//+ldmz57N448/zubNmxk0aBATJ06koKDA16X5nVWrVjFjxgzWrVvHJ598gtPpZMKECZSXl/u6NL+3YcMG/vjHPzJw4EBfl+LXjh07xqhRowgODubjjz9mx44d/OY3v6Fjx46+Ls0vPf/888ybN4/f//73fPvttzz//PO88MIL/O53v/N1aX6jvLycQYMG8eqrrzb4+AsvvMDcuXN57bXXWL9+PeHh4UycOJGqqqqzf1FDWsUFF1xgzJgxw3Pf5XIZKSkpRmZmpg+rCgwFBQUGYKxatcrXpfi10tJSo2fPnsYnn3xi/OhHPzJ+/vOf+7okv/V///d/xujRo31dRsC4/PLLjZ/97Gf1jl1zzTXG1KlTfVSRfwOMxYsXe+673W4jKSnJePHFFz3HiouLDZvNZrz33ntn/TpqAWoF1dXVbNq0ifHjx3uOmc1mxo8fz9q1a31YWWAoKSkBICYmxseV+LcZM2Zw+eWX1/ueS8v44IMPGDZsGD/5yU9ISEhgyJAhvPHGG74uy29ddNFFrFixgu+//x6Ar7/+mjVr1jBp0iQfVxYYcnJyyMvLq/d3S3R0NCNGjDin31BthtoKCgsLcblcJCYm1juemJjId99956OqAoPb7WbWrFmMGjWK/v37+7ocv7Vo0SI2b97Mhg0bfF1KQNizZw/z5s1j9uzZPPzww2zYsIF7770Xq9XKtGnTfF2e3/nFL36B3W6nT58+WCwWXC4Xzz77LFOnTvV1aQEhLy8PoMHf0LrHzoYCkPi1GTNmsH37dtasWePrUvzW/v37+fnPf84nn3xCSEiIr8sJCG63m2HDhvHcc88BMGTIELZv385rr72mANQC/va3v/Huu++ycOFCzjvvPLZs2cKsWbNISUnR592OqQusFcTFxWGxWMjPz693PD8/n6SkJB9V5f9mzpzJ0qVLWblyJZ07d/Z1OX5r06ZNFBQUcP755xMUFERQUBCrVq1i7ty5BAUF4XK5fF2i30lOTqZfv371jvXt25fc3FwfVeTfHnzwQX7xi19www03MGDAAG655Rbuu+8+MjMzfV1aQKj7nfT2b6gCUCuwWq0MHTqUFStWeI653W5WrFjByJEjfViZfzIMg5kzZ7J48WL++9//kp6e7uuS/Nq4cePYtm0bW7Zs8dyGDRvG1KlT2bJlCxaLxdcl+p1Ro0adsrTD999/T5cuXXxUkX+rqKjAbK7/c2mxWHC73T6qKLCkp6eTlJRU7zfUbrezfv36c/oNVRdYK5k9ezbTpk1j2LBhXHDBBbz88suUl5fz05/+1Nel+Z0ZM2awcOFClixZQmRkpKePODo6mtDQUB9X538iIyNPGV8VHh5ObGysxl21kPvuu4+LLrqI5557juuuu44vv/yS119/nddff93XpfmlK664gmeffZa0tDTOO+88vvrqK1566SV+9rOf+bo0v1FWVsauXbs893NyctiyZQsxMTGkpaUxa9YsnnnmGXr27El6ejqPPvooKSkpTJ48+exf9Bxmqkkz/e53vzPS0tIMq9VqXHDBBca6det8XZJfAhq8zZ8/39elBQxNg295H374odG/f3/DZrMZffr0MV5//XVfl+S37Ha78fOf/9xIS0szQkJCjG7duhm//OUvDYfD4evS/MbKlSsb/Ht72rRphmHUToV/9NFHjcTERMNmsxnjxo0zsrOzz+k1TYahpSxFREQksGgMkIiIiAQcBSAREREJOApAIiIiEnAUgERERCTgKACJiIhIwFEAEhERkYCjACQiIiIBRwFIREREAo4CkIi0aVlZWZhMJoqLi33y+itWrKBv375N2tR12bJlDB48WHtEibQDCkAi0maMHTuWWbNm1Tt20UUXcfjwYaKjo31S00MPPcQjjzzSpE1dMzIyCA4O5t13322FykTkXCgAiUibZrVaSUpKwmQytfprr1mzht27dzNlypQmP2f69OnMnTu3BasSEW9QABKRNmH69OmsWrWKV155BZPJhMlkYu/evad0gS1YsIAOHTqwdOlSevfuTVhYGNdeey0VFRW8/fbbdO3alY4dO3LvvffW67ZyOBw88MADdOrUifDwcEaMGEFWVtZpa1q0aBGXXXYZISEhnmNff/01l1xyCZGRkURFRTF06FA2btzoefyKK65g48aN7N6926ufj4h4V5CvCxARAXjllVf4/vvv6d+/P0899RQA8fHx7N2795RzKyoqmDt3LosWLaK0tJRrrrmGq6++mg4dOvDvf/+bPXv2MGXKFEaNGsX1118PwMyZM9mxYweLFi0iJSWFxYsXk5GRwbZt2+jZs2eDNX322WfcdNNN9Y5NnTqVIUOGMG/ePCwWC1u2bCE4ONjzeFpaGomJiXz22Wd0797dS5+OiHibApCItAnR0dFYrVbCwsJISko67blOp5N58+Z5Asa1117LX/7yF/Lz84mIiKBfv35ccsklrFy5kuuvv57c3Fzmz59Pbm4uKSkpADzwwAMsW7aM+fPn89xzzzX4Ovv27fOcXyc3N5cHH3yQPn36ADQYnlJSUti3b1+zPwMRaT0KQCLS7oSFhdVrXUlMTKRr165ERETUO1ZQUADAtm3bcLlc9OrVq951HA4HsbGxjb5OZWVlve4vgNmzZ3P77bfzl7/8hfHjx/OTn/zklJae0NBQKioqzvr9iUjLUwASkXbn5C4nAJPJ1OCxuunoZWVlWCwWNm3adMpsrpND0w/FxcVx7NixeseeeOIJbrrpJj766CM+/vhjHn/8cRYtWsTVV1/tOefo0aPEx8ef1XsTkdahACQibYbVam3SejvNNWTIEFwuFwUFBYwZM6ZZz9uxY8cpx3v16kWvXr247777uPHGG5k/f74nAFVVVbF7926GDBnitfpFxPs0C0xE2oyuXbuyfv169u7dS2FhodcWFOzVqxdTp07l1ltv5Z///Cc5OTl8+eWXZGZm8tFHHzX6vIkTJ7JmzRrP/crKSmbOnElWVhb79u3j888/Z8OGDfTt29dzzrp167DZbIwcOdIrtYtIy1AAEpE244EHHsBisdCvXz/i4+PJzc312rXnz5/Prbfeyv3330/v3r2ZPHkyGzZsIC0trdHnTJ06lW+++Ybs7GwALBYLRUVF3HrrrfTq1YvrrruOSZMm8eSTT3qe89577zF16lTCwsK8VruIeJ/JMAzD10WIiLRVDz74IHa7nT/+8Y9nPLewsJDevXuzceNG0tPTW6E6ETlbagESETmNX/7yl3Tp0qVJ3XF79+7lD3/4g8KPSDugFiAREREJOGoBEhERkYCjACQiIiIBRwFIREREAo4CkIiIiAQcBSAREREJOApAIiIiEnAUgERERCTgKACJiIhIwFEAEhERkYDz/wF+kHA+Qz61DwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(times, all_ssim)\n",
    "plt.xlabel(\"time (s)\")\n",
    "plt.ylabel(\"ssim\")\n",
    "plt.show()\n",
    "\n",
    "plt.plot(times, all_psnr)\n",
    "plt.xlabel(\"time (s)\")\n",
    "plt.ylabel(\"psnr (db)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f2f22cd-b524-43c1-ac3d-c14661282896",
   "metadata": {},
   "source": [
    "## Compare performances on syntetic data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "4841e24c-d794-4fe9-9a62-e5ebefefc3e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Functions as implemented in classical libraries.\n",
    "\n",
    "def naive_psnr(im1: np.ndarray, im2: np.ndarray, weights = None) -> float:\n",
    "    \"\"\"Compute the psnr with numpy.\"\"\"\n",
    "    if weights is None:\n",
    "      weights = [1.0 for _ in range(im1.shape[2])]\n",
    "    layers_mse = ((im1 - im2)**2).mean(axis=(0, 1)).tolist()\n",
    "    tot = sum(weights)\n",
    "    mse = sum(l*w/tot for w, l in zip(weights, layers_mse))\n",
    "    return -10.0*math.log10(mse) if mse > 1e-10 else 100.0\n",
    "\n",
    "def naive_ssim(\n",
    "  im1: np.ndarray, im2: np.ndarray, data_range : float = 1.0, weights = None, sigma: float = 1.5\n",
    ") -> float:\n",
    "    \"\"\"Compute the ssim with numpy and cv2.\"\"\"\n",
    "    # get gaussian window\n",
    "    r = int(3.5 * sigma + 0.5)  # same as skimage.metrics.structural_similarity\n",
    "    gauss = np.exp(-(np.arange(-r, r+1)**2) / (2.0 * sigma**2))\n",
    "    gauss_i, gauss_j = np.meshgrid(gauss, gauss, indexing=\"ij\")\n",
    "    gauss = gauss_i * gauss_j\n",
    "    gauss /= gauss.sum()\n",
    "    # compute statistics for all patches\n",
    "    mu1 = cv2.filter2D(im1, ddepth=-1, kernel=gauss)\n",
    "    mu2 = cv2.filter2D(im2, ddepth=-1, kernel=gauss)\n",
    "    mu11, mu22, mu12 = mu1 * mu1, mu2 * mu2, mu1 * mu2\n",
    "    s11 = cv2.filter2D(im1*im1, ddepth=-1, kernel=gauss) - mu11\n",
    "    s22 = cv2.filter2D(im2*im2, ddepth=-1, kernel=gauss) - mu22\n",
    "    s12 = cv2.filter2D(im1*im2, ddepth=-1, kernel=gauss) - mu12\n",
    "    # crop patches\n",
    "    mu11, mu22, mu12 = mu11[r:-r, r:-r], mu22[r:-r, r:-r], mu12[r:-r, r:-r]\n",
    "    s11, s22, s12 = s11[r:-r, r:-r], s22[r:-r, r:-r], s12[r:-r, r:-r]\n",
    "    # ssim formula\n",
    "    c1, c2 = (0.01 * data_range)**2, (0.03 * data_range)**2\n",
    "    ssim = ((2.0*mu12 + c1) * (2.0*s12 + c2)) / ((mu11 + mu22 + c1) * (s11 + s22 + c2))\n",
    "    # average\n",
    "    if weights is None:\n",
    "      weights = [1.0 for _ in range(im1.shape[2])]\n",
    "    weights = np.asarray(weights, dtype=im1.dtype)\n",
    "    return float((ssim.mean(axis=(0, 1)) * weights).sum() / weights.sum())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e61e8dd0-6a60-40d8-a451-987b5e1dcbc3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create syntetic 4k images.\n",
    "\n",
    "shape = (2160, 3840, 3)  # 4k\n",
    "# shape = (1080, 1920, 3)  # full hd\n",
    "im1 = np.random.random((2160, 3840, 3))\n",
    "im2 = 0.8*im1 + 0.2*np.random.random((2160, 3840, 3))\n",
    "im1, im2 = im1.astype(np.float32), im2.astype(np.float32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "af9f22f5-d4db-4865-878a-946c05b0f561",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PSNR: 120.8ms vs 4.9ms\n"
     ]
    }
   ],
   "source": [
    "time_naive_psnr = min(timeit.repeat(lambda: naive_psnr(im1, im2), repeat=5, number=10)) / 10\n",
    "time_c_psnr = min(timeit.repeat(lambda: compute_psnr(im1, im2), repeat=5, number=10)) / 10\n",
    "print(f\"PSNR: {1000*time_naive_psnr:.1f}ms vs {1000*time_c_psnr:.1f}ms\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "02e29fdc-1afa-47fa-8039-7c435cb59167",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SSIM: 1111.5ms vs 502.3ms\n"
     ]
    }
   ],
   "source": [
    "time_naive_ssim = min(timeit.repeat(lambda: naive_ssim(im1, im2, data_range=1.0), repeat=3, number=10)) / 10\n",
    "time_c_ssim = min(timeit.repeat(lambda: compute_ssim(im1, im2, data_range=1.0), repeat=3, number=10)) / 10\n",
    "print(f\"SSIM: {1000*time_naive_ssim:.1f}ms vs {1000*time_c_ssim:.1f}ms\")"
   ]
  }
 ],
 "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.13.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}