{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "626a5f2f",
   "metadata": {},
   "source": [
    "# Stand-alone tools\n",
    "Along with the core functionality of the [Lines](../api/sfrmaker.lines.rst#sfrmaker.lines.Lines) and [SFRData](../api/sfrmaker.sfrdata.rst#sfrmaker.sfrdata.SFRData) classes, SFRmaker also has stand-alone tools that may be useful in work related to the SFR package or other model components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9ff268d2",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-04-03T00:33:30.628999Z",
     "iopub.status.busy": "2025-04-03T00:33:30.628848Z",
     "iopub.status.idle": "2025-04-03T00:33:31.691922Z",
     "shell.execute_reply": "2025-04-03T00:33:31.691383Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import flopy\n",
    "import matplotlib.pyplot as plt\n",
    "import sfrmaker\n",
    "from sfrmaker.routing import find_path, get_upsegs, make_graph, make_reverse_graph"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2430b24",
   "metadata": {},
   "source": [
    "**Set up the Tyler Forks tests model SFR package**  \n",
    "See [Using SFRmaker in a scripting context](SFRmaker_demo.ipynb) for more details"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "1fde68c1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-04-03T00:33:31.693969Z",
     "iopub.status.busy": "2025-04-03T00:33:31.693613Z",
     "iopub.status.idle": "2025-04-03T00:33:33.598507Z",
     "shell.execute_reply": "2025-04-03T00:33:33.597949Z"
    },
    "tags": [
     "\"hide-output\""
    ]
   },
   "outputs": [],
   "source": [
    "%%capture\n",
    "lns = sfrmaker.Lines.from_nhdplus_v2(NHDPlus_paths='../tylerforks/NHDPlus/',\n",
    "                            bbox_filter='../tylerforks/grid.shp')\n",
    "m = flopy.modflow.Modflow.load('tf.nam', model_ws='../tylerforks/tylerforks')\n",
    "mg = flopy.discretization.StructuredGrid(delr=m.dis.delr.array * .3048,  # cell spacing along a row\n",
    "                                         delc=m.dis.delc.array * .3048,  # cell spacing along a column\n",
    "                                         xoff=682688, yoff=5139052,  # lower left corner of model grid\n",
    "                                         angrot=0,  # grid is unrotated\n",
    "                                         # projected coordinate system of model (UTM NAD27 zone 15 North)\n",
    "                                         crs=26715\n",
    "                                         )\n",
    "sfr = lns.to_sfr(grid=mg, model=m, model_length_units='feet')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54936b76",
   "metadata": {},
   "source": [
    "## Methods for working with the SFR routing network, or any other directed acyclic graph (DAG)\n",
    "\n",
    "* the `sfrmaker.routing.find_path` function to trace a routing path from a given reach or segment to the outlet (or for example in the context of GSFLOW, a cascade path from an HRU to an outlet)\n",
    "\n",
    "* the `sfrmaker.routing.get_upsegs` method to get a list of all segments or reaches upstream of a point within the SFR routing network\n",
    "\n",
    "\n",
    "### Create a routing dictionary\n",
    "* from vectors of to/from nodes, which in this case are the segments\n",
    "* if the `one_to_many=True` (the default), a set of one or more downstream connections is returned for each node\n",
    "* if `one_to_many=False`, a single integer representing the downstream connection is returned\n",
    "* SFRmaker's use of dictionaries to work with DAGs was inspired by [this essay](https://www.python.org/doc/essays/graphs/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f5867efe",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-04-03T00:33:33.600766Z",
     "iopub.status.busy": "2025-04-03T00:33:33.600308Z",
     "iopub.status.idle": "2025-04-03T00:33:33.605499Z",
     "shell.execute_reply": "2025-04-03T00:33:33.605025Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{4}"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "routing = make_graph(sfr.segment_data['nseg'], sfr.segment_data['outseg'])\n",
    "routing[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "edb2d15c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-04-03T00:33:33.606962Z",
     "iopub.status.busy": "2025-04-03T00:33:33.606703Z",
     "iopub.status.idle": "2025-04-03T00:33:33.610030Z",
     "shell.execute_reply": "2025-04-03T00:33:33.609669Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "routing = make_graph(sfr.segment_data['nseg'], sfr.segment_data['outseg'], one_to_many=False)\n",
    "routing[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c6f7e91",
   "metadata": {},
   "source": [
    "**The routing dictionary produced by ``make_graph`` can also be reversed:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d4fa3676",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-04-03T00:33:33.611697Z",
     "iopub.status.busy": "2025-04-03T00:33:33.611320Z",
     "iopub.status.idle": "2025-04-03T00:33:33.614428Z",
     "shell.execute_reply": "2025-04-03T00:33:33.614061Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{np.int64(1), np.int64(2)}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "routing_r = make_reverse_graph(routing)\n",
    "routing_r[4]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d89ce32",
   "metadata": {},
   "source": [
    "### Finding a path through the routing network\n",
    "The returned path includes all nodes between the starting point (in this case 1) and the outlet (0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a64e5605",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-04-03T00:33:33.615884Z",
     "iopub.status.busy": "2025-04-03T00:33:33.615737Z",
     "iopub.status.idle": "2025-04-03T00:33:33.618720Z",
     "shell.execute_reply": "2025-04-03T00:33:33.618369Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1, 4, 6, 10, 14, 26, 42, 0]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "path = find_path(routing, 1)\n",
    "path"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f07f41c",
   "metadata": {},
   "source": [
    "#### Plotting an elevation profile along a path\n",
    "This can be useful, for example, when an elevation profile along a stream is desired, or if one wants to map the SFR reach that each HRU ultimately discharges to in a GSFLOW model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b2250a66",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-04-03T00:33:33.620244Z",
     "iopub.status.busy": "2025-04-03T00:33:33.620056Z",
     "iopub.status.idle": "2025-04-03T00:33:33.718989Z",
     "shell.execute_reply": "2025-04-03T00:33:33.718494Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'Distance in River miles')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAGwCAYAAABIC3rIAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAXWhJREFUeJzt3Xd4VHWixvHvpPdJQkiDBIJIb4HQVUSRohQbiAhiWVzLiigWsIHuKqtr25VrwbWsrgULICqigEoRqSH0TiAJSUhCkkkjdeb+kXV2YwJmYMKZJO/neea5zDlnzrwz12Vezvmd8zPZbDYbIiIiIs2Ym9EBRERERIymQiQiIiLNngqRiIiINHsqRCIiItLsqRCJiIhIs6dCJCIiIs2eCpGIiIg0ex5GB2gsrFYr6enpBAYGYjKZjI4jIiIi9WCz2SgsLCQ6Oho3t9MfB1Ihqqf09HRiYmKMjiEiIiJnITU1ldatW592vQpRPQUGBgLVX2hQUJDBaURERKQ+CgoKiImJsf+On44KUT39eposKChIhUhERKSR+b3hLhpULSIiIs2eCpGIiIg0eypEIiIi0uypEImIiEizp0IkIiIizZ4KkYiIiDR7KkQiIiLS7KkQiYiISLOnQiQiIiLNngqRiIiINHsqRCIiItLsqRCJiIhIs6dCZLDSiipWH8g2OoaIiEizpkJkoIoqK3/6KJFb3t3EhxuPGR1HRESk2VIhMpC7yUR0sC82Gzy2eBdvrj5sdCQREZFmSYXIQG5uJp4a25W7L70AgHnf7mONTp+JiIicdypEBjOZTDw8shOTB8QC8Px3+7DZbAanEhERaV5UiFzE/cM64O/lzq7jBSzflWl0HBERkWZFhchFtAjw5vaL4gB44fv9VFl1lEhEROR8USFyIX+4pB3Bfp4czi5mUWKa0XFERESaDRUiFxLk48ldQ6oHWL+y8iBllVUGJxIREWkeVIhczNRBbYkI8uZ4/ik+2ZRqdBwREZFmQYXIxfh4unPvZRcC8OoPB8kvKTc4kYiISNOnQuSCbugbQ7uW/uQUlfPYkl26DF9ERKSBqRC5IE93N16e0AsPNxPf7MhgSdJxoyOJiIg0aSpELqpnTDD3XV596uzJJbtJzS0xOJGIiEjTpULkwu669AJ6xwZTWFbJzE+3695EIiIiDUSFyIV5uLvx8g298PdyZ9PRXNYfzjE6koiISJOkQuTi2rTwZ3SPaABW7c0yOI2IiEjTpELUCFzeORyARYlpHMkuMjiNiIhI06NC1AgM7RROr5hgCkorue29zeQV695EIiIizqRC1Ah4urvx1s0JtAr25ejJEuYs3W10JBERkSZFhaiRaBnozeuTe+NmgqXb01l3UAOsRUREnEWFqBHp0TqYKQPaAPDEl7s0+auIiIiTqBA1MjNHdKRloDfJOcV8mZRudBwREZEmQYWokQny8eTWwW0B+HhTirFhREREmggVokZofJ8YPNxMbEvJ59udGUbHERERafRUiBqhloHe9qNED3y6nd3pFmMDiYiINHIqRI3UIyM7cfGFYZyqqGLav7aQq3sTiYiInDUVokbKw92N+Tf2pl2YP+mWUl74fr/RkURERBotFaJGzOznyV+v6wHAJ5tS2JdZYHAiERGRxkmFqJHrFxfKld0jsdrgz1/vwWazGR1JRESk0VEhagJmj+qMl7sbPx86yaq9WUbHERERaXRUiJqAmFA/br84DoBnlu2lvNJqcCIREZHGRYWoibj70gsIC6i+g/UHG44ZHUdERKRRUSFqIgJ9PHloRAcAXll5gJyiMoMTiYiINB4qRE3I9X1i6BodRGFpJc99u8/oOCIiIo2GClET4u5m4s9XdwPgs61pbD2Wa3AiERGRxkGFqInpHRvCDQkxADyxZDeVVRpgLSIi8ntUiJqgh0d2xOzryZ6MAj7cmGJ0HBEREZdnaCFas2YNY8aMITo6GpPJxJIlS2qsnzt3Lp06dcLf35+QkBCGDRvGxo0ba2xz6aWXYjKZajwmTpxYY5u8vDymTJmC2WzGbDYzZcoU8vPzG/jTGadFgDcPjegIwAvf7ye7UAOsRUREzsTQQlRcXEzPnj2ZP39+nes7dOjA/Pnz2blzJ+vWraNt27YMHz6c7OzsGttNmzaNjIwM++PNN9+ssX7SpEkkJSWxfPlyli9fTlJSElOmTGmwz+UKbuwXS4/W5uoB1ss1wFpERORMTDYXmevBZDKxePFirr766tNuU1BQgNlsZuXKlVx++eVA9RGiXr168corr9T5mr1799KlSxc2bNhA//79AdiwYQMDBw5k3759dOzYsV75fn1vi8VCUFCQQ5/NKFuP5XHd6+vx8nAj6ckr8PPyMDqSiIjIeVXf3+9GM4aovLycBQsWYDab6dmzZ411H374IWFhYXTt2pUHH3yQwsJC+7pffvkFs9lsL0MAAwYMwGw2s379+tO+X1lZGQUFBTUejU3v2GBah/hSXmll7cEco+OIiIi4LJcvRF9//TUBAQH4+Pjw8ssvs2LFCsLCwuzrb7rpJj7++GN++uknnnjiCb744guuvfZa+/rMzEzCw8Nr7Tc8PJzMzMzTvu+8efPsY47MZjMxMTHO/WDngclk4oouEQC8svIgpRVVBicSERFxTS5fiIYOHUpSUhLr169n5MiRTJgwgays/05gOm3aNIYNG0a3bt2YOHEin3/+OStXriQxMdG+jclkqrVfm81W5/JfzZ49G4vFYn+kpqY694OdJ3cNuYBQfy/2ZhTwV92sUUREpE4uX4j8/f1p3749AwYM4O2338bDw4O33377tNv37t0bT09PDh48CEBkZCQnTpyotV12djYRERGn3Y+3tzdBQUE1Ho1ReJAPL46vPsX43vqjrNhT+7sQERFp7ly+EP2WzWajrOz0l5Hv3r2biooKoqKiABg4cCAWi4VNmzbZt9m4cSMWi4VBgwY1eF5XMLRTOLdfFAfAQ59vJ8NyyuBEIiIirsXQQlRUVERSUhJJSUkAJCcnk5SUREpKCsXFxTz66KNs2LCBY8eOkZiYyB/+8AfS0tIYP348AIcPH+bpp59my5YtHD16lGXLljF+/Hji4+MZPHgwAJ07d2bkyJFMmzaNDRs2sGHDBqZNm8bo0aPrfYVZU/DwyI50axVEfkkFc5fuNjqOiIiISzG0EG3ZsoX4+Hji4+MBeOCBB4iPj+fJJ5/E3d2dffv2cd1119GhQwdGjx5NdnY2a9eupWvXrgB4eXmxatUqRowYQceOHZk+fTrDhw9n5cqVuLu729/nww8/pHv37gwfPpzhw4fTo0cPPvjgA0M+s1G8Pdx5cXwv3Ezw3e4TbDmqec5ERER+5TL3IXJ1jfE+RHWZ9cUOPtmcSu/YYL64a9AZB5aLiIg0dk3uPkTiHA9c0QFfT3cSU/L5Mind6DgiIiIuQYWomQkP8uFPl7UH4C/f7MVyqsLgRCIiIsZTIWqGpl3cjgta+pNTVMZL3+83Oo6IiIjhVIiaIS8PN/48rhsAH2w4xs40i8GJREREjKVC1EwNah/GuF7RWG3w+JKdVFk1tl5ERJovFaJm7LErOxPo7cH2NAsfbUoxOo6IiIhhVIiasfAgH2YO7wDAC9/tJ7e43OBEIiIixlAhauamDGxL56ggLKcq+Nt3GmAtIiLNkwpRM+fuZuKpsdV3/v5kc4oGWIuISLOkQiT0iwtlXK9obDZ4cukurBpgLSIizYwKkQAwe1Rn/Lzc2ZaSz6Jtx42OIyIicl6pEAkAkWYf7r3sQgD++u0+Ckp1B2sREWk+VIjE7raL2hIXVn0H63+sPGh0HBERkfNGhUjsvD3cmTOmCwDvrT/KoaxCgxOJiIicHypEUsOlHcMZ1jmCSquNuUv3YLNpgLWIiDR9KkRSy5Oju+Dl4ca6Qzl8tzvT6DgiIiINToVIaolt4ccfL2kHwJ+/3sup8iqDE4mIiDQsFSKp092Xtifa7MPx/FO8sfqw0XFEREQalAqR1MnXy53HrqoeYP3G6sPkFJUZnEhERKThqBDJaV3ZPZILWvpTVmllR1q+0XFEREQajAqRnJbJZKJzVBAAH29KpaLKanAiERGRhqFCJGc0ZUAbPN1NrNhzgvsXJukyfBERaZJUiOSM+rdrwRuT++DpbuLrHRksStQ8ZyIi0vQ4XIjWrFlDZWVlreWVlZWsWbPGKaHEtVzeOYL7r+gAwDPL9pJXXG5wIhEREedyuBANHTqU3NzcWsstFgtDhw51SihxPdMubkeHiAByi8v567f7jI4jIiLiVA4XIpvNhslkqrX85MmT+Pv7OyWUuB5PdzeevaY7AAu3pLL5aO1SLCIi0lh51HfDa6+9Fqi+8uiWW27B29vbvq6qqoodO3YwaNAg5ycUl5HQNpSJfWP4ZHMqjy3eydf3XoyXh4ahiYhI41fvQmQ2m4HqI0SBgYH4+vra13l5eTFgwACmTZvm/ITiUmaN6sSKPSc4cKKIf647wt2Xtjc6koiIyDmrdyF69913AWjbti0PPvigTo81U8F+Xjx2VWce+HQ7/1h1kHG9WtEq2Pf3XygiIuLCHD7fMWfOHLy9vVm5ciVvvvkmhYWFAKSnp1NUVOT0gOJ6rolvRb+4UEorrMxbttfoOCIiIufM4UJ07Ngxunfvzrhx47jnnnvIzs4G4Pnnn+fBBx90ekBxPSaTibljuuJmgq93ZLDhyEmjI4mIiJwThwvRfffdR0JCAnl5eTXGEV1zzTWsWrXKqeHEdXWJDmJS/1gAnl+uy/BFRKRxc7gQrVu3jscffxwvL68ay9u0acPx47qLcXNy72UXArAtNR9LSYXBaURERM6ew4XIarVSVVVVa3laWhqBgYFOCSWNQ0SQDxe09Mdmg88T04yOIyIictYcLkRXXHEFr7zyiv25yWSiqKiIOXPmcOWVVzozmzQCExJiAPjLN3v4anu6wWlERETOjsnm4PTl6enpDB06FHd3dw4ePEhCQgIHDx4kLCyMNWvWEB4e3lBZDVVQUIDZbMZisRAUFGR0HJdhs9l4fMkuPtyYgoebiU/uGEBC21CjY4mIiAD1//2u932IfhUdHU1SUhIff/wxiYmJWK1Wbr/9dm666aYag6yleTCZTPx5XDfySspZtjOTN1Yf4Z8qRCIi0sg4fISoudIRojM7lFXIsJfW4GaCn2ddRpRZ5VhERIxX39/vs5qI6oMPPuCiiy4iOjqaY8eOAfDyyy/z5Zdfnl1aafTahwfSLy4Uqw3eW3/U6DgiIiIOcbgQvf766zzwwAOMGjWKvLw8+xVnISEhNQZbS/Nzx8XtAHj356Ok5ZUYnEZERKT+HC5Er776Km+99RaPPfYYHh7/HYKUkJDAzp07nRpOGpfLO4czoF0o5ZVW7vloG6UVtW/PICIi4oocLkTJycnEx8fXWu7t7U1xcbFTQknjZDKZeO66HgT7ebI9NZ+HP9+BhqiJiEhj4HAhiouLIykpqdbyb7/9li5dujgjkzRibVr48/pNffBwM7F0ezoL1hwxOpKIiMjvcrgQPfTQQ9xzzz0sXLgQm83Gpk2beOaZZ3j00Ud56KGHGiKjNDIDL2jBU+O6AvDyygNkWkoNTiQiInJmDt+H6NZbb6WyspKHH36YkpISJk2aRKtWrfj73//OxIkTGyKjNEKT+sWyOPE4W47l8eL3+/nb+J5GRxIRETmtet2HaOnSpYwaNQpPT88ay3NycrBarU327tT/S/chcty2lDyueW09JhN8fe9FdI02Gx1JRESaGafeh+iaa64hPz8fAHd3d7KysgAICwtrFmVIzk58bAhjekZjs8Gzy/ZqgLWIiLisehWili1bsmHDBqB67iqTydSgoaTpeHhER7zc3fj50El+3J9ldBwREZE61asQ3XnnnYwbNw53d3dMJhORkZG4u7vX+RD5XzGhftw6uC0Azy7bR2WV1dhAIiIidajXoOq5c+cyceJEDh06xNixY3n33XcJDg5u4GjSVNw9tD2fbknlUFYRn2xOZfKANkZHEhERqaHeV5l16tSJTp06MWfOHMaPH4+fn19D5pImxOzryYxhHZizdDcvrzjAuF7RBPp4/v4LRUREzhOH70M0Z84clSFx2KT+sbQL8+dkcTmv/3TY6DgiIiI1nNVs9yKO8nR345FRnQD44JdjWK264kxERFyHCpGcN5d1CsfT3URhWSV7MwuMjiMiImJnaCFas2YNY8aMITo6GpPJxJIlS2qsnzt3Lp06dcLf35+QkBCGDRvGxo0ba2xTVlbGvffeS1hYGP7+/owdO5a0tLQa2+Tl5TFlyhTMZjNms5kpU6bY76sk54+nuxv94kIBuP29LRzJLjI4kYiISDVDC1FxcTE9e/Zk/vz5da7v0KED8+fPZ+fOnaxbt462bdsyfPhwsrOz7dvMmDGDxYsX88knn7Bu3TqKiooYPXo0VVVV9m0mTZpEUlISy5cvZ/ny5SQlJTFlypQG/3xS28s39OLC8AAyC0q58a0N5JeUGx1JRESkflN3/NaqVatYtWoVWVlZWK017yvzzjvvnF0Qk4nFixdz9dVXn3abX2+/vXLlSi6//HIsFgstW7bkgw8+4IYbbgAgPT2dmJgYli1bxogRI9i7dy9dunRhw4YN9O/fH4ANGzYwcOBA9u3bR8eOHeuVT1N3OM/JojLGv/ELR3KKGd+nteY5ExGRBuPUqTv+11NPPcXw4cNZtWoVOTk55OXl1Xg0lPLychYsWIDZbKZnz+of0K1bt1JRUcHw4cPt20VHR9OtWzfWr18PwC+//ILZbLaXIYABAwZgNpvt29SlrKyMgoKCGg9xjhYB3vxtfA9MJvhsaxo/H8oxOpKIiDRzDs92/8Ybb/Dee++dt1NOX3/9NRMnTqSkpISoqChWrFhBWFgYAJmZmXh5eRESElLjNREREWRmZtq3qWu+tfDwcPs2dZk3bx5PPfWUEz+J/K8+bUKZMqAN7/9yjEcX7+S7GZfg46k7nYuIiDEcPkJUXl7OoEGDGiJLnYYOHUpSUhLr169n5MiRTJgwwT657On8dr61uuZe+7052WbPno3FYrE/UlNTz/5DSJ0eGtGR8EBvjp0sYc2B7N9/gYiISANxuBD94Q9/4KOPPmqILHXy9/enffv2DBgwgLfffhsPDw/efvttACIjIykvL691qi4rK4uIiAj7NidOnKi13+zsbPs2dfH29iYoKKjGQ5wr0MeTyzpVH71LTMk3NoyIiDRrDp8yKy0tZcGCBaxcuZIePXrg6VlzCoaXXnrJaeHqYrPZKCsrA6BPnz54enqyYsUKJkyYAEBGRga7du3i+eefB2DgwIFYLBY2bdpEv379ANi4cSMWi+W8HumSuvVuE8Inm1P5fGsaE/vG0DbM3+hIIiLSDDlciHbs2EGvXr0A2LVrV411ZzoFVZeioiIOHTpkf56cnExSUhKhoaG0aNGCZ555hrFjxxIVFcXJkyd57bXXSEtLY/z48QCYzWZuv/12Zs6cSYsWLQgNDeXBBx+ke/fuDBs2DIDOnTszcuRIpk2bxptvvgnAHXfcwejRo+t9hZk0nKu6R/HOumT2ZRYy6a0NLPzjQGJCNTWMiIicX2d12b2z/PTTTwwdOrTW8qlTp/LGG28wadIkNm7cSE5ODi1atKBv3748/vjj9O3b175taWkpDz30EB999BGnTp3i8ssv57XXXiMmJsa+TW5uLtOnT2fp0qUAjB07lvnz5xMcHFzvrLrsvuFkF5YxccEvHM4upkNEAMvvuwQ3N8fKtYiISF3q+/ttaCFqTFSIGlampZQB81YBsOXxYYQFeBucSEREmoL6/n7X65TZtddey3vvvUdQUBDXXnvtGbddtGiRY0lFgEizD2ZfTyynKjiQWUhYexUiERE5f+p1lZnZbLaPD/p1PrDTPUTOVt+21fOc3fHBVn45fNLgNCIi0pzolFk96ZRZwysorWDav7awMTkXLw83ltw9mC7R+q5FROTsNdjUHSINJcjHk3/d1o9LOrSkvNLK3KW7UV8XEZHzQYVIXIqPpzt/vbY7Pp5ubDqay1c7MoyOJCIizYAKkbic6GBf7r60PQDzlu2lpLzS4EQiItLUqRCJS7rjkna0DvElw1LKGz8dNjqOiIg0cSpE4pJ8PN15/KrOALyx5gipuSUGJxIRkabMqYXo/fff5/Bh/WtenGNE10gGt29BeaWVZ77Za3QcERFpwpxaiG655Ra6dOnCvffe68zdSjNlMpmYM6Yr7m4mlu/OZFtKntGRRESkiXJqIbJarezfv59u3bo5c7fSjHWICGRcr2gA/r0hxeA0IiLSVDl9DFHbtm354x//6OzdSjM2eUAbAL7ekU5+SbnBaUREpCmq11xmv2W1Wjl06BBZWVlYrdYa6y655BKnBBP5VXxMMF2igtiTUcBba4/w0IhORkcSEZEmxuFCtGHDBiZNmsSxY8dq3UXYZDJRVVXltHAiUP3f1X3DLuSPH2zl7XXJjO3Zio6RgUbHEhGRJsThU2Z33nknCQkJ7Nq1i9zcXPLy8uyP3NzchsgowvAuEVzUPozSCivT3t9CXrFOnYmIiPM4PLmrv78/27dvp3379g2VySVpclfj5RWXM/b/1pGae4ohHVry3q19MZlMRscSEREX1mCTu/bv359Dhw6dUziRsxHi78VbNyfg5eHG6gPZrD2YY3QkERFpIhweQ3Tvvfcyc+ZMMjMz6d69O56enjXW9+jRw2nhRH6rU2QQUwa04e11yfztu/1cfGGYjhKJiMg5c/iUmZtb7YNKJpMJm83WpAdV65SZ6zhZVMYlz/9IcXkVb09N4PLOEUZHEhERF1Xf32+HjxAlJyefUzCRc9UiwJsb+sbyzs/JfJmUrkIkIiLnzOFC1KZNm4bIIeKQ0T2jeOfnZH7Yl0VpRRU+nu5GRxIRkUbsrO5UffjwYe69916GDRvGFVdcwfTp0zWpq5xXvVoH0yrYl6KySv723X6j44iISCPncCH67rvv6NKlC5s2baJHjx5069aNjRs30rVrV1asWNEQGUVqcXMz8cTozgC8vS6Zz7akGpxIREQaM4cHVcfHxzNixAj++te/1lg+a9Ysvv/+exITE50a0FVoULVremnFAf6x6iBe7m58fEd/+rQJNTqSiIi4kAa7D9HevXu5/fbbay2/7bbb2LNnj6O7EzknMy6/kJFdIymvsvLHDxKxlFQYHUlERBohhwtRy5YtSUpKqrU8KSmJ8PBwZ2QSqTc3NxMvTuhJXJg/OUVlfLUj3ehIIiLSCDl8ldm0adO44447OHLkCIMGDcJkMrFu3Tqee+45Zs6c2RAZRc7I39uDG/vF8OyyfSzdns7kAboSUkREHOPwGCKbzcYrr7zCiy++SHp69b/Go6Ojeeihh5g+fXqTvWuwxhC5tvT8Uwz66w8ArJ91GdHBvgYnEhERV1Df32+HC9H/KiwsBCAwMPBsd9FoqBC5vglv/sKm5FwevbITd1xygdFxRETEBTTYoOr/FRgY2CzKkDQOY3tGA/BlksYRiYiIY+o1hqh3796sWrWKkJAQ4uPjz3harKledi+u78ruUTz11W52pxewPTWfnjHBRkcSEZFGol6FaNy4cXh7e9v/3FTHCUnjFurvxZge0Szadpy31h5h/qTeRkcSEZFG4pzGEDUnGkPUOOxJL+DKf6zF3c3ETw9eSkyon9GRRETEQA02hqhdu3acPHmy1vL8/HzatWvn6O5EnKpLdBAXtQ+jymrj3Z+PGh1HREQaCYcL0dGjR6mqqqq1vKysjLS0NKeEEjkXf7g4DoDPtqaiA6AiIlIf9b4x49KlS+1//u677zCbzfbnVVVVrFq1iri4OOemEzkLAy9oAUBhaSV5JRWE+nsZnEhERFxdvQvR1VdfDYDJZGLq1Kk11nl6etK2bVtefPFFp4YTORveHu60DPQmu7CM5JxiFSIREfld9T5lZrVasVqtxMbGkpWVZX9utVopKytj//79jB49uiGzitRbz9bVRzCf+mo3ZZW1T/GKiIj8L4fHECUnJxMWFtYQWUScZu7YrgT7ebIjzcLTX+0xOo6IiLg4hyd3BSguLmb16tWkpKRQXl5eY9306dOdEkzkXLQO8eOVG3px63ub+XBjCtf2bkWfNqFGxxIRERflcCHatm0bV155JSUlJRQXFxMaGkpOTg5+fn6Eh4erEInLuLRjONf3bs1nW9N4a00yfaaoEImISN0cPmV2//33M2bMGHJzc/H19WXDhg0cO3aMPn368MILLzRERpGzNu2S6ntjfbcnk2Mniw1OIyIirsrhQpSUlMTMmTNxd3fH3d2dsrIyYmJieP7553n00UcbIqPIWesQEciQDi2x2eCddclGxxERERflcCHy9PS0z2UWERFBSkoKAGaz2f5nEVcy7eLqo0SfbkmjsLTC4DQiIuKKHC5E8fHxbNmyBYChQ4fy5JNP8uGHHzJjxgy6d+/u9IAi52pw+xbEhPpyqqKKxJR8o+OIiIgLcrgQPfvss0RFRQHw5z//mRYtWnDXXXeRlZXFggULnB5Q5FyZTCYS/nOF2baUPIPTiIiIK3L4KrOEhAT7n1u2bMmyZcucGkikIcTHBrN423G26QiRiIjUweEjRE899RSHDx9uiCwiDaZXTDAAicfyyC0uP/PGIiLS7DhciL744gs6dOjAgAEDmD9/PtnZ2Q2RS8SpukQFERfmT2FZJXd/uJWKKqvRkURExIU4XIh27NjBjh07uOyyy3jppZdo1aoVV155JR999BElJSUNkVHknHm4u/HmlD74e7mz4UiupvMQEZEaTDabzXYuO/j555/56KOP+OyzzygtLaWgoMBZ2VxKQUEBZrMZi8VCUFCQ0XHkLK3cc4JpH2zBZoPFdw8iPjbE6EgiItKA6vv77fARot/y9/fH19cXLy8vKip0jxdxbcO6RHBF5wgANh/NNTiNiIi4irMqRMnJyTzzzDN06dKFhIQEEhMTmTt3LpmZmc7OJ+J0vWKDAdieZjE2iIiIuAyHL7sfOHAgmzZtonv37tx6661MmjSJVq1aNUQ2kQbRs3UwUH3FWZXVhrubydhAIiJiOIePEA0dOpQdO3aQlJTEQw89dE5laM2aNYwZM4bo6GhMJhNLliyxr6uoqOCRRx6he/fu+Pv7Ex0dzc0330x6enqNfVx66aWYTKYaj4kTJ9bYJi8vjylTpmA2mzGbzUyZMoX8/Pyzzi2NW6+YYAJ9PMiwlPLxJk03IyIiZ3mn6q5du1JeXs7+/fuprKw86zcvLi6mZ8+ezJ8/v9a6kpISEhMTeeKJJ0hMTGTRokUcOHCAsWPH1tp22rRpZGRk2B9vvvlmjfWTJk0iKSmJ5cuXs3z5cpKSkpgyZcpZ55bGzd/bg5lXdADgb9/t132JRETE8VNmp06d4k9/+hP/+te/ADhw4ADt2rVj+vTpREdHM2vWrHrva9SoUYwaNarOdWazmRUrVtRY9uqrr9KvXz9SUlKIjY21L/fz8yMyMrLO/ezdu5fly5ezYcMG+vfvD8Bbb73FwIED2b9/Px07dqx3Xmk6Jg9owyebU9mXWcjibce5/aI4oyOJiIiBHD5CNGvWLLZv385PP/2Ej4+PffmwYcNYuHChU8P9lsViwWQyERwcXGP5hx9+SFhYGF27duXBBx+ksLDQvu6XX37BbDbbyxDAgAEDMJvNrF+//rTvVVZWRkFBQY2HNB0e7m6M61V9undT8kmD04iIiNEcPkK0ZMkSFi5cyIABAzCZ/jsYtUuXLg06pUdpaSmzZs1i0qRJNe4jcNNNNxEXF0dkZCS7du1i9uzZbN++3X50KTMzk/Dw8Fr7Cw8PP+NVcfPmzeOpp55y/gcRl9EvrnrC13UHc0jPP0V0sK/BiURExCgOHyHKzs6us2AUFxfXKEjOVFFRwcSJE7Farbz22ms11k2bNo1hw4bRrVs3Jk6cyOeff87KlStJTEy0b1NXLpvNdsa8s2fPxmKx2B+pqanO+0DiEuJjgukVE0xxeRUPf74Dq/Wc7lEqIiKNmMOFqG/fvnzzzTf257+Wil/H5ThbRUUFEyZMIDk5mRUrVvzuXaJ79+6Np6cnBw8eBCAyMpITJ07U2i47O5uIiIjT7sfb25ugoKAaD2la3NxMvDihJz6ebqw7lMMHG44ZHUlERAzi8CmzefPmMXLkSPbs2UNlZSV///vf2b17N7/88gurV692arhfy9DBgwf58ccfadGixe++Zvfu3VRUVBAVFQVU3zfJYrGwadMm+vXrB8DGjRuxWCwMGjTIqXml8bmgZQCzR3VmztLdvLLyADf0jcHH093oWCIicp45fIRo0KBB/Pzzz5SUlHDBBRfw/fffExERwS+//EKfPn0c2ldRURFJSUkkJSUB1XfATkpKIiUlhcrKSq6//nq2bNnChx9+SFVVFZmZmWRmZlJeXn2Z9OHDh3n66afZsmULR48eZdmyZYwfP574+HgGDx4MQOfOnRk5ciTTpk1jw4YNbNiwgWnTpjF69GhdYSYA3NQ/ltYhvuSVVPBFYprRcURExADnPLnrufjpp58YOnRoreVTp05l7ty5xMXVfSn0jz/+yKWXXkpqaiqTJ09m165dFBUVERMTw1VXXcWcOXMIDQ21b5+bm8v06dNZunQpAGPHjmX+/Pm1rlY7E03u2rS9vS6ZP3+9h3Yt/Vl5/xDcdPdqEZEmob6/3/UqRI5cct5Uy4IKUdNWVFbJwGdXUVhWydtTE7i88+nHl4mISONR39/veo0hCg4O/t0ryH69aquqqsqxpCIuIMDbgxv7x7JgzRH+uTZZhUhEpJmpVyH68ccfGzqHiOFuGdSWt9cl88uRk+w6bqFbK7PRkURE5DypVyEaMmRIQ+cQMVx0sC9XdY9i6fZ03l6XzMs39DI6koiInCcOX2UGsHbtWiZPnsygQYM4fvw4AB988AHr1q1zajiR8+0PF1cP5P9qezoZllMGpxERkfPF4UL0xRdfMGLECHx9fUlMTKSsrAyAwsJCnn32WacHFDmferQOpl9cKJVWG/9arxs1iog0Fw4Xor/85S+88cYbvPXWW3h6etqXDxo0qMZ0GSKN1bSL2wHw0cZjFJdVGpxGRETOB4cL0f79+7nkkktqLQ8KCiI/P98ZmUQMdXmncOLC/CkoreSzLZrDTkSkOXC4EEVFRXHo0KFay9etW0e7du2cEkrESG5uJm67qHos0Ts/H6VKk76KiDR5DheiP/7xj9x3331s3LgRk8lEeno6H374IQ8++CB33313Q2QUOe+u792aYD9PUnJLWLEn0+g4IiLSwBye3PXhhx/GYrEwdOhQSktLueSSS/D29ubBBx/kT3/6U0NkFDnvfL3cmdy/DfN/PMQ/1yYzsluU0ZFERKQBnfVcZiUlJezZswer1UqXLl0ICAhwdjaXoqk7mp+sglIueu5HyqusLL57EPGxIUZHEhERB9X39/us7kME4OfnR0JCAv369WvyZUiap/AgH8b2igbgn+uSDU4jIiIN6awLkUhz8OuNGr/dmUFaXonBaUREpKGoEImcQafIIAa0C8Vqg2U7M4yOIyIiDUSFSOR3XNm9ekD1ij0nDE4iIiINRYVI5HcM6xwBwJZjeeQUlRmcRkREGkK9LrtfunRpvXc4duzYsw4j4oqig33p3srMzuMWFiWmccclFxgdSUREnKxehejqq6+u8dxkMvG/V+ubTCb7n6uqqpyTTMSFTBnYhoc/38Ebq48wqX8bArwdvoWXiIi4sHqdMrNarfbH999/T69evfj222/Jz8/HYrGwbNkyevfuzfLlyxs6r4ghro1vRbswf3KLy3lXl+CLiDQ5Do8hmjFjBn//+98ZMWIEQUFBBAYGMmLECF566SWmT5/eEBlFDOfh7sZ9wy4E4F+/aH4zEZGmxuFCdPjwYcxmc63lZrOZo0ePOiOTiEu6snsUwX6e5BSVszH5pNFxRETEiRwuRH379mXGjBlkZPz3niyZmZnMnDmTfv36OTWciCvxdHdjRJdIAL7anm5wGhERcSaHC9E777xDVlYWbdq0oX379rRv357Y2FgyMjJ4++23GyKjiMsY95+pPD7bksb+zEKD04iIiLOc1eSuNpuNFStWsG/fPmw2G126dGHYsGE1rjZrajS5q0D1f/t3fLCVFXtO0KdNCJ/9cSBubk33v3sRkcauvr/fZz3bPUBpaSne3t5Nugj9SoVIfpWef4orXlpNcXkVf7m6G5MHtDE6koiInEaDzXZvtVr585//TKtWrQgICCA5ufoS5CeeeEKnzKRZiA725cERHQH423f7KSmvNDiRiIicK4cL0V/+8hfee+89nn/+eby8vOzLu3fvzj//+U+nhhNxVTcPbEubFn5YTlWwZJsGWIuINHYOF6L333+fBQsWcNNNN+Hu7m5f3qNHD/bt2+fUcCKuyt3NxJT/nCp7/5ejnMOZZxERcQEOF6Ljx4/Tvn37WsutVisVFRVOCSXSGIzvE4Ovpzv7MgvZlJxrdBwRETkHDheirl27snbt2lrLP/vsM+Lj450SSqQxMPt5cnV8KwDe/+WYwWlERORcODxD5Zw5c5gyZQrHjx/HarWyaNEi9u/fz/vvv8/XX3/dEBlFXNbUQW34eFMKy3dnklVYSnigj9GRRETkLDh8hGjMmDEsXLiQZcuWYTKZePLJJ9m7dy9fffUVV1xxRUNkFHFZnSKD6B0bTJXVxudb04yOIyIiZ8nhI0QAI0aMYMSIEc7OItIoTewXS2JKPgs3p3LnJRfoRo0iIo2Qw0eIfrVlyxY++OAD/v3vf7N161ZnZhJpVEb3iCLQ24NjJ0vYoElfRUQaJYePEKWlpXHjjTfy888/ExwcDEB+fj6DBg3i448/JiYmxtkZRVyan5cHY3pF89HGFL7ekcGgC8KMjiQiIg5y+AjRbbfdRkVFBXv37iU3N5fc3Fz27t2LzWbj9ttvb4iMIi7vkgurS9D21Hxjg4iIyFlx+AjR2rVrWb9+PR07drQv69ixI6+++iqDBw92ajiRxqJH62AA9mUWkmE5RZTZ19hAIiLiEIePEMXGxtZ5A8bKykpatWrllFAijU2U2cd+tdkDC7dTZdWdq0VEGhOHC9Hzzz/Pvffey5YtW+zTFWzZsoX77ruPF154wekBRRoDk8nEixN64eflzi9HTrJgzRGjI4mIiANMtnpMwhQSEoLJ9N9LiYuLi6msrMTDo/qM269/9vf3Jze3aU5hUFBQgNlsxmKxEBQUZHQccVGfbk7l4S924OluYuOjwwj19/r9F4mISIOp7+93vcYQvfLKK87KJdKkjU9ozeurD5OcU8zO4xaGdGhpdCQREamHehWiqVOnNnQOkSbBZDLRJTqI5Jxi9mYUqBCJiDQSZ3WnaoCsrCyysrKwWq01lvfo0eOcQ4k0Zl2igvhmRwYr9pzgDxfF4eF+1vc/FRGR88ThQrR161amTp1qv/fQ/zKZTFRVVTktnEhjNLZnNH9feZCtx/K4+rWf+eu1PejWymx0LBEROQOH/+l666230qFDB9avX8+RI0dITk62P44c0ZU1IjGhfrwysRdBPh7sOl7Ata+vZ19mgdGxRETkDOp1ldn/CgwMZNu2bbRv376hMrkkXWUmjsoqLGX6x9vYcCSX/nGh/PsP/fHU6TMRkfOqvr/fDv/tfPnll7N9+/ZzCifSHIQH+vDC+J54e7ixMTmXO97fQlFZpdGxRESkDg4fIcrJyWHq1Kn069ePbt264enpWWP92LFjnRrQVegIkZytH/ad4O4PEymtsNI+PIAFU/rQrmWA0bFERJqF+v5+O1yIli5dypQpUygsLKy9syY8qFqFSM7FtpQ87vz3Vk4UlBHo48GX9wxWKRIROQ8a7JTZ9OnTmTJlChkZGVit1hqPplqGRM5VfGwIX917Ed1aBVFYWsk7PyfXukpTRESM43AhOnnyJPfffz8RERENkUekyQoP9OG+yzsA8O8NKUx7fwvZhWUGpxIRETiLQnTttdfy448/NkQWkSZvWOdw7h/WAU93Eyv3ZnHHB1uorLL+/gtFRKRBOXxjxg4dOjB79mzWrVtH9+7daw2qnj59utPCiTQ1JpOJ+4ZdyPCuEUx44xe2peTz9Nd7mDOmK+5upt/fgYiINAiHB1XHxcWdfmcmU5O9OaMGVYuzfZl0nPs+SQLgii4RzJ8Uj7eHu7GhRESamAYbVP2/d6b+7cPRMrRmzRrGjBlDdHQ0JpOJJUuW2NdVVFTwyCOP0L17d/z9/YmOjubmm28mPT29xj7Kysq49957CQsLw9/fn7Fjx5KWllZjm7y8PKZMmYLZbMZsNjNlyhTy8/Md/egiTjWuVyv+PrEXXh5urNhzgv/74ZDRkUREmi1Db5tbXFxMz549mT9/fq11JSUlJCYm8sQTT5CYmMiiRYs4cOBArfsczZgxg8WLF/PJJ5+wbt06ioqKGD16dI0r3iZNmkRSUhLLly9n+fLlJCUlMWXKlAb/fCK/Z1yvVrw0oScA8388xBurD+vqMxERAzh8ygwgLS2NpUuXkpKSQnl5eY11L7300tkFMZlYvHgxV1999Wm32bx5M/369ePYsWPExsZisVho2bIlH3zwATfccAMA6enpxMTEsGzZMkaMGMHevXvp0qULGzZsoH///gBs2LCBgQMHsm/fPjp27FivfDplJg3FZrPx2JJdfLQxBYDr+7Tmr9d2x0PTfIiInLP6/n47PKh61apVjB07lri4OPbv30+3bt04evQoNpuN3r17n1Po32OxWDCZTAQHBwOwdetWKioqGD58uH2b6OhounXrxvr16xkxYgS//PILZrPZXoYABgwYgNlsZv369actRGVlZZSV/feS6IICTc4pDcNkMvHM1d3oHBXE3KW7+XxrGkWllfz9xl4aUyQicp44/E/Q2bNnM3PmTHbt2oWPjw9ffPEFqampDBkyhPHjxzdERgBKS0uZNWsWkyZNsje8zMxMvLy8CAkJqbFtREQEmZmZ9m3Cw8Nr7S88PNy+TV3mzZtnH3NkNpuJiYlx4qcRqclkMjFlQBteu6k3Xu5uLN+dyZ8+2kaFLskXETkvHC5Ee/fuZerUqQB4eHhw6tQpAgICePrpp3nuueecHhCqB1hPnDgRq9XKa6+99rvb22w2TKb/XsL8v38+3Ta/NXv2bCwWi/2Rmpp6duFFHDCiayRv35JgH2g9Y2ESVqvGFImINDSHC5G/v7/9VFJ0dDSHDx+2r8vJyXFesv+oqKhgwoQJJCcns2LFihrn/yIjIykvLycvL6/Ga7Kysux30o6MjOTEiRO19pudnX3Gu217e3sTFBRU4yFyPlx8YUvenNwHT3cT3+zIYOEWlXERkYbmcCEaMGAAP//8MwBXXXUVM2fO5JlnnuG2225jwIABTg33axk6ePAgK1eupEWLFjXW9+nTB09PT1asWGFflpGRwa5duxg0aBAAAwcOxGKxsGnTJvs2GzduxGKx2LcRcTVDO4Uze1RnAP789R6+2p7+O68QEZFz4fCg6pdeeomioiIA5s6dS1FREQsXLqR9+/a8/PLLDu2rqKiIQ4f+e++V5ORkkpKSCA0NJTo6muuvv57ExES+/vprqqqq7GN+QkND8fLywmw2c/vttzNz5kxatGhBaGgoDz74IN27d2fYsGEAdO7cmZEjRzJt2jTefPNNAO644w5Gjx5d7yvMRIxw88A2/Lg/i7UHc7j3420czi5ixrAORscSEWmSzuqye2f56aefGDp0aK3lU6dOZe7cuae9K/aPP/7IpZdeClQPtn7ooYf46KOPOHXqFJdffjmvvfZajUHQubm5TJ8+naVLlwIwduxY5s+fb79arT502b0YobLKyosrDvD6T4fxcDOx8oEhtA3zNzqWiEijUd/f77MqRPn5+Xz++eccPnyYhx56iNDQUBITE4mIiKBVq1bnFNxVqRCJkSa8+QubknMJ9PHgmWu6M7ZntNGRREQahQabumPHjh106NCB5557jhdeeME+BcbixYuZPXv2WQcWkdN7cXxPesUEU1hayfSPt/H97tPfMkJERBzncCF64IEHuOWWWzh48CA+Pj725aNGjWLNmjVODSci1WJC/fj8zoHc2C8WgBkLk1h/2PlXdYqINFcOF6LNmzfzxz/+sdbyVq1anfFGhyJybjzc3Xh6XFcuvjCMkvIqpr6ziUWJab//QhER+V0OFyIfH586p7HYv38/LVu2dEooEambp7sbb92cwFU9oqiosvHAp9v5eFOK0bFERBo9hwvRuHHjePrpp6moqACq7wKdkpLCrFmzuO6665weUERq8vF059WJ8dw6uC0Asxft5BOVIhGRc+JwIXrhhRfIzs4mPDycU6dOMWTIENq3b09gYCDPPPNMQ2QUkd9wczPx5Ogu9lL02JJdbErONTaUiEgjdtb3Ifrhhx9ITEzEarXSu3dv+40Qmypddi+uyGazcf/CJJYkpRMR5M3H0wbQrmWA0bFERFxGg9yHqLKyEh8fH5KSkujWrZtTgjYWKkTiqorLKhn3fz9zKKuIYD9P3pjchwHtWvz+C0VEmoEGuQ+Rh4cHbdq0oaqq6pwDiohz+Ht78PG0AfSKCSa/pIJb3t3E/sxCo2OJiDQqDo8hevzxx5k9eza5uRqvIOIqWgZ688kdA7iofRilFVbu/nAr2YVlRscSEWk0HB5DFB8fz6FDh6ioqKBNmzb4+9ecVykxMdGpAV2FTplJY5BTVMZV/1jLiYIyWgX78s+pCXSO0n+vItJ81ff32+HZ7seNG4fJZDqncCLSMMICvPlo2gD+8K8tJOcUc/3r6/nXbf1IaBtqdDQREZdm6Gz3jYmOEEljkl9Szp3/3sqGI7n4ebnz3q396BenUiQizU+DTe7arl07Tp48WWt5fn4+7dq1c3R3ItIAgv28ePeWflzUvnqaj1ve3cSWoxr3JyJyOg4XoqNHj9Z5lVlZWRlpaZpXScRV+Hq588+pCfa5z259dzO7jluMjiUi4pLqPYZo6dKl9j9/9913mM1m+/OqqipWrVpFXFycc9OJyDnx8XRnwZQEpr6ziU1Hc5ny9kYW/nEgHSICjY4mIuJS6j2GyM2t+mCSyWTity/x9PSkbdu2vPjii4wePdr5KV2AxhBJY1ZYWsHkf25ke5qF8EBvPv3jQNqG+f/+C0VEGjmnjyGyWq1YrVZiY2PJysqyP7darZSVlbF///4mW4ZEGrtAH0/+dVs/OkUGklVYxs3vbKK0QjdYFRH5lcNjiJKTkwkLC2uILCLSgIL9vPjg9v5EBvmQklvCG6sP1zraKyLSXNW7EG3cuJFvv/22xrL333+fuLg4wsPDueOOOygr051xRVxZy0Bvpl9+IQCvrDzItPe3klVYanAqERHj1bsQzZ07lx07dtif79y5k9tvv51hw4Yxa9YsvvrqK+bNm9cgIUXEeW7sF8MjIzvh5e7Gyr0nuOb/1pOef8roWCIihqp3IUpKSuLyyy+3P//kk0/o378/b731Fg888AD/+Mc/+PTTTxskpIg4j8lk4q5LL2DpvYOJC/PneP4ppry9kZNFOsIrIs1XvQtRXl4eERER9uerV69m5MiR9ud9+/YlNTXVuelEpMF0igzi33/oT7TZh8PZxdzy7mYKSyuMjiUiYoh6F6KIiAiSk5MBKC8vJzExkYEDB9rXFxYW4unp6fyEItJgWgX78sEf+hPq78XO4xamvb9FV5+JSLNU70I0cuRIZs2axdq1a5k9ezZ+fn5cfPHF9vU7duzgggsuaJCQItJwLmgZwPu39SPA24MNR3K558NEyiutRscSETmv6l2I/vKXv+Du7s6QIUN46623eOutt/Dy8rKvf+eddxg+fHiDhBSRhtWtlZm3pybg7eHGqn1Z3L8wicoqlSIRaT4cnu3eYrEQEBCAu7t7jeW5ubkEBATUKElNie5ULc3B6gPZTPvXFsqrrFzbuxUvXN8TNzeT0bFERM5ag812bzaba5UhgNDQ0CZbhkSaiyEdWjJ/UjzubiYWJR7n8S936eaNItIsOFyIRKRpG941kpdv6IXJBB9tTOEv3+xVKRKRJk+FSERqGdszmueu6wHA2+uSeXnFAYMTiYg0LBUiEanThIQYnhrbFYB//HCI1386bHAiEZGGo0IkIqc1dVBbHhnZCYDnlu/jX+uPGhtIRKSBqBCJyBnddekFTL+sPQBzlu7m0y26I72IND0qRCLyu+6/ogO3XxQHwKwvdrBsZ4bBiUREnEuFSER+l8lk4vGrOnNjv1isNnjyy90UlVUaHUtExGlUiESkXkwmE0+P60qU2YecojJGvrKGNQeyjY4lIuIUKkQiUm+e7m7Mn9SbVsG+pOWd4uZ3NvHRxhSjY4mInDMVIhFxSJ82IXx//yXc2C8GgMeW7OTLpOMGpxIROTcqRCLiMH9vD569pjuTB8Ris8HMT7eTnn/K6FgiImdNhUhEzorJZOLpsd1o28KPSquNRxfvxGrVFB8i0jipEInIWXNzMzHtknYA/LQ/mz5/WcETS3axZNtxUnNLNAeaiDQaJpv+xqqXgoICzGYzFouFoKAgo+OIuJT5Pxzkhe9rz3fm5+VO6xBfOkUG0a1VEF2jzXSNDiLYz8uAlCLSHNX391uFqJ5UiETO7ERBKVuP5dkfu9MtVFTV/ddL6xBfukWb7SWpd5sQzL6e5zmxiDQHKkROpkIk4pjSiirS809x7GQJezIK2HXcwq50C6m5tQdf+3m5M+iCFrRrGUCnyEB6xgQT18IfNzeTAclFpClRIXIyFSIR57CUVLA7w8Lu4wXsSreQlJrPsZMltbYL9PGgR2szPVsH0zMmmF4xwUQE+RiQWEQaMxUiJ1MhEmkYNpuNLcfy2JdRwOHsYnYdt7DzuIWySmutbSOCvGsUpO6tzQT56FSbiJyeCpGTqRCJnD8VVVYOnChkR5qF7an5JKXmc+BEIXVd1d8uzJ9urcz0aG1mZLdIWgX7YjLpVJuIVFMhcjIVIhFjlZRXsut4ATvSqgvS9rT8OscjBXh70CkykMs6h9MpMpAOEYEqSSLNmAqRk6kQibiek0Vl7DhuYfdxCyv2ZrE9Nb/O7cICvBjcPoy2LfzpEh1Ez9bBRJo1HkmkOVAhcjIVIhHXV1pRRWpuCWsO5lSfZsss5EBWIXX9LdcxIpBOUYFcGB5gH5MUqPFIIk1OfX+/Pc5jJhGRBuXj6c6FEYFcGBFoX1ZcVsmaA9kcPVnC4ewidqcXsD+zgP0nCtl/otC+nckEHcID6d0mmPjYEHrHhnBBS3+dahNpJnSEqJ50hEik6ThZVMaWY3kcyS5mT0YB21LySMurPR7J7OtJfGwwvf9TkHrGmHUUSaSR0SkzJ1MhEmnasgpKSUzJZ1tKHokpeexIq33pv8lUfaotPjaELtFBxPxnWhKNRxJxXSpETqZCJNK8lFda2ZtRQGJKHokp+SQey+N4fu2jSADhgd70aB1Mj9Zm+80kQ/w1X5uIK2gUhWjNmjX87W9/Y+vWrWRkZLB48WKuvvpq+/pFixbx5ptvsnXrVk6ePMm2bdvo1atXjX1ceumlrF69usayG264gU8++cT+PC8vj+nTp7N06VIAxo4dy6uvvkpwcHC9s6oQiUj1UaQ8tqXkczi7mJTcYg5nF1NVxw2SYkJ96dE6mJ6tzXRvVX0TyQBvDdsUOd8axaDq4uJievbsya233sp1111X5/rBgwczfvx4pk2bdtr9TJs2jaefftr+3NfXt8b6SZMmkZaWxvLlywG44447mDJlCl999ZWTPomINAfhQT6M7BbFyG5R9mWnyqvYnW5he5qFHWn57EizkJxTTGruKVJzT/HNjgyg+nTbBS0D7EeQerQ20zkqCB9Pd6M+joj8D0ML0ahRoxg1atRp10+ZMgWAo0ePnnE/fn5+REZG1rlu7969LF++nA0bNtC/f38A3nrrLQYOHMj+/fvp2LHj2YUXEQF8vdxJaBtKQttQ+zLLqQp2HbewPS2fHanVRSndUsqhrCIOZRWxKPE4AB5uJjpGBtqPJPVoHUyHiAA83N2M+jgizVaTOH774Ycf8u9//5uIiAhGjRrFnDlzCAysvuz2l19+wWw228sQwIABAzCbzaxfv/60haisrIyysjL784KCgob9ECLSZJh9PRncPozB7cPsy7ILy9iRls/2NAs7/3Mk6WRxObvTC9idXsDHm6q38/Zwo2t0ED1aB5PQNoSRXSNVkETOg0ZfiG666Sbi4uKIjIxk165dzJ49m+3bt7NixQoAMjMzCQ8Pr/W68PBwMjMzT7vfefPm8dRTTzVYbhFpXloGenN55wgu7xwBVE9qezz/VPV8bf85krTruIXCssrqQdwp+by3/iiB3h7Etwmhe6sgurcy062VWVORiDSARl+I/ndsUbdu3bjwwgtJSEggMTGR3r17A9T5F4fNZjvjXyizZ8/mgQcesD8vKCggJibGiclFpDkzmUy0DvGjdYgfV3avHpNktdo4klPMzuP5bD2Wx5Jt6RT+58aSaw5k218bHuhNp6gg2oX5c2FEAMM6RxARpEv/Rc5Foy9Ev9W7d288PT05ePAgvXv3JjIykhMnTtTaLjs7m4iIiNPux9vbG29v74aMKiJSg5ubifbhAbQPD+Ca+NY8flUX9mUWsvO4hV1pFnYet3DgRCFZhWVkFf63JD22eBeRQT50jKyejmRAXAt6x4Zg9tNNJEXqq8kVot27d1NRUUFUVPW/uAYOHIjFYmHTpk3069cPgI0bN2KxWBg0aJCRUUVEzsjH051e/5ln7VelFVVsT83nSE4xR7KL2Hqs+j5JmQWlZBaUsvpANm+uPgJUj2Xq1iqIq3u1ol3LANq3DFBJEjkNQwtRUVERhw4dsj9PTk4mKSmJ0NBQYmNjyc3NJSUlhfT0dAD2798PQGRkJJGRkRw+fJgPP/yQK6+8krCwMPbs2cPMmTOJj49n8ODBAHTu3JmRI0cybdo03nzzTaD6svvRo0frCjMRaXR8PN3p364F/du1sC8rKqtkf2Yh+zIL2HXcwur92aRbSrGcquDnQyf5+dBJADzdTXSJNhPXwo9LOrSkU2QQ7Vr669J/EQy+MeNPP/3E0KFDay2fOnUq7733Hu+99x633nprrfVz5sxh7ty5pKamMnnyZHbt2kVRURExMTFcddVVzJkzh9DQ/14Cm5ubW+vGjPPnz9eNGUWkySopr+RoTgnLd2Xwy5GTHM87RbqltNZ23h5u9I4NoX14AIMuaEG3VmaizD66sk2ajEZxp+rGRIVIRBq7I9lFHDhRyI40C78cOUlyTjH5JRW1tvP1dKd7azNDOrSkU2QgcWH+tGsZYEBikXOnQuRkKkQi0tTYbDb2ZRayO72AnWn5rDuUQ2ruKcqrrLW2bdfSnwtaBtA5Koj42GDiY4IJ9tN8beL6VIicTIVIRJqDKquN5Jwi1h7MYcORk6TlnWJvRgF1TNdGu5b+xMeEVBek2GA6RQbh7qb7I4lrUSFyMhUiEWmusgpK2ZNRwNGcYnakWdiWmk9yTnGt7cICvOgcFURsqB89W1eXpAtaBuCmkiQGUiFyMhUiEZH/yi0uZ3tqPttS8tiWmk9SSj6FZZW1tgv08aBXTDD940LpEBFIu5YBtAvzV0mS80aFyMlUiERETq+iysrWY3mk5pZwKLuIpJTq+dpOVVTV2jYyyIeeMWYuaBlAt1ZmesYEE2320XQk0iBUiJxMhUhExDGVVVb2nyhk67E8Nh/NIzmniMNZxXWWpLAAL3q2DqZH62B6xpjp2TqYEH8N2pZzp0LkZCpEIiLnrrSiii1H8ziUVcj+E0XsPJ7PvoxCKusYtd2mhR99YkPo3SaEPm1C6BARqEHb4jAVIidTIRIRaRilFVXsyShge2o+21OrT7UdqWPQdqC3B71ig+nzn4LUKyaYQB9NRSJnpkLkZCpEIiLnj6WkgqS0fLYezWVrSh7bUvIpKa95qs1kgmizL7GhfvRpE0JC2xDatPAnLszfoNTiilSInEyFSETEOL+OR0o8lseWY3lsPZZHWt6pOrftEhVEt1ZBXBgeyIURAfSPa4Gvl+Zra65UiJxMhUhExLXkFJVxNKeYpdvTSc0tIcNSyr7Mwlrb+Xq60ykqkE6RQXSJCqR762B6tDLr0v9mQoXIyVSIRERcX1peCUmp+Rw8UWS//P94fu0jSWZfTzpHBdIlykznqEB6twmhXZi/Lv1vgur7++1xHjOJiIg0qNYhfrQO8bM/t1ptHM4uYm9mIfsyCtibUcDmo3lYTlWw4UguG47k2rcNC/AmoU0I3VoFMfCCMC5o6a/52poRHSGqJx0hEhFpGsorrRw4Ucie/xSk3ccLSErLp7yy9qS2Zl9PEtqE0DOmehqSdi39uTA8AA93NwOSy9nQKTMnUyESEWm6Siuq7Jf8rz2Uw/7MAk4UlNW5bai/F12igogL86dzVBBdooPoGBGogdsuSoXIyVSIRESal5LySg5nFfPz4RwOZRVVn3rLKKC0ovaRJDcTxIX50yXazMUXhtE5MoiYUF+dcnMBKkROpkIkIiKnyqvYlprH8bxTHMoqsp92yykqr3P7+Njg6qIUFUR8bAhdo4Pw8dSRpPNJhcjJVIhEROR0sgpL2ZNewKbkXH4+fJL0/FNkF9Y+5ebpbqJLtJkercy0DvElJtSPAe1aEKp52xqMCpGTqRCJiIgjjp0sZnuahWM51f83KTXvtEeSWof40ikyiM7/uV/SgHahtAjwPs+JmyYVIidTIRIRkXNhs9lIyztFYkoe+zILycg/xb7MwjpvJgnVV7i1aeHHpR1aEtfSn9hQP2JC/WgZ4K37JTlAhcjJVIhERKQh5BWXsy+zkP2ZBezLLCQxJY8DJ4pOu72PpxsxIdXlKDbUj85RgcTHhhAW4E2In6fK0m+oEDmZCpGIiJwvBaUVZOSXsin5JDuPW0jNPUVKbgkZllNYz/CrHWX2ISLIh+hgH9q3DOCC8ADahwfQLiyg2d4WQHeqFhERaaSCfDwJivSkY2RgjeXllVYyLNXl6NfHhsMnOZZbQn5JBRmWUjIspSSl1t5nZJAPbVr4ERfmT5sW/rRt4UfbMH/atPDDz0t1QN+AiIhII+Hl4UabFtWF5rcKSis4eKKQ7MJy0vJKOJRVVP3ILiK/pILMglIyC0rZmJxb67Xhgd60DasuSW1a+P+nNFX/OcC7eVSF5vEpRUREmrggH0/6tAmtc11+STnJOcUcO1nC0ZPFHM0p5ujJEo6dLCavpIKswjKyCsvYVEdZCgvwJi7Mz35UqX14AL3bhBDm742bW9MZr6QxRPWkMUQiItIUWUoqqkvSyWKO5lSXpOrnJeQW132bAKi+O3fHyCAuDA+gdYgvrUJ8iQmpPiUXHeyLu4uUJQ2qdjIVIhERaW4spypIOVlC8slijuUUk3yymI1Hcjmef+qMr/PycKPtf8YrxYUF0C7Mn3Ytq0/Fhfp7ndcr4VSInEyFSEREpFpFlZWswjK2p+aTmlvC8fxTpOWd4tjJYlJyS6ioOn21MPt60j48gA4RAYT6exEXFkCXqCBC/b1oEeCFp7ubU7PqKjMRERFpEJ7ubrQK9qVVsG+tdVVWG8fzTnEkp4ijOcUcySkmOaeYI9nFHM8/heVUBVuP5bH1WF6t1y6Y0ofhXSPPx0eoRYVIREREnMbdzURsCz9iW/hBx5rrTpVXcfRkMTvS8smwlJKcU8zejAJyiyvIKyk3dE43FSIRERE5L3y93OkcFUTnqNqnrqxnuuPkeaBCJCIiIoYz+hJ+545cEhEREWmEVIhERESk2VMhEhERkWZPhUhERESaPRUiERERafZUiERERKTZUyESERGRZk+FSERERJo9FSIRERFp9lSIREREpNlTIRIREZFmT4VIREREmj0VIhEREWn2NNt9PdlsNgAKCgoMTiIiIiL19evv9q+/46ejQlRPhYWFAMTExBicRERERBxVWFiI2Ww+7XqT7fcqkwBgtVpJT08nMDAQk8nktP0WFBQQExNDamoqQUFBTttvU6bvzDH6vhyn78wx+r4co+/LcefyndlsNgoLC4mOjsbN7fQjhXSEqJ7c3Nxo3bp1g+0/KChI/8NwkL4zx+j7cpy+M8fo+3KMvi/Hne13dqYjQ7/SoGoRERFp9lSIREREpNlTITKYt7c3c+bMwdvb2+gojYa+M8fo+3KcvjPH6PtyjL4vx52P70yDqkVERKTZ0xEiERERafZUiERERKTZUyESERGRZk+FSERERJo9FSKDvfbaa8TFxeHj40OfPn1Yu3at0ZFc0rx58+jbty+BgYGEh4dz9dVXs3//fqNjNRrz5s3DZDIxY8YMo6O4tOPHjzN58mRatGiBn58fvXr1YuvWrUbHclmVlZU8/vjjxMXF4evrS7t27Xj66aexWq1GR3MJa9asYcyYMURHR2MymViyZEmN9Tabjblz5xIdHY2vry+XXnopu3fvNiasizjTd1ZRUcEjjzxC9+7d8ff3Jzo6mptvvpn09HSnvLcKkYEWLlzIjBkzeOyxx9i2bRsXX3wxo0aNIiUlxehoLmf16tXcc889bNiwgRUrVlBZWcnw4cMpLi42OprL27x5MwsWLKBHjx5GR3FpeXl5DB48GE9PT7799lv27NnDiy++SHBwsNHRXNZzzz3HG2+8wfz589m7dy/PP/88f/vb33j11VeNjuYSiouL6dmzJ/Pnz69z/fPPP89LL73E/Pnz2bx5M5GRkVxxxRX2uTObozN9ZyUlJSQmJvLEE0+QmJjIokWLOHDgAGPHjnXOm9vEMP369bPdeeedNZZ16tTJNmvWLIMSNR5ZWVk2wLZ69Wqjo7i0wsJC24UXXmhbsWKFbciQIbb77rvP6Egu65FHHrFddNFFRsdoVK666irbbbfdVmPZtddea5s8ebJBiVwXYFu8eLH9udVqtUVGRtr++te/2peVlpbazGaz7Y033jAgoev57XdWl02bNtkA27Fjx875/XSEyCDl5eVs3bqV4cOH11g+fPhw1q9fb1CqxsNisQAQGhpqcBLXds8993DVVVcxbNgwo6O4vKVLl5KQkMD48eMJDw8nPj6et956y+hYLu2iiy5i1apVHDhwAIDt27ezbt06rrzySoOTub7k5GQyMzNr/AZ4e3szZMgQ/QY4wGKxYDKZnHIkV5O7GiQnJ4eqqioiIiJqLI+IiCAzM9OgVI2DzWbjgQce4KKLLqJbt25Gx3FZn3zyCYmJiWzevNnoKI3CkSNHeP3113nggQd49NFH2bRpE9OnT8fb25ubb77Z6Hgu6ZFHHsFisdCpUyfc3d2pqqrimWee4cYbbzQ6msv79e/5un4Djh07ZkSkRqe0tJRZs2YxadIkp0ySq0JkMJPJVOO5zWartUxq+tOf/sSOHTtYt26d0VFcVmpqKvfddx/ff/89Pj4+RsdpFKxWKwkJCTz77LMAxMfHs3v3bl5//XUVotNYuHAh//73v/noo4/o2rUrSUlJzJgxg+joaKZOnWp0vEZBvwFnp6KigokTJ2K1Wnnttdecsk8VIoOEhYXh7u5e62hQVlZWrX8xyH/de++9LF26lDVr1tC6dWuj47isrVu3kpWVRZ8+fezLqqqqWLNmDfPnz6esrAx3d3cDE7qeqKgounTpUmNZ586d+eKLLwxK5PoeeughZs2axcSJEwHo3r07x44dY968eSpEvyMyMhKoPlIUFRVlX67fgN9XUVHBhAkTSE5O5ocffnDK0SHQVWaG8fLyok+fPqxYsaLG8hUrVjBo0CCDUrkum83Gn/70JxYtWsQPP/xAXFyc0ZFc2uWXX87OnTtJSkqyPxISErjppptISkpSGarD4MGDa93K4cCBA7Rp08agRK6vpKQEN7eaPyPu7u667L4e4uLiiIyMrPEbUF5ezurVq/UbcAa/lqGDBw+ycuVKWrRo4bR96wiRgR544AGmTJlCQkICAwcOZMGCBaSkpHDnnXcaHc3l3HPPPXz00Ud8+eWXBAYG2o+smc1mfH19DU7negIDA2uNr/L396dFixYad3Ua999/P4MGDeLZZ59lwoQJbNq0iQULFrBgwQKjo7msMWPG8MwzzxAbG0vXrl3Ztm0bL730ErfddpvR0VxCUVERhw4dsj9PTk4mKSmJ0NBQYmNjmTFjBs8++ywXXnghF154Ic8++yx+fn5MmjTJwNTGOtN3Fh0dzfXXX09iYiJff/01VVVV9t+C0NBQvLy8zu3Nz/k6NTkn//d//2dr06aNzcvLy9a7d29dRn4aQJ2Pd9991+hojYYuu/99X331la1bt242b29vW6dOnWwLFiwwOpJLKygosN1333222NhYm4+Pj61du3a2xx57zFZWVmZ0NJfw448/1vn31tSpU202W/Wl93PmzLFFRkbavL29bZdccolt586dxoY22Jm+s+Tk5NP+Fvz444/n/N4mm81mO7dKJSIiItK4aQyRiIiINHsqRCIiItLsqRCJiIhIs6dCJCIiIs2eCpGIiIg0eypEIiIi0uypEImIiEizp0IkIiIizZ4KkUgzZTKZWLJkidExnO6WW27h6quvPi/v1bZtW1555ZXz8l7O8Nvv5tJLL2XGjBmG5RFxJSpEIk3ILbfcgslkwmQy4enpSUREBFdccQXvvPNOrQk3MzIyGDVqVL3225jK09///nfee++9c9rH3Llz7d+jm5sb0dHR3HTTTaSmptbYbvPmzdxxxx3n9F7nkzO+G5GmSoVIpIkZOXIkGRkZHD16lG+//ZahQ4dy3333MXr0aCorK+3bRUZG4u3tbWDShmE2mwkODj7n/XTt2pWMjAzS0tJYuHAhO3fuZMKECTW2admyJX5+fuf8Xqdjs9lq/P/sXDnruxFpilSIRJoYb29vIiMjadWqFb179+bRRx/lyy+/5Ntvv61xdOB/j/qUl5fzpz/9iaioKHx8fGjbti3z5s0Dqk8LAVxzzTWYTCb788OHDzNu3DgiIiIICAigb9++rFy5skaWtm3b8uyzz3LbbbcRGBhIbGxsrdnj09LSmDhxIqGhofj7+5OQkMDGjRvt67/66iv69OmDj48P7dq146mnnjpjSajrtND06dN5+OGHCQ0NJTIykrlz5/7u9+jh4UFkZCTR0dFcfPHFTJs2jQ0bNlBQUFDj8/16yuzGG29k4sSJNfZRUVFBWFgY7777LlBdcJ5//nnatWuHr68vPXv25PPPP7dv/9NPP2Eymfjuu+9ISEjA29ubtWvX1sp29OhRTCYTn376KRdffDG+vr707duXAwcOsHnzZhISEggICGDkyJFkZ2ef9rv5rfLych5++GFatWqFv78//fv356effrKvP3bsGGPGjCEkJAR/f3+6du3KsmXLfve7FGkMVIhEmoHLLruMnj17smjRojrX/+Mf/2Dp0qV8+umn7N+/n3//+9/24rN582YA3n33XTIyMuzPi4qKuPLKK1m5ciXbtm1jxIgRjBkzhpSUlBr7fvHFF0lISGDbtm3cfffd3HXXXezbt8++jyFDhpCens7SpUvZvn07Dz/8sP303nfffcfkyZOZPn06e/bs4c033+S9997jmWeecejz/+tf/8Lf35+NGzfy/PPP8/TTT7NixYp6vz4zM5NFixbh7u6Ou7t7ndvcdNNNLF26lKKiIvuy7777juLiYq677joAHn/8cd59911ef/11du/ezf3338/kyZNZvXp1jX09/PDDzJs3j71799KjR4/T5pozZw6PP/44iYmJeHh4cOONN/Lwww/z97//nbVr13L48GGefPLJen/OW2+9lZ9//plPPvmEHTt2MH78eEaOHMnBgwcBuOeeeygrK2PNmjXs3LmT5557joCAgHrvX8Sl2USkyZg6dapt3Lhxda674YYbbJ07d7Y/B2yLFy+22Ww227333mu77LLLbFartc7X/u+2Z9KlSxfbq6++an/epk0b2+TJk+3PrVarLTw83Pb666/bbDab7c0337QFBgbaTp48Wef+Lr74Ytuzzz5bY9kHH3xgi4qKOm2G334HQ4YMsV100UU1tunbt6/tkUceOe0+5syZY3Nzc7P5+/vbfH19bYANsE2fPr3Gdm3atLG9/PLLNpvNZisvL7eFhYXZ3n//ffv6G2+80TZ+/HibzWazFRUV2Xx8fGzr16+vsY/bb7/dduONN9psNpvtxx9/tAG2JUuWnDabzWazJScn2wDbP//5T/uyjz/+2AbYVq1aZV82b948W8eOHe3P6/pu7rvvPpvNZrMdOnTIZjKZbMePH6/xXpdffrlt9uzZNpvNZuvevbtt7ty5Z8wm0lh5GNjFROQ8stlsmEymOtfdcsstXHHFFXTs2JGRI0cyevRohg8ffsb9FRcX89RTT/H111+Tnp5OZWUlp06dqnWE6H+PcJhMJiIjI8nKygIgKSmJ+Ph4QkND63yPrVu3snnz5hpHhKqqqigtLaWkpKTe43d+e5QlKirKnuF0OnbsyNKlSykrK+PLL7/ks88+O+ORKU9PT8aPH8+HH37IlClTKC4u5ssvv+Sjjz4CYM+ePZSWlnLFFVfUeF15eTnx8fE1liUkJDj8uSIiIgDo3r17jWW/9zl/lZiYiM1mo0OHDjWWl5WV0aJFCwCmT5/OXXfdxffff8+wYcO47rrrzngES6QxUSESaSb27t1LXFxcnet69+5NcnIy3377LStXrmTChAkMGzasxviW33rooYf47rvveOGFF2jfvj2+vr5cf/31lJeX19jO09OzxnOTyWQ/Jebr63vGzFarlaeeeoprr7221jofH58zvra+GU7Hy8uL9u3bA9UDrA8ePMhdd93FBx98cNrX3HTTTQwZMoSsrCxWrFiBj4+P/Uq+X9/vm2++oVWrVjVe99vB7f7+/g5/rl/L7m+X/d7n/JXVasXd3Z2tW7fWOi3462mxP/zhD4wYMYJvvvmG77//nnnz5vHiiy9y77331us9RFyZCpFIM/DDDz+wc+dO7r///tNuExQUxA033MANN9zA9ddfz8iRI8nNzSU0NBRPT0+qqqpqbL927VpuueUWrrnmGqB6PNDRo0cdytWjRw/++c9/2t/nt3r37s3+/fvtxcRITzzxBB06dOD++++nd+/edW4zaNAgYmJiWLhwId9++y3jx4/Hy8sLgC5duuDt7U1KSgpDhgw5n9HrJT4+nqqqKrKysrj44otPu11MTAx33nknd955J7Nnz+att95SIZImQYVIpIkpKysjMzOTqqoqTpw4wfLly5k3bx6jR4/m5ptvrvM1L7/8MlFRUfTq1Qs3Nzc+++wzIiMj7Zdot23bllWrVjF48GC8vb0JCQmhffv2LFq0iDFjxmAymXjiiSfqfTTiVzfeeCPPPvssV199NfPmzSMqKopt27YRHR3NwIEDefLJJxk9ejQxMTGMHz8eNzc3duzYwc6dO/nLX/5yrl+VQ9q1a8e4ceN48skn+frrr+vcxmQyMWnSJN544w0OHDjAjz/+aF8XGBjIgw8+yP3334/VauWiiy6ioKCA9evXExAQwNSpU8/XR6lThw4duOmmm7j55pt58cUXiY+PJycnhx9++IHu3btz5ZVXMmPGDEaNGkWHDh3Iy8vjhx9+oHPnzobmFnEWXWUm0sQsX76cqKgo2rZty8iRI/nxxx/5xz/+wZdffnnaK6QCAgJ47rnnSEhIoG/fvhw9epRly5bh5lb9V8SLL77IihUriImJsY93efnllwkJCWHQoEGMGTOGESNGnPbIyel4eXnx/fffEx4ezpVXXkn37t3561//as85YsQIvv76a1asWEHfvn0ZMGAAL730Em3atDmHb+jszZw5k2+++abGbQF+66abbmLPnj20atWKwYMH11j35z//mSeffJJ58+bRuXNnRowYwVdffXXaU5nn27vvvsvNN9/MzJkz6dixI2PHjmXjxo3ExMQA1eO37rnnHjp37szIkSPp2LEjr732msGpRZzDZLPZbEaHEBERETGSjhCJiIhIs6dCJCIiIs2eCpGIiIg0eypEIiIi0uypEImIiEizp0IkIiIizZ4KkYiIiDR7KkQiIiLS7KkQiYiISLOnQiQiIiLNngqRiIiINHv/DzZ3qSO8LuPNAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "reach_data = pd.DataFrame(m.sfr.reach_data)\n",
    "reach_data.index = reach_data['reachID']\n",
    "\n",
    "reach_routing = make_graph(reach_data['reachID'], reach_data['outreach'], one_to_many=False)\n",
    "path = find_path(reach_routing, 10)\n",
    "\n",
    "# leave off the outlet segment (0; which doesn't exist) when getting the elevations\n",
    "path_reach_data = reach_data.loc[path[:-1]]\n",
    "elevations = path_reach_data['strtop']\n",
    "distances_mi = np.cumsum(path_reach_data['rchlen'])/5280  # distance along path in miles\n",
    "\n",
    "plt.plot(distances_mi, elevations)\n",
    "plt.ylabel('Streambed elevation, in feet')\n",
    "plt.xlabel('Distance in River miles')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a321216",
   "metadata": {},
   "source": [
    "### Getting the upstream contributing segments\n",
    "\n",
    "The ``get_upsegs`` method takes a reverse routing dictionary (one-to-many) and returns a set of all of the nodes that are upstream of that point. This function is integral to the [streambed elevation smoothing](Streambed_elevation_demo.ipynb) in SFRmaker, and can also be used to aggregate flow components (for example, net groundwater inflow) at a point on the stream network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ee910bc5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-04-03T00:33:33.720595Z",
     "iopub.status.busy": "2025-04-03T00:33:33.720426Z",
     "iopub.status.idle": "2025-04-03T00:33:33.724063Z",
     "shell.execute_reply": "2025-04-03T00:33:33.723603Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{np.int64(1),\n",
       " np.int64(2),\n",
       " np.int64(3),\n",
       " np.int64(4),\n",
       " np.int64(5),\n",
       " np.int64(6),\n",
       " np.int64(9),\n",
       " np.int64(10),\n",
       " np.int64(13),\n",
       " np.int64(14),\n",
       " np.int64(25),\n",
       " np.int64(26)}"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "get_upsegs(routing_r, 42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db61101f",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}