{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "324ab578",
   "metadata": {},
   "source": [
    "# `merchant-nemhauser` — Merchant & Nemhauser's (1978) exit-function SO-DTA\n",
    "\n",
    "**What.** `merchant-nemhauser` is the first network dynamic traffic assignment\n",
    "model: a discrete-time, single-destination system-optimal program that routes\n",
    "a time-varying demand through links whose outflow is governed by per-link\n",
    "`exit functions` `e_a(t) <= g_a(x_a(t))`. This benchmark ships the Carey (1987)\n",
    "LP relaxation (exit as an inequality, not the original nonconvex equality),\n",
    "solved to global optimality by `scipy.optimize.linprog` (HiGHS).\n",
    "\n",
    "**Why it is in the benchmark.** It founds the field of network DTA — every\n",
    "later analytical-DTA model in this track (`vi-due`, `pm-td`, `lp-so-dta`)\n",
    "answers a question this paper first posed. Its distinctive result is that\n",
    "holding traffic back below its exit-function bound (ramp metering) can be\n",
    "STRICTLY optimal, which this notebook demonstrates on a second anchor. See the\n",
    "[model compendium](../../docs/MODELS.md) (Merchant & Nemhauser 1978) and\n",
    "[docs/design/adr-020-merchant-nemhauser.md](../../docs/design/adr-020-merchant-nemhauser.md)\n",
    "(P1).\n",
    "\n",
    "**Scope.** This notebook solves the LP on two built-in anchors and certifies\n",
    "both; it does not benchmark DTA formulations against each other — for the\n",
    "spillback-capable LP, see [`lp-so-dta`](../07-dta/02-lp-so-dta.ipynb).\n",
    "\n",
    "**Canon.** `[merchant1978model]`, `[merchant1978optimality]`, [docs/REFERENCES.md](../../docs/REFERENCES.md) / [docs/references.bib](../../docs/references.bib)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dfb8ce66",
   "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 `SODTAEvaluator` from the\n",
    "emitted inflow/exit/occupancy arrays alone — conservation, node balance, the\n",
    "exit-function bound, and a weak-duality backstop against the harness's own\n",
    "freshly-solved LP optimum `Z*`. The solver's self-reported objective is never\n",
    "trusted ([README](../../README.md), *Certified, not self-reported*)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "50783d73",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T13:48:52.549070Z",
     "iopub.status.busy": "2026-07-21T13:48:52.548593Z",
     "iopub.status.idle": "2026-07-21T13:48:54.581117Z",
     "shell.execute_reply": "2026-07-21T13:48:54.579677Z"
    }
   },
   "outputs": [],
   "source": [
    "# Setup. `merchant-nemhauser` is a core model: a plain `pip install -e .`\n",
    "# suffices — no optional extra, so no guard cell. The inline backend is\n",
    "# Agg-based (headless CI renders into the notebook); NEVER\n",
    "# matplotlib.use(\"Agg\") in-kernel — it silently suppresses inline capture.\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "\n",
    "from tabench import (\n",
    "    SODTAEvaluator,\n",
    "    mn_metering_scenario,\n",
    "    mn_parallel_scenario,\n",
    "    solve_so_dta,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df955103",
   "metadata": {},
   "source": [
    "## Anchor A: parallel-route capacity metering\n",
    "\n",
    "6 vehicles at node 0 choose between a fast capacitated route (`g = min(x, 2)`,\n",
    "1 period free-flow) and a slow uncapacitated one (2 periods). `SODTAScenario`\n",
    "is frozen and content-hashed (P2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a78d3d11",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T13:48:54.586132Z",
     "iopub.status.busy": "2026-07-21T13:48:54.585819Z",
     "iopub.status.idle": "2026-07-21T13:48:54.591981Z",
     "shell.execute_reply": "2026-07-21T13:48:54.591081Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "scenario      : mn-parallel\n",
      "content hash  : 9bbd85f7fd84617f…\n",
      "links=3  nodes=3  periods=5  destination=2\n",
      "demand at node 0, period 0 : 6.0\n"
     ]
    }
   ],
   "source": [
    "scenario_a = mn_parallel_scenario()\n",
    "print(f\"scenario      : {scenario_a.name}\")\n",
    "print(f\"content hash  : {scenario_a.content_hash()[:16]}…\")\n",
    "print(f\"links={scenario_a.n_links}  nodes={scenario_a.n_nodes}  \"\n",
    "      f\"periods={scenario_a.n_periods}  destination={scenario_a.destination}\")\n",
    "print(f\"demand at node 0, period 0 : {scenario_a.demand[0, 0]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d47739a",
   "metadata": {},
   "source": [
    "## Solve (Anchor A)\n",
    "\n",
    "No `Budget`/`RngBundle`/`Trace` — `solve_so_dta` builds the canonical LP\n",
    "(`canonical_lp`) and solves it to global optimality with HiGHS, emitting a\n",
    "`DTATrajectory` (per-period inflow/exit/occupancy arrays, the P1-certifiable\n",
    "artifact) plus an LP dual certificate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "fba2dd96",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T13:48:54.595651Z",
     "iopub.status.busy": "2026-07-21T13:48:54.595389Z",
     "iopub.status.idle": "2026-07-21T13:48:54.611774Z",
     "shell.execute_reply": "2026-07-21T13:48:54.610915Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "self-reported objective : 10.0000  (provenance only)\n",
      "exits (fast, slow)      : [[-0.0, -0.0, 0.0], [2.0, 2.0, -0.0], [2.0, 0.0, 2.0], [0.0, 0.0, -0.0], [-0.0, -0.0, -0.0]]\n"
     ]
    }
   ],
   "source": [
    "trajectory_a = solve_so_dta(scenario_a)\n",
    "print(f\"self-reported objective : {trajectory_a.provenance['objective']:.4f}  (provenance only)\")\n",
    "print(f\"exits (fast, slow)      : {trajectory_a.exits.round(3).tolist()}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc324bed",
   "metadata": {},
   "source": [
    "## Certify (P1) — global optimum, via LP duality\n",
    "\n",
    "`SODTAEvaluator` independently RE-SOLVES the canonical LP at construction (its\n",
    "own fresh `Z*`), recomputes conservation/node-balance/exit-bound feasibility\n",
    "from the emitted arrays, and verifies the trajectory's own dual certificate by\n",
    "pure arithmetic (weak duality + zero gap = global optimality) — never trusting\n",
    "either the solver's claim or its dual."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "0588cfa2",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T13:48:54.614701Z",
     "iopub.status.busy": "2026-07-21T13:48:54.614151Z",
     "iopub.status.idle": "2026-07-21T13:48:54.620644Z",
     "shell.execute_reply": "2026-07-21T13:48:54.619780Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "feasible            : 1\n",
      "so_optimality_gap    : 0.000e+00\n",
      "total_cost           : 10.0000\n",
      "dual_gap             : 0.000e+00\n",
      "dual_infeasibility   : 0.000e+00\n"
     ]
    }
   ],
   "source": [
    "evaluator_a = SODTAEvaluator(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\"dual_gap             : {metrics_a['dual_gap']:.3e}\")\n",
    "print(f\"dual_infeasibility   : {metrics_a['dual_infeasibility']:.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\"], 10.0, atol=1e-6)\n",
    "assert metrics_a[\"dual_gap\"] < 1e-6\n",
    "assert metrics_a[\"dual_infeasibility\"] < 1e-6"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "404c982c",
   "metadata": {},
   "source": [
    "## Anchor B: the strict-holding result\n",
    "\n",
    "A series bottleneck `0 -A-> 1 -B-> 2` where the downstream street is twice as\n",
    "costly per vehicle-period (`weight 2` vs `1`). The distinctive Merchant-\n",
    "Nemhauser result: the RELAXED (Carey) LP optimum strictly PREFERS holding\n",
    "traffic back on link A below its exit-function bound — metering it at rate 1\n",
    "while `g_A(x_A(1)) = 2` — over the naive equality form `e = g(x)`, which is\n",
    "decision-free here and costs strictly more."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e2fe1f03",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T13:48:54.622819Z",
     "iopub.status.busy": "2026-07-21T13:48:54.622551Z",
     "iopub.status.idle": "2026-07-21T13:48:54.631572Z",
     "shell.execute_reply": "2026-07-21T13:48:54.630725Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "feasible             : 1\n",
      "so_optimality_gap     : 0.000e+00\n",
      "total_cost (relaxed)  : 18.0000\n",
      "exit_slack_max        : 1.0000  (> 0 means SOME period held traffic below its exit bound)\n",
      "link A exits (rate-1 metered) : [-0.0, 1.0, 1.0, 1.0, 1.0, -0.0, -0.0]\n",
      "naive equality-form cost (no holding, recomputed) : 22.0000\n"
     ]
    }
   ],
   "source": [
    "scenario_b = mn_metering_scenario()\n",
    "trajectory_b = solve_so_dta(scenario_b)\n",
    "evaluator_b = SODTAEvaluator(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 (relaxed)  : {metrics_b['total_cost']:.4f}\")\n",
    "print(f\"exit_slack_max        : {metrics_b['exit_slack_max']:.4f}  \"\n",
    "      f\"(> 0 means SOME period held traffic below its exit bound)\")\n",
    "assert metrics_b[\"feasible\"] == 1.0\n",
    "assert metrics_b[\"so_optimality_gap\"] < 1e-6\n",
    "assert np.isclose(metrics_b[\"total_cost\"], 18.0, atol=1e-6)\n",
    "# The distinctive result: holding is used (exit bound not tight everywhere).\n",
    "assert metrics_b[\"exit_slack_max\"] > 0.5\n",
    "print(f\"link A exits (rate-1 metered) : {trajectory_b.exits[:, 0].round(3).tolist()}\")\n",
    "\n",
    "# Recompute the naive equality-form cost DIRECTLY from the scenario's exit\n",
    "# functions (never quoted): forcing e_a(t) = g_a(x_a(t)) exactly (no holding)\n",
    "# on this instance is decision-free -- one feasible trajectory -- and it costs\n",
    "# strictly more than the relaxed optimum.\n",
    "occ, cost = np.zeros(2), 0.0\n",
    "x = np.array([0.0, 0.0])\n",
    "demand = scenario_b.demand[:, 0]\n",
    "for t in range(scenario_b.n_periods):\n",
    "    g_a = min(x[0], 2.0)   # exit function g_A(x) = min(x, 2)\n",
    "    g_b = min(x[1], 1.0)   # exit function g_B(x) = min(x, 1)\n",
    "    x_next_a = x[0] + demand[t] - g_a\n",
    "    x_next_b = x[1] + g_a - g_b\n",
    "    cost += scenario_b.cost_weights[0] * x[0] + scenario_b.cost_weights[1] * x[1]\n",
    "    x = np.array([x_next_a, x_next_b])\n",
    "print(f\"naive equality-form cost (no holding, recomputed) : {cost:.4f}\")\n",
    "assert cost > metrics_b[\"total_cost\"] + 1.0  # strictly worse than the relaxed optimum"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3d406bb",
   "metadata": {},
   "source": [
    "## Visualize\n",
    "\n",
    "`tabench.viz` is a road-network flow visualizer; the M-N artifact is a\n",
    "per-link OCCUPANCY trajectory over discrete periods, not a static link-flow\n",
    "vector, so this notebook plots the certified occupancy trajectories directly\n",
    "(a house time-series plot, not `tabench.viz`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "afe93d41",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T13:48:54.634301Z",
     "iopub.status.busy": "2026-07-21T13:48:54.633982Z",
     "iopub.status.idle": "2026-07-21T13:48:54.871286Z",
     "shell.execute_reply": "2026-07-21T13:48:54.870601Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6sAAAFeCAYAAABjIokGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABYh0lEQVR4nO3deXxM9/7H8fcksiIRhNipPbbYJUrSUimK3Cq63aDo1UtL05Zy21q6hLaKW2q5tl6qlJZWKY0l9BatLbWVi1pbQovEElnP7w8/c02ziDHJTGZez8djHjLf8z3nfM7JmE8+Z/kek2EYhgAAAAAAcCBu9g4AAAAAAIA/o1gFAAAAADgcilUAAAAAgMOhWAUAAAAAOByKVQAAAACAw6FYBQAAAAA4HIpVAAAAAIDDoVgFAAAAADgcilUAAAAAgMOhWAUKWUREhBo2bGjvMGzis88+U+nSpXX16lV7h2I3f/zxh4oXL641a9bYOxQATozc4brGjh0rk8mk33//3WbLjIiIUERExB37xcfHy2QyKT4+3tzWr18/Va9e3WaxFFULFiyQyWTSiRMnzG1t2rTRiBEj7BeUE6JYBVxYZmamKlasKJPJpG+++eau5x0zZoyef/55lShRooAivLN33nlHK1euzHf/GTNmqFevXqpatapMJpP69euXa9/Lly/r2WefVWBgoIoXL64HHnhAu3fvtuhTpkwZDRw4UK+//rqVWwAARUP16tVlMpnML29vb9WuXVuvvPKKLl68mK9luELu2LVrlx555BEFBQWpRIkSaty4sf75z38qMzPz3gOHQxs5cqSmT5+uc+fO2TsUp0GxCriwjRs36uzZs6pevbo++eSTu5p31apVOnz4sJ599tkCii5/7vYPjokTJ2rjxo1q0KCBihUrlmu/rKwsde3aVYsXL9bQoUP17rvv6vz584qIiNCRI0cs+g4ePFi7d+/Wxo0brd0MACgSQkJCtHDhQi1cuFDTpk1Tx44dNWXKFD388MP5mt/Zc8euXbsUFhamEydOaOTIkZo0aZLuu+8+DRs2TDExMTaI3P7+9a9/6fDhw/YOwyH16NFDfn5++uijj+wditPI/X8bgCIrIyNDWVlZ8vT0zLPfokWL1KxZM/Xt21ejR4/WtWvXVLx48XytY/78+Wrbtq0qVapki5ALzebNm81HxvM6qr98+XJt3bpVy5Yt02OPPSZJ6t27t+rUqaMxY8Zo8eLF5r7169dXw4YNtWDBAj344IMFvg0AUBDykzsqVaqkp59+2vx+4MCBKlGihN5//30dOXJEtWvXznMdzp47Zs2aJUnasmWLSpcuLUn629/+pvDwcC1YsEBTp04tlHgLkoeHh71DcFhubm567LHH9O9//1vjxo2TyWSyd0hFHmdW4bJOnjypv//976pbt658fHxUpkwZ9erVy+LeA+l/9yR8//33iomJMV8S+pe//EUXLlzIttxvvvlG4eHhKlmypPz8/NSyZUuLwuaWgwcP6oEHHpCvr68qVaqkd999N1uf8+fPa8CAASpfvry8vb3VpEkTffzxxxZ9Tpw4IZPJpPfff19TpkxRzZo15eXlpYMHD+a5/SkpKVqxYoUef/xx9e7dWykpKfryyy/zseekGzduaO3aterYsWOO0xctWqRWrVrJ19dXAQEBat++vb799luLPh999JEaNGggLy8vVaxYUUOGDNHly5ct+hw5ckQ9e/ZUUFCQvL29VblyZT3++ONKSkqSJJlMJl27dk0ff/yx+bK0vC7NkqRq1arlK3ksX75c5cuX16OPPmpuCwwMVO/evfXll18qNTXVov9DDz2kVatWyTCMOy4bQNHl6rkjJ0FBQZKU5xlHyTVyR3Jysry9vVWqVCmL9goVKsjHx+eO8+fl8uXL6tevn0qVKiV/f3/1799f169ft+iTkZGhN9980/z7rF69ukaPHp0tZ+XkzJkzioqKUvHixVWuXDm9+OKLOc7353tWb/8szZ4927zuli1baseOHdnmX7ZsmYKDg+Xt7a2GDRtqxYoV+b4P1mQyaezYsdnaq1evbvE7TE9P17hx41S7dm15e3urTJkyuv/++xUXF2cx36FDh/TYY4+pdOnS8vb2VosWLfTVV19lW/6BAwf04IMPysfHR5UrV9Zbb72lrKysHGN86KGHdPLkSSUkJNxxe3BnnFmFy9qxY4e2bt2qxx9/XJUrV9aJEyc0Y8YMRURE6ODBg/L19bXo//zzzysgIEBjxozRiRMnNGXKFA0dOlRLly4191mwYIGeeeYZNWjQQKNGjVKpUqW0Z88erV27Vk8++aS536VLl/Twww/r0UcfVe/evbV8+XKNHDlSjRo1UufOnSXdLCYjIiJ09OhRDR06VDVq1NCyZcvUr18/Xb58WcOGDbOIb/78+bpx44aeffZZeXl5mY/o5uarr77S1atX9fjjjysoKEgRERH65JNPLOLMza5du5SWlqZmzZplmzZu3DiNHTtWYWFhGj9+vDw9PfXDDz9o48aN6tSpk6Sbg0WMGzdOHTt21HPPPafDhw9rxowZ2rFjh77//nt5eHgoLS1NkZGRSk1N1fPPP6+goCD9+uuv+vrrr3X58mX5+/tr4cKFGjhwoFq1amW+pKxmzZp3jD8/9uzZo2bNmsnNzfKYXqtWrTR79mz997//VaNGjcztzZs31+TJk3XgwAGnGQQFQHaunjvS09PNA/3cuHFDe/bs0QcffKD27durRo0aec7rCrkjIiJCS5cu1d/+9jfFxMTI19dX33zzjb744gu9995797Ts3r17q0aNGoqNjdXu3bs1Z84clStXThMnTjT3GThwoD7++GM99thjeumll/TDDz8oNjZWP//8s1asWJHrslNSUtShQwedOnVKL7zwgipWrKiFCxfe1e0tixcv1pUrV/S3v/1NJpNJ7777rh599FH98ssv5rOxq1evVp8+fdSoUSPFxsbq0qVLGjBggM3PtI8dO1axsbHm33NycrJ27typ3bt366GHHpJ0swC9dZb/1VdfVfHixfXZZ58pKipKn3/+uf7yl79Iks6dO6cHHnhAGRkZ5n6zZ8/O9eBD8+bNJUnff/+9mjZtatPtckkG4KKuX7+erW3btm2GJOPf//63uW3+/PmGJKNjx45GVlaWuf3FF1803N3djcuXLxuGYRiXL182SpYsabRu3dpISUmxWO7t84WHh2dbR2pqqhEUFGT07NnT3DZlyhRDkrFo0SJzW1pamhEaGmqUKFHCSE5ONgzDMI4fP25IMvz8/Izz58/ne/sfeeQRo23btub3s2fPNooVK5avZcyZM8eQZOzbt8+i/ciRI4abm5vxl7/8xcjMzLSYdmsfnD9/3vD09DQ6depk0WfatGmGJGPevHmGYRjGnj17DEnGsmXL8oylePHiRt++fe8Y893OW7x4ceOZZ57J1r569WpDkrF27VqL9q1btxqSjKVLl1oVC4CiwZVzR7Vq1QxJ2V5t27Y1fv/99zvO7wq5IyMjwxg6dKjh4eFh3j/u7u7GjBkzrFqXYRjGmDFjDEnZctJf/vIXo0yZMub3CQkJhiRj4MCBFv1efvllQ5KxceNGc1t4eLgRHh5ufn/rc/PZZ5+Z265du2bUqlXLkGRs2rTJ3N63b1+jWrVq5ve3PktlypQxLl68aG7/8ssvDUnGqlWrzG2NGjUyKleubFy5csXcFh8fb0iyWGZuJBljxozJ1l6tWjWL30mTJk2Mrl275rmsDh06GI0aNTJu3LhhbsvKyjLCwsKM2rVrm9uGDx9uSDJ++OEHc9v58+cNf39/Q5Jx/PjxbMv29PQ0nnvuuTtuD+6My4Dhsm4/Ipaenq4//vhDtWrVUqlSpbKN+CpJzz77rMUlQO3atVNmZqZOnjwpSYqLi9OVK1f06quvytvb22LeP186VKJECYt7fjw9PdWqVSv98ssv5rY1a9YoKChITzzxhLnNw8NDL7zwgq5evarNmzdbLLNnz54KDAzM17b/8ccfWrduncWye/bsKZPJpM8++yxf80tSQECARfvKlSuVlZWlN954I9sZyVv7YP369UpLS9Pw4cMt+gwaNEh+fn5avXq1JMnf31+StG7dumyXORWGlJQUeXl5ZWu/9btNSUmxaL+1L2z5aAEAjseVc4cktW7dWnFxcYqLi9PXX3+tt99+WwcOHFD37t2zfS/+mSvkDnd3d9WsWVORkZH6+OOPtXTpUnXr1k3PP//8XQ3olJPBgwdbvG/Xrp3++OMPJScnS5L5EWp/HsjppZdekiTzPsrJmjVrVKFCBfMYDZLk6+t7VwNh9enTx+J3265dO0kyfz5/++037du3T9HR0Rb3/YaHh1tcqWQLpUqV0oEDB7INiHjLxYsXtXHjRvXu3VtXrlzR77//rt9//11//PGHIiMjdeTIEf3666+Sbu6bNm3aqFWrVub5AwMD9dRTT+W6/oCAAP4esBGKVbislJQUvfHGG6pSpYq8vLxUtmxZBQYG6vLly+b7Wm5XtWpVi/e3vpAvXbokSTp27Jgk5esS0MqVK2f7IyQgIMC8LOnmfVG1a9fOlrjr169vnn67O11+dbulS5cqPT1dTZs21dGjR3X06FFdvHhRrVu3vqtRgY0/3Z957Ngxubm5KTg4ONd5bsVdt25di3ZPT0/dd9995uk1atRQTEyM5syZo7JlyyoyMlLTp0/P8XdTEHx8fHK8V+fGjRvm6be7tS8YTAFwbq6cOySpbNmy6tixozp27KiuXbtq9OjRmjNnjrZu3ao5c+bkaxnOnDsmTJigiRMn6tNPP1V0dLR69+6tFStW6P7779eQIUOUkZFh9bLv9Fk6efKk3NzcVKtWLYt+QUFBKlWqVLbf/e1OnjypWrVqZft8/Xl/32t8krLFl1vbvRg/frwuX76sOnXqqFGjRnrllVe0d+9e8/SjR4/KMAy9/vrrCgwMtHiNGTNG0s17v2/FndPAYXntG8Mw+HvARihW4bKef/55vf322+rdu7c+++wzffvtt4qLi1OZMmVyvGne3d09x+X8Oenmhy2XdcvdDNxwqyBt27atateubX795z//0bZt2yyO0uekTJkykmTxB1JBmDRpkvbu3avRo0crJSVFL7zwgho0aKAzZ84U6Hqlm4NhnD17Nlv7rbaKFStatN/aF2XLli3w2ADYjyvnjtx06NBB0s0RcPPiCrnjo48+0oMPPphtxODu3bvrt99+yzYQ193I7+/fXkVSQXw+8+vPz7Bt3769jh07pnnz5qlhw4aaM2eOmjVrZj6gcuv/6ssvv2y+UuDPr3spoC9fvszfAzbCAEtwWcuXL1ffvn01adIkc9uNGzeyjSqYX7cGZ9i/f79NjhBWq1ZNe/fuVVZWlsUR8kOHDpmnW+P48ePaunWrhg4dqvDwcItpWVlZ+utf/6rFixfrtddey3UZ9erVMy/r9kt3atasqaysLB08eFAhISG5bpckHT58WPfdd5+5PS0tTcePH882SmSjRo3UqFEjvfbaa9q6davatm2rmTNn6q233pJUcEk5JCRE3333Xbb9/8MPP8jX11d16tSx6H/8+HFJ/zt7AcA5uWruyMuts4VXr17Ns58r5I7ExMRshZN085JxSfd0ZvVOqlWrpqysLB05csQiFyUmJury5ct5/u6rVaum/fv3ZzsjaMvnqd5a/9GjR7NNy6ktJwEBAdn+r6WlpeV4cLl06dLq37+/+vfvr6tXr6p9+/YaO3asBg4caP4MeXh45Do69e1x53Q5cW775tdff1VaWhp/D9gIZ1bhstzd3bMd7fvwww9zTDL50alTJ5UsWVKxsbHmS0VvseaoYpcuXXTu3DmLESMzMjL04YcfqkSJEtkKzfy6dVZ1xIgReuyxxyxevXv3Vnh4+B0vBW7evLk8PT21c+dOi/aoqCi5ublp/Pjx2c4w3NoHHTt2lKenp/75z39a7Je5c+cqKSlJXbt2lXRz+P8/J/VGjRrJzc3N4vLc4sWLW/1HYl4ee+wxJSYm6osvvjC3/f7771q2bJm6deuW7X7WXbt2yd/fXw0aNLB5LAAch6vmjrysWrVKktSkSZM8+7lC7qhTp47i4uLM9+dKN8/6ffbZZypZsqTNRh3OSZcuXSRJU6ZMsWj/4IMPJMm8j3Kb97ffftPy5cvNbdevX9fs2bNtFl/FihXVsGFD/fvf/7Y4sLF582bt27cvX8uoWbNmtjP4s2fPzvb/7/b9L92837tWrVrmz0C5cuUUERGhWbNm5Vjo3v54qS5dumj79u368ccfLabn9rfSrl27JElhYWH52ibkjTOrcFmPPPKIFi5cKH9/fwUHB2vbtm1av369+TKlu+Xn56fJkydr4MCBatmypZ588kkFBATop59+0vXr17M94+5Onn32Wc2aNUv9+vXTrl27VL16dS1fvlzff/+9pkyZopIlS1oV5yeffKKQkBBVqVIlx+ndu3fX888/r927d+f4eAHp5iBDnTp10vr16zV+/Hhze61atfSPf/xDb775ptq1a6dHH31UXl5e2rFjhypWrKjY2FgFBgZq1KhRGjdunB5++GF1795dhw8f1kcffaSWLVuaBw/ZuHGjhg4dql69eqlOnTrKyMjQwoUL5e7urp49e5rX2bx5c61fv14ffPCBKlasqBo1aqh169a5bv+qVav0008/Sbp5pHvv3r3mI+3du3dX48aNJd0sVtu0aaP+/fvr4MGDKlu2rD766CNlZmZq3Lhx2ZYbFxenbt26cY8K4ORcNXfc8uuvv2rRokWSbp7R+umnnzRr1iyVLVtWzz//fJ7zukLuePXVV/X000+rdevWevbZZ+Xj46NPP/1Uu3bt0ltvvWV+hIt083mlH3/8sY4fP56vZ4zeSZMmTdS3b1/Nnj1bly9fVnh4uH788Ud9/PHHioqK0gMPPJDrvIMGDdK0adMUHR2tXbt2qUKFClq4cGG2RzHdq3feeUc9evRQ27Zt1b9/f126dEnTpk1Tw4YN73hmXrr5aJ7BgwerZ8+eeuihh/TTTz9p3bp12S65DQ4OVkREhJo3b67SpUtr586dWr58uYYOHWruM336dN1///1q1KiRBg0apPvuu0+JiYnatm2bzpw5Y/59jxgxQgsXLtTDDz+sYcOGmR9dc+sqhj+Li4tT1apVeWyNrRTu4MOA47h06ZLRv39/o2zZskaJEiWMyMhI49ChQ9mGP7/1+IEdO3ZYzL9p06Zsw7kbhmF89dVXRlhYmOHj42P4+fkZrVq1Mj799FPz9PDwcKNBgwbZ4vnzUPCGYRiJiYnmGD09PY1GjRoZ8+fPt+hza8j49957747bvGvXLkOS8frrr+fa58SJE4Yk48UXX8xzWV988YVhMpmMU6dOZZs2b948o2nTpoaXl5cREBBghIeHG3FxcRZ9pk2bZtSrV8/w8PAwypcvbzz33HPGpUuXzNN/+eUX45lnnjFq1qxpeHt7G6VLlzYeeOABY/369RbLOXTokNG+fXvDx8fHkHTHRxH07ds3x0cvSMq2by9evGgMGDDAKFOmjOHr62uEh4dn+xwYhmH8/PPPhqRssQFwPq6YO27586Nr3NzcjHLlyhlPPPGEcfTo0XwtwxVyx9q1a43w8HCL/T9z5sxsy+zZs6fh4+NjEX9Obj265sKFCxbttz5jtz86JT093Rg3bpxRo0YNw8PDw6hSpYoxatQoi8ezGEb2R9cYhmGcPHnS6N69u+Hr62uULVvWGDZsmLF27dp8P7omp8+ScnjUzJIlS4x69eoZXl5eRsOGDY2vvvrK6Nmzp1GvXr0894NhGEZmZqYxcuRIo2zZsoavr68RGRlpHD16NNv/v7feesto1aqVUapUKcPHx8eoV6+e8fbbbxtpaWkWyzt27JgRHR1tBAUFGR4eHkalSpWMRx55xFi+fLlFv7179xrh4eGGt7e3UalSJePNN9805s6dm23/Z2ZmGhUqVDBee+21O24L8sdkGIVw1zMAp5OZmang4GD17t1bb775pr3Dsavhw4dry5Yt2rVrF2dWASAP5I7/KV++vKKjo/Xee+/ZOxS7CwkJUWBgoOLi4uwdyj1ZuXKlnnzySR07dkwVKlSwdzhOgXtWAVjF3d1d48eP1/Tp0/N16Y6z+uOPPzRnzhy99dZbFKoAcAfkjpsOHDiglJQUjRw50t6hFKr09PRs9xTHx8frp59+UkREhH2CsqGJEydq6NChFKo2xJlVAAAAAAXuxIkT6tixo55++mlVrFhRhw4d0syZM+Xv76/9+/dbfe83nBcDLAEAAAAocAEBAWrevLnmzJmjCxcuqHjx4uratasmTJhAoYoccWYVAAAAAOBwuGcVAAAAAOBwKFYBAAAAAA7H5e5ZzcrK0m+//aaSJUsycicAoEAYhqErV66oYsWKcnMrWseFyZMAgIKW3zzpcsXqb7/9pipVqtg7DACACzh9+rQqV65s7zDuCnkSAFBY7pQnXa5YLVmypKSbO8bPz8/O0QAAnFFycrKqVKlizjlFCXkSAFDQ8psnXa5YvXVJk5+fH0kYAFCgiuJltORJAEBhuVOeLFo30gAAAAAAXALFKgAAAADA4bjcZcAAAOtlZmYqPT3d3mE4DHd3dxUrVqxIXu4LAIWBvOHaPDw85O7ubvX8FKsAgHy5evWqzpw5I8Mw7B2KQ/H19VWFChXk6elp71AAwKGQN2AymVS5cmWVKFHCqvkpVgEAd5SZmakzZ87I19dXgYGBnEnUzWfEpaWl6cKFCzp+/Lhq165d5J6pCgAFhbwBwzB04cIFnTlzRrVr17bqDKvDFKsTJkzQqFGjNGzYME2ZMiXXfsuWLdPrr7+uEydOqHbt2po4caK6dOlSeIECgAtKT0+XYRgKDAyUj4+PvcNxGD4+PvLw8NDJkyeVlpYmb2/vAlsXeRJAUULegCQFBgbqxIkTSk9Pt6pYdYhDwDt27NCsWbPUuHHjPPtt3bpVTzzxhAYMGKA9e/YoKipKUVFR2r9/fyFFCgCujSPj2RXG2VTyJICiirzh2u7192/3M6tXr17VU089pX/9619666238uw7depUPfzww3rllVckSW+++abi4uI0bdo0zZw5szDChZUMw1BKeqa9wygSfDzc+WIH8slkMunSpUsqVaqUunTposmTJ6tu3bp5ztOvXz+FhIRo+PDhefbLysrSsGHDtGbNGplMJg0fPlxDhw61YfT5Q550LM6az8g9AByR3YvVIUOGqGvXrurYseMdk/C2bdsUExNj0RYZGamVK1fmOk9qaqpSU1PN75OTk+8pXtw9wzD02Mxt2nXykr1DKRJaVAvQssGh/NEA3KU1a9bYdHmLFi3SwYMH9d///ldJSUlq2rSpHnjgATVo0MCm67kT8qTjcOZ8Ru6BKyjIA5y39O3bVytWrNDZs2dVvHjxXPv99NNPGj16tFavXn03m5CnhIQEHTp0SI8//rjNlpmbFi1a6P3331dERIRefvllNWvWTE8++aTN12PXy4CXLFmi3bt3KzY2Nl/9z507p/Lly1u0lS9fXufOnct1ntjYWPn7+5tfVapUuaeYcfdS0jOdMrEXlJ0nLznlUXs4F8MwdD0to0Bfdzt6ZPXq1ZWQkCBJ5uTZrl071axZU4MHD85xnu+++07BwcHauXNntmlLly7VoEGD5O7urtKlS6tPnz769NNP73pf3QvypGNx5nxG7kFBcsScsWbNmjsWqncrOTlZq1atUpMmTbRs2bI8+44aNUqvvvqqTdefkJCgJUuWWDVvRkaG1esdMWKExo4dq8xM23+H2O3M6unTpzVs2DDFxcUV6IAUo0aNsjjKnJycTCK2o52vdZSvp/XPWnJm19My1eKt9fYOA8iXlPRMBb+xrkDXcXB8pHw9rU9Tx44d06ZNm5Senq7g4GBt27ZNoaGh5ulLly5VbGysVq9erRo1amSb/9SpU6pWrZr5ffXq1bV9+3ar47lb5EnH5iz5jNyDwuCIOaN69epauXKlQkJCFBERoRYtWuiHH37Qb7/9poceeijHWye+++47/e1vf9O///1vtWjRItv0Tz/9VB07dtQTTzyhDz74QP369ctx3adOndKBAwfUrl07SdKJEycUEhKi559/XqtXr9aVK1e0YMECLV++XJs2bVJGRoaWLFmihg0bSpIWLlyoadOmKT09XSVKlNCHH36oChUq6I033lBSUpJCQkLUpk0bzZw5Uzt27NDIkSOVnJyszMxMjR49Wr169TKv829/+5vi4uIUHR2txx9/XC+88IJOnDihlJQU9ejRw3xFz9atW/X3v/9dGRkZatmypUVxW65cOdWsWVPffvutOnfunO/fQX7YrVjdtWuXzp8/r2bNmpnbMjMztWXLFk2bNk2pqanZRowKCgpSYmKiRVtiYqKCgoJyXY+Xl5e8vLxsGzys5uvpfk9/fAJAfvXp00fFihVTsWLFFBISomPHjpmL1YULF8rd3V2bNm1SQECAnSPNGXnSsZHPAOdyrwc4JWnu3LkaP368OnTooOeee06HDx/O8ezt5s2b1bJlS4u2pKQkNW/eXG+++abmzp2ryMhIrVq1SpMnT9Z7772ncePGadmyZfr+++/16aefasuWLfLy8tJ3332nJ598UgcOHND48eO1cuVK860fly9f1rPPPqs1a9aoQoUK+v3339WsWTOFhYWZ19mgQQNNnDhR0s3bRkaPHq3w8HBlZGTokUce0bJly9SjRw/16dNH8+fPV8eOHfXtt99qwYIFFvGHhoZqw4YNzlOsdujQQfv27bNo69+/v+rVq6eRI0fmOLTxrZ1w+zXjcXFxFh8kAEDB8/Fw18HxkQW+jntx+9lId3d3i6PAjRs31nfffad9+/apffv2Oc5ftWpVnTx50pxjTpw4oapVq95TTHeDPAnAWRSFnHGvBzj37duns2fPqlOnTnJzc9PTTz+tefPmmQvB2505cybbLRve3t6KioqSdPN+0BIlSuiBBx6QJLVq1UqffPKJJOnLL7/UTz/9pNatW5vnvXjxolJSUrKtZ+vWrfrll1+yFZCHDx/WfffdJw8PDz399NOSpGvXrmnDhg0WBzyvXr2qw4cP69ChQypWrJg6duwoSerUqZPuu+8+i2UGBQXp4MGDOe6be2G3YrVkyZLmU9m3FC9eXGXKlDG3R0dHq1KlSuZ7dYYNG6bw8HBNmjRJXbt21ZIlS7Rz507Nnj270OMHAFdmMpmK9FmlJk2aaMSIEerWrZumTZumhx9+OFufXr166V//+pd69eqlpKQkLV26VF9//XWhxUieBOAsikLOuNcDnHPnztWVK1fMRVx6erqysrL09ttvq1gxy2339fXVjRs3LNpuv8LF3d0913gMw1Dfvn31zjvv3HGbDMNQgwYNtHXr1mzTTpw4IV9fX/Pj127d87t9+/Zst57s3bs32/x/Hoztxo0bBfI8XYd4zmpuTp06pbNnz5rfh4WFafHixZo9e7aaNGmi5cuXa+XKldmSOQAAd1K/fn2tW7dOw4YN0+eff55t+l//+lfVq1dPtWvXVsuWLRUTE6NGjRrZIdLckScBoOA1adJEq1at0jPPPKO1a9dmm56WlqZFixZp+/btOnHihE6cOKFff/1VVatWzXG038aNG+vw4cNWxdK9e3ctWrRIp06dknTzMWu3Bgn08/NTUlKSuW9YWJiOHz+u9ev/d196QkKC0tLSsi331pncCRMmmNt+++03nTlzRvXq1VNGRoY2bdokSVq/fr2OHTtmMf/PP/+sJk2aWLVNeXGoQxzx8fF5vpduHunu1atX4QQEAHBot4/+eOLECfPPf84fy5cvN/98+302NWvWzPUPBnd3d02fPt0mcdoKeRIA7OPWAc4uXbronXfeUc+ePc3TVq5cqWrVqqlevXoW8zz11FOaO3euevToYdF+//3368yZM7p48aJKly59V3G0a9dO7777rv7yl78oIyNDaWlp6tq1q1q0aKEOHTro/fffV+PGjRUWFqaZM2dq9erVevnll/XSSy8pPT1dVatWzfVxZp988oliYmLUsGFDmUwmFS9eXLNmzVLlypW1dOlS/f3vf1dmZqZatmxpUZgahqENGzbYfHRjSTIZdzvOcxGXnJwsf39/JSUlyc/Pz97huITraRnmEeDudXRPZ8Z+giO7ceOGjh8/rho1ahToyLRFUU77pijnmqIce0Fzxu9pZ9wmOAbyRt7ee+89SdIrr7xi50ju3dq1a7Vo0SItWrQo27TcPgf5zTUOfRkwAAAAADibYcOGqUSJEvYOwyaSkpL07rvvFsiyOXwGAAAAAIXI09NTzz33nL3DsIk+ffoU2LI5swoAAAAAcDgUqwAAAAAAh0OxCgAAAABwOBSrAAAAAACHQ7EKACiyTCaTLl++LEnq0qVLvh6y3q9fP02ZMuWO/VavXq3mzZvLy8tLw4cPv7dAAQB2V5A5Y+zYsQoMDFRISIiaNGmili1bauvWrfcYMRgNGADgFNasWWPT5dWuXVvz5s3TsmXLdPXqVZsuGwBgX7bOGZL01FNPmQvbJUuWaNiwYdqxY4fN1+NKOLMKALh7hiGlXSvYl2HcVUjVq1dXQkKCJCkiIkIvv/yy2rVrp5o1a2rw4ME5zvPdd98pODhYO3fuzDatTp06atKkiYoV47guANwTF8gZf5aUlKSAgIC7ignZkYEBAHcv/br0TsWCXcfo3yTP4lbPfuzYMW3atEnp6ekKDg7Wtm3bFBoaap6+dOlSxcbGavXq1apRo4YtIgYA5MRFcsYnn3yi+Ph4JSUlKTk5WevWrbM6HtxEsQoAcEp9+vRRsWLFVKxYMYWEhOjYsWPmPzwWLlwod3d3bdq0iSPfAACb5IzbLwPesGGDHn30UR0+fFg+Pj6FsQlOiWIVAHD3PHxvHsUu6HXcA29vb/PP7u7uysjIML9v3LixvvvuO+3bt0/t27e/p/UAAO7ABXNGhw4ddOPGDe3fv18tW7a8p9hcGcUqAODumUz3dLmVvTVp0kQjRoxQt27dNG3aND388MP2DgkAnJcL5oyffvpJV69eVfXq1Qs+QCfGAEsAAJdUv359rVu3TsOGDdPnn3+ebfqGDRtUuXJlffDBB5o7d64qV66sr776yg6RAgDs7U45Q7p5z+qtR9f07dtXCxcuVGBgYCFH6lw4swoAKLKM20Z/PHHihPnn+Ph4i37Lly83/7xgwQLzzzVr1sz1OXsdOnTQmTNnbBInAMD+CjJnjB07VmPHjrVFmLgNZ1YBAAAAAA6HYhUAAAAA4HAoVgEAAAAADseuxeqMGTPUuHFj+fn5yc/PT6Ghofrmm29y7b9gwQKZTCaL1+3DTAMA4CzIkQCcwe33icL13Ovv364DLFWuXFkTJkxQ7dq1ZRiGPv74Y/Xo0UN79uxRgwYNcpzHz8/P4sZmk8lUWOECAFBoyJEAijIPDw+ZTCZduHBBgYGBfB+5IMMwdOHCBZlMJnl4eFi1DLsWq926dbN4//bbb2vGjBnavn17ronYZDIpKCioMMIDAMBuyJEAijJ3d3dVrlxZZ86csRh5F67FZDKpcuXKcnd3t2p+h3l0TWZmppYtW6Zr164pNDQ0135Xr15VtWrVlJWVpWbNmumdd97JNWlLUmpqqlJTU83vk5OTbRo3AMB+TCaTLl26pFKlSqlLly6aPHmy6tatm+c8/fr1U0hIiIYPH55nv3/+85+aPXu2+ZLaESNG6Omnn7Zh9PlXUDlSIk8CKDglSpRQ7dq1lZ6ebu9QYCceHh5WF6qSAxSr+/btU2hoqG7cuKESJUpoxYoVCg4OzrFv3bp1NW/ePDVu3FhJSUl6//33FRYWpgMHDqhy5co5zhMbG6tx48YV5CYAABzAmjVrbLq8Bg0a6Pvvv5e/v79Onz6tpk2bKjQ0VDVr1rTpevJS0DlSIk8CKFju7u73VKzAtdl9NOC6desqISFBP/zwg5577jn17dtXBw8ezLFvaGiooqOjFRISovDwcH3xxRcKDAzUrFmzcl3+qFGjlJSUZH6dPn26oDYFAGBH1atXV0JCgiQpIiJCL7/8stq1a6eaNWtq8ODBOc7z3XffKTg4WDt37sw2rUOHDvL395ckValSRUFBQYWeQwo6R0rkSQCA47L7mVVPT0/VqlVLktS8eXPt2LFDU6dOvWNylW6eVm7atKmOHj2aax8vLy95eXnZLF4AwM1BE1IyUgp0HT7FfO5pQI5jx45p06ZNSk9PV3BwsLZt22ZxCe3SpUsVGxur1atXq0aNGnkua/369bp06ZJatmxpdTzWKOgcKZEnAQCOy+7F6p9lZWVZ3DuTl8zMTO3bt09dunQp4KgAALdLyUhR68WtC3QdPzz5g3w9fK2ev0+fPipWrJiKFSumkJAQHTt2zFysLly4UO7u7tq0aZMCAgLyXM6+ffvUv39/LV26VMWLF7c6HlsgRwIAXIldi9VRo0apc+fOqlq1qq5cuaLFixcrPj5e69atkyRFR0erUqVKio2NlSSNHz9ebdq0Ua1atXT58mW99957OnnypAYOHGjPzQAAOKDbnzHq7u6ujIwM8/vGjRvru+++0759+9S+fftcl3Hw4EE98sgjmjdvnu6///4CjffPyJEAAFdn12L1/Pnzio6O1tmzZ+Xv76/GjRtr3bp1euihhyRJp06dkpvb/26rvXTpkgYNGqRz584pICBAzZs319atW3MdbAIAUDB8ivnohyd/KPB1FJQmTZpoxIgR6tatm6ZNm6aHH344W5+ff/5ZXbp00ezZs815qTCRIwEArs6uxercuXPznB4fH2/xfvLkyZo8eXIBRgQAyA+TyXRPl+g6gvr162vdunXq0qWL3nnnHfXs2dNi+gsvvKCkpCSNHDlSI0eOlCRNnDhRkZGRhRIfORIA4Ooc7p5VAADyyzAM88+3P3T+z4Xc8uXLzT8vWLDA/HPNmjV1+PDhHJcdFxdnkxgBAIB17P7oGgAAAAAA/oxiFQAAAADgcChWAQAAAAAOh2IVAAAAAOBwKFYBAPl2+4BGuCkrK8veIQAA4JQYDRgAcEceHh4ymUy6cOGCAgMDZTKZ7B2S3RmGobS0NF24cEFubm7y9PS0d0gAADgVilUAwB25u7urcuXKOnPmjMUjYiD5+vqqatWqcnPjYiUAAGyJYhUAkC8lSpRQ7dq1lZ6ebu9QHIa7u7uKFSvGmWYAAAoAxSoAIN/c3d3l7u5u7zAAAIAL4JolAAAAAIDDoVgFAAAAADgcilUAAAAAgMOhWAUAAAAAOByKVQAAAACAw6FYBQAAAAA4HIpVAAAAAIDDoVgFAAAAADgcilUAAAAAgMOxa7E6Y8YMNW7cWH5+fvLz81NoaKi++eabPOdZtmyZ6tWrJ29vbzVq1Ehr1qwppGgBACg85EgAgKuza7FauXJlTZgwQbt27dLOnTv14IMPqkePHjpw4ECO/bdu3aonnnhCAwYM0J49exQVFaWoqCjt37+/kCMHAKBgkSMBAK7OZBiGYe8gble6dGm99957GjBgQLZpffr00bVr1/T111+b29q0aaOQkBDNnDkzX8tPTk6Wv7+/kpKS5OfnZ7O4kbvraRkKfmOdJOng+Ej5ehazc0SOif0EOI+CyjUFnSMl8mRenPF7+vZt2vlaR/l6uts5Itvy8XCXyWSydxgA/iS/ucZhvmUzMzO1bNkyXbt2TaGhoTn22bZtm2JiYizaIiMjtXLlylyXm5qaqtTUVPP75ORkm8QLAEBhKagcKZEn8T8t3lpv7xBsrkW1AC0bHErBChRRdh9gad++fSpRooS8vLw0ePBgrVixQsHBwTn2PXfunMqXL2/RVr58eZ07dy7X5cfGxsrf39/8qlKlik3jBwCgoBR0jpTIk67Ox8NdLaoF2DuMArPz5CWlpGfaOwwAVrL7mdW6desqISFBSUlJWr58ufr27avNmzfnmozv1qhRoyyONCcnJ5OIAQBFQkHnSIk86epMJpOWDQ51uoLuelqmU54pBlyN3YtVT09P1apVS5LUvHlz7dixQ1OnTtWsWbOy9Q0KClJiYqJFW2JiooKCgnJdvpeXl7y8vGwbNAAAhaCgc6REnsTNgtUZ7r8F4Hzsfhnwn2VlZVncO3O70NBQbdiwwaItLi4u1/t3AABwJuRIAIArsethtFGjRqlz586qWrWqrly5osWLFys+Pl7r1t0clS46OlqVKlVSbGysJGnYsGEKDw/XpEmT1LVrVy1ZskQ7d+7U7Nmz7bkZAADYHDkSAODq7Fqsnj9/XtHR0Tp79qz8/f3VuHFjrVu3Tg899JAk6dSpU3Jz+9/J37CwMC1evFivvfaaRo8erdq1a2vlypVq2LChvTYBAIACQY4EALg6uxarc+fOzXN6fHx8trZevXqpV69eBRQRAACOgRwJAHB1DnfPKgAAAAAAFKsAAAAAAIdDsQoAAAAAcDgUqwAAAAAAh0OxCgAAAABwOBSrAAAAAACHQ7EKAAAAAHA4FKsAAAAAAIdDsQoAAAAAcDgUqwAAAAAAh0OxCgAAAABwOBSrAAAAAACHQ7EKAAAAAHA4FKsAAAAAAIdDsQoAAAAAcDgUqwAAAAAAh0OxCgAAAABwOBSrAAAAAACHY1Wxeu3aNVvHAQAAAACAmVXFavny5fXMM8/oP//5zz2tPDY2Vi1btlTJkiVVrlw5RUVF6fDhw3nOs2DBAplMJouXt7f3PcUBAICjIUcCAFydVcXqokWLdPHiRT344IOqU6eOJkyYoN9+++2ul7N582YNGTJE27dvV1xcnNLT09WpU6c7nrn18/PT2bNnza+TJ09asxkAADgsciQAwNUVs2amqKgoRUVF6cKFC1q4cKEWLFig119/XZGRkXrmmWfUvXt3FSt250WvXbvW4v2CBQtUrlw57dq1S+3bt891PpPJpKCgIGtCBwCgSCBHAgBc3T0NsBQYGKiYmBjt3btXH3zwgdavX6/HHntMFStW1BtvvKHr16/f1fKSkpIkSaVLl86z39WrV1WtWjVVqVJFPXr00IEDB3Ltm5qaquTkZIsXAABFTUHkSIk8CQBwXPdUrCYmJurdd99VcHCwXn31VT322GPasGGDJk2apC+++EJRUVH5XlZWVpaGDx+utm3bqmHDhrn2q1u3rubNm6cvv/xSixYtUlZWlsLCwnTmzJkc+8fGxsrf39/8qlKlyt1uJgAAdlVQOVIiTwIAHJdVlwF/8cUXmj9/vtatW6fg4GD9/e9/19NPP61SpUqZ+4SFhal+/fr5XuaQIUO0f//+Ow7aFBoaqtDQ0GzrmTVrlt58881s/UeNGqWYmBjz++TkZBIxAKBIKagcKZEnAQCOy6pitX///nr88cf1/fffq2XLljn2qVixov7xj3/ka3lDhw7V119/rS1btqhy5cp3FYuHh4eaNm2qo0eP5jjdy8tLXl5ed7VMAAAcRUHmSIk8CQBwXFYVq2fPnpWvr2+efXx8fDRmzJg8+xiGoeeff14rVqxQfHy8atSocdexZGZmat++ferSpctdzwsAgKMiRwIAXJ1VxWp8fLzc3d0VGRlp0b5u3TplZWWpc+fO+VrOkCFDtHjxYn355ZcqWbKkzp07J0ny9/eXj4+PJCk6OlqVKlVSbGysJGn8+PFq06aNatWqpcuXL+u9997TyZMnNXDgQGs2BQAAh0SOBAC4OqsGWHr11VeVmZmZrd0wDL366qv5Xs6MGTOUlJSkiIgIVahQwfxaunSpuc+pU6d09uxZ8/tLly5p0KBBql+/vrp06aLk5GRt3bpVwcHB1mwKAAAOiRwJAHB1Vp1ZPXLkSI6Jr169enneF/NnhmHcsU98fLzF+8mTJ2vy5Mn5XgcAAEURORIA4OqsOrPq7++vX375JVv70aNHVbx48XsOCgAAAADg2qwqVnv06KHhw4fr2LFj5rajR4/qpZdeUvfu3W0WHAAAAADANVlVrL777rsqXry46tWrpxo1aqhGjRqqX7++ypQpo/fff9/WMQIAAAAAXIxV96z6+/tr69atiouL008//SQfHx81btxY7du3t3V8AAAAAAAXZFWxKkkmk0mdOnVSp06dbBkPAAAAAADWF6sbNmzQhg0bdP78eWVlZVlMmzdv3j0HBgAAAABwXVYVq+PGjdP48ePVokULVahQQSaTydZxAQAAAABcmFXF6syZM7VgwQL99a9/tXU8AAAAAABYNxpwWlqawsLCbB0LAAAAAACSrCxWBw4cqMWLF9s6FgAAAAAAJFl5GfCNGzc0e/ZsrV+/Xo0bN5aHh4fF9A8++MAmwQEAAAAAXJNVxerevXsVEhIiSdq/f7/FNAZbAgAAAADcK6uK1U2bNtk6DgAAAAAAzKy6ZxUAAAAAgIJk1ZlVSdq5c6c+++wznTp1SmlpaRbTvvjii3sODAAAAADguqw6s7pkyRKFhYXp559/1ooVK5Senq4DBw5o48aN8vf3t3WMAAAAAAAXY1Wx+s4772jy5MlatWqVPD09NXXqVB06dEi9e/dW1apVbR0jAAAAAMDFWFWsHjt2TF27dpUkeXp66tq1azKZTHrxxRc1e/ZsmwYIAAAAAHA9VhWrAQEBunLliiSpUqVK5sfXXL58WdevX7dddAAAAAAAl2RVsdq+fXvFxcVJknr16qVhw4Zp0KBBeuKJJ9ShQ4d8Lyc2NlYtW7ZUyZIlVa5cOUVFRenw4cN3nG/ZsmWqV6+evL291ahRI61Zs8aazQAAwGGRIwEArs6qYnXatGl6/PHHJUn/+Mc/FBMTo8TERPXs2VNz587N93I2b96sIUOGaPv27YqLi1N6ero6deqka9eu5TrP1q1b9cQTT2jAgAHas2ePoqKiFBUVZT67CwCAMyBHAgBcnckwDMPeQdxy4cIFlStXTps3b1b79u1z7NOnTx9du3ZNX3/9tbmtTZs2CgkJ0cyZM++4juTkZPn7+yspKUl+fn42ix25u56WoeA31kmSDo6PlK+n1U9McmrsJ8B5FESuKYwcWVCxOwu+p4sOfleAY8tvrrH6f25mZqZWrFihn3/+WZIUHBysHj16qFgx678MkpKSJEmlS5fOtc+2bdsUExNj0RYZGamVK1davV4UBkM+SpXSrukePnbOLS1DPkpRislNf1y/qpQMd3tH5PACvIvLzc2qC0SAIoccCVjnelqmvUOwKR8Pd5lMJnuHARQKq6qGAwcOqHv37jp37pzq1q0rSZo4caICAwO1atUqNWzY8K6XmZWVpeHDh6tt27Z5zn/u3DmVL1/eoq18+fI6d+5cjv1TU1OVmppqfp+cnHzXseEeGYaWe45TC7f/Su/bOxjH5SOpRY3ySvD2Upcvx9g7nCLBJ7Omtvf7goIVTq+gcqREnoTza/HWenuHYFMtqgVo2eBQCla4BKv+whs4cKAaNGigM2fOaPfu3dq9e7dOnz6txo0b69lnn7UqkCFDhmj//v1asmSJVfPnJjY2Vv7+/uZXlSpVbLp85EP69ZuFKvKUYjIpwdvL3mEUKSnux3TpRu737wHOoqBypESehHPy8XBXi2oB9g6jQOw8eUkp6c51thjIjVVnVhMSErRz504FBPzvSyAgIEBvv/22WrZsedfLGzp0qL7++mtt2bJFlStXzrNvUFCQEhMTLdoSExMVFBSUY/9Ro0ZZXBKVnJxMIraj68MOybc490Dl6Pof0pc3n1/8TddV8vEtY+eAHNfF61f16OpO9g4DKBQFmSMl8iSck8lk0rLBoU5V1F1Py3S6s8TAnVhVrNapU0eJiYlq0KCBRfv58+dVq1atfC/HMAw9//zzWrFiheLj41WjRo07zhMaGqoNGzZo+PDh5ra4uDiFhobm2N/Ly0teXpytchgevpJncXtH4ZgyUsw/lvYtIV/fknYMBoC9FUaOlMiTcF4mk4mBlYAizqr/wbGxsXrhhRc0duxYtWnTRpK0fft2jR8/XhMnTrS43yWv0Z2GDBmixYsX68svv1TJkiXN99T4+/vLx8dHkhQdHa1KlSopNjZWkjRs2DCFh4dr0qRJ6tq1q5YsWaKdO3dq9uzZ1mwKAAAOiRwJAHB1VhWrjzzyiCSpd+/e5pu7bz0Bp1u3bub3JpNJmZm5X34xY8YMSVJERIRF+/z589WvXz9J0qlTpywGTwkLC9PixYv12muvafTo0apdu7ZWrlxp1aBOAAA4KnIkAMDVWVWsbtq0ySYrz88jXuPj47O19erVS7169bJJDAAAOCJyJADA1VlVrIaHh9s6DgAAAAAAzKwqVrds2ZLn9Pbt21sVDAAAAAAAkpXF6p/vn5Fk8WDivO5TBQAAAADgTtzu3CW7S5cuWbzOnz+vtWvXqmXLlvr2229tHSMAAAAAwMVYdWbV398/W9tDDz0kT09PxcTEaNeuXfccGAAAAADAdVl1ZjU35cuX1+HDh225SAAAAACAC7LqzOrevXst3huGobNnz2rChAkKCQmxRVwAAAAAABdmVbEaEhIik8mU7Rlwbdq00bx582wSGAAAAADAdVlVrB4/ftzivZubmwIDA+Xt7W2ToAAAAAAArs2qYrVatWq2jgMAAAAAADOrBlh64YUX9M9//jNb+7Rp0zR8+PB7jQkAAAAA4OKsKlY///xztW3bNlt7WFiYli9ffs9BAQAAAABcm1XF6h9//JHjs1b9/Pz0+++/33NQAAAAAADXZlWxWqtWLa1duzZb+zfffKP77rvvnoMCAAAAALg2qwZYiomJ0dChQ3XhwgU9+OCDkqQNGzZo0qRJmjJlii3jAwAAAAC4IKuK1WeeeUapqal6++239eabb0qSqlevrhkzZig6OtqmAQIAAAAAXI9VxaokPffcc3ruued04cIF+fj4qESJEraMCwAAAADgwqwqVo8fP66MjAzVrl1bgYGB5vYjR47Iw8ND1atXt1V8AAAAAAAXZNUAS/369dPWrVuztf/www/q16/fvcYEAAAAAHBxVhWre/bsyfE5q23atFFCQkK+l7NlyxZ169ZNFStWlMlk0sqVK/PsHx8fL5PJlO117ty5u9wCAAAcH3kSAODKrCpWTSaTrly5kq09KSlJmZmZ+V7OtWvX1KRJE02fPv2u1n/48GGdPXvW/CpXrtxdzQ8AQFFAngQAuDKr7llt3769YmNj9emnn8rd3V2SlJmZqdjYWN1///35Xk7nzp3VuXPnu15/uXLlVKpUqbueDwCAooQ8CQBwZVYVqxMnTlT79u1Vt25dtWvXTpL03XffKTk5WRs3brRpgDkJCQlRamqqGjZsqLFjx+Z4STIAAK6KPAkAcAZWXQYcHBysvXv3qk+fPjp//ryuXLmi6OhoHTp0SA0bNrR1jGYVKlTQzJkz9fnnn+vzzz9XlSpVFBERod27d+c6T2pqqpKTky1eAAA4I/IkAMCZWP2cVV9fX5UuXVoVKlSQJJUoUcJ8SXBBqVu3rurWrWt+HxYWpmPHjmny5MlauHBhjvPExsZq3LhxBRoXAACOgDwJAHAmVp1Z3blzp2rWrKnJkyfr4sWLunjxoiZPnqyaNWvmefS2ILRq1UpHjx7NdfqoUaOUlJRkfp0+fboQowMAwL7IkwCAosqqM6svvviiunfvrn/9618qVuzmIjIyMjRw4EANHz5cW7ZssWmQeUlISDCf3c2Jl5eXvLy8Ci0eAAAcCXkSAFBUWVWs7ty506JQlaRixYppxIgRatGiRb6Xc/XqVYujvcePH1dCQoJKly6tqlWratSoUfr111/173//W5I0ZcoU1ahRQw0aNNCNGzc0Z84cbdy4Ud9++601mwEAgEMjTwIAXJlVxaqfn59OnTqlevXqWbSfPn1aJUuWzPdydu7cqQceeMD8PiYmRpLUt29fLViwQGfPntWpU6fM09PS0vTSSy/p119/la+vrxo3bqz169dbLAMAAGdBngQAuDKritU+ffpowIABev/99xUWFiZJ+v777/XKK6/oiSeeyPdyIiIiZBhGrtMXLFhg8X7EiBEaMWKENSEDAFDkkCcBAK7MqmL1/fffl8lkUnR0tDIyMiRJHh4eeu655zRhwgSbBggAAAAAcD1WFauenp6aOnWqYmNjdezYMUlSzZo15evra9PgAAAAAACuyernrEo3n7XaqFEjW8UCAAAAAIAkK5+zCgAAAABAQaJYBQAAAAA4HIpVAAAAAIDDoVgFAAAAADgcilUAAAAAgMOhWAUAAAAAOByKVQAAAACAw6FYBQAAAAA4HIpVAAAAAIDDoVgFAAAAADgcilUAAAAAgMOhWAUAAAAAOByKVQAAAACAw6FYBQAAAAA4HIpVAAAAAIDDoVgFAAAAADgcuxarW7ZsUbdu3VSxYkWZTCatXLnyjvPEx8erWbNm8vLyUq1atbRgwYICjxMAAHsgTwIAXJldi9Vr166pSZMmmj59er76Hz9+XF27dtUDDzyghIQEDR8+XAMHDtS6desKOFIAAAofeRIA4MqK2XPlnTt3VufOnfPdf+bMmapRo4YmTZokSapfv77+85//aPLkyYqMjCyoMAEAsAvyJADAldm1WL1b27ZtU8eOHS3aIiMjNXz4cPsEBACAAyFPwiqGIaVft3cUuJO0DPnohiTp+tVkydPdzgHZlo+Hu0wmk73DsC0PX8nZtqmQFali9dy5cypfvrxFW/ny5ZWcnKyUlBT5+Phkmyc1NVWpqanm98nJyQUeJwAA9kCexF0zDGlepHT6B3tHgjvwlfSz9/+/+ac9I0G+VWkjPbOWgvUeOP1owLGxsfL39ze/qlSpYu+QAABwGORJF5d+nUIVKCint3PVwj0qUmdWg4KClJiYaNGWmJgoPz+/HI8WS9KoUaMUExNjfp+cnEwiBgA4JfIk7snLRyVPX3tHgTwYhqGU9Ex7h2FT19My1e7dTZKkXa91lK9nkSpPcpZ2XXq/lr2jcApF6tMQGhqqNWvWWLTFxcUpNDQ013m8vLzk5eVV0KEBAGB35EncE09fybO4vaNAHkySfJ3tv2tahlL0/9c3exaXnKFYhc3Y9TLgq1evKiEhQQkJCZJuDrmfkJCgU6dOSbp5tDc6Otrcf/Dgwfrll180YsQIHTp0SB999JE+++wzvfjii/YIHwCAAkWeBAC4MrsWqzt37lTTpk3VtGlTSVJMTIyaNm2qN954Q5J09uxZc0KWpBo1amj16tWKi4tTkyZNNGnSJM2ZM4fh+AEATok8CQBwZXY9zx4RESHDMHKdvmDBghzn2bNnTwFGBQCAYyBPAgBcmdOPBgwAAAAAKHooVgEAAAAADodiFQAAAADgcChWAQAAAAAOh2IVAAAAAOBwKFYBAAAAAA6HYhUAAAAA4HAoVgEAAAAADodiFQAAAADgcChWAQAAAAAOh2IVAAAAAOBwKFYBAAAAAA6HYhUAAAAA4HAoVgEAAAAADodiFQAAAADgcChWAQAAAAAOh2IVAAAAAOBwKFYBAAAAAA6HYhUAAAAA4HAcolidPn26qlevLm9vb7Vu3Vo//vhjrn0XLFggk8lk8fL29i7EaAEAKDzkSACAq7J7sbp06VLFxMRozJgx2r17t5o0aaLIyEidP38+13n8/Px09uxZ8+vkyZOFGDEAAIWDHAkAcGV2L1Y/+OADDRo0SP3791dwcLBmzpwpX19fzZs3L9d5TCaTgoKCzK/y5csXYsQAABQOciQAwJXZtVhNS0vTrl271LFjR3Obm5ubOnbsqG3btuU639WrV1WtWjVVqVJFPXr00IEDBwojXAAACg05EgDg6uxarP7+++/KzMzMdtS3fPnyOnfuXI7z1K1bV/PmzdOXX36pRYsWKSsrS2FhYTpz5kyO/VNTU5WcnGzxAgDA0RVGjpTIkwAAx2X3y4DvVmhoqKKjoxUSEqLw8HB98cUXCgwM1KxZs3LsHxsbK39/f/OrSpUqhRwxAACF425zpESeBAA4LrsWq2XLlpW7u7sSExMt2hMTExUUFJSvZXh4eKhp06Y6evRojtNHjRqlpKQk8+v06dP3HDcAAAWtMHKkRJ4EADguuxarnp6eat68uTZs2GBuy8rK0oYNGxQaGpqvZWRmZmrfvn2qUKFCjtO9vLzk5+dn8QIAwNEVRo6UyJMAAMdVzN4BxMTEqG/fvmrRooVatWqlKVOm6Nq1a+rfv78kKTo6WpUqVVJsbKwkafz48WrTpo1q1aqly5cv67333tPJkyc1cOBAe24GAAA2R44EALgyuxerffr00YULF/TGG2/o3LlzCgkJ0dq1a80DSpw6dUpubv87AXzp0iUNGjRI586dU0BAgJo3b66tW7cqODjYXpsAAECBIEcCAFyZ3YtVSRo6dKiGDh2a47T4+HiL95MnT9bkyZMLISoAAOyPHAkAcFVFbjRgAAAAAIDzo1gFAAAAADgcilUAAAAAgMOhWAUAAAAAOByKVQAAAACAw6FYBQAAAAA4HIpVAAAAAIDDoVgFAAAAADgcilUAAAAAgMOhWAUAAAAAOByKVQAAAACAw6FYBQAAAAA4HIpVAAAAAIDDoVgFAAAAADgcilUAAAAAgMOhWAUAAAAAOByKVQAAAACAw6FYBQAAAAA4HIpVAAAAAIDDcYhidfr06apevbq8vb3VunVr/fjjj3n2X7ZsmerVqydvb281atRIa9asKaRIAQAoXORIAICrsnuxunTpUsXExGjMmDHavXu3mjRposjISJ0/fz7H/lu3btUTTzyhAQMGaM+ePYqKilJUVJT2799fyJEDAFCwyJEAAFdm92L1gw8+0KBBg9S/f38FBwdr5syZ8vX11bx583LsP3XqVD388MN65ZVXVL9+fb355ptq1qyZpk2bVsiRAwBQsMiRAABXVsyeK09LS9OuXbs0atQoc5ubm5s6duyobdu25TjPtm3bFBMTY9EWGRmplStXFmSo2WRlZurSlQuFus6iKuXaVclkuvlzRoqU7mHniBxTSkbK/96kXZeKXbNfMI4u7X/75mLSeYv3gC0ElAyUm7u7XWMoyjlSkoysLKVcv1Lo6y1I19My5aMbN9+kXZOd/4yyjbTr9o4AMLuelmnvEGwjLUO+///j9WvJUlqGXcOxNR/fkjK5Fc45T7t+y/7+++/KzMxU+fLlLdrLly+vQ4cO5TjPuXPncux/7ty5HPunpqYqNTXV/D45Ofkeo77p0pULivjyIZssyyVUr3LzX/ZZ/kxtLBmGvaNwWD4mk/kz9ei3UfYNBk4pvkecypQKsmsMhZEjpYLLkynXr8j3/ao2WZaj8JX0s/f/v3nfnpEAzqnFW+vtHYJN+OiG+bvCd2o9+wZTAK6/fEq+JfwLZV12vwy4oMXGxsrf39/8qlKlir1DAvLU9MYN+VCo5snHMNT0xg17hwE4BfIkJElV2kgevnfuB9iYj4e7WlQLsHcYNpUiL+3IqmPvMJyCXc+sli1bVu7u7kpMTLRoT0xMVFBQzke0g4KC7qr/qFGjLC6JSk5OtkkiDigZqPgecfe8HFfi7VOi0C4ZKLIMQz6GIdP/XzaNnJkkLcjK0qUbXP6LghFQMtDeIRRKjpQKLk/6+JbU9ZdP3fNyHJGPh7vzfU97+Jpv2QEKk8lk0rLBoUpJd5JLgG8xInU93Tkvs/fxLVlo67Jrserp6anmzZtrw4YNioqKkiRlZWVpw4YNGjp0aI7zhIaGasOGDRo+fLi5LS4uTqGhoTn29/LykpeXl61Dl5u7u90vEQNcmZukMt6F92UJFLbCyJFSweVJk5tboV0mBqBoM5lM8vV0gnvA/8yL78B7ZfdPRUxMjPr27asWLVqoVatWmjJliq5du6b+/ftLkqKjo1WpUiXFxsZKkoYNG6bw8HBNmjRJXbt21ZIlS7Rz507Nnj3bnpsBAIDNkSMBAK7M7sVqnz59dOHCBb3xxhs6d+6cQkJCtHbtWvMAEadOnZLbbZeOhoWFafHixXrttdc0evRo1a5dWytXrlTDhg3ttQkAABQIciQAwJWZDMO1RnJJTk6Wv7+/kpKS5OfnZ+9wAABOqCjnmqIcOwCgaMhvrmG0GwAAAACAw6FYBQAAAAA4HIpVAAAAAIDDoVgFAAAAADgcilUAAAAAgMOx+6NrCtutwY+Tk5PtHAkAwFndyjFFccB98iQAoKDlN0+6XLF65coVSVKVKlXsHAkAwNlduXJF/v7+9g7jrpAnAQCF5U550uWes5qVlaXffvtNJUuWlMlkuqdlJScnq0qVKjp9+jTPorsD9lX+sJ/yj32Vf+yr/LPVvjIMQ1euXFHFihXl5la07rghT+aNbSo6nHG7nHGbJOfcLrYpb/nNky53ZtXNzU2VK1e26TL9/Pyc5kNY0NhX+cN+yj/2Vf6xr/LPFvuqqJ1RvYU8mT9sU9HhjNvljNskOed2sU25y0+eLFqHewEAAAAALoFiFQAAAADgcChW74GXl5fGjBkjLy8ve4fi8NhX+cN+yj/2Vf6xr/KPfWVbzrg/2aaiwxm3yxm3SXLO7WKbbMPlBlgCAAAAADg+zqwCAAAAABwOxSoAAAAAwOFQrAIAAAAAHA7FqpWmT5+u6tWry9vbW61bt9aPP/5o75Ac0pYtW9StWzdVrFhRJpNJK1eutHdIDik2NlYtW7ZUyZIlVa5cOUVFRenw4cP2DsshzZgxQ40bNzY/4ys0NFTffPONvcNyeBMmTJDJZNLw4cPtHYrDGTt2rEwmk8WrXr169g6ryHO2POmM+cwZc4+r5Ahn+E535u/eX3/9VU8//bTKlCkjHx8fNWrUSDt37rR3WFarXr16tt+VyWTSkCFDCnzdFKtWWLp0qWJiYjRmzBjt3r1bTZo0UWRkpM6fP2/v0BzOtWvX1KRJE02fPt3eoTi0zZs3a8iQIdq+fbvi4uKUnp6uTp066dq1a/YOzeFUrlxZEyZM0K5du7Rz5049+OCD6tGjhw4cOGDv0BzWjh07NGvWLDVu3NjeoTisBg0a6OzZs+bXf/7zH3uHVKQ5Y550xnzmjLnHFXKEM32nO+N376VLl9S2bVt5eHjom2++0cGDBzVp0iQFBATYOzSr7dixw+L3FBcXJ0nq1atXwa/cwF1r1aqVMWTIEPP7zMxMo2LFikZsbKwdo3J8kowVK1bYO4wi4fz584YkY/PmzfYOpUgICAgw5syZY+8wHNKVK1eM2rVrG3FxcUZ4eLgxbNgwe4fkcMaMGWM0adLE3mE4FWfPk86az5w19zhTjnCm73Rn/e4dOXKkcf/999s7jAI1bNgwo2bNmkZWVlaBr4szq3cpLS1Nu3btUseOHc1tbm5u6tixo7Zt22bHyOBMkpKSJEmlS5e2cySOLTMzU0uWLNG1a9cUGhpq73Ac0pAhQ9S1a1eL7yxkd+TIEVWsWFH33XefnnrqKZ06dcreIRVZ5Mmiy9lyjzPmCGf7TnfG796vvvpKLVq0UK9evVSuXDk1bdpU//rXv+wdls2kpaVp0aJFeuaZZ2QymQp8fcUKfA1O5vfff1dmZqbKly9v0V6+fHkdOnTITlHBmWRlZWn48OFq27atGjZsaO9wHNK+ffsUGhqqGzduqESJElqxYoWCg4PtHZbDWbJkiXbv3q0dO3bYOxSH1rp1ay1YsEB169bV2bNnNW7cOLVr10779+9XyZIl7R1ekUOeLJqcKfc4a45wtu90Z/3u/eWXXzRjxgzFxMRo9OjR2rFjh1544QV5enqqb9++9g7vnq1cuVKXL19Wv379CmV9FKuAgxkyZIj279/vFPdtFJS6desqISFBSUlJWr58ufr27avNmzc7xR8jtnL69GkNGzZMcXFx8vb2tnc4Dq1z587mnxs3bqzWrVurWrVq+uyzzzRgwAA7RgYUHmfKPc6YI5zxO91Zv3uzsrLUokULvfPOO5Kkpk2bav/+/Zo5c6ZTFKtz585V586dVbFixUJZH8XqXSpbtqzc3d2VmJho0Z6YmKigoCA7RQVnMXToUH399dfasmWLKleubO9wHJanp6dq1aolSWrevLl27NihqVOnatasWXaOzHHs2rVL58+fV7NmzcxtmZmZ2rJli6ZNm6bU1FS5u7vbMULHVapUKdWpU0dHjx61dyhFEnmy6HG23OOMOcIVvtOd5bu3QoUK2Q6M1K9fX59//rmdIrKdkydPav369friiy8KbZ3cs3qXPD091bx5c23YsMHclpWVpQ0bNjjN/RAofIZhaOjQoVqxYoU2btyoGjVq2DukIiUrK0upqan2DsOhdOjQQfv27VNCQoL51aJFCz311FNKSEgo8n/UFKSrV6/q2LFjqlChgr1DKZLIk0WHq+QeZ8gRrvCd7izfvW3bts32CKj//ve/qlatmp0isp358+erXLly6tq1a6GtkzOrVoiJiVHfvn3VokULtWrVSlOmTNG1a9fUv39/e4fmcK5evWpxhOz48eNKSEhQ6dKlVbVqVTtG5liGDBmixYsX68svv1TJkiV17tw5SZK/v798fHzsHJ1jGTVqlDp37qyqVavqypUrWrx4seLj47Vu3Tp7h+ZQSpYsme2+s+LFi6tMmTJF/n40W3v55ZfVrVs3VatWTb/99pvGjBkjd3d3PfHEE/YOrchyxjzpjPnMGXOPs+YIZ/xOd9bv3hdffFFhYWF655131Lt3b/3444+aPXu2Zs+ebe/Q7klWVpbmz5+vvn37qlixQiwhC3y8YSf14YcfGlWrVjU8PT2NVq1aGdu3b7d3SA5p06ZNhqRsr759+9o7NIeS0z6SZMyfP9/eoTmcZ555xqhWrZrh6elpBAYGGh06dDC+/fZbe4dVJBT1xxwUlD59+hgVKlQwPD09jUqVKhl9+vQxjh49au+wijxny5POmM+cMfe4Uo4o6t/pzvzdu2rVKqNhw4aGl5eXUa9ePWP27Nn2DumerVu3zpBkHD58uFDXazIMwyi80hgAAAAAgDvjnlUAAAAAgMOhWAUAAAAAOByKVQAAAACAw6FYBQAAAAA4HIpVAAAAAIDDoVgFAAAAADgcilUAAAAAgMOhWAUAAAAAOByKVQD51q9fP0VFRd3TMuLj42UymXT58mWbxAQAgKMgTwK2VczeAQAoOqZOnSrDMOwdBgAADok8CdgWxSqAO8rMzJTJZJK/v7+9QwEAwOGQJ4GCwWXAgBOKiIjQ0KFDNXToUPn7+6ts2bJ6/fXXzUd7U1NT9fLLL6tSpUoqXry4Wrdurfj4ePP8CxYsUKlSpfTVV18pODhYXl5eOnXqVLbLm1JTU/XCCy+oXLly8vb21v33368dO3ZYxLJmzRrVqVNHPj4+euCBB3TixIlC2AMAAOSOPAkUDRSrgJP6+OOPVaxYMf3444+aOnWqPvjgA82ZM0eSNHToUG3btk1LlizR3r171atXLz388MM6cuSIef7r169r4sSJmjNnjg4cOKBy5cplW8eIESP0+eef6+OPP9bu3btVq1YtRUZG6uLFi5Kk06dP69FHH1W3bt2UkJCggQMH6tVXXy2cHQAAQB7Ik0ARYABwOuHh4Ub9+vWNrKwsc9vIkSON+vXrGydPnjTc3d2NX3/91WKeDh06GKNGjTIMwzDmz59vSDISEhIs+vTt29fo0aOHYRiGcfXqVcPDw8P45JNPzNPT0tKMihUrGu+++65hGIYxatQoIzg42GIZI0eONCQZly5dstXmAgBwV8iTQNHAPauAk2rTpo1MJpP5fWhoqCZNmqR9+/YpMzNTderUseifmpqqMmXKmN97enqqcePGuS7/2LFjSk9PV9u2bc1tHh4eatWqlX7++WdJ0s8//6zWrVtbzBcaGnpP2wUAgC2QJwHHR7EKuJirV6/K3d1du3btkru7u8W0EiVKmH/28fGxSOIAALgC8iTgOLhnFXBSP/zwg8X77du3q3bt2mratKkyMzN1/vx51apVy+IVFBSU7+XXrFlTnp6e+v77781t6enp2rFjh4KDgyVJ9evX148//pgtDgAA7I08CTg+ilXASZ06dUoxMTE6fPiwPv30U3344YcaNmyY6tSpo6eeekrR0dH64osvdPz4cf3444+KjY3V6tWr87384sWL67nnntMrr7yitWvX6uDBgxo0aJCuX7+uAQMGSJIGDx6sI0eO6JVXXtHhw4e1ePFiLViwoIC2GACA/CNPAo6Py4ABJxUdHa2UlBS1atVK7u7uGjZsmJ599llJ0vz58/XWW2/ppZde0q+//qqyZcuqTZs2euSRR+5qHRMmTFBWVpb++te/6sqVK2rRooXWrVungIAASVLVqlX1+eef68UXX9SHH36oVq1a6Z133tEzzzxj8+0FAOBukCcBx2cyjP9/oBQApxEREaGQkBBNmTLF3qEAAOBwyJNA0cBlwAAAAAAAh0OxCgAAAABwOFwGDAAAAABwOJxZBQAAAAA4HIpVAAAAAIDDoVgFAAAAADgcilUAAAAAgMOhWAUAAAAAOByKVQAAAACAw6FYBQAAAAA4HIpVAAAAAIDDoVgFAAAAADic/wOO+tAUzlnZoAAAAABJRU5ErkJggg==",
      "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",
    "for a in range(scenario_a.n_links):\n",
    "    axes[0].step(range(trajectory_a.occupancies.shape[0]), trajectory_a.occupancies[:, a],\n",
    "                 where=\"post\", label=f\"link {a}\")\n",
    "axes[0].set_title(f\"anchor A (cost {metrics_a['total_cost']:.0f})\")\n",
    "axes[0].set_xlabel(\"period\")\n",
    "axes[0].set_ylabel(\"occupancy\")\n",
    "axes[0].legend(fontsize=8)\n",
    "\n",
    "for a, label in enumerate([\"A (metered)\", \"B\"]):\n",
    "    axes[1].step(range(trajectory_b.occupancies.shape[0]), trajectory_b.occupancies[:, a],\n",
    "                 where=\"post\", label=f\"link {label}\")\n",
    "axes[1].set_title(f\"anchor B (cost {metrics_b['total_cost']:.0f}, holding used)\")\n",
    "axes[1].set_xlabel(\"period\")\n",
    "axes[1].legend(fontsize=8)\n",
    "fig.tight_layout()\n",
    "display(fig)\n",
    "plt.close(fig)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07d34f1a",
   "metadata": {},
   "source": [
    "## Takeaways & pointers\n",
    "\n",
    "- **Certified, not self-reported.** Both anchors' costs came from\n",
    "  `SODTAEvaluator`, which re-solves the LP itself and cross-checks the dual\n",
    "  certificate by pure arithmetic — never the solver's own objective claim.\n",
    "- **Holding is a real phenomenon**, not an LP artifact to explain away: anchor\n",
    "  B strictly prefers it, at a real cost gap versus the naive exit-equality\n",
    "  policy recomputed above.\n",
    "- **Where next.** the spillback-capable cell LP\n",
    "  [`lp-so-dta`](02-lp-so-dta.ipynb) (finite storage, which exit functions\n",
    "  cannot represent); the lineage in the\n",
    "  [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": "merchant-nemhauser"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}