{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "50e781db",
   "metadata": {},
   "source": [
    "# `od-kalman` — Davis & Nihan's (1993) linear-Gaussian OD estimation\n",
    "\n",
    "**What.** `od-kalman` (`DavisNihanKalmanEstimator`, T2, ADR-012) recovers OD\n",
    "demand from a TIME SERIES of link counts under Davis & Nihan's proven\n",
    "large-population limit: the day-to-day link-count process is a stationary\n",
    "linear-Gaussian VAR/VARMA around the equilibrium loading (their Prop. 2 for\n",
    "the mean, Prop. 3 for the covariance). `od-kalman` whitens its GLS solve by\n",
    "the EXACT multinomial spatial covariance across sensors (route-sharing links\n",
    "correlate positively, competing-route links negatively) plus an AR(1)\n",
    "effective-sample-size correction `tau` for day-to-day persistence — both\n",
    "strictly richer than `gls`'s IID-count assumption.\n",
    "\n",
    "**Why it is in the benchmark.** It needs its own observation model\n",
    "(`DayToDayCounts`, a VAR(1) count series centered on the UE loading, ADR-012)\n",
    "and is the benchmark's only T2 estimator that uses the CROSS-LINK spatial\n",
    "structure of counts rather than treating each sensor as independent — a\n",
    "mathematically distinct estimator from `gls`/`od-congested`, not a rename.\n",
    "See the [model compendium](../../docs/MODELS.md) (Davis & Nihan 1993) and\n",
    "[docs/design/adr-012-dn-kalman.md](../../docs/design/adr-012-dn-kalman.md)\n",
    "(P1).\n",
    "\n",
    "**Scope.** The exact multinomial spatial-covariance closed form, the DN-GLS\n",
    "scalar closed form, the off-diagonal whitening's load-bearing effect, the\n",
    "AR(1) `tau` persistence correction, and an end-to-end certified recovery on\n",
    "Braess under a day-to-day count series.\n",
    "\n",
    "**Canon.** `[davis1993large]`, [docs/REFERENCES.md](../../docs/REFERENCES.md) / [docs/references.bib](../../docs/references.bib)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35660148",
   "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 closed forms are recomputed\n",
    "algebraically in-cell (no trusted digits), and the end-to-end recovery is\n",
    "recomputed by the P1 `ODCertifier` (via `run_estimation_experiment`) from\n",
    "the emitted OD matrix against the harness's own pinned BFW assignment, never\n",
    "from `od-kalman`'s self-report\n",
    "([README](../../README.md), *Certified, not self-reported*)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "d8d71f2e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T17:05:21.009746Z",
     "iopub.status.busy": "2026-07-21T17:05:21.009296Z",
     "iopub.status.idle": "2026-07-21T17:05:23.365346Z",
     "shell.execute_reply": "2026-07-21T17:05:23.364331Z"
    }
   },
   "outputs": [],
   "source": [
    "# Setup. `od-kalman` is a core estimator: 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",
    "    Budget,\n",
    "    Demand,\n",
    "    RngBundle,\n",
    "    Trace,\n",
    "    braess_scenario,\n",
    "    run_estimation_experiment,\n",
    "    two_route_scenario,\n",
    "    viz,\n",
    ")\n",
    "from tabench.core.rng import SOURCE_OBSERVATION\n",
    "from tabench.estimation import ODTrace, ar1_tau, dn_gls_solve, gls_solve\n",
    "from tabench.estimation import DavisNihanKalmanEstimator, GLSEstimator\n",
    "from tabench.estimation.base import EstimationTask\n",
    "from tabench.models.frank_wolfe import BiconjugateFrankWolfeModel\n",
    "from tabench.observe._dn_process import dn_spatial_covariance\n",
    "from tabench.observe.levels import DayToDayCounts\n",
    "\n",
    "\n",
    "def _dn_task(scenario, truth, sensors, prior_d, n_periods, scale, rho, seed=0):\n",
    "    m = np.zeros((2, 2))\n",
    "    m[0, 1] = prior_d\n",
    "    ds = DayToDayCounts(np.asarray(sensors), n_periods, scale, rho, k_inner=80).observe(\n",
    "        scenario, truth, RngBundle(seed).generator(SOURCE_OBSERVATION)\n",
    "    )\n",
    "    return EstimationTask(\n",
    "        name=\"t\", network=scenario.network, prior=Demand(m), dataset=ds,\n",
    "        identifiability={}, scenario_hash=scenario.content_hash(), seed=seed,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c9641e7",
   "metadata": {},
   "source": [
    "## The DN spatial covariance: an exact multinomial closed form\n",
    "\n",
    "On the two-route corridor, routes `A={link 0, link 1}` and `B={link 2, link\n",
    "3}`: the DN link-count covariance is exactly the multinomial form `Var(link\n",
    "0) = (D^2/N) p_A p_B`, with same-route links POSITIVELY correlated and\n",
    "across-route links NEGATIVELY correlated — recomputed here directly, not\n",
    "quoted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "47f2dcc1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T17:05:23.369484Z",
     "iopub.status.busy": "2026-07-21T17:05:23.369107Z",
     "iopub.status.idle": "2026-07-21T17:05:23.401296Z",
     "shell.execute_reply": "2026-07-21T17:05:23.399896Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Var(link 0) DN            : 0.018750\n",
      "closed-form (D^2/N) p_A p_B : 0.018750\n",
      "Cov(link 0, link 1) [same route A]  : 0.018750  (should equal Var)\n",
      "Cov(link 0, link 2) [A vs B]        : -0.018750  (should equal -Var)\n"
     ]
    }
   ],
   "source": [
    "sc2 = two_route_scenario(sue_theta=None)\n",
    "truth_trace = Trace()\n",
    "BiconjugateFrankWolfeModel().solve(\n",
    "    sc2, Budget(iterations=5000, target_relative_gap=1e-10), RngBundle(0), truth_trace\n",
    ")\n",
    "truth2 = truth_trace.final.link_flows\n",
    "demand = float(sc2.demand.matrix[0, 1])\n",
    "p_a = truth2[0] / demand  # recomputed route-A proportion\n",
    "p_b = 1.0 - p_a\n",
    "n_trav = np.array([round(50.0 * demand)])\n",
    "q = dn_spatial_covariance(\n",
    "    sc2.network, sc2.demand, np.array([demand]), n_trav, 80, pairs=[(0, 1)]\n",
    ")\n",
    "var0 = (demand**2 / n_trav[0]) * p_a * p_b  # closed form, recomputed\n",
    "print(f\"Var(link 0) DN            : {q[0, 0]:.6f}\")\n",
    "print(f\"closed-form (D^2/N) p_A p_B : {var0:.6f}\")\n",
    "print(f\"Cov(link 0, link 1) [same route A]  : {q[0, 1]:.6f}  (should equal Var)\")\n",
    "print(f\"Cov(link 0, link 2) [A vs B]        : {q[0, 2]:.6f}  (should equal -Var)\")\n",
    "assert np.isclose(q[0, 0], var0, atol=1e-9)\n",
    "assert np.isclose(q[0, 1], var0, atol=1e-9)   # same route -> positively correlated\n",
    "assert np.isclose(q[0, 2], -var0, atol=1e-9)  # A vs B -> negatively correlated\n",
    "assert np.allclose(q, q.T)  # symmetric covariance"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd6bdb63",
   "metadata": {},
   "source": [
    "## The scalar closed form\n",
    "\n",
    "Single pair, single sensor: `g* = (g_pr/w^2 + p*c/s^2) / (1/w^2 + p^2/s^2)`\n",
    "— the DN whitening's GLS solve at explicit prior/count variances."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "8dc3bc70",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T17:05:23.404433Z",
     "iopub.status.busy": "2026-07-21T17:05:23.404196Z",
     "iopub.status.idle": "2026-07-21T17:05:23.410543Z",
     "shell.execute_reply": "2026-07-21T17:05:23.409747Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dn_gls_solve   : 3.757576\n",
      "closed form    : 3.757576\n"
     ]
    }
   ],
   "source": [
    "p, c, g_pr, s2, w2 = 0.625, 2.5, 3.0, 0.5, 4.0\n",
    "expected = (g_pr / w2 + p * c / s2) / (1.0 / w2 + p * p / s2)\n",
    "got = dn_gls_solve(\n",
    "    np.array([[p]]), np.array([c]), np.array([g_pr]), np.array([[s2]]), np.array([w2])\n",
    ")\n",
    "print(f\"dn_gls_solve   : {got[0]:.6f}\")\n",
    "print(f\"closed form    : {expected:.6f}\")\n",
    "assert np.isclose(got[0], expected, atol=1e-10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30981eed",
   "metadata": {},
   "source": [
    "## Off-diagonal whitening is load-bearing\n",
    "\n",
    "Two sensors on the SAME OD pair whose counts are correlated are NOT the two\n",
    "INDEPENDENT measurements `gls`'s diagonal covariance assumes: with a\n",
    "diagonal `Sigma`, `dn_gls_solve` reduces to plain `gls_solve` exactly, but a\n",
    "correlated (non-diagonal) `Sigma` moves the estimate materially — the\n",
    "spatial sense in which `od-kalman` is not a `gls` rename."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "dbb53754",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T17:05:23.413383Z",
     "iopub.status.busy": "2026-07-21T17:05:23.413190Z",
     "iopub.status.idle": "2026-07-21T17:05:23.419670Z",
     "shell.execute_reply": "2026-07-21T17:05:23.418888Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "diagonal Sigma   : 4.823944  (via plain gls_solve: 4.823944)\n",
      "correlated Sigma : 4.772727\n"
     ]
    }
   ],
   "source": [
    "p_obs = np.array([[0.6], [0.4]])  # two sensors, one pair\n",
    "counts = np.array([3.0, 2.0])\n",
    "prior = np.array([4.0])\n",
    "w_var = np.array([9.0])\n",
    "var = np.array([1.0, 1.0])\n",
    "diagonal = dn_gls_solve(p_obs, counts, prior, np.diag(var), w_var)\n",
    "correlated = dn_gls_solve(p_obs, counts, prior, np.array([[1.0, 0.8], [0.8, 1.0]]), w_var)\n",
    "via_gls = gls_solve(p_obs, counts, prior, w_var, var)\n",
    "print(f\"diagonal Sigma   : {diagonal[0]:.6f}  (via plain gls_solve: {via_gls[0]:.6f})\")\n",
    "print(f\"correlated Sigma : {correlated[0]:.6f}\")\n",
    "assert np.isclose(diagonal[0], via_gls[0], atol=1e-8)\n",
    "assert abs(correlated[0] - diagonal[0]) > 0.05"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a6fb6a9",
   "metadata": {},
   "source": [
    "## The AR(1) temporal correction `tau`\n",
    "\n",
    "Day-to-day persistence inflates the effective sampling variance: `tau`\n",
    "tracks the AR(1) variance-inflation factor `(1+rho)/(1-rho)` and collapses\n",
    "to `1` for an IID (`rho=0`) series — the effective-sample-size knob that\n",
    "distinguishes `od-kalman` from every mean-collapsing estimator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "2ec8892c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T17:05:23.422491Z",
     "iopub.status.busy": "2026-07-21T17:05:23.422114Z",
     "iopub.status.idle": "2026-07-21T17:05:23.512559Z",
     "shell.execute_reply": "2026-07-21T17:05:23.511810Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rho=0.0: tau=1.0000  target (1+rho)/(1-rho)=1.0000  rho_hat=0.0000\n",
      "rho=0.6: tau=3.8868  target (1+rho)/(1-rho)=4.0000  rho_hat=0.5907\n"
     ]
    }
   ],
   "source": [
    "for rho in (0.0, 0.6):\n",
    "    ds = DayToDayCounts(np.array([0]), 4000, 50.0, rho, k_inner=80).observe(\n",
    "        sc2, truth2, RngBundle(0).generator(SOURCE_OBSERVATION)\n",
    "    )\n",
    "    tau, rho_hat = ar1_tau(ds.payload[\"counts\"])\n",
    "    target = (1.0 + rho) / (1.0 - rho)\n",
    "    print(f\"rho={rho}: tau={tau:.4f}  target (1+rho)/(1-rho)={target:.4f}  rho_hat={rho_hat:.4f}\")\n",
    "    assert abs(tau - target) < 0.15 * max(target, 1.0)\n",
    "    if rho == 0.0:\n",
    "        assert abs(tau - 1.0) < 0.15"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57339aa7",
   "metadata": {},
   "source": [
    "## End to end: a certified Braess recovery under a day-to-day count series\n",
    "\n",
    "`braess_scenario(6.0)` is frozen and content-hashed (P2). Run through the\n",
    "public `run_estimation_experiment` API with `noise=\"day_to_day\"` (the DN\n",
    "VAR(1) observation level): a `cv=0` prior (=truth) certifies feasible and\n",
    "recovers to the finite-population sampling floor from four observed links,\n",
    "with the middle link `3->4` reserved as the held-out set (an empty held-out\n",
    "set is rejected on both T2 tracks, ADR-002/ADR-023)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "14e61838",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T17:05:23.515420Z",
     "iopub.status.busy": "2026-07-21T17:05:23.515218Z",
     "iopub.status.idle": "2026-07-21T17:05:23.820911Z",
     "shell.execute_reply": "2026-07-21T17:05:23.819866Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "scenario      : braess\n",
      "content hash  : cf00f411cdccec88…\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "od_feasible : 1\n",
      "od_rmse     : 0.0003\n",
      "self_obs_count_rmse (provenance) : 0.0098\n"
     ]
    }
   ],
   "source": [
    "sc = braess_scenario(6.0)\n",
    "print(f\"scenario      : {sc.name}\")\n",
    "print(f\"content hash  : {sc.content_hash()[:16]}…\")\n",
    "\n",
    "cfg = {\n",
    "    \"sensors\": {\"kind\": \"explicit\", \"links\": [0, 1, 3, 4]},\n",
    "    \"heldout\": {\"kind\": \"explicit\", \"links\": [2]},\n",
    "    \"n_periods\": 60,\n",
    "    \"noise\": \"day_to_day\",\n",
    "    \"population_scale\": 80.0,\n",
    "    \"rho\": 0.5,\n",
    "    \"prior\": {\"kind\": \"stale\", \"cv\": 0.0},\n",
    "    \"identifiability_k_inner\": 40,\n",
    "}\n",
    "result = run_estimation_experiment(\n",
    "    sc, [DavisNihanKalmanEstimator(k_inner=60, outer_iters=15)],\n",
    "    Budget(sp_calls=5000), seed=0, macroreps=1, estimation=cfg,\n",
    ")\n",
    "row = [r for r in result.rows if r[\"estimator\"] == \"od-kalman\"][-1]\n",
    "print(f\"od_feasible : {row['od_feasible']:.0f}\")\n",
    "print(f\"od_rmse     : {row['od_rmse']:.4f}\")\n",
    "print(f\"self_obs_count_rmse (provenance) : {row['self_obs_count_rmse']:.4f}\")\n",
    "assert row[\"od_feasible\"] == 1.0\n",
    "assert row[\"od_rmse\"] < 0.5  # recovers near the planted truth at the DN sampling floor\n",
    "assert np.isfinite(float(row[\"self_obs_count_rmse\"]))\n",
    "assert result.manifest[\"estimation\"][\"noise\"] == \"day_to_day\"\n",
    "assert result.manifest[\"estimation\"][\"rho\"] == 0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1f8eef8",
   "metadata": {},
   "source": [
    "## Visualize\n",
    "\n",
    "`tabench.viz` draws the road `Network`'s link flows for the two-route\n",
    "recovery-from-DN-series result vs the truth (recomputed here at full-sensor\n",
    "noiseless-limit settings to keep the plotted OD matrix exactly matched to a\n",
    "certified recovery)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "d844dbc3",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-21T17:05:23.824284Z",
     "iopub.status.busy": "2026-07-21T17:05:23.823786Z",
     "iopub.status.idle": "2026-07-21T17:05:25.948988Z",
     "shell.execute_reply": "2026-07-21T17:05:25.947955Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAesAAAG0CAYAAAASKw+DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA/dklEQVR4nO3deVwU9f8H8NfA7solIKiAmgl4o5m3ImBipuWNinlllwqZpt9vGV+/HnhbX7M8SvP6VWaHmamZZYpmeJt0aIUHoCUK3qCC7MH+/jAo5Bp2Zndnd15PH/sodmc+8+HaF+/P5zMzgtFoMIOIiIgUy8XeHSAiIqKKMayJiIgUjmFNRESkcAxrIiIihWNYExERKRzDmoiISOEY1kRERArHsCYiIlI4hjUREZHCMayJrGTBggWIiIgo/jg+Ph7Dhw+vUhs+Pr7Yvn17lfb5v/97D82bh8HXtwbeeeedUv0gIsfDsCbRevfujYSEBHt3Q1a2DLKFCxfinXfeseoxcnNz8corr2DSpJeQmvo7nn76aasej4hsQ2PvDhBVldlshslkgkbjWD++Pj4+Vj/GhQsXYDAY8NhjjyEwMNDqxyMi22BlTaLEx8dj//4DWLFiJXx8fOHj44vz58+ja9dHsHTpsuLthg8fDn//mrh9+zYAIDMzEz4+vkhLSwcA3LhxE+PGjUP9+g8iMDAIgwYNRlpaWoXHTk5Oho+PL3bt2oWoqK6oVas2Dh06hIKCAkyZMgWhoQ1Ru3YAevbshePHU4r327BhA+rXr1+ire3bt8PHx7f49YULX8OJEyeLP6cNGzYAAG7evIkXX5yAkJBQ1Kv3APr06YsTJ05I/hr+cxi8d+/emDJlCqZPn4EHH2yARo0aY8GCBRW2MX/+fDRu3AQnT54s9dqGDRvQuXM4AKBVq4eLv0f3KywsxGuvvYZmzZqjVq3aiIiIwO7du4tfHzXqKbz88ivFHyckJMDHxxenT58GAOj1egQF1cHevd8BALZs2YrOncMREBCIBg2C0a9ff9y5c0f8F4aIKsWwJlEWLlyIDh06YPTo0Th9+hROnz6FevXqISKiC/bv3w/gXsV78OAh+Pj44PDhwwCAAwcOoE6dOggNDQEAvPBCPH788Sd88snH2LXrW5jNZgwePAQGg6HSPiQmJiIxMRFHjx5FWFgLzJgxA9u2fYmVK1fg++/3ISQkGDExMbh+/YaozykmJgYvvvgimjVrVvw5xcTEAABGj34aV69ewaZNm7Bv33do1aoV+vXrX9z2+fPn4ePji+Tk5Cp/Lf/p448/gaenB/bsScLs2bPw2muvY8+evaW2M5vNeOWVV/Dxx5/g66+/RosWLcr8fLZu3QoA2LNnT/H36H4rVqzA8uVvY86cOTh48ACio7vjySeHFf/R9M/vKQDs338A/v7+SE6+91xKSgoMBgM6duyArKwsPPfccxg5cgSOHj2Cr77ajr59+8Js5s38iOTEsCZRfHx8oNNp4eHhjoCAAAQEBMDV1RURERE4fPgQTCYTTp48CZ1Oh9jYIcVv7MnJ+9GlSxcAQFpaGnbs+BrLli1FeHg4WrZsiTVrVuPSpUvYvv2rSvswdepUREd3Q0hIMKpV02Ht2nWYM2c2evTogaZNm2Lp0qVwd3fH+vXrRX1O7u7u8PLyhEbjWvw5ubu749ChQ0hJScH777+PNm1aIzQ0FPPmzYWPj09xGGq1WjRq1AgeHh4WfkXvCQsLQ0JCAkJDQzFs2DC0bt0a+/btK7GN0WjCmDFjsW/f99i585viP3zK+nz8/GoAAGrW9C/+Ht1v2bLleOmllzB48CA0atQIs2fPQsuWLfHOOysAABEREUhNTcXVq1dx48ZNnDp1CvHxccUBnpy8H23atIGHhweysrJgNBrRt29fPPjggwgLC8OYMc/Dy8tL0teFiEpyrEk/UpzOnTvj1q3b+PnnX3D06BF06dIFERERePPNtwDcq6wnTpwIADh16hQ0Gg3atWtXvL+fnx8aNmyI06dPAQAGDRqMQ4cOAQAeeOABHDlyuHjb1q1bF/9/RkbGX9Vdx+LntFot2rZtU9yWpU6ePInbt28jOLhkKObn5yMjIwMAUKdOHfzwwzFJxwHuhfU/BQYG4OrVKyWemzp1KnQ6HZKSdsPf31/S8XJzc3Hp0iV06tSxxPOdOnXEiRP3htabN2+OGjVqYP/+A9DptHjooYfQs2dPrF69BsC972nRoryWLVuia9euCA/vgujoaERHR6N///6oUcNXUj+JqCSGNUni6+uLFi1aYP/+/Th69Ci6deuG8PAueOaZZ3H27FmkpaUhIqKL6PaWLVuK/Py7AACttuSPZ1WrWBcXF9w/GmswGCvd7/btOwgMDCzzlClfX3kXid3/OQqCgMLCwhLPdev2CDZt+hxJSUmIjY2V9fhlEQQB4eHh2L9/P6pV0yEiIgItWrRAQUEBfvvtNxw9ehQTJkwAALi6umLr1i04cuQI9uzZg1Wr3sWcOXOQlLQbDRo0sHpfidSCw+Akmlarg8lkKvV8REQXJCcn4+DBg4iMjICfXw00adIYixYtQmBgIBo2bAgAaNKkCYxGI3744Yfifa9fv46zZ8+iSZOmAFA8vx0aGlJqcdg/BQcHQ6fT4ciRI8XPGQwGpKT8WNxWzZo1cevWrRKLne5fJHbvcyoZjq1atUJ2djY0GtfivhQ9pFa2lnj88cexZs1qTJgwEZs2fS6pLW9vbwQFBeHw4SMlnj98+AiaNm1a/HHRvHVy8n5ERkbAxcUFXbqEY+nSpSgoKChRmQuCgE6dOmHq1KlITk6GTqer8rnhRFQxhjWJVr9+ffzww3GcP38e165dK64AIyIikJSUBI1Gg8aNGxc/t3HjZ8Xz1QAQGhqK3r2fwMSJL+HQoUM4ceIExowZi6CgIPTu/USV+uLp6YnnnnsW06fPwO7du5GamoqJEyciLy8Po0aNAgC0bdsOHh4emD17NtLTM/DZZ5/ho48+KtHOgw/Wx/nz5/HLL7/g2rVrKCgoQLduj6BDhw4YMWIEkpL24Pz58zhy5Ahmz56DlJQfAQAXL15Eu3btcfz4cQu/mlXTt29fvPvuuxg/fjy2bNkqqa2JEydgyZIl+PzzzThz5gxmzkzEiRMnEB8fV7xN0bx1amoqOnXqVPzcxo2foXXr1vD09AQA/PDDD1i06A2kpPyIP//8E9u2fYmrV6+iSZMmkvpIRCUxrEm0iRMnwNXVFR07dkJISCj+/PNPAEDnzuEoLCwsEcwREREwmUylLjjy9tvv4OGHW2Ho0KHo0eMxmM1mbNr0GbRabZX7k5iYiH79+mLs2HGIiuqK9PQMbN68uXi+1M+vBlatehfffrsL4eHh2LRpU6mLuvTr1w/du3dHnz59ERISik2bNkEQBHz22UaEh4dj/PjxaNu2HZ599jn8+eefqF27FoB7VfyZM2eQl5dX5X5basCA/lixYgXGjRuHbdu2WdxOXFwcxo9/AdOmTUPnzuFIStqNTz75GKGhocXbhIWFwcfHBy1btixeLFbW97R69eo4ePAghgwZgrZt22Hu3LmYN28uevToYfknSkSlCEajgedYEBERKRgrayIiIoVjWBMRESkcw5qIiEjhGNZEREQKx7AmIiJSOIY1ERGRwjGsiYiIFI5hTUREpHAMayIiIoVjWBMRESkcw5qIiEjhGNZEREQKx7AmIiJSOIY1ERGRwjGsiYiIFI5hTUREpHAMayIiIoVjWBMRESkcw5qIiEjhGNZEREQKx7AmIiJSOIY1ERGRwjGsiYiIFI5hTUREpHAaWx7sXx9fwr7UO7Y8ZLGuTT2xeFiQXY5NREQkhU0r632pd2AwmW15SACAwWS22x8JREREUtm0sgYArauA5GmNbXrMyLmnbXo8IiJby83NRXp6Gi5fvgyjwQCNVovatWsjJCQU3t7e9u5eKXq9HkajUXI7Go0GOp1Ohh4pm83DmoiI5JGTk4OkpN34bu8eZGVllbtdYGAgukV3R3R0d/j4+Niwh2XT6/UYMmwUqrlID2tfX18sf3uF0wc2w5qIyMHo9Xps3Pgpvtq+HSaTEU2aNEW3bt0REhqKunXrQafVQm8wIDPzAtLT0pCSchwff7QBGz/9FL379EFs7FC7hpvRaEQ1FyO+u94YRrPls7EaoRCP4DSMRiPDmoiIlCM9PR3Llr6FzMxMRERGYeCAGDQIDi5z24CAALRp0xaDh8TiXEYGvtiyGdu2bsHxH45hwsRJCAkJsXHvSzKaXWAyu9q1D46Cp24RETmI3379FYkzp+POnTwk/Oe/mDz53+UG9f0aBAdj8uR/IyFhKu7cyUPizOn47ddfrdzjSggyPFSCYU1E5ADS09OxYMF81Kjhh4Wv/Q/t23ewqJ32HTpi4Wv/u9fOwvlIT0+XuafiCTL8UwuGNRGRwun1eixb+hbc3d0xM3E2atWqJam9WrVqYWbibLi5uWPZ0reg1+tl6mlVCYAg4cGwJiIipdi48VNkZmYiLv4FyUFdpFatWoiLi0dmZiY2bvxUljbJehjWREQKlpOTg6+2b0dEZJTFQ9/lad+hI7pERGLHV9uRk5Mja9vicNJaLMWH9d27d/H0qGHo0qE1uncNx9BB/ZGRnmbvbtmM2Wz7K74RkXIkJe2GyWTEwAExovdZv349vLy88OWXX1a67cCBMTAajdizJ0lKNy3DrBZN8WENAKOeehr7j6Qgad9B9Hz8Cfx70gR7d8nq7hQUIuOKASczDSgwMLCJ1Oq7vXvQpElT0au+z58/j/feew8dOoirwoODQ9C4SRPstUdYk2iKD2s3Nzd079ETgnDvT6g2bdvjzz//sHOvrKcopNOvGHG74F5I640MayI1ys3NRVZWFtq0aStq+8LCQowfPx6LFi2q0kVC2rRui6ysLNy6dcvSrlqEq8HFU3xY32/NqhXo+fgT9u6G7MoKaSJSt/S/pvxCQkNFbb9s2TJ06tQJrVu3rtJxitpPt/UUo5SV4MUrwtXBoa5gtuTNRTiXkY6Nb1Q+D+NITmXpoa/gErmZN41wzVXPDyWRLVV3ExDoo8y3wsuXLwMA6tatV+m2v/76K7Zu3YqdO3dW+ThF7WdnZ1d5X2mkTjyr531RmT+hZVixfCl2bP8SGzdvhYeHh727I69KCmmzGSjkQjMiqzCblfuGbzQYAAA6rbbSbQ8ePIjz58+jVatWAO4F74QJE5CVlYUxY8ZUuK9Opy1xPFIehwjrle8sxxebN2Hj5q3w8fG1d3dk1yRIh9t3C5Gda0KevnQo16uhQXV3h5uxICKJNH+FtF5EiI4ZM6ZEKPfq1Qvjx49H3759K91XrzeUOJ7NKPfvJMVRfFhfvJiJWTOm4sEGDTB4QB8AgE6nw45v99q5Z/LycnOBl5tLhaFNROpSu3ZtAEBm5gUEBARY7TiZmRcAwKrHKIvURWJqWmCm+LCuU6cuLl3NtXc3bKas0NZq1PMDSUR/Cwn5a+FXWproFeFFvvnmG9HbpqellTgeKQ/HVhXKy80FobW1aFFXCzctw5pIjby9vREYGIiUlONWPU7Kj8cRGBiI6tWrW/U4pXA1uGgMa4UTVPTDSESlPdItGqdOpeJcRoZV2s/ISMfpU6fQLbq7VdqvGC9hJhbDmohIwbp3fxSurhp8sWWzVdr/4ovN0Gg0iLZLWJNYDGsiIgXz8fFB7z59sD/5exw7dlTWto8dPYID+5PxRO8+8PHxkbVtkhfDmohI4WJjh6Ju3bpYueIdXLlyRZY2r1y5gpUr30HduvUQGztUljarShAEyQ+1YFgTESmcTqfDhImTcPduPmYlzpAc2FeuXEHizOm4e/cuJkx8qUrXEZcX56zFYlgTETmAkJAQJCRMxY0b15Hw6isWD4kfO3oECa++jJs3byAhYSpCQkJk7ilZA8OaiMhBNA8LQ+KsOfD09MDCBfOwePEi0avEMzLSsXjxIixcOB+enp5InDUHzcPCrNzjSgiQeOqWfbtvS4q/KAoREf0tJCQEr72+CBs3foodX23Hgf3JaNykCdq0bouQ0FDUrVsPOp0Wer0BmZkXkJ6WhpQfj+P0qVPQaDTo138AYmOH2nHomyzBsCYicjA6nQ4jR45C3779sGdPEvbuScInn3xU7vZBQUEYPnwEukV3V9Sqb0Hg5UbFsnlYG0xmRM49bfNjal3V800lInXw8fHBwIExGDgwBrm5ucjISEd2djaMBgM0Wi0CAgIQEhJq+yuTkexsGtZdm3piX+odWx4SAKB1FdC1qafNj0vOY8SIEUhO3o+uXbti/foPSr3eu3dvZGdfhptbNQDArl274O7ubutukop5e3ujVauH7d2NKuL9rMWyaVgvHhZky8MRySYuLg4jR47ERx99XO42H3zwPpo3b27DXhE5OPVkrWRcDU4kQmRkJLy8vOzdDSInw/OsxWJYE8nk+eefR0REJJYvX27vrhCRk+FqcCIZrF69GnXq1EFOTg6GDRuGRo0aoWfPnvbuFpGi2XI1uMFgwNq1a3Dil19w61Yu/Pz80K//gEpvYHLz5k1MnvQSatasif8tesPivkrFsCaSQZ06dQAUrc4diJSUFIY1UaVst8DMZDKhhq8vps+YiYCAAJw5cwYL5s+Fv79/hQvz1q1dg+DgYNy6dUtCP6XjMDiRREajEdeuXQMA6PV67Nq1G02bNrNzr4jon9zc3DD0yWEIDAyEIAho3LgxwsJaIPX338vd59ixo7h9+zaioqJs2NOyMayJROjXrz9Gj34au3btQrNmzXH06FEMHjwEly5dQkFBAQYOjEF4eDgiI6PQrFkzDBjQ395dJlI+mdaX5efnIy8vr/hhMBgqPbRer8fZs2fw4IMNynw9784dfPD+exgzdqyET1A+HAYnEmHbtq2lntu06bPi///++3227A6Rk5BnGDw+rmSgDh4SW+FtP81mM1auXIGgoCB06NixzG0+/HA9uj7SDUFBdXAqNVVCH+XBsCYiIrsQIM8CsxUrV5W4CJFWqy13H7PZjDWrV+HSxUxMnz4TLi6lB5h///03nDqVitde+5/FfZMbw5qIiByau7s7PDw8Kt3ObDZj7ZrVOHv2DKbPSISHZ9lXtjxx4gSys7MxbtwYAPdWkuv1ejz37NNY9MabqFGjhqz9F4NhTURE9iHY9nKja9euwalTqZgxc1aFFznq06cvund/tPjjQ4cOYk/Sbvx32gz4eHtb3FspGNZEROT0rly5jG93fgOtVosX4uOKn4+MisLYseMwf95cNG3WDDExg+Dh4VGiUvfy9ISrqwb+/v726DoAQDAaDWa7HZ2IiFQnLy8PT48eheSCLjBJqBldYURktQN47/31oobBHRkrayIisg8bD4M7MoY1ERHZhVyrwdWAYU1ERPbBylo0XsGMiIhI4RjWRERECsdhcCIisg8Og4vGypqIiEjhWFkTEZFdcDW4eAxrIiKyDw6Di8ZhcCIiIoVjWBMRESkch8GJiETQGwtx7XYhanm7QFPGPZDJAhwGF41hTURUDoPJjJy8QuTkm5Cn//v5IF+GtTwY1mIxrK3s64Op2Lz3BOrW8sHAR1qgVeM69u4SEVXg74AuRJ6+9E0J/b0Y1GR7DGsruXrzDjo+vRRXc/KKn3vz42Q0CKqB5wd0RP+oMNQP9LVfB4moFKPJjNNZBhSWc+NgDx2g0zCs5SIIrKvF4k+dlYQ/t7xEUBc5d+kGpq34Bi2HvYHo+HexbOMB/JF10/YdJKJSXF0AN235EeDj4WrD3qiBIMNDHVhZW8GRk38g+/rtSrc7nnoBx1MvYNqKb9C2aT0M7NaCFTeRHQmCgAY1NTh31VBijrqIjzvrG/mpJ3ClYFhbwbbkX6u8z/3BPbpPW4zs1QaurnxzILIlFwHQuAgASo6Fe+gEaF0ZLGQfTAIrCK3rL2n/46kXMHHRVvR4cTXu5Jfx5z0RWYXZbMaf103IvWuGt5sAD93fr/l48O1SdoIg/aES/Omzgr6RzWVp53jqBfx3xTeytEVEFSsK6pz8Qvi4u6C+vwYNamrhoRPgInAI3Bo4Yy0ef/qsoFYNL3QMe0CWtjbvOYHCwkJZ2iKist0f1A/4uUIQBLi6CAippUGTQC2HwMmuGNZWMmZAR1nayblzF9fKWFVORPIoL6iLCIIADYPaSqQOgavn+8KwtpJe4U1RTSt9/Z5vdXfU9PWUoUdEdL/KgpqsjQPhYjGsraS6RzU82qGh5HbiB3XimweRFTCoFYALzERjWFvRwEdaSNq/y0MN8K/hUTL1hoiKMKjJ0TCsrUjqUHhUm2DoZBhKJ6K/MaiVg4Pg4jGsrUjqUPiC9/Zi8Uffy9gjInVjUCsMh8FFY1hbmdSh8FmrdzGwiWTAoCZHxrC2MjlWhTOwiaRhUCsVB8LFYlhbmVyrwhnYRJZhUCsYh8FFY1jbwLP9OsjSDgObqGoY1MrGulo8hrUNPNqhEQZ0DSvztYhWDfDOqwPh6iLuW8HAJhKHQU3OhOcF2ci66bFo2TAZH+38EemZ11E/0Beje7fFS09GQOPqCjedFmPmbYJJxHXAZ63eBQA8B5uoHAxqByF5KFs931OGtY24urrg5ZFd8fLIrjAYTdBqXEu8Pii6JQAwsIkkYlA7EqmD2er5vnIY3A7uD+oig6JbYvV/B3NInMhCDGpyVqysFYYVNpFlGNSOiJW1WAxrBWJgE1UNg9oxSZ2yVtN3mGGtUAxsInEY1I6MlbVYnLNWMM5hE1WMQU1qwcpa4VhhE5WNQe0EeOqWaKysHQArbKKSGNSkNgxrB8HAJrqHQU1qxGFwB8IhcVI7BrVzEQRB0vdPUNEwOMPawTCwSa0Y1M6Iq8HFYlg7IAY2qQ2D2klxgZlonLN2UJzDJrVgUBMxrB0aA5ucHYOa6B4Ogzs4DomTs2JQOz8uMBOPlbUTYIVNzoZBTVQSK2snwQqbnAWDWk24GlwshrUTYWCTo2NQqwxXg4vGYXAnwyFxclQMaqLyMaydEAObHA2DmqhiHAZ3UhwSJ0fBoFY3Na3oloJh7cQY2KR0DGqV45y1aAxrJ8fAJqViUBNXg4vHOWsV4Bw2KQ2DmqhqWFmrBCtsUgoGNRWTWlirCMNaRRjYZG8Mavon4a9/UvZXCw6DqwyHxMleGNRElmNlrUKssMnWGNRUJq4GF41hrVIMbLIVBjWVz3arwQ0GA9auXYMTv/yCW7dy4efnh379ByA6unupbXNycvD+e/+H3377Ffn5+QgICERs7FC0a99eQl+lYVirGAObrI1BTRWy4QIzk8mEGr6+mD5jJgICAnDmzBksmD8X/v7+aNXq4RLb3r2bjwbBwRgxchRq1KiBlJTjWPLWm1iw4DXUe+AB23T4PpyzVjnOYZO1MKhJSdzc3DD0yWEIDAyEIAho3LgxwsJaIPX330ttGxAQiH79+sPf3x8uLi5o16496tSpg9NnTtuh5/ewsiZW2CQ7BjWJI88weH5+folntVottFpthXvq9XqcPXsGERGRlR4lJycHFy5k4sEHH7S8qxIxrAkAA5vkw6AmsQRBkPSzUXTqVnzc2BLPDx4Si9jYoeXuZzabsXLlCgQFBaFDx44VHsNoMOCtNxejc3g4QkMbWtxXqRjWVIyBTVIxqMkeVqxcBXd39+KPK6qqzWYz1qxehUsXMzF9+ky4VDAFaDQY8MYbi1CtWjXEjYuTtc9VxbCmEhjYZCkGNVWdPMPg7u7u8PDwqHRrs9mMtWtW4+zZM5g+IxEenp7lbms0GLB48RswGo2Y8moCNJUMq1sbw5pKYWBTVTGoySI2vtzo2rVrcOpUKmbMnAUvL69ytzMajVj85hsoKLiLVxOmVjr/bQsMayoTA5vEYlCT5Wx3nvWVK5fx7c5voNVq8UL830PakVFRGDt2HObPm4umzZohJmYQTp86hR+OHYNWq8Nzzz5TvO3AmBjExAyS0F/LCUajwWyXI5ND+HzPCdGBDQAzx/RgYKsIg5oskZeXh6dHj0KK72gUCjqL23Ex69Hm5vt47/31oobBHRkra6oQK2wqD4OapOKNPMRjWFOlGNh0PwY1yYK3yBSNVzAjUXilMyrCoCayPVbWJBorbGJQk7xst8DM0TGsqUoY2OrFoCb5SbxFplk9P38Ma6oyBrb6MKjJGrjATDzOWZNFOIetHgxqIvtjZU0WY4Xt/BjUZFWcshaNYU2SMLCdF4OarI9pLRaHwUkyDok7HwY1kbKwsiZZsMJ2HgxqshlB4mpwFf1cMqxJNgxsx8egJlvianDxGNYkKwa242JQk81xylo0zlmT7DiH7XgY1ETKxsqarIIVtuNgUJP9sLQWi2FNVsPAVj4GNdmVAIkLzGTrieJxGJysikPiysWgJnIcrKzJ6lhhKw+DmpSAq8HFY2VNNsEKWzkY1ESOh5U12QwrbPtjUJOi8KIoorGyJptihW0/DGoix8XKmmyOFbbtMahJmXjqllgMa7ILBrbtMKhJqQRBkPSzqKafYw6Dk91wSNz6GNREzoFhTXbFwLYeBjWR8+AwONkdh8Tlx6Amh8DV4KIxrEkRGNjyYVCT4+ACM7E4DE6KwSFx6RjURM6JlTUpCitsyzGoydFwFFw8hjUpDgO76hjU5Jg4DC4Wh8FJkTgkLh6Dmsj5sbImxWKFXTkGNTk0joOLxrAmRWNgl49BTY6Pw+BiMaxJ8RjYpTGoyRmwsBaPc9bkEDiH/TcGNZH6sLImh8EKm0FNzobD4GIxrMmhqDmwGdTkdDgOLhqHwcnhqHFInEFNpG6srMkhqanCZlCT8+IwuFgMa3JYaghsBjU5M0a1eAxrcmjOHNgManJ6nLMWjXPW5PCccQ6bQU3kuL788kvcvHlT1jYZ1uQUnCmwGdSkHoIMD+VZtOgNNGrUGFFRXTFt2nR8++23uH37tqQ2GdbkNJwhsBnUpCpFw+BSHgq0b993OH36NKZMeQV6fQFmzkxEcHAIHnusp8Vtcs6anIojz2EzqImcR40avmjcuDEuXcpCVlY2Ll68iEIR70nlYViT03HEwGZQEzmP5557HgcOHIC/vx+6du2KYcOexPLly+Dt7W1xmwxrckqOFNgMalIrQRAk/awr9fdk79698Pb2xqOP9kBkZCTCwzvDw8NDUpucsyanVeYcdtpXwM+rgPQdpbaftXoXHmoXgXbt2iMiIgIRERHIz8+3ah8Z1KRuzrnALD09DR9+uB4BAbWxZs1qtGz5EHr27IX58+db3CYra3JqpSrs2q0A/+bA9d/L3P78pRuIn5SAha+MsnrfGNREzqtFixZo0KABGjZsiODgYGzYsAHHjh3D1KlTLWqPlTU5vRIVdvV6gKu2wu1XfH7I6qvEGdRE+Ks4lrIa3N6fQNkSE2fh0Ud7IDg4BImJiSgsLMTbb7+NtLQ0i9tkZU2qUKLCrmzjjJ2Y9a9vcWh3P3y2brHsfWFQEzm33NxcjB8/HpGREahZs6YsbTKsSTWKAvv5V99CuUvOGjwG6LwAUwG+3bkdcf+pi5UL/i1bHxjURM5v8eI3iv//2rVr8Pf3l9wmh8FJVQZFt8TLo7qWfzEFnde9/7pWA2o0wsebd8o2JM6gJiqpaDW4lIcS5efnY9KkyQgMDELDho0QGBiESZMm486dOxa3ybAm1YlqHYy2TeuWvtKZuRAw/rX6u9AE5JwH3P1kudIZg5qoLM65Gnzq1P/i7Nkz2LZtK06dSsWXX25DWloapk2bbnGbHAYnVenXrz9OnjyJvLw8eHv8hpxa3VB48QhQP/peNX1m673QhhnwaQD4NgQg7TxsBjVROZz0Hplff/01Dh48CD+/GgCA2rVr4/3330PnzuF4803L1sEwrElVtm3bWuLjz/ecwJh5QX9fOKXZk+Xua0lgM6iJ1MdsNsPFpeTvuSC4wGw2W9wmh8FJ1ax58w8GNVFlnHMYvGfPnnjqqdFISfkRV69exfHjKXjmmWfQq1cvi9tkWJPqWSOwGdRElRNk+KdE8+fPwwMP1EOvXr3QqFFjPPHEE6hbtw7mzZtrcZscBieCvNcSZ1ATiSX1NpfK/L3y8vLC22+/jeXLl+Pq1auoWbOm5PcAhjXRX+QIbAY1ERURBAG1atWSpS2GNdE/SAlsBjWRetWv/6Co3/fz589Z1D7Dmug+lgS2IAgY9FhnBjVRVQgSh8EV9Dv20UcbrNo+w5qoDFUJbFdXF7joqjOoiRTMYDBg7do1OPHLL7h1Kxd+fn7o138AoqO7l7l9Xl4eVq96Fykpx6HT6dCz1+MYPHhIue3PnJmIpKTdAICFCxciISFB1v5zNThROcSsEnd1dcHM+IGI7tgce478hs93HmRQE4lky9XgJpMJNXx9MX3GTLz/wYd4YfwErP/gffz8809lbr9u3Vrcvn0b76x4F7Nmz0XS7t3Yt++7cts/c+YMjEYjAGD58rer8mUQhZU1UQUqqrDvD+pZ72yGqdAMMyy70hmR6sg0DJ6fn1/iaa1WC6225K1w3dzcMPTJYcUfN27cGGFhLZD6++9o1erhEtsWFBTg4IH9mDN3Hjw9PeHp6YnHH38ce5KS0LXrI2V2JTIyAhERkQgNDUV+fj5GjBhZ5nYbNnxYxU/yHoY1USXKCuzyghqQdmlSIjWR62qj8XFjSzw/eEgsYmOHVrivXq/H2bNnEBERWeq1ixczYTQa0aBBcPFzDRoE44svNpfb3rp167B161acP38e3377LVq2bCH+ExGBYU0kwj8DGwLKDeoiDGwi21mxchXc3d2LP76/qr6f2WzGypUrEBQUhA4dO5Z6/e7du6hWzQ2urq7Fz3l4epaq4P+pWrVqiI2NBQDcvJkj+5w1w5pIpKLAzrxpQrcO5Qd1EQY2UcXkWgzu7u4ODw8PUfuYzWasWb0Kly5mYvr0mXApY02Km5sb9PoCmEym4sDOy7tT4g+Ciki5Ull5GNZEIpnNZrRv1QyN8wux92jFQV2EgU1UPqn3pK7qvmazGWvXrMbZs2cwfUYiPDw9y9yuTp26cHV1xflz5xASGgoAOHfuHOrXr29xX6XianAiEe6/4EldX1eInW2T437YRCTd2rVrcOpUKqZNnwkvL69yt6tWrRrCw7vg008/Rt6dO7h06SK++XoHors/asPelsTKmqgSZV2ZrL6M1xInUitBAFxsdE2UK1cu49ud30Cr1eKF+Lji5yOjojB27DjMnzcXTZs1Q0zMIADAs889j1WrViIubix0Oh169Xq83JXgtsCwJqpARZcQlfPmH0RqJEhcDl6VsK5VqzY2fvZ5ua9P/e+0Eh97eHhg0qR/WdSv1avXYMyY50s9/9JLk7BkyVsWtclhcKJyiLnWtzXvh03k7JzzbtbA22+/ja1bt5Z47l//+jd+/fVXi9tkZU1UhqrclIMVNhH906ZNn6Ffv/7w9/dHREQEpkyZgpSUFGzdusXiNhnWRPex5O5ZDGyiqhMknrul1Ev7NmzYEOvXf4ARI0YiKioSqamnsG3bVvj4+FjcJsOaVOPChQsYN24crly5Co3GFa+8MgUDBw4osU1cXDy+338Q7h5ecHUR8MlHH0DwDxHVPgObqGpcBMDsHDfdwsmTJ0t8XK1aNYwbNw4rV67EmjWrceHCBVy4cAEtWlh2ZTOGNamGRqPBggUL8NBDDyE7Oxtduz6Cxx7rAc+/zrU0m824ozdjYsI8PPF4L4vunsXAJlKniIhICIIAs7n0tRf69u0H4N5IwI0b1y1qn2FNqhEYGIjAwEAAQEBAAPz9/XDjxg14enoWD30bjGZ46gRJt7lkYBOJY8vV4NZ28+YNq7bP1eCkSj/++BNMpkLUq1evxBy1ViNgyeuJiIiIQGLiLJhMJova5ypxosoVTVlLeagFK2tSnevXbyAuLg5Lly4ptZjs9XkzERgYiIKCAsTFxWPdunUYM2aMRcdhhU1UMRdBgNkJF5hdvHgR8+bNw08//YRbt26XeO2XX362qE2GNalKQUEBRowYjsmTJ6FDhw5lrPoOAnDvQv7Dhj2JLVu2SDoeA5tIfcaOHQt3dw9MmjRJ9A1GKsOwJtUwm82Ij38BUVFRGDp0aJmnZ2VlZSEwMBCFhYXYseNrNG3aTPJxGdhEZZNaFyuzrgZ++ulnpKenQafTydYm56xJNQ4fPozNmzdj+/av0LFzJAY+3hXZ51Px+szJ+PHHnwAAzz8/BuHh4QgP7wKTyYS4uHGyHJtz2ESlOeucddOmTZGdnS1rm6ysSTU6d+6MGzeul6qoH+2yrHib7du/tNrxWWETqUPfvn0xbNgwPP/8GNSuXavEa0888YRFbTKsSTUsuTKZ3BjYRH9zplO3/mnNmjUAgDfeeKPE84IgMKyJKqKEoC7CwCa6x1lXg5848YvsbTKsyekpKaiLMLCJnLeytgaGNTk1JQZ1EQY2kfPo2bMXdu78BsDflx4tS3KyZQtHGdbktJQc1EUY2KRmznTq1vPPP1f8/y+8EC97+wxrckqOENRFGNikVpJPv1LQr/SQIUOK/3/48OGyt8+wJqfjSEFdhIFN5Nh27NghajuuBieCYwZ1EQY2qY3UBWZKqqxffTWh0m146hYRHDuoizCwSU0EOM84uDVO1/onhjU5BWcI6iIMbFILFyeqrK2N1wYnh+dMQV2E1xInon9iZU0OzRmDuggrbHJ2zjRnbW0Ma3JYzhzURRjY5Myc67fVujgMTg5JDUFdhEPiRMTKmhyOmoK6CCtsckaC1KuiOPnv/T8xrMmhqDGoizCwydlwNbh4DGtyGGoO6iIMbHImXGAmHuesySEwqP/GOWwi9WFlTYrHoC6NFTY5A1bW4jGsSdEY1OVjYJOjcwHTWiwOg5NiMagrxyFxInVgZU2KxKAWjxU2OSoOg4vHsCbFYVBXHQObHBHDWjyGNSkKg9pyDGxyNAxr8ThnTYrBoJaOc9hEzomVNSkCg1o+rLDJUQiCIOn33Kyi9wiGNdkdg1p+DGxyBFIvDQ4BMMvWG2XjMDjZFYPaejgkTuQ8WFmT3TCorY8VNimZ1PVlgHoqa4Y12QWD2nYY2KRULhKHwc0CUPlPtHNgWJPNMahtj4FNSiR1gZma7mfNOWuyKQa1/XAOm8hxsbImm2FQ2x8rbFISOVaDqwXDmmyCQa0cDGxSCjnmrNWCw+BkdQxq5eGQOJFjYWVNVsWgVi5W2GRvcpy6pRYMa7IaBrXyMbDJnu7NWUtZDS5fX5SOYU1WwaB2HAxsshcuMBOPc9YkOwa14+EcNpGysbImWTGoHRcrbLI1F0hcDS5bT5SPYU2yYVA7PgY22RKHwcVjWJMsGNTOg4FNtiL89U9KC2rBOWuSjEHtfDiHTaQsrKxJEga182KFTdbGK5iJx7AmizGonR8Dm6yKYS0ah8HJIgxq9eCQOJH9sbKmKmNQqw8rbLIGqavB1fS2w7CmKmFQqxcDm+QmCIKk9w81vfcwrEk0BjUxsElOLsK9h8VU9PbDsCZRGNRUhIFNjuqbr3fgu+++wx9/nMfDrVtjypSEcre98OefWLduLTIy0qHRaNGuXTs8/cyzqFatmg17/DcuMKNKMajpflx0RnIQZHhURQ0/P8QMGoTu3R+tdNslS95CnTp1sHr1WrzxxmKcP38en2/6rIpHlA/DmirEoKbyMLBJqqIFZlIeVdGxYyd06NAR1b29K9328uVsREZFQaPVwtvHB+3atcMff/xh4WcqHcOaysWgpsowsEkJ8vPzkZeXV/wwGAyS2+zbtx/27dsHfUEBbt64gaNHj6Jtu3Yy9NYynLOmMjGoSSzOYZOlXAQBLjKcuxUfN7bE04OHxCI2dqiUruHh1m2w4p3leOqpkSgsLET79h3QrVu0pDalYFhTKQxqqioGNllCrvOsV6xcBXd39+LntVqtpH7dvn0bc2bPwtChQ/HYYz1xt6AA69atxbKlSzD5X/+W1LalOAxOJTCoyVIcEqeqkmvO2t3dHR4eHsUPqWGdnZ0FvV6Px5/oDY1WCy8vL/To0QMpKSkyfNaWYVhTMQY1ScXAJiUzmUzQ6/UoNJlgLjRDr9fDWMb8dt06deHm5oadO7+ByWRCfn4+knbvRnBwsB16fQ+HwQkAg5rkwyFxEssFtq0YP/98EzZ9trH445EjhqF58zAkzpqN+fPmommzZoiJGQQ3d3e8mvAfbPhwPT75+CO4uLigSZOmGP/iizbsbUmC0Wgw2+3opAgMarKGz/ecEB3YADBzTA8Gtkrk5eXh6dGjUL//Irho3SvfoRyFhnz8sfVlvPf+enh4eMjYQ+XhMLjKMajJWjgkTiQfDoOrGIOarI1D4lQR3shDPIa1SjGoyVYY2FQe3iJTPIa1CjGoydYY2FQW3nVLPM5ZqwyDmuyFc9hElmNlrSIMarI3VthUgsRhcDVV1gxrlWBQk1IwsKmI8Nc/KfurBYfBVYBBTUrDIXGiqmFl7eQY1KRUrLBJ6gIzs4reyhjWToxBTUrHwFY3nrolHsPaSTGoyVEwsNWLYS0e56ydEIOaHA3nsIkqxsrayTCoyVGxwlYfARIvN6qi1eAMayfCoCZHx8BWF6m3yFTTLSM5DO4kGNTkLDgkTlQaK2snwKAmZ8MKWx24wEw8hrWDY1CTs2JgOz+GtXgMawfGoCZnx8B2bi6CABcJ71lmFb3fcc7aQTGoSS04h03EytohMahJbVhhOycOg4vHsHYwDGpSKwa282FYi8dhcAfCoCa145A4qRUrawfBoCa6hxW28+BFUcRjWDsABjVRSQxs5yAIEi83qqL3QYa1wjGoicrGwHYCEuesVXRpcM5ZKxmDmqhinMMmtWBlrVAMaiJxWGE7Lq4GF49hrUAMaqKqYWA7Jhfh3kPK/mrBYXCFYVATWYZD4uTMWFkrCIOaSBpW2I5F+OuflP3VgmGtEAxqInkwsB0H56zFY1grAIOaSF4MbMfAOWvxOGdtZwxqIuvgHDY5E1bWdsSgJrIuVtjKxmFw8RjWdsKgJrINBrZy8XKj4nEY3A4Y1ES2xSFxcnSsrG2MQU1kH6ywlUeAtMt7q+mdk2FtQwxqIvtiYCsLV4OLx7C2EQY1kTIwsJWDC8zE45y1DTCoiZSFc9jkaFhZWxmDmkiZWGHbH1eDi8ewtiIGNZGyMbDti3PW4nEY3EoY1ESOgUPi5AhYWVsBg5rIsbDCtg8BEheYydYT5WNYy4xBTeSYLAlsQQAmD2NgS8F3R3E4DC4jBjWRY6vqkHjiql04+Ms563bKibkIguSHWjCsZcKgJnIOVQ3sGe9+a+UeEXEYXDaXbxUyqImcRFWGxE+fv2KLLjklXhRFPIa1TLyqCYC3K2pXd2FQEzkBsYEtuAgoLCyEi8hKnP7GsBaPYS0Tz2ou8Kxm714QkZzEBPaArmEMagsxrMXjTxgRUQUGRbfEmmmDUU1burZp+IA/EkZ3s0OvSG1YWRMRVSKmW0u0alQHM1btxG/p2TCaCtEvsjkmDI1AoH91e3fPYblAWsWopmqTYU1EJEJoPX9smD3c3t1wKrw2uHhq+sOEiIjIIbGyJiIiu+ACM/EY1kREZBe865Z4DGsiIrILVtbicc6aiIhI4VhZExGRXXA1uHgMayIisgsB0m6RqZ6o5jA4ERGR4rGyJiIiu+BqcPEY1kREZBe2Xg3+zdc78N133+GPP87j4datMWVKQoXbJyXtxratW3H9+jV4e3vj6WeeRfv2HSzvsAQMayIRLly4gHHjxuHKlavQaFzxyitTMHDggBLbxMfH48CBg/D2vnet6A8+WI+QkGA79JbIMdh6gVkNPz/EDBqEE7/8gmvXr1W47e5d3+Krr7Zj0uTJaNAgGDk5OSgouGtxX6ViWBOJoNFosGDBAjz00EPIzs5G166P4LHHesDT07PEdq+//hp69eplp14SqVN+fn6Jj7VaLbRabantOnbsBAA4d+5chWFdaDLh008/xYsTJiA4OAQA4OvrK1+HLcCwJhIhMDAQgYGBAICAgAD4+/vhxo0bpcKaiMQTJM5ZFxXW8XFjSzw/eEgsYmOHWtzuxYsXkZNzExnp6Vj17kqYTCY83LoNnnpqNDw8PCzvsAQMa6Iq+vHHn2AyFaJevXqlXps2bTrmzJmDHj0ew/Tp0+Dq6mqHHhI5BrlO3VqxchXc3d2Lny+rqq6K27dvAwBOnPgFCxa+DgBY8tZivP/e/yH+hfGS2rYUT90iqoLr128gLi4OS5a8Veq1mTNn4tixo0hKSsK5c+ewbt0623eQSIXc3d3h4eFR/JAa1m5ubgCAAQNj4O3tDW9vbwwYGIPjx3+Qo7sWYVgTiVRQUIARI4Zj8uRJ6NixY6nXAwMDIQgC3NzcMGzYk0hJSbFDL4kcR9FqcCkPa6hTpw60Wp11GrcQw5pIBLPZjPj4FxAVFYUnn3yyzG2ysrIAAIWFhdix42s0bdrMll0kcjgugiD5URUmkwl6vR6FJhPMhWbo9XoYDYZS2+mqVUNkVBS2btmC27dv486dO9i6ZQva2em0LYBz1kSiHD58GJs3b0aLFmH46quvAADvvvsuVqxYiWeffRZt2rTG88+PwfXr11BYaEa7du0QFzfOzr0mUjZbn2f9+eebsOmzjcUfjxwxDM2bhyFx1mzMnzcXTZs1Q0zMIADA008/g7VrVuPF8fHQarVo2649Ro9+2vLOSiQYjQaz3Y5ORESqk5eXh6dHj8LIV9+Frpp75TuUQ1+Qjw9fG4f33l9vt1XatmLTyvpfH1/CvtQ7tjxksa5NPbF4WJBdjk1ERKXxftbi2TSs96XegcFkhtbVtl9hg8lstz8SiIiobAIkhrVsPVE+m89Za10FJE9rbNNjRs49bdPjERHZWm5uLtLT03D58mUYDQZotFrUrl0bISGh8Pb2tnf3SCIuMCMiclA5OTlIStqN7/buKT4boSyBgYHoFt0d0dHd4ePjY8MeVswFAlwk1MdS9nU0DGsiIgej1+uxceOn+Gr7dphMRjRp0hTdunVHSGgo6tatB51WC73BgMzMC0hPS0NKynF8/NEGbPz0U/Tu0wexsUOh09n/PGLOWYvHsCYiciDp6elYtvQtZGZmIiIyCgMHxKBBcNl3dwsICECbNm0xeEgszmVk4Istm7Ft6xYc/+EYJkychJCQEBv3viSGtXgMayJyCGkXriH5pwxcy7mDzi0fRIewB6BR2bXXf/v1VyxcOB9ubu5I+M9/q3Rv5QbBwZg8+d+I6BKBlStXIHHmdCQkTEXzsDAr9pjkwrAmIsX734ff4bX3v4PBaCp+ruED/pgy6hEMim6pitBOT0/HggXz4efnh5mJs1GrVi2L2mnfoSMaBIdgVuIMLFw4H4mz5titwnaReNctKfs6Gl5ulIgUbeu+XzF3bVKJoAaAs39ew9j5n6PjM8vw6a6fYDSZymnB8en1eixb+hbc3d0lBXWRWrVqYWbibLi5uWPZ0reg1+tl6mnV3BsGFyQ87NJtu2BYE5Givf3ZwQpfV0Nob9z4KTIzMxEX/4LkoC5Sq1YtxMXFIzMzExs3fipLm2Q9DGsiUrRfzl4StZ2zhnZOTg6+2r4dEZFRVZqjFqN9h47oEhGJHV9tR05OjqxtiyHI8FALh5iznvafV7Dzm69x4c8/sGvvfrRo+ZC9uySZ2WxG+sXruJBt+18QIkeSX1D6rkgVKQrtWat3YfLwKDzTt51Dz2knJe2GyWTEwAExlW7br18/ZGdnw8XFBV5eXli0aBFatWpV4T4DB8bgwP5k7NmThIEDKz+GnLgaXDyHCOvefQfghQmT0L93T3t3RRa7j57Bi//7Apeu3rJ3V4icVuaVXLy8ZDumrfwGn8wbiW5tQ+3dJYt8t3cPmjRpWu7pWf/0wQcfwNfXFwCwbds2jBs3DocPH65wn+DgEDRu0gR77RTWUhaJqSmsHWIYvHN4F9SpU9fe3ZDFzkOnMDhhPYOayEbuFhgx8OX38NOpTHt3pcpyc3ORlZWFNm3aitq+KKiL9hVEplmb1m2RlZWFW7f4vqRUDhHWzsJsNmPe/+2B2cy7khLZkhnAxDe22rsbVZaengYACAkVPyowZswYNGnSBHPmzMHq1atF7VPUftHxbEXaSnBB9B8jzsAhhsGdRYHBiJ/PXLR3N4hU6cyfV+3dhSq7fPkyAKBu3Xqi9ykK6A0bNmDGjBnYvHlzpfsUtZ+dnW1BLy3HOWvxWFnbkKCqtYtEyuKIv31Gw73FdTqttsr7jhgxAt9//z2uXbtW6bY6nbbE8Uh5GNY2VE2nwaMdGtm7G0Sq1OmhB+3dhSrT/BXSehEhevPmTVy69Pdpbl9++SX8/Pzg5+dX6b56vaHE8WzFRYaHWjjEMPgr/3oJSbt24vLlbAyLHQgvLy8cOvazvbtlkWnPdsfRX/9A7p0Ce3eFSDV0Glcs/Xd/e3ejymrXrg0AyMy8gICAgAq3zc3NxahRo5Cfnw8XFxfUrFkTmzZtEjWvm5l5AQAqPYbcOAwunkOE9f8WL7F3F2TTukld7F0Rh1eX78CRk3/gVh5Dm8iaGgTVwI4lz6FuLeXcx1mskJC/Fn6lpVW6Irx+/frYt2+fRcdJT0srcTxbkbpIjAvMyKoaPlATn7/2FEymQty8nW/v7hApWsiAhRbt17NTY/z32e5o1aiOzD2yHW9vbwQGBiIl5TgGD4m12nFSfjyOwMBAVK9e3WrHIGkY1nbk6uoCfx9Pe3eDSNF8vNyQc/uu6O37RDbDq091w0MNg6zYK9t5pFs0Pvn4I5zLyBB1YZSqyshIx+lTpzBs+AjZ264M77olnprm54nIAbVpKu6CSH0imyF59QvYMHu40wQ1AHTv/ihcXTX4Ykvlp2BZ4osvNkOj0SA6urtV2q9I0Zy1lIdaMKyJSNGmPdMdri7lv1U5a0gX8fHxQe8+fbA/+XscO3ZU1raPHT2CA/uT8UTvPvDxcbw5fTVhWBORorVr/gDemxmLAD+v4ucEQXD6kP6n2NihqFu3LlaueAdXrlyRpc0rV65g5cp3ULduPcTGDpWlzapiZS0e56yJSPH6RYWhe/tGSD13Gddy8tA+7AHUqO5u727ZjE6nw4SJk5A4czpmJc7AzMTZku5rfeXKFSTOnI67d+8i4T//hU6nk7G34gl//ZOyv1qwsiYih+DprkPbZvXwWKfGqgrqIiEhIUhImIobN64j4dVXLB4SP3b0CBJefRk3b95AQsJUhISEyNxT8YoWmEl5qAXDmojIQTQPC0PirDnw9PTAwgXzsHjxIpzLyBC1b0ZGOhYvXoSFC+fD09MTibPmoHlYmJV7THLhMDgRkQMJCQnBa68vwsaNn2LHV9txYH8yGjdpgjat2yIkNBR169aDTqeFXm9AZuYFpKelIeXH4zh96hQ0Gg369R+A2Nihdhv6/id9Qb6keWd9gXquU8GwJiJyMDqdDiNHjkLfvv2wZ08S9u5JwieffFTu9kFBQRg+fAS6RXdXxKpvjUYDX19fLJn1ouS2fH19odE4f5QJRqPBZjdXbjvzLAwmM7Sutp1oKDrm8VkNbXpcIiJbyc3NRUZGOrKzs2E0GKDRahEQEICQkFBFXplMr9fDaDRKbkej0ShilMDabPrnSNemntiXeseWhwQAaF0FdG3KK4URkfPy9vZGq1YP27sboul0OlWErFxsWlkTERFR1XE1OBERkcIxrImIiBSOYU1ERKRwDGsiIiKFY1gTEREpHMOaiIhI4RjWRERECsewJiIiUjiGNRERkcIxrImIiBSOYU1ERKRwDGsiIiKFY1gTEREpHMOaiIhI4RjWRERECsewJiIiUjiGNRERkcIxrImIiBSOYU1ERKRwDGsiIiKFY1gTEREpHMOaiIhI4RjWRERECsewJiIiUjiGNRERkcIxrImIiBTu/wEMjwDkS8pSygAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 480x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAc0AAAHXCAYAAADeCAHgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABJkElEQVR4nO3dd3hUVf4G8PdOn0lPJoVAIgHF/SGsqARC6BjK0sFFQXbdIrAUu2BHXV1AAesCIrIqCEoTlqIUZWlBQAQBC9hAIAJJJmXSJpl2fn+EDAkpzCST3MnM+3mePA+Z3Nz5chjy5pzznXslu90mQERERNekkLsAIiKi5oKhSURE5CaGJhERkZsYmkRERG5iaBIREbmJoUlEROQmhiYREZGbGJpERERuYmgSERG5iaFJ1ASeeOIJTJkyxa1jV65ciR49ejRyRfXz/vvL0K7djYiPb4njx4/LXQ5Rk1PJXQARNQ82mw2PP/44NmxYj9TUVLnLIZIFZ5pEBLvdDiHqvgx1ZmYmSktL0b59+3o9h81mq9f3EfkShiYRgI4dO+KVV15Fnz590aJFPO6444/Izc3DI488isTERNxyy604dOiQ6/jCwkI88MCDaNfuRrRrdyMeeuhhFBcXu76+f/9+dOuWivj4lhg//k8oLCyq8nynT5/BXXfdhTZt2qJDhw6YN28enE7nNetcuHAhhg4dVuWxjz9ej86dkwEAx44dw+23p6FVqwQkJbXBXXfdVeu5wsLCsWTJEqSkdEOLFvEoKiqqta7jx48jObkLAKB9+5tw882dAABFRUWYPn0GbrqpA9q2vR7/+Mc/YDabAQBnz55FWFg4VqxYgU6dbsH//V97V41Dhw7Fdde1RqdOt+D995e5apozZw7uuusuTJ8+A4mJibjppg74+OP1rq87nU4sXrwYnTsno2XLVrjlllvx+eefAwCEEK6vJSYmYsiQIfjhhx+uOaZEnmBoEl22YcN6rFjxAU6dOonffvsNaWlp6NOnN86cOYMxY/6Ihx9+2HXsE088gdOnT+PgwQM4cOAL/PTTj3jyyacAAHl5+Rg3bhwmTpyIc+fO4k9/Go81a9a4vrekpAQjRgxH7969cerUSWzduhUff7weK1asuGaNY8aMwcGDB5GRkeF6bPXq1a5wnDHjMQwaNAjnzp3FqVMn8cADD9R5vrVr12HDhvXIyDgPpVJZa10333wzDh48AAD4/vvvcPz4MQDAtGn3IS8vD/v3p+PEieOw2eyYMWNGlefYunUrdu/ehRMnjiMzMxMjR47C3/9+L06f/gUffrgSc+bMwe7de1zH79z5P6SmpuLMmTN45pmn8cADD6CwsBAAsGTJEixa9BbeeecdZGScx6ZNG5GQkAAAWLr0P/jggw+wevUqnD59GsOGDcNdd42F1Wq95rgSuYuhSXTZ3/9+L1q1aoWwsDD0798fkZGRGD58OJRKJUaPHo3vvz8Jq9UKp9OJNWvW4vnnn0NkZCSioqLw7LPPYtWqVXA6ndi+fRvi4lrg73//G1QqFf7whz+gV69erufZvn0HwsLCMXXqVGg0GiQkJGDy5MlYu3bdNWuMiYlBnz59sGbNWgBAdnY2du3ahbFjy0NTrVbh/PnzuHjxIrRaLbp3717n+R588AG0aNECWq3W47pMJhM2bdqE+fPnIzw8HEFBQXj66aewfv0GOBwO13GPP/44wsPDYTAYsGrVanTvnorRo0dBqVSiffv2GD9+PNauXes6/uabb3Z9fezY8tD7+edfAAD/+c+7ePLJJ3DLLZ0gSRISEhJw4403AgCWLl2Kp556Cm3btoVKpcLkyZNRWlqKr7766prjSuQuNgIRXRYTE+36s8Ggr/K5Xq+HEAIlJSWwWq2wWq1ITEx0fb1169YoKytDTk4OLl685Jr9VEhISEBZWSkA4Ny5czh58mSV73c6BVq2bOlWnWPHjsW8efPwyCMPY926dejatYvr+RYsWIiXX34JvXv3QXh4OCZNmohJkybVeq5WrVq5/uxpXWfPnoPT6cTNN/++yuMKhQKZmZmVnuPKWJw7dw47dnxW5TkcDie6devm+jw2Nsb1Z0mSoNfrUFRUPtM8f/482rZtW2M9586dw6RJ/4BSeWUuYLXacOHChZr/8kT1wNAk8pDRaIRGo8G5c+cQE1P+A/7cuXPQarWIiopCixZxOH/+fJXvycjIQHS0EQDQsmVLdOrUCTt3fl6v5x8yZDAefvhhfP31MaxatRoTJtzr+lqbNkl4++23IYTAwYMHMWLESCQnd8Ett3Sq8VwKxZWA8bSuVq1aQqFQ4NSpUzAYDNW+fvbs2cvPIVV5jqFDh+K999516zmulpCQgNOnT6NLly7VvtayZUu89NIcpKWl1evcRO7g8iyRhxQKBcaM+SNeeOFF5ObmITc3F//85wu46667oFAoMGDAQFy8eBHvv78Mdrsd27dvx969e13fP2jQQGRlZeGdd5aitLQUDocDP/30E/bt2+fW8+v1egwfPhwvvvgifvjhB4wcOdL1tY8++ghZWVmQJAlhYWFQKBRVZl518bSu2NhYDBkyBDNmzEBOTg6A8g7bzZs31/ocY8fehb1792Ljxo2w2Wyw2Ww4ceIEjhw56laNf/vbX/HSSy/jxIkTEELg/PnzrmafiRMnYNas2fjpp58AAAUFBfjkk09c+6FE3sDQJKqHl156CYmJiejatSu6dk1BmzZtMHv2LABAZGQEPvxwJRYvXozExOuwfPlyjBkzxvW9wcHB2LhxI/bs2YOOHX+PpKQk3HvvBGRmZrn9/OPGjcXOnTsxZMgQhISEuB7fvXs3unfvgfj4lhg37m68+OIL+P3vf1/Hma6oT11vvbUIYWFh6NOnL1q1SsCgQX/AsWO1X/QgPj4e69d/jPfeex/t2t2I66+/AdOnz3A72CZPnox77/07/vrXv6Fly1YYMWIkzp8vb4qaNGkS7r77bvzpT39Gq1YJ6NKlq1v7xESekOx2W91vziIiIiIAnGkSERG5jaFJRETkJoYmERGRmxiaREREbmJoEhERuYmhSURE5Ca/Ds2Ky55d65ZHRERE7vDr0LRYLPjrX/4Mi8XS4HPl5mR7oSL/xfGpG8endhybunF86taQ8RFCoLDQ7NZt+Sr4dWh6kyeDGog4PnXj+NSOY1M3jk/d6js+Qgjk5+egpLgIdrv7N0hnaBIRUUCpCExrWRnCI6Kg0Wjd/l6GJhERBYyrA1Or1Xn0/bw1GBERBQxJkqDV6GAwBHscmABnmkREFACEECgtLW8KNQTVLzABhiYREfm5iiVZc34eHA5Hg87F5dk6OJwCX54uwaFfLDDlW2AMN6FrWz26tDFAWelu9ERE5Juu3sNUKpUNOh9DsxYHfi7Gvz/LQWaBHUB5W7NCUYDNxwoQG6rC/f2j0O36IJmrJCKi2jS06acmXJ6twa6TRZi5PguZBXa0jdHg/rQoPN5fj/vTotAmWoPMAjtmrs/CrpNFjVbDkCFD8P3331d7fMqUKdi2bZtb5zh79ix69+7j5cr8z5w5c7BkyRK5yyAiLxNCQAjhtcAEONOsJqvAjpc/yYYQAtNuj8Ko20IhSRJM2aUwRodhxK2h2HCkAAt35uDlT7LRoZUO0SEcxqs5HI4GL4P4Qw1E1PSEEHA6HVAqVYiIMEKSvLedxpnmVbYcL4DNITCoYwhGdw6rNtiSJGF05zAM7BgCm0Ngy7FCj5/jtddeR0pKN3Trloo1a9YAKF/+ffDBh9C5czLGjh0Li6W0znMIITB9+gz8858vAADuvPMu9OrVGykp3VznrGzlypX485/vwbBhw9GhQwd8+OGHmDNnDlJTUzFy5CiUlZUBAGbPno0+ffoiJaUbnnzySdf3d+zYEXPmzEGPHj3Qp09fXLp0qdpzzJkzB5MnT0b//gPw2GOP4/TpMxg1ajR69+6DoUOH4ezZswCAn3/+GUOHDkX37t3Rp09fmM1mWCwWTJo0Campqejbtx9OnDgBp9OJm2/uhOLiYgDll0Xs0KED7HZ7receMmQInnjiCfTu3QerVq3C55/vRFpaf/To0RMTJ06C1WoFALz33vu45ZZb0b//APz4408e/fsRkXzsdieW78/D0Nd+RdcXfsYfFhag30unMWdLFvJL7K4l2bzcHAghvBqYAEOzmu3flC+5ju4cWudxd1z++rZvPAvNI0eOYsOGDdi9exc++eQTzJo1GxcvXsSmTZuRlZWJw4e/xMyZM3Hs2LFazyGEwCOPPIrQ0FA899yzAIDFixdj79492Lnzc8yf/4orBCv74YcfsHr1Kmzbtg3Tp8/A7373f/jiiy8QGRmJHTt2AAAmT56C3bt34cCBL3D+fAYOHjzo+v74+JZIT09H//5pWL58eY21nT59Bp98sgWvvDIfjz76KF5//TXs2bMbM2ZMx8yZ5bVOnDgJDz/8CPbv348tWzbDYDDgnXeWIjg4BF988QXmzn0ZU6ZMgUKhQP/+/bFt23YAwI4dO9C3bz+oVKpazw0AKpUae/bsxqBBg/Dmm29i8+ZNSE/fh9atr8OyZctw8eJFvPHGG9i1639Yv/5jfP311x79GxKRPH68WIp+L5/Bq9tMyMi1ocwmYHUAucUOrD5kRtrLv2LlvguwlpUhJLT6pMcbuK5YicMpYCq0w6BRoG1M3ZdVahujhUGjgKnQDodTuN1Ne+jQQQwfPhw6nQ46nQ69e/fG0aNHcfDgAYwePRqSJOGmm27CTTfdVOs5XnzxRfTt2w/PPjvT9diiRQuxdetWAEBGRgYyMjKgUlX95+3duxcMBgMMBgPUajUGD/4DAKB9+/Y4d+4cAGDPnj148803UVZWiuxsE9LS0pCSkgIAGDZsKACgU6dO+PTTrTXWNnjwYGg0GhQVFeHAgQMYP348gPKgNxiCUFBQgLy8PNx+ez8AQHBwMADg4MEDePDBBwEAycnJsFhKYTabMXLkCLz99hLcccdobNy4CePHj6/13BVGjRoJADh8+DC+++479O/fHwBQVmbFgAEDcOTIEfTq1Qvh4eGXa/5DrWNNRL7hfI4V97yTgVKbgF4jYWCHEIxPDUdZUS4OnNdg9ZdmZBc68MpnJQgxRGJEnHf2MK8mS2harVa8/vqr+C0jAxqNBqGhYZg4cRLiWrSoclxWVhbuv28aEhMTXY89On0G4uLiGqUuCeXLrzaHuGYQOpwCNoeAQpLgrXef1PRb0TvvvINly5YBAD777DMAwG23dcaXX36JoqIiBAcHY+/evTh48BB27tx5OYj7oKysrFpoVr6+okKhgFardf3Z4XCitLQUTz75JHbv3oW4uDg8/fQzsFrLKn2/BgCgVCrhdNb8XieDQQ+gfLk5Ojoa6enpVb5eUFDg0ZikpqZi2rT7kJubi8OHD2PJkrdRUlJS47kr6PVXahg4cAAWLVpU5etbtmxplN9AiajxPPVxJkptApFBSqy9LxFRweU/30wKBSb2icI9KUG4993f8O1FJ+Zty8ewWyOgUHh/MVW25dm0tP54/Y1/Y978V5GcnIzFi9+q8Ti9Xod5819xfTRWYAKAQiGhbbQGNofA4TN1307s8GkLbA6BtjEaj34Ap6R0w+bNm1FWVoa8vHzs3bsXt912G1JSumHDhg0QQuDkyZP47rvvAAATJ05Eeno60tPTXWEwZMhg/PWvf8H48X9CWVkZCgsLERkZCZ1OhxMnTuDbb7+t19+/tLQUkiQhMjISZrMZn3zySb3OAwChoaGIiYlxzX4dDge+//57hIaGIiIiAv/73y4AQFFREWw2G1JSumHt2nUAgCNHjsBg0CMsLMy1RDtjxmPo06cPVCpVree+WpcuXbBv3z7XLLqgoAC//vorbrvtNuzduxdmsxlFRUXYutW9bmQikkdOkR3fZZT3ebx2dwtXYAJXOmS1Oh2WTrgOWrWEojInNtej38QdsoSmRqPBrbfe5gqbG9q1Q3Z2lhylVDPslhAAwEcH82Fz1HzzaptD4KND+eXHdwrx6Py33noLRo4cid69+2Dw4MF46qknERcXh+HDh8FojEZyche88MIL6NSpU53nGT9+PPr3T8O9905AWloaioqK0KVLV8yf/8o1v7c24eHhGDduHLp06YqxY8ciOTm5XuepsHTpUrz99hJ0794d3bqlYs+ePQCAJUvexvz585Camorhw0egpKQEEydOgNlsRmpqKqZPn4GFCxe6zjNy5AisW7cOI0eOvOa5KzMajXjjjTfx5z/fg9TUVAwePBjnz59HixYt8MADD6Bv334YNWp0vceLiJrG2sNmOAUQF6bCzYl61+PlXbJOFBaaAQA6jQrdL79/fsNRz1a13CXZ7baak6EJ/fvNNxAcHIy//f3eKo9nZWXhwQfuQ+vWreF0OpGc3AWjR98BRS1vI7DZbLDZrtwXzWKxYMrkSXh/2QcwGAxu1WKxOnHvfzKQWWBH9xuC8OCAKEQFq2DKzoQxOhY5RXa8sSMH+38qRmyoCv+5txX0GvZTVYwP1YzjUzuOTd04PsCszZlY+2UBOl2nw/sTEgDUfuGCxbtysPh/uWgbo8HH91/n9VpkbwRav/5jXLp0Cc8+93y1r0VERGDx2+8gLCwMRYWFeO21V7F5y2aMGDGyxnNt2LAe69ZWf7tFjikLJXp9Dd9Rsxm3q/DCpzbs+6EQX/xYiFsSVIg0COSW/Iqj5+1wCiBcL2HG7SoUm7NR7PaZ/ZfDYYcpO1PuMnwWx6d2HJu6cXwAhaMUAkBuoRWm7EzXDBMon/MVFphRiPLZ5rlMCwQAFTwbN3d/MZE1NDdt2ogvDx3CzGefczWlVKZWqxEWFgYACA4JQd9+/ZCevq/W0Bw1ajSGDh3m+rxiphlljHF7pgkAxmhgcZwdy9Lz8L/vi3A0w3n5MnoKaNRK9Pu/YNzTIwKxobL/zuEz+Ntw3Tg+tePY1I3jA4zqWorVR87jfJ4TDk0kglWlKCosQHhEFAoLzFXG54szpyEB6P67MBijjV6vRbaf+ls2b8L+9HTMfPY5BAXVfA1Xs9mMoKAgqFQq2Gw2fHnoEJJaJ9V6TrVaDbVa7ZX6YkNVeGxwNKb0i8Sxs6W4ZMpDnDECna7TIUTHq8wQETWV37XQISFSjfO5NkxffQnv3xsPjUYLtVrjmmECwJufmZBf4oBaKWFC78hGqUWW0MzJycHy5csQGxuLfz7/HIDywJs95yWsXvURIiIjMWDAQJw6dRJrVq+6/JYIBzp06IjRd/yxSWsN0SnR88YgmCKLYIzmBdqJiOTw6CAjHv7oIr45X4oxCzPw8CAjerYrfxvc6ewyvLbNhH0/lgAAxiSHwtBIvSayhGZUVBTWrP24xq/dNXac689du6aga9eUpiqLiIh8kBACN8eVYmpPNd7aZ8PpbCvu/+AC9BoJEgQs1gJUdLQO6hiCx4bENFot3JQjIiKfVblL9p5ecejeHnh9hwlHfi2FxSogUP7eyRtiNZjQOxIDO3r2NkBPMTSJiMhnFZjzq7ytpH1LYMnfWqHE6sSv2VZkm3Jw8w2xCDc0TZwxNImIyGcFBQdDp9dXux+mQaNA+5Y6mDSqJgtMgHc5ISIiHyOEQFFhAYTTCZVK7bUbSHsDQ5OIiHxGxR5mcXEhbHa73OVUw9AkIiKfcPWl8SrurORLGJpERCS72q4l62vYCERERLKTJAkatRYGQ7DPBibA0CQiIhkJIWC1lkGr1SEouHHfY+kNXJ4lIiJZVCzJ5ufnwuFwyF2OWxiaRETU5KrsYYZHQlnLfZJ9DUOTiIiaVHNp+qkJQ5OIiJqUEALC6Wx2gQmwEYiIiJqIEAJOpxNKpRIRkdGQJEnukjzGmSYRETU6V9NPnglCiGYZmABDk4iIGlnlPczgkLBmG5gAQ5OIiBpRc276qQlDk4iIGo3NZoXNavWLwATYCERERI1ACAEA0Gi0MEbHQqFoHu/DvBbONImIyKsqlmSLigoAwG8CE2BoEhGRF1Xew9RotHKX43UMTSIi8gp/a/qpCUOTiIi8oqS4yK8DE2AjEBEReYkhKBgarRZqtUbuUhoNZ5pERFRvQgiY83Nhs1ohSZJfBybA0CQionqq2MMsLbXAKZxyl9MkGJpEROSxQGj6qQlDk4iIPFZgzgu4wATYCERERPVgCAqGTm8IqMAEONMkIiI3CSFQVFQAIQTUak3ABSbA0CQiIjdU7GEWFxXCbrPJXY5sGJpERFSnq5t+1Br/fltJXRiaRERUq0Dtkq0NG4GIiKhOapUGBkNwwAcmwNAkIqIaCCFgs1qh0WoRHBIqdzk+g8uzRERURcWSbH5+DpxOh9zl+BSGJhERuVTewwwLj/SrG0h7A0OTiIgAsOnHHQxNIiICADidTjgdTgZmHdgIREQU4IQQEMIJpVKJyKhoSJIkd0k+izNNIqIAVrEkm5ebAyEEA/MaGJpERAGq8h5mcEgoA9MNDE0iogDEpp/6YWgSEQUgq7UMVquVgekhNgIREQUQIQQAQKvVIdoYC4WS78P0BGeaREQBwnV7r+JCAGBg1gNDk4goAFTew1SrA/fWXg3F0CQi8nNs+vEehiYRkZ8rLi5kYHoJG4GIiPxcUFAItBod1BouyzYUZ5pERH5ICAFzfi5sNhskSWJgeglDk4jIz1TsYZaWWng/TC9jaBIR+RE2/TQuhiYRkR8xm/MYmI2IjUBERH7EYAiCXm9gYDYSzjSJiJo5IQSKiwshhIBGo2VgNiKGJhFRM1axh1lUWAC73SZ3OX6PoUlE1Exd3fTDy+M1PoYmEVEzxC5ZebARiIiomVIpVTBEBDMwmxBDk4ioGRFCwGazQqPRIiQ0XO5yAg6XZ4mImomKJdn8vBw4nU65ywlIDE0iomag8h5mWHgkFAr++JaDLMuzVqsVr7/+Kn7LyIBGo0FoaBgmTpyEuBYtqh175MhX+GD5MjidTiQmXoep0+6DwWCQoWoiInmw6cd3yParSlpaf7z+xr8xb/6rSE5OxuLFb1U7ptRiweK3FmHGY4/jzX8vREREBD5et1aGaomI5ON0OuGwOxiYPkCW0NRoNLj11tsgSRIA4IZ27ZCdnVXtuK+PfY3WrZPQsmUrAMDAgYOwf396ree12WwoKSlxfVgslsb5CxARNQEhBIQQUCqViDLGMDB9gE90z376ySfo3Dm52uMmkwnR0dGuz6NjYpCXlw+HwwGlUlnt+A0b1mPd2jXVHs8xZaFEr29QjQ6HHabszAadw59xfOrG8akdx6ZmQojLzT4C2VmXXJMMqspbrx9jdKxbx8kemuvXf4xLly7h2eeeb/C5Ro0ajaFDh7k+t1gsmDJ5EqKMMQ3eBzVlZ7o9qIGI41M3jk/tODbVVd7DVCgUiI6Jk7skn9XUrx9Z2682bdqILw8dwlNPPwOtVlvt60ajEdnZ2a7Ps7OyEBERXuMsEwDUajUMBoPrQ9/A2SURUVO7uulHktgl60tk+9fYsnkT9qen45mZzyIoKKjGYzp1ugVnzpzGb79lAAC2b9+G1O49mrJMIqImZS0rhbXMyqYfHyXL8mxOTg6WL1+G2NhY/PP55wCUzxJnz3kJq1d9hIjISAwYMBB6vR6TJ0/FvLkvw+FwIiExAfdNu1+OkomIGpUQApIkQavTwxitqXVFjeQlS2hGRUVhzdqPa/zaXWPHVfm8c3IyOidXbxIiIvIXFUuyGo0WQUEhDEwfxsVyIiIZVd7DVKnUcpdD18DQJCKSCa/00/wwNImIZFJUVMDAbGZkf58mEVGgCgoKgVarg0ZT/S135Js40yQiakJCCJjNebDbbVAoFAzMZoahSUTURCr2MEstJXA4HHKXQ/XA0CQiagJs+vEPDE0ioiZgzs9lYPoBNgIRETUBvSEIekMQA7OZ40yTiKiRCCFQUlwEIQS0Wh0D0w8wNImIGkHFHmZhoRl2u13ucshLGJpERF52ddOPWs3L4/kLhiYRkRexS9a/MTSJiLxMqVAyMP0Uu2eJiLxACAG73Qa1WoPQsAi5y6FGwpkmEVEDVSzJ5uWa4HQ65S6HGhFDk4ioASrvYYaFR0Kh4I9Vf8Z/XSKiemLTT+BhaBIR1ZPT4YDDbmdgBhA2AhEReUgIASEElCoVooyxkCRJ7pKoiXCmSUTkgYolWXN+LoQQDMwAw9AkInJT5T1MQ1AwAzMAMTSJiNzAph8CGJpERG4pKytlYBIbgYiI6lKxb6nT6aGOjoVSyR+bgYwzTSKiWpQvyeaipKQYABiYxNAkIqrJlT3MUiiVSrnLIR/B0CQiugqbfqg2DE0ioqsUFRYwMKlGXKAnIrpKUHAItDodNBqt3KWQj+FMk4gI5UuyBeY8OOx2KBQKBibViKFJRAGvYg/TYimBw+GQuxzyYQxNIgpoVzf9aLScYVLtGJpEFLCEEDDn57Lph9zGRiAiCliSJEGnN0BvCGJgkls40ySigCOEQElJMYQQ0On0DExyG2eaRBRQKu9hajQaqFRquUuiZoQzTSIKGFc3/TAwyVMMTSIKCLw0HnkDQ5OIAoZCUjAwqUG4p0lEfk0IAbvdDrVajbDwSLnLoWaOM00i8lsVS7L5eSYIp1PucsgPMDSJyC9V3sMMDYuApOCPO2o4voqIyO+w6YcaC0OTiPyOw2GH3WZnYJLXsRGIiPyGEAIAoFKpYYyOhSRJMldE/oYzTSLyCxVLsub8XABgYFKjYGgSUbNXeQ9TbwiSuxzyYwxNImrW2PRDTYmhSUTNWmmphYFJTYaNQETULAkhyu+HqdNDrdZApeKPM2p8nGkSUbNTviSbC4ulBJIkMTCpyTA0iahZubKHWQoFr/JDTYyvOCJqNtj0Q3JjaBJRs1FYaGZgkqy4EUBEzUZQUAh0Oj00Gq3cpVCA8nimuXnzZuTn5zdCKURE1QkhUFCQD4fDAaVSycAkWXkcmvPnv4IbbmiHXr1645lnZmLHjh0oKipqjNqIKMBV7GFaSorhcNjlLofI89Dcs2c3fvzxRzz22AxYrWV47rnnkZTUBgMGDGyM+ogoQF3d9MMZJvmCeu1pRkSEo127drh48RIuXcrEhQsX4ORd0YnIA6VWJw6dtuDcJStatyhGt7Z6qFTlv8cLIWDOz2XTD/kcj0Pz3nsnYP/+/YiKikTv3r0xbtxYLFjwb4SGhjZGfUTkZzLNNszZko0vfi6B1S4gAEi4AJ1aQq8bg/Dk0GhEBKmg0+mhNwQxMMmneByau3btQmhoKNLS+qNnz55ITe0Gg8HQGLURkZ/5/rdSTHg3AyXW8vteBmkVMKgFiqyAxSqw49siHPi5GB9MSkDraP5cId/j8Z7m6dO/YMWKDxAbG4OlS99Bx46/x8CBgzB79uzGqI+I/ERRqQMT3/sNJVaBcIMSL42Jw/5n2uKjv4dg/9NtMHNENII0EgpKBf66NANWO7d8yPfUa0+zQ4cOaN26Na6//nokJSVh5cqVOHz4MJ566im3z/Huu//Bka8OIzs7G3PnzkfrpKRqx3z33beYPWsW4uPjXY/NmjUbGi0bAoiamwU7c1Bc5kSIToENDyQiIujKjx9JktCvrR03jtNiwodlyC9xYnl6Hib0iZKxYqLqPA7N55//J9LT03H8+HHccMP16NmzJxYuXIgePXp6dJ6UlBSMGDESz858us7j4uPjMW/+K56WSUQ+ZuuJQgDA3d3CqwRm5S7ZGxKiMewWM9Z+WYC1XxUwNMnneByaBQUFmDZtGnr27AGj0VjvJ27f/qZ6f29tbDYbbDab63OLxeL15yAiz5VanTCXOKGQgL/1iKj2dQmSq0t2aj8V1n5ZgOwCvi+TfI/Hofnqq1dmfTk5OYiKatzfBDMzL+Hxx6ZDoVCgT99+GDhwUK3HbtiwHuvWrqn2eI4pCyV6fYPqcDjsMGVnNugc/ozjU7dAH5+8EicEAIUEFJmzUYTyGSYAOJ0O2O0SCgvMKITZ9fY1p0BAj1mFQH/tXIu3xscYHevWcR6HpsViwZNPPoVVq1ahrKwMWq0WY8eOxaxZ/0JQUJDHhdYlKakNFi9eAkNQEHJycjBn9r8QEhKC1NTuNR4/atRoDB06rEqtUyZPQpQxpsEdvqbsTLcHNRBxfOoW6OMT6XRCKRXB7gRy7WG4IU6L/Pwc2G02KBTKKmNz4KdiAEXQqKSAHrMKgf7auZamHh+Pu2efeupp/PzzT9i0aSN++OEUNm/ehF9++QXPPDPT68UZDAYYLgdxVFQUuvfoiVMnT9Z6vFqtLv+eyx/6Bs4uicg7FAoF2sWVN/C98VmOaw8zNCwCkiRVOXbRrlwAwM0JfH8m+R6PQ3Pr1q1YvvwDdOnSBTExMUhOTsayZe/j008/9XpxeXl5rqUai8WCo0e+qrHLloh838Q+5XuZ+38qweovi2q80s+inTn45nwpAGDq7WwCIt/j8fKsEAIKRdXfDCVJ4dqfcNeStxfj6NEjyM/Px6xZL0Kn0+PfCxZi8VuL0LlzMjonJ+PQwQPYsWM7lEolHA4HUrqlom/ffp6WTEQ+4Pb2IehxvRn7frZgcboN23/IxJ1dwhCltuPCqTx8/JUZ53PLG/mG3ByCW67jShH5Hslut3mUdg888CB+/fVXPP/880hMTMDZs+fw4osvIjExEW+++UZj1VkvJSUl+Otf/oz3l33APc1GxvGpW6CPT8Uv1UIIPL4mE59/V4SKHzzll9ErJ0nAyFtD8dzIwB2rqwX6a+damnp8PJ5pzp49C48//jgGDRoEm80GjUaDP/7xDsya9a/GqI+ImrmK92FKkgLh4ZGYN7YFMnKteOOzHHx91oJSqwMGjQpd2+pxX1oUYsPUcpdMVCuPQzM4OBgLFy7EggULYDKZYDQaq23kExEB1W/vVaFVpAbz7moBgDMpal7qdRk9oPyyV9HR0d6shYj8yNWBybuVkD9wKzQTE69zazZ59uyvDa2HiPyExVLCwCS/41Zofvjhysaug4j8hBACkiRBrzdAo9FApeIeJfkPt0Lzueeex86dnwMAXnrpJTzxxBONWhQRNU9CCJjzc6HTG6DT6RmY5HfcurjBTz/9BLu9/OLJCxYsbNSCiKh5qtjDLCsrZXMg+S23Zpo9e/ZAjx490bZtW1gsFowf/6caj1u5coVXiyOi5oFNPxQo3ArNd999Fxs3bsTZs2exY8cOdOzYobHrIqJmpLDAzMCkgOBWaGq1Wtx5550AgPx8M/c0iaiKoOAQ6HR6aLRauUshalQeX7CdV/4hIqB8SbawIB9OpwNKpZKBSQHB49AkIqrYwywpKXY1CRIFAoYmEXnk6qYfjYYzTAocDE0ichu7ZCnQeRya77yztMbHH3zwoYbWQkQ+TpIk6LR6BiYFLI9Dc+HChdi4cWOVxx555FF89913XiuKiHyLEAKllhIAgN4QxMCkgOXxXU7WrVuL4cNHICoqCj169MBjjz2Go0ePYuPG/zZCeUQkt8pLsmq1BkpVvW+ORNTsefzqv/766/HBB8sxfvyf0KtXT5w69QM2bdqIsLCwxqiPiGR09R4mA5MCnVv/A7799tsqn2u1WvzjH//A4sWLsXTpO8jIyEBGRgY6dOCVgoj8BZt+iKpzKzR79OgJSZIghKj2tWHDhgMobxDIy8v1bnVEJB8hAAEGJlElboVmfn5eY9dBRD5CCAGHwwGVSoXwiCjesYSoEr5Pk4hcKpZk83JNrptJE9EVHu/qX7hwAbNmzcKxY8dQWFhU5WsnThz3WmFE1LSu3sNkYBJV53FoTpo0CXq9AQ899BAMBkNj1ERETYxNP0Tu8Tg0jx07jtOnf4FGo2mMeohIBna7DTarjYFJdA0e72n+7ne/Q2ZmZmPUQkRNTAgBIQTUag2M0bEMTKJr8HimOWzYMIwbNw4TJkxETEx0la8NHjzYa4URUeOqWJJVKJQIC4uAQsG+QKJr8Tg0ly4tv2D7K6+8UuVxSZIYmkTNxNV7mETkHo9D85tvTjRGHUTURNj0Q1R/XI8hCjAWSzEDk6ie3JppDhw4CNu3bwNw5ZJ6Ndm3b6/3KiOiRqHXB0Gj1kKlVstdClGz41ZoTphwr+vPU6dOabRiiKhxCCFgzs913QuTgUlUP26F5pgxY1x/vvvuuxutGCLyvsp7mHpDkNzlEDVrboXmp59+6tbJ2D1L5FvY9EPkXW6F5uOPP3HNY/iWEyLfU1iQz8Ak8iK3QpNvMyFqngxBIdDpDNBotXKXQuQX+JYTIj8jhEBhoRlOpxMqlYqBSeRFDE0iP1Kxh1lSXAS73SZ3OUR+h6FJ5CeubvrRaDjDJPI2hiaRH2CXLFHT8Pjas0TkeyRJglajg8EQzMAkakRuhWZdl86rjJfRI2paQgiUlZVCp9PDEBQsdzlEfs+t0OSl84h8T+UlWXV0LJRKLhwRNTa3/pfx0nlEvuXqPUwGJlHTqFcj0IoVKzB8+AikpqYCANLT07F+/QavFkZENWPTD5F8PA7NefPmY9GiRbjjjjuQkZEBAIiLi8Obb77p9eKIqDohBIQQDEwiGXgcmsuXL8fatWvxl7/cA6C8OahNmzY4c+aMt2sjokqEEHA47FAoFIiIMDIwiWTgcWiWlJQgLi4OAFwdtTabDVpeqouo0VQsyeblmiCEcKubnYi8z+PQTE7ujKVLl1Z57IMPVqBr165eK4qIrqi8hxkSGs7AJJKRxy13c+a8hOHDh2Plyg9RXFyM/v0HICsrCxs3/rcRyiMKbGz6IfItHodmUlJrfPnlIWzfvgPnzp1Dy5YtMWjQQAQF8Y7wRN5mt9lgs1oZmEQ+ol5v7tLr9Rg5coS3ayGiy4QQAAC1RgNjdBwUCl4mmsgXuBWa06ZNc+tkCxcubFAxRHRlSVapVCE0NJyBSeRD3PrfGBoa6vqQJAXWrl2HzMwsaDRaZGVlY926j6FQKBu7ViK/V3kPk8uxRL7HrZnmnDlzXH/+85/vwfLlyzBo0CDXY9u3b8cHH6zwfnVEAYRNP0S+z+N1n127dmHAgAFVHktLS8Pu3bu9VRNRQCopKWJgEvk4j0MzMTGh2qxy5cqVSEhI8FpRRIHIYAhGZFQ0A5PIh3ncPTtv3jyMG3c33nrrLSQkJOD8+fO4cOECPvrow8aoj8ivCSFgNufBYAiCRqOFWq2RuyQiqoPHodm9e3ecOHEc27Ztw6VLmWjRIg4DBgxERER4I5RH5L8q72Hq9Qa5yyEiN9TrfZrh4eEYO3YscnJyEBUV5e2aiPwem36Imqd6XbD9wQcfQlxcC1x//Q2Ii2uBhx56GMXFxY1RH5FfKijIZ2ASNUMeh+bTTz+DX375GZs2bcQPP5zC5s2b8Msvv+CZZ2Y2Rn1EfikoKJiBSdQMebw8u3XrVnzxxReIjIwAAMTExGDZsvfRrVsqXnvtVa8XSOQvhBAoLipEUFAwVCo1VCq13CURkYc8Dk0hBBSKqrcmkiSF61qZ7nr33f/gyFeHkZ2djblz56N1UlKNx/1v5+f47383QAiBmzp0xIQJE6FS1Wsrlkg2lfcwNVotNBref5aoOfJ4eXbgwIG4556/4OjRr2EymXDkyFH87W9/q3KFIHekpKTghRdnITo6utZjsjIzsXr1Krzwwr/w5r8Xwpyfj88//8zTkolkJYSA0+l07WEyMImaL49Dc/bsWUhIaIVBgwbhhhvaYfDgwWjZMh6zZv3Lo/O0b3/TNTtvDx48gNs6JyM8IgKSJKH/gAHYn55e6/E2mw0lJSWuD4vF4lFNRN5WMcMEBPcwifyAx+ucwcHBWLhwIRYsWACTyQSj0dhod5I3mUxVZqIx0TEwmUy1Hr9hw3qsW7um2uM5piyU6PUNqsXhsMOUndmgc/gzjk/tnE4nAKCwwIxCmGWuxvfwtVM3jk/dvDU+xuhYt46r9+agw+GAVqtFYWGh67HQ0ND6ns4rRo0ajaFDh7k+t1gsmDJ5EqKMMTAYGvbmcVN2ptuDGog4PlUJIWC1XrlTCcendhybunF86tbU4+NxaB4+fBgPPfQQTp485Wr+EUJAkiTk5eV6tTij0YhLmVd+g8jKzoLRaKz1eLVaDbWaHYkkL1fTj9UKozEWSiVvm0fkLzze05w8eQqGDh2KAwe+wPHjx3D8+DGcOHEcx48f83pxXVNScOSrw8jPy4MQAp/t2IHu3bt7/XmIvKXKlX7CIxmYRH7G45lmdnY2nnjiiQbvYy55ezGOHj2C/Px8zJr1InQ6Pf69YCEWv7UInTsno3NyMmJj4zDmzrswc+bTAMqbh9L6D7jGmYnkwUvjEfk/j0NzzJgx+PTTTzFkyJAGPfGkf0yu8fHJU6ZW+TwtrT/S0vo36LmImoIQAsLpZGAS+TGPQ/OZZ55BWloa3njjzWrvsVy5ckUt30Xkvyreh6lUKhERGd1o3eREJD+PQ3PSpEnQaDRISUmBwdCwt3EQNXcVS7JOhwORUTEMTCI/53Fo7t+/Hz/8cAohISGNUQ9Rs3H1HiYDk8j/edw9e+ONN6KoqKgxaiFqNtj0QxSYPJ5pDhs2DHfeeRfuvfdexMRU3dMcPHiw1woj8mU2mxU2q5WBSRRgPA7N9957DwDwyiuvVHlckiSGJvm9igt6aDRaGKNjoVDwfZhEgcTj0PzmmxONUQeRz6tYklWp1AgJCWNgEgUgj/c0iQJRlfth8tZeRAGLoUl0DWz6IaIKDE2iaygpLmJgEhGABtwajChQGIKCodFqoVZr5C6FiGTGmSZRDcqXZHNhs1ohSRIDk4gAMDSJqqnYwywrtcApnHKXQ0Q+hKFJVAmbfoioLgxNokoKzHkMTCKqFRuBiCoxBAVDpzcwMImoRpxpUsATQqCoqABCCKjVGgYmEdWKoUkBrWIPs7ioEHa7Te5yiMjHMTQpYF3d9MO3lRDRtTA0KSCxS5aI6oONQBSw1CoNDIZgBiYRuY2hSQFFCAGb1QqNVovgkFC5yyGiZobLsxQwKpZk8/Nz4HQ65C6HiJohhiYFhMp7mGHhkbyBNBHVC0OT/B6bfojIWxia5PecTiecDicDk4gajI1A5LeEEBDCCaVSicioaEiSJHdJRNTMcaZJfqliSTYvNwdCCAYmEXkFQ5P8TuU9zOCQUAYmEXkNQ5P8Cpt+iKgxMTTJr1itZbBZrQxMImoUbAQivyCEAABotToYjbFQKPk+TCLyPs40qdlz3d6ruBAAGJhE1GgYmtSsVd7D5K29iKixMTSp2WLTDxE1NYYmNVvFRYUMTCJqUmwEomYrKDgEWq0Oag2XZYmoaXCmSc2KEALm/FzYbDZIksTAJKImxdCkZqNiD7O01ML7YRKRLBia1Cyw6YeIfAFDk5oFszmPgUlEsmMjEDULBkMQ9HoDA5OIZMWZJvksIQSKiwohhIBGo2VgEpHsGJrkkyr2MIuKCmC32+Quh4gIAEOTfNDVTT+8PB4R+QqGJvkUdskSkS9jIxD5HJVSBUNEMAOTiHwOQ5N8ghACNpsVGo0WIaHhcpdDRFQjLs+S7CqWZPPzcuB0OuUuh4ioVgxNklXlPcyw8EgoFHxJEpHv4k8okg2bfoiouWFokmycTiccdgcDk4iaDTYCUZMTQkAIAaVSiShjDCRJkrskIiK3cKZJTapy048QgoFJRM0KQ5OaTOU9zKDgEAYmETU7DE1qEmz6ISJ/wNCkJmEtK4W1zMrAJKJmjY1A1Kgq9i21Oj2M0RoolUq5SyIiqjfONKnRVCzJFhcXAgADk4iaPYYmNYrKe5gqlVrucoiIvIKhSV7Hph8i8lcMTfK6oqICBiYR+SU2ApHXBQWFQKvVQaPRyl0KEZFXyRaaFy9ewMIFC1BYWACDwYCp0+5DQkJilWO+++5bzJ41C/Hx8a7HZs2aDY2WP4x9jRACZnMegoKCoVKpGZhE5JdkC80lb7+NtLQ09OnbDwcPHMCihQsw56W51Y6Lj4/HvPmvyFAhuUsIAafTiVJLCXQ6PRt/iMhvybKnaTabcfr0L+jZqzcAoGtKCkymHFy6eLFB57XZbCgpKXF9WCwWb5RLdaho+gEE9zCJyO/JMtPMMZkQHh7het+eJEkwGo0wmUyIa9GiyrGZmZfw+GPToVAo0KdvPwwcOKjW827YsB7r1q6p4fmyUKLXN6hmh8MOU3Zmg87hjxwOBwABACgsMKMQZnkL8lF8/dSOY1M3jk/dvDU+xuhYt47z6UagpKQ2WLx4CQxBQcjJycGc2f9CSEgIUlO713j8qFGjMXToMNfnFosFUyZPQpQxBgaDoUG1mLIz3R7UQFJWVgqgPDA5PrXj66d2HJu6cXzq1tTjI8vybJTRiPz8vMuzlPIlPpPJBKPRWOU4g8EAQ1BQ+fdERaF7j544dfJkredVq9Xl33P5Q9/A2SXVTAiBkuIiCCGg1eq4JEtEAUOW0AwLC0NSUhvs27sHAHDo4EFERUVVW5rNy8uD0+kEUD5rPHrkK7ROSmryeumKij3MwkIz7Ha73OUQETUp2ZZnJ036BxYuXIANG9ZDrzdg6tRpAIDFby1C587J6JycjEMHD2DHju1QKpVwOBxI6ZaKvn37yVVywLv6Sj9qNbtkiSiwyBaa8S1bYtbsOdUenzxlquvPg/4wGIP+MLgpy6Ja8NJ4RES8jB55QKlQMjCJKKD5dPcsyU8IAbvdBrVag9CwCLnLISKSFWeaVKuKJdm8XJOrIYuIKJAxNKlGlfcww8IjoVDwpUJExJ+EVA2bfoiIasbQpGqcDgccdjsDk4joKmwEIhchBIQQUKpUiDLGQpIkuUsiIvIpnGkSgCtLsub8XAghGJhERDVgaFKVPUxDUDADk4ioFgzNAMemHyIi9zE0A1xZWSkDk4jITWwEClAV+5Y6nR7q6FgolXwpEBFdC2eaAahiSbakpBgAGJhERG5iaAaYynuYSqVS7nKIiJoVhmYAYdMPEVHDMDQDSFFhAQOTiKgBuJkVQIKCQ6DV6aDRaOUuhYioWeJM088JIVBgzoPDbodCoWBgEhE1AEPTj1XsYVosJXA4HHKXQ0TU7DE0/dTVTT8aLWeYREQNxdD0Q0IImPNz2fRDRORlbATyQ5IkQa83QG8IYmASEXkRZ5p+RAiBkpJiCCGg1ekZmEREXsaZpp+ovIep0WigUqnlLomIyO9wpukHrm76YWASETUOhmYzx0vjERE1HYamH1BICgYmEVET4J5mMyWEgN1uh1qtRlh4pNzlEBEFBM40m6GKJdm8PBOE0yl3OUREAYOh2cxU3sMMC4uApOA/IRFRU+FP3GaETT9ERPJiaDYjDocDdpudgUlEJBM2AjUDQggAgEqlgjE6FpIkyVwREVFg4kzTx1Usyebn5wIAA5OISEYMTR9WeQ/TYAiSuxwiooDH0PRRbPohIvI9DE0fVVpqYWASEfkYNgL5GCEEJEmCTqeHWq2BSsV/IiIiX8GZpg8pX5LNhcVSAkmSGJhERD6GoekjruxhlkLBq/wQEfkk/nT2AWz6ISJqHhiaPqCw0MzAJCJqBrhp5gOCg0Kg0+mh0WjlLoWIiOrAmaZMhBAoKMiHw+GAQqlkYBIRNQMMTRlU7GFaSorhsNvlLoeIiNzE0GxiVzf9aLScYRIRNRcMzSZU8T5MNv0QETVPbARqQpIkQa/Tw2AIYmASETVDDM0mIIRAqaUEOr0BOr1B7nKIiKieGJqNrPIeplqjgUqllrskIiKqJ+5pNqKrm34YmEREzRtDs5Hw0nhERP6HodmIJEgMTCIiP8I9TS8TQsDhsEOlUiM8IkrucoiIyIs40/SiiiXZvFwThBByl0NERF7G0PSSynuYoWERkCRJ7pKIiMjLGJpewKYfIqLAwND0AofdDrvNxsAkIvJzbARqgIp9S5VaDWN0HJdkiYj8HGea9VSxJGs25wEAA5OIKAAwNOuh8h6mnteSJSIKGAxND7Hph4gocMm2p3nx4gUsXLAAhYUFMBgMmDrtPiQkJFY77n87P8d//7sBQgjc1KEjJkyYCJVKvq3YUksJA5OIKEDJNtNc8vbbSEtLwxtvLsCIEaOwaOGCasdkZWZi9epVeOGFf+HNfy+EOT8fn3/+mQzVXmn60ekNiDLGMDCJiAKQLKFpNptx+vQv6NmrNwCga0oKTKYcXLp4scpxBw8ewG2dkxEeUX6xgP4DBmB/enqT1yuEgNPpRGmpBZIk8W4lREQBSpZ1zhyTCeHhEVAqlQDKO0+NRiNMJhPiWrRwHWcymRAdHe36PCY6BiaTqdbz2mw22Gw21+cWi6XBtVbsYQKCHbJERAHOr96nuWHDeqxbu6ba4zmmLJTo9R6fr2KGCZQvzRYWmFEIc0PL9EsOhx2m7Ey5y/BZHJ/acWzqxvGpm7fGxxgd69ZxsoRmlNGI/Pw8OBwOKJVKCCFgMplgNBqrHGc0GnEp88pgZGVnVTumslGjRmPo0GGuzy0WC6ZMnoQoYwwMBs/fGlJgzofFUozwiCgUFpjdHtRAZMrO5PjUgeNTO45N3Tg+dWvq8ZFlTzMsLAxJSW2wb+8eAMChgwcRFRVVZWkWKN/rPPLVYeTn5UEIgc927ED37t1rPa9arYbBYHB96Osxu6wsKDgEERFGNv0QEREAGZdnJ036BxYuXIANG9ZDrzdg6tRpAIDFby1C587J6JycjNjYOIy58y7MnPk0AKB9+5uQ1n9Ao9YlhEBRoRlBwSFQKpWufVciIiLZQjO+ZUvMmj2n2uOTp0yt8nlaWn+kpfVvkpoqX7hAq9NDo2FgEhHRFbwi0GVXX+lHo9HKXRIREfkYhiYqAjOXV/ohIqI6+dVbTupLkiTotDoYDEEMTCIiqlVAzzSFELBYSgAAegYmERFdQ8DONCvvYWrUGihlvAg8ERE1DwE507y66YeBSURE7gi40OT9MImIqL4CLjQhBCDAwCQiIo8FzLqkEAIOhwMqlQrhEVG8YwkREXksIGaaFUuyebkmCMFbfBERUf0ExEzTbM6FUqHgDJOIiBokIGaabPohIiJv8OuZphDlN4/WavVwOJwoKSmp97ksFkuDvt/fcXzqxvGpHcembhyfunlzfPR6/TVXI/06NEtLSwEADz74gMyVEBGRr3t/2QcwGAx1HiPZ7TbRRPU0OafTiby8POh0ugbtZVosFkyZPAlvLV7S4Btb+yOOT904PrXj2NSN41M3b49PwM80FQoFoqKivHY+vV5/zd9CAhnHp24cn9pxbOrG8albU45PQDQCEREReQNDk4iIyE0MTTeo1Wr8ccydUKvVcpfikzg+deP41I5jUzeOT93kGB+/bgQiIiLyJs40iYiI3MTQJCIichNDk4iIyE1+/T5NT128eAELFyxAYWEBDAYDpk67DwkJidWO+9/Oz/Hf/26AEAI3deiICRMmQqXy/6F0Z3y+++5bzJ41C/Hx8a7HZs2aDY1W29TlNql33/0Pjnx1GNnZ2Zg7dz5aJyXVeFygvnbcGZ9Afe1YrVa8/vqr+C0jAxqNBqGhYZg4cRLiWrSoduyRI1/hg+XL4HQ6kZh4HaZOu8/v37/p7vhkZWXh/vumITHxys+kR6fPQFxcnFfr8f//rR5Y8vbbSEtLQ5++/XDwwAEsWrgAc16aW+WYrMxMrF69Ci+/PA9h4eGY+/JL+PzzzzBo0B9kqrrpuDM+ABAfH49581+RoUL5pKSkYMSIkXh25tO1HhPIrx13xgcIzNcOAKSl9cctt9wKSZKwbeunWLz4LTz/zxeqHFNqsWDxW4vw/D9fQMuWrfCfpe/g43Vr8ed7/iJT1U3HnfEBAL1e1+ivHy7PXmY2m3H69C/o2as3AKBrSgpMphxcunixynEHDx7AbZ2TER4RAUmS0H/AAOxPT5ej5Cbl7vgEqvbtb7rm1acC9bUDuDc+gUqj0eDWW29zXb7thnbtkJ2dVe24r499jdatk9CyZSsAwMCBg7B/v/+/ftwdn6bCmeZlOSYTwsMjoFQqAQCSJMFoNMJkMlVZBjCZTIiOjnZ9HhMdA5PJ1OT1NjV3xwcAMjMv4fHHpkOhUKBP334YOHCQHCX7nEB97XiCrx3g008+QefOydUev/r1Ex0Tg7y8fDgcDtf/y0BQ2/gAQFlZGZ584jE4nU4kJ3fB6NF3QOHlsWFoklclJbXB4sVLYAgKQk5ODubM/hdCQkKQmtpd7tLIx/G1A6xf/zEuXbqEZ597Xu5SfFJd4xMREYHFb7+DsLAwFBUW4rXXXsXmLZsxYsRIr9bA5dnLooxG5OfnweFwACi/F6fJZILRaKxynNFoRHZ2tuvzrOysasf4I3fHx2AwwBAUVP49UVHo3qMnTp082eT1+qJAfe24K9BfO5s2bcSXhw7hqaefgbaG5qerXz/ZWVmIiAgPmFnmtcZHrVYjLCwMABAcEoK+/frh5MnvvV4HQ/OysLAwJCW1wb69ewAAhw4eRFRUVLWlx64pKTjy1WHk5+VBCIHPduxA9+7+/5uwu+OTl5cHp9MJoPy2PUePfFVrJ2mgCdTXjrsC+bWzZfMm7E9PxzMzn0XQ5V8crtap0y04c+Y0fvstAwCwffs2pHbv0ZRlysad8TGbzbDb7QAAm82GLw8dQlJr779+eBm9Si789hsWLlyAoqJC6PUGTJ06DYnXXYfFby1C587J6Jxcvo7++eefYeN/NwAob3CYOOkfAfG2AXfGZ9vWT7Fjx3YolUo4HA6kdEvFmDF3Nuh+ps3BkrcX4+jRI8jPz0dISAh0Oj3+vWAhXzuXuTM+gfraycnJwZTJkxAbGwudrvyekGq1GrPnvITVqz5CRGQkBgwYCAD46vBhrFixHA6HEwmJCbhv2v2u2bm/cnd8Dh06iDWrV0GhUMDhcKBDh4748z1/8fp1aRmaREREbuLyLBERkZsYmkRERG5iaBIREbmJoUlEROQmhiYREZGbGJpERERuYmgSERG5iaFJ5EVnz55FWFg48vPzPfq+0aPvwI4dOxqnqGs4e/YsOndORllZWa3HrFy5Ej16XLn6TNeuKdi2bZtb558zZw7uvvtut+vZunUrOnbsiPj4ltiyZQuGDBmCRYsWuf39RI2JoUkBqWPHjtiyZUuDzxMWFo4TJ0406Bx79+6FyWTCgAEDAAD79u1DWFg44uNbolWrBFx//Q24444/4pNPPmlwvd9//z2MxugqIXbdddehS5dkvPvuu26f59Chgxg0qHHuQPLkk0/h6aefxoULv2Ho0KGN8hxE9cXQJKqB3W6HEE1zsax33lmKP/1pfJXHwsJCceHCb8jIOI+vvz6KsWPvwn333Y/5DbjBrtPpxAMPPIiUlK7VvjZu3DgsWfJOvc/tTWfPnkX79u3lLoOoRgxNCjj33PMXnD+fgXvvnYD4+JZ46KGHAZTPGpcsWYKUlG5o0SIeRUVF1WaSixYtwpAhQwAAffv2AwAMGDAQ8fEtqwTatm3b0KnTLUhMTMSUKVNgs9lqrMVms2Hnzp3o1atXrfWGhIRgzJgxmDdvHubOnYvc3Lx6/b0XL16MG29sV+NF4lNSUnDhwgX88MMPbp2r8ky9Yul27ty5aNv2elx//Q11Lqe+8MKL6N69Oy5dulTl8dzcXMTHt4TT6XSNaU1Lxjt3/g89evREQkIievbshV27dgMAMjMzYTRGo6ioCADw9ttvIywsHD/++COA8mXfbt1S3fr7EdWGoUkBZ/nyZUhIaIX//GcpLlz4Da+//prra2vXrsOGDeuRkXG+1rspVNi1638AgB07tuPChd8wffqjrq999tnn2LdvLw4dOoQ9e/ZizZo1NZ7jl19+QUlJCW644YZr1j18+DDYbDYcOfIVAODVV19DYmJirR9r1651fe+5c+fw1luL8eKLL9Z4brVajTZt2uCbb765Zh01OXnyFPR6PU6dOon33nsXM2c+i9Onz1Q5xm63Y9q0+3Do0EF8+umniIuLq/L1yMhIXLjwG4ArY3r1LaB++eU07r77bjz22AycOXMajz76KMaNG4dff/0VsbGxaNOmDQ4cOACgfNk7KSkJe/fuc33eq1fPev39iCowNIkqefDBB9CiRQtotVooFPX/7/H4448hJCQELVq0wO23345jx47VeFx+fj4MBoNb90TUaDSIiopCXl75TPORRx7GuXPnav0YM2aM63sfeuhhPP30U4iMjKz1/CEhIcjLy/fo71khKioK999/P9RqNXr27InExER8882VGXpJiQXjx49HQUEB1q9f77rvoafWr1+PHj16YPjw4VCpVBg5cgRSUlKwbt3HAICePXti7959cDqdOHjwEKZPfxT79lUOzdpn9ETuYGgSVdKqVSuvnCcmJsb156Agg2vJ8Grh4eEoKSlx3dy7LlarFTk5OYiIiPColtWrV8Nut2Ps2LF1HldYWIiIiHCPzl0hJia6yudX/52/+eYb7Nq1G08++USNNxB214ULF5CYmFjlsdatW+PChQsAykNz3759OH78BK677joMHjwEX3zxBUwmE06d+oH3L6UGY2hSQJKkml/6V88ug4KCYLFYXJ9fupR51Xkadq/Htm3bwmAw4KeffrrmsZs2bYZGo0HnzuX35pw//xXEx7es9aNiSXj37t04cuQIkpLaICmpDd5440189tnnuOGGdq5z22w2nD59Gh07dmzQ36c2Xbt2wfz58zBy5CicPHmy3ueJj4/HuXPnqjx27tw5xMfHAwB69uyBb775Blu2bEGvXr0QGRmBuLg4LFmyBB06dEB4eHhD/hpEDE0KTDExMThz5sw1j7v55t9j1arymdqJEyewevXqep2nNmq1Gv369XMtIdakqKgIH3+8Ho899hhmzJjhmg1On/4oLlz4rdaPO++8EwAwe/YcfPnlIaSn70N6+j78/e9/Q8+ePbFnz27Xcxw6dAgtWrTAjTfeWO+/y7Xcc889eO65ZzF8+Ah8++239TrH6NGjkZ6ejk8++QR2ux2bNm3CF198gTvuuANA+TLxjTe2w5IlS9CzZ/n+Za9evfDWW4u5NEtewdCkgPToo49gyZJ3kJiYiEceebTW4+bOnYvDh79EYuJ1eO655zFu3LgqX3/66afx+ONPIDHxOrz66mu1nKVuEydOwMqVH1Z5zGwucL1Ps1OnW7By5Uq8+eYbVZqN3BUREY6WLVu6PkJCQqDTaV2zMwD46KNVmDhxQr3q98T48ePx4osvYOTIUfV6f2vbtm2wYsUHmDNnDlq3TsLcuXOxYsUKJCW1dh3Ts2dPlJaWolu3FABAnz69UVBQgN69GZrUcJLdbmuaN6MRUa1GjRqNqVOnoH///k3+3OfOncMdd/wR6en7GrTfSBQIGJpERERu4vIsERGRmxiaREREbmJoEhERuYmhSURE5CaGJhERkZsYmkRERG5iaBIREbmJoUlEROQmhiYREZGbGJpERERu+n8kVZJZSbKbdgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 480x460 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "task_full = _dn_task(sc2, truth2, np.arange(4), prior_d=3.0, n_periods=400, scale=50.0, rho=0.5)\n",
    "trace_full = ODTrace()\n",
    "DavisNihanKalmanEstimator(k_inner=120, outer_iters=60).estimate(\n",
    "    task_full, Budget(sp_calls=10**9, iterations=200), RngBundle(0), trace_full\n",
    ")\n",
    "recovered_trace = Trace()\n",
    "from tabench import Scenario\n",
    "BiconjugateFrankWolfeModel().solve(\n",
    "    Scenario(\"recovered\", sc2.network, Demand(trace_full.final.od_matrix), family=sc2.family),\n",
    "    Budget(iterations=5000, target_relative_gap=1e-10), RngBundle(0), recovered_trace,\n",
    ")\n",
    "display(viz.plot_network_flows(sc2.network, truth2))\n",
    "display(viz.plot_flow_scatter(\n",
    "    (\"truth (D=4)\", truth2),\n",
    "    {\"od-kalman recovered\": recovered_trace.final.link_flows},\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f75f2551",
   "metadata": {},
   "source": [
    "## Takeaways & pointers\n",
    "\n",
    "- **Certified, not self-reported.** Both the closed forms and the Braess\n",
    "  `od_rmse` above were recomputed independently — the closed forms\n",
    "  algebraically, the Braess recovery by `ODCertifier` inside\n",
    "  `run_estimation_experiment`.\n",
    "- **The spatial covariance is load-bearing, not decorative.** Off-diagonal\n",
    "  whitening moved the estimate materially away from `gls`'s diagonal-V\n",
    "  assumption — verified directly, not asserted from the paper.\n",
    "- **Where next.** the per-cell-covariance sibling [`gls`](01-gls.ipynb)\n",
    "  (diagonal-Sigma special case, verified exactly above); the within-day\n",
    "  dynamic cousin [`od-dynamic`](07-od-dynamic.ipynb); 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": "estimation",
   "unit": "od-kalman"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}