{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c46346c2",
   "metadata": {},
   "source": [
    "# `sue-msa` — Logit stochastic user equilibrium by MSA on the two-route network\n",
    "\n",
    "**What.** Stochastic user equilibrium (SUE) replaces Wardrop's *every used route is\n",
    "shortest* with *every traveller minimises PERCEIVED cost*: route choice is a logit\n",
    "model over the true costs, so flow spreads across routes by a dispersion parameter\n",
    "θ. `sue-msa` reaches the logit-SUE fixed point with the predetermined `1/k`\n",
    "successive-averages step over Dial's efficient-path STOCH loading — no path\n",
    "enumeration (`[sheffi1985urban]`, [docs/REFERENCES.md](../../docs/REFERENCES.md)).\n",
    "\n",
    "**Why it is in the benchmark.** It is the deterministic-UE ladder's stochastic\n",
    "counterpart and the anchor for the SUE task (ADR-001). See the\n",
    "[model compendium](../../docs/MODELS.md) and\n",
    "[docs/ARCHITECTURE.md](../../docs/ARCHITECTURE.md) (P1).\n",
    "\n",
    "**Scope.** Runs on the built-in two-route logit-SUE anchor (θ = 0.5, demand 4) and\n",
    "certifies the SUE fixed-point residual. It does not treat probit SUE — see\n",
    "[`sue-probit-msa`](11-sue-probit-msa.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f14a13a5",
   "metadata": {},
   "source": [
    "## How this notebook is graded\n",
    "\n",
    "**A notebook never claims a number it does not compute in that cell.** Every scored\n",
    "quantity below is recomputed live by the P1 `Evaluator` from the flows the model\n",
    "emitted, in the cell where it is claimed. Model self-reports are shown only as\n",
    "provenance and diffed against the certificate as an honesty check, exactly as the\n",
    "harness treats them ([README](../../README.md), *Certified, not self-reported*)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8749cfc1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T13:45:22.895025Z",
     "iopub.status.busy": "2026-07-21T13:45:22.894848Z",
     "iopub.status.idle": "2026-07-21T13:45:24.888356Z",
     "shell.execute_reply": "2026-07-21T13:45:24.887165Z"
    }
   },
   "outputs": [],
   "source": [
    "# Setup. `sue-msa` is a core model: a plain `pip install -e .` suffices — no\n",
    "# optional extra, so no guard cell. The inline backend is Agg-based (headless CI\n",
    "# renders into the notebook); NEVER matplotlib.use(\"Agg\") in-kernel — it silently\n",
    "# suppresses inline figure capture.\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "\n",
    "from tabench import (\n",
    "    Budget,\n",
    "    DialSUEModel,\n",
    "    Evaluator,\n",
    "    RngBundle,\n",
    "    Trace,\n",
    "    two_route_scenario,\n",
    "    viz,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a967f41c",
   "metadata": {},
   "source": [
    "## The scenario\n",
    "\n",
    "Two disjoint 2-link routes — A = 1→3→2 (cost `2 + f_A`) and B = 1→4→2 (cost\n",
    "`1.5 + 2 f_B`) — with demand 4 and logit dispersion θ = 0.5, the analytic logit-SUE\n",
    "anchor. Scenarios are frozen and content-hashed (P2); the hash is the instance's\n",
    "identity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "498a0162",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T13:45:24.894219Z",
     "iopub.status.busy": "2026-07-21T13:45:24.893562Z",
     "iopub.status.idle": "2026-07-21T13:45:24.899822Z",
     "shell.execute_reply": "2026-07-21T13:45:24.899095Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "scenario      : tworoute\n",
      "content hash  : 9b98e58b339b702f…\n",
      "routes        : 1->3->2 (A) and 1->4->2 (B), disjoint\n",
      "total demand  : 4.0\n",
      "SUE family    : logit (theta = 0.5)\n"
     ]
    }
   ],
   "source": [
    "scenario = two_route_scenario()\n",
    "net = scenario.network\n",
    "\n",
    "print(f\"scenario      : {scenario.name}\")\n",
    "print(f\"content hash  : {scenario.content_hash()[:16]}…\")\n",
    "print(f\"routes        : 1->3->2 (A) and 1->4->2 (B), disjoint\")\n",
    "print(f\"total demand  : {scenario.demand.total}\")\n",
    "print(f\"SUE family    : {scenario.sue_family} (theta = {scenario.sue_theta})\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29a7e463",
   "metadata": {},
   "source": [
    "## Solve\n",
    "\n",
    "The model contract ([CONTRIBUTING.md](../../CONTRIBUTING.md)): a model receives\n",
    "`(scenario, budget, rng, trace)`, records checkpoints, and respects the budget.\n",
    "Budgets are hardware-free (iterations / shortest-path calls; wall-clock is recorded\n",
    "but never the ranking axis, P7). Whatever the model writes into `self_report` is\n",
    "provenance, not a score."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "91b8a19b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T13:45:24.903640Z",
     "iopub.status.busy": "2026-07-21T13:45:24.903037Z",
     "iopub.status.idle": "2026-07-21T13:45:25.004275Z",
     "shell.execute_reply": "2026-07-21T13:45:25.003505Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "model              : sue-msa\n",
      "budget spent       : 300 iterations, 301 shortest-path calls\n",
      "emitted flows      : [2.299096 2.299096 1.700904 1.700904]\n",
      "self-reported resid: 3.990e-07  (provenance only)\n"
     ]
    }
   ],
   "source": [
    "model = DialSUEModel()\n",
    "bundle = model.solve(scenario, Budget(iterations=300), RngBundle(0), Trace())\n",
    "\n",
    "final = bundle.final\n",
    "print(f\"model              : {model.name}\")\n",
    "print(f\"budget spent       : {final.coords.iterations} iterations, \"\n",
    "      f\"{final.coords.sp_calls} shortest-path calls\")\n",
    "print(f\"emitted flows      : {np.round(final.link_flows, 6)}\")\n",
    "print(f\"self-reported resid: {final.self_report['sue_fixed_point_residual']:.3e}  (provenance only)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd77cd41",
   "metadata": {},
   "source": [
    "## Certify (P1)\n",
    "\n",
    "The scored quantity for an SUE task is the **fixed-point residual** — how far the\n",
    "loaded flows are from reproducing themselves under the logit loading — recomputed by\n",
    "the harness `Evaluator`. The logit loading is closed-form, so the model's self-report\n",
    "should match the certificate exactly. We recompute the analytic fixed point in-cell\n",
    "with `brentq` rather than quoting it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "68ac3cc4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T13:45:25.008728Z",
     "iopub.status.busy": "2026-07-21T13:45:25.008512Z",
     "iopub.status.idle": "2026-07-21T13:45:25.015212Z",
     "shell.execute_reply": "2026-07-21T13:45:25.014531Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "certified SUE residual : 3.990e-07\n",
      "feasible               : 1\n",
      "logit fixed point f_A  : 2.299096  (recomputed via brentq)\n"
     ]
    }
   ],
   "source": [
    "from scipy.optimize import brentq\n",
    "\n",
    "evaluator = Evaluator(scenario)\n",
    "metrics = evaluator.evaluate(final.link_flows)\n",
    "residual = metrics[\"sue_fixed_point_residual\"]\n",
    "print(f\"certified SUE residual : {residual:.3e}\")\n",
    "print(f\"feasible               : {metrics['feasible']:.0f}\")\n",
    "\n",
    "assert metrics[\"feasible\"] == 1.0\n",
    "assert residual < 1e-4\n",
    "\n",
    "# Honesty diff (P1): the logit loading is closed-form, so self-report == certificate.\n",
    "assert np.isclose(final.self_report[\"sue_fixed_point_residual\"], residual, rtol=1e-6, atol=1e-12)\n",
    "\n",
    "# Analytic anchor RECOMPUTED in-cell: the binary-logit fixed point via brentq.\n",
    "# Route costs c_A = 2 + f_A, c_B = 1.5 + 2 f_B; f_A = D / (1 + exp(theta (c_A - c_B))).\n",
    "demand = scenario.demand.total\n",
    "theta = scenario.sue_theta\n",
    "\n",
    "def _logit_residual(f_a):\n",
    "    c_a = 2.0 + f_a\n",
    "    c_b = 1.5 + 2.0 * (demand - f_a)\n",
    "    return f_a - demand / (1.0 + np.exp(theta * (c_a - c_b)))\n",
    "\n",
    "f_a = brentq(_logit_residual, 1e-9, demand - 1e-9)\n",
    "ref_flows = np.array([f_a, f_a, demand - f_a, demand - f_a])\n",
    "print(f\"logit fixed point f_A  : {f_a:.6f}  (recomputed via brentq)\")\n",
    "assert np.allclose(final.link_flows, ref_flows, atol=1e-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0aa81d9a",
   "metadata": {},
   "source": [
    "## Visualize\n",
    "\n",
    "Both figures come from `tabench.viz`. Left/top: the SUE link flows on the two-route\n",
    "network — flow spreads across BOTH routes (the stochastic spreading θ controls).\n",
    "Right/bottom: the emitted flows against the `brentq`-recomputed logit fixed point;\n",
    "on-diagonal means the MSA loading reached the analytic SUE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "646e8bda",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T13:45:25.019944Z",
     "iopub.status.busy": "2026-07-21T13:45:25.019484Z",
     "iopub.status.idle": "2026-07-21T13:45:25.279326Z",
     "shell.execute_reply": "2026-07-21T13:45:25.278546Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAesAAAG0CAYAAAASKw+DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA/O0lEQVR4nO3dd3wU1d4G8Gc2uwspJCEhJBBAkkASCIgUqQlIkCK9BmkWlCaC6LXkIiV0fC+gAgJSvCoXCyIKAoIQikGqxAIqAZLQQq+hhGzJvn9gVkLaZGfLzM7z5bMf3c2ZMydtn/zOOTsrmExGC4iIiEi2NK4eABEREZWMYU1ERCRzDGsiIiKZY1gTERHJHMOaiIhI5hjWREREMsewJiIikjmGNRERkcwxrImIiGSOYU3kILNmzUJsbKz1/qhRozBw4MAy9eHn548NGzaU6Zj//vdj1K0bA3//ili0aFGhcRCR8jCsSbQuXbogMTHR1cOwK2cG2ezZs7Fo0SKHniM7OxtvvPEGxo17BUeP/oXnnnvOoecjIufQunoARGVlsVhgNpuh1Srrx9fPz8/h5zh79iyMRiM6dOiAkJAQh5+PiJyDlTWJMmrUKOze/RMWL14CPz9/+Pn549SpU2jT5gnMn7/A2m7gwIEIDKyE27dvAwCysrLg5+eP9PQMAMD16zcwYsQI1KjxCEJCqqBPn75IT08v8dwpKSnw8/PH1q1b0bp1GwQFVcbevXuRm5uLN998ExERtVC5cjA6duyEQ4dSrcetWrUKNWrUKNDXhg0b4Ofnb/347Nnv4PDhI9bPadWqVQCAGzdu4OWXxyA8PALVqlVH167dcPjwYclfwwenwbt06YI333wTEydOwiOP1ETt2pGYNWtWiX3MnDkTkZFROHLkSKGPrVq1Ci1atAQANGjwmPV79LC8vDy88847qFOnLoKCKiM2Nhbbtm2zfnzIkGfw+utvWO8nJibCz88fx44dAwAYDAZUqVIVO3bsBAB8++06tGjREsHBIahZMwzdu/fAnTt3xH9hiKhUDGsSZfbs2WjatCmeffZZHDuWhmPH0lCtWjXExrbC7t27AdyvePfs2Qs/Pz/s27cPAPDTTz+hatWqiIgIBwC89NIo/PLLr/jii8+xdesPsFgs6Nu3H4xGY6ljSEpKQlJSEg4cOICYmHqYNGkS1q//DkuWLMaPP+5CeHgYevfujWvXrov6nHr37o2XX34ZderUsX5OvXv3BgA8++xzuHLlMtasWYNdu3aiQYMG6N69h7XvU6dOwc/PHykpKWX+Wj7o88+/gLe3F7ZvT8bUqVPwzjv/h+3bdxRqZ7FY8MYbb+Dzz7/A999/j3r16hX5+axbtw4AsH37duv36GGLFy/GwoUfYNq0adiz5yfEx7fD008PsP7R9OD3FAB27/4JgYGBSEm5/1hqaiqMRiOaNWuKCxcu4IUXXsDgwYNw4MB+bNy4Ad26dYPFwjfzI7InhjWJ4ufnB71eBy8vTwQHByM4OBgeHh6IjY3Fvn17YTabceTIEej1eiQk9LM+saek7EarVq0AAOnp6di06XssWDAfLVu2RP369bF8+TKcP38eGzZsLHUM48ePR3x8W4SHh6FcOT1WrPgI06ZNRfv27REdHY358+fD09MTK1euFPU5eXp6wsfHG1qth/Vz8vT0xN69e5GamopPPvkEjRo1REREBGbMmA4/Pz9rGOp0OtSuXRteXl42fkXvi4mJQWJiIiIiIjBgwAA0bNgQu3btKtDGZDJj2LDh2LXrR2zZstn6h09Rn09AQEUAQKVKgdbv0cMWLFiIV155BX379kHt2rUxdeoU1K9fH4sWLQYAxMbG4ujRo7hy5QquX7+BtLQ0jBo10hrgKSm70ahRI3h5eeHChQswmUzo1q0bHnnkEcTExGDYsBfh4+Mj6etCRAUpa9GPZKdFixa4des2fvvtdxw4sB+tWrVCbGws3n33PQD3K+uxY8cCANLS0qDVatGkSRPr8QEBAahVqxaOHUsDAPTp0xd79+4FAFSvXh379++ztm3YsKH1/zMzM/+u7ppZH9PpdGjcuJG1L1sdOXIEt2/fRlhYwVDMyclBZmYmAKBq1ar4+eeDks4D3A/rB4WEBOPKlcsFHhs/fjz0ej2Sk7chMDBQ0vmys7Nx/vx5NG/erMDjzZs3w+HD96fW69ati4oVK2L37p+g1+vw6KOPomPHjli2bDmA+9/T/E159evXR5s2bdCyZSvEx8cjPj4ePXr0QMWK/pLGSUQFMaxJEn9/f9SrVw+7d+/GgQMH0LZtW7Rs2QrPPz8UJ06cQHp6OmJjW4nub8GC+cjJuQcA0OkK/niWtYrVaDR4eDbWaDSVetzt23cQEhJS5Eum/P3tu0ns4c9REATk5eUVeKxt2yewZs3XSE5ORkJCgl3PXxRBENCyZUvs3r0b5crpERsbi3r16iE3Nxd//vknDhw4gDFjxgAAPDw8sG7dt9i/fz+2b9+OpUs/xLRp05CcvA01a9Z0+FiJ1ILT4CSaTqeH2Wwu9HhsbCukpKRgz549iIuLRUBARURFRWLOnDkICQlBrVq1AABRUVEwmUz4+eefrcdeu3YNJ06cQFRUNABY17cjIsILbQ57UFhYGPR6Pfbv3299zGg0IjX1F2tflSpVwq1btwpsdnp4k9j9z6lgODZo0AAXL16EVuthHUv+TWpla4unnnoKy5cvw5gxY7FmzdeS+vL19UWVKlWwb9/+Ao/v27cf0dHR1vv569YpKbsRFxcLjUaDVq1aYv78+cjNzS1QmQuCgObNm2P8+PFISUmBXq8v82vDiahkDGsSrUaNGvj550M4deoUrl69aq0AY2NjkZycDK1Wi8jISOtjq1d/ZV2vBoCIiAh06dIZY8e+gr179+Lw4cMYNmw4qlSpgi5dOpdpLN7e3njhhaGYOHEStm3bhqNHj2Ls2LG4e/cuhgwZAgBo3LgJvLy8MHXqVGRkZOKrr77CZ599VqCfRx6pgVOnTuH333/H1atXkZubi7Ztn0DTpk0xaNAgJCdvx6lTp7B//35MnToNqam/AADOnTuHJk0ex6FDh2z8apZNt27d8OGHH2L06NH49tt1kvoaO3YM3n//fXz99VocP34ckycn4fDhwxg1aqS1Tf669dGjR9G8eXPrY6tXf4WGDRvC29sbAPDzzz9jzpy5SE39BWfOnMH69d/hypUriIqKkjRGIiqIYU2ijR07Bh4eHmjWrDnCwyNw5swZAECLFi2Rl5dXIJhjY2NhNpsLXXDkgw8W4bHHGqB///5o374DLBYL1qz5CjqdrszjSUpKQvfu3TB8+Ai0bt0GGRmZWLt2rXW9NCCgIpYu/RA//LAVLVu2xJo1awpd1KV79+5o164dunbthvDwCKxZswaCIOCrr1ajZcuWGD16NBo3boKhQ1/AmTNnULlyEID7Vfzx48dx9+7dMo/bVj179sDixYsxYsQIrF+/3uZ+Ro4cidGjX8KECRPQokVLJCdvwxdffI6IiAhrm5iYGPj5+aF+/frWzWJFfU8rVKiAPXv2oF+/fmjcuAmmT5+OGTOmo3379rZ/okRUiGAyGfkaCyIiIhljZU1ERCRzDGsiIiKZY1gTERHJHMOaiIhI5hjWREREMsewJiIikjmGNRERkcwxrImIiGSOYU1ERCRzDGsiIiKZY1gTERHJHMOaiIhI5hjWREREMsewJiIikjmGNRERkcwxrImIiGSOYU1ERCRzDGsiIiKZY1gTERHJHMOaiIhI5hjWREREMsewJiIikjmGNRERkcwxrImIiGRO68yTvfb5eew6eseZp7RqE+2NeQOquOTcREREUji1st519A6MZoszTwkAMJotLvsjgYiISCqnVtYAoPMQkDIh0qnnjJt+zKnnIyJytuzsbGRkpOPSpUswGY3Q6nSoXLkywsMj4Ovr6+rhFWIwGGAymST3o9Vqodfr7TAieXN6WBMRkX3cvHkTycnbsHPHdly4cKHYdiEhIWgb3w7x8e3g5+fnxBEWzWAwoN+AISinkR7W/v7+WPjBYrcPbIY1EZHCGAwGrF79JTZu2ACz2YSoqGi0bdsO4RERCA2tBr1OB4PRiKyss8hIT0dq6iF8/tkqrP7yS3Tp2hUJCf1dGm4mkwnlNCbsvBYJk8X21VitkIcncAwmk4lhTURE8pGRkYEF899DVlYWYuNao1fP3qgZFlZk2+DgYDRq1Bh9+yXgZGYmvvl2Ldav+xaHfj6IMWPHITw83MmjL8hk0cBs8XDpGJSCL90iIlKIP//4A0mTJ+LOnbtI/PfbePXVfxUb1A+rGRaGV1/9FxITx+POnbtImjwRf/7xh4NHXArBDjeVYFgTESlARkYGZs2aiYoVAzD7nf/g8ceb2tTP402bYfY7/7nfz+yZyMjIsPNIxRPs8E8tGNZERDJnMBiwYP578PT0xOSkqQgKCpLUX1BQECYnTUX58p5YMP89GAwGO420rARAkHBjWBMRkVysXv0lsrKyMHLUS5KDOl9QUBBGjhyFrKwsrF79pV36JMdhWBMRydjNmzexccMGxMa1tnnquziPN22GVrFx2LRxA27evGnXvsXhorVYsg/re/fu4bkhA9CqaUO0a9MS/fv0QGZGuquH5TQWi/Ov+EZE8pGcvA1mswm9evYWfczKlSvh4+OD7777rtS2vXr1hslkwvbtyVKGaRtmtWiyD2sAGPLMc9i9PxXJu/ag41Od8a9xY1w9JIe7k5uHzMtGHMkyItfIwCZSq507tiMqKlr0ru9Tp07h448/RtOm4qrwsLBwREZFYYcrwppEk31Yly9fHu3ad4Qg3P8TqlHjx3HmzGkXj8px8kM647IJt3Pvh7TBxLAmUqPs7GxcuHABjRo1FtU+Ly8Po0ePxpw5c8p0kZBGDRvjwoULuHXrlq1DtQl3g4sn+7B+2PKli9Hxqc6uHobdFRXSRKRuGX8v+YVHRIhqv2DBAjRv3hwNGzYs03ny+89w9hKjlJ3g1h3h6qCoK5i9/+4cnMzMwOq5pa/DKEnaBQMMJVwiN+uGCR7Z6vmhJHKmCuUFhPjJ86nw0qVLAIDQ0Gqltv3jjz+wbt06bNmypcznye//4sWLZT5WGqkLz+p5XpTnT2gRFi+cj00bvsPqtevg5eXl6uHYVymFtMUC5HGjGZFDWCzyfcI3GY0AAL1OV2rbPXv24NSpU2jQoAGA+8E7ZswYXLhwAcOGDSvxWL1eV+B8JD+KCOslixbim7VrsHrtOvj5+bt6OHYXVUWP2/fycDHbjLuGwqFcraIWFTwVt2JBRBJp/w5pg4gQHTZsWIFQ7tSpE0aPHo1u3bqVeqzBYCxwPqeR799JsiP7sD53LgtTJo3HIzVrom/PrgAAvV6PTT/scPHI7MunvAY+5TUlhjYRqUvlypUBAFlZZxEcHOyw82RlnQUAh56jKFI3ialpg5nsw7pq1VCcv5Lt6mE4TVGhrdOq5weSiP4RHv73xq/0dNE7wvNt3rxZdNuM9PQC5yP54dyqTPmU1yCisg71QnUor2NYE6mRr68vQkJCkJp6yKHnSf3lEEJCQlChQgWHnqcQ7gYXjWEtc4KKfhiJqLAn2sYjLe0oTmZmOqT/zMwMHEtLQ9v4dg7pv2S8hJlYDGsiIhlr1+5JeHho8c23ax3S/zffrIVWq0W8S8KaxGJYExHJmJ+fH7p07YrdKT/i4MEDdu374IH9+Gl3Cjp36Qo/Pz+79k32xbAmIpK5hIT+CA0NxZLFi3D58mW79Hn58mUsWbIIoaHVkJDQ3y59lpUgCJJvasGwJiKSOb1ejzFjx+HevRxMSZokObAvX76MpMkTce/ePYwZ+0qZriNuX1yzFothTUSkAOHh4UhMHI/r168h8a03bJ4SP3hgPxLfeh03blxHYuJ4hIeH23mk5AgMayIihagbE4OkKdPg7e2F2bNmYN68OaJ3iWdmZmDevDmYPXsmvL29kTRlGurGxDh4xKUQIPGlW64dvjPJ/qIoRET0j/DwcLzzf3OwevWX2LRxA37anYLIqCg0atgY4RERCA2tBr1eB4PBiKyss8hIT0fqL4dwLC0NWq0W3Xv0REJCfxdOfZMtGNZERAqj1+sxePAQdOvWHdu3J2PH9mR88cVnxbavUqUKBg4chLbx7WS161sQeLlRsZwe1kazBXHTjzn9nDoP9XxTiUgd/Pz80KtXb/Tq1RvZ2dnIzMzAxYsXYTIaodXpEBwcjPDwCOdfmYzszqlh3SbaG7uO3nHmKQEAOg8BbaK9nX5ech+DBg1CSsputGnTBitXflrgY7du3cJTTz1lvX/q1Gn8+9+JeOmll5w9TFIxX19fNGjwmKuHUUZ8P2uxnBrW8wZUcebpiOxm5MiRGDx4MD777PNCH6tQoQJ2794NALBYLKhfvz46d+7s7CESKY96slYy7gYnEiEuLg4+Pj6ltjtw4AAqVw5GzZo1HT8oIsXj66zFYlgT2dE333yD3r17uXoYRORmuBucyE4sFgvWr/8OW7f+4OqhECkCd4OLx7AmspO9e/eievXqCA0NdfVQiBSCG8zE4jQ4kZ1wCpyIHIVhTSRC9+498Oyzz2Hr1q2oU6cuDhw4gL59++H8+fMAgLy8PGzYsBE9evRw8UiJFIT7y0TjNDiRCOvXryv02Jo1X1n/X6PR4K+//nTmkIjcAKfBxWJYExGRSwjgBjOxOA1OREQkc6ysiYjINQROg4vFypqIiEjmGNZEREQyx2lwIiJyDU6Di8awJiIil+BucPEY1kRE5BqsrEXjmjUREZHMMayJiIhkjtPgRETkGpwGF42VNRERkcyxsiYiIpfgbnDxGNZEROQanAYXjdPgREREMsewJiIikjlOgxMRiWAw5eHq7TwE+Wqg1bDOsQtOg4vGsCYiKobRbMHNu3m4mWPGXcM/j1fxZ1jbB8NaLIa1g32/5yjW7jiM0CA/9HqiHhpEVnX1kIioBP8EdB7uGiyFPh7ow6BWIqPRiBUrluPw77/j1q1sBAQEoHuPnoiPb1eo7c2bN/HJx//Fn3/+gZycHAQHhyAhoT+aPP64C0Z+H8PaQa7cuINmz83HlZt3rY+9+3kKalapiBd7NkOP1jGoEeLvugESUSEmswXHLhiRVzijAQBeekCvZVjbiyA4r642m82o6O+PiZMmIzg4GMePH8esmdMRGBiIBg0eK9D23r0c1AwLw6DBQ1CxYkWkph7C+++9i1mz3kG16tUljNh2/KlzkJYvLCwQ1PlOnr+OCYs3o/6AuYgf9SEWrP4Jpy/ccP4AiagQDw1QXld8BPh5eThxNGog2OEmTvny5dH/6QEICQmBIAiIjIxETEw9HP3rr0Jtg4ND0L17DwQGBkKj0aBJk8dRtWpVHDt+TMonKwkrawfYf+Q0Ll67XWq7Q0fP4tDRs5iweDMaR1dDr7b1WHETuZAgCKhZSYuTV4wF1qjz+XmyvrE/6evOOTk5Be7rdDrodLoSjzEYDDhx4jhiY+NK7f/mzZs4ezYLjzzyiKRxSsGwdoD1KX+U+ZiHg/vZro0xuFMjeHjwyYHImTQCoNUIAArOhXvpBeg81LOhSUlGjRxe4H7ffglISOhfbHuLxYIlSxajSpUqaNqsWYl9m4xGvPfuPLRo2RIREbXsMl5bMKwdICI0UNLx+cH9yYZD+G7e8/D21NtpZERUEovFgjPXzMi+Z4FveQGmPIu1wvbz4h/Odmenl24tXrIUnp6e1kdLqqotFguWL1uK8+eyMHHiZGhKeBmeyWjE3LlzUK5cOYwcMVLCOKXjT58DdIura5d+Dh09i7cXb7ZLX0RUsvygvpmTBz9PDWoEalGzkg5eegEagVPgjmCvFWtPT094eXlZb8WFtcViwYrly3DixHG8PWESvLy9ix2byWjEvHlzYTKZ8K/X34C2lGl1R+NPnwMEVfRBsxj77Bhcu/0w8vLy7NIXERXt4aCuHuABQRDgoREQHqRFVIiOU+BuYMWK5UhLO4oJEyfDx8en2HYmkwnz3p2L3Nx7eOPNt0pd/3YGToM7yLCezbD/jzOS+7l55x6u3ryLoIrF/2ARke2KC+p8giBAy03gDiL8PRVuI4v4Yy9fvoQftmyGTqfDS6P+mdKOa90aw4ePwMwZ0xFdpw569+6DY2lp+PngQeh0erww9Hlr2169e6N37z62j1cChrWDdGoZjXI6LXKNJkn9+FfwRCX/4qdqiMh2pQU1OZrzrmAWFFQZq7/6utiPj397gvX/68bElNjWFTgN7iAVvMrhyabSdw6O6tOcTx5EDsCglgFBkH5TCYa1A/V6op6k41s9WhOvDWxtp9EQUT4GNSkNw9qB8qfCbdW6URj0Eo4nosIY1PLhvOuXKR/D2oGkToXP+ngH5n32ox1HRKRuDGqZ4TS4aAxrB5M6FT5l2VYGNpEdMKhJyRjWDiZ1KhxgYBNJxaCWK06Ei8WwdjB77QpnYBPZhkEtY5wGF41h7QRDuze1Sz8MbKKyYVDLG+tq8RjWTvBk09ro2SamyI/FNqiJRW/1gkcJF5N/EAObSBwGNbkTvi7IST6amID6tVLw2ZZfkJF1DTVC/PFsl8Z45elYaD08UF6vw7AZa2AWcR3wKcu2AgBfg01UDAa1QkieylbP95Rh7SQeHhq8PrgNXh/cBkaTGbqHLjbcJ74+ADCwiSRiUCuJ8y43qnScBneBh4M6X5/4+lj2dl9OiRPZiEFN7oqVtcywwiayDYNaiVhZi8WwliEGNlHZMKiVSeqStZq+wwxrmWJgE4nDoFYyVtZicc1axriGTVQyBjWpBStrmWOFTVQ0BrUb4Eu3RGNlrQCssIkKYlCT2jCsFYKBTXQfg5rUiNPgCsIpcVI7BrV7EQRB0vdPUNE0OMNaYRjYpFYManfE3eBiMawViIFNasOgdlPcYCYa16wVimvYpBYMaiKGtaIxsMndMaiJ7uM0uMJxSpzcFYPa/XGDmXisrN0AK2xyNwxqooJYWbsJVtjkLhjUasLd4GIxrN0IA5uUjkGtMtwNLhqnwd0Mp8RJqRjURMVjWLshBjYpDYOaqGScBndTnBInpWBQq5uadnRLwbB2YwxskjsGtcpxzVo0hrWbY2CTXDGoibvBxeOatQpwDZvkhkFNVDasrFWCFTbJBYOarKQW1irCsFYRBja5GoOaHiT8/U/K8WrBaXCV4ZQ4uQqDmsh2rKxViBU2ORuDmorE3eCiMaxVioFNzsKgpuJxN7hYDGsVY2CTozGoqUTcYCYa16xVjmvY5CgMaiL7YWVNrLDJ7hjUJA6nwcViWBMABjbZD4OaxBIEQdLPBl+6RarEKXGSikFN5BisrKkAVthkKwY1lR2nwcViWFMhDGwqKwY12YS7wUVjWFORGNgkFoOabMfKWiyuWVOxuIZNpWFQEzkHK2sqEStsKg6DmqTiG3mIx7CmUjGw6WEMarILrlmLxrAmURjYlI9BTUpkNBqxYsVyHP79d9y6lY2AgAB079ET8fHtimz/xRef4+CBA8jKOotOnZ7Cc88PdfKIC2JYk2gMbGJQk305b4OZ2WxGRX9/TJw0GcHBwTh+/DhmzZyOwMBANGjwWKH2ISEhGDxkCJK3bZMwPvvhBjMqE246Uy8GNdmf8M/bZNpyK0NYly9fHv2fHoCQkBAIgoDIyEjExNTD0b/+KrL9E0+0RcOGjeDp6Wmnz1UaVtZUZqyw1YdBTY5grw1mOTk5BR7X6XTQ6XQlHmswGHDixHHExsbZfH5nYliTTRjY6sGgJrkbNXJ4gft9+yUgIaF/se0tFguWLFmMKlWqoGmzZo4enl0wrMlmDGz3x6Amh7LTkvXiJUsLTFeXVFVbLBYsX7YU589lYeLEydCIXNJzNYY1ScLAdl8ManI8+6S1p6cnvLy8Sm1tsViwYvkynDhxHBMnJcHL21vCuZ2LYU2SMbDdD4Oa3NGKFcuRlnYUkyZPgY+PT4ltTSYT8vLyrDeDwQCNRgOt1jWxybAmu2Bguw8GNTmNdVe3hONFunz5En7Yshk6nQ4vjRppfTyudWsMHz4CM2dMR3SdOujduw8A4MMli7Fr105ru82bv0ebNk9g9MtjbB+vBILJZLS45Mzklr7eflh0YAPA5GHtGdgywqAmZ7h79y6ee3YIfgscjjyN3uZ+NHkGNLi6FB9/slLUNLiSsbImu2KFrVwManI6vumWaMrYBkeKwgunKA+DmkjeWFmTQ7DCVg4GNbkOS2uxGNbkMAxs+WNQk0sJkLjBzG4jkT1Og5NDcUpcvhjURMrBypocjhW2/DCoSQ7sdW1wNWBlTU7BCls+GNREysPKmpyGFbbrMahJVpx4URSlY2VNTsUK23UY1ETKxcqanI4VtvMxqEme+NItsRjW5BIMbOdhUJNcCYIg6WdRTT/HnAYnl+GUuOMxqIncA8OaXIqB7TgMaiL3wWlwcjlOidsfg5oUgbvBRWNYkywwsO2HQU3KwQ1mYnEanGSDU+LSMaiJ3BMra5IVVti2Y1CT0nAWXDyGNckOA7vsGNSkTJwGF4vT4CRLnBIXj0FN5P5YWZNsscIuHYOaFI3z4KIxrEnWGNjFY1CT8nEaXCyGNckeA7swBjW5AxbW4nHNmhSBa9j/YFATqQ8ra1IMVtgManI3nAYXi2FNiqLmwGZQk9vhPLhonAYnxVHjlDiDmkjdWFmTIqmpwmZQk/viNLhYDGtSLDUENoOa3BmjWjyGNSmaOwc2g5rcHtesReOaNSmeO65hM6iJlOu7777DjRs37Nonw5rcgjsFNoOa1EOww01+5syZi9q1I9G6dRtMmDARP/zwA27fvi2pT4Y1uQ13CGwGNalK/jS4lJsM7dq1E8eOHcObb74BgyEXkycnISwsHB06dLS5T65Zk1tR8ho2g5rIfVSs6I/IyEicP38BFy5cxLlz55An4jmpOAxrcjtKDGwGNZH7eOGFF/HTTz8hMDAAbdq0wYABT2PhwgXw9fW1uU+GNbklJQU2g5rUShAEST/rcv092bFjB3x9ffHkk+0RFxeHli1bwMvLS1KfXLMmt1XkGnb6RuC3pUDGpkLtpyzZiMiYhoiNjUVsbCyqV6+BRYsWOXSMDGpSN/fcYJaRkY7//W8lgoMrY/nyZahf/1F07NgJM2fOtLlPVtbk1gpV2JUbAIF1gWt/FW7socfFyl0xclh7vDogDvXr10fnzp0dNjYGNZH7qlevHmrWrIlatWohLCwMq1atwsGDBzF+/Hib+mNYk9srENgVqgG3zpbYfsqyrTiT/hcqVw5GzZo1HTImBjUR/i6OpVwUxW4jsaukpCnYvXs3fvvtN9SuXQtxcXH44IMPEBsbZ3OfDGtShQKBLaL9R598hg6tGztkLAxqIveWnZ2N0aNHIy4uFpUqVbJLn1yzJtXIX8PWCKX82FsswI10/JBmsfvrsBnURO5v3ry56NWrJypVqoSrV6/apU+GNalKn/j6eH1Im5Kn3u6cA/QVAL2PXS+cwqAmKih/N7iUmxzl5ORg3LhXERJSBbVq1UZISBWMG/cq7ty5Y3OfDGtSndYNw9A4OrT4K51dPwFUrG29a4/AZlATFcU9d4OPH/82Tpw4jvXr1yEt7Si++2490tPTMWHCRJv75Jo1qUr37j1w5MgR3L17F75ef+JmUFvkndsP1IgH9D7WKXBE9y9wnJTXYTOoiYrhpu+R+f3332PPnj0ICKgIAKhcuTI++eRjtGjREu++O8+mPhnWpCrr168rcP/r7YcxbEaVfy6cIghA/aFFHmtLYDOoidTHYrFAoyn4ey4IGlgsFpv75DQ4qZoj3/yDQU1UGvecBu/YsSOeeeZZpKb+gitXruDQoVQ8//zz6NSpk819MqxJ9RwR2AxqotIJdvgnRzNnzkD16tXQqVMn1K4dic6dOyM0tCpmzJhuc5+cBieCfa8lzqAmEkvq21zK8/fKx8cHH3zwARYuXIgrV66gUqVKkp8DGNZEf7NHYDOoiSifIAgICgqyS18Ma6IHSAlsBjWRetWo8Yio3/dTp07a1D/DmughtgS2IAjo06EFg5qoLASJ0+Ay+h377LNVDu2fYU1UhLIEtoeHBhp9BQY1kYwZjUasWLEch3//HbduZSMgIADde/REfHy7ItvfvXsXy5Z+iNTUQ9Dr9ejY6Sn07duv2P4nT05CcvI2AMDs2bORmJho1/FzNzhRMcTsEvfw0GDyqF6Ib1YX2/f/ia+37GFQE4nkzN3gZrMZFf39MXHSZHzy6f/w0ugxWPnpJ/jtt1+LbP/RRytw+/ZtLFr8IaZMnY7kbduwa9fOYvs/fvw4TCYTAGDhwg/K8mUQhZU1UQlKqrAfDuopi9bCnGeBBbZd6YxIdew0DZ6Tk1PgYZ1OB51OV+Cx8uXLo//TA6z3IyMjERNTD0f/+gsNGjxWoG1ubi72/LQb06bPgLe3N7y9vfHUU09he3Iy2rR5osihxMXFIjY2DhEREcjJycGgQYOLbLdq1f/K+Enex7AmKkVRgV1cUAPSLk1KpCb2utroqJHDCzzet18CEhL6Fz7gAQaDASdOHC/yPabPncuCyWRCzZph1sdq1gzDN9+sLba/jz76COvWrcOpU6fwww8/oH79euI/EREY1kQiPBjYEFBsUOdjYBM5z+IlS+Hp6Wm9/3BV/TCLxYIlSxajSpUqaNqsWaGP37t3D+XKlYeHh4f1MS9v70IV/IPKlSuHhIQEAMCNGzftvmbNsCYSKT+ws26Y0bZp8UGdj4FNVDJ7bQb39PSEl5eXqGMsFguWL1uK8+eyMHHiZGiK2JNSvnx5GAy5MJvN1sC+e/dOgT8ISiLlSmXFYVgTiWSxWPB4gzqIzMnDjgMlB3U+BjZR8aS+J3VZj7VYLFixfBlOnDiOiZOS4OXtXWS7qlVD4eHhgVMnTyI8IgIAcPLkSdSoUcPmsUrF3eBEIjx8wZNQfw+IXW2zx/thE5F0K1YsR1raUUyYOBk+Pj7FtitXrhxatmyFL7/8HHfv3MH58+ew+ftNiG/3pBNHWxAra6JSFHVlshp2vJY4kVoJAqBx0jVRLl++hB+2bIZOp8NLo0ZaH49r3RrDh4/AzBnTEV2nDnr37gMAGPrCi1i6dAlGjhwOvV6PTp2eKnYnuDMwrIlKUNIlRO355h9EaiRI3A5elrAOCqqM1V99XezHx789ocB9Ly8vjBv3mk3jWrZsOYYNe7HQ46+8Mg7vv/+eTX1yGpyoGGKu9e3I98Mmcnfu+W7WwAcffIB169YVeOy11/6FP/74w+Y+WVkTFaEsb8rBCpuIHrRmzVfo3r0HAgMDERsbizfffBOpqalYt+5bm/tkWBM9xJZ3z2JgE5WdIPG1W3K9tG+tWrWwcuWnGDRoMFq3jsPRo2lYv34d/Pz8bO6TYU2qcvbsWYwYMQKXL1+BVuuBN954E7169bR+3GKxYMDgoTh2LA0C8tA6tgXmzp0r6kmBgU1UNhoBsLjHm27hyJEjBe6XK1cOI0aMwJIlS7B8+TKcPXsWZ8+eRb16tl3ZjGFNqqLVajFr1iw8+uijuHjxItq0eQIdOrSHt7e3taIeN3EOqgb5oVpFDZ577nls3LgR3bp1E9U/A5tInWJj4yAIAiyWwtde6NatO4D7MwHXr1+zqX+GNalKSEgIQkJCAADBwcEIDAzA9evX4eXlZZ36rhrkh+oBHjCbzTAYcss81cbAJhLHmbvBHe3GjesO7Z+7wUm1fvnlV5jNeQgNDS20Rv3MM8+iVq3a8Pb2QefOncvcN3eJE5Uuf8layk0tWFmTKl27dh0jR47E+++/V+RmspUrP4XBYMDIkaOwc+cuxMe3LfM5WGETlUwjCLC44Qazc+fOYcaMGfj1119x69btAh/7/fffbOqTYU2qk5ubi0GDBmLcuFdQtXaTYnd96/V6dO3aFZs2bbIprAEGNpEaDR8+HJ6eXhg3bpzoNxgpDcOaVMVisWDUqJcQF9carTr0KxTURqMR58+fR40aNWA2m7Fly2Y0btxY0jkZ2ERFk1oXy7OuBn799TdkZKRDr9fbrU+uWZOq7Nu3D2vXrsW36zei11Nt8ELfeGSfP4oxY8YiNfUXGI1GDB36Alq0aIlWrWJRoYIvhg4dKvm8XMMmKsxd16yjo6Nx8eJFu/bJyppUpXnz5vg9/XKhinrhwgXWNtu2bXXIuVlhE6lDt27dMGDAALz44jBUrhxU4GO2bFgFGNakIrZcmczeGNhE/3Cnl249aPny5QCAuXPnFnhcEASGNVFJ5BDU+RjYRPe5627ww4d/t3ufDGtye3IK6nwMbCL3rawdgWFNbk2OQZ2PgU3kPjp27IQtWzYD+OfSo0VJSbFt4yjDmtyWnIM6HwOb1MydXrr14osvWP//pZdG2b1/hjW5JSUEdT4GNqmV5JdfyehXul+/ftb/HzhwoN37Z1iT21FSUOdjYBMp26ZNm0S1425wIigzqPMxsEltpG4wk1Nl/dZbiaW24Uu3iKDsoM7HwCY1EeA+8+COeLnWgxjW5BbcIajzMbBJLTRuVFk7Gq8NTornTkGdj9cSJ6IHsbImRXPHoM7HCpvcnTutWTsaw5oUy52DOh8Dm9yZe/22OhanwUmR1BDU+TglTkSsrElx1BTU+VhhkzsSpF4Vxc1/7x/EsCZFUWNQ52Ngk7vhbnDxGNakGGoO6nwMbHIn3GAmHtesSREY1P/gGjaR+rCyJtljUBfGCpvcAStr8RjWJGsM6uIxsEnpNGBai8VpcJItBnXpOCVOpA6srEmWGNTiscImpeI0uHgMa5IdBnXZMbBJiRjW4jGsSVYY1LZjYJPSMKzF45o1yQaDWjquYRO5J1bWJAsMavthhU1KIQiCpN9zi4qeIxjW5HIMavtjYJMSSL00OATAYrfRyBunwcmlGNSOwylxIvfByppchkHteKywSc6k7i8D1FNZM6zJJRjUzsPAJrnSSJwGtwhA6T/R7oFhTU7HoHY+BjbJkdQNZmp6P2uuWZNTMahdh2vYRMrFypqchkHteqywSU7ssRtcLRjW5BQMavlgYJNc2GPNWi04DU4Ox6CWH06JEykLK2tyKAa1fLHCJlezx0u31IJhTQ7DoJY/Bja50v01aym7we03FrljWJNDMKiVg4FNrsINZuIxrMnuGNTKw8AmNdj8/Sbs3LkTp0+fwmMNG+LNNxOLbZuRno7//vcjnD59ChUqVEC/hP5o0+YJ5w32IQxrsisGtXIxsMnZNJC4G7yM7SsGBKB3nz44/PvvuHrtarHt7ty5g1mzZqBfQn882e5JpGekY/q0aQiuHIzoOnVsH7AEDGuyGwa18jGwyZnsNQ2ek5NT4GGdTgedTleoebNmzQEAJ0+eLDGs09KOQqvVoUOHjgCA2rUj0axZMyQnb2NYk7IxqN0HA5ucRfj7n5QeAGDUyOEFHu3bLwEJCf1t7tViseDhuj3PYsGZ06ds7lMqhjVJxqB2PwxsUpLFS5bC09PTer+oqrosIiOjcO9eLjZ/vwlPtu+AEydO4OCB/fD19ZM6VJsxrEkSBrX7YmCTo9nrCmaenp7w8vKyz6AAVKhQAW8l/hv/W/kpVq/+EtWqVccTT7TF8ePH7XaOsmJYk80Y1O6PgU0OJePLjUZHR2P6jJnW++/Om4u6des67oSl4OVGySYMavXgpUnJXZjNZhgMBuSZzbDkWWAwGGAyGotsm5mZAaPRCENuLrZt24o///wDnbt0dfKI/8HKmsqMQa0+rLDJEaTuBi/rsV9/vQZrvlptvT940ADUrRuDpClTMXPGdETXqYPevfsAAL7ftAkHDuyH2ZyHqKgoTJqchICAANsHK5FgMhnL+lI1UjEGtbp9vf2w6MAGgMnD2jOwqZC7d+/iuWeHIPCp/0Cj8yz9gGLkGXNw9fs38PEnK+26Zi1HrKxJNAY1scIme9II9282U9HTD9esSRQGNeXjGjaR87GyplIxqOlhrLDJHqS+RaaanoUY1lQiBjUVh4FNUjl7g5mScRqcisWgptJwSpzIOVhZU5EY1CQWK2yylUYQoGFpLQrDmgphUFNZMbDJFpwGF49hTQUwqMlWDGwqK4a1eFyzJisGNUnFNWwix2BlTQAY1GQ/rLBJLA1YMYrFsCYGNdkdA5vE4DS4ePyjRuUY1OQonBInsh9W1irGoCZHY4VNJREEQdJzjpqerxjWKsWgJmdhYFNxOA0uHsNahRjU5GwMbCoK33VLPK5ZqwyDmlyFa9hEtmNlrSIManI1VthUgMRpcDVV1gxrlWBQk1wwsCmf8Pc/KcerBafBVYBBTXLDKXGismFl7eYY1CRXrLBJ6gYzi4qeyhjWboxBTXLHwFY3vnRLPIa1m2JQk1IwsNWLYS0e16zdEIOalIZr2EQlY2XtZhjUpFSssNVHgMTLjapoNzjD2o0wqEnpGNjqIvUtMi32GogCcBrcTTCoyV1wSpyoMFbWboBBTe6GFbY6cIOZeAxrhWNQk7tiYLs/hrV4DGsFY1CTu2NguzeNIEAj4TnLoqLnO65ZKxSDmtSCa9hErKwViUFNasMK2z1xGlw8hrXCMKhJrRjY7odhLR6nwRWEQU1qxylxUitW1grBoCa6jxW2++BFUcRjWCsAg5qoIAa2exAEiZcbVdHzIMNa5hjUREVjYLsBiWvWKro0ONes5YxBTVQyrmGTWrCylikGNZE4rLCVi7vBxWNYyxCDmqhsGNjKpBHu36QcrxacBpcZBjWRbTglTu6MlbWMMKiJpGGFrSzC3/+kHK8WDGuZYFAT2QcDWzm4Zi0ew1oGGNRE9sXAVgauWYvHNWsXY1ATOQbXsMmdsLJ2IQY1kWOxwpY3ToOLx7B2EQY1kXMwsOWLlxsVj9PgLsCgJnIuTomT0rGydjIGNZFrsMKWHwHSLu+tpmdOhrUTMaiJXIuBLS/cDS4ew9pJGNRE8sDAlg9nbzDb/P0m7Ny5E6dPn8JjDRvizTcTi2179swZfPTRCmRmZkCr1aFJkyZ47vmhKFeunO0DloBr1k7AoCaSF65hq1PFgAD07tMH7do9WWrb999/D1WrVsWyZSswd+48nDp1Cl+v+coJoywaw9rBGNRE8sTAdr383eBSbmXRrFlzNG3aDBV8fUtte+nSRcS1bg2tTgdfPz80adIEp0+ftvVTlYzT4A7EoCaSN06Ju5a91qxzcnIKPK7T6aDT6SSMDOjWrTt27dqFsJphuHv3Lg4cOIB2T5ZekTsKw9pBGNREysDAVr5RI4cXuN+3XwISEvpL6vOxho2weNFCPPPMYOTl5eHxx5uibdt4SX1KwbB2AAY1kbIwsF1DgMQNZn//d/GSpfD09LQ+LrWqvn37NqZNnYL+/fujQ4eOuJebi48+WoEF89/Hq6/9S1LftuKatZ0xqImUyZY17Hc/5xq2VIKEWz5PT094eXlZb1LD+uLFCzAYDHiqcxdodTr4+Pigffv2SE1NldSvFAxrO2JQEylbWQM7aelW7Pn9pGMH5cY0giD5VhZmsxkGgwF5ZjMseRYYDAaYjMZC7UKrhqJ8+fLYsmUzzGYzcnJykLxtG8LCwuz1qZcZp8HthEFN5B7KOiU+6cMfsO2D4aW2I9f7+us1WPPVauv9wYMGoG7dGCRNmYqZM6Yjuk4d9O7dB+U9PfFW4r+x6n8r8cXnn0Gj0SAqKhqjX37ZZWNnWNvJpVt5DGoiN1GWwD526rIzhuSWnH1RlISE/sVuPBv/9oQC96OjozFt+gxbh2Z3DGs78SknAL4eqFxBw6AmcgNiA1vQCMjLy4NG5NQ5/YNvkSkew9pOvMtp4O2aq9ARkYOICeyebWIY1DZiWIvHnzAiohL0ia+P5RP6opyucG1Tq3ogEp9t64JRkdqwsiYiKkXvtvXRoHZVTFq6BX9mXITJnIfucXUxpn8sQgIruHp4iqWBtIpRTdUmw5qISISIaoFYNXWgq4fhVmy5vvfDx6uFmv4wISIiUiRW1kRE5BLcYCYew5qIiFzCXu+6pQYMayIicglW1uJxzZqIiEjmWFkTEZFLcDe4eAxrIiJyiYff6tKW49WC0+BEREQyx8qaiIhcgrvBxWNYExGRS3A3uHicBicS4ezZs+jSpQuaNm2Gli1b4ptvvi3UZujQF9CqVSs0b94Cr776GvJKeR9kIrXL32Am5aYWrKyJRNBqtZg1axYeffRRXLx4EW3aPIEOHdrD29vb2ua9996Fr68vLBYLnn32OWzcuBHdunVz4aiJyF0wrIlECAkJQUhICAAgODgYgYEBuH79eoGw9vX1BQCYzWYYDLmq+qufyBaCxDVrNf2KcRqcqIx++eVXmM15qFatWqGPDRnyDGrVqg1vbx907tzZBaMjUg7BDje1YFgTlcG1a9cxcuRIvP/+e0V+fOXKT3HsWBosFgt27tzl3MERkdtiWBOJlJubi0GDBuLVV8ehWbNmxbbT6/Xo2rUrNm3a5MTRESlP/m5wKTe1YFgTiWCxWDBq1Eto3bo1nn766UIfNxqNOH36NID7a9ZbtmxGZGRtZw+TSFE0giD5phbcYEYkwr59+7B27VrUqxeDjRs3AgA+/PBDLF68BEOHDkV0dBSGDn0Bd+7cgcViQWxsLIYOHeriURPJG19nLR7DmkiEFi1a4MaN64UeX7hwgfX/t23b6swhEZGKODWsX/v8PHYdvePMU1q1ifbGvAFVXHJuIiIqjJW1eE4N611H78BotkDn4dyvsNFscdkfCUREVDQBEsPabiORP6dPg+s8BKRMiHTqOeOmH3Pq+YiInC07OxsZGem4dOkSTEYjtDodKleujPDwCOsFe0i5uGZNRKRQN2/eRHLyNuzcsR0XLlwotl1ISAjaxrdDfHw7+Pn5OXGEJdNAgEZCfSzlWKVhWBMRKYzBYMDq1V9i44YNMJtNiIqKRtu27RAeEYHQ0GrQ63QwGI3IyjqLjPR0pKYewuefrcLqL79El65dkZDQH3q93tWfBtesy4BhTUSkIBkZGVgw/z1kZWUhNq41evXsjZphYUW2DQ4ORqNGjdG3XwJOZmbim2/XYv26b3Ho54MYM3YcwsPDnTz6ghjW4jGsiUgR0s9eRcqvmbh68w5a1H8ETWOqQ+vh4ephOdWff/yB2bNnonx5TyT++208/nhT0cfWDAvDq6/+C7GtYrFkyWIkTZ6IxMTxqBsT48ARk70wrIlI9v7zv51455OdMJrM1sdqVQ/Em0OeQJ/4+qoI7YyMDMyaNRMBAQGYnDQVQUFBNvXzeNNmqBkWjilJkzB79kwkTZnmsgpbI/Fdt6QcqzS83CgRydq6XX9g+orkAkENACfOXMXwmV+j2fML8OXWX2Eym4vpQfkMBgMWzH8Pnp6ekoI6X1BQECYnTUX58p5YMP89GAwGO420bO5PgwsSbi4ZtkswrIlI1j74ak+JH1dDaK9e/SWysrIwctRLkoM6X1BQEEaOHIWsrCysXv2lXfokx2FYE5Gs/X7ivKh27hraN2/exMYNGxAb17pMa9RiPN60GVrFxmHTxg24efOmXfsWg+9nLZ4i1qwn/PsNbNn8Pc6eOY2tO3ajXv1HXT0kySwWCzLOXcPZi87/BSFSkpxcY5na54f2lGVb8erA1ni+WxNFr2knJ2+D2WxCr569S23bvXt3XLx4ERqNBj4+PpgzZw4aNGhQ4jG9evXGT7tTsH17Mnr1Kv0c9sTd4OIpIqy7dOuJl8aMQ48uHV09FLvYduA4Xv7PNzh/5Zarh0LktrIuZ+P19zdgwpLN+GLGYLRtHOHqIdlk547tiIqKLvblWQ/69NNP4e/vDwBYv349RowYgX379pV4TFhYOCKjorDDRWEtZZOYmsJaEdPgLVq2QtWqoa4ehl1s2ZuGvokrGdRETnIv14Rer3+MX9OyXD2UMsvOzsaFCxfQqFFjUe3zgzr/WEFkmjVq2BgXLlzArVt8XpIrRYS1u7BYLJjx3+2wWCyuHgqRqlgAjJ27ztXDKLOMjHQAQHiE+FmBYcOGISoqCtOmTcOyZctEHZPff/75nEXaTnBB9B8j7kAR0+DuItdowm/Hz7l6GESqdPzMFVcPocwuXboEAAgNrSb6mPyAXrVqFSZNmoS1a9eWekx+/xcvXrRhlLbjmrV4rKydSFDV3kUieVHib5/JeH9znV6nK/OxgwYNwo8//oirV6+W2lav1xU4H8kPw9qJyum1eLJpbVcPg0iVmj/6iKuHUGbav0PaICJEb9y4gfPn/3mZ23fffYeAgAAEBASUeqzBYCxwPmfR2OGmFoqYBn/jtVeQvHULLl26iAEJveDj44O9B39z9bBsMmFoOxz44zSy7+S6eihEqqHXemD+v3q4ehhlVrlyZQBAVtZZBAcHl9g2OzsbQ4YMQU5ODjQaDSpVqoQ1a9aIWtfNyjoLAKWew944DS6eIsL6P/Ped/UQ7KZhVCh2LB6JtxZuwv4jp3HrLkObyJFqVqmITe+/gNAg+byPs1jh4X9v/EpPL3VHeI0aNbBr1y6bzpORnl7gfM4idZMYN5iRQ9WqXglfv/MMzOY83Lid4+rhEMlaeM/ZNh3XsXkk3h7aDg1qV7XziJzH19cXISEhSE09hL79Ehx2ntRfDiEkJAQVKlRw2DlIGoa1C3l4aBDo5+3qYRDJmp9Pedy8fU90+65xdfDWM23xaK0qDhyV8zzRNh5ffP4ZTmZmirowSlllZmbgWFoaBgwcZPe+S8N33RJPTevzRKRAjaLFXRCpa1wdpCx7CaumDnSboAaAdu2ehIeHFt98W/pLsGzxzTdrodVqER/fziH9lyR/zVrKTS0Y1kQkaxOebwcPTfFPVe4a0vn8/PzQpWtX7E75EQcPHrBr3wcP7MdPu1PQuUtX+Pkpb01fTRjWRCRrTepWx8eTExAc4GN9TBAEtw/pByUk9EdoaCiWLF6Ey5cv26XPy5cvY8mSRQgNrYaEhP526bOsWFmLxzVrIpK97q1j0O7x2jh68hKu3ryLx2Oqo2IFT1cPy2n0ej3GjB2HpMkTMSVpEiYnTZX0vtaXL19G0uSJuHfvHhL//Tb0er0dRyue8Pc/KcerBStrIlIEb089Gtephg7NI1UV1PnCw8ORmDge169fQ+Jbb9g8JX7wwH4kvvU6bty4jsTE8QgPD7fzSMXL32Am5aYWDGsiIoWoGxODpCnT4O3thdmzZmDevDk4mZkp6tjMzAzMmzcHs2fPhLe3N5KmTEPdmBgHj5jshdPgREQKEh4ejnf+bw5Wr/4SmzZuwE+7UxAZFYVGDRsjPCICoaHVoNfrYDAYkZV1Fhnp6Uj95RCOpaVBq9Wie4+eSEjo77Kp7wcZcnMkrTsbctVznQqGNRGRwuj1egwePATdunXH9u3J2LE9GV988Vmx7atUqYKBAwehbXw7Wez61mq18Pf3x/tTXpbcl7+/P7Ra948ywWQyOu3NlRtPPgGj2QKdh3MXGvLPeWhKLaeel4jIWbKzs5GZmYGLFy/CZDRCq9MhODgY4eERsrwymcFggMlkktyPVquVxSyBozn1z5E20d7YdfSOM08JANB5CGgTzSuFEZH78vX1RYMGj7l6GKLp9XpVhKy9OLWyJiIiorLjbnAiIiKZY1gTERHJHMOaiIhI5hjWREREMsewJiIikjmGNRERkcwxrImIiGSOYU1ERCRzDGsiIiKZY1gTERHJHMOaiIhI5hjWREREMsewJiIikjmGNRERkcwxrImIiGSOYU1ERCRzDGsiIiKZY1gTERHJHMOaiIhI5hjWREREMsewJiIikjmGNRERkcwxrImIiGSOYU1ERCRzDGsiIiKZY1gTERHJ3P8DdFHjhBqOVb0AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 480x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAc0AAAHXCAYAAADeCAHgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABHAklEQVR4nO3dd3hUZdoG8PtMzUx6MikkJBCa36KsoHQIoAaI9CJFXHVVQIp1bVhQV6UouuuuBAGRtcAK0gQbIFgoEgQsoIKrgEAIhEzqJJkkU97vj5ghFWYmk5wp9++6cl3MzJkzzxyG3Lzv+8w5ktVqESAiIqLLUshdABERka9gaBIRETmJoUlEROQkhiYREZGTGJpEREROYmgSERE5iaFJRETkJIYmERGRkxiaRERETmJoErWAOXPmYObMmU5tu3r1avTv37+ZK3LPW2+9jU6drkBCQiJ++OEHucshanEquQsgIt9gsVjw2GOPYdOmjejbt6/c5RDJgiNNIoLVaoUQlz4NdU5ODsrLy9G5c2e3XsNisbj1PCJvwtAkAtClSxe88so/MGjQdWjVKgHjx9+E/PwC/O1vDyE5ORndul2D/fv3O7Y3mUy477770anTFejU6Qo88MCDKC0tdTy+d+9e9OnTFwkJibjllr/AZCqp9XonTpzEpEmT0K5de1x11VVYtGgR7Hb7ZevMyMjAiBEja923YcNGdO/eAwDw/fff44Yb0tC6dRJSUtph0qRJje4rPDwCy5cvR+/efdCqVQJKSkoareuHH35Ajx49AQCdO1+Jq6/uCgAoKSnBww8/giuvvArt23fA3XffjaKiIgDAqVOnEB4egVWrVqFr12740586O2ocMWIE2rRpi65du+Gtt9521LRgwQJMmjQJDz/8CJKTk3HllVdhw4aNjsftdjuWLl2K7t17IDGxNbp1uwY7duwAAAghHI8lJydj+PDh+OWXXy57TIlcwdAk+sOmTRuxatW7OHbsKM6ePYu0tDQMGjQQJ0+exIQJN+HBBx90bDtnzhycOHECmZn7sG/f1/j11//h8cefAAAUFBTi5ptvxrRp03D69Cn85S+34P3333c8t6ysDKNHj8LAgQNx7NhRfPrpp9iwYSNWrVp12RonTJiAzMxMZGVlOe5bu3atIxwfeeRRpKen4/TpUzh27Cjuu+++S+5v3br12LRpI7KyzkCpVDZa19VXX43MzH0AgJ9//gk//PA9AGD27HtQUFCAvXv34PDhH2CxWPHII4/Ueo1PP/0UX375BQ4f/gE5OTkYM2Ys7rzzLpw4cRz//e9qLFiwAF9++ZVj+507P0ffvn1x8uRJPPXUk7jvvvtgMpkAAMuXL8eSJa/jjTfeQFbWGWzZshlJSUkAgBUr3sS7776LtWvX4MSJExg5ciQmTZqMysrKyx5XImcxNIn+cOedd6F169YIDw/H4MGDERUVhVGjRkGpVGLcuHH4+eejqKyshN1ux/vvr8Ozzz6DqKgoREdH4+mnn8aaNWtgt9uxbdtWxMe3wp133gGVSoUbb7wRAwYMcLzOtm3bER4egVmzZkGj0SApKQkzZszAunXrL1tjbGwsBg0ahPffXwcAyM3NxRdffIHJk6tCU61W4cyZMzh37hy0Wi369et3yf3df/99aNWqFbRarct1GY1GbNmyBS+//DIiIiIQHByMJ598Ahs3boLNZnNs99hjjyEiIgJ6vR5r1qxFv359MW7cWCiVSnTu3Bm33HIL1q1b59j+6quvdjw+eXJV6P3223EAwJtvrsTjj89Bt25dIUkSkpKScMUVVwAAVqxYgSeeeALt27eHSqXCjBkzUF5ejoMHD172uBI5i41ARH+IjY1x/Fmv19W6rdPpIIRAWVkZKisrUVlZieTkZMfjbdu2RUVFBfLy8nDu3HnH6KdaUlISKirKAQCnT5/G0aNHaz3fbhdITEx0qs7Jkydj0aJF+NvfHsT69evRq1dPx+stXpyBF19ciIEDByEiIgLTp0/D9OnTG91X69atHX92ta5Tp07Dbrfj6qv/XOt+hUKBnJycGq9x8VicPn0a27d/Vus1bDY7+vTp47gdFxfr+LMkSdDpglBSUjXSPHPmDNq3b99gPadPn8b06XdDqbw4FqistCA7O7vhN0/kBoYmkYsMBgM0Gg1Onz6N2NiqX/CnT5+GVqtFdHQ0WrWKx5kzZ2o9JysrCzExBgBAYmIiunbtip07d7j1+sOHD8ODDz6I7777HmvWrMXUqXc5HmvXLgXLli2DEAKZmZkYPXoMevToiW7duja4L4XiYsC4Wlfr1olQKBQ4duwY9Hp9vcdPnTr1x2tItV5jxIgR+M9/Vjr1GnUlJSXhxIkT6NmzZ73HEhMTsXDhAqSlpbm1byJncHqWyEUKhQITJtyE5557Hvn5BcjPz8ff//4cJk2aBIVCgSFDhuLcuXN46623YbVasW3bNuzatcvx/PT0obhw4QLeeGMFysvLYbPZ8Ouvv2L37t1Ovb5Op8OoUaPw/PPP45dffsGYMWMcj7333nu4cOECJElCeHg4FApFrZHXpbhaV1xcHIYPH45HHnkEeXl5AKo6bD/88MNGX2Py5EnYtWsXNm/eDIvFAovFgsOHD+PQoW+dqvGOO/6KhQtfxOHDhyGEwJkzZxzNPtOmTcW8efPx66+/AgCKi4vx8ccfO9ZDiTyBoUnkhoULFyI5ORm9evVCr1690a5dO8yfPw8AEBUVif/+dzWWLl2K5OQ2eOeddzBhwgTHc0NCQrB582Z89dVX6NLlz0hJScFdd01FTs4Fp1//5psnY+fOnRg+fDhCQ0Md93/55Zfo168/EhIScfPNU/D888/hz3/+8yX2dJE7db3++hKEh4dj0KDr0Lp1EtLTb8T33zd+0oOEhARs3LgB//nPW+jU6Qp06NARDz/8iNPBNmPGDNx1153461/vQGJia4wePQZnzlQ1RU2fPh1TpkzBX/5yK1q3TkLPnr2cWicmcoVktVou/eUsIiIiAsCRJhERkdMYmkRERE5iaBIRETmJoUlEROQkhiYREZGTGJpERERO8uvQrD7t2eUueUREROQMvw5Ns9mMv95+K8xmc5P3lZ+X64GKAg+Pm3t43NzD4+aeQD1uQgiYTEVOXZavml+Hpie5clDpIh439/C4uYfHzT2BeNyEECgszENZaQmsVucvkM7QJCKigFIdmJUVFYiIjIZGo3X6uQxNIiIKGHUDU6sNcun5vDQYEREFDEmSoNUEQa8PcTkwAY40iYgoAAghUF5e1RSqD3YvMAGGJhER+bnqKdmiwgLYbLYm7YvTs5dgswt8c6IM+4+bYSw0wxBhRK/2OvRsp4eyxtXoiYjIO9Vdw1QqlU3aH0OzEft+K8Vrn+Uhp9gKoKolW6EoxoffFyMuTIV7B0ejT4dgmaskIqLGNLXppyEMzQZ8cbQE8z7MhRAC7WM1GPbnUASJEpRLIfj4BxNO5FZi7sYLeHJkDK77U4jc5RIRUQOEEBBCeCwwAYZmPReKrXjx46rAnH1DNMZeGwZJkmDMLYchJhyjrwnDpkPFyNiZhxc/zsVVrYMQE8rDSETkLYQQsNttUCpViIw0QJI8t5zGRqA6PvqhGBabQHqXUIzrHl7vYEuShHHdwzG0SygsNoGPvje5/BolJSUYN248+vTpiz59+mLHjp1ISWnneHz58uVYsGABAODEiZMYO3YcBg4chBEjRuLUqVP19rd7926MGjUakyZNQpcuf8bLL7+CN99cidTUAbj++huQl5cHAFiyZAmuvbY7+vXrhwcf/FvV+/3oI1x33fXo3z8VEyZMRGFhocvvh4iopVitdryztwAj/vk7ej33G3r+/Tdcv/AEFnx0AYVlVseUbEF+HoQQHg1MgKFZz7YjJQCAcd3DLrnd+D8e33rE9dDcuXMnoqKisG/f1/j6673o2bNHo9s+9NBDePXVf+Krr77EI488jLlzn25wuyNHjiAjIwOZmfvw+uuvo6KiHLt378KAAQOwZs0aAMCLL76EXbu+wt69e/Hss88AAPr164/PP9+JPXt244Ybrscbb6xw+f0QEbWE/50rx/UvncQ/thqRlW9BhUWg0iqQX2rD2v1FSHvxd6zenY3KigqEhtUf9HgC5xVrsNkFjCYr9BoF2sde+rRK7WO10GsUMJqssNmFS920V155JR5//AnMnfs0Ro4cgZ49eza4XUlJCfbt24dbbrkFQNWUg17fcPNRz549YDAYAAAJCQkYPHjwH6/VGQcPHgQAXHvttZg+fTrGjBmD4cOHAwCyss7g9ttvR25uLsxmM7p37+70+yAiailZ+ZW4/Y0smC0COrWEoV1CcUvfCOg1Ej75wYS13xQh12TDK5+VIVQfhdHxnlnDrIuhWYOEqulXi01cNghtdgGLTUAhSXD12ycdOnTAnj27sXXrVjzxxJOYOHFCrf8RVVRUAKjq2I2JicGePXtqPT8rKwuTJ08GADz44N8QGxtT69yJCoUCGo2m6j1JCthsVSdjXrfufezevQcff/wRMjKW4Msvv8Bjjz2Ghx56GDfccD22bt2K1av/69qbISJqAU+sz4HZIhAVrMS6e5IRHXIxvqYNisZtvYNx18qz+PGcHYu2FmLkNZFQKDw/mcrp2RoUCgntYzSw2AQOnLz05cQOnDDDYqvqrnV1CuDcuXPQ6/WYMmUKZs6cicOHjyAsLAynT5+GxWLBtm3bAABhYWGIjY3Fp59+CgCw2Wz4+eef0bp1a+zZswd79uzB+PHjnHpNu92OrKwsDBo0EC+88AKysrJgs9lQXGxCQkIrCCHw3ntrXHofREQtIa/Eih/PlgMA/jmlVa3ArO6Q1QYFYeXUNghSSyipsGPLd64vnTmDI806RnYLxT+3VeC9zEJc21YHtbJ+IFpsAu/tL6zavmuoy6/x008/Y+7cuVAqFQgK0mHx4tcwaNAgjBw5CnFxcejYsaNj2xUrVuDBBx/ECy+8AIvFijvu+Cs6d+7s8mvabDZMmzYdJSUmCAE89tijUCqVmDPnMUyaNBmRkZHo378fTp8+4/K+iYia0/oDRbDbgfhwFa5O1jnur276USpVCAuLgEajQr+Owdj5cwk++LYYY64N93gtktVqER7fq5coKyvDX2+/FW+9/S70er1TzzFX2nHXm1nIKbaiX8dg3D8kGtEhKhhzc2CIiUNeiRX/2p6Hvb+WIi5MhTfvag2dhgP2xlQfN3INj5t7eNzc4+3Hbf6HF/D+N0XomhyEt6YlAWj8xAXLvsjD65/no32sBhvubePxWjjSrEOnUWD+hHg8suYc9v5aiszjZejZTodQVQVM1vPYf9wMu6iaV58/IZ6BSUTUzEK0Vb9ni8xV/RmXOtPPhT/O4qZVN8+pTvkbvwFtDRosvi0RQ7uEQiEB+34rw9afLdj3WxlUSiC9SygW35aItgaN3KUSEfm99C5VZ1773ViJnCILyspKGj013udHSwEAvds7N7voKo40GxEXpsKjw2Iw8/oofH+qHOeNBYg3RKJrmyCEBjXthL9EROS8Tq2CkBSlxpl8Cx5ecx5vTU2ARqOFWl174PLaDiMKSm1QKyVMHRjVLLVwpHkZoUFKpF4RjOs6qZF6RTADk4hIBg+lGyBJwJGsckzIyELmSYvjsRO5Fbj33bN486sCAMD4HmHQN9PSGUeaRETk1YQQuDq+HPcMUCNjlwUncitx77vZ0GkkKCUJpRV2VHe0Du0SgjnDY5utFoYmERF5rZpNP39JjUefPwH/+iwPh343w1wpAAhIADrEaTB1QBTS/+z61wBdwdAkIiKvVVxUWKvpp3MisOyviSivtOP3vEpUWATaGNSI0LdMnDE0iYjIawWHhCBIp6vXJRukUeD/WjXP+WUvhY1ARETkVYQQKDEVQ9jtUKnUHruAtCcwNImIyGtUr2GWlppgsVrlLqcehiYREXmFumf6qb5akzdhaBIRkewudWo8b8JGICIikp0kSdCotdDrQ7w2MAGGJhERyUgIgcrKCmi1QQgOad7vWHoCp2eJiEgW1VOyhYX5sNlscpfjFIYmERG1uFprmBFRUCp947zeDE0iImpRvtL00xCGJhERtSghBITd7nOBCbARiIiIWogQAna7HUqlEpFRMZAkSe6SXMaRJhERNTtH00+BEUIInwxMgKFJRETNrOYaZkhouM8GJsDQJCKiZuTLTT8NYWgSEVGzsVgqYams9IvABNgIREREzUAIAQDQaLQwxMRBofCN72FeDkeaRETkUdVTsiUlxQDgN4EJMDSJiMiDaq5hajRaucvxOIYmERF5hL81/TSEoUlERB5RVlri14EJsBGIiIg8RB8cAo1WC7VaI3cpzYYjTSIicpsQAkWF+bBUVkKSJL8OTIChSUREbqpewywvN8Mu7HKX0yIYmkRE5LJAaPppCEOTiIhcVlxUEHCBCbARiIiI3KAPDkGQTh9QgQlwpElERE4SQqCkpBhCCKjVmoALTIChSURETqhewywtMcFqschdjmwYmkREdEl1m37UGv/+WsmlMDSJiKhRgdol2xg2AhER0SWpVRro9SEBH5gAQ5OIiBoghIClshIarRYhoWFyl+M1OD1LRES1VE/JFhbmwW63yV2OV2FoEhGRQ801zPCIKL+6gLQnMDSJiAgAm36cIcuaZmVlJV599R84m5UFjUaDsLBwTJs2HfGtWtXb9tChg3j3nbdht9uRnNwGs2bfA71eL0PVRET+zW63w26zMzAvQbaRZlraYLz6r9ew6OV/oEePHli69PV625SbzVj6+hI88uhj+PdrGYiMjMSG9etkqJaIyH8JIWC326BUKhEVHcPAvARZQlOj0eCaa66FJEkAgI6dOiE390K97b77/ju0bZuCxMTWAIChQ9Oxd++eRvdrsVhQVlbm+DGbzc3zBoiI/ERVYNpRkJ8HIYTj9zI1zCu+cvLJxx+je/ce9e43Go2IiYlx3I6JjUVBQSFstqr/EdW1adNGrF/3fr3784wXUKbTNalGm80KY25Ok/YRiHjc3MPj5h4eN9dUByZQNdLMM9YfvAQKQ0ycU9vJHpobN27A+fPn8fQzzzZ5X2PHjsOIESMdt81mM2bOmI5oQ2yT10GNuTlOH1S6iMfNPTxu7uFxc17Nph+FQoGY2Po9JVSfrN2zW7Zsxjf79+OJJ5+CVqut97jBYEBubq7jdu6FC4iMjGhwlAkAarUaer3e8aNr4uiSiMhfVVZWoLKyEhGR0ZAkfpHCWbIdqY8+3IK9e/bgqblPIzg4uMFtunbthpMnT+Ds2SwAwLZtW9G3X/+WLJOIyK8IISCEgFYbhBhDHJt+XCTL9GxeXh7eeedtxMXF4e/PPgOgapQ4f8FCrF3zHiKjojBkyFDodDrMmDELi156ETabHUnJSbhn9r1ylExE5POqp2TVag1CQsKgaGTWjhonS2hGR0fj/XUbGnxs0uSba93u3qMHuveo3yRERETOq7mGqdeHyF2Oz+JENhGRn+OZfjyHoUlE5OdKS00MTA+R/SsnRETUvIKDQ6HVBEGt0chdis/jSJOIyA8JIVBUmA+LxQJJkhiYHsLQJCLyM9VrmOXlZl4P08MYmkREfoRNP82LoUlE5EeKigoYmM2IjUBERH5Erw+GTqdnYDYTjjSJiHycEAKlpSYIIaDRaBmYzYihSUTkw6rXMEtMxbBaLXKX4/cYmkREPqpu049aza+VNDeGJhGRD2KXrDzYCERE5KNUShX0kSEMzBbE0CQi8iFCCFgsldBotAgNi5C7nIDD6VkiIh9RPSVbWJAHu90udzkBiaFJROQDaq5hhkdEQaHgr2858KgTEXk5Nv14D4YmEZGXs9vtsFltDEwvwEYgIiIvJYSAEAJKpRLRhlhIkiR3SQGPI00iIi9Us+lHCMHA9BIMTSIiL1NzDTM4JJSB6UUYmkREXoRNP96NoUlE5EUqK8pRWVHJwPRSbAQiIvIC1euW2iAdDDEaKJVKuUuiBnCkSUQks+op2dJSEwAwML0YQ5OISEY11zBVKrXc5dBlMDSJiGTCph/fw9AkIpJJSUkxA9PHsBGIiEgmwcGh0GqDoNFo5S6FnMSRJhFRCxJCoKioAFarBQqFgoHpYxiaREQtpHoNs9xcBpvNJnc55AaGJhFRC2DTj39gaBIRtYCiwnwGph9gIxARUQvQ6YOh0wczMH0cQ5OIqAlyiixYta8QucU2BKkl3NA5BKlXBAOompI1l5UyLP0IQ5OIyA1n8irx6PvncexcBYS4eP8H3xYjKliJmddHIq2jDZUVFVBrtFCrebYff8DQJCJy0a85Fbh12RmUW6rSMilKjeRoNYrL7TiaXYH8UhvmfWjEb9eq8LcbWzEw/QhDk4jIBXa7HdNWZqHcIhAVrMSS2xPwf60uTr2Wltvw6Joz2HPcgrWHrOj/JxtSr5CxYPIods8SEblg46FiFJbZoVVLeH92cq3ABAC9VoEFY8PRLUkDAMjYmSdHmdRMGJpERC74774iAMB1/xcCQ+jFyTohBCyWSkiShLDwSDwxMg4A8Mv5ChhNVllqJc9jaBIRueDCHwE4vnuY477qExcU5Btht9sBAJ1aBSFcr4AQwNHsCllqJc9jaBIRueKPTlmVUqq6WeNMP+ERUVAoLv5aVUhV29jsot5uyDcxNImIXBCur/q1+cnh4kueGu9coQWFZVXnl+0Yr5GlVvI8hiYRkQvG/TEt++kPJSg1W2CzWhs8Nd6iT4wQAmgTrUZiJEPTXzA0iYhccGvfSOjUEkoq7PjLinNQB9cOTLvdjn9vN+LzoyUAgDsHRMpVKjUDfk+TiMgFaqWEZ0bo8eTmUvyeW4nrXzyFrslB6BivQWGpHXt/LYWpvKoZaMhVIRh9TbjMFZMnMTSJiJxUvYbZLUHg31Ni8OzmAuSarDj0uxmHfjc7ttOoJEzoGY5HboyRsVpqDgxNIiIn1G366RcfhM8ejcC+X0vx3v4iFJTZoFVJ6N1Oh9v6R0Kj4uqXP2JoEhE5oaKivMEu2T4dg9GnY7CMlVFLYmgSEV2CEAKSJCEoSAd1TByUSv7aDGScPyAiakTVlGw+yspKAYCBSQxNIqKGXFzDLIdSqZS7HPISDE0iojoudaYfCmwMTSKiOkpMxQxMahAn6ImI6ggOCYU2KAgajVbuUsjLcKRJRISqKdniogLYrFYoFAoGJjWIoUlEAa96DdNsLoPNZpO7HPJiDE0iCmh1m340Wo4wqXEMTSIKWEIIFBXms+mHnMZGICIKWJIkIUinh04fzMAkp3CkSUQBRwiBsrJSCCEQFKRjYJLTONIkooBScw1To9FApVLLXRL5EI40iShg1G36YWCSqxiaRBQQeGo88gSGJhEFDIWkYGBSk3BNk4j8mhACVqsVarUa4RFRcpdDPo4jTSLyW9VTsoUFRgi7Xe5yyA8wNInIL9VcwwwLj4Sk4K87ajp+iojI77Dph5qLbGuaK1e+iUMHDyA3NxcvvfQy2qak1Nvmp59+xPx585CQkOC4b968+Tw3JBFdks1mhdViZWCSx8kWmr1798bo0WPw9NwnL7ldQkICFr38SgtVRUS+TAgBIQRUKjUMMXGQJEnuksjPyBaanTtf6fF9WiwWWCwWx22z2ezx1yAi7ySEgN1uR1FhPiIioxmY1Cy8/isnOTnn8dijD0OhUGDQdddj6ND0RrfdtGkj1q97v979ecYLKNPpmlSHzWaFMTenSfsIRDxu7uFxc011YAICFkslj52L+HkDDDFxTm3n1aGZktIOS5cuhz44GHl5eVgw/wWEhoaib99+DW4/duw4jBgx0nHbbDZj5ozpiDbEQq/XN6kWY26O0weVLuJxcw+Pm/NqNv0oFArExLaSuySfw8+b87y6e1av10MfHAwAiI6ORr/+qTh29Gij26vV6qrn/PGja+Lokoi8X3m52dElK0le/SuN/IBXf8IKCgr+mHKpGjV+e+hgg122RBR4hBAAgKAgHaINceySpRYh2/Ts8mVL8e23h1BYWIh5855HUJAOry3OwNLXl6B79x7o3qMH9mfuw/bt26BUKmGz2dC7T19cd931cpVMRF6iako2H0FBOuh0eqhUXr3SRH5Etk/a9LtnNHj/jJmzHH9Ov3EY0m8c1lIlEZEPqLmGqdcHy10OBRivnp4lIqqJZ/ohuTE0ichnmExFDEySFRcCiMhnBAeHIihIB42Gp9Ikebg80vzwww9RWFjYDKUQEdUnhEBxcSFsNhuUSiUDk2Tlcmi+/PIr6NixEwYMGIinnpqL7du3o6SkpDlqI6IAV72GaS4rhc1mlbscItdD86uvvsT//vc/PProI6isrMAzzzyLlJR2GDJkaHPUR0QBqm7TD0eY5A3cWtOMjIxAp06dcO7ceZw/n4Ps7GzHSQiIiJpKCIGiwnw2/ZDXcTk077prKvbu3Yvo6CgMHDgQN988GYsXv4awsLDmqI+IApAkSVUnLtAHMzDJq7gcml988QXCwsKQljYYqamp6Nu3T5NPhk5EBFSNMMvNZQjS6RGk4+8V8j4ur2meOHEcq1a9i7i4WKxY8Qa6dPkzhg5Nx/z585ujPiIKENVrmFWdsmz6Ie/k1prmVVddhbZt26JDhw5ISUnB6tWrceDAATzxxBOero+IAkDdph+VSi13SUQNcjk0n33279izZw9++OEHdOzYAampqcjIyED//qnNUR8R+TmeGo98icuhWVxcjNmzZyM1tT8MBkNz1EREAUaCxMAkn+ByaP7jH684/pyXl4fo6GiPFkREgUEIAZvNCpVKjYhI/h4h3+ByI5DZbMYDDzyI+PhW6NChI+LjW+GBBx5EaWlpc9RHRH6oekq2IN/ouJg0kS9wOTSfeOJJ/Pbbr9iyZTN++eUYPvxwC44fP46nnprbHPURkZ+puYYZFh4JSZLkLonIaS5Pz3766af4+uuvERUVCQCIjY3F22+/hT59+uKf//yHxwskIv/Bph/ydS6PNIUQUChq/89QkhScYiGiy7JZrbBaLAxM8lkuh+bQoUNx222349tvv4PRaMShQ9/ijjvuQHp6enPUR0R+QAgBIQRUajUMMfEMTPJZLofm/PnzkJTUGunp6ejYsROGDRuGxMQEzJv3QnPUR0Q+rnpKtqioAAC4hkk+zeU1zZCQEGRkZGDx4sUwGo0wGAz8R0BEDaq7hknk69w6jR5Q9b/FmJgYT9ZCRH6ETT/kj5wKzeTkNk6NJk+d+r2p9RCRnzCbyxiY5HecCs3//nd1c9dBRH5CCAFJkqDT6aHRaHjydfIrToXmM888i507dwAAFi5ciDlz5jRrUUTkm4QQKCrMr7oeZpCOgUl+x6nu2V9//RVWa9X17RYvzmjWgojIN1WvYVZUlLM5kPyWUyPN1NT+6N8/Fe3bt4fZbMYtt/ylwe1Wr17l0eKIyDew6YcChVOhuXLlSmzevBmnTp3C9u3b0aXLVc1dFxH5EFNxEQOTAoJToanVajFx4kQAQGFhEdc0iaiW4JBQBAXpoNFq5S6FqFm5fEYgnvmHiICqKVlTcSHsdhuUSiUDkwKCy6FJRFS9hllWVupoEiQKBAxNInJJ3aYfjYYjTAocDE0ichq7ZCnQuRyab7yxosH777//gabWQkReTpIkBGl1DEwKWC6HZkZGBjZv3lzrvr/97SH89NNPHiuKiLyLEALl5jIAgE4fzMCkgOXyVU7Wr1+HUaNGIzo6Gv3798ejjz6Kb7/9Fps3f9AM5RGR3GpOyarVGihVbl8cicjnufzp79ChA9599x3ccstfMGBAKo4d+wVbtmxGeHh4c9RHRDKqu4bJwKRA59S/gB9//LHWba1Wi7vvvhtLly7FihVvICsrC1lZWbjqKp4piMhfsOmHqD6nQrN//1RIkgQhRL3HRo4cBaCqQaCgIN+z1RGRfIQABBiYRDU4FZqFhQXNXQcReQkhBGw2G1QqFSIio3nFEqIa+D1NInKonpItyDc6LiZNRBe5vKqfnZ2NefPm4fvvv4fJVFLrscOHf/BYYUTUsuquYTIwiepzOTSnT58OnU6PBx54AHq9vjlqIqIWxqYfIue4HJrff/8DTpw4Do1G0xz1EJEMrFYLLJUWBibRZbi8pvl///d/yMnJaY5aiKiFCSEghIBarYEhJo6BSXQZLo80R44ciZtvvhlTp05DbGxMrceGDRvmscKIqHlVT8kqFEqEh0dCoWBfINHluByaK1ZUnbD9lVdeqXW/JEkMTSIfUXcNk4ic43JoHjlyuDnqIKIWwqYfIvdxPoYowJjNpQxMIjc5NdIcOjQd27ZtBXDxlHoN2b17l+cqI6JmodMFQ6PWQqVWy10Kkc9xKjSnTr3L8edZs2Y2WzFE1DyEECgqzHdcC5OBSeQep0JzwoQJjj9PmTKl2YohIs+ruYap0wfLXQ6RT3MqND/55BOndsbuWSLvwqYfIs9yKjQfe2zOZbfhV06IvI+puJCBSeRBToUmv2ZC5Jv0waEICtJDo9XKXQqRX+BXToj8jBACJlMR7HY7VCoVA5PIgxiaRH6keg2zrLQEVqtF7nKI/A5Dk8hP1G360Wg4wiTyNIYmkR9glyxRy3D53LNE5H0kSYJWEwS9PoSBSdSMnArNS506ryaeRo+oZQkhUFFRjqAgHfTBIXKXQ+T3nApNnjqPyPvUnJJVx8RBqeTEEVFzc+pfGU+dR+Rd6q5hMjCJWoZbjUCrVq3CqFGj0bdvXwDAnj17sHHjJo8WRkQNY9MPkXxcDs1Fi17GkiVLMH78eGRlZQEA4uPj8e9//9vjxRFRfUIICCEYmEQycDk033nnHaxbtw63334bgKrmoHbt2uHkyZOero2IahBCwGazQqFQIDLSwMAkkoHLoVlWVob4+HgAcHTUWiwWaHmqLqJmUz0lW5BvhBDCqW52IvI8l0OzR4/uWLFiRa373n13FXr16uWxoojoopprmKFhEQxMIhm53HK3YMFCjBo1CqtX/xelpaUYPHgILly4gM2bP2iG8ogCG5t+iLyLy6GZktIW33yzH9u2bcfp06eRmJiI9PShCA527YrwK1e+iUMHDyA3NxcvvfQy2qakNLjd5zt34IMPNkEIgSuv6oKpU6dBpWJ7PQUGq8UCS2UlA5PIS7iVPjqdDmPGjG7SC/fu3RujR4/B03OfbHSbCzk5WLt2DV58cRHCIyLw0osLsWPHZ0hPv7FJr03k7ao7ZNUaDQwx8VAoeJpoIm/gVGjOnj3bqZ1lZGQ4/cKdO1952W0yM/fh2u49EBEZCQAYPGQINm3cyNAkvyaEgN1uh8lUhLCwCAYmkRdx6l9jWFiY40eSFFi3bj1yci5Ao9HiwoVcrF+/AQqF0uPFGY1GxMTEOG7HxsTCaDQ2ur3FYkFZWZnjx2w2e7wmouZUvYYJCE7HEnkhp0aaCxYscPz51ltvwzvvvI309HTHfdu2bcO7767yfHUu2rRpI9ave7/e/XnGCyjT6Zq0b5vNCmNuTpP2EYh43JxXPcIEBADAVFwEE4rkLcrH8PPmHh43wBAT59R2Lq9pfvHFF3j77bdq3ZeWloa77prq6q4uy2Aw4HzOxb/IC7kXYDAYGt1+7NhxGDFipOO22WzGzBnTEW2IhV6vb1Itxtwcpw8qXcTj5rzSUhNKTMWIiIyGqbiIx80N/Ly5h8fNeS4vliQnJ9UbVa5evRpJSUkeK6par969cejgARQWFEAIgc+2b0e/fv0a3V6tVkOv1zt+dE0cXRK1JL0+BFHRMZyWJfJiLo80Fy1ahJtvnoLXX38dSUlJOHPmDLKzs/Hee/91aT/Lly3Ft98eQmFhIebNex5BQTq8tjgDS19fgu7de6B7jx6Ii4vHhImTMPePDtvOna9E2uAhrpZM5LWEECgqKoBeHwyNRgu1WiN3SUR0CS6HZr9+/XD48A/YunUrzp/PQatW8RgyZCgiIyNc2s/0u2c0eP+MmbNq3U5LG4y0tMGulknk9WqeuECna9ryARG1DLe+pxkREYHJkycjLy8P0dHRnq6JyO/xTD9EvsmtE7bff/8DiI9vhQ4dOiI+vhUeeOBBlJaWNkd9RH6puLiQgUnkg1wOzSeffArHj/+GLVs245dfjuHDD7fg+PHjeOqpuc1RH5FfCg4OYWAS+SCXp2c//fRTfP3114iKqjpLT2xsLN5++y306dMX//znPzxeIJG/EEKgtMSE4OAQqFRqqFRquUsiIhe5PNIUQkChqH1pIklSQAjhsaKI/E31GmZpqQkWq0XucojITS6H5tChQ3Hbbbfj22+/g9FoxKFD3+KOO+6odYYgIrqobtOPRsMLthP5KpdDc/78eUhKao309HR07NgJw4YNQ2JiAubNe6E56iPyaeySJfIvLq9phoSEICMjA4sXL4bRaITBYOCV5IkaIUkSNGot9PoQBiaRH3D7as42mw1arRYmk8lxX1hYmEeKIvJ1QghUVlZAqw1CcEio3OUQkYe4HJoHDhzAAw88gKNHjzmaf4QQkCQJBQX5Hi+QyNc4pmQrK2EwxEGp9Pxl84hIHi6H5owZM3HTTeOxcuVKnhCdqI66a5gMTCL/4nJo5ubmYs6cOVzHJKqDTT9E/s/l7tkJEybgk08+aY5aiHyaEALCbmdgEvkxl0eaTz31FNLS0vCvf/0bMTExtR5bvXpVI88i8l9CCNjtdiiVSkRGxXAWhsiPuRya06dPh0ajQe/evaHXc02TAlv1lKzdZkNUdCwDk8jPuRyae/fuxS+/HENoKNvoKbDVXcNkYBL5P5fXNK+44gqUlJQ0Ry1EPoNNP0SByeWR5siRIzFx4iTcddddiI2tvaY5bNgwjxVG5M0slkpYKisZmEQBxuXQ/M9//gMAeOWVV2rdL0kSQ5P8XvUJPTQaLQwxcVAo+D1MokDicmgeOXK4Oeog8nrVU7IqlRqhoeEMTKIA5PKaJlEgqrmGyUt7EQUuhibRZbDph4iqMTSJLqOstISBSUQAmnBpMKJAoQ8OgUarhVqtkbsUIpIZR5pEDaiaks2HpbISkiQxMIkIAEOTqJ7qNcyKcjPswi53OUTkRRiaRDWw6YeILoWhSVRDcVEBA5OIGsVGIKIa9MEhCNLpGZhE1CCONCngCSFQUlIMIQTUag0Dk4gaxdCkgFa9hllaYoLVapG7HCLycgxNClh1m374tRIiuhyGJgUkdskSkTvYCEQBS63SQK8PYWASkdMYmhRQhBCwVFZCo9UiJDRM7nKIyMdwepYCRvWUbGFhHux2m9zlEJEPYmhSQKi5hhkeEcULSBORWxia5PfY9ENEnsLQJL9nt9tht9kZmETUZGwEIr8lhIAQdiiVSkRFx0CSJLlLIiIfx5Em+aXqKdmC/DwIIRiYROQRDE3yOzXXMENCwxiYROQxDE3yK2z6IaLmxNAkv1JZWQFLZSUDk4iaBRuByC8IIQAAWm0QDIY4KJT8HiYReR5HmuTzHJf3KjUBAAOTiJoNQ5N8Ws01TF7ai4iaG0OTfBabfoiopTE0yWeVlpgYmETUotgIRD4rOCQUWm0Q1BpOyxJRy+BIk3yKEAJFhfmwWCyQJImBSUQtiqFJPqN6DbO83MzrYRKRLBia5BPY9ENE3oChST6hqKiAgUlEsmMjEPkEvT4YOp2egUlEsuJIk7yWEAKlJSYIIaDRaBmYRCQ7hiZ5peo1zJKSYlitFrnLISICwNAkL1S36YenxyMib8HQJK/CLlki8mZsBCKvo1KqoI8MYWASkddhaJJXEELAYqmERqNFaFiE3OUQETWI07Mku+op2cKCPNjtdrnLISJqFEOTZFVzDTM8IgoKBT+SROS9+BuKZMOmHyLyNQxNko3dbofNamNgEpHPYCMQtTghBIQQUCqViDbEQpIkuUsiInIKR5rUomo2/QghGJhE5FMYmtRiaq5hBoeEMjCJyOcwNKlFsOmHiPyBbGua585lI2PxYphMxdDr9Zg1+x4kJSXX2uann37E/HnzkJCQ4Lhv3rz50Gi1LV0uNVFlRTkqKyoZmETk02QLzeXLliEtLQ2Drrsemfv2YUnGYixY+FK97RISErDo5VdkqJA8QQgBANAG6WCI0UCpVMpcERGR+2SZni0qKsKJE8eROmAgAKBX794wGvNw/tw5OcqhZiKEgN1uR2mpCQAYmETk82QZaeYZjYiIiHT8EpUkCQaDAUajEfGtWtXaNifnPB579GEoFAoMuu56DB2a3uh+LRYLLJaL1140m83N8wbosqrXMAEBlUotdzlERB7h1d/TTElph6VLl0MfHIy8vDwsmP8CQkND0bdvvwa337RpI9ave7/e/XnGCyjT6ZpUi81mhTE3p0n7CBTVI0ygamrWVFwEE4rkLcrH8PPmHh439/C4AYaYOKe2kyU0ow0GFBYWwGazQalUQggBo9EIg8FQazu9Xn/xOdHR6Nc/FceOHm00NMeOHYcRI0Y6bpvNZsycMR3Rhtha+3KHMTfH6YMa6EymIpSVliAiMhqm4iIeNzfw8+YeHjf38Lg5T5Y1zfDwcKSktMPuXV8BAPZnZiI6Orre1GxBQYHjqhdmsxnfHjqItikpje5XrVZDr9c7fnRNHF2Se4KDQxEZZWCXLBH5HdmmZ6dPvxsZGYuxadNG6HR6zJo1GwCw9PUl6N69B7r36IH9mfuwffs2KJVK2Gw29O7TF9ddd71cJdMlCCFQXFyI4OAQqFRqaDT8WhAR+R/ZQjMhMRHz5i+od/+MmbMcf06/cRjSbxzWkmWRG2qeuCAoSMfGHyLyWzwjEDUJz/RDRIGEoUlNUlSYz8AkooDh1V85Ie+n0wdDpw9mYBJRQOBIk1wmhEBZaQmEENBqgxiYRBQwGJrkkuo1TJOpCFarVe5yiIhaFEOTnFa36UetZpcsEQUWhiY5hV2yREQMTXKBUqFkYBJRQGP3LF2SEAJWqwVqtQZh4ZFyl0NEJCuONKlR1VOyBflGxzmAiYgCGUOTGlRzDTM8IgoKBT8qRET8TUj1sOmHiKhhDE2qx26zwWa1MjCJiOpgIxA5CCEghIBSpUK0IQ6SJMldEhGRV+FIkwBcnJItKsyHEIKBSUTUAIYm1VrD1AeHMDCJiBrB0AxwbPohInIeQzPAVVSUMzCJiJzERqAAVb1uGRSkgzomDkolPwpERJfDkWYAqp6SLSsrBQAGJhGRkxiaAabmGqZSqZS7HCIin8LQDCBs+iEiahqGZgApMRUzMImImoCLWQEkOCQU2qAgaDRauUshIvJJHGn6OSEEiosKYLNaoVAoGJhERE3A0PRj1WuYZnMZbDab3OUQEfk8hqafqtv0o9FyhElE1FQMTT8khEBRYT6bfoiIPIyNQH5IkiTodHro9MEMTCIiD+JI048IIVBWVgohBLRBOgYmEZGHcaTpJ2quYWo0GqhUarlLIiLyOxxp+oG6TT8MTCKi5sHQ9HE8NR4RUcthaPoBhaRgYBIRtQCuafooIQSsVivUajXCI6LkLoeIKCBwpOmDqqdkCwqMEHa73OUQEQUMhqaPqbmGGR4eCUnBv0IiopbC37g+hE0/RETyYmj6EJvNBqvFysAkIpIJG4F8gBACAKBSqWCIiYMkSTJXREQUmDjS9HLVU7KFhfkAwMAkIpIRQ9OL1VzD1OuD5S6HiCjgMTS9FJt+iIi8D0PTS5WXmxmYREReho1AXkYIAUmSEBSkg1qtgUrFvyIiIm/BkaYXqZqSzYfZXAZJkhiYRERehqHpJS6uYZZDwbP8EBF5Jf529gJs+iEi8g0MTS9gMhUxMImIfAAXzbxASHAogoJ00Gi0cpdCRESXwJGmTIQQKC4uhM1mg0KpZGASEfkAhqYMqtcwzWWlsFmtcpdDREROYmi2sLpNPxotR5hERL6CodmCqr+HyaYfIiLfxEagFiRJEnRBOuj1wQxMIiIfxNBsAUIIlJvLEKTTI0inl7scIiJyE0OzmdVcw1RrNFCp1HKXREREbuKaZjOq2/TDwCQi8m0MzWbCU+MREfkfhmYzkiAxMImI/AjXND1MCAGbzQqVSo2IyGi5yyEiIg/iSNODqqdkC/KNEELIXQ4REXkYQ9NDaq5hhoVHQpIkuUsiIiIPY2h6AJt+iIgCA0PTA2xWK6wWCwOTiMjPsRGoCarXLVVqNQwx8ZySJSLycxxpuql6SraoqAAAGJhERAGAoemGmmuYOp5LlogoYDA0XcSmHyKiwCXbmua5c9nIWLwYJlMx9Ho9Zs2+B0lJyfW2+3znDnzwwSYIIXDlVV0wdeo0qFTyLcWWm8sYmEREAUq2kebyZcuQlpaGf/17MUaPHoslGYvrbXMhJwdr167Bc8+9gH+/loGiwkLs2PGZDNVebPoJ0ukRbYhlYBIRBSBZQrOoqAgnThxH6oCBAIBevXvDaMzD+XPnam2XmbkP13bvgYjIqpMFDB4yBHv37GnxeoUQsNvtKC83Q5IkXq2EiChAyTLPmWc0IiIiEkqlEkBV56nBYIDRaER8q1aO7YxGI2JiYhy3Y2NiYTQaG92vxWKBxWJx3DabzU2utXoNExDskCUiCnB+9T3NTZs2Yv269+vdn2e8gDKdzuX9VY8wgaqpWVNxEUwoamqZAcVms8KYmyN3GT6Hx809PG7u4XEDDDFxTm0nS2hGGwwoLCyAzWaDUqmEEAJGoxEGg6HWdgaDAedzLv5FXsi9UG+bmsaOHYcRI0Y6bpvNZsycMR3Rhljo9a5/NaS4qBBmcykiIqNhKi5y+qDSRcbcHB43N/C4uYfHzT08bs6TZU0zPDwcKSntsHvXVwCA/ZmZiI6OrjU1C1StdR46eACFBQUQQuCz7dvRr1+/RverVquh1+sdPzo3Rpc1BYeEIjLSwKYfIiICIOP07PTpdyMjYzE2bdoInU6PWbNmAwCWvr4E3bv3QPcePRAXF48JEydh7twnAQCdO1+JtMFDmrUuIQRKTEUIDgmFUql0rLsSERHJFpoJiYmYN39BvftnzJxV63Za2mCkpQ1ukZpqnrhAG6SDRsPAJCKii3hGoD/UPdOPRqOVuyQiIvIyDE1UB2Y+z/RDRESX5FdfOXGXJEkI0gZBrw9mYBIRUaMCeqQphIDZXAYA0DEwiYjoMgJ2pFlzDVOj1kAp40ngiYjINwTkSLNu0w8Dk4iInBFwocnrYRIRkbsCLjQhBCDAwCQiIpcFzLykEAI2mw0qlQoRkdG8YgkREbksIEaa1VOyBflGCMFLfBERkXsCYqRZVJQPpULBESYRETVJQIw02fRDRESe4NcjTSGqLh6t1epgs9lRVlbm9r7MZnOTnh+oeNzcw+PmHh439/C4VdHpdJedjfTr0CwvLwcA3H//fTJXQkRE3u6tt9+FXq+/5DaS1WoRLVRPi7Pb7SgoKEBQUFCT1jLNZjNmzpiO15cub/KFrQMJj5t7eNzcw+PmHh63iwJ+pKlQKBAdHe2x/el0usv+L4Tq43FzD4+be3jc3MPj5pyAaAQiIiLyBIYmERGRkxiaTlCr1bhpwkSo1Wq5S/EpPG7u4XFzD4+be3jcXOPXjUBERESexJEmERGRkxiaRERETmJoEhEROcmvv6fpqnPnspGxeDFMpmLo9XrMmn0PkpKS6233+c4d+OCDTRBC4MqrumDq1GlQqQL3UDpz3H766UfMnzcPCQkJjvvmzZsPjVbb0uV6jZUr38ShgweQm5uLl156GW1TUhrcjp+32pw5bvy81VZZWYlXX/0HzmZlQaPRICwsHNOmTUd8q1b1tj106CDefedt2O12JCe3wazZ9/D7mzUE7r+8BixftgxpaWkYdN31yNy3D0syFmPBwpdqbXMhJwdr167Biy8uQnhEBF56cSF27PgM6ek3ylS1/Jw5bgCQkJCARS+/IkOF3ql3794YPXoMnp77ZKPb8PNWnzPHDeDnra60tMHo1u0aSJKErZ9+gqVLX8ezf3+u1jblZjOWvr4Ez/79OSQmtsabK97AhvXrcOttt8tUtffh9OwfioqKcOLEcaQOGAgA6NW7N4zGPJw/d67WdpmZ+3Bt9x6IiIyEJEkYPGQI9u7ZI0fJXsHZ40b1de585WXPWMXPW33OHDeqTaPR4JprrnWcIq5jp07Izb1Qb7vvvv8ObdumIDGxNQBg6NB07N0b2J+3ujjS/EOe0YiIiEgolUoAgCRJMBgMMBqNtaYwjEYjYmJiHLdjY2JhNBpbvF5v4exxA4CcnPN47NGHoVAoMOi66zF0aLocJfsUft7cx89b4z75+GN0796j3v11P28xsbEoKCiEzWZz/BsPdAxNahEpKe2wdOly6IODkZeXhwXzX0BoaCj69u0nd2nkh/h5a9zGjRtw/vx5PP3Ms3KX4pM4PfuHaIMBhYUFsNlsAKquxWk0GmEwGGptZzAYkJub67h9IfdCvW0CibPHTa/XQx8cXPWc6Gj065+KY0ePtni9voafN/fw89awLVs245v9+/HEk09B20BTVN3PW+6FC4iMjOAoswaG5h/Cw8ORktIOu3d9BQDYn5mJ6OjoelOMvXr3xqGDB1BYUAAhBD7bvh39+gXu/16dPW4FBQWw2+0Aqi5F9O2hg412i9JF/Ly5h5+3+j76cAv27tmDp+Y+jeA//kNRV9eu3XDy5AmcPZsFANi2bSv69uvfkmV6PZ5Gr4bss2eRkbEYJSUm6HR6zJo1G8lt2mDp60vQvXsPdO9RtQawY8dn2PzBJgBVTQnTpt8d0F8BcOa4bf30E2zfvg1KpRI2mw29+/TFhAkTm3SdU1+3fNlSfPvtIRQWFiI0NBRBQTq8tjiDn7fLcOa48fNWW15eHmbOmI64uDgEBVVdM1OtVmP+goVYu+Y9REZFYciQoQCAgwcOYNWqd2Cz2ZGUnIR7Zt/rGLUTQ5OIiMhpnJ4lIiJyEkOTiIjISQxNIiIiJzE0iYiInMTQJCIichJDk4iIyEkMTSIiIicxNMlndOnSBR999JFH9jV+/E14440VHtlXc7LZbOjbty9+/vlnAMCCBQswZcoUmau6tFOnTqF79x6oqKhodJvVq1ejf/+LZ5rp1as3tm7d6tT+XT0Gn376Kbp06YKEhER89NFHGD58OJYsWeL084lqYmhSQNqwYT2mTZsKANi9ezeSk+tfbLyu1157Dddccy0SE1ujXbv2GDVqNE6dOgWgfghUmzlzJubMmQOgKkzCwyOQkJBY66f68Ya89957aNeuPTp37uzO22yy6poLCwudfk6bNm3Qs2cPrFy50unn7N+fifT05rkKyeOPP4Enn3wS2dlnMWLEiGZ5DQocgXsuLiIXrF27FsuXL8fatWvRuXNnFBYW4vPPv3DrtGw///wTIiIinNr2jTdW4PHHGw/Vy7FYLFCr1W4/310333wz7rvvfsycObPFX7uuU6dOyfafDvI/HGmSz1q7di169OiJ5ORkDB2aju+//97xWGFhIW677XYkJyeje/ceWLZsGcLDIxyPV0/R5efn46abJqCoqNgx8vv666/rvdaBAwcxYMBAxy/fiIgIjBs31qkRqrvOnTuHw4cP1ztBu9Vqw+zZ96B16yR063YNPvzwQ8djM2fOxOzZ9+D22/+K1q2TsHLlSlgsFsybNw9XX90VbdumYPLkyThX4yLh4eERePPNlejduw9at07C5MmTUVRUBAC4/vobAFSd8zYhIRHvv/8+AGDz5s3o2rUbkpKSce+992HixElYsGCBY5+9e/dGdnY2fvnlF6fea82p9+pR+0svvYT27TugQ4eOl5xOfe6559GvXz+cP3++1v35+flISEiE3W7HkCFDkZCQ2OCU8c6dn6N//1QkJSUjNXUAvvjiSwBATk4ODIYYlJSUAIDjM/S///0PQNW0b58+fZ16f+Q/GJrkk/bu3Yu//e0hvPrqP3H8+HGMHj0K48ff5Phl/+ijj6KsrAxHjhzBRx99iDVr1ja4n6ioKKxfvw7h4WHIzj6L7Oyz6Nu3/i/C3r174YMPPsCiRS8jMzMT5eXlzfr+AODIkSNISEhAaGhorft37NiBa6+9Fr//fhLz58/DXXdNxYkTJx2Pb9iwAbfeeitOnz6FW2+9Fc899zwyM/dj27at+N//fkH79h1w55131trnpk2b8OGHW/Djj0dw9my2I6Q+/3wngKrRcXb2WUycOBG//fYbpk6dhgUL5uPkyRPo2rUrduzYUWt/arUa7dq1w5EjR9x670ePHoNOp8OxY0fxn/+sxNy5T9d6jwBgtVoxe/Y92L8/E5988gni4+NrPR4VFYXs7LMAgO3btyE7+2y9y2EdP34CU6ZMwaOPPoKTJ0/goYcews0334zff/8dcXFxaNeuHfbt2wcA2LVrF1JSUrBr127H7QEDUt16f+S7GJrkk9asWYuJEyeiX79+UKvVmDVrFiIiIrB9+3bYbDZs3LgJTzzxBMLDwxEfH4/77ruvSa930003ISMjA9988w0mTpyIlJR2uPfe+1BaWuryvrp0+TOSk5MdP6tWrWpwu+qreNTVoUMH3HnnHVCpVLjxxhuRmpqKDRvWOx6//vrrkJZ2AxQKBXQ6Hd58803Mnz8P8fHx0Gg0mDv3KWRm7kdWVpbjOQ88cD9iYmIQERGBUaNG1Rq117Vx40YMHDgQN954I1QqFe666060b9++3nahoaEoKCh0/sDUEB0djXvvvRdqtRqpqalITk7GkSOHHY+XlZlxyy23oLi4GBs3bkR4eLhbr7Nx40b0798fo0aNgkqlwpgxo9G7d2+sX78BAJCamopdu3bDbrcjM3M/Hn74IezeXTM0B7j1uuS7uKZJPik7O7te402bNm1w9mw28vLyYLFY0Lp1ouOxpKTWTX7NMWNGY8yY0RBCIDMzE9OmTcPLL7+CZ555Gmq1GhaLtd5zLBYrVKraa4pHjhx2ak0zIiICJpOp3v1JSUn1bmdnX5xubd364nvNy8tDaWkphg0bjprLrxqNBmfPnnVsGxsb63gsOFjvmJJsyLlz5xusoS6TyYTIyIhG93MpsbExtW7XrenIkSMwmUz48ssvGryYsrOys7PrTbG3bdsW2dnZAKpC89VXX8UPPxxGmzZtMGzYcPz978/BaDTi2LFfeG3TAMSRJvmkhIQEnD59utZ9p0+fRmJiAqKjo6FWq5GVddbx2JkzWXV34aBQuPbPQJIk9OnTB6NHj8bPP/8EoDq4smG11g7OkydPok0b99Y9u3Tpguzs7HoBdubMmVq3s7KykJBw8aLfNd9PVFQU9Ho9du7cgdOnTzt+cnLOo1evXpetoaFj06pVfIM11GSxWHDixAl06dLlsq/hjl69euLllxdhzJixOHr0qNv7aexzlJCQAABITe3/xxT/RxgwYACioiIRHx+P5cuX46qrrnK6oYv8B0OTfNKkSROxbt06ZGZmwmq1YtmyZcjPz8fgwYOhVCoxduxYLFy4EEVFRcjJycHixYsb3VdsbCxMphLk5uY2us2qVavw8ccfO7568fPPP+Pjjz9Bz55VwXPttdciJiYGzz33PEpLS1FZWYlVq1bh2LFjGDZsmFvvsVWrVujSpQv27t1b6/7ffvsNb731NqxWK7Zt24Zdu3Zh3LhxDe5DoVDgzjvvwJNPPuUItvz8fGzYsNGpGgwGAxQKBU6e/N1x39ixY/HVV19h27ZtsFqteOutt/Hbb7/Vet7+/fvRqlUrXHHFFS68Y9fcdttteOaZpzFq1Gj8+OOPbu1j3Lhx2LNnDz7++GNYrVZs2bIFX3/9NcaPHw+gapr4iis6Yfny5UhNrVq/HDBgAF5/fSmnZgMUQ5N8UlV35Yu45557kZLSDhs2bMT69esd//N/6aWXoNFocNVVXTB8+AiMHTsGGo2mwX117NgRt956K3r27IXk5GRH40dN4eHhWLx4Ma6+uisSE1tjypRbcNNN4/HAA/cDqJruXLduHX777VdcfXVXdOzYCatXr8aGDeuRmJhYa1/VnajVP7fddnuj73PatKlYtWp1rfvS0tJw4MABtG2bgscem4Ply5c3uKZY7ZlnnkHPnj0wcuQoJCa2xsCBA/H55583un1NOp0Oc+Y8hptuugnJyclYt24dOnbsiGXLluHRRx9DSko7HDx4EGlpabWe9957axzfg21Ot9xyC55//jmMGTMWhw8fvvwT6mjfvh1WrXoXCxYsQNu2KXjppZewatUqpKS0dWyTmpqK8vJy9OnTGwAwaNBAFBcXY+BAhmYgkqxWi5C7CKLmtn79esybNx/fffet3KW4xGazITU1FW+++Sb+9Kc/yV1Oo6ZMmYIuXbrg8ccfx+nTpzF+/E3Ys2d3k9YbibwRG4HILx0/fhxFRcXo1q0rTpw4gUWLXsaYMWPkLstlSqWywe+NerPk5GQcOPCN3GUQNQuGJvml0tJSTJ8+HWfPnkVYWBhGjhyJRx55WO6yiMjHcXqWiIjISWwEIiIichJDk4iIyEkMTSIiIicxNImIiJzE0CQiInISQ5OIiMhJDE0iIiInMTSJiIic9P947qaIGHBfGwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 480x460 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(viz.plot_network_flows(net, final.link_flows))\n",
    "display(viz.plot_flow_scatter((\"logit SUE (brentq)\", ref_flows), {\"sue-msa\": final.link_flows}))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47d5dccc",
   "metadata": {},
   "source": [
    "## Takeaways & pointers\n",
    "\n",
    "- **Certified, not self-reported.** The SUE residual came from `Evaluator`; the\n",
    "  closed-form logit self-report was only diffed against it.\n",
    "- **Stochastic spreading.** Unlike a deterministic UE, flow uses both routes by the\n",
    "  logit dispersion θ — the fixed point recomputed here via `brentq`.\n",
    "- **Where next.** The Monte-Carlo probit variant: [`sue-probit-msa`](11-sue-probit-msa.ipynb); the deterministic-UE baseline it mirrors: [`msa`](02-msa.ipynb);\n",
    "  the ADR-001 SUE design in the [model compendium](../../docs/MODELS.md)."
   ]
  }
 ],
 "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.10.12"
  },
  "tabench": {
   "covers": [],
   "requires_extra": null,
   "track": "static",
   "unit": "sue-msa"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}