{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6423a47e",
   "metadata": {},
   "source": [
    "# `od-dynamic` — Cascetta, Inaudi & Marquis's (1993) within-day dynamic OD estimation\n",
    "\n",
    "**What.** `od-dynamic` (T2d, ADR-023) estimates a SEQUENCE of time-sliced OD\n",
    "matrices `d_h` from time-sliced link counts, linked by an EXOGENOUS lagged\n",
    "assignment map (frozen free-flow travel times — no iterative equilibrium).\n",
    "The paper offers two GLS variants: SIMULTANEOUS (all slices jointly,\n",
    "efficient but needs the whole horizon before solving) and SEQUENTIAL (slice\n",
    "by slice, earlier estimates frozen once made — online-capable, but provably\n",
    "less efficient). This notebook covers all three registered T2d units:\n",
    "`od-dynamic-sim`, `od-dynamic-seq`, and the `prior-profile` baseline they\n",
    "must both beat.\n",
    "\n",
    "**Why it is in the benchmark.** It is the benchmark's first WITHIN-DAY\n",
    "estimation task — the estimand is a `(H, Z, Z)` demand PROFILE, not a single\n",
    "`(Z, Z)` matrix, and the exogenous lag tensor is a genuinely new artifact\n",
    "(no assignment equilibrium is solved; certification is exact linear\n",
    "algebra). The simultaneous/sequential efficiency gap is the paper's\n",
    "headline result, made an executable, hand-verified fact below. See the\n",
    "[model compendium](../../docs/MODELS.md) (Cascetta, Inaudi & Marquis 1993)\n",
    "and\n",
    "[docs/design/adr-023-od-dynamic.md](../../docs/design/adr-023-od-dynamic.md)\n",
    "(P1).\n",
    "\n",
    "**Scope.** Anchor A1 (integer lag: both estimators reduce EXACTLY to the\n",
    "static `gls` closed form, per slice); Anchor A2 (fractional lag — \"the\"\n",
    "anchor: simultaneous strictly dominates sequential, exact rationals); and an\n",
    "end-to-end certified run on the two-route corridor via the public\n",
    "`run_dynamic_estimation_experiment` API.\n",
    "\n",
    "**Canon.** `[cascetta1993dynamic]`, [docs/REFERENCES.md](../../docs/REFERENCES.md) / [docs/references.bib](../../docs/references.bib)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa4a7de5",
   "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 — the anchors are recomputed as\n",
    "exact closed forms in-cell (no trusted digits, exact rationals via\n",
    "`fractions.Fraction` where the paper's own algebra is being verified), and\n",
    "the end-to-end recovery is recomputed by the P1 `DynamicODCertifier` (via\n",
    "`run_dynamic_estimation_experiment`) from the emitted OD profile, never from\n",
    "the estimator's self-report\n",
    "([README](../../README.md), *Certified, not self-reported*)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "61a61659",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T13:47:53.311443Z",
     "iopub.status.busy": "2026-07-21T13:47:53.311184Z",
     "iopub.status.idle": "2026-07-21T13:47:55.251676Z",
     "shell.execute_reply": "2026-07-21T13:47:55.250536Z"
    }
   },
   "outputs": [],
   "source": [
    "# Setup. `od-dynamic` is a core estimator family: 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",
    "from fractions import Fraction\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from tabench import (\n",
    "    Budget,\n",
    "    DynamicPriorBaseline,\n",
    "    RngBundle,\n",
    "    SequentialDynamicGLSEstimator,\n",
    "    SimultaneousDynamicGLSEstimator,\n",
    "    run_dynamic_estimation_experiment,\n",
    "    two_route_scenario,\n",
    ")\n",
    "from tabench.estimation import (\n",
    "    dynamic_gls_sequential,\n",
    "    dynamic_gls_simultaneous,\n",
    "    predict_interval_counts,\n",
    "    stacked_tensor_map,\n",
    "    tensor_blocks,\n",
    ")\n",
    "\n",
    "\n",
    "def _tensor(m0: float, m1: float) -> np.ndarray:\n",
    "    # Single-pair single-sensor lag tensor with lag-0/lag-1 fractions (m0, m1).\n",
    "    m = np.zeros((2, 1, 1), dtype=np.float64)\n",
    "    m[0, 0, 0] = m0\n",
    "    m[1, 0, 0] = m1\n",
    "    return m\n",
    "\n",
    "\n",
    "def _sim(m, counts, prior, w_prior, v_count, n_intervals):\n",
    "    n_slices, n_pairs = prior.shape\n",
    "    a = stacked_tensor_map(m, n_slices, n_intervals)\n",
    "    return dynamic_gls_simultaneous(\n",
    "        a, counts.reshape(-1), prior.reshape(-1), w_prior.reshape(-1), v_count.reshape(-1)\n",
    "    ).reshape(n_slices, n_pairs)\n",
    "\n",
    "\n",
    "def _seq(m, counts, prior, w_prior, v_count, n_intervals):\n",
    "    blocks = tensor_blocks(m, prior.shape[0], n_intervals)\n",
    "    return dynamic_gls_sequential(blocks, counts, prior, w_prior, v_count, n_intervals)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ee0383a",
   "metadata": {},
   "source": [
    "## Anchor A1: integer lag reduces EXACTLY to static `gls`\n",
    "\n",
    "With `tau = Delta` exactly (lag fractions `M0=0, M1=p=1`), each slice couples\n",
    "to exactly one interval — the stacked GLS block-diagonalizes and BOTH\n",
    "estimators reduce, per slice, to the static `gls` closed form\n",
    "`d_h = (z_h + p*c_{h+1}) / (1 + p^2)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0b563cb4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T13:47:55.256332Z",
     "iopub.status.busy": "2026-07-21T13:47:55.255942Z",
     "iopub.status.idle": "2026-07-21T13:47:55.265039Z",
     "shell.execute_reply": "2026-07-21T13:47:55.264369Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "predicted interval counts : [0. 4. 6. 5.]\n",
      "simultaneous : [3.5 4.5 4. ]\n",
      "sequential   : [3.5 4.5 4. ]\n",
      "static gls   : [3.5 4.5 4. ]\n"
     ]
    }
   ],
   "source": [
    "p = 1.0\n",
    "m = _tensor(0.0, p)\n",
    "n_slices, n_intervals = 3, 4  # T = H + 1: the last slice is observed at lag 1\n",
    "truth_a1 = np.array([[4.0], [6.0], [5.0]])\n",
    "counts_a1 = predict_interval_counts(m, truth_a1, n_intervals)\n",
    "print(f\"predicted interval counts : {counts_a1.reshape(-1)}\")\n",
    "assert np.allclose(counts_a1.reshape(-1), [0.0, 4.0, 6.0, 5.0])\n",
    "\n",
    "prior_a1 = np.full((n_slices, 1), 3.0)\n",
    "w, v = np.ones((n_slices, 1)), np.ones((n_intervals, 1))\n",
    "sim_a1 = _sim(m, counts_a1, prior_a1, w, v, n_intervals)\n",
    "seq_a1 = _seq(m, counts_a1, prior_a1, w, v, n_intervals)\n",
    "closed_a1 = np.array(\n",
    "    [[(prior_a1[h, 0] + p * counts_a1[h + 1, 0]) / (1.0 + p * p)] for h in range(n_slices)]\n",
    ")\n",
    "print(f\"simultaneous : {sim_a1.reshape(-1)}\")\n",
    "print(f\"sequential   : {seq_a1.reshape(-1)}\")\n",
    "print(f\"static gls   : {closed_a1.reshape(-1)}\")\n",
    "assert np.allclose(sim_a1, closed_a1, atol=1e-12)\n",
    "assert np.allclose(seq_a1, closed_a1, atol=1e-12)\n",
    "assert np.allclose(sim_a1, seq_a1, atol=1e-12)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8d2d4bf",
   "metadata": {},
   "source": [
    "## Anchor A2 (\"the\" anchor): fractional lag, simultaneous strictly dominates\n",
    "\n",
    "`tau = Delta/2` (`M0=M1=1/2`), `H=2` slices, `T=3` intervals, truth `(4,6)`,\n",
    "prior `(3,3)`, `V=W=I`. The simultaneous solve lets later counts `c_2, c_3`\n",
    "REVISE slice 1 through the lag entry; the sequential solve freezes slice 1\n",
    "from `c_1` alone once made. Exact rationals (Cascetta et al.'s own algebra,\n",
    "recomputed via `fractions.Fraction`, not quoted from the paper)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e305abdd",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T13:47:55.269042Z",
     "iopub.status.busy": "2026-07-21T13:47:55.268587Z",
     "iopub.status.idle": "2026-07-21T13:47:55.279327Z",
     "shell.execute_reply": "2026-07-21T13:47:55.278665Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "predicted interval counts : [2. 5. 3.]\n",
      "simultaneous (exact) : [Fraction(128, 35), Fraction(142, 35)]  = [3.657142857142857, 4.057142857142857]\n",
      "sequential (exact)   : [Fraction(16, 5), Fraction(94, 25)]  = [3.2, 3.76]\n",
      "|error| simultaneous : [0.34285714 1.94285714]\n",
      "|error| sequential   : [0.8  2.24]\n"
     ]
    }
   ],
   "source": [
    "m2 = _tensor(0.5, 0.5)\n",
    "n_intervals2 = 3\n",
    "truth_a2 = np.array([[4.0], [6.0]])\n",
    "counts_a2 = predict_interval_counts(m2, truth_a2, n_intervals2)\n",
    "print(f\"predicted interval counts : {counts_a2.reshape(-1)}\")\n",
    "assert np.allclose(counts_a2.reshape(-1), [2.0, 5.0, 3.0])\n",
    "\n",
    "prior_a2 = np.array([[3.0], [3.0]])\n",
    "w2, v2 = np.ones((2, 1)), np.ones((3, 1))\n",
    "sim_a2 = _sim(m2, counts_a2, prior_a2, w2, v2, n_intervals2).reshape(-1)\n",
    "seq_a2 = _seq(m2, counts_a2, prior_a2, w2, v2, n_intervals2).reshape(-1)\n",
    "\n",
    "# Exact simultaneous closed form: (I + A'A) d = z + A'c, solved over the rationals.\n",
    "half = Fraction(1, 2)\n",
    "a_rows = [[half, Fraction(0)], [half, half], [Fraction(0), half]]\n",
    "ztc = [Fraction(3), Fraction(3)]\n",
    "for t in range(3):\n",
    "    for h in range(2):\n",
    "        ztc[h] += a_rows[t][h] * Fraction(int(counts_a2[t, 0]))\n",
    "ata = [[sum(a_rows[t][i] * a_rows[t][j] for t in range(3)) for j in range(2)] for i in range(2)]\n",
    "hmat = [[ata[i][j] + (1 if i == j else 0) for j in range(2)] for i in range(2)]\n",
    "det = hmat[0][0] * hmat[1][1] - hmat[0][1] * hmat[1][0]\n",
    "sim_exact = [\n",
    "    (hmat[1][1] * ztc[0] - hmat[0][1] * ztc[1]) / det,\n",
    "    (hmat[0][0] * ztc[1] - hmat[1][0] * ztc[0]) / det,\n",
    "]\n",
    "print(f\"simultaneous (exact) : {sim_exact}  = {[float(x) for x in sim_exact]}\")\n",
    "assert sim_exact == [Fraction(128, 35), Fraction(142, 35)]\n",
    "assert np.allclose(sim_a2, [float(x) for x in sim_exact], atol=1e-12)\n",
    "\n",
    "# Exact sequential: d_1 from c_1 alone (lag 0), then d_2 from c_2 minus frozen d_1 (lag 1).\n",
    "d1 = (Fraction(3) + half * Fraction(2)) / (1 + half * half)\n",
    "resid2 = Fraction(5) - half * d1\n",
    "d2 = (Fraction(3) + half * resid2) / (1 + half * half)\n",
    "print(f\"sequential (exact)   : {[d1, d2]}  = {[float(d1), float(d2)]}\")\n",
    "assert [d1, d2] == [Fraction(16, 5), Fraction(94, 25)]\n",
    "assert np.allclose(seq_a2, [float(d1), float(d2)], atol=1e-12)\n",
    "\n",
    "# The headline: simultaneous strictly dominates sequential, componentwise vs truth.\n",
    "err_sim = np.abs(sim_a2 - truth_a2.reshape(-1))\n",
    "err_seq = np.abs(seq_a2 - truth_a2.reshape(-1))\n",
    "print(f\"|error| simultaneous : {err_sim}\")\n",
    "print(f\"|error| sequential   : {err_seq}\")\n",
    "assert np.all(err_sim < err_seq)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f87debc8",
   "metadata": {},
   "source": [
    "## End to end: two-route, all three T2d units, certified\n",
    "\n",
    "The public `run_dynamic_estimation_experiment` API runs `prior-profile`,\n",
    "`od-dynamic-sim`, and `od-dynamic-seq` on the SAME within-day task and\n",
    "certifies every checkpoint through `DynamicODCertifier`: `od-dynamic-sim`\n",
    "must strictly beat both `od-dynamic-seq` and `prior-profile` on the\n",
    "descriptive `od_rmse` — the efficiency ordering Anchor A2 established above,\n",
    "now confirmed on a real assignment-derived task, not a hand-built tensor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b6c8e769",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T13:47:55.282945Z",
     "iopub.status.busy": "2026-07-21T13:47:55.282635Z",
     "iopub.status.idle": "2026-07-21T13:47:55.302439Z",
     "shell.execute_reply": "2026-07-21T13:47:55.301935Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "scenario      : tworoute\n",
      "content hash  : 1ac534dfed2698fb…\n",
      "prior-profile    od_feasible=1  od_rmse=0.7693  heldout_count_rmse=1.9246\n",
      "od-dynamic-sim   od_feasible=1  od_rmse=0.2389  heldout_count_rmse=1.7798\n",
      "od-dynamic-seq   od_feasible=1  od_rmse=0.5522  heldout_count_rmse=1.8667\n"
     ]
    }
   ],
   "source": [
    "sc = two_route_scenario(sue_theta=None)\n",
    "print(f\"scenario      : {sc.name}\")\n",
    "print(f\"content hash  : {sc.content_hash()[:16]}…\")\n",
    "card = {\n",
    "    \"sensors\": {\"kind\": \"explicit\", \"links\": [3]},\n",
    "    \"heldout\": {\"kind\": \"explicit\", \"links\": [2]},\n",
    "    \"n_slices\": 3, \"slice_length\": 2.0, \"n_days\": 12, \"noise\": \"poisson\",\n",
    "    \"prior\": {\"kind\": \"stale\", \"cv\": 0.3},\n",
    "}\n",
    "result = run_dynamic_estimation_experiment(\n",
    "    sc,\n",
    "    [DynamicPriorBaseline(), SimultaneousDynamicGLSEstimator(), SequentialDynamicGLSEstimator()],\n",
    "    Budget(sp_calls=2000), seed=1, macroreps=4, estimation=card,\n",
    ")\n",
    "finals = {r[\"estimator\"]: r for r in result.rows if r[\"macrorep\"] == 0}\n",
    "for name in (\"prior-profile\", \"od-dynamic-sim\", \"od-dynamic-seq\"):\n",
    "    row = finals[name]\n",
    "    print(f\"{name:16s} od_feasible={row['od_feasible']:.0f}  od_rmse={row['od_rmse']:.4f}  \"\n",
    "          f\"heldout_count_rmse={row['heldout_count_rmse']:.4f}\")\n",
    "    assert row[\"od_feasible\"] == 1.0\n",
    "    assert np.isfinite(float(row[\"heldout_count_rmse\"]))\n",
    "assert result.manifest[\"identifiability\"][\"linear_identifiable\"] is True\n",
    "assert finals[\"od-dynamic-sim\"][\"od_rmse\"] < finals[\"od-dynamic-seq\"][\"od_rmse\"]\n",
    "assert finals[\"od-dynamic-sim\"][\"od_rmse\"] < finals[\"prior-profile\"][\"od_rmse\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f5720e5",
   "metadata": {},
   "source": [
    "## Visualize\n",
    "\n",
    "`tabench.viz` draws a road `Network`'s per-LINK flows; the `od-dynamic`\n",
    "artifact is a per-SLICE demand PROFILE for one OD pair, a different shape\n",
    "entirely, so this notebook plots the certified per-slice profile directly\n",
    "(a house profile plot, not `tabench.viz`) — Anchor A2's exact simultaneous\n",
    "and sequential estimates against the truth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "2067f3ca",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T13:47:55.305899Z",
     "iopub.status.busy": "2026-07-21T13:47:55.305416Z",
     "iopub.status.idle": "2026-07-21T13:47:55.416674Z",
     "shell.execute_reply": "2026-07-21T13:47:55.415826Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhwAAAFJCAYAAADQVfSlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABDSElEQVR4nO3de1zO9/8/8Mclnc9RKqUSIkXOQggTwjDnbWJz+jjEms8wG2WmmfMOZPZRzWFsjhszx8KcSVYjp5U5rtBRFPX6/eHX++tyXem60rsr9bjfbt1urtf79X6/n+/r/b4uj+t9VAghBIiIiIhkVE3XBRAREVHlx8BBREREsmPgICIiItkxcBAREZHsGDiIiIhIdgwcREREJDsGDiIiIpIdAwcRERHJjoGDiIiIZMfAoSOhoaFQKBS4d++erkup8FxdXTFy5EidzLtoPZHude7cGZ07dy6z6ZXXutXl9qvLeVP5GjlyJFxdXUs1bnltJwwc/9+KFSugUCjQpk0bXZeiE4MHD4ZCocD06dPVDk9KSsJHH30EHx8fmJubw8HBAYGBgThz5kw5V6p78+fPx/bt23Vdxmvr9u3bCA0NRXx8vOzzys3NRWhoKGJjY2WfV1W1YcMGLFu2TNdlVAnl+dmRAwPH/7d+/Xq4urri1KlTuHr1qq7LKVdZWVn49ddf4erqih9//BHqHq/z/fffY/Xq1WjZsiUWL16MkJAQXLp0CW3btsX+/ftlre/SpUtYvXq1rPPQBgPHq7l9+zbCwsK0/tLcu3cv9u7dq9U4ubm5CAsLq7KBozw+Owwc5edln53Vq1fj0qVL5V+UFhg4ACQnJ+PYsWNYsmQJbG1tsX79el2XVCaEEHj06FGJ/bZs2YKCggKsWbMGN27cwOHDh1X6DBs2DDdu3MD333+PsWPH4r///S9OnjwJGxsbhIaGylD9/zE0NIS+vr6s86CKKzc3FwBgYGAAAwMDHVfzeuFnp+rQ19eHoaGhrst4KQYOPNu7YW1tjcDAQAwcOFBt4EhJSYFCocCiRYvw3Xffwd3dHYaGhmjVqhVOnz6t0j8pKQmDBw+Gra0tjI2N4eHhgVmzZqn0y8jIwMiRI2FlZQVLS0uMGjVK+oIt8vTpU3z22WfSPF1dXfHxxx8jLy9PqZ+rqyt69+6NPXv2oGXLljA2NsaqVas0Wv433ngD/v7+aNSokdrlb9GiBczMzJTaatSoAT8/P1y8eFGpPTc3F0lJSRqdn3LlyhW89dZbsLe3h5GREZycnDB06FBkZmYqLdfzxxejoqKgUCjwxx9/IDg4GLa2trCyssK4ceOQn5+PjIwMjBgxAtbW1rC2tsZHH32ktNcmNjYWCoVC5Vdv0TqOiooqtl6FQoGHDx8iOjoaCoUCCoVCqu369euYMGECPDw8YGxsjBo1amDQoEFISUlRmkZR/UePHkVISAhsbW1hamqK/v37Iy0tTWWeu3fvhp+fH0xNTWFubo7AwED89ddfKv0OHjwo9bOyssKbb76psm6KO86r7nyGffv2oUOHDrCysoKZmRk8PDzw8ccfF/veaDJebGwsWrVqBQAYNWqU9B4WveedO3eGl5cXzp49i44dO8LExEQaV905HI8fP0ZoaCgaNGgAIyMjODg4YMCAAbh27RpSUlJga2sLAAgLC5PmVVxA7tSpE5o2bap2mIeHBwICAl663EIIzJs3D05OTjAxMYG/v7/a9QQAf//9NwYNGgQbGxuYmJigbdu22LVrl1Kfou30p59+QlhYGGrXrg1zc3MMHDgQmZmZyMvLw9SpU2FnZwczMzOMGjVK7XeCus+OJtvejh07EBgYCEdHRxgaGsLd3R2fffYZCgoKpD6dO3fGrl27cP36den9fX77ysvLw5w5c1CvXj0YGhrC2dkZH330kUqdcmxr2taQl5eHDz74ALa2tjA3N0ffvn1x8+ZNlW1Gm88QAKxbtw4tWrSAsbExbGxsMHToUNy4cUOpT9F2f+HCBfj7+8PExAS1a9fGl19+KfUp6bOjrq5FixahXbt2qFGjBoyNjdGiRQts3ry5pLdVNtV1NucKZP369RgwYAAMDAwwbNgwrFy5EqdPn5ZW7vM2bNiA7OxsjBs3DgqFAl9++SUGDBiAv//+W/ol8eeff8LPzw/6+voYO3YsXF1dce3aNfz666/4/PPPlaY3ePBguLm5ITw8HHFxcfj+++9hZ2eHBQsWSH1Gjx6N6OhoDBw4EB9++CFOnjyJ8PBwXLx4Edu2bVOa3qVLlzBs2DCMGzcOY8aMgYeHx0uX/fbt24iJiUF0dDSAZ3syli5dim+++UajX5N3795FzZo1ldpOnToFf39/zJkz56V7P/Lz8xEQEIC8vDxMnjwZ9vb2uHXrFnbu3ImMjAxYWlq+dN5F44SFheHEiRP47rvvYGVlhWPHjqFOnTqYP38+fvvtNyxcuBBeXl4YMWJEictTkrVr12L06NFo3bo1xo4dCwBwd3cHAJw+fRrHjh3D0KFD4eTkhJSUFKxcuRKdO3fGhQsXYGJiolK/tbU15syZg5SUFCxbtgyTJk3Cpk2blOYXFBSEgIAALFiwALm5uVi5ciU6dOiAc+fOSV8w+/fvR8+ePVG3bl2Ehobi0aNH+Prrr9G+fXvExcVpfTLZX3/9hd69e6NJkyaYO3cuDA0NcfXqVRw9evSVxmvUqBHmzp2L2bNnY+zYsfDz8wMAtGvXTprG/fv30bNnTwwdOhTvvPMOatWqpXZeBQUF6N27Nw4cOIChQ4diypQpyM7Oxr59+5CYmIhu3bph5cqV+M9//oP+/ftjwIABAIAmTZqond67776LMWPGIDExEV5eXlL76dOncfnyZXzyyScvXfbZs2dj3rx56NWrF3r16oW4uDh0794d+fn5Sv3+/fdftGvXDrm5uQgODkaNGjUQHR2Nvn37YvPmzejfv79S//DwcBgbG2PGjBm4evUqvv76a+jr66NatWpIT09HaGgoTpw4gaioKLi5uWH27NkvrRPQbNuLioqCmZkZQkJCYGZmhoMHD2L27NnIysrCwoULAQCzZs1CZmYmbt68iaVLlwKA9MOksLAQffv2xR9//IGxY8eiUaNGSEhIwNKlS3H58mXpsKRc25o2NQDPvmfXrVuH4cOHo127djh48CACAwNLfC9f5vPPP8enn36KwYMHY/To0UhLS8PXX3+Njh074ty5c7CyspL6pqeno0ePHhgwYAAGDx6MzZs3Y/r06fD29kbPnj01+uy8aPny5ejbty/efvtt5OfnY+PGjRg0aBB27tz5ystWKqKKO3PmjAAg9u3bJ4QQorCwUDg5OYkpU6Yo9UtOThYARI0aNcSDBw+k9h07dggA4tdff5XaOnbsKMzNzcX169eVplFYWCj9e86cOQKAeO+995T69O/fX9SoUUN6HR8fLwCI0aNHK/WbNm2aACAOHjwotbm4uAgA4vfff9d4+RctWiSMjY1FVlaWEEKIy5cvCwBi27ZtJY57+PBhoVAoxKeffqrUHhMTIwCIOXPmvHT8c+fOCQDi559/fmk/FxcXERQUJL2OjIwUAERAQIDSe+rr6ysUCoUYP3681Pb06VPh5OQkOnXqpFJfTEyM0nyK1nFkZKTUVrSenmdqaqpUT5Hc3FyVtuPHjwsA4ocfflCpv1u3bkr1f/DBB0JPT09kZGQIIYTIzs4WVlZWYsyYMUrTvHv3rrC0tFRq9/HxEXZ2duL+/ftS2/nz50W1atXEiBEjpLagoCDh4uKiUueLy7l06VIBQKSlpan0fRlNxjt9+rTK+1ykU6dOAoCIiIhQO+z59bhmzRoBQCxZskSlb9H7mpaWVuy2+OIyZ2RkCCMjIzF9+nSlfsHBwcLU1FTk5OQUu0ypqanCwMBABAYGKq3Tjz/+WABQ2l6mTp0qAIgjR45IbdnZ2cLNzU24urqKgoICIcT/badeXl4iPz9f6jts2DChUChEz549lWrw9fVVWbfFfXZK2vaEUL89jxs3TpiYmIjHjx9LbYGBgWq3qbVr14pq1aopLacQQkRERAgA4ujRo0IIebc1TWso+p6dMGGCUr/hw4erbD+afoZSUlKEnp6e+Pzzz5X6JSQkiOrVqyu1F233z39P5OXlCXt7e/HWW29JbS/77Kir68V1mJ+fL7y8vESXLl2U2l/cTuRS5Q+prF+/HrVq1YK/vz+AZ7vMhwwZgo0bNyrtOiwyZMgQWFtbS6+LUubff/8NAEhLS8Phw4fx3nvvoU6dOkrjqtvdNn78eKXXfn5+uH//PrKysgAAv/32GwAgJCREqd+HH34IACq7Yd3c3Erc9fu89evXIzAwEObm5gCA+vXro0WLFiWex5Kamorhw4fDzc0NH330kdKwzp07QwhR4rkdRXsw9uzZo3IYSRPvv/++0nvapk0bCCHw/vvvS216enpo2bKltH7kZGxsLP37yZMnuH//PurVqwcrKyvExcWp9B87dqxS/X5+figoKMD169cBPNtdnJGRgWHDhuHevXvSn56eHtq0aYOYmBgAwJ07dxAfH4+RI0fCxsZGml6TJk3wxhtvSNuQNop+ee3YsQOFhYWyj/c8Q0NDjBo1qsR+W7ZsQc2aNTF58mSVYaW53NXS0hJvvvmm0onTBQUF2LRpE/r16wdTU9Nix92/fz/y8/MxefJkpXlPnTpVpe9vv/2G1q1bo0OHDlKbmZkZxo4di5SUFFy4cEGp/4gRI5TOwyjazt977z2lfm3atMGNGzfw9OnTEpe1pG0PUN6es7Ozce/ePfj5+UmHTEvy888/o1GjRmjYsKHS9tulSxcAkLZfObc1TWso+owEBwcrja9u/Wlq69atKCwsxODBg5XmbW9vj/r160vzLmJmZoZ33nlHem1gYIDWrVu/0nfX8+swPT0dmZmZ8PPzU/t9VB6qdOAoKCjAxo0b4e/vj+TkZFy9ehVXr15FmzZt8O+//+LAgQMq47wYIorCR3p6OoD/Cx7P75J9mZKmd/36dVSrVg316tVT6mdvbw8rKyulLwjgWeDQ1MWLF3Hu3Dm0b99eWvarV6+ic+fO2LlzpxR6XvTw4UP07t0b2dnZ2LFjh8q5HZpyc3NDSEgIvv/+e9SsWRMBAQH49ttvlc7feJkX37uiAOPs7KzSXvR+yunRo0eYPXs2nJ2dYWhoiJo1a8LW1hYZGRlql6mkdX/lyhUAQJcuXWBra6v0t3fvXqSmpgKAtA2oO3zWqFEj3Lt3Dw8fPtRqWYYMGYL27dtj9OjRqFWrFoYOHYqffvqpxP8QSjve82rXrq3R4bxr167Bw8MD1auX3ZHhESNG4J9//sGRI0cAPAsS//77L959992Xjle0DurXr6/Ubmtrq/QDpahvcevq+WkV0WY7Lyws1OjzU9K2Bzw7ZNG/f39YWlrCwsICtra20n+ImszjypUr+Ouvv1S23QYNGgCAtP3Kua1pWkPR92zR4dEiJR2SLmn5hRCoX7++yvwvXrwozbuIk5OTSlC2trZ+pe+unTt3om3btjAyMoKNjQ1sbW2xcuVKjb9jy1qVPofj4MGDuHPnDjZu3IiNGzeqDF+/fj26d++u1Kanp6d2WkLNpaSa0HR6mv5iez7RlmTdunUAgA8++AAffPCByvAtW7ao/NLMz8/HgAED8Oeff2LPnj0aB6viLF68GCNHjsSOHTuwd+9eBAcHIzw8HCdOnICTk9NLxy3uvVPX/vz7Wdx7qW6PljYmT56MyMhITJ06Fb6+vrC0tIRCocDQoUPVfnmWtO6Lxlm7di3s7e1V+pXmP1pNl93Y2BiHDx9GTEwMdu3ahd9//x2bNm1Cly5dsHfv3mJrL+14L05DVwICAlCrVi2sW7cOHTt2xLp162Bvb49u3brprCZttnNAs++iksbNyMhAp06dYGFhgblz58Ld3R1GRkaIi4vD9OnTNQqQhYWF8Pb2xpIlS9QOLwpMcm5rmtagDU0/Q4WFhVAoFNi9e7faZXjxh1pZ/99y5MgR9O3bFx07dsSKFSvg4OAAfX19REZGYsOGDaWa5quq0oFj/fr1sLOzw7fffqsybOvWrdi2bRsiIiK0+gKsW7cuACAxMbFManRxcUFhYSGuXLki/QICnp14lpGRARcXl1JNVwiBDRs2wN/fHxMmTFAZ/tlnn2H9+vVKgaOwsBAjRozAgQMH8NNPP6FTp06lmveLvL294e3tjU8++QTHjh1D+/btERERgXnz5pXJ9F9U9GsuIyNDqf3FX5bFKe4LZ/PmzQgKCsLixYultsePH6vMR1NFv7bs7Oxe+h9e0Tag7hr8pKQk1KxZUzocYG1trbYedcterVo1dO3aFV27dsWSJUswf/58zJo1CzExMS+tp6Txyurunu7u7jh58iSePHlS7KWf2s5LT08Pw4cPR1RUFBYsWIDt27djzJgxJQalonVw5coV6TsAeHaI9cVfqC4uLsWuq+enpUuxsbG4f/8+tm7dio4dO0rtycnJKn2Le4/d3d1x/vx5dO3atcT1INe2pmkNRd+zRXvNiqhbT5p+htzd3SGEgJubm7RH5VVpsz1v2bIFRkZG2LNnj9LlspGRkWVSS2lU2UMqjx49wtatW9G7d28MHDhQ5W/SpEnIzs7GL7/8otV0bW1t0bFjR6xZswb//POP0rDSJNVevXoBgMqNdYoSe2nPND569ChSUlIwatQotcs/ZMgQxMTE4Pbt29I4kydPxqZNm7BixQrpjH91NL0sNisrS+V4s7e3N6pVq6ZyyVpZcnFxgZ6ensr9RlasWKHR+Kampmq/cPT09FTW8ddff13qPScBAQGwsLDA/Pnz8eTJE5XhRZcxOjg4wMfHB9HR0Up1JSYmYu/evdI2BDz7EszMzMSff/4ptd25c0flaqcHDx6ozM/HxwcAXrpuNBmvKPyUNogVeeutt3Dv3j188803KsOK1kPRlUHazOvdd99Feno6xo0bh5ycHKXj6sXp1q0b9PX18fXXXyttA+puiNWrVy+cOnUKx48fl9oePnyI7777Dq6urvD09NS4VrkUBaznlyU/P1/tZ8TU1FTtLvrBgwfj1q1bam889ujRI+kwn5zbmqY19OzZEwDw1VdfKfVRt/40/QwNGDAAenp6CAsLU/leEELg/v37xS5bcbT57Ojp6UGhUCh9/6SkpOj0poVVdg/HL7/8guzsbPTt21ft8LZt20o3ARsyZIhW0/7qq6/QoUMHNG/eHGPHjoWbmxtSUlKwa9cure+u2LRpUwQFBeG7776TdnOeOnUK0dHR6Nevn3Syq7bWr18PPT29YgNL3759MWvWLGzcuBEhISFYtmwZVqxYAV9fX5iYmEiHY4r0799f+jBoelnswYMHMWnSJAwaNAgNGjTA06dPsXbtWujp6eGtt94q1XJpwtLSEoMGDcLXX38NhUIBd3d37Ny5U+WYanFatGiB/fv3Y8mSJXB0dISbmxvatGmD3r17Y+3atbC0tISnpyeOHz+O/fv3o0aNGqWq08LCAitXrsS7776L5s2bY+jQobC1tcU///yDXbt2oX379tJ/tgsXLkTPnj3h6+uL999/X7os1tLSUmkdDB06FNOnT0f//v0RHBwsXWbboEEDpRPJ5s6di8OHDyMwMBAuLi5ITU3FihUr4OTkpHSy44s0Gc/d3R1WVlaIiIiAubk5TE1N0aZNG63OPwKenW/xww8/ICQkBKdOnYKfnx8ePnyI/fv3Y8KECXjzzTdhbGwMT09PbNq0CQ0aNICNjQ28vLxeeiiwWbNm8PLykk44bN68eYm12NraYtq0aQgPD0fv3r3Rq1cvnDt3Drt371a5bHzGjBn48ccf0bNnTwQHB8PGxgbR0dFITk7Gli1bUK2a7n8HtmvXDtbW1ggKCkJwcDAUCgXWrl2r9kdTixYtsGnTJoSEhKBVq1YwMzNDnz598O677+Knn37C+PHjERMTg/bt26OgoABJSUn46aefpPsFybmtaVqDj48Phg0bhhUrViAzMxPt2rXDgQMH1N51WtPPkLu7O+bNm4eZM2ciJSUF/fr1g7m5OZKTk7Ft2zaMHTsW06ZN02q9aPPZCQwMxJIlS9CjRw8MHz4cqamp+Pbbb1GvXj2lsFSuZL8OpoLq06ePMDIyEg8fPiy2z8iRI4W+vr64d++edMnkwoULVfpBzWV3iYmJon///sLKykoYGRkJDw8PpctHiy6hevGSrqLL1pKTk6W2J0+eiLCwMOHm5ib09fWFs7OzmDlzptKlaUI8u7QpMDCwxGXPz88XNWrUEH5+fi/t5+bmJpo1ayaEeHbJFYBi/56vV9PLYv/++2/x3nvvCXd3d2FkZCRsbGyEv7+/2L9/v8pyqbu07/Tp00r9intPg4KChKmpqVJbWlqaeOutt4SJiYmwtrYW48aNE4mJiRpdFpuUlCQ6duwojI2NlS55TE9PF6NGjRI1a9YUZmZmIiAgQCQlJWlcf3GX68bExIiAgABhaWkpjIyMhLu7uxg5cqQ4c+aMUr/9+/eL9u3bC2NjY2FhYSH69OkjLly4IF60d+9e4eXlJQwMDISHh4dYt26dynIeOHBAvPnmm8LR0VEYGBgIR0dHMWzYMHH58mWV6T1P0/F27NghPD09RfXq1ZXe806dOonGjRurnfaLl8UK8eyyv1mzZkmfDXt7ezFw4EBx7do1qc+xY8dEixYthIGBgdJ2qW7dFvnyyy8FADF//vyXLu/zCgoKRFhYmHBwcBDGxsaic+fOIjExUe0lh9euXRMDBw6Uvh9at24tdu7cqdSnaHt48bJxbbb/V9n2jh49Ktq2bSuMjY2Fo6Oj+Oijj8SePXtU+uXk5Ijhw4cLKysrAUDp0sz8/HyxYMEC0bhxY2FoaCisra1FixYtRFhYmMjMzBRCyL+taVKDEEI8evRIBAcHixo1aghTU1PRp08fcePGDbXfZZp8hops2bJFdOjQQZiamgpTU1PRsGFDMXHiRHHp0iWpT3HbvbpLXYv77Kjr+7///U/Ur19fGBoaioYNG4rIyEi1dZbXZbEKIUp5RgoRUSW1fPlyfPDBB0hJSVG5ooOqFoVCUeLeWtKM7vfdERFVIEII/O9//0OnTp0YNojKUJU9h4OI6HkPHz7EL7/8gpiYGCQkJGDHjh26LomoUmHgICLCs6t+hg8fDisrK3z88cfFnlBORKXDcziIiIhIdjyHg4iIiGT3Wh9SKSwsxO3bt2Fubl5mdy8kIiIizQghkJ2dDUdHxxLvIfNaB47bt2+X6l74REREVHZu3LhR4vOvXuvAUfRI9Rs3bsDCwkLH1RAREVUtWVlZcHZ2lv4/fpnXOnAUHUaxsLBg4CAiItIRTU5r4EmjREREJDsGDiIiIpLda31IpSQFBQVqH+tNVJHo6+tLjwMnIqqsKm3gyMnJwc2bN9U+TpmoIlEoFHBycoKZmZmuSyEikk2lDBwFBQW4efMmTExMYGtry3t0UIUlhEBaWhpu3ryJ+vXrc08HEVVaOg8ct27dwvTp07F7927k5uaiXr16iIyMRMuWLUs9zSdPnkAIAVtbWxgbG5dhtURlz9bWFikpKXjy5AkDBxFVWjoNHOnp6Wjfvj38/f2xe/du2Nra4sqVK7C2ti6T6XPPBr0OuJ0SUVWg08CxYMECODs7IzIyUmpzc3PTYUVEREQkB50Gjl9++QUBAQEYNGgQDh06hNq1a2PChAkYM2aM2v55eXnIy8uTXmdlZWk8L9cZu165XnVSvgjUuG9oaChmzJgBIyMjreezbNkyDB06FPb29tK0MjIysGzZMq2nRUREVN50Gjj+/vtvrFy5EiEhIfj4449x+vRpBAcHw8DAAEFBQSr9w8PDERYWpoNKy0ZYWBimTp2qEjiePn2K6tVfviqWLVuGzp07S4GDiF6dXD9EXgfa/FgiKgs6DRyFhYVo2bIl5s+fDwBo1qwZEhMTERERoTZwzJw5EyEhIdLronu4vw7Gjx8PAPDz84Oenh4cHR1hb2+Pq1evIjU1FUlJSVAoFEhPT4eVlRUAoGbNmjhz5gx++OEH3L59G0OGDIGxsTGioqIAAHfu3EGfPn1w7do12NvbY/PmzbCxsdHREhIRERVPp3cadXBwgKenp1Jbo0aN8M8//6jtb2hoKD035XV7fkpERAQA4MiRI4iPj4ednR3Onj2LXbt2ISkp6aXjzp49G46Ojti0aRPi4+Ph4+MDADh58iSioqJw4cIF2NnZYdWqVXIvBhERUanoNHC0b98ely5dUmq7fPkyXFxcdFRR+Ro0aJBGT9grTo8ePVCjRg0AgK+vL65du1ZWpREREZUpnQaODz74ACdOnMD8+fNx9epVbNiwAd999x0mTpyoy7LKzYt3ltTT00NBQYH0+vHjxy8d//lzQfT09PD06dOyLZCIiKiM6DRwtGrVCtu2bcOPP/4ILy8vfPbZZ1i2bBnefvttXZYlG3Nzc2RmZhY7vF69ejh58iQAYOvWrXj48KE0zMLC4qXjEhERVWQ6v9No79690bt3b9nnUxHOyP7www/xxhtvwMTEBI6OjirDly5diuDgYHzyyScIDAyUDpcAQHBwMMaMGQMTExPppFEiIqLXhUK8xk83y8rKgqWlJTIzM5VOIH38+DGSk5Ph5uZWqnteEJUnbq+6w8tiiV5Ncf8Pq6PTQypERERUNTBwEBERkewYOIiIiEh2DBxEREQkOwYOIiIikh0DBxEREclO5/fhKDehljJNlzfjIiIiKgn3cFQwiYmJcHV11ahvTk4OFAqFvAVpKCIiAgsXLqxw0xo0aBCOHz9eJtPSVHx8PDZu3KjU5ufnh+Tk5HKtg4ioIqk6ezhIVuPHj69w0zp16hQePHgAX1/fMpmepuLj47F9+3YMHTpUavvwww8xZ84c/PDDD+VaCxFRRcE9HOVkz549aN68OZo0aYJOnTrhwoUL0rDQ0FDUr18fLVq0UPll/KJVq1ahfv36aNasGZYuXSq1L1q0CGPHjpVeZ2RkoGbNmnjw4AGioqLQrVs3DBs2DN7e3mjZsiX+/vtvAMDdu3fh7++PFi1aoHHjxpg0aRIKCwsBQGk8T09PtGvXDhcuXED//v3RqFEjdO/eHTk5OdIyTJ06VZr/ggUL4O3tjaZNm6Jt27bIzc1VWZYrV66gffv2aNq0Kby9vfHJJ5+oTEubGtS9V8OHD5deZ2dnY8yYMWjdujWaNGmCsWPHIj8/H5cuXYKTk5P0nixatAg9evRAYWEhEhIS0KFDBzRv3hyenp6YN2+eNL38/Hz897//hZeXF5o2bYoePXogNTUVs2fPRkxMDHx8fKTwFBgYiN27d/N5OERUZTFwlIPU1FQMHz4c0dHR+PPPPzF27FgMHDgQQgjs2rULP//8M86ePYszZ84gJSWl2OkkJiZizpw5OHz4MM6dO4dHjx5Jw0aPHo3t27cjIyMDABAZGYk333wTNjY2AIDTp09j/vz5SEhIQLdu3bBgwQIAgJWVFX799VecPXsWf/75J1JSUvDTTz9J0z19+jQWLFiACxcuwN3dHX369EFERAQuXrwIAwMDREdHq9QZHR2NLVu24I8//sD58+exe/duGBoaqvT75ptv0Lt3b5w/fx4JCQkICQlRu9ylqQEAYmNj0aZNG+n1hx9+CD8/P5w6dQrnz59HYWEhli9fDg8PDyxcuBCDBw9GbGwsvv32W6xduxbVqlWDq6srDhw4gLi4OJw9exZbtmzBiRMnAADh4eG4fPkyzp49i/Pnz2Pt2rWws7PD3Llz4e/vj/j4eERERAAA9PX14e3tjSNHjhS7fomIKjMGjnJw8uRJeHt7w9vbGwDw9ttv4/bt27h16xYOHDiAwYMHw8LCAgqFAuPGjSt2OgcPHkTPnj3h4OAAAPjPf/4jDbOyssLAgQOxZs0aCCGwcuVKTJo0SRru6+sLNzc36d/Xrl0DABQWFmL69Olo2rQpmjVrhjNnziA+Pl5pvDp16gAAWrZsiVatWqFWrVoAnj3t98qVKyp17ty5E+PHj4el5bMTda2traGnp6fSr2PHjli9ejVmzZqFvXv3wsrKSu1yl6YGALh586bUDwC2b9+OhQsXwsfHB82aNcORI0dw9epVAMCwYcPQvHlzBAQEYO3atbC1tQUAPHr0CKNHj4a3tzfatm2L69evS+/Pzp07MWXKFClMFY1THHt7e9y8efOlfYiIKiuew1HBPH8S6BdffCEdYinaI1FcX+DZE2X79u2LRo0awdbWFs2aNZOGPf9QMD09PTx9+hQAsGTJEqSmpuLkyZMwMjJCSEgIHj9+XOx4xU1HE8HBwTh8+DAAYO3atXjrrbfQrl077Nu3D9988w2WLVuG3377TWW80tZgYmKitCxCCGzZsgUNGjRQ6fv06VMkJibCxsYGt27dkto//vhj1KxZE+fOnUP16tUxYMAApWlq4/HjxzA2Ni7VuERErzvu4SgHbdu2RUJCAhITEwEAGzduRO3atVG7dm1069YNP//8M7KzsyGEwHfffSeNN2PGDMTHxyM+Ph4BAQHo0qULfv/9d9y9excApN31RRo2bIi6deti7NixSns3XiY9PR329vYwMjLC3bt38fPPP7/y8vbt2xcRERHS+QoZGRkoKCjAV199JS2Pt7c3rly5glq1amHEiBH48ssvpUMVZaVJkya4dOmS9Lpfv35YsGCBFFDS09OlPRwzZsyAh4cHjhw5gmnTpknt6enpcHJyQvXq1XHp0iXs27dPaTmXL1+OvLw8AEBaWhoAwMLCQu25GhcvXkTTpk3LdBmJiF4XVWcPhw7vl2Fra4v169djxIgRePr0KaytrfHzzz9DoVCgV69eOHXqFJo3bw4LCwv07Nmz2Ol4eXkhNDQUfn5+MDMzw4ABA1T6jBkzBpMmTcLAgQM1qm3KlCkYOHAgGjduDEdHR3Tr1q3Uy1nk3Xffxe3bt9GuXTtUr14dpqam2L9/P0xMTJT6bd68GevWrYOBgQEKCwtVAtSrGjhwIPbs2SMt09KlSzFjxgz4+PigWrVqqF69Or788kskJSXh999/x6lTp2BiYoIlS5Zg8ODBOHbsGD755BO8++67iI6Ohru7O7p06SJNf/r06Zg1axaaN28OfX19ODo64rfffkPXrl2xaNEiNGnSBO3atUNERARSUlJQUFDAwEFEVZZCCCF0XURpZWVlwdLSEpmZmbCwsJDaHz9+jOTkZLi5uSntfq8KJk2ahFq1auHTTz/VdSk6l5OTg3bt2uH48eMwNTXVaS0zZsxAvXr1MHr0aJVhVXl71TXXGbt0XYLOpHwRqOsSqBIo7v9hdarOHo5K7vbt2+jSpQtsbGywZ88eXZdTIZiZmWHp0qVITk6Gl5eXTmtxdHTEe++9p9MaiIh0iYGjknB0dERSUpKuy6hwunbtqusSADw7YZaIqCrjSaNEREQkOwYOIiIikh0DBxEREcmOgYOIiIhkV2VOGvWO9pZluglBCWU6vcTERPTu3fulz1QpkpOTA3Nzc1SEK5sjIiKQnZ2N//73v7ouRcmgQYMQEhJSrk+MjY+PR1JSktLTYv38/PDDDz9It5cnIqpquIeDysT48eMrXNjQ5ePpX3zqb9Hj6YmIqioGjnLCx9Nr/nj6J0+eYMaMGWjdujV8fHwwePBgpKenAwDu3LmDgIAAeHp6olu3bhg6dChCQ0OLfa/4eHoiooqBgaMc8PH02j2efuHChTA1NcWpU6ek564UhZHg4GC0bt0aFy5cQHR0NA4cOFDs+8XH0xMRVRwMHOWAj6fX7vH027dvx7p16+Dj4wMfHx/8+OOPSE5OBgAcOHBAuj147dq10bdv32LfLz6enoio4tDpSaOhoaEICwtTavPw8KjSd8zk4+l/gxACX3/9Nbp3717iNF98D57Hx9MTEVUcOt/D0bhxY9y5c0f6++OPP3RdUpnj4+m1ezx9v379sHTpUum8j9zcXPz1118AgG7dumHNmjUAnp3P8csvvxRbBx9PT0RUcej8stjq1avD3t5e9vmU9eWr2uDj6bV7PP306dORl5eHNm3aSHswpk+fjsaNG2P58uUYOXIkPD09Ubt2baXHxb+Ij6cnIqo4dPp4+tDQUCxcuBCWlpYwMjKCr68vwsPDpXMGXpSXlyf9mgSePRbX2dmZj6d/TlV7PP20adNgZmam9koVPp6eSsLH0xO9Gm0eT6/TQypt2rRBVFQUfv/9d6xcuRLJycnw8/NDdna22v7h4eGwtLSU/pydncu54orr9u3baNiwIeLi4pQuT63Knn88va7x8fREVNXpdA/HizIyMuDi4oIlS5bg/fffVxnOPRxUGXF71R3u4SB6Ndrs4dD5ORzPs7KyQoMGDaQT9l5kaGio9n4OREREVLHp/CqV5+Xk5ODatWvSfSZeVQXaeUNULG6nRFQV6HQPx7Rp09CnTx+4uLjg9u3bmDNnDvT09DBs2LBXmq6+vj4UCgXS0tJga2v70ns1EOmSEAJpaWlQKBTQ19fXdTlERLLRaeC4efMmhg0bhvv378PW1hYdOnTAiRMnSrxjY0n09PTg5OSEmzdvavTUVSJdUigUcHJyUns3ViKiykKngaOkB5W9CjMzM9SvXx9PnjyRbR5EZUFfX59hg4gqvQp10mhZ09PT4xc5ERFRBVChTholIiKiyomBg4iIiGTHwEFERESyY+AgIiIi2TFwEBERkewYOIiIiEh2DBxEREQkOwYOIiIikh0DBxEREcmOgYOIiIhkx8BBREREsmPgICIiItkxcBAREZHsGDiIiIhIdgwcREREJDsGDiIiIpIdAwcRERHJjoGDiIiIZMfAQURERLJj4CAiIiLZMXAQERGR7Bg4iIiISHYMHERERCQ7Bg4iIiKSHQMHERERya66Jp2aNWsGhUKh0QTj4uJeqSAiIiKqfDQKHP369ZP+/fjxY6xYsQKenp7w9fUFAJw4cQJ//fUXJkyYIEuRRERE9HrTKHDMmTNH+vfo0aMRHByMzz77TKXPjRs3Sl3IF198gZkzZ2LKlClYtmxZqadDREREFY9GgeN5P//8M86cOaPS/s4776Bly5ZYs2aN1kWcPn0aq1atQpMmTbQel4iISiHUUtcV6E5opq4rqJK0PmnU2NgYR48eVWk/evQojIyMtC4gJycHb7/9NlavXg1ra2utxyciIqKKT+s9HFOnTsV//vMfxMXFoXXr1gCAkydPYs2aNfj000+1LmDixIkIDAxEt27dMG/evJf2zcvLQ15envQ6KytL6/kRERFR+dM6cMyYMQN169bF8uXLsW7dOgBAo0aNEBkZicGDB2s1rY0bNyIuLg6nT5/WqH94eDjCwsK0LblUXGfsKpf5VEQpXwTqugQiIqpktA4cADB48GCtw8WLbty4gSlTpmDfvn0aH4qZOXMmQkJCpNdZWVlwdnZ+pTqIiIhIfqUKHACQn5+P1NRUFBYWKrXXqVNHo/HPnj2L1NRUNG/eXGorKCjA4cOH8c033yAvLw96enpK4xgaGsLQ0LC0JRMREZGOaB04rly5gvfeew/Hjh1TahdCQKFQoKCgQKPpdO3aFQkJCUpto0aNQsOGDTF9+nSVsEFERESvL60Dx8iRI1G9enXs3LkTDg4OGt+B9EXm5ubw8vJSajM1NUWNGjVU2omIiOj1pnXgiI+Px9mzZ9GwYUM56iEiIqJKSOvA4enpiXv37slRC2JjY2WZLhEREemW1jf+WrBgAT766CPExsbi/v37yMrKUvojIiIiepHWezi6desG4NlJn8/T9qRRIiIiqjq0DhwxMTFy1EFERFQuvKO9dV2CTiQEJZTcSUZaB45OnTrJUQcRERFVYqW+8Vdubi7++ecf5OfnK7Xzia9ERET0Iq0DR1paGkaNGoXdu3erHc5zOIiIiOhFpXpabEZGBk6ePInOnTtj27Zt+PfffzFv3jwsXrxYjhqpvIVa6roC3QjN1HUFRESVltaB4+DBg9ixYwdatmyJatWqwcXFBW+88QYsLCwQHh6OwEA+aZSIiIiUaX0fjocPH8LOzg4AYG1tjbS0NACAt7c34uLiyrY6IiIiqhS0DhweHh64dOkSAKBp06ZYtWoVbt26hYiICDg4OJR5gURERPT60/qQypQpU3Dnzh0AwJw5c9CjRw+sX78eBgYGiIqKKuv6iIiIqBLQOnC888470r9btGiB69evIykpCXXq1EHNmjXLtDgiIiKqHEp9H44iJiYmaN68eVnUQkRERJWU1oFDCIHNmzcjJiYGqampKCwsVBq+devWMiuOiIiIKodS3Ydj1apV8Pf3R61ataBQKOSoi4iIiCoRrQPH2rVrsXXrVvTq1UuOeoiIiKgS0vqyWEtLS9StW1eOWoiIiKiS0jpwhIaGIiwsDI8ePZKjHiIiIqqEtD6kMnjwYPz444+ws7ODq6sr9PX1lYbzbqNERET0Iq0DR1BQEM6ePYt33nmHJ40SERGRRrQOHLt27cKePXvQoUMHOeohIiKiSkjrczicnZ1hYWEhRy1ERERUSWkdOBYvXoyPPvoIKSkpMpRDRERElVGpnqWSm5sLd3d3mJiYqJw0+uDBgzIrjoiIiCoHrQPHsmXLZCiDSPe8o711XYLOJAQl6LoEIqrkSnWVChEREZE2tD6HAwCuXbuGTz75BMOGDUNqaioAYPfu3fjrr7/KtDgiIiKqHLQOHIcOHYK3tzdOnjyJrVu3IicnBwBw/vx5zJkzp8wLJCIiotef1oFjxowZmDdvHvbt2wcDAwOpvUuXLjhx4oRW01q5ciWaNGkCCwsLWFhYwNfXF7t379a2JCIiIqrgtA4cCQkJ6N+/v0q7nZ0d7t27p9W0nJyc8MUXX+Ds2bM4c+YMunTpgjfffJOHZoiIiCoZrU8atbKywp07d+Dm5qbUfu7cOdSuXVurafXp00fp9eeff46VK1fixIkTaNy4sUr/vLw85OXlSa+zsrK0mh8RERHphtZ7OIYOHYrp06fj7t27UCgUKCwsxNGjRzFt2jSMGDGi1IUUFBRg48aNePjwIXx9fdX2CQ8Ph6WlpfTn7Oxc6vkRERFR+dE6cMyfPx8NGzaEs7MzcnJy4OnpiY4dO6Jdu3b45JNPtC4gISEBZmZmMDQ0xPjx47Ft2zZ4enqq7Ttz5kxkZmZKfzdu3NB6fkRERFT+tD6kYmBggNWrV+PTTz9FYmIicnJy0KxZM9SvX79UBXh4eCA+Ph6ZmZnYvHkzgoKCcOjQIbWhw9DQEIaGhqWaDxEREemO1oGjSJ06dVCnTp1XLsDAwAD16tUDALRo0QKnT5/G8uXLsWrVqleeNhEREVUMGgWOkJAQjSe4ZMmSUhcDAIWFhUonhhIREdHrT6PAce7cOaXXcXFxePr0KTw8PAAAly9fhp6eHlq0aKHVzGfOnImePXuiTp06yM7OxoYNGxAbG4s9e/ZoNR0iIiKq2DQKHDExMdK/lyxZAnNzc0RHR8Pa2hoAkJ6ejlGjRsHPz0+rmaempmLEiBG4c+cOLC0t0aRJE+zZswdvvPGGVtMhIiKiik3rczgWL16MvXv3SmEDAKytrTFv3jx0794dH374ocbT+t///qft7ImIiOg1pPVlsVlZWUhLS1NpT0tLQ3Z2dpkURURERJWL1oGjf//+GDVqFLZu3YqbN2/i5s2b2LJlC95//30MGDBAjhqJiIjoNaf1IZWIiAhMmzYNw4cPx5MnT55NpHp1vP/++1i4cGGZF0hERESvP60Dh4mJCVasWIGFCxfi2rVrAAB3d3eYmpqWeXFERERUOZT6xl+mpqZo0qRJWdZCRERElZTW53AQERERaYuBg4iIiGTHwEFERESyY+AgIiIi2Wl90ujBgwexdetWpKSkQKFQwM3NDQMHDkTHjh3lqI+IiIgqAa32cIwfPx7dunXDjz/+iPv37yMtLQ3r16+Hv78/Jk+eLFeNRERE9JrTOHBs27YNkZGRWLNmDe7du4fjx4/jxIkTSEtLw+rVq/Hdd9/hl19+kbNWIiIiek1pHDgiIyMREhKCkSNHQqFQ/N8EqlXDe++9h6lTp/JhbERERKSWxoEjLi4O/fv3L3b4gAEDcPbs2TIpioiIiCoXjQPHvXv34OTkVOxwJycn3L9/v0yKIiIiospF48CRn58PfX39YodXr14d+fn5ZVIUERERVS5aXRb76aefwsTERO2w3NzcMimIiIiIKh+NA0fHjh1x6dKlEvsQERERvUjjwBEbGytjGURERFSZlfrW5vfu3cO9e/fKshYiIiKqpLQKHBkZGZg4cSJq1qyJWrVqoVatWqhZsyYmTZqEjIwMmUokIiKi153Gh1QePHgAX19f3Lp1C2+//TYaNWoEALhw4QKioqJw4MABHDt2DNbW1rIVS0RERK8njQPH3LlzYWBggGvXrqFWrVoqw7p37465c+di6dKlZV4kERERvd40PqSyfft2LFq0SCVsAIC9vT2+/PJLbNu2rUyLIyIiospB48Bx584dNG7cuNjhXl5euHv3bpkURURERJWLxoGjZs2aSElJKXZ4cnIybGxsyqImIiIiqmQ0DhwBAQGYNWuW2tuX5+Xl4dNPP0WPHj3KtDgiIiKqHLQ6abRly5aoX78+Jk6ciIYNG0IIgYsXL2LFihXIy8vD2rVr5ayViIiIXlMaBw4nJyccP34cEyZMwMyZMyGEAAAoFAq88cYb+Oabb+Ds7KzVzMPDw7F161YkJSXB2NgY7dq1w4IFC+Dh4aHdUhAREVGFptXD29zc3LB7926kp6fjypUrAIB69eqV+tyNQ4cOYeLEiWjVqhWePn2Kjz/+GN27d8eFCxdgampaqmkSERFRxaNV4ChibW2N1q1bv/LMf//9d6XXUVFRsLOzw9mzZ9U+CC4vLw95eXnS66ysrFeugYiIiORX6mepyCEzMxMAit1jEh4eDktLS+lP20M4REREpBsVJnAUFhZi6tSpaN++Pby8vNT2mTlzJjIzM6W/GzdulHOVREREVBqlOqQih4kTJyIxMRF//PFHsX0MDQ1haGhYjlURERFRWagQgWPSpEnYuXMnDh8+DCcnJ12XQ0RERGVMp4FDCIHJkydj27ZtiI2NhZubmy7LISIiIpnoNHBMnDgRGzZswI4dO2Bubi49i8XS0hLGxsa6LI2IiIjKkE5PGl25ciUyMzPRuXNnODg4SH+bNm3SZVlERERUxnR+SIWIiIgqvwpzWSwRERFVXgwcREREJDsGDiIiIpIdAwcRERHJjoGDiIiIZMfAQURERLJj4CAiIiLZMXAQERGR7Bg4iIiISHYMHERERCQ7Bg4iIiKSHQMHERERyY6Bg4iIiGTHwEFERESyY+AgIiIi2TFwEBERkewYOIiIiEh2DBxEREQkOwYOIiIikh0DBxEREcmOgYOIiIhkx8BBREREsmPgICIiItkxcBAREZHsGDiIiIhIdgwcREREJDsGDiIiIpKdTgPH4cOH0adPHzg6OkKhUGD79u26LIeIiIhkotPA8fDhQzRt2hTffvutLssgIiIimVXX5cx79uyJnj176rIEIiIiKgc6DRzaysvLQ15envQ6KytLh9UQERGRpl6rk0bDw8NhaWkp/Tk7O+u6JCIiItLAaxU4Zs6ciczMTOnvxo0bui6JiIiINPBaHVIxNDSEoaGhrssgIiIiLb1WeziIiIjo9aTTPRw5OTm4evWq9Do5ORnx8fGwsbFBnTp1dFgZERERlSWdBo4zZ87A399feh0SEgIACAoKQlRUlI6qIiIiorKm08DRuXNnCCF0WQIRERGVA57DQURERLJj4CAiIiLZMXAQERGR7Bg4iIiISHYMHERERCQ7Bg4iIiKSHQMHERERyY6Bg4iIiGTHwEFERESyY+AgIiIi2TFwEBERkewYOIiIiEh2DBxEREQkOwYOIiIikh0DBxEREcmOgYOIiIhkx8BBREREsmPgICIiItkxcBAREZHsGDiIiIhIdgwcREREJDsGDiIiIpIdAwcRERHJjoGDiIiIZMfAQURERLJj4CAiIiLZMXAQERGR7Bg4iIiISHYVInB8++23cHV1hZGREdq0aYNTp07puiQiIiIqQzoPHJs2bUJISAjmzJmDuLg4NG3aFAEBAUhNTdV1aURERFRGdB44lixZgjFjxmDUqFHw9PREREQETExMsGbNGl2XRkRERGWkui5nnp+fj7Nnz2LmzJlSW7Vq1dCtWzccP35cpX9eXh7y8vKk15mZmQCArKysMq+tMC+3zKf5ushSCF2XoBMFjwp0XYLOyPEZeh3wc141VdXPuhyf86JpClHy9qTTwHHv3j0UFBSgVq1aSu21atVCUlKSSv/w8HCEhYWptDs7O8tWY1VkqesCdOairgvQGcv/VN21XlVV7TVeNT/rcn7Os7OzYWn58unrNHBoa+bMmQgJCZFeFxYW4sGDB6hRowYUCoUOK6OykJWVBWdnZ9y4cQMWFha6LoeIZMDPeeUihEB2djYcHR1L7KvTwFGzZk3o6enh33//VWr/999/YW9vr9Lf0NAQhoaGSm1WVlZylkg6YGFhwS8iokqOn/PKo6Q9G0V0etKogYEBWrRogQMHDkhthYWFOHDgAHx9fXVYGREREZUlnR9SCQkJQVBQEFq2bInWrVtj2bJlePjwIUaNGqXr0oiIiKiM6DxwDBkyBGlpaZg9ezbu3r0LHx8f/P777yonklLlZ2hoiDlz5qgcNiOiyoOf86pLITS5loWIiIjoFej8xl9ERERU+TFwEBERkewYOIiIiEh2DBxEREQkOwYO0rnDhw+jT58+cHR0hEKhwPbt23VdEhGVsfDwcLRq1Qrm5uaws7NDv379cOnSJV2XReWIgYN07uHDh2jatCm+/fZbXZdCRDI5dOgQJk6ciBMnTmDfvn148uQJunfvjocPH+q6NConvCyWKhSFQoFt27ahX79+ui6FiGSUlpYGOzs7HDp0CB07dtR1OVQOuIeDiIjKXWZmJgDAxsZGx5VQeWHgICKiclVYWIipU6eiffv28PLy0nU5VE50fmtzIiKqWiZOnIjExET88ccfui6FyhEDBxERlZtJkyZh586dOHz4MJycnHRdDpUjBg4iIpKdEAKTJ0/Gtm3bEBsbCzc3N12XROWMgYN0LicnB1evXpVeJycnIz4+HjY2NqhTp44OKyOisjJx4kRs2LABO3bsgLm5Oe7evQsAsLS0hLGxsY6ro/LAy2JJ52JjY+Hv76/SHhQUhKioqPIviIjKnEKhUNseGRmJkSNHlm8xpBMMHERERCQ7XhZLREREsmPgICIiItkxcBAREZHsGDiIiIhIdgwcREREJDsGDiIiIpIdAwcRERHJjoGDiIiIZMfAQVRFxcbGQqFQICMjQ9elKHF1dcWyZcuk1wqFAtu3b5dtfiNHjkS/fv1kmz4RPcNnqRBVAZ07d4aPj4/Sf+Tt2rXDnTt3YGlpqbvCNHDnzh1YW1vrugwiekUMHERVlIGBAezt7XVdRolehxqJqGQ8pEJUyY0cORKHDh3C8uXLoVAooFAokJKSonJIJSoqClZWVti5cyc8PDxgYmKCgQMHIjc3F9HR0XB1dYW1tTWCg4NRUFAgTT8vLw/Tpk1D7dq1YWpqijZt2iA2NrbYeoQQCA0NRZ06dWBoaAhHR0cEBwcX2//FQyo3b97EsGHDYGNjA1NTU7Rs2RInT56Uhu/YsQPNmzeHkZER6tati7CwMDx9+rTE92nRokVwcHBAjRo1MHHiRDx58qTEcYhIc9zDQVTJLV++HJcvX4aXlxfmzp0LALC1tUVKSopK39zcXHz11VfYuHEjsrOzMWDAAPTv3x9WVlb47bff8Pfff+Ott95C+/btMWTIEADApEmTcOHCBWzcuBGOjo7Ytm0bevTogYSEBNSvX19lHlu2bMHSpUuxceNGNG7cGHfv3sX58+c1WpacnBx06tQJtWvXxi+//AJ7e3vExcWhsLAQAHDkyBGMGDECX331Ffz8/HDt2jWMHTsWADBnzpxipxsTEwMHBwfExMTg6tWrGDJkCHx8fDBmzBiN6iIiDQgiqvQ6deokpkyZotQWExMjAIj09HQhhBCRkZECgLh69arUZ9y4ccLExERkZ2dLbQEBAWLcuHFCCCGuX78u9PT0xK1bt5Sm3bVrVzFz5ky1tSxevFg0aNBA5Ofnqx3u4uIili5dKr0GILZt2yaEEGLVqlXC3Nxc3L9/X+24Xbt2FfPnz1dqW7t2rXBwcFDbXwghgoKChIuLi3j69KnUNmjQIDFkyJBixyEi7XEPBxFJTExM4O7uLr2uVasWXF1dYWZmptSWmpoKAEhISEBBQQEaNGigNJ28vDzUqFFD7TwGDRqEZcuWoW7duujRowd69eqFPn36oHr1kr+O4uPj0axZM9jY2Kgdfv78eRw9ehSff/651FZQUIDHjx8jNzcXJiYmasdr3Lgx9PT0pNcODg5ISEgosR4i0hwDBxFJ9PX1lV4rFAq1bUWHMHJycqCnp4ezZ88q/YcNQCmkPM/Z2RmXLl3C/v37sW/fPkyYMAELFy7EoUOHVOb1ImNj45cOz8nJQVhYGAYMGKAyzMjIqNjxXraMRFQ2GDiIqgADAwOlEz3LSrNmzVBQUIDU1FT4+flpPJ6xsTH69OmDPn36YOLEiWjYsCESEhLQvHnzl47XpEkTfP/993jw4IHavRzNmzfHpUuXUK9ePa2XhYjkxcBBVAW4urri5MmTSElJgZmZWbGHJLTVoEEDvP322xgxYgQWL16MZs2aIS0tDQcOHECTJk0QGBioMk5UVBQKCgrQpk0bmJiYYN26dTA2NoaLi0uJ8xs2bBjmz5+Pfv36ITw8HA4ODjh37hwcHR3h6+uL2bNno3fv3qhTpw4GDhyIatWq4fz580hMTMS8efPKZJmJqHR4WSxRFTBt2jTo6enB09MTtra2+Oeff8ps2pGRkRgxYgQ+/PBDeHh4oF+/fjh9+jTq1Kmjtr+VlRVWr16N9u3bo0mTJti/fz9+/fXXYs/5eJ6BgQH27t0LOzs79OrVC97e3vjiiy+kwzkBAQHYuXMn9u7di1atWqFt27ZYunSpRmGGiOSlEEIIXRdBRERElRv3cBAREZHsGDiIiIhIdgwcREREJDsGDiIiIpIdAwcRERHJjoGDiIiIZMfAQURERLJj4CAiIiLZMXAQERGR7Bg4iIiISHYMHERERCS7/weFrTu71EE4YgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 550x340 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(5.5, 3.4))\n",
    "slices = np.arange(1, 3)\n",
    "width = 0.25\n",
    "ax.bar(slices - width, truth_a2.reshape(-1), width, label=\"truth\")\n",
    "ax.bar(slices, sim_a2, width, label=\"od-dynamic-sim (exact)\")\n",
    "ax.bar(slices + width, seq_a2, width, label=\"od-dynamic-seq (exact)\")\n",
    "ax.set_xticks(slices)\n",
    "ax.set_xlabel(\"time slice h\")\n",
    "ax.set_ylabel(\"OD demand\")\n",
    "ax.set_title(\"Anchor A2: simultaneous strictly dominates sequential\")\n",
    "ax.legend(fontsize=8)\n",
    "fig.tight_layout()\n",
    "display(fig)\n",
    "plt.close(fig)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c9feef6",
   "metadata": {},
   "source": [
    "## Takeaways & pointers\n",
    "\n",
    "- **Certified, not self-reported.** Anchors A1/A2 were recomputed as exact\n",
    "  rationals in-cell; the end-to-end two-route run was certified by\n",
    "  `DynamicODCertifier` inside `run_dynamic_estimation_experiment`, never\n",
    "  from a self-report.\n",
    "- **The efficiency ordering is real, not asymptotic folklore.** Anchor A2\n",
    "  proves it exactly on a hand-built tensor; the end-to-end run confirms the\n",
    "  SAME ordering (`sim < seq < prior`) on a real assignment-derived task.\n",
    "- **`prior-profile` is a first-class competitor.** It is the baseline both\n",
    "  GLS variants must beat (`ADR-023` `covers: [od-dynamic-sim, od-dynamic-seq,\n",
    "  prior-profile]`).\n",
    "- **Where next.** the static-T2 sibling this reduces to under integer lag\n",
    "  [`gls`](01-gls.ipynb) (Anchor A1, verified exactly above); the\n",
    "  spatial-covariance estimator [`od-kalman`](06-od-kalman.ipynb); the\n",
    "  lineage 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": [
    "od-dynamic-sim",
    "od-dynamic-seq",
    "prior-profile"
   ],
   "requires_extra": null,
   "track": "estimation",
   "unit": "od-dynamic"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}