{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2d5f4bca",
   "metadata": {},
   "source": [
    "# `lp-so-dta` — Ziliaskopoulos's (2000) LP single-destination SO-DTA on CTM cells\n",
    "\n",
    "**What.** `lp-so-dta` formulates single-destination system-optimal DTA as ONE\n",
    "linear program by embedding Daganzo's Cell Transmission Model as linear cell-\n",
    "flow constraints: the CTM's `min(sending, receiving)` Godunov flux equality is\n",
    "relaxed to its four linear `<=` bounds (conservation stays an equality). Unlike\n",
    "`merchant-nemhauser`'s exit functions, cells carry FINITE storage, so this LP\n",
    "can represent queue spillback.\n",
    "\n",
    "**Why it is in the benchmark.** It is the strongest analytic-DTA anchor in the\n",
    "benchmark: a certifiable GLOBAL optimum via LP duality, reproducible with\n",
    "`scipy.optimize.linprog` (no external solver) — no other track in this family\n",
    "gives up so little to get an exact answer. See the\n",
    "[model compendium](../../docs/MODELS.md) (Ziliaskopoulos 2000) and\n",
    "[docs/design/adr-021-lp-so-dta.md](../../docs/design/adr-021-lp-so-dta.md)\n",
    "(P1).\n",
    "\n",
    "**Scope.** This notebook solves the LP on two built-in anchors — a control-\n",
    "free corridor and a diverge with a tiny-storage bottleneck cell (the\n",
    "spillback effect `merchant-nemhauser`'s exit functions cannot represent) —\n",
    "and certifies both.\n",
    "\n",
    "**Canon.** `[ziliaskopoulos2000linear]`, [docs/REFERENCES.md](../../docs/REFERENCES.md) / [docs/references.bib](../../docs/references.bib)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcfa3f00",
   "metadata": {},
   "source": [
    "## How this notebook is graded\n",
    "\n",
    "**A notebook never claims a number it does not compute in that cell.** Every\n",
    "scored quantity below is recomputed live by the P1 `CellSODTAEvaluator` from\n",
    "the emitted occupancy/flow arrays alone — conservation (with exogenous demand\n",
    "and the absorbing sink), the four CTM bound families, the storage envelope, and\n",
    "a weak-duality backstop against the harness's own freshly-solved LP optimum\n",
    "`Z*`. The solver's self-reported objective is never trusted\n",
    "([README](../../README.md), *Certified, not self-reported*)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "73d4cd2f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T13:48:57.848481Z",
     "iopub.status.busy": "2026-07-21T13:48:57.847844Z",
     "iopub.status.idle": "2026-07-21T13:48:59.907173Z",
     "shell.execute_reply": "2026-07-21T13:48:59.905835Z"
    }
   },
   "outputs": [],
   "source": [
    "# Setup. `lp-so-dta` 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\n",
    "# CI renders into the notebook); NEVER matplotlib.use(\"Agg\") in-kernel — it\n",
    "# silently suppresses inline figure capture.\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "\n",
    "from tabench import (\n",
    "    CellSODTAEvaluator,\n",
    "    solve_cell_so_dta,\n",
    "    zil_corridor_scenario,\n",
    "    zil_diverge_spillback_scenario,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbf02cde",
   "metadata": {},
   "source": [
    "## Anchor A: a control-free corridor\n",
    "\n",
    "`R -> A -> B(Q=1, N=2) -> S`, 6 vehicles queued at the source. The LP has no\n",
    "diverge to control (a single path) — no decision, just the physics.\n",
    "`CellSODTAScenario` is frozen and content-hashed (P2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d15f6d7c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T13:48:59.912156Z",
     "iopub.status.busy": "2026-07-21T13:48:59.911700Z",
     "iopub.status.idle": "2026-07-21T13:48:59.919630Z",
     "shell.execute_reply": "2026-07-21T13:48:59.918833Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "scenario      : zil-corridor\n",
      "content hash  : 80ca151d03f194f2…\n",
      "cells=4  sink=3  periods=10\n",
      "capacity=[inf 10.  1. inf]  storage=[inf 20.  2. inf]\n",
      "initial occupancy at source : 6.0\n"
     ]
    }
   ],
   "source": [
    "scenario_a = zil_corridor_scenario()\n",
    "print(f\"scenario      : {scenario_a.name}\")\n",
    "print(f\"content hash  : {scenario_a.content_hash()[:16]}…\")\n",
    "print(f\"cells={scenario_a.n_cells}  sink={scenario_a.sink}  periods={scenario_a.n_periods}\")\n",
    "print(f\"capacity={scenario_a.capacity}  storage={scenario_a.storage}\")\n",
    "print(f\"initial occupancy at source : {scenario_a.initial_occupancy[0]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5a0ce66",
   "metadata": {},
   "source": [
    "## Solve + certify (Anchor A)\n",
    "\n",
    "No `Budget`/`RngBundle`/`Trace` — `solve_cell_so_dta` builds the canonical LP\n",
    "(`cell_canonical_lp`) and solves it with HiGHS to global optimality, emitting a\n",
    "`CellTrajectory`. `CellSODTAEvaluator` independently re-solves the LP and\n",
    "certifies feasibility + optimality by pure arithmetic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2744dcab",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T13:48:59.924228Z",
     "iopub.status.busy": "2026-07-21T13:48:59.923566Z",
     "iopub.status.idle": "2026-07-21T13:48:59.941915Z",
     "shell.execute_reply": "2026-07-21T13:48:59.941174Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "feasible           : 1\n",
      "so_optimality_gap   : 0.000e+00\n",
      "total_cost          : 33.0000\n",
      "holding_max         : 1.0000  (Tier-B diagnostic, not an error)\n",
      "dual_gap            : 0.000e+00\n"
     ]
    }
   ],
   "source": [
    "trajectory_a = solve_cell_so_dta(scenario_a)\n",
    "evaluator_a = CellSODTAEvaluator(scenario_a)\n",
    "metrics_a = evaluator_a.certify(trajectory_a)\n",
    "print(f\"feasible           : {metrics_a['feasible']:.0f}\")\n",
    "print(f\"so_optimality_gap   : {metrics_a['so_optimality_gap']:.3e}\")\n",
    "print(f\"total_cost          : {metrics_a['total_cost']:.4f}\")\n",
    "print(f\"holding_max         : {metrics_a['holding_max']:.4f}  (Tier-B diagnostic, not an error)\")\n",
    "print(f\"dual_gap            : {metrics_a['dual_gap']:.3e}\")\n",
    "assert metrics_a[\"feasible\"] == 1.0\n",
    "assert metrics_a[\"so_optimality_gap\"] < 1e-6\n",
    "assert np.isclose(metrics_a[\"total_cost\"], 33.0, atol=1e-6)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37502d18",
   "metadata": {},
   "source": [
    "## Anchor B: a finite-storage effect exit functions cannot represent\n",
    "\n",
    "`R -> A -> {B(Q=1, N=1, tiny storage), C -> D(Q=2 each)} -> S`: 6 vehicles at\n",
    "`R` face a diverge between a short route through a ONE-vehicle-storage\n",
    "bottleneck cell `B` and a longer route `C -> D`. Certified below: `B` sits\n",
    "jam-full at `t=2` in the optimum (the spillback pair lemma binds)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ee0db070",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T13:48:59.943979Z",
     "iopub.status.busy": "2026-07-21T13:48:59.943693Z",
     "iopub.status.idle": "2026-07-21T13:48:59.958189Z",
     "shell.execute_reply": "2026-07-21T13:48:59.957596Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "feasible            : 1\n",
      "so_optimality_gap    : 0.000e+00\n",
      "total_cost           : 26.0000\n",
      "holding_max          : 0.0000\n",
      "cell B occupancy at t=2 : 1.0000  (storage cap: 1.0)\n",
      "relaxed (N_B=2) total_cost : 25.0000\n",
      "finite storage on cell B strictly costs 1 extra veh-interval (26 vs 25) —\n",
      "an effect the exit-function merchant-nemhauser model cannot represent.\n"
     ]
    }
   ],
   "source": [
    "scenario_b = zil_diverge_spillback_scenario()\n",
    "trajectory_b = solve_cell_so_dta(scenario_b)\n",
    "evaluator_b = CellSODTAEvaluator(scenario_b)\n",
    "metrics_b = evaluator_b.certify(trajectory_b)\n",
    "print(f\"feasible            : {metrics_b['feasible']:.0f}\")\n",
    "print(f\"so_optimality_gap    : {metrics_b['so_optimality_gap']:.3e}\")\n",
    "print(f\"total_cost           : {metrics_b['total_cost']:.4f}\")\n",
    "print(f\"holding_max          : {metrics_b['holding_max']:.4f}\")\n",
    "assert metrics_b[\"feasible\"] == 1.0\n",
    "assert metrics_b[\"so_optimality_gap\"] < 1e-6\n",
    "assert np.isclose(metrics_b[\"total_cost\"], 26.0, atol=1e-6)\n",
    "# Cell B (index 2) is jam-full (storage 1) at t=2 in every optimum.\n",
    "b_occ_at_2 = trajectory_b.occupancies[2, 2]\n",
    "print(f\"cell B occupancy at t=2 : {b_occ_at_2:.4f}  (storage cap: {scenario_b.storage[2]})\")\n",
    "assert np.isclose(b_occ_at_2, scenario_b.storage[2], atol=1e-6)\n",
    "\n",
    "# The distinctive point: relaxing B's storage from 1 to 2 strictly LOWERS the\n",
    "# optimum -- a genuine finite-storage effect, recomputed here by rebuilding the\n",
    "# scenario with N_B=2 and re-solving, not quoted.\n",
    "from dataclasses import replace\n",
    "relaxed_storage = list(scenario_b.storage)\n",
    "relaxed_storage[2] = 2.0\n",
    "relaxed = replace(scenario_b, name=\"zil-diverge-relaxed-storage\", storage=relaxed_storage)\n",
    "relaxed_traj = solve_cell_so_dta(relaxed)\n",
    "relaxed_metrics = CellSODTAEvaluator(relaxed).certify(relaxed_traj)\n",
    "print(f\"relaxed (N_B=2) total_cost : {relaxed_metrics['total_cost']:.4f}\")\n",
    "assert relaxed_metrics[\"feasible\"] == 1.0\n",
    "assert np.isclose(relaxed_metrics[\"total_cost\"], 25.0, atol=1e-6)\n",
    "assert relaxed_metrics[\"total_cost\"] < metrics_b[\"total_cost\"]\n",
    "print(\"finite storage on cell B strictly costs 1 extra veh-interval (26 vs 25) —\")\n",
    "print(\"an effect the exit-function merchant-nemhauser model cannot represent.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd126baa",
   "metadata": {},
   "source": [
    "## Visualize\n",
    "\n",
    "`tabench.viz` is a road-network flow visualizer; the LP artifact is a per-cell\n",
    "OCCUPANCY trajectory, not a static link-flow vector, so this notebook plots\n",
    "the certified occupancy trajectories directly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "734d6197",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T13:48:59.960281Z",
     "iopub.status.busy": "2026-07-21T13:48:59.959818Z",
     "iopub.status.idle": "2026-07-21T13:49:00.222239Z",
     "shell.execute_reply": "2026-07-21T13:49:00.221179Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6sAAAFeCAYAAABjIokGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABkSklEQVR4nO3deVwU9R8G8GeBZVluVOQIPPBCBEXBE88UbzzKIzXFo8wrNTrUn5lHKWpqmlepeZRpKB6VmYpXpqbiQUre5oUmnlxyL9/fH8TmuqiwrMwsPO/Xa1+xs7Mzzw7GZz4735lRCCEEiIiIiIiIiGTETOoARERERERERE9js0pERERERESyw2aViIiIiIiIZIfNKhEREREREckOm1UiIiIiIiKSHTarREREREREJDtsVomIiIiIiEh22KwSERERERGR7LBZJSIiIiIiItlhs0okgYEDB6JSpUovnO/atWtQKBRYvXr1S8+UZ8OGDShTpgxSUlKKbZ1yk5WVBU9PTyxZskTqKERUCkyZMgUKhUJnWqVKlTBw4EBpApVQ+dVehUKBKVOmaJ+vXr0aCoUC165d006rVKkSOnfuXDwh/5VX/+fMmWPwMm7evAkrKyscOnTIiMlMT6NGjfDRRx9JHYMMxGaViLQ0Gg0mT56Md999F7a2tpLlmDFjBrZu3VqgedPS0jBkyBD4+vrCwcEBtra2qFOnDhYsWICsrCydeQ8cOIAuXbrA09MTVlZWcHV1Rfv27fUKuVKpRFhYGKZPn4709HRjfSwiIqJiM23aNDRs2BBBQUGSZVi3bh3mz59foHlTU1OxePFitG3bFm5ubrCzs0PdunWxdOlSaDSafN9z5coV9O3bF+XLl4darUa1atUwceJEnXnGjRuHxYsX486dO0X9OCQBC6kDEJVGy5cvR05OjtQx9Pz888+4cOEChg4dKmmOGTNmoEePHujWrdsL501LS8Nff/2Fjh07olKlSjAzM8Phw4fx3nvv4ejRo1i3bp123osXL8LMzAzDhg2Dq6srHj16hLVr16J58+b45Zdf0L59e+28gwYNwvjx47Fu3ToMHjz4ZXxMIqJnunDhAszMeEzBmORae1+Ge/fuYc2aNVizZo2kOdatW4fY2FiMHTv2hfP+/fffePfdd9G6dWuEhYXB3t4eO3fuxIgRI3DkyBG9zxITE4OWLVvilVdewfvvv4+yZcvixo0buHnzps58Xbt2hb29PZYsWYJp06YZ8+NRMWCzSmRk6enpsLS0zHcn4/Hjx7CxsYFSqZQg2X/rf5ZVq1YhKCgIr7zySjGmKpoyZcrgyJEjOtOGDRsGBwcHLFq0CPPmzYOrqysA4K233sJbb72lM++IESPg5eWF+fPn6zSrjo6OaNu2LVavXs1mlYiKnUqlKvZ1Pq9+Gdu1a9dQuXJl7Nu3Dy1btnzp6wMgWe2Vwtq1a2FhYYGQkBCpoxSYq6srzpw5g1q1ammnvfPOOxg8eDBWrVqFSZMmoWrVqgCAnJwc9O/fH97e3ti3bx/UavUzl2tmZoYePXrg22+/xdSpU/WG3JO88Ss7KpVu3bqFIUOGwN3dHSqVCpUrV8bw4cORmZmpnefvv/9Gz549UaZMGVhbW6NRo0b45ZdfdJazf/9+KBQK/PDDD/j444/xyiuvwNraGklJSRg4cCBsbW1x5coVdOzYEXZ2dujXrx+A/M+bSUhIwMCBA+Hg4ABHR0eEhoYiISEh3/x79+5Fs2bNYGNjA0dHR3Tt2hXnzp3TmSfvHKizZ8+ib9++cHJyQtOmTZ+5TdLT07Fjxw60adMm39fXrl2LBg0awNraGk5OTmjevDl27dqlM8+SJUtQq1YtqFQquLu7Y+TIkXqf4dKlS3j99dfh6uoKKysreHh44I033kBiYiKA3POHHj9+jDVr1kChUEChUBh03lbe9n3WNsxjbW0NZ2fnfOcLDg7GwYMH8fDhw0Kvn4goPwcPHkT9+vVhZWWFKlWq4Ouvv853vifPWT1+/DgUCkW+R8l27twJhUKBbdu2aafdunULgwcPhouLC1QqFWrVqoWVK1fqvO959QsANm7cCB8fH1hZWcHX1xdbtmzJt3bl5ORg/vz5qFWrFqysrODi4oJ33nkHjx49KsJW+k9UVBSaNm0KR0dH2NraokaNGvjf//6n9zkiIiLwv//9D66urrCxsUGXLl30jrAV9HoRz7Jr1y74+/vDysoKPj4+2Lx5s87rDx8+xAcffAA/Pz/Y2trC3t4eHTp0wJ9//qm3rPT0dEyZMgXVq1eHlZUV3Nzc8Nprr+HKlSvPXL8QAkOHDoWlpaXeup+2detWNGzYMN9Teo4ePYqOHTvCyckJNjY2qF27NhYsWKAzT0H2M5KTkzF27FhUqlQJKpUK5cuXR3BwME6ePAkAaNmyJX755Rdcv35dW8+ft/3LlSun06jm6d69OwDorH/Xrl2IjY3F5MmToVarkZqa+syhwkBuPb9+/TpiYmKeOQ/JE4+sUqlz+/ZtNGjQAAkJCRg6dCi8vb1x69YtREZGIjU1FZaWloiPj0eTJk2QmpqK0aNHo2zZslizZg26dOmCyMhI7R/OPJ9++iksLS3xwQcfICMjA5aWlgCA7OxstGvXDk2bNsWcOXNgbW2dbyYhBLp27YqDBw9i2LBhqFmzJrZs2YLQ0FC9eXfv3o0OHTrAy8sLU6ZMQVpaGhYuXIigoCCcPHlSrxD07NkT1apVw4wZMyCEeOZ2OXHiBDIzM1GvXj2916ZOnYopU6agSZMmmDZtGiwtLXH06FHs3bsXbdu2BZDbHE+dOhVt2rTB8OHDceHCBSxduhTR0dE4dOgQlEolMjMz0a5dO2RkZODdd9+Fq6srbt26hW3btiEhIQEODg747rvv8NZbb6FBgwba4chVqlR59i/0X5mZmUhKSkJaWhqOHz+OOXPmoGLFitpvYZ+UlJSEzMxM3L9/H99++y1iY2N1dn7yBAQEQAiBw4cPF/vFNYio5Dlz5gzatm0LZ2dnTJkyBdnZ2Zg8eTJcXFye+77AwEB4eXlhw4YNenUhIiICTk5OaNeuHQAgPj4ejRo1gkKhwKhRo+Ds7Ixff/0VQ4YMQVJSkt5wzPzq1y+//ILevXvDz88P4eHhePToEYYMGZLvqJt33nkHq1evxqBBgzB69GhcvXoVixYtwqlTp7R/+w31119/oXPnzqhduzamTZsGlUqFy5cv53vBoOnTp0OhUGDcuHG4e/cu5s+fjzZt2iAmJua5R90K6tKlS+jduzeGDRuG0NBQrFq1Cj179sSOHTsQHBwMIPdL7q1bt6Jnz56oXLky4uPj8fXXX6NFixY4e/Ys3N3dAeReH6Jz587Ys2cP3njjDYwZMwbJycmIiopCbGxsvjVPo9Fg8ODBiIiIwJYtW9CpU6dnZs3KykJ0dDSGDx+u91pUVBQ6d+4MNzc3jBkzBq6urjh37hy2bduGMWPGACj4fsawYcMQGRmJUaNGwcfHBw8ePMDBgwdx7tw51KtXDxMnTkRiYiLi4uLwxRdfAIBB18PIO9e0XLly2mm7d+8GkDsCITAwECdOnIClpSW6d++OJUuWoEyZMjrLCAgIAAAcOnQIdevWLXQGkpAgKmUGDBggzMzMRHR0tN5rOTk5Qgghxo4dKwCI33//XftacnKyqFy5sqhUqZLQaDRCCCH27dsnAAgvLy+Rmpqqs6zQ0FABQIwfP15vPaGhoaJixYra51u3bhUAxOzZs7XTsrOzRbNmzQQAsWrVKu10f39/Ub58efHgwQPttD///FOYmZmJAQMGaKdNnjxZABB9+vQp0HZZsWKFACDOnDmjM/3SpUvCzMxMdO/eXfu58+Rtr7t37wpLS0vRtm1bnXkWLVokAIiVK1cKIYQ4deqUACA2btz43Cw2NjYiNDS0QLnzrF+/XgDQPgIDA8Xp06fznbddu3ba+SwtLcU777wj0tLS9Oa7ffu2ACBmzZpVqCxERPnp1q2bsLKyEtevX9dOO3v2rDA3NxdP75JVrFhR5+/ghAkThFKpFA8fPtROy8jIEI6OjmLw4MHaaUOGDBFubm7i/v37Ost74403hIODg7ZWPa9++fn5CQ8PD5GcnKydtn//fgFAp3b9/vvvAoD4/vvvdd6/Y8eOfKc/6erVqwKA2Ldv3zPn+eKLLwQAce/evWfOk/c5XnnlFZGUlKSdvmHDBgFALFiwQDvt6dorhBAAxOTJk7XPV61aJQCIq1evaqdVrFhRABCbNm3STktMTBRubm6ibt262mnp6el6dfLq1atCpVKJadOmaaetXLlSABDz5s3T+zx5dTVv+3z++eciKytL9O7dW6jVarFz585nbos8ly9fFgDEwoULdaZnZ2eLypUri4oVK4pHjx7lu14hCr6f4eDgIEaOHPncLJ06ddLb5oWRkZEhfHx8ROXKlUVWVpZ2epcuXQQAUbZsWdGvXz8RGRkpJk2aJCwsLESTJk10Pk8eS0tLMXz4cIOzkDQ4DJhKlZycHGzduhUhISEIDAzUez3vPIbt27ejQYMGOsNmbW1tMXToUFy7dg1nz57VeV9oaOgzv7nN75vNp23fvh0WFhY685qbm+Pdd9/Vme+ff/5BTEwMBg4cqPOtYe3atREcHIzt27frLXvYsGEvXD8APHjwAADg5OSkM33r1q3IycnBJ598onceU9722r17NzIzMzF27Fided5++23Y29trh087ODgAyB22lpqaWqBcBdWqVStERUVh48aNGDZsGJRKJR4/fpzvvDNnzsSuXbvwzTffoFGjRsjMzER2drbefHnb4v79+0bNSkSlj0ajwc6dO9GtWzdUqFBBO71mzZrao6LP07t3b2RlZekM/9y1axcSEhLQu3dvALmjdDZt2oSQkBAIIXD//n3to127dkhMTNQO0czzdP26ffs2zpw5gwEDBugcBWvRogX8/Px03rtx40Y4ODggODhYZ10BAQGwtbXFvn37tPOmpKTozJM3TDgxMVFnet4pIUDutQMA4Mcff3zhhZEGDBgAOzs77fMePXrAzc0t37poCHd3d51RVfb29hgwYABOnTqlPfKnUqm0NVCj0eDBgwfaoctPbvdNmzahXLlyejUegN75lJmZmejZsye2bduG7du3a0czPc+z6vmpU6dw9epVjB07Vrttn15vYfYzHB0dcfToUdy+ffuFmQw1atQonD17FosWLYKFxX8DQvNur1e/fn2sXbsWr7/+OqZNm4ZPP/0Uhw8fxp49e/SW5eTkxHpugtisUqly7949JCUlwdfX97nzXb9+HTVq1NCbXrNmTe3rT6pcuXK+y7GwsICHh8cLc12/fh1ubm56w2OezpC33mdlu3//vl6D9qxszyKeGip85coVmJmZwcfH57n588tlaWkJLy8v7euVK1dGWFgYVqxYgXLlyqFdu3ZYvHixzs6JoVxcXNCmTRv06NEDS5cuRefOnREcHJzvper9/f0RHByMwYMHIyoqCseOHcv3vNi8bcGLMRBRUd27dw9paWmoVq2a3mv5/U1/Wp06deDt7Y2IiAjttIiICJQrVw6vvvqqdh0JCQlYtmwZnJ2ddR6DBg0CANy9e1dnuU/XiLy/1/mdQvH0tEuXLiExMRHly5fXW19KSorOuvKGJOc98k456datm870rl27at/Tu3dvBAUF4a233oKLiwveeOMNbNiwId/G9entqlAoULVqVZ37pRZF1apV9WpB9erVAUC7jpycHHzxxReoVq0aVCoVypUrB2dnZ5w+fVqnzl25cgU1atTQab6eJTw8HFu3bkVkZGShL0SVXz0H8Nx9oMLsZ8yePRuxsbHw9PREgwYNMGXKFPz999+Fyvg8n3/+OZYvX45PP/0UHTt21Hkt7wuWPn366Ezv27cvAODw4cN6yxNCsJ6bIJ6zSmQEzzqq+uS3rFIp6Lk6ZcuWBQA8evSoQA22oebOnYuBAwfixx9/xK5duzB69GiEh4fjyJEjRl1vjx49MHHiRPz444945513njmfpaUlunTpgpkzZyItLU1ne+V98//keTJERFLp3bs3pk+fjvv378POzg4//fQT+vTpo2168pq4N998M99rHgC5R8ieVJTzOXNyclC+fHl8//33+b7u7Oys/fmjjz7Cm2++qX0eHx+PN998E3PmzEGdOnW00588GqhWq3HgwAHs27cPv/zyC3bs2IGIiAi8+uqr2LVrF8zNzQ3O/jLMmDEDkyZNwuDBg/Hpp5+iTJkyMDMzw9ixYw2+ZU67du2wY8cOzJ49Gy1btoSVldUL3/NkPX+ZevXqhWbNmmHLli3YtWsXPv/8c8yaNQubN29Ghw4dirTs1atXY9y4cRg2bBg+/vhjvdfzzv99+nzv8uXLA8j/syckJLCemyA2q1SqODs7w97eHrGxsc+dr2LFirhw4YLe9PPnz2tfN6aKFStiz549SElJ0Tm6+nSGvPU+K1u5cuWee2ua5/H29gYAXL16VWeoV5UqVZCTk4OzZ8/C39//mfnzcnl5eWmnZ2Zm4urVq3pXGPbz84Ofnx8+/vhjHD58GEFBQfjqq6/w2WefATDOkcy0tDQAKNBR27S0NAghkJycrLPjdvXqVQD/HVEnIjKUs7Mz1Go1Ll26pPdafn/T89O7d29MnToVmzZtgouLC5KSkvDGG2/orMPOzg4ajeaZV3Z/kby/55cvX9Z77elpVapUwe7duxEUFPTCptfHx0dnhE7e0ciAgIDnHjE0MzND69at0bp1a8ybNw8zZszAxIkTsW/fPp3P+PR2FULg8uXLes25oS5fvqx3ZO7ixYsA/rv6fGRkJFq1aoVvvvlG571PN0lVqlTB0aNHkZWV9cILUDVq1AjDhg1D586d0bNnT2zZsuWFR2QrVKgAtVqtrWFPrhcAYmNjn/nvo7D7GW5ubhgxYgRGjBiBu3fvol69epg+fbq2WTWknv/4449466238Nprr2Hx4sX5zhMQEIDly5fj1q1bOtPzhiQ/+UUJkHuF7MzMTNZzE8RhwFSqmJmZoVu3bvj5559x/Phxvdfzhsx07NgRx44dwx9//KF97fHjx1i2bBkqVar03CGxhujYsSOys7OxdOlS7TSNRoOFCxfqzOfm5gZ/f3+sWbNG51YrsbGx2LVrl94wmcIICAiApaWl3nbp1q0bzMzMMG3aNL1vhvO2V5s2bWBpaYkvv/xSZ9jRN998g8TERO1VC5OSkvTODfXz84OZmRkyMjK002xsbF54y5k89+/fz/cqxytWrAAAnXOTnx7+BuTuRGzatAmenp7ab2TznDhxAgqFAo0bNy5QFiKiZzE3N0e7du2wdetW3LhxQzv93Llz2LlzZ4GWUbNmTfj5+SEiIgIRERFwc3ND8+bNddbx+uuvY9OmTfl+KXvv3r0XrsPd3R2+vr749ttvtecFAsBvv/2GM2fO6Mzbq1cvaDQafPrpp3rLyc7OLvDf8WfJ77ZheV+aPlkzAODbb79FcnKy9nlkZCT++eefIh/hy3P79m1s2bJF+zwpKQnffvst/P39tffyNjc316tHGzdu1GuoXn/9ddy/fx+LFi3SW09+9axNmzb44YcfsGPHDvTv3/+FR2mVSiUCAwP16nm9evVQuXJlzJ8/X+93k7fegu5naDQavS+Dy5cvD3d3d716XphTfQ4cOIA33ngDzZs3x/fff//M0Wldu3aFSqXCqlWrdLZHXu3Pu0JznhMnTgAAmjRpUuAsJA88skqlzowZM7Br1y60aNECQ4cORc2aNfHPP/9g48aNOHjwIBwdHTF+/HisX78eHTp0wOjRo1GmTBmsWbMGV69exaZNm4w+tDckJARBQUEYP348rl27pr1/W35/4D///HN06NABjRs3xpAhQ7SXlHdwcMCUKVMMzmBlZYW2bdti9+7dmDZtmnZ61apVMXHiRHz66ado1qwZXnvtNahUKkRHR8Pd3R3h4eFwdnbGhAkTMHXqVLRv3x5dunTBhQsXsGTJEtSvX1879Gvv3r0YNWoUevbsierVqyM7OxvfffeddgcrT0BAAHbv3o158+bB3d0dlStXRsOGDfPNvXbtWnz11Vfo1q0bvLy8kJycjJ07dyIqKgohISHac7kAoEOHDvDw8EDDhg1Rvnx53LhxA6tWrcLt27d1zgPLExUVhaCgIO2QKiKiopg6dSp27NiBZs2aYcSIEcjOzsbChQtRq1YtnD59ukDL6N27Nz755BNYWVlhyJAhevVo5syZ2LdvHxo2bIi3334bPj4+ePjwIU6ePIndu3cX6L7RM2bMQNeuXREUFIRBgwbh0aNHWLRoEXx9fXUa2BYtWuCdd95BeHg4YmJi0LZtWyiVSly6dAkbN27EggUL0KNHj8JtpCdMmzYNBw4cQKdOnVCxYkXcvXsXS5YsgYeHh959w8uUKYOmTZti0KBBiI+Px/z581G1alW8/fbbBq//SdWrV8eQIUMQHR0NFxcXrFy5EvHx8Vi1apV2ns6dO2PatGkYNGgQmjRpgjNnzuD777/XGXEE5F4M6ttvv0VYWBiOHTuGZs2a4fHjx9i9ezdGjBihc95unm7dumHVqlUYMGAA7O3tn3l/3jxdu3bFxIkTkZSUBHt7ewC5X9gvXboUISEh8Pf3x6BBg+Dm5obz58/jr7/+0n5pUpD9jOTkZHh4eKBHjx6oU6cObG1tsXv3bkRHR2Pu3LnaHAEBAYiIiEBYWBjq168PW1tbhISE5Jv5+vXr6NKlCxQKBXr06IGNGzfqvF67dm3tkXJXV1dMnDgRn3zyCdq3b49u3brhzz//xPLly9GnTx/Ur19f571RUVGoUKECb1tjior/AsRE0rt+/boYMGCAcHZ2FiqVSnh5eYmRI0eKjIwM7TxXrlwRPXr0EI6OjsLKyko0aNBAbNu2TWc5eZfMz+9WLKGhocLGxibf9ed3+fwHDx6I/v37C3t7e+Hg4CD69++vvdXLk7euEUKI3bt3i6CgIKFWq4W9vb0ICQkRZ8+e1Zkn79Y1z7vk/9M2b94sFAqFuHHjht5rK1euFHXr1hUqlUo4OTmJFi1aiKioKJ15Fi1aJLy9vYVSqRQuLi5i+PDhOpfH//vvv8XgwYNFlSpVhJWVlShTpoxo1aqV2L17t85yzp8/L5o3by7UarUA8Nzb2ERHR4uePXuKChUqCJVKJWxsbES9evXEvHnzdC5zn5evadOmoly5csLCwkI4OzuLkJAQceDAAb3lJiQkCEtLS7FixYoCbDkiooL57bffREBAgLC0tBReXl7iq6++0v69ftLTt67Jc+nSJe2ttw4ePJjvOuLj48XIkSOFp6enUCqVwtXVVbRu3VosW7ZMO8/z6pcQQvzwww/C29tbqFQq4evrK3766Sfx+uuvC29vb715ly1bJgICAoRarRZ2dnbCz89PfPTRR+L27dvP3A4FuXXNnj17RNeuXYW7u7uwtLQU7u7uok+fPuLixYt6n2P9+vViwoQJonz58kKtVotOnTrp3CJIiKLduqZTp05i586donbt2kKlUglvb2+9bZeeni7ef/994ebmJtRqtQgKChJ//PGHaNGihWjRooXOvKmpqWLixImicuXK2t9Rjx49xJUrV3S2z+eff67zviVLlggA4oMPPnjmdhMi99+AhYWF+O677/ReO3jwoAgODhZ2dnbCxsZG1K5dW+82Ny/az8jIyBAffvihqFOnjnY5derUEUuWLNFZTkpKiujbt69wdHTUu/XR0/J+l896PPl7EiL3djsLFy4U1atXF0qlUnh6eoqPP/5YZGZm6syn0WiEm5ub+Pjjj5+7zUieFELkM96AiEoljUYDHx8f9OrVK99hXaXJ/PnzMXv2bFy5csUoN5QnIjJ1/v7+cHZ2RlRUlNRRtPbv349WrVph48aNRTqKWxINGTIEFy9exO+//y51FElt3boVffv2xZUrV+Dm5iZ1HCoknrNKRFrm5uaYNm0aFi9erDPUq7TJysrCvHnz8PHHH7NRJaJSJysrS+/6Avv378eff/5Z6NunkHQmT56M6OhoHDp0SOookpo1axZGjRrFRtVE8cgqEREREWldu3YNbdq0wZtvvgl3d3ecP38eX331FRwcHBAbGyur8/h5ZJWoZOMFloiIiIhIy8nJCQEBAVixYgXu3bsHGxsbdOrUCTNnzpRVo0pEJR+PrBIREREREZHs8JxVIiIiIiIikh02q0RERERERCQ7Jn3Oak5ODm7fvg07OzsoFAqp4xAREQEAhBBITk6Gu7s7zMyk+16YdZKIiOSooHXSpJvV27dvw9PTU+oYRERE+bp58yY8PDwkWz/rJBERydmL6qRJN6t2dnYAcj+kvb29xGmIiIhyJSUlwdPTU1unpMI6SUREclTQOmnSzWrekCZ7e3sWYSIikh2ph96yThIRkZy9qE7yAktEREREREQkO2xWiYiIiIiISHbYrBIREREREZHssFklIiIiIiIi2WGzSkRERERERLIjebN669YtvPnmmyhbtizUajX8/Pxw/PhxqWMRERFJjjWSiIhKM0lvXfPo0SMEBQWhVatW+PXXX+Hs7IxLly7ByclJylhERESSY40kIqLSTtJmddasWfD09MSqVau00ypXrlzsOYQQSMvSFPt6C0OtNJf8fn1EZDqEEBBpaVLHMEkKtVoWf2/lUiMB06iTLyQE1MiQxe+2SJTWgKl/BpKEEALZGRlSx6ASwEKlKra/pZI2qz/99BPatWuHnj174rfffsMrr7yCESNG4O233853/oyMDGQ88T9ZUlKSUXKkZWng88lOoyzrZQms6ISNwxqbfpElopdOCIHrffsh7dQpqaOYpBonT0BhbS11jELXSKB018nnE4i0nIpAs4tSByk6z0bA4B1sWKlQhBD44ZOPcPviOamjUAkwek0klFZWxbIuSc9Z/fvvv7F06VJUq1YNO3fuxPDhwzF69GisWbMm3/nDw8Ph4OCgfXh6ehZzYukcv/7I9L/VJqJiIdLS2KiWAIWtkUDprpPPo0ZGyWhUAeDmESArVeoUZGKyMzLYqJJJUgghhFQrt7S0RGBgIA4fPqydNnr0aERHR+OPP/7Qmz+/b4w9PT2RmJgIe3t7g3PIeXhTaqYGgZ/tBgCcndYO1paSHgwnIhOQk5qKC/UCAADVDh2EmVotcSLTYoxhwElJSXBwcChSfSpsjQRKZ50siNSUJJT7MncIdeqY87C2MXxbSCYzFZhTNffn/90GLG2kzUMmJSs9HV+G9gAADF+2FkpV8RwVo5LJGMOAC1onJe183Nzc4OPjozOtZs2a2LRpU77zq1QqqFQqo+dQKBRsAomoRDJTq2EmgyGtVHiFrZEA6+QzWZr/97PSmo0elWpKlVWxDeEkKipJhwEHBQXhwoULOtMuXryIihUrSpSIiIhIHlgjiYiotJO0WX3vvfdw5MgRzJgxA5cvX8a6deuwbNkyjBw5UspYREREkmONJCKi0k7SZrV+/frYsmUL1q9fD19fX3z66aeYP38++vXrJ2UsIiIiybFGEhFRaSf5CSidO3dG586dpY5BREQkO6yRRERUmkl6ZJWIiIiIiIgoP2xWiYiIiIiISHbYrBIREREREZHssFklIiIiIiIi2WGzSkRERERERLLDZpWIiIiIiIhkh80qERERERERyQ6bVSIiIiIiIpIdNqtEREREREQkO2xWiYiIiIiISHbYrBIREREREZHssFklIiIiIiIi2WGzSkRERERERLLDZpWIiIiIiIhkh80qERERERERyQ6bVSIiIiIiIpIdNqtEREREREQkO2xWiYiIiIiISHbYrBIREREREZHssFklIiIiIiIi2WGzSkRERERERLLDZpWIiIiIiIhkh80qERERERERyQ6bVSIiIiIiIpIdNqtEREREREQkO2xWiYiIiIiISHbYrBIREREREZHssFklIiIiIiIi2WGzSkRERERERLIjabM6ZcoUKBQKnYe3t7eUkYiIiGSBNZKIiEo7C6kD1KpVC7t379Y+t7CQPBIREZEssEYSEVFpJnnVs7CwgKurq9QxZE5AjQwg8zFk8CvTp7QGFAqpUxAVOyEERFqa1DH05MgwExmGNZKIiEozyTufS5cuwd3dHVZWVmjcuDHCw8NRoUKFfOfNyMhARkaG9nlSUlJxxZSOEIi0nIpAs4vAHKnDPINnI2DwDjasVKoIIXC9bz+knToldRQqwQpTI4FSWieJiKjEkvSc1YYNG2L16tXYsWMHli5diqtXr6JZs2ZITk7Od/7w8HA4ODhoH56ensWcWAJZqbmNqpzdPAJkpUqdgqhYibQ02Teq6nr1oFCrpY5BBipsjQRKaZ0kIqISS9Ijqx06dND+XLt2bTRs2BAVK1bEhg0bMGTIEL35J0yYgLCwMO3zpKSkUlWIU8ech7WNvdQx/pOZCsypKnUKIslVO3QQZjJsChVqNRQc8WCyClsjAdZJIiIqWSQfBvwkR0dHVK9eHZcvX873dZVKBZVKVcypZERpDVjaSJ2CiJ5iplbDzNpa6hhUwr2oRgKsk0REVLLI6j6rKSkpuHLlCtzc3KSOQkREJCuskUREVNpI2qx+8MEH+O2333Dt2jUcPnwY3bt3h7m5Ofr06SNlLCIiIsmxRhIRUWkn6TDguLg49OnTBw8ePICzszOaNm2KI0eOwNnZWcpYREREkmONJCKi0k7SZvWHH36QcvVERESyxRpJRESlnazOWSUiIiIiIiIC2KwSERERERGRDLFZJSIiIiIiItlhs0pERERERESyw2aViIiIiIiIZIfNKhEREREREckOm1UiIiIiIiKSHTarREREREREJDtsVomIiIiIiEh22KwSERERERGR7LBZJSIiIiIiItlhs0pERERERESyw2aViIiIiIiIZIfNKhEREREREckOm1UiIiIiIiKSHTarREREREREJDtsVomIiIiIiEh22KwSERERERGR7LBZJSIiIiIiItlhs0pERERERESyw2aViIiIiIiIZIfNKhEREREREckOm1UiIiIiIiKSHTarREREREREJDtsVomIiIiIiEh22KwSERERERGR7FhIHYCIiIiIiKgk0Gg0yMrKkjqGrCiVSpibmxv0XjarRERERERERZSSkoK4uDgIIaSOIisKhQIeHh6wtbUt9HvZrBIRERERERWBRqNBXFwcrK2t4ezsDIVCIXUkWRBC4N69e4iLi0O1atUKfYRVNs3qzJkzMWHCBIwZMwbz58+XOg4REZGssE4SEclXVlYWhBBwdnaGWq2WOo6sODs749q1a8jKyip0syqLCyxFR0fj66+/Ru3ataWOQkREJDusk0REpoFHVPUVZZtIfmQ1JSUF/fr1w/Lly/HZZ59JHYcMlZkqdQJ9SmuAfzBKBCEERFqa1DF05MgsD5VcrJNERFQQFhYW8PX1RVZWFry8vPDdd9/B0dFR6lhFInmzOnLkSHTq1Alt2rR5YRHOyMhARkaG9nlSUtLLjkcFNaeq1An0eTYCBu9gw2rihBC43rcf0k6dkjoKkSRYJ4mIqCAcHR0RExMDAOjfvz8WL16MiRMnShuqiCQdBvzDDz/g5MmTCA8PL9D84eHhcHBw0D48PT1fckJ6LqV1bkMoVzePAFkyPOJLhSLS0mTdqKrr1YOC56bQS8I6SUREhggKCkJcXJzUMYpMsiOrN2/exJgxYxAVFQUrK6sCvWfChAkICwvTPk9KSmIhlpJCkXvkUm4NYWaqPI/0UpFVO3QQZjJrDBVqNc9PoZeCdZKIyDQJIZCWpTH6ctVK8wLtc2g0GkRFRWHw4MFGz1DcJGtWT5w4gbt376JevXraaRqNBgcOHMCiRYuQkZGhd7UolUoFlUpV3FHpeRQKwNJG6hRUSpip1TCztpY6BlGxYJ0kIjJNaVka+Hyy0+jLPTutHawtn92+JSQkwN/fX3ubmHbt2hk9Q3GTbBhw69atcebMGcTExGgfgYGB6NevH2JiYgp9WWMiIqKShHWSiIgKI++c1evXr0OhUGDJkiVSRyoyyY6s2tnZwdfXV2eajY0NypYtqzediIiotGGdJCIyTWqlOc5OM/5RTbWyYF9S2tjY4Msvv8Trr7+OESNGwMJC8mvqGsyg5I8fP4aNDYd+EhERERERPUmhUDx3uG5xCAwMhJ+fHzZs2IC+fftKmqUoDNqKLi4u6NWrFwYPHoymTZsaLcz+/fuNtiwiIqKShnWSiIie5f79+zrPt23bJlES4zHonNW1a9fi4cOHePXVV1G9enXMnDkTt2/fNnY2IiIiIiIiKqUMOrLarVs3dOvWDffu3cN3332H1atXY9KkSWjXrh0GDx6MLl26mPTYaCIiOdJoNMjKypI6Bj1FqVTyYkdEREQvQZE6SmdnZ4SFhSEsLAwLFy7Ehx9+iO3bt6NcuXIYNmwYxo8fD2veZoKIqMhSUlIQFxcHIYTUUegpCoUCHh4esLW1lToKERFRiVKkZjU+Ph5r1qzB6tWrcf36dfTo0QNDhgxBXFwcZs2ahSNHjmDXrl3GykpEVCppNBrExcXB2toazs7OBbohOBUPIQTu3bunvacdj7ASEREZj0HN6ubNm7Fq1Srs3LkTPj4+GDFiBN588004Ojpq52nSpAlq1qxprJxERKVWVlYWhBBwdnaGWq2WOg49xdnZGdeuXUNWVhabVSIiIiMyqFkdNGgQ3njjDRw6dAj169fPdx53d3dMnDixSOGIiOg/PKIqT/y9EBERvRwGNav//PPPC89FVavVmDx5skGhiIiIiIiIqHDWrFmDt99+G/Hx8XBycpI6TpEZdOua/fv3Y+fOnXrTd+7ciV9//bXIoYiISF4sLCzg7++PWrVqISQkBAkJCfnOl5KSgjZt2sjmQlCdOnXCo0ePpI5BRERULCIiIlC/fn1s2bJF6ihGYVCzOn78eGg0Gr3pQgiMHz++yKGIiEheHB0dERMTg7/++guOjo5YvHhxvvOtWLECvXr1emlDY3Nycgo1f79+/fDVV1+9lCxERERy8vDhQ1y8eBGzZ89GRESE1HGMwqBm9dKlS/Dx8dGb7u3tjcuXLxc5FBER5U8IgdTMbKM/CnMkNCgoCHFxcfm+tm7dOnTt2hUAcPv2bQQFBaFOnTqoXbs2Tp8+DQCYNWsWfH194efnh++//x5A7oidHj16aJfTo0cP7N+/HwBQtmxZjBo1Cn5+frh48SKmT58OPz8/1K5dG1988QUA4Pjx42jRogUCAgIQEhKChw8fAgA6d+5cYgo2ERGZCCGAzMfGf7ygVm/evBldu3ZFkyZNcOnSJdy/f7+YPvDLY9A5qw4ODvj7779RqVIlnemXL1+GjY2NMXIREVE+0rI08PlE/zSMojo7rR2sLV9cEjQaDaKiojB48GC91zIyMhAfHw8XFxcAwPr169GyZUtMnz4d2dnZyMzMRHR0NDZs2IDjx48jNTUV9evXR6tWrZ67zocPH6JDhw5YtGgRtm/fjr179+L48eNQqVR4+PAhsrKy8P7772PLli0oU6YMVq5cifDwcHz++eewt7dHeno6kpKSYG9vb9jGISIiKoysVGCGu/GX+7/bgOWze62IiAh89tlnUCgU6N69OzZt2oR33nnH+DmKkUHNateuXTF27Fhs2bIFVapUAZDbqL7//vvo0qWLUQMSEZH0EhIS4O/vr72faLt27fTmefDggc7FHOrXr48BAwbAwsICPXr0gJ+fHw4dOoTXX38dVlZWsLKyQuvWrREdHQ0HB4dnrlutVqNTp04AgN27d2PQoEFQqVQAgDJlyiA2NhZ//vknXn31VQBAdnY2atWqpX1/2bJlER8fz2aViIhKrLt37+LgwYPo3bs3ACAzMxPe3t6ls1mdPXs22rdvD29vb3h4eAAA4uLi0KxZM8yZM8eoAYmI6D9qpTnOTtNvFI2x3OfJO2f18ePHCA4OxpIlSzB69GideaysrJCenq593rx5cxw6dAjbtm1Dnz59MGPGjGcu38LCQud81IyMDO3PL7r6fE5ODurWrYt9+/bl+3p6ejrvT0tERMVHaZ17FPRlLPcZNm3ahGHDhmlPjwEALy8v3LlzB66ursbPUkwMOmfVwcEBhw8fxi+//IIRI0bg/fffx549e7B37144OjoaOSIREeVRKBSwtrQw+qOgF0SysbHBl19+iblz5yI7O1vntTJlyiAtLU07/fr163B1dcU777yD/v374/Tp02jatCk2b96MjIwMPHr0CHv37kWDBg1QoUIFnD17FtnZ2YiPj8fhw4fzXX+bNm2watUqbTP78OFDeHt74+bNmzhx4gSA3Eb3/Pnz2vc8ePAA7u4vYTgWERFRfhSK3OG6xn48p1ZHRESgW7duOtNCQkIQGRn5kj/sy2XQkVUgd4epbdu2aNu2rTHzEBGRzAUGBsLPzw8bNmxA3759dV5r0aIFjh49iqCgIOzfvx+ff/45lEolHB0dsX79eri6uqJnz54ICAiAQqHA1KlT4ebmBgDo2LEjfHx8UKNGDdStWzffdXfs2BEnTpxAvXr1oFQqMWjQIIwZMwYREREYM2YMkpOTodFoMGnSJHh7eyMmJgYNGjSAmZlB380SERGZhLyLEj5pwYIFxR/EyAxuVvfs2YM9e/bg7t27ercSWLlyZZGDERGRfDx9RcFt27blO9+IESOwevVqBAUFITQ0FKGhoXrzjBs3DuPGjdObPm/ePMybN++F6540aRImTZqkMy0gIAAHDx7Ue+/atWsxbNiwfLMSERGRvBnUrE6dOhXTpk1DYGAg3NzcXtr99IiIyLQ0atQIZ8+ehRBCFrXBx8dHe+ElIiIiMi0GNatfffUVVq9ejf79+xs7DxERmbj8bmsjFTllISIiosIx6CSezMxMNGnSxNhZiIiIiIiIiAAY2Ky+9dZbWLdunbGzEBEREREREQEwcBhweno6li1bht27d6N27dpQKpU6r+d3gQwiIiIiIiKigjLoyOrp06fh7+8PMzMzxMbG4tSpU9pHTEyMkSMSEZEcrFmzBpaWlnj06JHUUYiIiOgpFhYW8Pf3h7+/P+rXr//MvuzmzZvo0aMHACAmJga7du3SvvbVV18hIiLCoPWHhobi0qVLBr33WQw6srpv3z6jhiAiIvmLiIhA/fr1sWXLFl64iIiISGYcHR21DeqmTZswbdo0bN68WW++uXPnYujQoQBym9XY2Fi0bdsWAIp0u7d33nkHc+bMwddff23wMp7Gu6QTEdELPXz4EBcvXsTs2bMN/saViIiIikdSUhIcHR3zfW3btm149dVXodFo8Mknn+Dbb7+Fv78/tm/fjilTpmDRokUAgJYtW2LcuHEIDAyEr68v/vrrL+Tk5MDb2xuJiYkAgOTkZHh5eSE7OxuNGzfGvn37oNFojPY5DDqyCgDHjx/Hhg0bcOPGDWRmZuq8ll8HT0RERiAEkJVq/OUqrYHn3Bd18+bN6Nq1K5o0aYJLly7h/v37KFeunPFzEBERmTghBNKy04y+XLWF+rn3ME9ISIC/vz9SU1Px4MEDHD58WG+ev//+Gy4uLrCwyG0Dp02bhtjYWMyZMwcAcOzYMZ35lUoljh8/jpUrV2LevHn45ptv0KtXL2zYsAFvv/02IiMj0b17d+3yKlWqhHPnzsHX19con9mgZvWHH37AgAED0K5dO+zatQtt27bFxYsXER8fj+7duxslGBER5SMrFZjhbvzl/u82YGnzzJcjIiLw2WefQaFQoHv37ti0aRPeeecd4+cgIiIycWnZaWi4rqHRl3u071FYK62f+fqTw4AjIyMxcuRI7N69W2eeO3fuwNnZucDrzOvtAgIC8P333wMABg4ciNDQULz99tv47rvvMH/+fO38zs7O+Oeff4zWrBo0DHjGjBn44osv8PPPP8PS0hILFizA+fPn0atXL1SoUMEowYiISB7u3r2LgwcPonfv3qhUqRLWr1/PocBEREQy1rlz53yPrFpZWSE9Pb3Ay1GpVAAAc3Nz7fBeLy8vWFhYYO/evUhMTETt2rW186enp0OtVhcx/X8MOrJ65coVdOrUCQBgaWmJx48fQ6FQ4L333sOrr76KqVOnGi0gERE9QWmdexT0ZSz3GTZt2oRhw4bhiy++0E7z8vLCnTt34OrqavwsREREJkxtocbRvkdfynIL6vDhw/Dy8tKbXq1aNVy9elX73M7ODsnJyYXOMnDgQLz55pv46KOPdKZfuXIFNWvWLPTynsWgZtXJyUn7oV555RXExsbCz88PCQkJSE19CedSERFRLoXiucN1X4aIiAi9LyFDQkIQGRmJUaNGFWsWIiIiuVMoFM8drvuy5J2zKoSAhYUFli1bpjePnZ0dXF1dERcXBw8PD7Rq1QozZ85E3bp1MX369AKvq0ePHhg2bBj69u2rnfbgwQNYW1ujbNmyRvk8gIHNavPmzREVFQU/Pz/07NkTY8aMwd69exEVFYXWrVsXeDlLly7F0qVLce3aNQBArVq18Mknn6BDhw6GxCIiopdg//79etMWLFhQ/EFKGdZIIiIqjOzs7ALNN3z4cKxduxbjx49HmTJlEB0drX2tY8eO2p+frP++vr46z6Ojo9G+fXuUL19eO239+vV46623DP8A+TCoWV20aJF2rPPEiROhVCpx+PBhvP766/j4448LvBwPDw/MnDkT1apVgxACa9asQdeuXXHq1CnUqlXLkGhEREQlAmskERG9DL1798bq1asNfv/UqVOxevVqbN26VWe6g4MD+vTpU7RwTzGoWS1Tpoz2ZzMzM4wfP96glYeEhOg8nz59OpYuXYojR46wEBMRUanGGklERC+DQqHAoEGDDH7/5MmTMXnyZL3p/fv3L0qsfBl8n1WNRoMtW7bg3LlzAAAfHx907dpVe48dQ5a3ceNGPH78GI0bNzY0FpGuTJmeQ/2Ce1pKQQgBkWb8e4IVVY4MMxEVN9ZI0iHX2lYQMqx/BSWEQHZGhtQxDJKVUfCrvxLJiUGd5V9//YUuXbrgzp07qFGjBgBg1qxZcHZ2xs8//1yo++qcOXMGjRs3Rnp6OmxtbbFlyxb4+PjkO29GRgYynvgjkZSUZEh8Kk3mVJU6Qf48GwGDd8imYAshcL1vP6SdOiV1FCJ6QmFqJMA6WWrItbYVhMzqX0EJIfDDJx/h9sVzUkchKlUMus/qW2+9hVq1aiEuLg4nT57EyZMncfPmTdSuXRtDhw4t1LJq1KiBmJgYHD16FMOHD0doaCjOnj2b77zh4eFwcHDQPjw9PQ2JTyWd0jq3GMrZzSNAlny+GRdpabJvVNX16kFhxPt2EZmCwtRIgHWyRDOF2lYQMqt/BZWdkVEiGlX3Gj6w+Pe+mUSmwKAjqzExMTh+/DicnJy005ycnDB9+nTUr1+/UMuytLRE1aq53xAGBAQgOjoaCxYswNdff60374QJExAWFqZ9npSUxEJM+hSK3G9t5VgMM1Nl/414tUMHYSbDplChVkNhYt/EExVVYWokwDpZosm5thWECdS/ghq+bC2UKiupYxjEQqViLSWTYlCzWr16dcTHx+td4OHu3bvaomqonJwcnSFMT1KpVFDx2yAqCAnuRVlSmKnVMLMu/nuDkbxZWFhoT/FQKpVYvnw5/P399ea7efMm3nvvPURGRmL16tWIjY3FnDlzCrSOa9eu4dixY+jVqxeA3C9G7969i7Zt2wIApkyZgnLlyr30e7teu3YNPXr0wPHjx3WmL1++HObm5hg8ePBLXf/zPK9GAqyTJR5rmywoVVZQWplms0olX1xcHEaPHo0///wTTk5OqFy5MhYtWgQXFxepoxnEoGHA4eHhGD16NCIjIxEXF4e4uDhERkZi7NixmDVrFpKSkrSP55kwYQIOHDiAa9eu4cyZM5gwYQL279+Pfv36GfRhiIjo5XB0dERMTAxiYmIwfvx4TJs2Ld/55s6dW+jTQfJcu3YNGzZs0D6PiYnBrl27DFrWy9C/f38sXbq02NbHGklERIUhhEDXrl3RqVMnXLlyBcePH8fo0aNx7949qaMZzKBmtXPnzjh79ix69eqFihUromLFiujVqxdiY2MREhICJycnODo66gwTzs/du3cxYMAA1KhRA61bt0Z0dDR27tyJ4OBggz4MERG9fElJSXB0dMz3tW3btuHVV1/VPr969SqaN2+O6tWrY/78+drps2bNgq+vL/z8/PD9998DyL1v9+7du+Hv74+vv/4an3zyCb799lv4+/tj+/btOuu5cuUK2rVrh8DAQLz66qu4du0aAKBly5YYN24cAgMD4evri7/++gsA8PjxYwwcOBD169dHQEAAoqKiAADJycno378/ateujTp16uD333/XWU9sbCwCAwNx5coVWFlZoWLFijh58mRRNl+BsUYSEVFh7NmzB7a2thgyZIh2WrNmzQp18Vu5MWgY8L59+4yy8m+++cYoyyEiKi2EEEjLNv7tfNQWzz8nOCEhAf7+/khNTcWDBw9w+PBhvXn+/vtvuLi46NzCLDo6GqdPn4aFhQUCAwMREhKChw8fYsOGDTh+/DhSU1NRv359tGrVCtOnT8eiRYsQGRkJIHdI65PDiI8dO6Zd7ogRI/D111+jUqVK2Lt3Lz788ENs3LgRQO4w5ePHj2PlypWYN28evvnmG0yfPh2dO3fG6tWrcf/+fTRt2hTnzp3Dp59+igoVKuC7775DTk4OkpOT8ejRIwDA6dOnMWjQIGzYsAFVqlQBANSrVw+HDx9GvXr1irjFX4w1kojINL2s2wG+6PodZ8+eLZb6VJwMalZbtGhh7BxERFQAadlpaLiuodGXe7TvUVgrn32uct4wYACIjIzEyJEjsXv3bp157ty5A2dnZ51p7du31x6F7dixI/744w/cv38fr7/+OqysrGBlZaU9aujg4FCgrCkpKfj999/RrVs3ALk7BTY2/53H1717dwC5FyTKO2q7a9cubNu2DZ999hmA3COt8fHx2L17N3766ScAgJmZGRwcHPDo0SPcvn0bb7zxBrZt2wYvLy/tsp2dnbVHcYmIiPIj0tJwoV6A0Zdb4+QJKErZdUUMalYPHDjw3NebN29uUBgiIpK/zp07Y8CAAXrTrayskJ6ue+P5J78BVigURrkKZU5ODlxcXLTN89PyLjBkbm4OjUajfc/PP/+MihUrFmgdTk5OKFeuHI4dO6bTrKanp0Mtw6tlExER1axZE5s3b5Y6hlEZ1Ky2bNlSb9qTOyB5OwdERGRcags1jvY9+lKWW1CHDx/WaeDyVKtWDVevXtWZtmPHDiQmJsLCwgK//vorRowYgUePHmHYsGF4//33kZqair1792Lq1Km4ffs2kpOTte+1s7PTeZ7H3t4eLi4u+PnnnxESEgKNRoNz584995yctm3b4ssvv8TcuXMB5F68yd/fH23atMHSpUsxffp07TBgAFCr1fjpp58QHBwMR0dHtG/fHgBw+fJlfiFLRETPpVCrUePkiZey3Odp06YNxo0bh9WrV2PgwIEAgIMHD8LR0dFkz1s16AJLjx490nncvXsXO3bsQP369WV15UYiopJGoVDAWmlt9MeLjnjmnbNap04dfPjhh1i2bJnePHZ2dnB1dUVcXJx2Wv369RESEoK6deti6NChqFKlCgIDA9GzZ08EBASgefPmmDp1Ktzc3FC7dm1kZWXB398fK1asQKtWrXDy5EnUrVtX7wJL69atw8KFC1GnTh34+flhz549z80/adIkJCYmonbt2vDx8dGeBztp0iRcu3YNfn5+qFevHs6cOaN9j4ODA37++WeMGzcOhw4dAgAcPXpU5wJSRERET1MoFDCztjb640W1WqFQYOvWrdi6dSuqVKmCWrVqYeHChXqn6JgSg46s5ndeUXBwMCwtLREWFoYTJ4z/TQIREUknOzu7QPMNHz4ca9euxfjx4zFw4EDtN7tPGzduHMaNG6czTalUYu/evTrToqOjtT937NhR+7OXl1e+X47u379f+7Ovr6/2uY2NDVasWKE3v52dnfa81ifl3WPVxcUFf/75J4DcKwNXr179hVe6JyIikkqFChWwdetWqWMYjUFHVp/FxcUFFy5cMOYiiYjIhPTu3dtkbzz+Ivfu3Xvm/WWJiIjI+Aw6snr69Gmd50II/PPPP5g5cyb8/f2NkYuIiEyQQqHAoEGDpI7xUrRq1UrqCERERKWKQc2qv78/FAoFhBA60xs1aoSVK1caJRgRERERERGVXgY1q09f7dHMzAzOzs6wsrIySigiIiIiIiIq3QxqVgt6nzoiIiIiIiIiQxh0gaXRo0fjyy+/1Ju+aNEijB07tqiZiIiIiIiIqJQzqFndtGkTgoKC9KY3adIEkZGRRQ5FRETyMm3aNNSqVQt+fn4IDAzUOx0kT/fu3XH79u1nLqdjx45IS0t77roqVaqElJQUnWkXLlx45m1wiIiICLCwsEDdunXh4+ODgIAALF++XOpIRWbQMOAHDx7ke69Ve3t73L9/v8ihiIhIPg4fPox9+/YhJiYGSqUScXFxsLGx0ZsvJiYGarUa7u7uz1zW9u3bDcpQo0YNxMfHIy4uDh4eHgYtg4iIqCRzdHTEqVOnAAA3btxAt27dIITA0KFDJU5mOIOOrFatWhU7duzQm/7rr7/Cy8uryKGIiCh/QgjkpKYa/fH01d2fdOfOHZQrVw5KpRIA4OHhAScnJ7351q1bh65duwIANBoN3nzzTfj4+MDPzw+rVq0C8N9R02vXrqFOnToIDQ1FzZo10bt3b70MSUlJaNmyJX7++WcAQKdOnbBhwwajbEciIqKXRQiBrAyN0R/Pq9VPq1ChAubOnYslS5a8xE/68hl0ZDUsLAyjRo3CvXv38OqrrwIA9uzZg7lz52L+/PnGzEdERE8QaWm4UC/A6MutcfIEFNbW+b4WHByMyZMnw8fHB8HBwejfvz8CAwP15jty5AjefvttALlHWa9evYqzZ88CABITE/XmP3fuHNavX4+aNWuiVatWOHjwIJo1a6adv0+fPvjoo4/QuXNnAEC9evUwb948hIWFGeUzExERvQzZmTlYNuY3oy936IIWUKrMCzx/vXr1cOHCBaPnKE4GHVkdPHgw5s6di2+++QatWrVCq1atsHbtWixdulS7o0JERCWDnZ0dTp06hQULFkCtViM4OBhRUVF68925cwfOzs4AAC8vL9y+fRsjR47Erl278j11pEaNGvDx8YFCoUDdunVx7do17WsdO3bEhx9+qG1UAcDZ2Rn//POP8T8gERFRCVSYI7FyZdCRVQAYPnw4hg8fjnv37kGtVsPW1taYuYiIKB8KtRo1Tp54Kct9HgsLCwQHByM4OBjlypXDjz/+iODgYJ15rKyskJ6eDgBwcnLCmTNnsH37dnzxxRfYtWsX5syZozO/SqXS/mxubg6NRqN93qRJE/z6668ICQnRTktPT4f6BTmJiIikZmFphqELWryU5RZGTEwMvL29jZ6jOBnUrF69ehXZ2dmoVq2a9lt0ALh06RKUSiUqVapkrHxERPQEhULxzOG6L8uFCxdgYWGBKlWqQAiB2NhY+Pj46M3n7e2Ny5cvw9XVFffv34elpSV69eqFihUrYuLEiYVa5+zZszFmzBhMmjQJn376KQDg8uXLqFmzplE+ExER0cuiUCgKNVz3Zbh58yY++OADjBo1StIcRWXQMOCBAwfi8OHDetOPHj3KWwsQEZUwKSkpePPNN1GrVi34+voiJycH7777rt587du3x2+/5Z6jc+vWLbRo0QJ16tTBiBEjMHny5EKtU6FQYPny5YiNjdVeC+G3335Dhw4divx5iIiISqKEhAT4+/vDx8cH3bp1w7BhwzBkyBCpYxWJQUdWT506le99Vhs1amTy3TsREekKCAjAH3/88cL5evXqhQ4dOuB///sf6tSpo718/pPyzku1tbXF8ePHtdOfHCL85LmrW7ZsAQBkZmbi+PHjmDt3roGfgoiIqGTLzs6WOoLRGXRkVaFQIDk5WW96YmKizjlHRERUetja2mLChAmIj483+rJv3bqFzz77DObm0g6rIiIiouJjULPavHlzhIeH6zSmGo0G4eHhaNq0qdHCERGRaenYsSNcXV2NvtzKlSujZcuWRl8uERERyZdBw4BnzZqF5s2bo0aNGtp74v3+++9ISkrC3r17jRqQiIiIiIiISh+Djqz6+Pjg9OnT6N27N+7evYvk5GQMGDAA58+fh6+vr7EzEhERSsb90koi/l6IiIheDoPvs2ptbY0yZcrAzc0NQO65SjyXiIjI+JRKJRQKBe7duwdnZ2coFAqpI9G/hBC4d+9e7m0KlEqp4xAREZUoBjWrx48fR7t27aBWq9GgQQMAwBdffIEZM2Zg165dqFevnlFDEhGVZubm5vDw8EBcXJzOlXJJHhQKBTw8PPiFLRERkZEZ1Ky+99576NKlC5YvXw4Li9xFZGdn46233sLYsWNx4MABo4YkIirtbG1tUa1aNWRlZUkdhZ6iVCrZqBIRkSxMmzYNERERMDMzg0qlwsaNG1G5cmWpYxnM4COrTzaqAGBhYYGPPvoIgYGBRgtHRET/MTc3Z1NERERE+Tp8+DD27duHmJgYKJVKxMXFwcbGRupYRWJQs2pvb48bN27A29tbZ/rNmzdhZ2dnlGBERERERESmRgiB7IwMoy/XQqV67nUr7ty5g3LlymmvoeDh4WH0DMXNoGa1d+/eGDJkCObMmYMmTZoAAA4dOoQPP/wQffr0KfBywsPDsXnzZpw/fx5qtRpNmjTBrFmzUKNGDUNiERERlRiskUREpik7IwNfhvYw+nJHr4mE0srqma8HBwdj8uTJ8PHxQXBwMPr372/yo14NunXNnDlz8Nprr2HAgAGoVKkSKlWqhIEDB6JHjx6YNWtWgZfz22+/YeTIkThy5AiioqKQlZWFtm3b4vHjx4bEIiIiKjFYI4mIqDDs7Oxw6tQpLFiwAGq1GsHBwYiKipI6VpEYdGTV0tISCxYsQHh4OK5cuQIAqFKlCqytrQu1nB07dug8X716NcqXL48TJ06gefPmhkQjIiIqEVgjiYhMk4VKhdFrIl/Kcl84j4UFgoODERwcjHLlyuHHH39EcHCw0bMUF4Pvswrk3mvVz8/PWFmQmJgIAChTpozRlkkkN0IAQqMAEh8AyjSp4wAActLkkYOIno01koQQEKb69zozDSJLgWyogITHgDJH6kSFkpWR/t/PmRpAoZEwjeEsLM14r+5ioFAonjtc92W5cOECLCwsUKVKFQghEBsbCx8fn2LPYUxFalaNKScnB2PHjkVQUBB8fX3znScjIwMZT5ysnJSUVFzxiIxCCIHre8oh7b4lENlO6jj5E0LqBET0lILUSIB1siQTQuB6335IO3VK6igGEQBO1v0ciQ5VgE9ipY5TaEL8d9uwlR8ehEKhlDCN4dyqOKD7B/XYsJZQKSkpGDVqlPZvf0BAAN59912JUxWNbJrVkSNHIjY2FgcPHnzmPOHh4Zg6dWoxpiIyLpGtyG1UZUpdLgMKCzarRHJTkBoJsE6WZCItzWQbVQDIMbPMbVRJUv9cSUR2Zg6UKt4GrSQKCAjAH3/8IXUMo5JFszpq1Chs27YNBw4ceO4llidMmICwsDDt86SkJHh6ehZHRCLjeOKbzGr7dsFMrZYwzBOyUoEFtaEwF/y2lUhmClojAdbJ0qLaoYPyqR8FlJWSgt/+PaIa+rEPlCZ278eMxAQsH5f788Cp9aByMq3h+FkZGqz66PlfdhHJkaTNqhAC7777LrZs2YL9+/ejcuXKz51fpVJBVYATi4lMgZlDWZgV8qJkL03mY4BHVIlkpbA1EmCdLC3M1Gr51I8CMtNka39WOVhDaWcnYRpD/DcMWKky55FJomIiabM6cuRIrFu3Dj/++CPs7Oxw584dAICDgwPUJvaNIRERkTGxRhIRUWln0H1WjWXp0qVITExEy5Yt4ebmpn1ERERIGYuIiEhyrJFERKZH8EKVeoqyTSQfBkxERET6WCOJiEyHUqmEQqHAvXv34OzszGuA/EsIgXv37uXezkdZ+Ktoy+ICS0RERERERKbK3NwcHh4eiIuLw7Vr16SOIysKhQIeHh4wNy/8ud5sVomIiIiIiIrI1tYW1apVQ1ZW1otnLkWUSqVBjSrAZpWIiIiIiMgozM3NDW7MSJ+kF1giIiIiIiIiyg+bVSIiIiIiIpIdNqtEREREREQkO2xWiYiIiIiISHbYrBIREREREZHssFklIiIiIiIi2WGzSkRERERERLLDZpWIiIiIiIhkh80qERERERERyQ6bVSIiIiIiIpIdNqtEREREREQkO2xWiYiIiIiISHbYrBIREREREZHssFklIiIiIiIi2WGzSkRERERERLLDZpWIiIiIiIhkh80qERERERERyQ6bVSIiIiIiIpIdNqtEREREREQkO2xWiYiIiIiISHbYrBIREREREZHssFklIiIiIiIi2WGzSkRERERERLLDZpWIiIiIiIhkh80qERERERERyQ6bVSIiIiIiIpIdSZvVAwcOICQkBO7u7lAoFNi6dauUcYiIiGSFdZKIiEozSZvVx48fo06dOli8eLGUMYiIiGSJdZKIiEozCylX3qFDB3To0EHKCERERLLFOklERKWZpM0qFU5qpgbIzJY6hg4hBNI16VLH0COEANJlmCstTfvzw7QUKKCRMM0TMh9DrVAAANIS4wGltcSB9FkpzaH4N6OsWKgBOeYCoLZQy3ObEVGRpGalwSxL6hSFk5WdnlubkY2s5ARAmNYHyEpO1P6cmpUOTVaqhGkKLytLJvsbRIVkUs1qRkYGMjIytM+TkpIkTFP8ms3ehzRYSR3jCQLWFb+CufV1qYPoEgLTvtPA+5bUQZ6v/eZXkWEpo0aikmfuf7d3lTYHGU3d8nWxpv0aNqylSGmvkyVZbqOXq+WGFvKqHwVgkW2J/skBEJrbWPqe1GmKpv3mdkhVS52icCw0lngLnwPQ/bdEJHcmdTXg8PBwODg4aB+enp5SR3rp1EpzqSM8myJLfo0qAFUWZN+onvcAMpRSp6CS7tTdU0jLTnvxjFRilMY6WVqkZctvtFBhmGsUEJrbUscoMqfHadCYmXazx7pApsSkjqxOmDABYWFh2udJSUklvhA/eUTkxMdtAEsbCdPoSstOQ8uNnwAAfu2+B2oLeXzNKNLScG9uCwBAud07oFDLI9eTnK2ssF9uR7uEAGQ4rCk1U4O2XxwAABwc1wpqS5n82cpMBRbUzv35w8uyGjqdlp2GlhtaSh2DJFAa62RptOP1HbC2KyN1jEJ59PABftgzHADQb/Z8lHVykjhR4aQ8vI24jr1hniPQpeuPsC7jIXWkQkl+nILIY7FSxyAqNJns9RWMSqWCSqWSOoZkrC0tALnsqAOA4r+jvmXUtrCWyc56Dsxx79+fy5Zxhpm1PHKZBnupA+hRZ2YjTeTmUluXy/3/QA4sHuc2+EDueasy+fdPpVtpr5OlhZW5WjY1t6BSLR5rf7awdYDSvqyEaQpPmZUKi5zcv/lWFqa3/bMseM4qmSZJ9/pSUlJw+fJl7fOrV68iJiYGZcqUQYUKFSRMRkREJD3WSSIiKs0kbVaPHz+OVq1aaZ/nDV0KDQ3F6tWrJUpFREQkD6yTRERUmknarLZs2ZJXJCMiInoG1kkiIirNTOpqwERERERERFQ6sFklIiIiIiIi2WGzSkRERERERLLDZpWIiIiIiIhkh80qERERERERyQ6bVSIiIiIiIpIdNqtEREREREQkO2xWiYiIiIiISHbYrBIREREREZHssFklIiIiIiIi2WGzSkRERERERLLDZpWIiIiIiIhkh80qERERERERyQ6bVSIiIiIiIpIdNqtEREREREQkO2xWiYiIiIiISHbYrBIREREREZHssFklIiIiIiIi2WGzSkRERERERLLDZpWIiIiIiIhkh80qERERERERyQ6bVSIiIiIiIpIdNqtEREREREQkO2xWiYiIiIiISHbYrBIREREREZHssFklIiIiIiIi2WGzSkRERERERLLDZpWIiIiIiIhkh80qERERERERyY4smtXFixejUqVKsLKyQsOGDXHs2DGpIxEREckCayQREZVWkjerERERCAsLw+TJk3Hy5EnUqVMH7dq1w927d6WORkREJCnWSCIiKs0kb1bnzZuHt99+G4MGDYKPjw+++uorWFtbY+XKlVJHIyIikhRrJBERlWYWUq48MzMTJ06cwIQJE7TTzMzM0KZNG/zxxx/FliNHo0Fqwu1iW1+hZKYBmn+/U0h+CFikS5vnCemaNKgyBQAgJzUNOUqJA/0rJy1N6gj0kqRmaqSO8J/MbFj/+2Pq4yQgM1vSOE9Ky/7v/4GHiXeRZmElYRrT42TnDDNzc6ljyKZGAkB2VhYe3rparOs0pszUFKRk2gAAEq7fgLWNvcSJCi81JQEZ5ioAQMKjFKRnSboLV2iJCY+1Pyekp0KRmixhmsJLTf0vf9qDe5DNTk8BpT2R/17cDaSpWRfIcOXcK8BCWTz/D0j6l+7+/fvQaDRwcXHRme7i4oLz58/rzZ+RkYGMjAzt86SkJKPkSE24jZtBbY2yrJfDNfc/G+WX8bt//3tzblNJc1DpEPjZbqkjaKmRjnP/1nrrBd7ShnmaQgFU8gQAdPglROIwpmd/1yiUdXSVOkahayTw8urkw1tX8d24MKMsSzr1cv8zbcLzZ5MzX4/c/44bKW2OInrjl9eRocp48YwyosoU2n2e+136SZrFEBozS6D5FwCAnfPjJU5Dpu71CYBrxSrFsi7JhwEXRnh4OBwcHLQPT09PqSORjKnr1YNCrZY6BhWRWmmOwIpOUsfQkwYVonOqSx0jX2ohUDddPqMwqPiwTpLcZavKIENpWo0qAGQogfMeUqcwnFlOJhwSr0gdg6jQFEIIIdXKMzMzYW1tjcjISHTr1k07PTQ0FAkJCfjxxx915s/vG2NPT08kJibC3t7wIT2yHgacx1yde7REhtQWVlDIMJtCrZZlLio8IQTSsmQ0BDiPEEBWqtQp8iWEQLqGDashjDEMOCkpCQ4ODkWqT4WtkcDLq5OmPgwYyP1/Igsq2dbSgrK0toeZmUkda9ASQgBqhenmz8mBKjXVZPctcnJykJrOU6Wo6IwxDLigdVLSYcCWlpYICAjAnj17tIU4JycHe/bswahRo/TmV6lUUKlURs9hZm4O27L89plIrhQKBawtZXp+lspB6gTPZCN1ACqSwtZI4OXVSQulEuUryXMkAVGxspXv3/yCML2ztam0k3zvLywsDKGhoQgMDESDBg0wf/58PH78GIMGDZI6GhERkaRYI4mIqDSTvFnt3bs37t27h08++QR37tyBv78/duzYoXdBCSIiotKGNZKIiEozSc9ZLSpjnBNERERkbHKpT3LJQURE9KSC1ifTPMOdiIiIiIiISjQ2q0RERERERCQ7bFaJiIiIiIhIdtisEhERERERkeywWSUiIiIiIiLZkfzWNUWRdyHjpKQkiZMQERH9J68uSX3BfdZJIiKSo4LWSZNuVpOTkwEAnp6eEichIiLSl5ycDAcHB0nXD7BOEhGRPL2oTpr0fVZzcnJw+/Zt2NnZQaFQFGlZSUlJ8PT0xM2bN3kvukLgdis8brPC4zYrPG6zwjPmNhNCIDk5Ge7u7jAzk+6MG9bJ/5h6fsD0PwPzS4v5pWfqn0GKOmnSR1bNzMzg4eFh1GXa29ub5D8eqXG7FR63WeFxmxUet1nhGWubSXlENQ/rpD5Tzw+Y/mdgfmkxv/RM/TMUZ53kBZaIiIiIiIhIdtisEhERERERkeywWf2XSqXC5MmToVKppI5iUrjdCo/brPC4zQqP26zwuM2ez9S3j6nnB0z/MzC/tJhfeqb+GaTIb9IXWCIiIiIiIqKSiUdWiYiIiIiISHbYrBIREREREZHssFklIiIiIiIi2WGz+q/FixejUqVKsLKyQsOGDXHs2DGpI8lWeHg46tevDzs7O5QvXx7dunXDhQsXpI5lUmbOnAmFQoGxY8dKHUXWbt26hTfffBNly5aFWq2Gn58fjh8/LnUs2dJoNJg0aRIqV64MtVqNKlWq4NNPPwUvTaDrwIEDCAkJgbu7OxQKBbZu3arzuhACn3zyCdzc3KBWq9GmTRtcunRJmrAyYqp18kW/b7kz9Zq7dOlS1K5dW3tfxsaNG+PXX3+VOpbBTLF+T5kyBQqFQufh7e0tdaxCMeX9gUqVKultf4VCgZEjR0odrUCk3rdgswogIiICYWFhmDx5Mk6ePIk6deqgXbt2uHv3rtTRZOm3337DyJEjceTIEURFRSErKwtt27bF48ePpY5mEqKjo/H111+jdu3aUkeRtUePHiEoKAhKpRK//vorzp49i7lz58LJyUnqaLI1a9YsLF26FIsWLcK5c+cwa9YszJ49GwsXLpQ6mqw8fvwYderUweLFi/N9ffbs2fjyyy/x1Vdf4ejRo7CxsUG7du2Qnp5ezEnlw5Tr5It+33Jn6jXXw8MDM2fOxIkTJ3D8+HG8+uqr6Nq1K/766y+poxWaKdfvWrVq4Z9//tE+Dh48KHWkAjP1/YHo6GidbR8VFQUA6Nmzp8TJCkbyfQtBokGDBmLkyJHa5xqNRri7u4vw8HAJU5mOu3fvCgDit99+kzqK7CUnJ4tq1aqJqKgo0aJFCzFmzBipI8nWuHHjRNOmTaWOYVI6deokBg8erDPttddeE/369ZMokfwBEFu2bNE+z8nJEa6uruLzzz/XTktISBAqlUqsX79egoTyUFLq5NO/b1NUEmquk5OTWLFihdQxCsWU6/fkyZNFnTp1pI5hsJK2PzBmzBhRpUoVkZOTI3WUApF636LUH1nNzMzEiRMn0KZNG+00MzMztGnTBn/88YeEyUxHYmIiAKBMmTISJ5G/kSNHolOnTjr/3ih/P/30EwIDA9GzZ0+UL18edevWxfLly6WOJWtNmjTBnj17cPHiRQDAn3/+iYMHD6JDhw4SJzMdV69exZ07d3T+H3VwcEDDhg1LbU1gnZQXU665Go0GP/zwAx4/fozGjRtLHadQTL1+X7p0Ce7u7vDy8kK/fv1w48YNqSMVWEnaH8jMzMTatWsxePBgKBQKqeMUiNT7FhbFshYZu3//PjQaDVxcXHSmu7i44Pz58xKlMh05OTkYO3YsgoKC4OvrK3UcWfvhhx9w8uRJREdHSx3FJPz9999YunQpwsLC8L///Q/R0dEYPXo0LC0tERoaKnU8WRo/fjySkpLg7e0Nc3NzaDQaTJ8+Hf369ZM6msm4c+cOAORbE/JeK21YJ+XDVGvumTNn0LhxY6Snp8PW1hZbtmyBj4+P1LEKzNTrd8OGDbF69WrUqFED//zzD6ZOnYpmzZohNjYWdnZ2Usd7oZK0P7B161YkJCRg4MCBUkcpMKn3LUp9s0pFM3LkSMTGxprUuQ9SuHnzJsaMGYOoqChYWVlJHcck5OTkIDAwEDNmzAAA1K1bF7Gxsfjqq69MrjgVlw0bNuD777/HunXrUKtWLcTExGDs2LFwd3fnNiMqAUy15taoUQMxMTFITExEZGQkQkND8dtvv5lEw1oS6veTR8Bq166Nhg0bomLFitiwYQOGDBkiYbKCKUn7A9988w06dOgAd3d3qaMUmNT7FqW+WS1XrhzMzc0RHx+vMz0+Ph6urq4SpTINo0aNwrZt23DgwAF4eHhIHUfWTpw4gbt376JevXraaRqNBgcOHMCiRYuQkZEBc3NzCRPKj5ubm96OTM2aNbFp0yaJEsnfhx9+iPHjx+ONN94AAPj5+eH69esIDw83uYIulby/+/Hx8XBzc9NOj4+Ph7+/v0SppMU6KQ+mXHMtLS1RtWpVAEBAQACio6OxYMECfP311xIne7GSWL8dHR1RvXp1XL58WeooBVJS9geuX7+O3bt3Y/PmzVJHKRSp9y1K/TmrlpaWCAgIwJ49e7TTcnJysGfPHpM7n6K4CCEwatQobNmyBXv37kXlypWljiR7rVu3xpkzZxATE6N9BAYGol+/foiJiTG5QlccgoKC9G7PcPHiRVSsWFGiRPKXmpoKMzPdP+vm5ubIycmRKJHpqVy5MlxdXXVqQlJSEo4ePVpqawLrpLRKYs3NyclBRkaG1DEKpCTW75SUFFy5ckXnCzk5Kyn7A6tWrUL58uXRqVMnqaMUitT7FqX+yCoAhIWFITQ0FIGBgWjQoAHmz5+Px48fY9CgQVJHk6WRI0di3bp1+PHHH2FnZ6c9j8vBwQFqtVridPJkZ2end36RjY0NypYta1LnHRWn9957D02aNMGMGTPQq1cvHDt2DMuWLcOyZcukjiZbISEhmD59OipUqIBatWrh1KlTmDdvHgYPHix1NFlJSUnROaJw9epVxMTEoEyZMqhQoQLGjh2Lzz77DNWqVUPlypUxadIkuLu7o1u3btKFlpgp18kX/b7lztRr7oQJE9ChQwdUqFABycnJWLduHfbv34+dO3dKHa1ASkL9/uCDDxASEoKKFSvi9u3bmDx5MszNzdGnTx+poxVISdgfyMnJwapVqxAaGgoLC9NqvyTftyiWaw6bgIULF4oKFSoIS0tL0aBBA3HkyBGpI8kWgHwfq1atkjqaSTG1S99L4eeffxa+vr5CpVIJb29vsWzZMqkjyVpSUpIYM2aMqFChgrCyshJeXl5i4sSJIiMjQ+posrJv3758/4aFhoYKIXJvXzNp0iTh4uIiVCqVaN26tbhw4YK0oWXAVOvki37fcmfqNXfw4MGiYsWKwtLSUjg7O4vWrVuLXbt2SR2rSEytfvfu3Vu4ubkJS0tL8corr4jevXuLy5cvSx2rUEx9f2Dnzp0CgEnWEqn3LRRCCFE8bTERERERERFRwZT6c1aJiIiIiIhIftisEhERERERkeywWSUiIiIiIiLZYbNKREREREREssNmlYiIiIiIiGSHzSoRERERERHJDptVIiIiIiIikh02q0RERERERCQ7bFaJSMfAgQPRrVu3Ii1j//79UCgUSEhIMEomIiIiuWCdJCo+FlIHICJ5WbBgAYQQUscgIiKSJdZJouLDZpWIAAAajQYKhQIODg5SRyEiIpId1kmi4sdhwEQmqmXLlhg1ahRGjRoFBwcHlCtXDpMmTdJ+25uRkYEPPvgAr7zyCmxsbNCwYUPs379f+/7Vq1fD0dERP/30E3x8fKBSqXDjxg294U0ZGRkYPXo0ypcvDysrKzRt2hTR0dE6WbZv347q1atDrVajVatWuHbtWjFsASIiomdjnSQyfWxWiUzYmjVrYGFhgWPHjmHBggWYN28eVqxYAQAYNWoU/vjjD/zwww84ffo0evbsifbt2+PSpUva96empmLWrFlYsWIF/vrrL5QvX15vHR999BE2bdqENWvW4OTJk6hatSratWuHhw8fAgBu3ryJ1157DSEhIYiJicFbb72F8ePHF88GICIieg7WSSITJ4jIJLVo0ULUrFlT5OTkaKeNGzdO1KxZU1y/fl2Ym5uLW7du6byndevWYsKECUIIIVatWiUAiJiYGJ15QkNDRdeuXYUQQqSkpAilUim+//577euZmZnC3d1dzJ49WwghxIQJE4SPj4/OMsaNGycAiEePHhnr4xIRERUK6ySR6eM5q0QmrFGjRlAoFNrnjRs3xty5c3HmzBloNBpUr15dZ/6MjAyULVtW+9zS0hK1a9d+5vKvXLmCrKwsBAUFaacplUo0aNAA586dAwCcO3cODRs21Hlf48aNi/S5iIiIjIF1ksi0sVklKoFSUlJgbm6OEydOwNzcXOc1W1tb7c9qtVqniBMREZUGrJNEpoHnrBKZsKNHj+o8P3LkCKpVq4a6detCo9Hg7t27qFq1qs7D1dW1wMuvUqUKLC0tcejQIe20rKwsREdHw8fHBwBQs2ZNHDt2TC8HERGR1FgniUwbm1UiE3bjxg2EhYXhwoULWL9+PRYuXIgxY8agevXq6NevHwYMGIDNmzfj6tWrOHbsGMLDw/HLL78UePk2NjYYPnw4PvzwQ+zYsQNnz57F22+/jdTUVAwZMgQAMGzYMFy6dAkffvghLly4gHXr1mH16tUv6RMTEREVHOskkWnjMGAiEzZgwACkpaWhQYMGMDc3x5gxYzB06FAAwKpVq/DZZ5/h/fffx61bt1CuXDk0atQInTt3LtQ6Zs6ciZycHPTv3x/JyckIDAzEzp074eTkBACoUKECNm3ahPfeew8LFy5EgwYNMGPGDAwePNjon5eIiKgwWCeJTJtCiH9vNkVEJqVly5bw9/fH/PnzpY5CREQkO6yTRKaPw4CJiIiIiIhIdtisEhERERERkexwGDARERERERHJDo+sEhERERERkeywWSUiIiIiIiLZYbNKREREREREssNmlYiIiIiIiGSHzSoRERERERHJDptVIiIiIiIikh02q0RERERERCQ7bFaJiIiIiIhIdtisEhERERERkez8H3QBPXeb3+nvAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 950x360 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(9.5, 3.6))\n",
    "labels_a = [\"R (source)\", \"A\", \"B (bottleneck)\", \"S (sink)\"]\n",
    "for i, label in enumerate(labels_a):\n",
    "    axes[0].step(range(trajectory_a.occupancies.shape[0]), trajectory_a.occupancies[:, i],\n",
    "                 where=\"post\", label=label)\n",
    "axes[0].set_title(f\"corridor (cost {metrics_a['total_cost']:.0f})\")\n",
    "axes[0].set_xlabel(\"period\")\n",
    "axes[0].set_ylabel(\"occupancy\")\n",
    "axes[0].legend(fontsize=7)\n",
    "\n",
    "labels_b = [\"R\", \"A\", \"B (tiny)\", \"C\", \"D\", \"S\"]\n",
    "for i, label in enumerate(labels_b):\n",
    "    axes[1].step(range(trajectory_b.occupancies.shape[0]), trajectory_b.occupancies[:, i],\n",
    "                 where=\"post\", label=label)\n",
    "axes[1].set_title(f\"diverge+spillback (cost {metrics_b['total_cost']:.0f})\")\n",
    "axes[1].set_xlabel(\"period\")\n",
    "axes[1].legend(fontsize=7)\n",
    "fig.tight_layout()\n",
    "display(fig)\n",
    "plt.close(fig)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c8db3f3",
   "metadata": {},
   "source": [
    "## Takeaways & pointers\n",
    "\n",
    "- **Certified, not self-reported.** Both anchors' costs came from\n",
    "  `CellSODTAEvaluator`, which re-solves the LP itself and cross-checks the\n",
    "  dual certificate — never the solver's own objective claim.\n",
    "- **Finite storage is the point.** `merchant-nemhauser`'s exit functions have\n",
    "  no analogue of `N_B`; the 26-vs-25 gap recomputed above is real physics this\n",
    "  LP can represent that the exit-function model cannot.\n",
    "- **Where next.** the exit-function twin\n",
    "  [`merchant-nemhauser`](01-merchant-nemhauser.ipynb) (no finite storage); the\n",
    "  CTM cell scheme this LP relaxes [`ctm`](../05-dnl/01-ctm.ipynb); the lineage\n",
    "  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": "dta",
   "unit": "lp-so-dta"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}