{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Day 11, possibly dynamic programming?\n",
    "\n",
    "- [Day 11](https://adventofcode.com/2018/day/11)\n",
    "\n",
    "This is a challenge that, to me, sounds like we need to use dynamic programming. For a large problem set, you'd only want to keep a running total. Since we want the top-left grid coordinate, you'd start at `(max(x), max(y))` and work your way backwards, taking advantage of the calculations already done for the `x + 1, y`, `x + 2, y`, `x, y + 1`, ..., `x + 2, y + 2` positions. For a 300 x 300 grid that would mean you only need to keep the last 600 calculation results in memory and let you use `max()` on the running calculation.\n",
    "\n",
    "However, for a 300x300 grid it is simpler to just vectorise the hell out of the grid, using `numpy`.\n",
    "\n",
    "The formula for any given grid coordinate power value is\n",
    "\n",
    "$$\\lfloor\\frac{((x + 10)y + serial) \\times (x + 10)}{100}\\rfloor \\bmod 10 - 5$$\n",
    "\n",
    "but I must note that subtracting 5 at the end doesn't actually matter to the outcome. Either the cell score falls in the range $[0, 10)$ or $[-5, 5)$, with the 3 x 3 grid score in the range $[0, 81]$ or $[-45, 36]$. Not that `numpy` much cares.\n",
    "\n",
    "Summing the sliding 3x3 windows is a little more interesting here. There are 298 x 298 complete 3 x 3 sub-grids that need to be considered here (from `((1, 1) ... (3, 3))` all the way to `((298, 298) ... (300, 300))`), and we need to create sums for all those sub windows. I'm using the [`numpy.lib.stride_tricks.as_strided()` function](https://docs.scipy.org/doc/numpy/reference/generated/numpy.lib.stride_tricks.as_strided.html) here to step over the whole matrix in those 3x3 windows, so we can sum them all and produce a new `(298 x 298)` matrix of sums at coordinates that match the top-level corner of each sub-matrix.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.lib.stride_tricks import as_strided\n",
    "\n",
    "# These values never change, so can be made globals\n",
    "GRIDSIZE = 300\n",
    "XX, YY = np.meshgrid(np.arange(1, GRIDSIZE + 1), np.arange(1, GRIDSIZE + 1))\n",
    "RACK_ID = XX + 10\n",
    "\n",
    "\n",
    "def grid_powers(serial):\n",
    "    # calculate power levels\n",
    "    return (RACK_ID * YY + serial) * RACK_ID // 100 % 10 - 5\n",
    "\n",
    "\n",
    "def summed_grid_powers(power_levels, window_size=3):\n",
    "    # sum levels for 3 x 3 subgrids; substitute edges for zeros\n",
    "\n",
    "    window_count = GRIDSIZE - window_size + 1\n",
    "    # output shape, 2d grid of 2d windows\n",
    "    shape = (window_count, window_count, window_size, window_size)\n",
    "    # per shape axis, the stride across power_levels matches up to the\n",
    "    # same axes.\n",
    "    strides = power_levels.strides * 2\n",
    "\n",
    "    # we want to sum every subwindow, so it is time to start striding\n",
    "    # we need to produce a (window_count, window_count, window_size, window_size)\n",
    "    # matrix that then is summed on the last 2 axes.\n",
    "    return as_strided(power_levels, shape, strides).sum(axis=(2, 3))\n",
    "\n",
    "\n",
    "def max_grid(serial):\n",
    "    summed = summed_grid_powers(grid_powers(serial))\n",
    "\n",
    "    # produce the (x, y) coordinates for the largest 3x3 grid top-left coordinate\n",
    "    # argmax() flattens the array and gives us an index based on that, so we need\n",
    "    # numpy.unravel to give back the original y, x coordinates.\n",
    "    y, x = np.unravel_index(summed.argmax(), summed.shape)\n",
    "    # Translate from zero to one-based indexing\n",
    "    return x + 1, y + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "power_tests = {\n",
    "    # serial, x, y: power level\n",
    "    (8, 3, 5): 4,\n",
    "    (57, 122, 79): -5,\n",
    "    (39, 217, 196): 0,\n",
    "    (71, 101, 153): 4,\n",
    "}\n",
    "\n",
    "for (tserial, x, y), expected in power_tests.items():\n",
    "    # indexing a [y, x] arranged matrix with 0-based offsets\n",
    "    assert grid_powers(tserial)[y - 1, x - 1] == expected\n",
    "\n",
    "max_tests = {\n",
    "    18: (33, 45),\n",
    "    42: (21, 61),\n",
    "}\n",
    "\n",
    "for tserial, expected in max_tests.items():\n",
    "    assert max_grid(tserial) == expected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import aocd\n",
    "\n",
    "data = aocd.get_data(day=11, year=2018)\n",
    "serial = int(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Part 1: 44,37\n"
     ]
    }
   ],
   "source": [
    "x, y = max_grid(serial)\n",
    "print(f\"Part 1: {x},{y}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "367 µs ± 11.7 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)\n",
      "2 ms ± 10.9 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit grid_powers(serial)\n",
    "%timeit max_grid(serial)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2, variable window size\n",
    "\n",
    "Now I'm really glad I worked out how to do striding here. We can now produce separate summed matrices for each size. Vectorised sums over striding views does slow down dramatically in the middle somewhere, see timings below.\n",
    "\n",
    "We still may want to explore dynamic programming here, however, as that would let us calculate values for all possible sizes as you traverse from 300, 300 down to 1,1 in one sweep. That's for later however.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "357 µs ± 4.49 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)\n",
      "1.64 ms ± 23.9 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)\n",
      "57.3 ms ± 147 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "12.6 µs ± 113 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit grid_powers(serial)\n",
    "power_levels = grid_powers(serial)\n",
    "%timeit summed_grid_powers(power_levels)\n",
    "%timeit summed_grid_powers(power_levels, 150)  # performance degrades\n",
    "%timeit summed_grid_powers(power_levels, 300)  # performance degrades"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def optimal_window_size(serial):\n",
    "    power_levels = grid_powers(serial)\n",
    "    # big matrix with 300 flattened power level arrays,\n",
    "    # all padded back to 300 x 300 size so we can determine which one\n",
    "    # has the most power output, then extract the size and grid position\n",
    "    by_size = np.stack(\n",
    "        [\n",
    "            np.pad(summed_grid_powers(power_levels, i + 1), (0, i), \"constant\")\n",
    "            for i in range(GRIDSIZE)\n",
    "        ]\n",
    "    ).reshape(GRIDSIZE, -1)\n",
    "    size = by_size.max(axis=1).argmax()\n",
    "    y, x = np.unravel_index(by_size[size].argmax(), power_levels.shape)\n",
    "    return x + 1, y + 1, size + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "optimal_tests = {\n",
    "    18: (90, 269, 16),\n",
    "    42: (232, 251, 12),\n",
    "}\n",
    "for tserial, expected in optimal_tests.items():\n",
    "    assert optimal_window_size(tserial) == expected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Part 2: 235,87,13\n"
     ]
    }
   ],
   "source": [
    "x, y, s = optimal_window_size(serial)\n",
    "print(f\"Part 2: {x},{y},{s}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAGzCAYAAAA1yP25AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABnRUlEQVR4nO3deXhTVf4/8PdNmqRrWrpTaCmUpaxFC61FlqJArQ4K6IjijAW+4shQR6eov3ZmBKrjMOqojILijqOgCDPAqCOCBWQHCxQVaNnKTne6pSVNk/P7o+RC6JZC06TJ+/U8fSA3N/eefLgh795zzr2SEEKAiIiIqBNS2LsBRERERDeKQYaIiIg6LQYZIiIi6rQYZIiIiKjTYpAhIiKiTotBhoiIiDotBhkiIiLqtBhkiIiIqNNikCEiIqJOi0GGyAFIkoTU1FR7NwMLFiyAJEn2bsYNmz59OiIjI61e19vb27YNIgANx/eCBQscbluOrLN/FjsSg0wntWzZMkiSJP+4u7ujb9++SE1NRWFhob2b53AuXLiABQsWICcnp0P2t2LFCixatKhD9kXNq6mpwYIFC7BlyxZ7N8Up7Ny5EwsWLEB5ebm9m0IkY5Dp5F544QV8+umnWLx4MUaMGIF33nkHCQkJqKmpsXfTHMqFCxeQmZnJINOKv/zlL6itrbV3M27Y+++/j7y8PPlxTU0NMjMzGWTayc6dO5GZmdnmIFNbW4u//OUvtmkUuTw3ezeAbk5ycjKGDRsGAHjssccQEBCA119/HevWrcPDDz9s59Y1z2Qyoa6uDu7u7vZuCl3Dzc0Nbm6d978FlUpl7ybQFdd+xvk5J1viGRknc8cddwAA8vPzAQD19fV48cUXERUVBY1Gg8jISPzpT3+CXq+XX5OWloaAgABceyP0J598EpIk4c0335SXFRYWQpIkvPPOO/IyvV6P+fPno3fv3tBoNAgPD8dzzz1nsX3g6hiQ5cuXY+DAgdBoNFi/fn2L7+Xtt9+W1w0LC8OcOXMa/SYYGRmJ6dOnN3ptYmIiEhMTAQBbtmzB8OHDAQAzZsyQu+OWLVsmrzto0CDs27cPI0aMgIeHB3r27ImlS5dabNPcnXfq1CmL5Vu2bIEkSfJv/YmJifjmm29w+vRpeV/WjttYvnw5+vXrB3d3d8TGxmLr1q3yc5s3b4YkSVizZk2j161YsQKSJGHXrl3NbttgMCAzMxN9+vSBu7s7AgICMHLkSGzcuFFe5/p++enTp1t0YV77c+04BWuPA2uUl5dDqVRaHHslJSVQKBSNjtPZs2cjNDTUor3mWp86dQpBQUEAgMzMzCbbDQDnz5/HpEmT4O3tjaCgIDzzzDMwGo2ttjM7OxtJSUkIDAyUj5mZM2cCAIQQiIyMxH333dfodZcvX4avry9+97vfAbh6/Hz55ZfIzMxEt27d4OPjgwceeAAVFRXQ6/V4+umnERwcDG9vb8yYMaPZz9eqVaswYMAAeHh4ICEhAT///DMA4N1330Xv3r3h7u6OxMTERscwAOzZswd33XUXfH194enpiTFjxmDHjh3y8wsWLMCzzz4LAOjZs6dcT/O2WvqMN1f3//u//0NYWBg0Gg169uyJ2bNno66urtXaX+/AgQNITk6GVquFt7c37rzzTuzevbvRej/99BPGjBkDDw8PdO/eHX/961/x8ccfN/m5vl5BQQFmzJiB7t27Q6PRoGvXrrjvvvsave7bb7/FmDFj4OPjA61Wi+HDh2PFihXy89u2bcOvf/1rREREyJ+VP/7xj1afCf3ss88QGxsLDw8P+Pv746GHHsLZs2eteq2z6ry/elGTTpw4AQAICAgA0HCW5pNPPsEDDzyAuXPnYs+ePVi4cCGOHDkifyGOGjUKb7zxBg4dOoRBgwYBaPiwKRQKbNu2DX/4wx/kZQAwevRoAA2/cd17773Yvn07Hn/8cfTv3x8///wz3njjDRw9ehRr1661aNumTZvw5ZdfIjU1FYGBgS1+uS9YsACZmZkYN24cZs+ejby8PLzzzjv48ccfsWPHjjb95t2/f3+88MILmDdvHh5//HGMGjUKADBixAh5nUuXLuHuu+/Ggw8+iIcffhhffvklZs+eDbVaLX85WevPf/4zKioqcO7cObzxxhsAYNWg0h9++AErV67EH/7wB2g0Grz99tu46667sHfvXgwaNAiJiYkIDw/H8uXLMXnyZIvXLl++HFFRUUhISGh2+wsWLMDChQvx2GOPIS4uDpWVlcjOzsb+/fsxfvz4Jl/zu9/9DuPGjbNYtn79eixfvhzBwcEA2n4ctMbPzw+DBg3C1q1b5WNv+/btkCQJZWVlOHz4MAYOHAig4Zg0/3teLygoCO+88w5mz56NyZMnY8qUKQCAIUOGyOsYjUYkJSUhPj4e//jHP/D999/jtddeQ1RUFGbPnt1sG4uKijBhwgQEBQUhPT0dfn5+OHXqFP7zn/8AaPji/s1vfoNXXnkFZWVl8Pf3l1/71VdfobKyEr/5zW8strlw4UJ4eHggPT0dx48fx1tvvQWVSgWFQoFLly5hwYIF2L17N5YtW4aePXti3rx5Fq/ftm0b/vvf/2LOnDny9n71q1/hueeew9tvv43f//73uHTpEl555RXMnDkTmzZtkl+7adMmJCcnIzY2FvPnz4dCocDHH3+MO+64A9u2bUNcXBymTJmCo0eP4vPPP8cbb7yBwMBAuc7Xbseaz/iFCxcQFxeH8vJyPP7444iOjsb58+exevVq1NTUQK1WN1v76x06dAijRo2CVqvFc889B5VKhXfffReJiYn44YcfEB8fD6AhOI0dOxaSJCEjIwNeXl744IMPoNForNrP/fffj0OHDuHJJ59EZGQkioqKsHHjRpw5c0Z+n8uWLcPMmTMxcOBAZGRkwM/PDwcOHMD69esxbdo0AMCqVatQU1OD2bNnIyAgAHv37sVbb72Fc+fOYdWqVS224aWXXsLzzz+PBx98EI899hiKi4vx1ltvYfTo0Thw4AD8/PysrptTEdQpffzxxwKA+P7770VxcbE4e/as+OKLL0RAQIDw8PAQ586dEzk5OQKAeOyxxyxe+8wzzwgAYtOmTUIIIYqKigQA8fbbbwshhCgvLxcKhUL8+te/FiEhIfLr/vCHPwh/f39hMpmEEEJ8+umnQqFQiG3btllsf+nSpQKA2LFjh7wMgFAoFOLQoUOtvreioiKhVqvFhAkThNFolJcvXrxYABAfffSRvKxHjx4iJSWl0TbGjBkjxowZIz/+8ccfBQDx8ccfN7kuAPHaa6/Jy/R6vRg6dKgIDg4WdXV1QoirNc/Pz7d4/ebNmwUAsXnzZnnZPffcI3r06NHqezUDIACI7Oxsednp06eFu7u7mDx5srwsIyNDaDQaUV5eLi8rKioSbm5uYv78+S3uIyYmRtxzzz0trjN//nzR0n8Lx44dE76+vmL8+PGivr5eCNG248Bac+bMsTj20tLSxOjRo0VwcLB45513hBBClJaWCkmSxD//+U95vZSUFIu6FxcXCwBN1iYlJUUAEC+88ILF8ltuuUXExsa22L41a9YIAOLHH39sdp28vDwBQG6v2b333isiIyPlz5H5+Bk0aJB8rAkhxMMPPywkSRLJyckWr09ISGh0bAEQGo3G4th89913BQARGhoqKisr5eUZGRkWx7HJZBJ9+vQRSUlJcpuEEKKmpkb07NlTjB8/Xl726quvNvkZMLehuc/49f8Gjz76qFAoFE3W79o2NOX6bU2aNEmo1Wpx4sQJedmFCxeEj4+PGD16tLzsySefFJIkiQMHDsjLSktLhb+/f7PvyezSpUsCgHj11VebXae8vFz4+PiI+Ph4UVtb2+x7qqmpafTahQsXCkmSxOnTp+Vl138WT506JZRKpXjppZcsXvvzzz8LNze3RstdCbuWOrlx48YhKCgI4eHheOihh+Dt7Y01a9agW7du+N///gegoevoWnPnzgUAfPPNNwAafqOKjo6WuzF27NgBpVKJZ599FoWFhTh27BiAht/4Ro4cKXc9rFq1Cv3790d0dDRKSkrkH3P31ubNmy32O2bMGAwYMKDV9/T999+jrq4OTz/9NBSKq4forFmzoNVq5Xa3Jzc3N/lUPwCo1Wr87ne/Q1FREfbt29fu+2tKQkICYmNj5ccRERG477778N1338ldHY8++ij0ej1Wr14tr7dy5UrU19c3+g3/en5+fjh06JD879lWOp0OkydPRpcuXfD5559DqVQCaPtxYI1Ro0ahsLBQHri7bds2jB49GqNGjZLPDG7fvh1CiGbPyFjriSeeaLTvkydPtvga82++X3/9NQwGQ5Pr9O3bF/Hx8Vi+fLm8rKysDN9++y0eeeSRRlNrH330UYszjfHx8RBCNDojGB8fj7Nnz6K+vt5i+Z133mlxBsR8JuL++++Hj49Po+Xm95iTk4Njx45h2rRpKC0tlf/9dDod7rzzTmzduhUmk6nFephZ8xk3mUxYu3YtJk6cKI/vu1ZbphwbjUZs2LABkyZNQq9eveTlXbt2xbRp07B9+3ZUVlYCaDiTmJCQgKFDh8rr+fv745FHHml1Px4eHlCr1diyZQsuXbrU5DobN25EVVUV0tPTG40JuvY9eXh4yH/X6XQoKSnBiBEjIITAgQMHmm3Df/7zH5hMJjz44IMWn7PQ0FD06dPnhj5nzoJBppNbsmQJNm7ciM2bN+Pw4cM4efIkkpKSAACnT5+GQqFA7969LV4TGhoKPz8/nD59Wl527RfEtm3bMGzYMAwbNgz+/v7Ytm0bKisrcfDgQYsvjWPHjuHQoUMICgqy+Onbty+AhtPv1+rZs6dV78ncrn79+lksV6vV6NWrl0W720tYWBi8vLwslpnfR2t95+2lT58+jZb17dsXNTU1KC4uBgBER0dj+PDhFl+Oy5cvx2233dbo3/l6L7zwAsrLy9G3b18MHjwYzz77LH766Ser2zdr1iycOHECa9askbsugbYfB9YwH2fbtm2DTqfDgQMHMGrUKIwePdriONVqtYiJiWnz9s3c3d0tukYAoEuXLs1+WZmNGTMG999/PzIzMxEYGIj77rsPH3/8caOxK48++ih27NghH7OrVq2CwWDAb3/720bbjIiIsHjs6+sLAAgPD2+03GQyoaKi4oZfD0B+j+Zgm5KS0ujf8IMPPoBer2+0r+ZY8xkvLi5GZWWl3I19M4qLi1FTU9Po/wqgoUvZZDLJ40dOnz7d5Gektc8NAGg0Grz88sv49ttvERISgtGjR+OVV15BQUGBvI65W7+193XmzBlMnz4d/v7+8risMWPGAECLdT527BiEEOjTp0+jf6cjR47c0OfMWXCMTCcXFxfX5G8117LmN5yRI0fi/fffx8mTJ+VxB5IkYeTIkdi2bRvCwsJgMpksgozJZMLgwYPx+uuvN7nN6/8DvfY3kfbS3HszGo3yGYOO2FdHevTRR/HUU0/h3Llz0Ov12L17NxYvXtzq60aPHo0TJ05g3bp12LBhAz744AO88cYbWLp0KR577LEWX/vPf/4Tn3/+OT777DOL32iBth8H1ggLC0PPnj2xdetWREZGQgiBhIQEBAUF4amnnsLp06exbds2jBgxwuKsXVvd6DEiSRJWr16N3bt346uvvsJ3332HmTNn4rXXXsPu3bvlMVEPPfQQ/vjHP2L58uX405/+hM8++wzDhg1r8ou3ubY0t1xcM+j5Zl5vPtvy6quvNvq3NbP2woG2+Iw7iqeffhoTJ07E2rVr8d133+H555/HwoULsWnTJtxyyy1WbcNoNGL8+PEoKyvD//t//w/R0dHw8vLC+fPnMX369BbPfJlMJkiShG+//bbJf1NXvrgjg4wT69GjB0wmE44dO4b+/fvLywsLC1FeXo4ePXrIy8wBZePGjfjxxx+Rnp4OoOHL75133pHPWFzb9REVFYWDBw/izjvvbNcrUJrblZeXZ3G6uK6uDvn5+RaDT7t06dLkNS1Onz5t8drW2nfhwgXodDqLszJHjx4FAPl0fZcuXQCg0f6aOkN0I/Voqsvn6NGj8PT0tDhr8NBDDyEtLQ2ff/45amtroVKpMHXqVKv24e/vjxkzZmDGjBmorq7G6NGjsWDBghaDzLZt2/DMM8/g6aefbvI0vK2Og1GjRmHr1q3o2bMnhg4dCh8fH8TExMDX1xfr16/H/v37kZmZ2eI2bH1l1Ntuuw233XYbXnrpJaxYsQKPPPIIvvjiC7me/v7+uOeee7B8+XI88sgj2LFjh8NdXygqKgoAoNVqGw3svl571DMoKAharRa//PJLu2zL09PT4tpBZrm5uVAoFHKQ7tGjB44fP95ovaaWNScqKgpz587F3LlzcezYMQwdOhSvvfYaPvvsM7mOv/zyS7NneX7++WccPXoUn3zyCR599FF5+bUzB1vatxACPXv2lM92UgN2LTmxu+++GwAa/cdp/s35nnvukZf17NkT3bp1wxtvvAGDwYDbb78dQMOXyYkTJ7B69WrcdtttFtcYefDBB3H+/Hm8//77jfZdW1sLnU53Q+0eN24c1Go13nzzTYvfOj/88ENUVFRYtDsqKgq7d++2mLL59ddfN5qOaA4ozV3Iq76+Hu+++678uK6uDu+++y6CgoLk8Gb+j+raKdFGoxHvvfdeo+15eXlZfTrebNeuXdi/f7/8+OzZs1i3bh0mTJhg8RtYYGAgkpOT8dlnn2H58uW466675BkkLSktLbV47O3tjd69e7c4RfrixYt48MEHMXLkSLz66qtNrmOr42DUqFE4deoUVq5cKQdthUKBESNG4PXXX4fBYGh1fIynpyeA5v/db9SlS5canRExn824vp6//e1vcfjwYTz77LNQKpV46KGH2rUtNys2NhZRUVH4xz/+gerq6kbPm7s1gdY/R9ZQKBSYNGkSvvrqK2RnZzd6/vq6tkSpVGLChAlYt26dRRdwYWEhVqxYgZEjR0Kr1QIAkpKSsGvXLouLYpaVlVl00zanpqYGly9ftlgWFRUFHx8f+d97woQJ8PHxwcKFCxuta35P5s/xte9RCIF//vOfrbZhypQpUCqVyMzMbFQjIUSjz7cr4RkZJxYTE4OUlBS89957KC8vx5gxY7B371588sknmDRpEsaOHWux/qhRo/DFF19g8ODB8tmHW2+9FV5eXjh69Kg8fdDst7/9Lb788ks88cQT2Lx5M26//XYYjUbk5ubiyy+/xHfffddqt1dTgoKCkJGRgczMTNx111249957kZeXh7fffhvDhw+3GNT62GOPYfXq1bjrrrvw4IMP4sSJExa/HZlFRUXBz88PS5cuhY+PD7y8vBAfHy/36YeFheHll1/GqVOn0LdvX6xcuRI5OTl477335AGYAwcOxG233YaMjAx5Su0XX3zRaNAl0PDlsHLlSqSlpWH48OHw9vbGxIkTW3zfgwYNQlJSksX0awBNnnV49NFH8cADDwAAXnzxRavqOmDAACQmJiI2Nhb+/v7Izs7G6tWrW7zH0x/+8AcUFxfjueeewxdffGHx3JAhQzBkyJA2HQfmafWbN2+Wr/PTHHNIycvLw9/+9jd5+ejRo/Htt99Co9HI1wdqjoeHBwYMGICVK1eib9++8Pf3x6BBg256fMYnn3yCt99+G5MnT0ZUVBSqqqrw/vvvQ6vVyr9AmN1zzz0ICAjAqlWrkJycLE9bdxQKhQIffPABkpOTMXDgQMyYMQPdunXD+fPnsXnzZmi1Wnz11VcAIIf6P//5z3jooYegUqkwceLERuPLWvO3v/0NGzZswJgxY+Qp+xcvXsSqVauwffv2Nk0j/utf/4qNGzdi5MiR+P3vfw83Nze8++670Ov1eOWVV+T1nnvuOXz22WcYP348nnzySXn6dUREBMrKylo823T06FHceeedePDBBzFgwAC4ublhzZo1KCwslIOpVqvFG2+8gcceewzDhw/HtGnT0KVLFxw8eBA1NTX45JNPEB0djaioKDzzzDM4f/48tFot/v3vf7c6Jgto+D/sr3/9KzIyMnDq1ClMmjQJPj4+yM/Px5o1a/D444/jmWeesbpuTqXjJ0pRezBPBW5p+qcQQhgMBpGZmSl69uwpVCqVCA8PFxkZGeLy5cuN1l2yZIkAIGbPnm2xfNy4cQKAyMrKavSauro68fLLL4uBAwcKjUYjunTpImJjY0VmZqaoqKiQ1wMg5syZ06b3uHjxYhEdHS1UKpUICQkRs2fPFpcuXWq03muvvSa6desmNBqNuP3220V2dnaj6ddCCLFu3ToxYMAA4ebmZjEVe8yYMWLgwIEiOztbJCQkCHd3d9GjRw+xePHiRvs6ceKEGDdunNBoNCIkJET86U9/Ehs3bmw0/bq6ulpMmzZN+Pn5CQCtTsU21+ezzz4Tffr0ERqNRtxyyy0W27yWXq8XXbp0Eb6+vo2mejbnr3/9q4iLixN+fn7Cw8NDREdHi5deesliyu/1Uz7NU9Ob+rl2Cqy1x8HcuXOFJEniyJEjVrU5ODhYABCFhYXysu3btwsAYtSoUY3Wv376tRBC7Ny5U8TGxgq1Wm3R7pSUFOHl5dVoG61NQRdCiP3794uHH35YRERECI1GI4KDg8WvfvUri+nz1/r9738vAIgVK1Y0es48/XrVqlUWy5v7jJvbV1xcLC9r6vOVn5/f5JTh5vZ34MABMWXKFBEQECA0Go3o0aOHePDBBxt97l988UXRrVs3oVAoLKYtt/QZv/54EaLh8gKPPvqoCAoKEhqNRvTq1UvMmTNH6PX6JrfR0rb2798vkpKShLe3t/D09BRjx44VO3fubPTaAwcOiFGjRgmNRiO6d+8uFi5cKN58800BQBQUFDS7z5KSEjFnzhwRHR0tvLy8hK+vr4iPjxdffvllo3X/+9//ihEjRggPDw+h1WpFXFyc+Pzzz+XnDx8+LMaNGye8vb1FYGCgmDVrljh48GCjy0M0dxz++9//FiNHjhReXl7Cy8tLREdHizlz5oi8vLwW6+bMJCHacB6PyAklJiaipKSkXfrsO0p9fT3CwsIwceJEfPjhh/ZujtXi4uLQo0ePVi/85Wz++Mc/4sMPP0RBQYHc3UWO4emnn8a7776L6urqdp8gQB2DXUtEndDatWtRXFxsMWDQ0Zmn8H/yySf2bkqHunz5Mj777DPcf//9DDF2VltbazGzqrS0FJ9++ilGjhzJENOJMcgQdSJ79uzBTz/9hBdffBG33HKLfP2JzkCr1d7QvZc6q6KiInz//fdYvXo1SktL8dRTT9m7SS4vISEBiYmJ6N+/PwoLC/Hhhx+isrISzz//vL2bRjeBQYaoE3nnnXfka7mYb3pJjunw4cN45JFHEBwcjDfffLPZa7RQx7n77ruxevVqvPfee5AkCbfeeis+/PBD+f5x1DlxjAwRERF1WryODBEREXVaDDJERETUaTn9GBmTyYQLFy7Ax8fH5pcrJyIiovYhhEBVVRXCwsJavKea0weZCxcu3NBN64iIiMj+zp49i+7duzf7vNMHGR8fHwANhTDfc+NmGQwGbNiwARMmTJAvX09NY63ahvWyHmtlPdaqbVgv69myVpWVlQgPD5e/x5vj9EHG3J2k1WrbNch4enpCq9XyIG8Fa9U2rJf1WCvrsVZtw3pZryNq1dqwEA72JSIiok6LQYaIiIg6LQYZIiIi6rQYZIiIiKjTYpAhIiKiTotBhoiIiDotBhkiIiLqtBhkiIiIqNNikCEiIqJOi0GGiIiIOi0GGSIiIuq0HD7InD17FomJiRgwYACGDBmCVatW2btJRERE5CAcPsi4ublh0aJFOHz4MDZs2ICnn34aOp3O3s1qN2dKa/D+1pMor6mzd1OIiIg6HYe/+3XXrl3RtWtXAEBoaCgCAwNRVlYGLy8vO7fs5tUbTZj5yY84XlSNf+0+hVfuj0F8T38oFC3f6ZOIiIga2PyMzNatWzFx4kSEhYVBkiSsXbu20TpLlixBZGQk3N3dER8fj7179za5rX379sFoNCI8PNzGrbatbceK8YfPD+AfG47ieFE1AOBsWS0efn83Rr68CfvPXLJzC4mIiDoHmwcZnU6HmJgYLFmypMnnV65cibS0NMyfPx/79+9HTEwMkpKSUFRUZLFeWVkZHn30Ubz33nu2brLN/f3bXPz34AUs/eEEAODJO3rjgdju8Na44ULFZcxc9iOOF1XZuZVERESOz+ZdS8nJyUhOTm72+ddffx2zZs3CjBkzAABLly7FN998g48++gjp6ekAAL1ej0mTJiE9PR0jRoxocX96vR56vV5+XFlZCQAwGAwwGAw3+3bkbV37Z1uUVOtx6EKl/Lh7Fw/8bmQPaFRK/CW5L1KW7cPBcxWY8fGP+GpOArw0Dt/716KbqZUrYr2sx1pZj7VqG9bLeraslbXblIQQot333tzOJAlr1qzBpEmTAAB1dXXw9PTE6tWr5WUAkJKSgvLycqxbtw5CCEybNg39+vXDggULWt3HggULkJmZ2Wj5ihUr4Onp2U7vpO1+uSRhX7GE7l4C/z2jRHcvgYejjNCqAK366no6A/DqT0pcqpMwOtSE+3ua7NZmIiIie6mpqcG0adNQUVEBrVbb7Hp2/XW/pKQERqMRISEhFstDQkKQm5sLANixYwdWrlyJIUOGyONrPv30UwwePLjJbWZkZCAtLU1+XFlZifDwcEyYMKHFQrSFwWDAxo0bMX78eKhUqlbX19eb8OJrW1FSXYcDZQ3LfhXbC4+P79Pk+qEDSzDzk/3YWqDAE/fEY3hkl3Zptz20tVaujvWyHmtlPdaqbVgv69myVuYeldY4fL/FyJEjYTJZf1ZCo9FAo9E0Wq5Sqdq9yNZuc03OWZRUN0yvNp//SowOafa1d/TviqnDwrEy+yz+vO4w/veHURAQ8FQ7/D9Xs2xRf2fGelmPtbIea9U2rJf1bPUdaw27XkcmMDAQSqUShYWFFssLCwsRGhpqp1a1LyEE3t92EgAQ4d/QteWlVuLWiJbPsvz5V/0RqnVHfokOcS99j8ELNuCz3adt3l4iIqLOxK5BRq1WIzY2FllZWfIyk8mErKwsJCQk2LFl7Wf1vnM4VlQNb40bVj2RgIfjwjH/3oFQu7Vceq27CgunNHSfVenrYTQJPL/uF6z88Qw6cFgTERGRQ7N5X0V1dTWOHz8uP87Pz0dOTg78/f0RERGBtLQ0pKSkYNiwYYiLi8OiRYug0+nkWUyd2fnyWrzw1WEAwO/HRiFE646FU4ZY/fqx0cH47P/iUWc0YlNuET7bfQb/798/49Pdp/FhynCEaN1t1XQiIqJOweZBJjs7G2PHjpUfmwfipqSkYNmyZZg6dSqKi4sxb948FBQUYOjQoVi/fn2jAcCdjRACf/rPz6jS1+PWCD/8bnTUDW1nZJ9AAEBi32D4e2nwwbaT+OV8JZZsPo4X7hvUnk0mIiLqdGweZBITE1vtCklNTUVqaqqtm9KhvjtUiB+OFkOtVOAfv46B8iZvO6BQSEgb3xe3Rvhh+sc/4pufLuL5Xw2ASunwt8siIiKyGX4L2kBNXT1e/LqhS+nx0b3QK8i73bY9sncgArzUKNXVYcfxknbbLhERUWfEINPOhBBI//fPOF9ei25+Hpgztne7bt9NqcCvhjTcRPO/ORfaddtERESdDYNMO/t092n89+AFuCkkLHpoKDzUynbfx71DuwEAvjtUgLp6XvmXiIhcF4NMO/tgWz4AID05GsMj/W2yj1vC/eDvpYauzoiD58ptsg8iIqLOgEGmHZ0u1eFMWQ1USgkPx0XYbD8KhYQRUQEAgO3HOE6GiIhcF4NMO9p6tBgAENuji83vWn1774Zp2TtPMMgQEZHrYpBpR1uvnB0Z1SfI5vsaeSXIHDhTjmp9vc33R0RE5IgYZNqJwWjCrhOlAIDRHRBkwv09EeHviXqTwN78Upvvj4iIyBExyLST/acvoVpfjy6eKgwM03bIPs3dS1uPsnuJiIhcE4NMO9lwuOEO3mP7BUNxk1fxtVZiv4YzP5tyi3gjSSIickkMMu1ACIH1vxQAAJIGhXbYfkf2DoRKKeFMWQ1Olug6bL9ERESOgkGmHRy6UInz5bVwVyk6ZHyMmZfGDfE9G6Zhb84t6rD9EhEROQoGmXbw3aGGszGJfYNtciXfloyNDgYAbMkr7tD9EhEROQIGmXaw7cq06/EDQjp832P6Ngz43XuqjONkiIjI5TDI3CQhBE4WVwMABnXz7fD9h/l5AADq6k2oqTN2+P6JiIjsiUHmJl2qMaDycsMF6XoEeHb4/j1USqiUDbOkKmoNHb5/IiIie2KQuUn5V2YLhfm6w13VseNjAECSJPh6qAAwyBARkethkLlJp64EmchAL7u1QXslyJTXMMgQEZFrYZC5SadK7R9keEaGiIhcFYPMTTJ3LfUMsF+Q8bsSZCoZZIiIyMUwyNwknpEhIiKyHwaZmyCEwKmSGgBAz8COn7FkxiBDRESuikHmJpRU16FaXw+FBIT7M8gQERF1NAaZm3C6rOFsTJifBzRuHT/12kzLIENERC6KQeYmlFTXAQBCtO52bQfPyBARkatikLkJVVeu6GsOEvbCIENERK6KQeYmVF1uCA5adze7tsOX06+JiMhFdYogM3nyZHTp0gUPPPCAvZtiwXyPJa29z8h48owMERG5pk4RZJ566in861//snczGpGDjLtjdC2V1xoghLBrW4iIiDpSpwgyiYmJ8PHxsXczGqm6cgZE6+EYXUtGk4CuzmjXthAREXUkmweZrVu3YuLEiQgLC4MkSVi7dm2jdZYsWYLIyEi4u7sjPj4ee/futXWz2oWjnJHxUCmhVjb8U7J7iYiIXInNTyXodDrExMRg5syZmDJlSqPnV65cibS0NCxduhTx8fFYtGgRkpKSkJeXh+Dg4DbvT6/XQ6/Xy48rKysBAAaDAQZD+3zJm7dTUdsw/dpTJbXbtm+U1sMNJdV1KK2sRbCXfc8QXctcF3vXp7NgvazHWlmPtWob1st6tqyVtduURAcOqpAkCWvWrMGkSZPkZfHx8Rg+fDgWL14MADCZTAgPD8eTTz6J9PR0eb0tW7Zg8eLFWL16dYv7WLBgATIzMxstX7FiBTw92/fqu38/qMTFGgmz+xsR7WffsSl/y1GisFZC6gAj+vhynAwREXVuNTU1mDZtGioqKqDVaptdz66/utfV1WHfvn3IyMiQlykUCowbNw67du26oW1mZGQgLS1NflxZWYnw8HBMmDChxUK0hcFgwMaNGyHc3AHoMW70CAzp7tsu275Ry87tQeHZCvSPuRUTBoTYtS3XMtdq/PjxUKns2wXXGbBe1mOtrMdatQ3rZT1b1srco9IauwaZkpISGI1GhIRYfvGGhIQgNzdXfjxu3DgcPHgQOp0O3bt3x6pVq5CQkNDkNjUaDTQaTaPlKpWq3YtcpW8YI+Pv42H3g93PUw0A0NUJu7elKbaovzNjvazHWlmPtWob1st6tqiVtdtznMEULfj+++/t3YRGTALQ6RtmCNn7gngAr+5LRESuya7TrwMDA6FUKlFYWGixvLCwEKGhoXZqlXUuXzPL2cfOs5YAwEPdEKZqDZx+TURErsOuQUatViM2NhZZWVnyMpPJhKysrGa7jhxFbUOvUsPUZzf7X47nyuxrmHhBPCIiciE27xOprq7G8ePH5cf5+fnIycmBv78/IiIikJaWhpSUFAwbNgxxcXFYtGgRdDodZsyYYeum3ZTaKyc+7H0xPDOFJAEATCYGGSIich02/xbOzs7G2LFj5cfmGUUpKSlYtmwZpk6diuLiYsybNw8FBQUYOnQo1q9f32gAsKOprW8IDva+GJ6ZHGSYY4iIyIXYPMgkJia2ev+f1NRUpKam2rop7armSteSvW8YaWYOMkZ2LRERkQux/+COTsrcteTjADOWAEDRkGM4RoaIiFwKg8wNMg/2dZSuJeWVJMMcQ0REroRB5gbVGq+MkXGQwb6SuWuJg2SIiMiFMMjcoMsOd0am4U92LRERkSthkLlBV6dfO0aQ4fRrIiJyRQwyN8jRxshInH5NREQuiEHmBjnaGBklp18TEZELYpC5QY52RsY8/bq1a/YQERE5E8c4ndAJPdrHiIGxt2FQd197NwUAoFCYx8jYuSFEREQdiEHmBoV6AvE9/aFSOcoZGXYtERGR62HXkpPg9GsiInJFDDJOgtOviYjIFTHIOAlOvyYiIlfEIOMklFdmLXGMDBERuRIGGSehkG8aySBDRESug0HGSVwdI2PnhhAREXUgBhknwenXRETkihhknIR5+jW7loiIyJUwyDgJ86wlI6ctERGRC2GQcRIKTr8mIiIXxCDjJHhlXyIickUMMk7i6hkZBhkiInIdDDJOgtOviYjIFTHIOAlOvyYiIlfEIOMkrlzYl9OviYjIpTDIOAnzLQo4/ZqIiFwJg4yT4PRrIiJyRQwyToJX9iUiIlfEIOMkJA72JSIiF8Qg4ySUnH5NREQuiEHGSfCCeERE5IoYZJyEefo1gwwREbkSBhknwenXRETkihhknIS5a4knZIiIyJUwyDgJ3v2aiIhcEYOMk+D0ayIickWdIsh8/fXX6NevH/r06YMPPvjA3s1xSJx+TURErsjN3g1oTX19PdLS0rB582b4+voiNjYWkydPRkBAgL2b5lA4/ZqIiFyRw5+R2bt3LwYOHIhu3brB29sbycnJ2LBhg72b5XAkTr8mIiIXZPMgs3XrVkycOBFhYWGQJAlr165ttM6SJUsQGRkJd3d3xMfHY+/evfJzFy5cQLdu3eTH3bp1w/nz523d7E5HKU+/tnNDiIiIOpDNg4xOp0NMTAyWLFnS5PMrV65EWloa5s+fj/379yMmJgZJSUkoKiqyddOcytXp1zwjQ0RErsPmY2SSk5ORnJzc7POvv/46Zs2ahRkzZgAAli5dim+++QYfffQR0tPTERYWZnEG5vz584iLi2t2e3q9Hnq9Xn5cWVkJADAYDDAYDDf7duRtXfunIzAZ6xv+FMKh2uWItXJkrJf1WCvrsVZtw3pZz5a1snabkujAX+ElScKaNWswadIkAEBdXR08PT2xevVqeRkApKSkoLy8HOvWrUN9fT369++PLVu2yIN9d+7c2exg3wULFiAzM7PR8hUrVsDT09MWb8shFNYCf8txg4dS4O9xRns3h4iI6KbU1NRg2rRpqKiogFarbXY9u85aKikpgdFoREhIiMXykJAQ5ObmAgDc3Nzw2muvYezYsTCZTHjuuedanLGUkZGBtLQ0+XFlZSXCw8MxYcKEFgvRFgaDARs3bsT48eOhUqnaZZs361SpDn/L2QGlmwp3351k7+bIHLFWjoz1sh5rZT3Wqm1YL+vZslbmHpXWOPz0awC49957ce+991q1rkajgUajabRcpVK1e5Ftsc0bpVGpATR0LTlKm67lSLXqDFgv67FW1mOt2ob1sp6tvmOtYdfp14GBgVAqlSgsLLRYXlhYiNDQUDu1qnPi9GsiInJFdg0yarUasbGxyMrKkpeZTCZkZWUhISHBji3rfMzTr3nzayIiciU271qqrq7G8ePH5cf5+fnIycmBv78/IiIikJaWhpSUFAwbNgxxcXFYtGgRdDqdPIuJrCNf2ZdJhoiIXIjNg0x2djbGjh0rPzYPxE1JScGyZcswdepUFBcXY968eSgoKMDQoUOxfv36RgOAqWUK3v2aiIhckM2DTGJiYqsXaUtNTUVqaqqtm+LUrt5rqeGieOa7YRMRETkzh7/XEllHcU1w4UkZIiJyFQwyTkJ5TZAxMskQEZGLYJBxEtI1/5IcJ0NERK6CQcZJKNm1RERELohBxklcO0bGyCnYRETkIhhknISCXUtEROSCGGScxLVnZEwmOzaEiIioAzHIOAmLIMMzMkRE5CIYZJyE4prr33H6NRERuQoGGSchSRLvgE1ERC7H4YPMkiVLEBkZCXd3d8THx2Pv3r32bpLDMk/BZo4hIiJX4dBBZuXKlUhLS8P8+fOxf/9+xMTEICkpCUVFRfZumkMyj5Ph9GsiInIVDh1kXn/9dcyaNQszZszAgAEDsHTpUnh6euKjjz6yd9McEruWiIjI1dj87tc3qq6uDvv27UNGRoa8TKFQYNy4cdi1a1ezr9Pr9dDr9fLjyspKAIDBYIDBYGiXtpm3017bay/KKyN+9XXt915vlqPWylGxXtZjrazHWrUN62U9W9bK2m06bJApKSmB0WhESEiIxfKQkBDk5uY2+7qFCxciMzOz0fINGzbA09OzXdu4cePGdt3ezTIalQAkbN68BUEe9m6NJUerlaNjvazHWlmPtWob1st6tqhVTU2NVes5bJC5URkZGUhLS5MfV1ZWIjw8HBMmTIBWq22XfRgMBmzcuBHjx4+HSqVql222h+cPbILeWI9Ro8egV5CXvZsDwHFr5ahYL+uxVtZjrdqG9bKeLWtl7lFpjcMGmcDAQCiVShQWFlosLywsRGhoaLOv02g00Gg0jZarVKp2L7IttnkzFFe6lpRuSodqF+B4tXJ0rJf1WCvrsVZtw3pZz1bfsdZw2MG+arUasbGxyMrKkpeZTCZkZWUhISHBji1zXObp15y0RERErsJhz8gAQFpaGlJSUjBs2DDExcVh0aJF0Ol0mDFjhr2b5pAkTr8mIiIX49BBZurUqSguLsa8efNQUFCAoUOHYv369Y0GAFMDBadfExGRi3HoIAMAqampSE1NtXczOgXz9Gve/ZqIiFyFw46RobZTyGNkeEaGiIhcA4OME1Fc+ddkkCEiIlfBIONEeEaGiIhcDYOME+H0ayIicjUMMk7EfNNITr8mIiJXwSDjRNi1REREroZBxolw+jUREbkaBhknIvGMDBERuRgGGSei5PRrIiJyMQwyToRjZIiIyNV0iiDz9ddfo1+/fujTpw8++OADezfHYcldSxwjQ0RELsLh77VUX1+PtLQ0bN68Gb6+voiNjcXkyZMREBBg76Y5HKV5+jXPyBARkYtw+DMye/fuxcCBA9GtWzd4e3sjOTkZGzZssHezHJK5a0kwyBARkYuwaZDZunUrJk6ciLCwMEiShLVr1za53pIlSxAZGQl3d3fEx8dj79698nMXLlxAt27d5MfdunXD+fPnbdnsTktxZfq1kV1LRETkImwaZHQ6HWJiYrBkyZJm11m5ciXS0tIwf/587N+/HzExMUhKSkJRUZEtm+aUruQYDvYlIiKXYdMxMsnJyUhOTm5xnddffx2zZs3CjBkzAABLly7FN998g48++gjp6ekICwuzOANz/vx5xMXFNbs9vV4PvV4vP66srAQAGAwGGAyGm3k7MvN22mt77cWcSg319Q7TNketlaNivazHWlmPtWob1st6tqyVtduURAcNqJAkCWvWrMGkSZPkZXV1dfD09MTq1astlqekpKC8vBzr1q1DfX09+vfvjy1btsiDfXfu3NnsYN8FCxYgMzOz0fIVK1bA09Ozvd+WQ3n7sAJ5FQr8prcRw4N4VoaIiDqvmpoaTJs2DRUVFdBqtc2uZ9dZSyUlJTAajQgJCbFYHhISgtzcXACAm5sbXnvtNYwdOxYmkwnPPfdcizOWMjIykJaWJj+urKxEeHg4JkyY0GIh2sJgMGDjxo0YP348VCpVu2yzPawu3oe8ilIMGRKDu28Js3dzADhurRwV62U91sp6rFXbsF7Ws2WtzD0qrWlzkElPT8fLL7/c4jpHjhxBdHR0WzfdrHvvvRf33nuvVetqNBpoNJpGy1UqVbsX2RbbvBlu5kv7KhQO1S7A8Wrl6Fgv67FW1mOt2ob1sp6tvmOt0eYgM3fuXEyfPr3FdXr16mXVtgIDA6FUKlFYWGixvLCwEKGhoW1tmsvj9GsiInI1bQ4yQUFBCAoKapedq9VqxMbGIisrSx4jYzKZkJWVhdTU1HbZhyvh9GsiInI1Nh0jU11djePHj8uP8/PzkZOTA39/f0RERAAA0tLSkJKSgmHDhiEuLg6LFi2CTqeTZzGR9Tj9moiIXI1Ng0x2djbGjh0rPzYPwk1JScGyZcsAAFOnTkVxcTHmzZuHgoICDB06FOvXr280AJhap1Swa4mIiFyLTYNMYmKiVV+qqamp7EpqB+abRhpNDDJEROQaHP5eS2Q982Bf5hgiInIVDDJORMkxMkRE5GIYZJzI1TMyDDJEROQaGGSciHn6NbuWiIjIVTDIOBHz9GsO9iUiIlfBIONEeGVfIiJyNQwyToRX9iUiIlfDIONEeGVfIiJyNQwyTkTJWUtERORiGGSciMQgQ0RELoZBxokoOf2aiIhcDIOME5HHyDDJEBGRi2CQcSK8si8REbkaBhknwunXRETkahhknAinXxMRkathkHEinH5NRESuhkHGiXD6NRERuRoGGSfC6ddERORqGGScCKdfExGRq2GQcSLsWiIiIlfDIONElJx+TURELoZBxomYu5YEz8gQEZGLYJBxIuYr+xoZZIiIyEUwyDiRq7cosHNDiIiIOgiDjBPhlX2JiMjVMMg4Efk6MjwlQ0RELoJBxolw+jUREbkaBhknwunXRETkahhknAinXxMRkathkHEiCnYtERGRi2GQcSJXryNj54YQERF1EAYZJ6K48q/JriUiInIVDDJORD4jw+nXRETkIhhknAjHyBARkathkHEiVy+IZ+eGEBERdRAGGSfCWxQQEZGrYZBxIuxaIiIiV9MpgszkyZPRpUsXPPDAA/ZuikPj9GsiInI1nSLIPPXUU/jXv/5l72Y4PE6/JiIiV9MpgkxiYiJ8fHzs3QyHx+nXRETkamwaZLZu3YqJEyciLCwMkiRh7dq1Ta63ZMkSREZGwt3dHfHx8di7d68tm+W0ro6RsXNDiIiIOoibLTeu0+kQExODmTNnYsqUKU2us3LlSqSlpWHp0qWIj4/HokWLkJSUhLy8PAQHB7d5n3q9Hnq9Xn5cWVkJADAYDDAYDDf2Rq5j3k57ba+9CJMRAGA0mhymbY5aK0fFelmPtbIea9U2rJf1bFkra7cpiQ4aUCFJEtasWYNJkyZZLI+Pj8fw4cOxePFiAIDJZEJ4eDiefPJJpKeny+tt2bIFixcvxurVq1vcz4IFC5CZmdlo+YoVK+Dp6Xnzb8SBHauQsPiwEqEeAhlDjfZuDhER0Q2rqanBtGnTUFFRAa1W2+x6Nj0j05q6ujrs27cPGRkZ8jKFQoFx48Zh165dN7TNjIwMpKWlyY8rKysRHh6OCRMmtFiItjAYDNi4cSPGjx8PlUrVLttsD3vyy7D4cDY8vbxx992327s5ABy3Vo6K9bIea2U91qptWC/r2bJW5h6V1tg1yJSUlMBoNCIkJMRieUhICHJzc+XH48aNw8GDB6HT6dC9e3esWrUKCQkJTW5To9FAo9E0Wq5Sqdq9yLbY5s3QqBvaIgCHahfgeLVydKyX9Vgr67FWbcN6Wc9W37HWaPNg3/T0dEiS1OLPtSGkPXz//fcoLi5GTU0Nzp0712yIcXW8si8REbmaNp+RmTt3LqZPn97iOr169bJqW4GBgVAqlSgsLLRYXlhYiNDQ0LY2zeVx+jUREbmaNgeZoKAgBAUFtcvO1Wo1YmNjkZWVJQ8CNplMyMrKQmpqarvsw5WYgwxPyBARkauw6RiZ6upqHD9+XH6cn5+PnJwc+Pv7IyIiAgCQlpaGlJQUDBs2DHFxcVi0aBF0Oh1mzJhhy6Y5JfPdr3lGhoiIXIVNg0x2djbGjh0rPzbPJkpJScGyZcsAAFOnTkVxcTHmzZuHgoICDB06FOvXr280AJhaJ3GMDBERuRibBpnExESr7vuTmprKrqR2wCv7EhGRq+kU91oi65i7lnhGhoiIXAWDjBPh9GsiInI1DDJOhNOviYjI1TDIOBFOvyYiIlfDIONEOEaGiIhcDYOMEzFPv2bXEhERuQoGGSdybdeSySSQ/u+f8NevD9u5VURERLbDIONE5Cv7CoEDZy/hix/P4oPt+Siu0tu5ZURERLbBIONENG4N/5xGk8Dy3Wfk5UcuVtqrSURERDbFIONE/DzVGB7ZBQDwnwPn5eWHGWSIiMhJMcg4mZQRkY2W8YwMERE5KwYZJ5M0MBShWncAgI+m4VZahy8wyBARkXNikHEyKqUCfxzfB0qFhPS7owEAJ4qrcdlgtHPLiIiI2h+DjBOaOjwCx19KxrS4CAR4qWESQF5Blb2bRURE1O4YZJyUJEmQJAn9u2oBcJwMERE5JwYZJzcgrCHIcOYSERE5IwYZJzfgyhkZDvglIiJnxCDj5MxnZHILqmDiPZiIiMjJMMg4uV6BXlC7KVCtr8fZSzX2bg4REVG7YpBxcm5KBfqF+ABg9xIRETkfBhkXII+T4YBfIiJyMgwyLqB/14YzMpyCTUREzoZBxgUMCPMFwK4lIiJyPgwyLqB/Vx8oJOBCxWVcKK+1d3OIiIjaDYOMC/BxV2FouB8A4IejxfZtDBERUTtikHERif2CAQBb8ors3BIiIqL24/BB5uzZs0hMTMSAAQMwZMgQrFq1yt5N6pQS+wUBAHYcL0VdvcnOrSEiImofDh9k3NzcsGjRIhw+fBgbNmzA008/DZ1OZ+9mdTqDwnwR6K1Gtb4e2afL7N0cIiKiduHwQaZr164YOnQoACA0NBSBgYEoK+MXcVspFBJG9204K/NDHsfJEBGRc7BpkNm6dSsmTpyIsLAwSJKEtWvXNrnekiVLEBkZCXd3d8THx2Pv3r1Nrrdv3z4YjUaEh4fbsNXO6+o4GQYZIiJyDm623LhOp0NMTAxmzpyJKVOmNLnOypUrkZaWhqVLlyI+Ph6LFi1CUlIS8vLyEBwcLK9XVlaGRx99FO+//36L+9Tr9dDr9fLjysqGa6cYDAYYDIZ2eFeQt9Ne2+soCZF+UEhAXmEVzpRUoauvu8332VlrZS+sl/VYK+uxVm3DelnPlrWydpuSEKJDboksSRLWrFmDSZMmWSyPj4/H8OHDsXjxYgCAyWRCeHg4nnzySaSnpwNoCCfjx4/HrFmz8Nvf/rbF/SxYsACZmZmNlq9YsQKenp7t82Y6sTd+VuJUtYSpvYwYEcK7YRMRkWOqqanBtGnTUFFRAa1W2+x6Nj0j05q6ujrs27cPGRkZ8jKFQoFx48Zh165dAAAhBKZPn4477rij1RADABkZGUhLS5MfV1ZWIjw8HBMmTGixEG1hMBiwceNGjB8/HiqVql222VHyPU9iUdZxlGm64u67h9p8f525VvbAelmPtbIea9U2rJf1bFkrc49Ka+waZEpKSmA0GhESEmKxPCQkBLm5uQCAHTt2YOXKlRgyZIg8xubTTz/F4MGDm9ymRqOBRqNptFylUrV7kW2xTVu7s38oFmUdx84TpagXCniolR2y385YK3tivazHWlmPtWob1st6tvqOtUabB/ump6dDkqQWf8whpD2MHDkSJpMJOTk58k9zIYZaNzBMi3B/D+jqjPjq4AV7N4eIiOimtPmMzNy5czF9+vQW1+nVq5dV2woMDIRSqURhYaHF8sLCQoSGhra1aWQFhULCb+J7YOG3uVi28xR+Paw7JEmyd7OIiIhuSJuDTFBQEIKCgtpl52q1GrGxscjKypIHAZtMJmRlZSE1NbVd9kGNPTgsHK9vPIrDFyux/8wlxPbwt3eTiIiIbohNryNTXV0tdwcBQH5+PnJycnDmzBl5nbS0NLz//vv45JNPcOTIEcyePRs6nQ4zZsywZdNcWhcvNe4bGgYAWL7nTCtrExEROS6bDvbNzs7G2LFj5cfm2UQpKSlYtmwZAGDq1KkoLi7GvHnzUFBQgKFDh2L9+vWNBgBT+5o6PBxfZp/D+l8K8OJ99fDS2HXcNxER0Q2x6bdXYmIirLlMTWpqKruSOtitEV0QGeCJU6U1WP9LAe6P7W7vJhEREbWZw99riWxDkiRMubUhvPx7/zk7t4aIiOjGMMi4sMm3dAMA7DxRij0nS+3cGiIiorZjkHFh4f6eeGh4ww04n/v3T6itM9q5RURERG3DIOPi/nRPf3T1dcfp0hr8Y0OevZtDRETUJgwyLk7rrsLfpjRcKfmjHfnIPlVm5xYRERFZj0GGMLZfMB6I7Q4hgOdW/4TLBnYxERFR58AgQwCA5+8ZgBCtBidLdHhj41F7N4eIiMgqDDIEAPD1VOFvkxu6mN7fdhIHzlyyc4uIiIhaxyBDsjv7h2DyLd1gEsCLXx+26mKGRERE9sQgQxYykqPhrlJg/5lybMotsndziIiIWsQgQxaCte6YPqInAODV7/Kgr+fAXyIiclwMMtTIE2N6QevuhtyCKsz98iBMJnYxERGRY2KQoUb8PNV4+5FYqJQSvv7pIt7cdMzeTSIiImoSgww1aWSfQCycMgQA8PaWEzhdqrNzi4iIiBpjkKFm3X9rN4zqE4i6ehMyv+IsJiIicjwMMtQsSZKw4N6BUCklbMotwmd7zti7SURERBYYZKhFUUHeeDapHwDgha8OYc/JUju3iIiI6CoGGWrVrFG9kDwoFAajwG8/2ot1Oeft3SQiIiIADDJkBUmS8NqDMRjXPwR19SY89UUO/vFdHqdlExGR3THIkFU81W5477exeGJMFABg8ebj+P3y/aipq7dzy4iIyJUxyJDVFAoJ6cnR+MevY6BSSlh/qABT3t6JbceKOaOJiIjsgkGG2uyB2O5YMes2+HupkVtQhd9+uBe/+3QfSqr19m4aERG5GAYZuiHDI/3x3dOjMeP2SKiUEjYcLkTSG1ux4VCBvZtGREQuhEGGbliQjwbzJw7EujkjER3qg1JdHR7/dB+eWXUQlZcN9m4eERG5AAYZumkDwrRYl3o7nhgTBYUErN53DsmLtuHIxUp7N42IiJycm70bQM5B46ZEenI0xvUPRtqXB3GmrAYPvLMTE2O6QpRKmGA0QaWydyuJiMjZOPwZmfLycgwbNgxDhw7FoEGD8P7779u7SdSCYZH++Cp1JBJ6BUBXZ8QXP57DypNKTPvwR6zLOY+TxdX2biIRETkRhz8j4+Pjg61bt8LT0xM6nQ6DBg3ClClTEBAQYO+mUTN8PVX4ZGYcvv3lIo5cqMCyHSdx4GwFDnyRA0kCnr9nAGaO7GnvZhIRkRNw+CCjVCrh6ekJANDr9RBC8JolnYDaTYH7hnbD3QODEVx1DHnKHsgt1OHg2XK88PVhnL1Ug7/cMwBKhWTvphIRUSdm066lrVu3YuLEiQgLC4MkSVi7dm2T6y1ZsgSRkZFwd3dHfHw89u7da/F8eXk5YmJi0L17dzz77LMIDAy0ZbOpnQW4A3+9byDW/n4E0pOjAQAf7ziFxz75EZtyC1FbZ7RzC4mIqLOyaZDR6XSIiYnBkiVLml1n5cqVSEtLw/z587F//37ExMQgKSkJRUVF8jp+fn44ePAg8vPzsWLFChQWFtqy2WQjkiThiTFRePPhW6BWKrA5rxgzl2Uj5oUNeOyTbOQWcJYTERG1jU27lpKTk5GcnNziOq+//jpmzZqFGTNmAACWLl2Kb775Bh999BHS09Mt1g0JCUFMTAy2bduGBx54oMnt6fV66PVXrzBbWdnw5WgwGGAwtM+1Tczbaa/tObOmapU8IAgRj8dhZfY5/HC0BBcqLuP7I4XYlFuIaXHheOqO3vDzdM0pTjy2rMdaWY+1ahvWy3q2rJW125REBw04kSQJa9aswaRJk+RldXV18PT0xOrVqy2Wp6SkoLy8HOvWrUNhYSE8PT3h4+ODiooK3H777fj8888xePDgJvezYMECZGZmNlq+YsUKeawNOQ4hgIu1wHdnFcgpazhB6KEUuD1UYGxXE7xdM88QEbm8mpoaTJs2DRUVFdBqtc2uZ9fBviUlJTAajQgJCbFYHhISgtzcXADA6dOn8fjjj8uDfJ988slmQwwAZGRkIC0tTX5cWVmJ8PBwTJgwocVCtIXBYMDGjRsxfvx4qHhxlBZZW6vHAOw+WYYXv8nF0aJqfH9ewk+VHvjn1CEY1qNLxzXYznhsWY+1sh5r1Tasl/VsWStzj0pr2hxk0tPT8fLLL7e4zpEjRxAdHd3WTTcpLi4OOTk5Vq+v0Wig0WgaLVepVO1eZFts01lZU6tR/ULwbZ9gfH+kEK9+l4fjRdWY9uGPGN8/BL9N6IERUYEuM8uJx5b1WCvrsVZtw3pZz1bfsdZoc5CZO3cupk+f3uI6vXr1smpbgYGBUCqVjQbvFhYWIjQ0tK1NIyegVEhIGhiKkb0D8Ze1v2DNgfPYcLgQGw4XIszXHVNu7Y6H4yPQzc/D3k0lIiIH0OYgExQUhKCgoHbZuVqtRmxsLLKysuQxMiaTCVlZWUhNTW2XfVDn5KVxwxtTh+L3iVH4167TWJdzHhcqLmPx5uN4d+sJTB0ejmeTouHrwd+WiIhcmU3HyFRXV+P48ePy4/z8fOTk5MDf3x8REREAgLS0NKSkpGDYsGGIi4vDokWLoNPp5FlM5Nr6hPjgxUmD8Od7+uP7I4VYvvsMdp0sxWe7z2BzbjH+ck9/3Nk/BGo3h7/bBhER2YBNg0x2djbGjh0rPzYPwk1JScGyZcsAAFOnTkVxcTHmzZuHgoICDB06FOvXr280AJhcm7tKiV8NCcOvhoRh14lSpP/nJ5wurcHs5fsR4KXGowmRSBnRA36eans3lYiIOpBNg0xiYqJVtxNITU1lVxJZLSEqAN/8YRTe3nwcq/edQ1GVHm98fxQfbD+JOWN74+HhEfB10evQEBG5Gp6Pp07JW+OG5+6Kxs70O/DWw7cgOtQHVZfr8fdvcxH3t+/x1BcHsPNECe/LRUTk5Bz+ppFELXFTKjAxJgx3D+6K/+w/hw+35yO3oArrci5gXc4FDO7mi4fjInBbL3/0CvK2d3OJiKidMciQU1AqJPx6WDgeiO2On85VYGX2Wfxn/zn8fL4CP6/5GQBwa4QfHonvgXuGdIW7SmnnFhMRUXtgkCGnIkkSYsL9EBPuh7nj+2LFnjPYcaIE2acuYf+Zcuw/U44Xvj6MB2K74ze39UDPQC97N5mIiG4Cgww5rQBvDZ68sw+evLMPiqouY1X2OazYcwbny2vx4fZ8fLQjH+P7h+BXMWEY2y8IPu4cIExE1NkwyJBLCPZxx5yxvfHEmChsPVqMf+06hc15xfJVg91VCtw1MBQTBoZiZJ9AaBlqiIg6BQYZcilKhYSx0cEYGx2Mo4VV+Pf+c9h4uBAni3VYm3MBa3MuwE0hYVhkF4zt17Ben2BvSJJr3OOJiKizYZAhl9U3xAcZyf2Rflc0cs6W46uDF7HlaBFOFuuw+2QZdp8sw8Jvc9HNzwNj+gVhbL9g3N47AJ5qfmyIiBwF/0cmlydJEm6J6IJbIrpgHgbgdKkOW/KKsTmvCLtOlOJ8eS1W7DmDFXvOQO2mwG29AjC2XxDuiA5GjwAOFiYisicGGaLr9AjwQsoIL6SMiERtnRG7T5ZiU24RNucV4dylWmw9WoytR4uR+dVh9Ar0wtjoYNwRHYzhkf685xMRUQdjkCFqgYdaKY+pEULgRHE1NuUWYVNuEbJPXcLJEh1Obs/Hh9vz4aVWYsLAUMwa1QsDwrT2bjoRkUtgkCGykiRJ6B3sg97BPnh8dBQqLxuw/VgJNuUWYUteMUqq9Vhz4DzWHDiP0X2DMH1ED4yICuTF94iIbIhBhugGad1VuHtwV9w9uCtMJoGcc+X4eMcpfPPTBbn7Se2mQHxPf4zpG4TRfYM4A4qIqJ0xyBC1A4VCwq0RXXBrRBc8M6EvPt5xChsOFeBCxWVsO1aCbcdKgG+OIMhHgxFRAVd+AhHu72nvphMRdWoMMkTtrEeAFxbcOxDzJw7A8aJq/HC0GD8cLcbe/DIUV+nlG1oCQPcuHhjZOxAjegdiRFQAfDUcLExE1BYMMkQ2IkkS+oT4oE+IDx4b1QuXDUYcOFOOXSdKsPNEKXLOluPcpVp88eNZfPHjWQBAvxBvBEAB3b7ziI0MQJ9gbygU7IoiImoOgwxRB3FXKZEQFYCEqACkAajW1+PH/DLsOF6CHSdKceRiJfIKqwEosHPtIQCAr4cK8T39cVuvANzWKwDRoT4MNkRE12CQIbITb42bPLUbAEqq9dh5rAhrt+VApwnAz+crUVFrkO8HBTDYEBFdj0GGyEEEemuQPCgU4owJd989HJJCiV8uVGL3yVLsPlmKH/PLGGyIiK7DIEPkoNyUCgwN98PQcD88MSYK9UZTm4LNsMgu6N9VC5WSA4iJyHkxyBB1EjcSbNxVCgzp7ocRUQEY1ScIA8O0vEAfETkVBhmiTqqlYLPnZCn2nylHRa0Be/PLsDe/DIu+PwY3RcNMqiHdfDG4uy+GdPdFv1AfaNwYboioc2KQIXIS1wcbk0ngZIkOP54qww95xfjxVBlKdXU4crESRy5WYmV2w5RvlVJCv1Af9AvRIirYC1FB3hjUzRdhvu68CjEROTwGGSInpVBI6B3sjd7B3ng4LgJCCFysuIyfzlXg5/Pl+Pl8JX4+V45LNQb8cr4Sv5yvtHh9iFaDW8K74NYefrglogsGd/NltxQRORwGGSIXIUkSwvw8EObngbsGhQIAhBA4d6kWv5yvwPGiahwvrsbRwmocLaxCYaUe6w8VYP2hAgCQu6X6h/oguqsPokO1iO7qg2Afd3u+LSJycQwyRC5MkiSE+3s2uudTTV09fj5Xgf1nyrH/zCUcOHMJJdVXu6Vw4Oq6AV5q9O+qxcAwLQaEaTGgqxa9gryh5DRwIuoADDJE1Iin2g3xvQIQ3ysAwNUzN0cuViK3oAq5BQ1/nirRoVRXh+3HS7D9eIn8eneVAv1CG0LNgLCGkBMd6gNPNf/LIaL21Wn+V6mpqUH//v3x61//Gv/4xz/s3Rwil3LtmZsJA0Pl5bV1RhwtrMKRi5U4fLEShy40nLGpqTPi4NlyHDxbfs02gJ6BXhjQVYuBYb7y2ZsgH40d3hEROYtOE2Reeukl3HbbbfZuBhFdw0OtREy4H2LC/eRlJpPA6bIaHL5QiUMXKnD4YiUOX6hEUZUeJ4t1OFmsw9c/XZTXD/LRNHRLdb3aNRUZ4MUrFBORVTpFkDl27Bhyc3MxceJE/PLLL/ZuDhG1QKGQ0DPQCz0DvXDPkK7y8uIqvRxqGs7eVCC/RIfiKj225BVjS16xvK6nWon+XS27pvqG+HDWFBE1YtMgs3XrVrz66qvYt28fLl68iDVr1mDSpEmN1luyZAleffVVFBQUICYmBm+99Rbi4uLk55955hm8+uqr2Llzpy2bS0Q2FOSjwRifIIzpGyQvq6mrR25BFQ5duBpwcq90Te07fQn7Tl+S11UqJEQFNXRNRQZ6oaCiFhfPKlCz/zwGd++C3sHeDDpELsimQUan0yEmJgYzZ87ElClTmlxn5cqVSEtLw9KlSxEfH49FixYhKSkJeXl5CA4Oxrp169C3b1/07dvXqiCj1+uh1+vlx5WVDdfGMBgMMBgM7fK+zNtpr+05M9aqbVytXioJGNzVG4O7egOxYQCAeqMJ+aU1OHKxCkcKqnDkYhUOX6zEpRrDlanh1ddsQYEf1hxq+JsEdO/igaighov69Qr0Qu8gL0QFeUHrobLDu3McrnZc3SzWy3q2rJW125SEEKLd997UjiSpyTMy8fHxGD58OBYvXgwAMJlMCA8Px5NPPon09HRkZGTgs88+g1KpRHV1NQwGA+bOnYt58+Y1uZ8FCxYgMzOz0fIVK1bA09OziVcQkaMTAqioA87XSDivA0ouS/BTA3oTcF4HnNdJqDE2P6bGRyUQ4iEQ7AEEuQsEuTf8GegOuPGemkQOqaamBtOmTUNFRQW0Wm2z69k1yNTV1cHT0xOrV6+2WJ6SkoLy8nKsW7fOYhvLli3DL7/80uKspabOyISHh6OkpKTFQrSFwWDAxo0bMX78eKhUrv2bXmtYq7Zhvax3ba3c3NxQXF2Hk8U6nCiuxoliXcNPiQ6Flfpmt6GQgDBfd/QI8EJkgCd6BHgi8spP9y4eTnPncB5XbcN6Wc+WtaqsrERgYGCrQcaug31LSkpgNBoREhJisTwkJAS5ubk3tE2NRgONpvF0TpVK1e5FtsU2nRVr1Tasl/XMtermr0Y3f2+M6mf5/0nVZQNOFOtwvKga+SXVOFVSg/wSHU6V6lBTZ8S58ss4V34ZO06UWrxOqZDQw98T0V190D9Ui+iuDdfC6d7Fo9Peg4rHVduwXtaz1XesNdocZNLT0/Hyyy+3uM6RI0cQHR3d1k23avr06e2+TSJybj7uKvlmmtcSQqC4Si+HmvySGpy68vdTpTpcNphwskSHkyU6/O/ngqvb07jJt2joFeSFHgGeiPD3Qri/B+8iTmQHbQ4yc+fObTVQ9OrVy6ptBQYGQqlUorCw0GJ5YWEhQkNDm3kVEdHNkyQJwVp3BGvd5SsYm5lMAoVVl3GssLrhKsZXBhyfKK5Glb4eP566hB9PXbpue0CYrwd6XOmi6hHghR7+V/4M8ISXplNc7YKo02nzJysoKAhBQUGtr2gFtVqN2NhYZGVlyWNkTCYTsrKykJqa2i77ICJqK4VCQldfD3T19cDoa6aL19WbcLKkGrkXq3CkoBKnS2pwuqwGp690U50vr8X58lrsvK6bCgACvdUW4abnlRlVvQK94aHmmRyiG2XTXxGqq6tx/Phx+XF+fj5ycnLg7++PiIgIAEBaWhpSUlIwbNgwxMXFYdGiRdDpdJgxY4Ytm0ZE1GZqN0XDXb9DtZiEbvJyIQRKqutwulSH06UNwaYh4DT8/VKNASXVdSiprrO4No5ZV193RAZ4ITLQCz0DPREZ0HBBwXB/T14bh6gVNg0y2dnZGDt2rPw4LS0NQMOspGXLlgEApk6diuLiYsybNw8FBQUYOnQo1q9f32gAMBGRo5IkCUE+GgT5aDAs0r/R8xW1BpwprcHpsoagc6pEh/wSHY4XV6O8xoCLFZdxseIydp0svW67Dd1VkYGeGBjmi1vC/RAR0HDPK607B6ESATYOMomJibBmdndqaiq7kojIafl6qDC4uy8Gd/dt9FxptR6nSq8ONM4vuRp2qvT1cnfVjuOljbYZ7u+B8C4NU8XD/T0R3sUT4f4e6Obnye4qchkcfUZEZEcB3hoEeGsQ26OLxXIhBEp1Dd1VJ4p0OHC2HIcvVODspVqU6epQUWtAxXkDfjlf2eR2u3iqEOzjjiAfDSpr63CuRImPz+1BNz9PdPV1R5ifB8L8zH96IMBL3WmnlZNrY5AhInJAkiQh0FuDQG8NYnv448Hh4fJzOn09zl2qxdmyGpy9VIOzZbU4e6kG5y7V4lxZDar09bhUY8ClGgPyCqvMW0TZ2QrknK1ocn9qNwXCfK8GG/Pf/b3U8HZ3Q3c/T4T5ucPNSS4SSM6DQYaIqJPx0rihX6gP+oX6NHpOCIGKWgMKKi+jqFKP4io93N2Aoz/tQ+/BsSiqNuBieS0uVNTiQvllXKyoRVGVHnX1poYurtKaZvfrppDQvYsHIuTZV57o6uuBUF93dPV1R7CPhkGHOhyDDBGRE5EkCX6eavh5qhF95XJcBoMB9aeAuwaGNHm11Lp6EworL+N8eS0uXgk458trcaG8tqELq9aAc5dqWw07CgkI9Nagq6/7lXBzNeSEahseB2s1nIlF7YpBhojIxandFA2Dhf2bv7GuySRQUHkZp0trcObK7Kuzl2pRUFGLixWXUVh5GQajQFGVHkVVehw813QXFgAEeKnlgBOivRJ0fD2uCUDu8FTz64mswyOFiIhapVBI8viZhKiARs+bTA2DkwsqLuNCRS0KKxumlBeYfyov40J5LfT1JpTq6lCqq8OhC00PVAYArbub5RmdK2d1rj3To3V34wBlYpAhIqKbp1BcvZZOU9PMgavjd8wBp+HPhjM6BVeCz8XyWujqjKi8XI/Ky1XXDFZuzFOtvKbr6uoZnV6BXhgW6Q+1G8fruAIGGSIi6hDXjt/p31Xb7HpVlw3XBJ3LctC5NvSU1xhQU2fEyWIdThbrGm3DR+OG3iHe8lmca/8M1roj0FsNbw3P6DgDBhkiInIoPu4q+Lir0Cek8awss9o645WzOLXXhZ5a5JytQEm1HgfOlLe4H3eVAoHeDWeRgrw1CLzyZ5C3CmcuSeh5sQph/l7w91RDoWDgcVQMMkRE1Ol4qJXoGdhwT6rrmUwChy9W4tyVwcgFlforY3ZqUXhlSnq1vh6XDaaGa+9cqm1iD0q8l7ur4W8KCYHeavkCg0Hm8HPtz5VlvMt5x2PFiYjIqSgUEgZ188Wgbk2P1QEazuiUVDfMsCqu0qO4Wo+SKj2KqhoGJR8/V4xaSY0ynQFGk0BhpR6FlfpW9+2hUjYKN009DvBWQ+PGaejtgUGGiIhcjoda2eyUc4PBgP/973+4++6xgEKJ0uo6FFU1XGCwpPpq8CmuuuZxlR66OiNqDUacKavBmbLmLyxo5uuhsgg4YX4e6ObXMCsr0EeDAC81ArzVnIreClaHiIioGSqlomGQsK97q+vq9PUWwaapsGNebjAK+WKDx4uqW9yuh0qJAG81Arw1CLwSbvy9NAj0bvh7gFfDGZ5Abw26eKpdbrYWgwwREVE78NK4wUvjhh4BjcftXMs8Df3aYFNUqb9y24iGmVml1XUoqdZDX29CrcHYwliexrTubgj01liEnADvhuDj79WwLPDKMj8PVacfyMwgQ0RE1IGunYbe0swsIQRq6owNoUanR2l1HUqr9Q0XFKyuQ+mVZSVXlpXp6mA0iSvX4KnHyZLG09Kvp5Agn93x92oINwFeajnomLu3zIHIEaesM8gQERE5IEmS5LM8EQHN3z7CzGRqONNjDjgNgUePkmtCz7WhqKLWAJMASqobur+soXZTXOneagg2XTzcUFGowOBLNegV3PzgaltikCEiInICCoWELl5qdPFSo3dw6+sbjCZc0tVZBJ0S+YyPHmXXPVdTZ0RdvQkXKi7jQsXla/eM1BqDzd5XaxhkiIiIXJBKqUDwlSsdW6Omrt7iTE+prg5FFbXI/iUPYX4eNm5t8xhkiIiIqFWeajd4+rtZTFk3GAwIrz6CAC+13drlWnO0iIiIyKkwyBAREVGnxSBDREREnRaDDBEREXVaDDJERETUaTHIEBERUafFIENERESdVqe4jkxkZCS0Wi0UCgW6dOmCzZs327tJRERE5AA6RZABgJ07d8Lb29vezSAiIiIHwq4lIiIi6rRsGmS2bt2KiRMnIiwsDJIkYe3atU2ut2TJEkRGRsLd3R3x8fHYu3evxfOSJGHMmDEYPnw4li9fbssmExERUSdi064lnU6HmJgYzJw5E1OmTGlynZUrVyItLQ1Lly5FfHw8Fi1ahKSkJOTl5SE4uOH2ndu3b0e3bt1w8eJFjBs3DoMHD8aQIUOa3J5er4def/V25JWVlQAa7gdhMLTP3TnN22mv7Tkz1qptWC/rsVbWY63ahvWyni1rZe02JSGEaPe9N7UjScKaNWswadIki+Xx8fEYPnw4Fi9eDAAwmUwIDw/Hk08+ifT09EbbefbZZzFw4EBMnz69yf0sWLAAmZmZjZavWLECnp6eTbyCiIiIHE1NTQ2mTZuGiooKaLXaZtez62Dfuro67Nu3DxkZGfIyhUKBcePGYdeuXQAazuqYTCb4+PiguroamzZtwoMPPtjsNjMyMpCWliY/rqioQEREBBISEuDj49Mu7TYYDNi8eTPGjh0LlUrVLtt0VqxV27Be1mOtrMdatQ3rZT1b1qqqqgoA0Nr5FrsGmZKSEhiNRoSEhFgsDwkJQW5uLgCgsLAQkydPBgAYjUbMmjULw4cPb3abGo0GGo1GfmzuWurZs2d7N5+IiIhsrKqqCr6+vs0+3+Ygk56ejpdffrnFdY4cOYLo6Oi2brpJvXr1wsGDB2/49WFhYTh79ix8fHwgSVK7tKmyshLh4eE4e/Zsi6e7iLVqK9bLeqyV9VirtmG9rGfLWgkhUFVVhbCwsBbXa3OQmTt3brPjU8x69epl1bYCAwOhVCpRWFhosbywsBChoaFtbVqTFAoFunfv3i7bup5Wq+VBbiXWqm1YL+uxVtZjrdqG9bKerWrV0pkYszYHmaCgIAQFBd1Qg66nVqsRGxuLrKwseRCwyWRCVlYWUlNT22UfRERE5LxsOkamuroax48flx/n5+cjJycH/v7+iIiIAACkpaUhJSUFw4YNQ1xcHBYtWgSdTocZM2bYsmlERETkBGwaZLKzszF27Fj5sXk2UUpKCpYtWwYAmDp1KoqLizFv3jwUFBRg6NChWL9+faMBwI5Eo9Fg/vz5FoOKqWmsVduwXtZjrazHWrUN62U9R6hVh11HhoiIiKi98V5LRERE1GkxyBAREVGnxSBDREREnRaDDBEREXVaDDJERETUaTHItNGSJUsQGRkJd3d3xMfHY+/evfZukt0tWLAAkiRZ/Fx7i4rLly9jzpw5CAgIgLe3N+6///5GV3N2Zlu3bsXEiRMRFhYGSZKwdu1ai+eFEJg3bx66du0KDw8PjBs3DseOHbNYp6ysDI888gi0Wi38/Pzwf//3f6iuru7Ad9ExWqvV9OnTGx1rd911l8U6rlKrhQsXYvjw4fDx8UFwcDAmTZqEvLw8i3Ws+eydOXMG99xzDzw9PREcHIxnn30W9fX1HflWOoQ19UpMTGx0fD3xxBMW67hCvd555x0MGTJEvlpvQkICvv32W/l5RzuuGGTaYOXKlUhLS8P8+fOxf/9+xMTEICkpCUVFRfZumt0NHDgQFy9elH+2b98uP/fHP/4RX331FVatWoUffvgBFy5cwJQpU+zY2o6l0+kQExODJUuWNPn8K6+8gjfffBNLly7Fnj174OXlhaSkJFy+fFle55FHHsGhQ4ewceNGfP3119i6dSsef/zxjnoLHaa1WgHAXXfdZXGsff755xbPu0qtfvjhB8yZMwe7d+/Gxo0bYTAYMGHCBOh0Onmd1j57RqMR99xzD+rq6rBz50588sknWLZsGebNm2ePt2RT1tQLAGbNmmVxfL3yyivyc65Sr+7du+Pvf/879u3bh+zsbNxxxx247777cOjQIQAOeFwJslpcXJyYM2eO/NhoNIqwsDCxcOFCO7bK/ubPny9iYmKafK68vFyoVCqxatUqedmRI0cEALFr164OaqHjACDWrFkjPzaZTCI0NFS8+uqr8rLy8nKh0WjE559/LoQQ4vDhwwKA+PHHH+V1vv32WyFJkjh//nyHtb2jXV8rIYRISUkR9913X7OvcdVaCSFEUVGRACB++OEHIYR1n73//e9/QqFQiIKCAnmdd955R2i1WqHX6zv2DXSw6+slhBBjxowRTz31VLOvceV6denSRXzwwQcOeVzxjIyV6urqsG/fPowbN05eplAoMG7cOOzatcuOLXMMx44dQ1hYGHr16oVHHnkEZ86cAQDs27cPBoPBom7R0dGIiIhg3dBw246CggKL+vj6+iI+Pl6uz65du+Dn54dhw4bJ64wbNw4KhQJ79uzp8Dbb25YtWxAcHIx+/fph9uzZKC0tlZ9z5VpVVFQAAPz9/QFY99nbtWsXBg8ebHEl9aSkJFRWVsq/fTur6+tltnz5cgQGBmLQoEHIyMhATU2N/Jwr1stoNOKLL76ATqdDQkKCQx5XNr1FgTMpKSmB0WhsdOuEkJAQ5Obm2qlVjiE+Ph7Lli1Dv379cPHiRWRmZmLUqFH45ZdfUFBQALVaDT8/P4vXhISEoKCgwD4NdiDmGjR1XJmfKygoQHBwsMXzbm5u8Pf3d7ka3nXXXZgyZQp69uyJEydO4E9/+hOSk5Oxa9cuKJVKl62VyWTC008/jdtvvx2DBg0CAKs+ewUFBU0ee+bnnFVT9QKAadOmoUePHggLC8NPP/2E//f//h/y8vLwn//8B4Br1evnn39GQkICLl++DG9vb6xZswYDBgxATk6Owx1XDDJ005KTk+W/DxkyBPHx8ejRowe+/PJLeHh42LFl5Gweeugh+e+DBw/GkCFDEBUVhS1btuDOO++0Y8vsa86cOfjll18sxqZR85qr17VjqQYPHoyuXbvizjvvxIkTJxAVFdXRzbSrfv36IScnBxUVFVi9ejVSUlLwww8/2LtZTWLXkpUCAwOhVCobjcwuLCxEaGionVrlmPz8/NC3b18cP34coaGhqKurQ3l5ucU6rFsDcw1aOq5CQ0MbDSivr69HWVmZy9ewV69eCAwMxPHjxwG4Zq1SU1Px9ddfY/Pmzejevbu83JrPXmhoaJPHnvk5Z9RcvZoSHx8PABbHl6vUS61Wo3fv3oiNjcXChQsRExODf/7znw55XDHIWEmtViM2NhZZWVnyMpPJhKysLCQkJNixZY6nuroaJ06cQNeuXREbGwuVSmVRt7y8PJw5c4Z1A9CzZ0+EhoZa1KeyshJ79uyR65OQkIDy8nLs27dPXmfTpk0wmUzyf7Su6ty5cygtLUXXrl0BuFathBBITU3FmjVrsGnTJvTs2dPieWs+ewkJCfj5558twt/GjRuh1WoxYMCAjnkjHaS1ejUlJycHACyOL1ep1/VMJhP0er1jHlftPnzYiX3xxRdCo9GIZcuWicOHD4vHH39c+Pn5WYzMdkVz584VW7ZsEfn5+WLHjh1i3LhxIjAwUBQVFQkhhHjiiSdERESE2LRpk8jOzhYJCQkiISHBzq3uOFVVVeLAgQPiwIEDAoB4/fXXxYEDB8Tp06eFEEL8/e9/F35+fmLdunXip59+Evfdd5/o2bOnqK2tlbdx1113iVtuuUXs2bNHbN++XfTp00c8/PDD9npLNtNSraqqqsQzzzwjdu3aJfLz88X3338vbr31VtGnTx9x+fJleRuuUqvZs2cLX19fsWXLFnHx4kX5p6amRl6ntc9efX29GDRokJgwYYLIyckR69evF0FBQSIjI8Meb8mmWqvX8ePHxQsvvCCys7NFfn6+WLdunejVq5cYPXq0vA1XqVd6err44YcfRH5+vvjpp59Eenq6kCRJbNiwQQjheMcVg0wbvfXWWyIiIkKo1WoRFxcndu/ebe8m2d3UqVNF165dhVqtFt26dRNTp04Vx48fl5+vra0Vv//970WXLl2Ep6enmDx5srh48aIdW9yxNm/eLAA0+klJSRFCNEzBfv7550VISIjQaDTizjvvFHl5eRbbKC0tFQ8//LDw9vYWWq1WzJgxQ1RVVdnh3dhWS7WqqakREyZMEEFBQUKlUokePXqIWbNmNfpFwlVq1VSdAIiPP/5YXseaz96pU6dEcnKy8PDwEIGBgWLu3LnCYDB08LuxvdbqdebMGTF69Gjh7+8vNBqN6N27t3j22WdFRUWFxXZcoV4zZ84UPXr0EGq1WgQFBYk777xTDjFCON5xJQkhRPuf5yEiIiKyPY6RISIiok6LQYaIiIg6LQYZIiIi6rQYZIiIiKjTYpAhIiKiTotBhoiIiDotBhkiIiLqtBhkiIiIqNNikCEiIqJOi0GGiIiIOi0GGSIiIuq0/j/l27PNrTpbTwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.yscale(\"symlog\")\n",
    "plt.grid(True)\n",
    "plt.title(\"Power output by size, with symmetric log scale\")\n",
    "power_levels = grid_powers(serial)\n",
    "maxima = [summed_grid_powers(power_levels, i + 1).max() for i in range(300)]\n",
    "_ = plt.plot(np.arange(1, len(maxima) + 1), maxima)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Avoiding the performance issues\n",
    "\n",
    "We are now constantly summing and re-summing values for the larger and larger window sizes. But for any given `size`, we can reuse `size - 1` if we can add the missing column and row values. We can keep a running total\n",
    "of those on hand, with a running sum of `x` and `y`, where `y` is offset by one to avoid summing the bottom-right-hand corner twice.\n",
    "\n",
    "We are essentially adding columns and rows into the original cell in the top left corner, expanding out to the next row and column each step. The cells here are relative to a given top-left corner cell:\n",
    "\n",
    "| ⋱   | 0                                        | 1                                        | 2                                         | 3                                        |\n",
    "| --- | ---------------------------------------- | ---------------------------------------- | ----------------------------------------- | ---------------------------------------- |\n",
    "| 0   | A                                        | <span style=\"color: firebrick\">B</span>  | <span style=\"color: forestgreen\">C</span> | <span style=\"color: steelblue\">D</span>  |\n",
    "| 1   | <span style=\"color: red\">E</span>        | <span style=\"color: red\">F</span>        | <span style=\"color: forestgreen\">G</span> | <span style=\"color: steelblue\">H</span>  |\n",
    "| 2   | <span style=\"color: limegreen\">I</span>  | <span style=\"color: limegreen\">J</span>  | <span style=\"color: limegreen\">K</span>   | <span style=\"color: steelblue\">L</span>  |\n",
    "| 3   | <span style=\"color: dodgerblue\">M</span> | <span style=\"color: dodgerblue\">N</span> | <span style=\"color: dodgerblue\">O</span>  | <span style=\"color: dodgerblue\">P</span> |\n",
    "\n",
    "1. For window-size 1, only `A` is needed\n",
    "2. For window-size 2 we add in the row cells `E` and `F` and column cell `B`\n",
    "3. For window-size 3, add in row cells `I`, `J` and `K`, and column cells `C` and `G`\n",
    "4. For window-size 3, add in row cells `M`, `N`, `O` and `P`, and column cells `D`, `H` and `L`\n",
    "\n",
    "Note that these are all for just one cell. But for the cell _below_ this one in the table, there are corresponding surrounding cells to add, and what is cell `E` for the one in the top row, is cell `A` in that very row, with `F` being the same as `B` in the row above. The same applies to the top cell in the column to the right; there `A` is the equivalent of `B` for the cell in the top left corner, etc.\n",
    "\n",
    "We can take advantage of that by keeping a running sum for rows and columns where we roll the running sums up (for rows) and left (for columns) before adding their numbers to our summed matrix. We use the original 1x1 window power matrix to supply the values. This matrix is rolled up and left a step each time, so first `F`, then `K` then `P` is in the top-left corner. We then roll the column sums up and add the rolled power values to the column sums and add these to the windowed sums. Then we roll the column sums to the left, add them to the windowed sums, and only then add the rolled power values to the column sums.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def optimal_window_size(serial):\n",
    "    power = grid_powers(serial)\n",
    "    summed = np.copy(power)\n",
    "    best_size = 1\n",
    "    highest_power = summed.max()\n",
    "    best_pos = summed.argmax()\n",
    "    # updated sums for rows and columns to update the current size with\n",
    "    power_rolled = rowsums = colsums = power\n",
    "\n",
    "    for size in range(2, GRIDSIZE + 1):\n",
    "        # matrix to update our running row and column sums;\n",
    "        # roll up and left, and clear the last row and column\n",
    "        power_rolled = np.roll(power_rolled, (-1, -1), (0, 1))\n",
    "        power_rolled[-1, :] = 0\n",
    "        power_rolled[:, -1] = 0\n",
    "        # roll rowsums up, clear last row and add power_rolled values\n",
    "        # Then add these to our summed windows\n",
    "        rowsums = np.roll(rowsums, -1, 0)\n",
    "        rowsums[-1, :] = 0\n",
    "        rowsums += power_rolled\n",
    "        summed += rowsums\n",
    "        # roll the column sums to the left, add these to the summed\n",
    "        # windows first, then add power_rolled to the column sums.\n",
    "        colsums = np.roll(colsums, -1, 1)\n",
    "        colsums[:, -1] = 0\n",
    "        summed += colsums\n",
    "        colsums += power_rolled\n",
    "        # finally, clear the next inner row and column of summed\n",
    "        summed[-(size - 1) :, :] = 0\n",
    "        summed[:, -(size - 1) :] = 0\n",
    "\n",
    "        # now our summed matrix is up to date for this window size. Check if it\n",
    "        # has a larger power output\n",
    "        power_output = summed.max()\n",
    "        if power_output > highest_power:\n",
    "            highest_power = power_output\n",
    "            best_size = size\n",
    "            best_pos = summed.argmax()\n",
    "\n",
    "    y, x = np.unravel_index(best_pos, summed.shape)\n",
    "    return x + 1, y + 1, best_size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(235, 87, 13)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "optimal_window_size(serial)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 ms ± 1.71 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit optimal_window_size(serial)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}