{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c0",
   "metadata": {},
   "source": [
    "# Competition indices on a simulated stand\n",
    "\n",
    "This notebook walks a full individual-tree competition workflow:\n",
    "\n",
    "1. generate a randomly located stand on two circular plots,\n",
    "2. fit the Näslund height curve and impute the heights that were never measured,\n",
    "3. compute several competition indices and compare competitor-selection rules.\n",
    "\n",
    "The stand is synthetic so the notebook is self-contained, but nothing here is\n",
    "specific to simulated data -- swap in a real stem map and the same calls apply."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1",
   "metadata": {},
   "source": [
    "## Notebook objectives\n",
    "\n",
    "- Show that a *partially measured* inventory is the normal case, and how\n",
    "  imputation fills the gaps without losing track of what was measured.\n",
    "- Show that competition indices are chosen **together with** a competitor-selection\n",
    "  rule, and that changing the rule changes the answer.\n",
    "- Show where each number comes from: every index and every published selection\n",
    "  rule carries its own citation.\n",
    "\n",
    "## Prerequisites\n",
    "\n",
    "- `pyforestry` installed from this repository.\n",
    "- Run the cells in order. All randomness is seeded, so the numbers below are\n",
    "  reproducible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c2",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-26T12:01:42.065054Z",
     "iopub.status.busy": "2026-07-26T12:01:42.065054Z",
     "iopub.status.idle": "2026-07-26T12:01:49.886004Z",
     "shell.execute_reply": "2026-07-26T12:01:49.886004Z"
    }
   },
   "outputs": [],
   "source": [
    "import random\n",
    "from math import cos, pi, sin, sqrt\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "\n",
    "from pyforestry.base.competition import (\n",
    "    INDEX_REGISTRY,\n",
    "    BitterlichBAF,\n",
    "    FixedRadius,\n",
    "    LeeGadowRadius,\n",
    "    MeanHeightRadius,\n",
    "    NearestNeighbours,\n",
    "    SearchCone,\n",
    "    competition_indices,\n",
    "    index_source,\n",
    "    selector_source,\n",
    ")\n",
    "from pyforestry.base.helpers import CircularPlot, Stand, Tree\n",
    "from pyforestry.base.helpers.primitives import Position\n",
    "from pyforestry.base.helpers.tree_species import TreeSpecies\n",
    "\n",
    "RNG = random.Random(20260726)  # seeded: every number below is reproducible"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3",
   "metadata": {},
   "source": [
    "## 1. Generate a random stand\n",
    "\n",
    "Two circular plots of 10 m radius. Stems are placed uniformly over each plot --\n",
    "note the `sqrt` on the radius, without which points would bunch toward the\n",
    "centre -- and diameters are drawn from a lognormal, which gives the usual\n",
    "right-skewed diameter distribution.\n",
    "\n",
    "Heights are the realistic part: in a real inventory only a **sample** of trees is\n",
    "measured for height. Here roughly one tree in four gets a measured height, and\n",
    "the rest are left blank for step 2 to fill in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-26T12:01:49.889268Z",
     "iopub.status.busy": "2026-07-26T12:01:49.889268Z",
     "iopub.status.idle": "2026-07-26T12:01:49.908506Z",
     "shell.execute_reply": "2026-07-26T12:01:49.908506Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "90 trees on 2 plots\n",
      "21 have a measured height, 69 do not\n"
     ]
    }
   ],
   "source": [
    "PLOT_RADIUS_M = 10.0\n",
    "TREES_PER_PLOT = 45\n",
    "HEIGHT_SAMPLE_RATE = 0.25\n",
    "\n",
    "SPECIES = TreeSpecies.Sweden.picea_abies\n",
    "\n",
    "\n",
    "def true_height(diameter_cm: float) -> float:\n",
    "    \"\"\"A 'true' height-diameter relationship used only to generate the sample.\"\"\"\n",
    "    return 1.3 + diameter_cm**2 / (1.15 + 0.20 * diameter_cm) ** 2\n",
    "\n",
    "\n",
    "def make_plot(plot_id: int, centre: Position) -> CircularPlot:\n",
    "    \"\"\"Place TREES_PER_PLOT stems uniformly on a circular plot.\"\"\"\n",
    "    trees = []\n",
    "    for k in range(TREES_PER_PLOT):\n",
    "        # Uniform over the disc: radius must be scaled by sqrt(u).\n",
    "        r = PLOT_RADIUS_M * sqrt(RNG.random())\n",
    "        theta = RNG.uniform(0.0, 2.0 * pi)\n",
    "        # Two cohorts: a main canopy plus a suppressed understory. Real stands\n",
    "        # have both, and without the understory the 0.3 * d_i size screen has\n",
    "        # nothing to remove.\n",
    "        if RNG.random() < 0.30:\n",
    "            diameter = min(45.0, max(4.0, RNG.lognormvariate(2.05, 0.30)))\n",
    "        else:\n",
    "            diameter = min(45.0, max(4.0, RNG.lognormvariate(3.05, 0.28)))\n",
    "\n",
    "        measured = None\n",
    "        if RNG.random() < HEIGHT_SAMPLE_RATE:\n",
    "            # Measured heights carry a little observation noise.\n",
    "            measured = round(true_height(diameter) + RNG.gauss(0.0, 0.4), 2)\n",
    "\n",
    "        trees.append(\n",
    "            Tree(\n",
    "                position=(centre.X + r * cos(theta), centre.Y + r * sin(theta)),\n",
    "                species=SPECIES,\n",
    "                diameter_cm=round(diameter, 1),\n",
    "                height_m=measured,\n",
    "                uid=f\"p{plot_id}t{k:02d}\",\n",
    "            )\n",
    "        )\n",
    "    return CircularPlot(id=plot_id, position=centre, radius_m=PLOT_RADIUS_M, trees=trees)\n",
    "\n",
    "\n",
    "stand = Stand(plots=[make_plot(1, Position(0.0, 0.0)), make_plot(2, Position(40.0, 0.0))])\n",
    "all_trees = [t for p in stand.plots for t in p.trees]\n",
    "\n",
    "n_measured = sum(1 for t in all_trees if t.height_m is not None)\n",
    "print(f\"{len(all_trees)} trees on {len(stand.plots)} plots\")\n",
    "print(f\"{n_measured} have a measured height, {len(all_trees) - n_measured} do not\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5",
   "metadata": {},
   "source": [
    "## 2. Fit and apply the Näslund height imputer\n",
    "\n",
    "`Stand.impute` fits Näslund's height-diameter curve to the stand's *own*\n",
    "measured pairs and writes a modelled height for every tree that lacks one.\n",
    "\n",
    "The measured heights are never touched. A modelled value goes into\n",
    "`Tree.imputed` alongside the imputer that produced it and its citation, and\n",
    "`Tree.value_of` reads the two together:\n",
    "\n",
    "- `tree.height_m` -- the measurement, or `None`\n",
    "- `tree.value_of(\"height_m\")` -- measured if present, else imputed\n",
    "- `tree.provenance(\"height_m\")` -- `\"measured\"` or `\"imputed\"`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c6",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-26T12:01:49.910231Z",
     "iopub.status.busy": "2026-07-26T12:01:49.910231Z",
     "iopub.status.idle": "2026-07-26T12:01:49.926112Z",
     "shell.execute_reply": "2026-07-26T12:01:49.925611Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "imputed a height for 69 trees\n",
      "\n",
      "example tree p1t02: dbh 5.6 cm\n",
      "  height_m (measured) : None\n",
      "  value_of            : 7.47 m\n",
      "  provenance          : imputed\n",
      "  produced by         : naslund_height_imputer\n",
      "  citation            : Näslund, M. (1936)\n",
      "                        Skogsförsöksanstaltens gallringsförsök i tallskog\n"
     ]
    }
   ],
   "source": [
    "assigned = stand.impute(\"height_m\")\n",
    "print(f\"imputed a height for {assigned} trees\")\n",
    "\n",
    "example = next(t for t in all_trees if t.provenance(\"height_m\") == \"imputed\")\n",
    "source = example.imputed_source(\"height_m\")\n",
    "print(f\"\\nexample tree {example.uid}: dbh {example.diameter_cm} cm\")\n",
    "print(f\"  height_m (measured) : {example.height_m}\")\n",
    "print(f\"  value_of            : {example.value_of('height_m'):.2f} m\")\n",
    "print(f\"  provenance          : {example.provenance('height_m')}\")\n",
    "print(f\"  produced by         : {example.imputed['height_m'].imputer_id}\")\n",
    "print(f\"  citation            : {source.author} ({source.year})\")\n",
    "print(f\"                        {source.title}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7",
   "metadata": {},
   "source": [
    "Every tree now has a usable height, and each one still knows which kind it is.\n",
    "That distinction matters downstream: a model fitted on measured heights should\n",
    "not silently be fed interpolated ones without the user knowing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c8",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-26T12:01:49.926112Z",
     "iopub.status.busy": "2026-07-26T12:01:49.926112Z",
     "iopub.status.idle": "2026-07-26T12:01:49.941817Z",
     "shell.execute_reply": "2026-07-26T12:01:49.941817Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            count   mean   min    max\n",
      "provenance                           \n",
      "imputed        69  15.00  6.43  20.94\n",
      "measured       21  13.83  8.15  17.78\n"
     ]
    }
   ],
   "source": [
    "summary = pd.DataFrame(\n",
    "    {\n",
    "        \"uid\": [t.uid for t in all_trees],\n",
    "        \"dbh_cm\": [float(t.diameter_cm) for t in all_trees],\n",
    "        \"height_m\": [t.value_of(\"height_m\") for t in all_trees],\n",
    "        \"provenance\": [t.provenance(\"height_m\") for t in all_trees],\n",
    "    }\n",
    ")\n",
    "print(summary.groupby(\"provenance\")[\"height_m\"].agg([\"count\", \"mean\", \"min\", \"max\"]).round(2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-26T12:01:49.941817Z",
     "iopub.status.busy": "2026-07-26T12:01:49.941817Z",
     "iopub.status.idle": "2026-07-26T12:01:50.240637Z",
     "shell.execute_reply": "2026-07-26T12:01:50.240637Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArEAAAG4CAYAAABSPb94AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAACEmElEQVR4nO3deVxUVf8H8M+wzbCOgCCgCLih5J6ZmOa+Z9qqWS7ZYqaZWm5lipZ7aYtPmvVLMy3tyTRNH3ekUlyJzCVXFFMQRNkGhmXm/P64zcAwMzDIwAzweb9e89K5c+/Mmct1+Hjme86RCSEEiIiIiIiqEQdbN4CIiIiIqLwYYomIiIio2mGIJSIiIqJqhyGWiIiIiKodhlgiIiIiqnYYYomIiIio2mGIJSIiIqJqhyGWiIiIiKodhlgiIiIiqnYYYmuIdevWQSaTQSaT4dChQ0aPCyHQpEkTyGQydO/evcrbZ4/GjBmD0NDQMvcLDQ3FY489ZtXXlslkiIqKuq9jLW3PuXPnEBUVhWvXrt3X69iaJefo1q1biIqKQnx8vNFjY8aMgYeHR+U0rpioqCjIZDL4+/sjKyvL6HFTP69Dhw5BJpPhxx9/BCC914kTJ97X61fkWqoo3edOdb3GKmL27Nlo2LAhnJycUKdOHQBA9+7dDT5fc3JyEBUVZfIz+ciRI4iKikJ6errV22bpZ9vnn3+OdevWGW0veX1WN9euXYNMJsOHH354X8eX9tnZvXt3tGzZsoItJGthiK1hPD098X//939G22NiYnDlyhV4enraoFVUUmxsLF5++eVKfY1z585h3rx5NTpg3Lp1C/PmzTMZYqtaamoqli5datG+7du3R2xsLHr27AlAuh7efvvtymweWdHPP/+MBQsWYNSoUYiJicH+/fsBSKHw888/1++Xk5ODefPmmQ2x8+bNq5QQaylzIba2qw2fnTUFQ2wNM2zYMGzZsgWZmZkG2//v//4PkZGRaNiwoY1aZj25ubkQQti6GRXSqVMnNGjQwNbNICvq378/VqxYgeTk5DL39fLyQqdOneDj4wNAuh4s6Tkj04QQyM3NrbLXO3PmDABg0qRJeOSRR9ChQwcAQEREBCIiIqqsHfYqJyfH1k2gWoIhtoZ57rnnAADff/+9fltGRga2bNmCsWPHmjwmPz8fH3zwAZo3bw65XA4/Pz+8+OKLSE1NNdhv8+bN6Nu3LwIDA+Hq6ooWLVpg5syZUKlUBvtdvXoVw4cPR1BQEORyOerVq4devXoZ9JaZ+wo0NDQUY8aM0d/XfV25d+9ejB07Fn5+fnBzc0NeXp6+TZGRkXB3d4eHhwf69euHP/74w+h5161bh/DwcMjlcrRo0QLr168v9Tyasnv3brRv3x6urq5o3rw5vv76a6N9kpOTMW7cODRo0AAuLi4ICwvDvHnzUFhYaLCfqff/+++/IzIyEgqFAvXr18d7772Hr776yuzXtaW1Z926dXjmmWcAAD169NCXmuh6Xf744w889thj8Pf3h1wuR1BQEAYNGoR//vmn1HOwb98+DBkyBA0aNIBCoUCTJk0wbtw43Llzx2A/3VfsZ8+exXPPPQelUol69eph7NixyMjIMNg3MzMTr7zyCnx9feHh4YH+/fvj4sWLpbYDkL7yfOihhwAAL774ov49ljyvly9fxsCBA+Hh4YHg4GC89dZb+utHx9J/A6X54IMPUFhYaNFX+3PmzEHHjh3h4+OjvyaXL19u9J+zgwcPonv37vD19YWrqysaNmyIp556qtSQoDv3JZn66l9X6mDJtX306FE88sgjUCgUCAoKwqxZs1BQUFDme9U5duwYBg8eDF9fXygUCjRu3BiTJ0/WP27uK3BT70dXfrF69Wq0aNECcrkcX331Ffz9/TFy5Eij50hPT4erqyumTp2q35aZmYm3334bYWFhcHFxQf369TF58mSjz7OSQkNDMXv2bABAvXr1DK654uUE165dg5+fHwBg3rx5+utzzJgxiIqKwrRp0wAAYWFhJkvBKvuzLTQ0FGfPnkVMTIz+9Uue/4KCArz77rsICgqCl5cXevfujQsXLhjso/t6/ddff0Xnzp3h5uam/12TmJiIF154Qf8506JFC3z00UfQarX643WlCyV7q3UlASV7ir/88ks0a9YMcrkcERER+O6770otn1i+fDnCwsLg4eGByMhIHD16tNTzUtZnp86JEyfQtWtXuLm5oVGjRli8eLHB+wIsv8Z01/PatWsRHh4OV1dXdOjQAUePHoUQAsuWLdO/h549e+Ly5culvodaRVCNsHbtWgFAnDhxQowcOVJ07NhR/9iqVauEu7u7yMzMFA888IDo1q2b/jGNRiP69+8v3N3dxbx588S+ffvEV199JerXry8iIiJETk6Oft/3339frFixQuzcuVMcOnRIrF69WoSFhYkePXoYtCU8PFw0adJEfPvttyImJkZs2bJFvPXWWyI6Olq/DwAxd+5co/cREhIiRo8ebfS+6tevL1599VXxv//9T/z444+isLBQLFiwQMhkMjF27Fjxyy+/iJ9++klERkYKd3d3cfbsWaPnGDJkiNixY4fYsGGDaNKkiQgODhYhISFlntuQkBDRoEEDERERIdavXy/27NkjnnnmGQFAxMTE6PdLSkrSP+cXX3wh9u/fL95//30hl8vFmDFjDJ6z5Pv/888/hUKhEK1btxabNm0S27dvFwMHDhShoaECgEhISChXe1JSUsTChQsFAPGf//xHxMbGitjYWJGSkiKys7OFr6+v6NChg/jhhx9ETEyM2Lx5s3jttdfEuXPnSj0Xq1atEosWLRLbt28XMTEx4ptvvhFt2rQR4eHhIj8/X7/f3LlzBQARHh4u5syZI/bt2yeWL18u5HK5ePHFF/X7abVa0aNHDyGXy8WCBQvE3r17xdy5c0WjRo3MXiM6GRkZ+p/t7Nmz9e/xxo0bQgghRo8eLVxcXESLFi3Ehx9+KPbv3y/mzJkjZDKZmDdvnv55yvNvwBTde01NTRVTpkwRTk5O4sKFCwY/r0GDBhkcM2rUKPH111+Lffv2ib1794r3339fuLq6GrQrISFBKBQK0adPH7Ft2zZx6NAhsXHjRjFy5Ehx7949/X4lz5OuPSXpzlV5ryUhhDh79qxwc3MTERER4vvvvxc///yz6Nevn2jYsKHRc5qye/du4ezsLFq3bi3WrVsnDh48KL7++msxfPhw/T6jR482+e/R1PvRfSa0bt1afPfdd+LgwYPizJkzYsqUKcLV1VVkZGQY7P/5558LAOL06dNCCCFUKpVo27atqFu3rli+fLnYv3+/+OSTT4RSqRQ9e/YUWq3W7HuJi4sTL730kgAgdu/ebXDNdevWTf/5qlarxe7duwUA8dJLL+mvz8uXL4sbN26IN954QwAQP/30k/4xXbur4rMtLi5ONGrUSLRr107/+nFxcUIIIaKjowUAERoaKp5//nmxc+dO8f3334uGDRuKpk2bisLCQv3zdOvWTfj4+Ijg4GDx2WefiejoaBETEyNSUlJE/fr1hZ+fn1i9erXYvXu3mDhxogAgxo8frz9e91rFfz8IIV3/AMTatWv127744gsBQDz11FPil19+ERs3bhTNmjUTISEhBu9Xd2xoaKjo37+/2LZtm9i2bZto1aqV8Pb2Funp6WbPS2mfnbr36+vrK5o2bSpWr14t9u3bJ15//XUBQHzzzTf65ynPNQZAhISEiM6dO4uffvpJbN26VTRr1kz4+PiIKVOmiCFDhujfb7169UTr1q1LvUZrE4bYGqJ4iNV9KJw5c0YIIcRDDz2kD1ElQ+z3338vAIgtW7YYPN+JEycEAPH555+bfD2tVisKCgpETEyMACD+/PNPIYQQd+7cEQDExx9/XGp7yxtiR40aZbBfYmKicHJyEm+88YbB9qysLBEQECCeffZZIYQUUIKCgkT79u0N/tFfu3ZNODs7WxxiFQqFuH79un5bbm6u8PHxEePGjdNvGzdunPDw8DDYTwghPvzwQwHA4JdPyff/zDPPCHd3d5GamqrfptFoREREhMngYUl7/vvf/5r85XDy5EkBQGzbtq3M914a3TVw/fp1AUD8/PPP+sd0wWPp0qUGx7z++utCoVDofxb/+9//BADxySefGOy3YMGCMkOsEEXXafFfdDqjR48WAMQPP/xgsH3gwIEiPDxcf/9+/w2UfK+pqanizp07QqlUiqeeekr/uKkQW5xGoxEFBQVi/vz5wtfXV39ufvzxRwFAxMfHl/r6FQ2xllxLw4YNE66uriI5OVm/rbCwUDRv3tyiENu4cWPRuHFjkZuba3af8oZYpVIp7t69a7D99OnTAoBYs2aNwfaOHTuKBx98UH9/0aJFwsHBQZw4ccJgP90537VrV6nvp/jPvLjiIVYIIVJTU81ex8uWLTN57qrys63k7wMd3e+QgQMHGmz/4YcfBAARGxtr8J4BiAMHDhjsO3PmTAFAHDt2zGD7+PHjhUwm0/9Hz9IQq9FoREBAgHj44YcN9rt+/brR+9Ud26pVK4PAffz4cQFAfP/996WeF3OfncXfb8n3FRERIfr166e/X55rDIAICAgQ2dnZ+m3btm0TAETbtm0Nfr4ff/yxwX/IajuWE9RA3bp1Q+PGjfH111/jr7/+wokTJ8yWEvzyyy+oU6cOBg8ejMLCQv2tbdu2CAgIMPiK5+rVqxgxYgQCAgLg6OgIZ2dndOvWDQBw/vx5AICPjw8aN26MZcuWYfny5fjjjz+MvmK5H0899ZTB/T179qCwsBCjRo0yaLdCoUC3bt307b5w4QJu3bqFESNGGHwlGRISgs6dO1v8+m3btjWoJ1YoFGjWrBmuX7+u3/bLL7+gR48eCAoKMmjTgAEDAEiD68yJiYlBz549UbduXf02BwcHPPvss/fdHnOaNGkCb29vzJgxA6tXr8a5c+fKPEYnJSUFr732GoKDg+Hk5ARnZ2eEhIQAKLoGinv88ccN7rdu3RpqtRopKSkAgOjoaADA888/b7DfiBEjLG5TaWQyGQYPHmzUhpI/N0v/DZTF19cXM2bMwJYtW3Ds2DGz+x08eBC9e/eGUqnU/1uaM2cO0tLS9Oembdu2cHFxwauvvopvvvkGV69eLd+bt5Al11J0dDR69eqFevXq6bc5Ojpi2LBhZT7/xYsXceXKFbz00ktQKBRWa3fPnj3h7e1tsK1Vq1Z48MEHsXbtWv228+fP4/jx4wafgb/88gtatmyJtm3bGvzM+/XrZ3aGl6pS1Z9tpTH17xeA0eeMt7e3fpCizsGDBxEREYGOHTsabB8zZgyEEDh48GC52nLhwgUkJycbfSY2bNgQjzzyiMljBg0aBEdHxzLbX14BAQFG78vU50p5rrEePXrA3d1df79FixYAgAEDBhj8fHXbK/oeagqG2BpIJpPhxRdfxIYNG7B69Wo0a9YMXbt2Nbnv7du3kZ6eDhcXFzg7OxvckpOT9bWO2dnZ6Nq1K44dO4YPPvgAhw4dwokTJ/DTTz8BgH5QhUwmw4EDB9CvXz8sXboU7du3h5+fHyZNmmRy+iFLBQYGGrUbAB566CGjdm/evFnf7rS0NADSh05JpraZ4+vra7RNLpcbDCa5ffs2duzYYdSeBx54AACM6kaLS0tLMwgIOqa2Wdoec5RKJWJiYtC2bVu88847eOCBBxAUFIS5c+eWWuOo1WrRt29f/PTTT5g+fToOHDiA48eP62vMTL12yXbK5XKDfdPS0uDk5GS0X3l+NqVxc3MzCk5yuRxqtVp/39J/A5aaPHkygoKCMH36dJOPHz9+HH379gUg1fcdPnwYJ06cwLvvvgug6Nw0btwY+/fvh7+/PyZMmIDGjRujcePG+OSTT8rVnrJYci2lpaXd978hXV2xtQcylvxM0Bk7dixiY2Px999/AwDWrl0LuVyuHy8ASD/z06dPG/28PT09IYQo98/cmqr6s600Zf371TH1s0hLSzO5PSgoSP94eej2r8jnpLn2l5elvw/Kc43pBnnquLi4lLq9+GdYbeZk6wZQ5RgzZgzmzJmD1atXY8GCBWb3q1u3Lnx9fbF7926Tj+um5Dp48CBu3bqFQ4cO6XtfAZicHiYkJEQ/zdfFixfxww8/ICoqCvn5+Vi9ejUA6R98ycE1gPkPtpIDO3Q9lj/++KO+J9AU3YeNqRHjlowiL4+6deuidevWZs+37sPbFF9fX/0vr+Ks3UadVq1aYdOmTRBC4PTp01i3bh3mz58PV1dXzJw50+QxZ86cwZ9//ol169Zh9OjR+u0VGWTg6+uLwsJCpKWlGfxiqKz3bYql/wYs5erqiqioKLz66qvYuXOn0eObNm2Cs7MzfvnlF4OAvW3bNqN9u3btiq5du0Kj0eDkyZP47LPPMHnyZNSrVw/Dhw83+fq658zLy9P/0gZK/09UWXx9fe/735BucFNZgwYVCoXJzwRz7TY1eA2QBrdOnToV69atw4IFC/Dtt99i6NChBr22devWhaurq8kBbLrHbcUeP9vKYupn4evri6SkJKPtt27dAlD0Potfr8WV/Lnr3m9Vfk5WhD1fYzUJQ2wNVb9+fUybNg1///23QeAo6bHHHsOmTZug0Wjw8MMPm91P9yFV/JciAHzxxReltqNZs2aYPXs2tmzZgri4OP320NBQnD592mDfgwcPIjs7u9Tn0+nXrx+cnJxw5coVo1KD4sLDwxEYGIjvv/8eU6dO1b+P69ev48iRI6UGy/J67LHHsGvXLjRu3Njoa86ydOvWDbt27cKdO3f0H25arRb//e9/77s9lvQ6yGQytGnTBitWrMC6desMfkam9i3+vDplXQOl6dGjB5YuXYqNGzdi0qRJ+u3fffedRcdbo2fF0n8D5TF27FisWLECM2fONCqnkclkcHJyMviaMzc3F99++63Z53N0dMTDDz+M5s2bY+PGjYiLizMbYnWjtE+fPq2fvQEAduzYcd/vp0ePHti+fTtu376t7/XSaDTYvHlzmcc2a9ZMX940depUo+uneLtTUlIMXiM/Px979uwpV1u9vb0xdOhQrF+/HpGRkUhOTjYqp3rsscewcOFC+Pr6IiwsrFzPXx6lXZ/mHqvKzzZLv725H7169cKiRYsQFxeH9u3b67evX78eMpkMPXr0AGB4vfbr10+/3/bt2w2eLzw8HAEBAfjhhx8MZplITEy0+me5tT5XquIaq+0YYmuwxYsXl7nP8OHDsXHjRgwcOBBvvvkmOnbsCGdnZ/zzzz+Ijo7GkCFD8MQTT6Bz587w9vbGa6+9hrlz58LZ2RkbN27En3/+afB8p0+fxsSJE/HMM8+gadOmcHFxwcGDB3H69GmDHr6RI0fivffew5w5c9CtWzecO3cOK1euhFKptOi9hYaGYv78+Xj33Xdx9epV9O/fH97e3rh9+zaOHz8Od3d3zJs3Dw4ODnj//ffx8ssv44knnsArr7yC9PR0REVFWe0rN5358+dj37596Ny5MyZNmoTw8HCo1Wpcu3YNu3btwurVq81+pfruu+9ix44d6NWrF9599124urpi9erV+qlYHBzKX/mjW1VmzZo18PT0hEKhQFhYGGJjY/H5559j6NChaNSoEYQQ+Omnn5Ceno4+ffqYfb7mzZujcePGmDlzJoQQ8PHxwY4dO7Bv375yt02nb9++ePTRRzF9+nSoVCp06NABhw8fLjXQFde4cWO4urpi48aNaNGiBTw8PBAUFFSuX2iW/hsoD0dHRyxcuFB/nK4WD5Dq9JYvX44RI0bg1VdfRVpaGj788EOjcLd69WocPHgQgwYNQsOGDaFWq/W9Or179zb72gMHDoSPjw9eeuklzJ8/H05OTli3bh1u3LhRrvdQ3OzZs7F9+3b07NkTc+bMgZubG/7zn/+UOR2Vzn/+8x8MHjwYnTp1wpQpU9CwYUMkJiZiz5492LhxIwBpjus5c+Zg+PDhmDZtGtRqNT799FNoNJpyt3fs2LHYvHkzJk6ciAYNGhidr8mTJ2PLli149NFHMWXKFLRu3RparRaJiYnYu3cv3nrrLav8h8bT0xMhISH4+eef0atXL/j4+KBu3boIDQ1Fq1atAACffPIJRo8eDWdnZ4SHh1fpZ5vuG5nNmzejUaNGUCgU+nZV1JQpU7B+/XoMGjQI8+fPR0hICHbu3InPP/8c48ePR7NmzQBIpQ+9e/fGokWL4O3tjZCQEBw4cEBfqqbj4OCAefPmYdy4cXj66acxduxYpKenY968eQgMDLyvz0hzzH12miojMKeqrrFaz5ajysh6is9OUBpTo1ELCgrEhx9+KNq0aSMUCoXw8PAQzZs3F+PGjROXLl3S73fkyBERGRkp3NzchJ+fn3j55ZdFXFycwQjS27dvizFjxojmzZsLd3d34eHhIVq3bi1WrFhhMEo0Ly9PTJ8+XQQHBwtXV1fRrVs3ER8fb3Z2AnPva9u2baJHjx7Cy8tLyOVyERISIp5++mmxf/9+g/2++uor0bRpU+Hi4iKaNWsmvv76a7OjoUsyN7q85EhkIaTRyJMmTRJhYWHC2dlZ+Pj4iAcffFC8++67BiNPYWLE8m+//SYefvhhIZfLRUBAgJg2bZpYsmSJAGAwJUx52vPxxx+LsLAw4ejoqP85/f333+K5554TjRs3Fq6urkKpVIqOHTuKdevWlXkuzp07J/r06SM8PT2Ft7e3eOaZZ0RiYqLZEfIlR2+bGiGfnp4uxo4dK+rUqSPc3NxEnz59xN9//23R7ARCSLMLNG/eXDg7OxscM3r0aOHu7m60v6nR7pb+GzDF3HsVQojOnTsLAEY/r6+//lqEh4cLuVwuGjVqJBYtWiT+7//+z+DcxMbGiieeeEKEhIQIuVwufH19Rbdu3cT27dsNnsvUeTp+/Ljo3LmzcHd3F/Xr1xdz584VX331lcnZCSy9lg4fPiw6depkcH2uWbPGotkJdO9nwIABQqlUCrlcLho3biymTJlisM+uXbtE27Zthaurq2jUqJFYuXKl2dkJJkyYYPa1NBqNCA4OFgDEu+++a3Kf7OxsMXv2bBEeHi5cXFyEUqkUrVq1ElOmTDGYhcEUS2cnEEKI/fv3i3bt2gm5XC4AGHy+zZo1SwQFBQkHBwej0fBV8dl27do10bdvX+Hp6amf5kmIohkD/vvf/xrsb2raq27duokHHnjA5PNfv35djBgxQvj6+gpnZ2cRHh4uli1bJjQajcF+SUlJ4umnnxY+Pj5CqVSKF154QT+LSsmZR9asWSOaNGli8H6HDBki2rVrZ9TOZcuWGbXJ0s8VU5+dpb1fU+fc0mvM1PVs7j2Y+9nUVjIhqvnSR0Q1WN++fXHt2jWLJv8nIqpt0tPT0axZMwwdOhRr1qyxdXOoirGcgMhOTJ06Fe3atUNwcDDu3r2LjRs3Yt++ffpBckREtVlycjIWLFiAHj16wNfXF9evX8eKFSuQlZWFN99809bNIxtgiCWyExqNBnPmzEFycjJkMhkiIiLw7bff4oUXXrB104iIbE4ul+PatWt4/fXXcffuXbi5uaFTp05YvXq1fipDql1YTkBERERE1Q4XOyAiIiKiaochloiIiIiqHYZYIiIiIqp2avzALq1Wi1u3bsHT09PsMoVEREREZB+EEMjKykJQUFCpC1nU+BB769YtBAcH27oZRERERFQON27cMLvSJVALQqynpycA6UR4eXnZuDVEREREVJrMzEwEBwfrM5w5NT7E6koIvLy8GGKJiIiIqomyykA5sIuIiIiIqh2GWCIiIiKqdhhiiYiIiKjaqfE1sZbSaDQoKCiwdTPISlxcXEqdloOIiIiqt1ofYoUQSE5ORnp6uq2bQlbk4OCAsLAwuLi42LopREREVAlqfYjVBVh/f3+4ublxQYQaQLfARVJSEho2bMifKRERUQ1Uq0OsRqPRB1hfX19bN4esyM/PD7du3UJhYSGcnZ1t3RwiIiKyslpdNKirgXVzc7NxS8jadGUEGo3Gxi0hIiKiylCrQ6wOv26uefgzJSIiqtkYYomIiIjIUL4KUKWZfkyVJj1uYwyxVC0dOnQIMpmMs0oQERFZW74KiF4IHIgyDrKqNGl79EKbB1mGWCIiIiIqUqAG8jKB7BTDIKsLsNkp0uMFalu2kiHWGjRagdgrafg5/iZir6RBoxW2blK1IIRAYWGhrZtBRERExbn7Ar2iAA//oiCbeqEowHr4S4+723ZmJ4bYCtp9JgldlhzEc18exZub4vHcl0fRZclB7D6TVKmv2717d7zxxhuYPHkyvL29Ua9ePaxZswYqlQovvvgiPD090bhxY/zvf//TH3Pu3DkMHDgQHh4eqFevHkaOHIk7d+4UvZfdu9GlSxfUqVMHvr6+eOyxx3DlyhX94/n5+Zg4cSICAwOhUCgQGhqKRYsWAQCuXbsGmUyG+Ph4/f7p6emQyWQ4dOgQgKISgD179qBDhw6Qy+X47bffIITA0qVL0ahRI7i6uqJNmzb48ccfDd7vrl270KxZM7i6uqJHjx64du2a9U8qERERSUoG2X1z7CrAAgyxFbL7TBLGb4hDUoZhd3pyhhrjN8RVepD95ptvULduXRw/fhxvvPEGxo8fj2eeeQadO3dGXFwc+vXrh5EjRyInJwdJSUno1q0b2rZti5MnT2L37t24ffs2nn32Wf3zqVQqTJ06FSdOnMCBAwfg4OCAJ554AlqtFgDw6aefYvv27fjhhx9w4cIFbNiwAaGhoeVu9/Tp07Fo0SKcP38erVu3xuzZs7F27VqsWrUKZ8+exZQpU/DCCy8gJiYGAHDjxg08+eSTGDhwIOLj4/Hyyy9j5syZVjmHREREZIa7LxA50XBb5ES7CLAAIBNC1OjvvjMzM6FUKpGRkQEvLy+Dx9RqNRISEhAWFgaFQlGu59VoBbosOWgUYHVkAAKUCvw+oyccHaw/3VP37t2h0Wjw22+/Se3RaKBUKvHkk09i/fr1AKTVyAIDAxEbG4tdu3bh2LFj2LNnj/45/vnnHwQHB+PChQto1qyZ0WukpqbC398ff/31F1q2bIlJkybh7Nmz2L9/v9EUVteuXUNYWBj++OMPtG3bFoDUE+vt7Y3o6Gh0794dhw4dQo8ePbBt2zYMGTIEgBSc69ati4MHDyIyMlL/fC+//DJycnLw3Xff4Z133sG2bdtw9uxZ/evOnDkTS5Yswb1791CnTh2jtlfkZ0tEREQwrIHVqYKe2NKyW3Hsib1PxxPumg2wACAAJGWocTzhbqW1oXXr1vq/Ozo6wtfXF61atdJvq1evHgAgJSUFp06dQnR0NDw8PPS35s2bA4C+ZODKlSsYMWIEGjVqBC8vL4SFhQEAEhMTAQBjxoxBfHw8wsPDMWnSJOzdu/e+2t2hQwf938+dOwe1Wo0+ffoYtG39+vX6dp0/fx6dOnUyCM7FAy8RERFZWfEA6+EP9JlvWCNrbvqtKlSrl52tiJQsy0bkWbrf/Si5nKpMJjPYpgt9Wq0WWq0WgwcPxpIlS4yeJzAwEAAwePBgBAcH48svv0RQUBC0Wi1atmyJ/Px8AED79u2RkJCA//3vf9i/fz+effZZ9O7dGz/++CMcHKT/DxXv2NetiFaSu7u7/u+6UoWdO3eifv36BvvJ5XKj5yQiIqJKVjLA6npee0UVbT8QZfPaWIbY++TvadlX1JbuV9nat2+PLVu2IDQ0FE5Oxj/2tLQ0nD9/Hl988QW6du0KAPj999+N9vPy8sKwYcMwbNgwPP300+jfvz/u3r0LPz8/AEBSUhLatWsHAAaDvMyJiIiAXC5HYmIiunXrZnafbdu2GWw7evRomc9NRERE98FZAcj//Rq/eFAtHmTlXtJ+NsQQe586hvkgUKlAcoYapvoJdTWxHcN8qrppJk2YMAFffvklnnvuOUybNg1169bF5cuXsWnTJnz55Zfw9vaGr68v1qxZg8DAQCQmJhoNnlqxYgUCAwPRtm1bODg44L///S8CAgJQp04dODg4oFOnTli8eDFCQ0Nx584dzJ49u8x2eXp64u2338aUKVOg1WrRpUsXZGZm4siRI/Dw8MDo0aPx2muv4aOPPsLUqVMxbtw4nDp1CuvWraukM0VERFTLubgDPd6R5oEt2dOqC7LOCmk/G2JN7H1ydJBh7uAIAFJgLU53f+7giEoZ1HU/goKCcPjwYWg0GvTr1w8tW7bEm2++CaVSCQcHBzg4OGDTpk04deoUWrZsiSlTpmDZsmUGz+Hh4YElS5agQ4cOeOihh3Dt2jXs2rVLX0rw9ddfo6CgAB06dMCbb76JDz74wKK2vf/++5gzZw4WLVqEFi1aoF+/ftixY4e+Jrdhw4bYsmULduzYgTZt2mD16tVYuHChdU8QERERFXFxN18q4O5r8wALcHaCCo9g330mCfN2nDMY5BWoVGDu4Aj0bxlYobbT/ePsBERERNWTpbMTsJyggvq3DESfiAAcT7iLlCw1/D2lEgJ76YElIiIiqokYYq3A0UGGyMb2MfEvERERUW3AmlgiIiIiqnYYYomIiIio2mGIJSIiIqJqhyGWiIiIiKodhlgiIiIiqnYYYomIiMi68lWAKs30Y6o06XGiCrJpiF20aBEeeugheHp6wt/fH0OHDsWFCxcM9hFCICoqCkFBQXB1dUX37t1x9uxZG7WYiIiISpWvAqIXAgeijIOsKk3aHr2QQZYqzKYhNiYmBhMmTMDRo0exb98+FBYWom/fvlCpii7spUuXYvny5Vi5ciVOnDiBgIAA9OnTB1lZWTZsue11794dkydPtnUzrKomviciolqnQA3kZQLZKYZBVhdgs1OkxwvUpT0LUZlsutjB7t27De6vXbsW/v7+OHXqFB599FEIIfDxxx/j3XffxZNPPgkA+Oabb1CvXj189913GDdunC2abRd++uknODs727oZ6N69O9q2bYuPP/7Y1k0hIiJ74O4L9IoqCqwHooDIiUDsSum+h7/0uDsXCaKKsaua2IyMDACAj48PACAhIQHJycno27evfh+5XI5u3brhyJEjJp8jLy8PmZmZBrdKY8OaHx8fH3h6elba8xMRUS1ljd9tuiDr4S8F131zGGDJ6uwmxAohMHXqVHTp0gUtW7YEACQnJwMA6tWrZ7BvvXr19I+VtGjRIiiVSv0tODi4chps45qf4l+9h4aG4oMPPsCoUaPg4eGBkJAQ/Pzzz0hNTcWQIUPg4eGBVq1a4eTJk/rj161bhzp16mDbtm1o1qwZFAoF+vTpgxs3buj3GTNmDIYOHWrwupMnT0b37t31j8fExOCTTz6BTCaDTCbDtWvXAADnzp3DwIED4eHhgXr16mHkyJG4c+eO/nlUKpW+vYGBgfjoo48q5TwREVE5WPN3m7uv1ANbXOREBliyGrsJsRMnTsTp06fx/fffGz0mk8kM7gshjLbpzJo1CxkZGfpb8VBmVXZW87NixQo88sgj+OOPPzBo0CCMHDkSo0aNwgsvvIC4uDg0adIEo0aNghBCf0xOTg4WLFiAb775BocPH0ZmZiaGDx9u8Wt+8skniIyMxCuvvIKkpCQkJSUhODgYSUlJ6NatG9q2bYuTJ09i9+7duH37Np599ln9sdOmTUN0dDS2bt2KvXv34tChQzh16pRVzwkREZWTNX+3qdKkEoLiYlea7+UlKie7CLFvvPEGtm/fjujoaDRo0EC/PSAgAACMel1TUlKMemd15HI5vLy8DG6VouRXJQeigNQLRf/Iq/grk4EDB2LcuHFo2rQp5syZg6ysLDz00EN45pln0KxZM8yYMQPnz5/H7du39ccUFBRg5cqViIyMxIMPPohvvvkGR44cwfHjxy16TaVSCRcXF7i5uSEgIAABAQFwdHTEqlWr0L59eyxcuBDNmzdHu3bt8PXXXyM6OhoXL15EdnY2/u///g8ffvgh+vTpg1atWuGbb76BRqOprNNDRESWsNbvtuKh18Mf6DPf8DkZZMkKbBpihRCYOHEifvrpJxw8eBBhYWEGj4eFhSEgIAD79u3Tb8vPz0dMTAw6d+5c1c01Zkc1P61bt9b/XRfwW7VqZbQtJSVFv83JyQkdOnTQ32/evDnq1KmD8+fPV6gtp06dQnR0NDw8PPS35s2bAwCuXLmCK1euID8/H5GRkfpjfHx8EB4eXqHXJSIiVLymtaK/20oG2F5RgF+4cThmkKUKsmmInTBhAjZs2IDvvvsOnp6eSE5ORnJyMnJzcwFIZQSTJ0/GwoULsXXrVpw5cwZjxoyBm5sbRowYYcumF7GTmp/iMxXoSi1MbdNqtQbHmSrL0G1zcHAwKD8ApN7bsmi1WgwePBjx8fEGt0uXLulnnSAiokpgrZrWivxuc1YAci/j0Fs8HMu9pP2IKsCmIXbVqlXIyMhA9+7dERgYqL9t3rxZv8/06dMxefJkvP766+jQoQNu3ryJvXv32s/I/Gpc81NYWGgw2OvChQtIT0/X95r6+fkhKSnJ4Jj4+HiD+y4uLkZlAO3bt8fZs2cRGhqKJk2aGNzc3d3RpEkTODs74+jRo/pj7t27h4sXL1r5HRIR1TLWqmmtyO82F3egxzume211QbbHO9J+RBVg83ICU7cxY8bo95HJZIiKikJSUhLUajViYmL0sxfYXDWv+XF2dsYbb7yBY8eOIS4uDi+++CI6deqEjh07AgB69uyJkydPYv369bh06RLmzp2LM2fOGDxHaGgojh07hmvXruHOnTvQarWYMGEC7t69i+eeew7Hjx/H1atXsXfvXowdOxYajQYeHh546aWXMG3aNBw4cEDfw+7gYBcl2kRE1Zc1alqt8bvNxd38a7j7MsCSVTA13K8aUPPj5uaGGTNmYMSIEYiMjISrqys2bdqkf7xfv3547733MH36dDz00EPIysrCqFGjDJ7j7bffhqOjIyIiIuDn54fExEQEBQXh8OHD0Gg06NevH1q2bIk333wTSqVSH1SXLVuGRx99FI8//jh69+6NLl264MEHH6zS909EZBesPed4RWpaa8DvNqo9ZKKGFyhmZmZCqVQiIyPDaKYCtVqNhIQEhIWFQaEoZ22Oru4oL9P4Q0H3ISD3stuvTNatW4fJkycjPT3d1k2pFBX62RIRVZXK/F2SekEKsDp95kuB1FbtIbJQadmtOJsuO1ut6Wp+CtTma36cFfxHTkREknyV8e8MXQ1rxj/AvtlAnw+kx4v3iOr2K8/vE3M1rWX1xPJ3G1UjLCeoCNb8EBGRJczNGuDuC0ROAtITgetHpCBb0TnHK1rTyt9tVE0wxNZSY8aMqbGlBEREdqe0WQNiPwXqNAQcnICs5IrNOc6aVqpFGGKJiIgqW1mzBigbAEM+Bxxdio65nznHOUcr1SIc2MXBPzUSf7ZEZJdK1roCUrCMnCT1yJbcfj+rP5qqvS3++qxpJTtn6cAu9sTCeBUrqv5q+P/NiKi6MrUSVtsXigKsNeYcZ00r1RK1enYCFxcXODg44NatW/Dz84OLi4vJZVipehFCIDU1FTKZzGDpXSIimys5a4AmH/j5dakmVtmgqOe1V1RRj+2BqPvrkSWq4Wp1iHVwcEBYWBiSkpJw69YtWzeHrEgmk6FBgwZwdHS0dVOIiCQlB11FTgR+XwHcvSrNTtB/qXENq25eVtawEhmp1TWxOkIIFBYWQqPRVHHrqLI4OzszwBKR/TA1a4BuPti97wLZtw17YosfxxpWqmW42EE56L525lfPREQ1lK0HO+lmDQCMZw3ou8B8jytLCIjMYk8sERHVbPaylKqtgzRRNcHZCYiIqHbIV5kfwa9KA3LSzC80oPuKPy9TCpiVibMGEFkVQywREVVf5pZzBYpC6rEvgEenmV9o4H7nYyUim2KIJSKi6qV4z2vJ5VxTLxU9XryX1dndcMWsiiztSkR2gSGWiIjsk6kyAV3P656ZwL0bhsupZvwDbBoO7Hwb2DfbOKSaWmjgfpZ2JSK7wBBLRET2x1yZQIEaUKUA144A3z9bFGQjJ0lzraozgIRfgaxk417WkgsNANL98q6IRUR2gSGWiIjsgyVlAgCAf1dWVGcC0fOl+tbYTwGPAMBRDvg1AxxdDHtZS87TWtGlXYnI5jjFFhER2Z6pabB0wTPjH6mXNfhhQGiB3HuAQglASEEWkJZvTU+Ulm91dJG26XpiAfMLDXBwF5Hd4RRbRERUfZTseVWlGZcJXI2RVrby8Af6LQK6vi0dq8kHUi8CngHSqlcle1nzc6R5YEsG1eL1tFzalajaYU8sERFVLksn+S/ZMxo5UapZvXsVuHutqEygz3zArW5RL23qRUCTJ/XODt8E+DU1fq6u0wEXNy40QFQNWJrdGGKJiKjylHe1rOLhEzBdJlC8lMDVG5A5ADeOSfsoGxiXC1TFalxEZDUsJyAiItszVSYAmF8tq/g0WLoyAY9iZQIKJXDtsDQ7gcIL6LsAGLhM6oFVNjAuR+gVxQBLVEMxxBIRUeUpXndqyWpZummwdAFWkwdkJ0u1sX7hQM/ZUngFoJ+lwMVdKiEwVd/K5VyJaiyWExARUeUrWSYAmA6wun3cfAAB4J/jxmUC6YnAwQ8Adz/jXlbWtxJVe6yJ/RdDLBGRnUi9IC33qtNnvtS7Cpie7spZAWTckuaANTU9FsMqUY3EmlgiIrIfZa2W5awwngaLZQJEVAr2xBIRUeUyN3VWyd5VS6fiIqIajT2xRERkmeLLvZakSiu23Ot9MFUm4BduPNhLlSYFVHMrZrHnlYhKYIglIqrNslOBve8ZTn+lk3oJ2PuuNM/r/QZZU2UCAFfLIqIKc7J1A4iIyEbyVcChxUBCjDQDwIGooqCZegnYNBzQFgINI6Wv+e+nJ9TFXZpBwFSZgC7IskyAiO4De2KJiGqrAjWgzZcCbHqitITrgSgg8ZgUYNUZgIMT0OUt81/zW4JlAkRUCRhiiYhqK11PqLJBUZBNuwr8OFYKsAqltBKWX1Nbt5SIyIhNQ+yvv/6KwYMHIygoCDKZDNu2bTN4PDs7GxMnTkSDBg3g6uqKFi1aYNWqVbZpLBFRTVQ8yHoGAHcuSKtkOcqBIZ8zwBKR3bJpiFWpVGjTpg1Wrlxp8vEpU6Zg9+7d2LBhA86fP48pU6bgjTfewM8//1zFLSUiqsHcfYG2LwB3rxVt8w4F4jeYn7WAiMjGbBpiBwwYgA8++ABPPvmkycdjY2MxevRodO/eHaGhoXj11VfRpk0bnDx5sopbSkRUg6VeAn5+vagH1i8cyE4uqpFlkCUiO2TXNbFdunTB9u3bcfPmTQghEB0djYsXL6Jfv35mj8nLy0NmZqbBjYiIzNDNQqCrgX36a8CnkfFgLwZZIrIzdh1iP/30U0RERKBBgwZwcXFB//798fnnn6NLly5mj1m0aBGUSqX+FhwcXIUtJiK6PxqtQOyVNPwcfxOxV9Kg0VbBYoqqNODwcmkaLd0groYPGw/2cnDmPK5EZHfsep7YTz/9FEePHsX27dsREhKCX3/9Fa+//joCAwPRu3dvk8fMmjULU6dO1d/PzMxkkCUiu7b7TBLm7TiHpAy1flugUoH5A0LRp2mdyluG1VkBuPoAIZ2BR6YWDeLSDfY6EAV4NwK6z+Q0WERkd2RCiCr4737ZZDIZtm7diqFDhwIAcnNzoVQqsXXrVgwaNEi/38svv4x//vkHu3fvtuh5LV1/l4jIFnafScL4DXEo/kHsCjW8kYUXnfaif1M3BD+9pCjIqtKAAhXw6zJppase71QsYOarTC9EoHstLkRARFXM0uxmtz2xBQUFKCgogIODYcWDo6MjtFqtjVpFRGQ9Gq3AvB3nICAFVwXyoYYLpjhtgS/SIQNw/fpt1N89Ew693gOc3YE9M4HkM4BXEOCF+19JS8fF3fzxFVnggIioktk0xGZnZ+Py5cv6+wkJCYiPj4ePjw8aNmyIbt26Ydq0aXB1dUVISAhiYmKwfv16LF++3IatJiKyjuMJd5GUoYYr1JjitAWeyMGawkHwRA6UshxkCDek5zui4MrvkCc9Dfg2BpJOSwcHtCpaIpaIqBayaYg9efIkevToob+vq2UdPXo01q1bh02bNmHWrFl4/vnncffuXYSEhGDBggV47bXXbNVkIiKrScmSamAVyIcncuAnS8erTjuxpnAQXnXaiSDZHQTLUqDVaICMG0DOHcDJFQh9BOi3iAGWiGo1u6mJrSysiSUiu2Ci9jT2Shqe+/IovJEJBfIwyWkb/GTpSBV18F/No1jqvAZKmQperi6QOwgpwPo1A/ovluZyJSKqgSzNbnY9xRYRUY2QrwKiFxrNt9oxzAfNvfIxw2kzXnTai08LhyJV1IGfLB2vO21HqqgDrcwJLs7OgMwB8AkFHF2A2JWct5WIaj2GWCKiylagBvIygewUgyDrmHsXX9bfCT9ZOryQAzXk+EozEADggkKEypKhlDtApi2UVtJy9ZHmcy3xPEREtRFDLBGRteSrTAdLd18gchLg5lMUQFMvAAeiEOychWZNmmCd22gAwMuOu+CMQoQ73oS/Uw4U2hygfgcgtDOgyQcgGGSJiGDHU2wREVUrupKBvEzjWQNUaUDsp1JNqy7I7psjPebhj+BeUdgpBO5seweyjHx45KngqmwGWerfgLMroPAEuk6XniM7BVB4SUFW7sWVtIio1mJPLBFRRWSnAqmXTJcMpF4CUv+W7menANoCoP0Yw+MjJwIAHA/OQz2HDPjXD4Fbsx6Q+TUDRu0AwroC6kwpwEZOAjz8AXd/oOfsii90QERUjXF2AiKi+5WdCmx8GsjPBoZvknpZdYHVSQHcigPysoHANkCdYCmE6npTdTz8gUenAce+KOrFdVYUzWSgSpOeU7c6V4Gaq2gRUY1maXZjiCUiul+pl4BNwwF1hvT1vi7I7pwKXNon9bw6OANN+wA9ZhcFWA9/qQc2dmXR/UenSStycflXIqrlOMUWEVFl82sqBVeFUgqym4YD/5yQemB1AbZ+O6DT64YBtleUNM9rryjpfnYK8Osy86/j7ssAS0RUAkMsEZEp5mYaAKTt+Srp7yWD7M+vS6FUF2BdPIBTawEHl6IAq+ttdfctCrIcpEVEVC4MsUREJZlZnABAUY1q9ELDINt3ASC00jahBRo+DAxaLgXU3HtAgUqqiS1ZLqALshykRURULgyxRESAYc9ryZkGUi8VPa4buJWXKe0HSI/vmVkUYGUOwL0EAA5FPa2596SSAnPzyDLAEhGVC0MsEVHJntfiX/Nn/CPVuu6aBux917Cu1d232OCuLGlGgtAugFegNCvBpuFAzl2WDBARVQLOTkBElUajFTiecBcpWWr4eyrQMcwHjg4yWzfLWPEeVpMBNUNa9tWvGaBsYPpxhRJ4ai2grC8FV1OzFnCGASKiMlma3bhiFxFVit1nkjBvxzkkZaj12wKVCswdHIH+LQNt06h8FaC6C7i4Gdam6paFjVkEZN6SAq1uCizPAKlswK8Z4Ogibdcd61pHGrgFSEHVr2nR8w3fJAVZFw9pP1NTZxER0X1jTywRWd3uM0kYvyEOJT9cdH2wq15oXyVBtnhPcIBCi4eufAaH678BAa2AfouKgqUqTappvRUv1bMqGwAOToAmH0hPBOo0lAIsYDzDQHYqkJteFGCLS70kBVgPv0p/r0RENQV7YonIJjRagXk7zhkFWAD6bTN/+guecmd0auxbaeUFJXuCvZGJha7n0V1xD67qw8CeWVKQBaQAe+2I9PeAVoAmTwqwqRcB71Ap1BZfnOBAVFGQ9fAzH1JNBVsiIrIK9sQSkVXFXknDc18etWjfyiovMNcT7INMzHb6FgM8r8LV2RFo8CAgBHAzTtqhQQfASQ6oUqUAq8krqmn1a2q+dpaIiKyGK3YRkU2kZKnL3ulfyRlqjN8Qh91nkqz2+qX1BN+FFz4oHIlD6qYQEMD1I0BirPSgLsCqM6Q62LBHpQBbp2HR1FhcnICIyG4wxBKRVfl7Wh7sdEFz3o5z0Git86XQ8YS7BoPJSroLL3yoHowMRYOijcr6UmvUGVJA7fMBMOhDqQdW2aCohKB4kOXiBERENsUQS0RW1THMB4FKBSytdBUAkjLUOJ5w1yqvX1ZPsDcy8brjNrhkJRZtTL8BpPwt9bzqSgRc3KUSAlM9rxVYnECjFYi9koaf428i9kqa1cI7EVFtw4FdRGRVjg4yzB0cgfEb4iADTH6tb0p5yhAASNNlFaiNalL9PRXwRibUcEEuDHuFvf+tie3k8DccZXIgpHNRTWxeFlCYZ/w6up5XK8zxapfTjhERVVPsiSUiq+vfMhCrXmiPAKXlpQUWlyHkq6Rpr4qvsKWjSkNHbxWi3H7EVKctcEVRWDQIsA4yODfuAgxaATz2MRDaWZpa65+T0qwFJZeGtcKysLrBZiVLHSqjLpiIqDZgiCWiStG/ZSB+n9ETG19+GHVcnc3uJ4PUG9kxzKf0J8xXAfduSOH14AeAKqWoVjX1khRs98yC46ZhiPTLgydy4Ip8/eFquCATHsiCGzQNO8Oh/2IpnLr7Av0WS0FW4SXN62rlAVuWTDtmzbpgIqLagOUERFRpHB1keKRJXSx+qhXGb5CmsSoe03R1s3MHR5Q+X+zd60DMUkDkS0fplnNVeAGpF4ANTwKOzkBhLgAH+Ie2hF/HKZDvTQL+7fnMhQLr3UYjtFcAmrcMNV6xq99ioEAFuFW817WksgabFa8LjmzMKbuIiCzBEEtElU5XXlCyHjSgtHpQXc2r0AKbR0g9rYGtAa8gwEkB5KRJ9az/nAKERgqxbr5AWGeg32L0cvdF93Yt9Ct2+XtKvb1mw7K7L4DKCZCW1vuWuy6YiKgWY4gloirRv2Ug+kQEWBYq714HYpYAjo5A21FAoRpwdIZIOg1Vxj1AfQ+ODjIocpIhE1rpGIUS8G8BdH1b38vq6CCzi55NS+t9yzM9GRFRbccQS0RVpsxQma8C0v8BtrwEpF+Xel7j1wNDViN781ho81LhoT4L7b/zHggZAMggU3gBdcMBRxdpaVg7W0lLN+1YcobaZF2sDFKvdJl1wUREpMeBXURkW9mp0sCsfJU0aOvQQiA/WyoPuBUPpPyNm/s/xblMBTxELmQAHCGkeloBZApX5MrcpNCq8DJcmMBO6KYdA2A0f67FdcFERGSAIZaIbCc7Fdj4NLBpOJB8DsjLlEoHAtsATq5AQQ7EPyfgnngQ7R0u6w/TQvrwKoAjMoQbjqhDIQrUAGR2G2TNTTsWoFRg1QvtOU8sEVE5sZyAiKpWdiqQmy6thpWbLvW6qjOAn18Dur0LXPofkJUE5N77dyoDAS+RDQEZZJACLCCD5t/7ycIb6QUOuFPoCr/isxYUX2HLTpSrLpiIiErFEEtEVUfX85qfDQzfJAXZ4ZuknticNGDbq9LsAzlpgLYAcHQCNAUQABwg9AH2lLYp/GXp8JLlINzhHxRonZBb4AN4KgF3P6Dja4C7j9WnyrIGexlsRkRU3bGcgIgqX75K+mq/eM/rpuFSLaybD9BrHpCXDWgLgYx/pNWz4ABoNQAA2b/DoWQAzmkb4B/hhzfzJ+Cu8EQhnNDS4RrkcjnQczbQ4x3AO9guAywREVmPTAhRo5eIyczMhFKpREZGBry8vGzdHCK7pdGKyvmaWzdgKy9TmjUg564UYNUZgIsH4BMG3IyTAmtBDiD3kP7UagAhIGQOyBHOcBb5yIUCgMAZbSgSRQA2FXbHcpfV0DjI0fiNn+HoG1rx9hIRkU1Zmt1s2hP766+/YvDgwQgKCoJMJsO2bduM9jl//jwef/xxKJVKeHp6olOnTkhMTKz6xhLVYLvPJKHLkoN47sujeHNTPJ778ii6LDmI3WeSKv7kBWopwOoGW7n5SCUELh5S7eu136XQWpD7b4DNBTQF0iIHzq6QNe2LHP8HcVjT0qDnVQYtLiEYr+S/hRuPbWSAJSKqZWwaYlUqFdq0aYOVK1eafPzKlSvo0qULmjdvjkOHDuHPP//Ee++9B4XCvgZrEFVnu88kYfyGOKNlUZMz1Bi/Ia7iQdbdV+qB9fAvCrLQSj2wuoUKNIXSICwhAL/m/5YTyACFN9DlLfg1aIwHgn1w1bERMoQb/hF++LTwCSiVdTD9hcfQq0PLirWRiIiqHbspJ5DJZNi6dSuGDh2q3zZ8+HA4Ozvj22+/ve/nZTkBkXkarUCXJQeNAqyObhL+32f0rHhpgSpNCrDZKVKJwc1TUg2sphCAAGSOQPhAwK0OED4Y2PUWUJgrzTYwZDUQvx5aFy/85dMPt7R1UMevPkf2ExHVQNWinKA0Wq0WO3fuRLNmzdCvXz/4+/vj4YcfNllyQET353jCXbMBFpBmuErKUON4wt2Kv5i7LxA50TDAuvsB3iGAgxOg8ARSzgBtR0HTtA/+6PY1smQeyBau0NQJAXpFwaHnO2jT8VEM6NQakY19GWCJiGoxuw2xKSkpyM7OxuLFi9G/f3/s3bsXTzzxBJ588knExMSYPS4vLw+ZmZkGNyIyLSXLfIC9n/1KpUqTBnjpAqyDE9DgIeC5zVKNrJsvkJeF7B9ewXOLvsUTP6Tg8TsT8Oit19HlP39hd0I+ZxwgIiI9uw2xWq1UKzdkyBBMmTIFbdu2xcyZM/HYY49h9erVZo9btGgRlEql/hYcHFxVTSaqdvw9Lasvv5OVB422WOWRbsosSCUJsVfS8HP8TcReSYMm+470eHG6UgJ1BuCkADwCgKZ9pdW5jq4E6rcDhm9CtswDVzOBy1nSFNYJCMJdKK1Xn0tERDWG3S52ULduXTg5OSEiIsJge4sWLfD777+bPW7WrFmYOnWq/n5mZiaDLJEZHcN8EKhUIDlDjdKK49/feR5f/Z6A+QNC0SfEBTi+GsjLxIH64zF7bxKSMtTwRiYUyMNMtx1o3TgYYc8skHpOi9fC1gkG+i+RBm65+RRtPxAFTc+5GJs3FZfznXAXSoPXF5Dqc+ftOIc+EQEsIyAiIvvtiXVxccFDDz2ECxcuGGy/ePEiQkJCzB4nl8vh5eVlcCMi0xwdZJg7WPqPYlmxMCMjHdd+fA83t74DqFKRkpSI1K3vQJ2RAm9kYrbTt/ja5UN4FaTi2PlrOPDXdelAZ4W0BKyHvzRLgX+4tFJX8VkL5F44+U8ujmf5GgVYHavW5xIRUbVn057Y7OxsXL58WX8/ISEB8fHx8PHxQcOGDTFt2jQMGzYMjz76KHr06IHdu3djx44dOHTokO0aTVTD9G8ZiFUvtMe8HedMDvJyhRoK5AMAPJGDhH+yEODVCKduaxEou4cop3VwghZtHa4AAM6LhlhaOAzyvUno3q4FHF3cpVW0CtRScC1OF2SdFUg+l25Re61Sn0tERNWeTXtiT548iXbt2qFdu3YAgKlTp6Jdu3aYM2cOAOCJJ57A6tWrsXTpUrRq1QpfffUVtmzZgi5dutiy2UQ1Tv+Wgfh9Rk+8N6iFwXZXqDHFaQumO20GACwpHIYb+Z64mZyMwnw1wmTJ6OkYj46O5wEAx7TN8X7hSNyFl2GvqYu7cYDVcfcFXNwtrs+1dD8iIqrZbNoT2717d5Q1Te3YsWMxduzYKmoRUe3l6CBDXU+5wTYF8uGJHPjJ0jHdaTOWFg7D0sJhWKv5Ce0dzsJVlgcnaFAoXHFZBOE/mqG4h6ISnvL0mpZVn6ubs7ZjmM99vkMiIqpJ7LYmloiqXoBCi0CkwhvS1HT34IWlhcOQKuogSHYHUU7fIEh2B17OWihk+ZCjACrIcU0EoABOeNlxl/5YoHy9pqXV5+ruzx0cwUFdREQEgCGWiHTyVXjoymf4VvER3nP61iDIrikchIayFPRy/ANfyT9GWM4ZuMrykQtn5Ak57gpPZAo3fY+tDzIReB+9prr63AClYfgNUCqw6oX26N8y0Gpvl4iIqje7nWKLiKpYgRoOefdQ37UArqq/8Z7Tt3i/cCQAYKLTVrhBDSdo4OeQDVmugIPcE8dzGkEtnOEly4FaOCNTuMH/3yDr13fhffWa9m8ZiD4RATiecBcpWWr4eyq4vCwRERlhiCUiibsv0G8xXDETPpd+Q+fcC/hA9jUggPYOl+DpoIaToyOcXNwATT4U7l7w6DYPnx66jdE538BPlo4M4YZ8Z0883DgUYa3MT4VXFkcHGSIbmxkIRkREBEAmyhpZVc1lZmZCqVQiIyODc8YSWUKVBuyZCXHtMAry1dBqtHDW5sHB1QsyhRLwCwccnAH1PUDZAJqecxF3/R78ji+Gi3sd1Ov9Jhw9fLlELBER3RdLsxt7Yolqm3yV6TlbASnAOiuArm9DlnELLqnngcJcwMUVaPgw0PcDwPnfcPrvaluOB+fhoV5RQOhS6ViGVyIiqgIMsUS1Sb4KiF4I5GVKiwwUD7K65WEdXAD1XeDeNUjzAsiAwnzpT+di8732ipL2l3sxvBIRUZVjOQFRbaDrfQX0Paj6ZWABoEAF/LoMyPgHuJcA5KsBB0egwYOApgC4+Yd0P/QRoN+ioiCr67llgCUiIiuxNLsxxBLVdCV7X4GiIKtQAgW5QNplwN0PyEgsCrC6wAoAe2YC145Ify8ZZImIiKzI0uzGeWKJaroCtRRgs1Ok8ApIYVahBK4dBhJjAXUm4F4XCOkKuHkbBtV/Zy1AaGdA4QW41pF6X4mIiGyINbFENZ27b1H9qi7IRk6UemA1eYCjHPBrBnSbASgbAKq7gIubYU+rLsgWqAA3zjxARES2V65ygoyMDGzduhW//fYbrl27hpycHPj5+aFdu3bo168fOnfuXJltvS8sJyD6l27gVnYKoMkHUi9K2/2aAY4uRTWyLBMgIiIbsmo5QVJSEl555RUEBgZi/vz5UKlUaNu2LXr16oUGDRogOjoaffr0QUREBDZv3my1N0FE5ZCvkoKqKboBWJETiwKsJk8auNV/sRRgdb205p6DiIjIjlhUTtCmTRuMGjUKx48fR8uWLU3uk5ubi23btmH58uW4ceMG3n77bas2lIhKYdHUWf8uUKALsI5ywEkOuNU1LjdgjywREdk5i8oJUlNT4efnZ/GTlnf/ysRyAqoVipcKFC8L0G3P+AdITwS86gOqVKBuUynAqjMNp9rSzfva4x3WvRIRkU1wiq1/McRSrVEyyEZOBGJXFgXYOg2lgVtdp0sDtwDTc8Zy3lciIrKhSg2xN2/exOHDh5GSkgKtVmvw2KRJk8rf2krEEEs1kUYrcDzhLlKy1PD3VKBjmA8cHWSGQVbH1VtaaUubb77UgL2vRERkJyotxK5duxavvfYaXFxc4OvrC5lMVvRkMhmuXr16/62uBAyxVNPsPpOEeTvOISlDrd8WqFRg7uAI9G8ZCKReAPbNKTqgz3ypB7ZAbbrOlatuERGRHam0EBscHIzXXnsNs2bNgoOD/a+VwBBLNcnuM0kYvyEOJf/R6v4r+dUzYeh1c5VhTyynziIiomqk0lbsysnJwfDhw6tFgCWqEf6dOkujFZi345xBgPVGJlyhhvj379k734PIui0F1z7zOXUWERHVWOVOoi+99BL++9//VkZbiKgk3dRZB6IQd/6yQQmBNzIx3WkzpjhtQSBSMc1pM9wL0pAi6kg9r37h0p8MskREVAOVe9nZRYsW4bHHHsPu3bvRqlUrODs7Gzy+fPlyqzWOyJ6ZHVxlTQVqae7X7BT4HV8Mb/TGPXjpA6yfLB26rtksuAECiGs6CQN0pQPFl5yVe0m1r0RERDVAuUPswoULsWfPHoSHhwOA0cAuotqgzMFV1lIshHqp/sF0p834SjMQLzvugp8sHamiDpYWDsM9eGFF4VNQIB+f1w00/RylDN6qkkBORERkReUe2OXt7Y0VK1ZgzJgxldQk6+LALrK2sgZXrXqhvXWDLACo0qDdH4Ujf55DXoEGAjAIsLrXD1Aq8PuMnuUKoFUWyImIiCxQaQO75HI5HnnkkQo1jqi6MjW4Ske3bd6Oc9BorbyGiLsvHDpPRLN6HgCkwPqVZqBBgAWAuYMjyh1gx2+IMwiwAJCcocb4DXHYfSbJGq0nIiKyunKH2DfffBOfffZZZbSFyO4dT7hrFPiKEwCSMtQ4nnDXui+sSgNiV8LfU4FWDZSQOzviZcdd8EYmAKkHtrw9wDYL5ERERFZQ7prY48eP4+DBg/jll1/wwAMPGA3s+umnn6zWOCJ7k5JlPsDez34WKbGcrH+fiah7ZCXSU/5BpNN+pHacifYtmpS7hrU8gTyyMeeYJSIi+1LuEFunTh08+eSTldEWIrvn72nZ6H5L9ytTiQCrW7TAoXcUfA5EwSc7BaFXPgVCo8q9mIFNAjkREZGVlDvErl27tjLaQVQtdAzzQaBSgeQMtcmv4XWDqzqG+VjnBZ0V0tRYgOGqW1aYOqvKAzkREZEVlTvEEtVmjg4yzB0cgfEb4iADDILs/Q6uKpWLO9DjHWm+2JI9rRZMnVWaKg/kREREVmTRwK7+/fvjyJEjZe6XlZWFJUuW4D//+U+FG0Zkr/q3DMSqF9ojQGnYQ3k/g6ss4uJuvlTA3fe+AixQFMiBogCuUymBnIiIyIos6ol95pln8Oyzz8LT0xOPP/44OnTogKCgICgUCty7dw/nzp3D77//jl27duGxxx7DsmXLKrvdRDbVv2Ug+kQEVPsFAnSBvOQ8sQGcJ5aIiOycxYsd5Ofn48cff8TmzZvx22+/IT09XXoCmQwRERHo168fXnnlFf1KXvaCix0QlY0rdhERkb2wNLuVe8UunYyMDOTm5sLX19domi17whBLNpOvMl3LCkizDtxnLSsREVFNVmkrdukolUoEBARUKMD++uuvGDx4MIKCgiCTybBt2zaz+44bNw4ymQwff/zxfb8eUZXJVwHRC6XZA1Rpho/pps2KXijtR0REROV23yHWGlQqFdq0aYOVK1eWut+2bdtw7NgxBAUFVVHLiCqoQA3kZUrzuxYPssXnfc3LlPYjIiKicrPpFFsDBgzAgAEDSt3n5s2bmDhxIvbs2YNBgwZVUcuIKqj4PK66IBs5EYhdabRwAREREZWfTXtiy6LVajFy5EhMmzYNDzzwgK2bQ1Q+uiDr4S8F131zGGCJiIisxK5D7JIlS+Dk5IRJkyZZfExeXh4yMzMNbkQ24+4r9cAWFzmRAZaIiKiCyh1iGzVqhLS0NKPt6enpaNSokVUaBQCnTp3CJ598gnXr1kEms3yqn0WLFkGpVOpvwcHBVmsTUbmp0qQSguJiVxoP9iIiIqJyKXeIvXbtGjQajdH2vLw83Lx50yqNAoDffvsNKSkpaNiwIZycnODk5ITr16/jrbfeQmhoqNnjZs2ahYyMDP3txo0bVmsTUbkUH8Tl4Q/0mV9UWmBq1gIiIiKymMUDu7Zv367/+549e6BUKvX3NRoNDhw4UGq4LK+RI0eid+/eBtv69euHkSNH4sUXXzR7nFwuh1wut1o7iIxYMv9rgdowwOpqYEsO9mJtLBER0X2xOMQOHToUgLRC1+jRow0ec3Z2RmhoKD766KNyvXh2djYuX76sv5+QkID4+Hj4+PigYcOG8PU1/OXu7OyMgIAAu1sVjGoR3fyveZnGAVTX8yr3Ah6ZLP0JGO5XPMjKvaTAS0REROVmcYjVarUAgLCwMJw4cQJ169at8IufPHkSPXr00N+fOnUqAGD06NFYt25dhZ+fyOpKzv+qC6jFSwcAQOYA9HjHdI+tLshyxS4iIqL7dt/LzlYXXHaWrK5krSvnfyUiIrIaS7PbfYXYAwcO4MCBA0hJSdH30Op8/fXX5W9tJWKIpUpRsucVYIAlIiKyAkuzW7lnJ5g3bx769u2LAwcO4M6dO7h3757BjahW4PyvRERENlXuZWdXr16NdevWYeTIkZXRHqLqwdz8r+yJJSIiqhLl7onNz89H586dK6MtRNUD538lIiKyuXKH2JdffhnfffddZbSFyP6VDLC9ogC/cOlPBlkiIqIqY1E5gW7qK0CaamvNmjXYv38/WrduDWdnZ4N9ly9fbt0WEtkTZwXnfyUiIrIDFs1OUHwu11KfTCbDwYMHK9woa+LsBGR1lqzYxflfiYiI7oul2c2intjo6GirNYzIrtxPIHVxNx9SOaiLiIioSpS7JpaoxtAtIWuqhlVX+xq9UNqPiIiI7Eq5p9h64oknIJPJjLbLZDIoFAo0adIEI0aMQHh4uFUaSFRpLF1CtkDN8gAiIiI7U+6eWKVSiYMHDyIuLk4fZv/44w8cPHgQhYWF2Lx5M9q0aYPDhw9bvbFEVqUbjFV8VoHUC8azD7BEgIiIyO6Ue9nZmTNnIjMzEytXroSDg5SBtVot3nzzTXh6emLBggV47bXXcPbsWfz++++V0ujy4MAuKhOXkCUiIrIblma3codYPz8/HD58GM2aNTPYfvHiRXTu3Bl37tzBX3/9ha5duyI9Pf2+Gm9NDLFkkdQLwL45Rff7zJfmfyUiIqIqZWl2K3c5QWFhIf7++2+j7X///Tc0Gg0AQKFQmKybJbJL5paQ5YIFREREdqvcIXbkyJF46aWXsGLFCvz+++84fPgwVqxYgZdeegmjRo0CAMTExOCBBx6wemOJrI5LyBIREVVL5S4n0Gg0WLx4MVauXInbt28DAOrVq4c33ngDM2bMgKOjIxITE+Hg4IAGDRpUSqPLg+UEZJapJWRLzk7A2lgiIqIqVWk1sSVfBIBdh0OGWDJLN09sXqZxUNUFWbkX0OMdTrFFRERURaokxFYHDLFUKi4hS0REZFesuuxs+/btceDAAXh7e6Ndu3alDtqKi4srf2uJbIVLyBIREVVLFoXYIUOGQC6XAwCGDh1ame0hIiIiIioTywmIiIiIyG5U2jyxAJCeno6vvvoKs2bNwt27dwFIZQQ3b968v9YSEREREZWDReUExZ0+fRq9e/eGUqnEtWvX8Morr8DHxwdbt27F9evXsX79+spoJxERERGRXrl7YqdOnYoxY8bg0qVLUCgU+u0DBgzAr7/+atXGERERERGZUu4Qe+LECYwbN85oe/369ZGcnGyVRhFZLF9lflUtVZr0OBEREdU45Q6xCoVCv8hBcRcuXICfn59VGkVkEd1iBaaWh9UtVhC9kEGWiIioBip3iB0yZAjmz5+PgoICAIBMJkNiYiJmzpyJp556yuoNJDKrQC2ttpWdYhhkiy8bm5cp7UdEREQ1SrlD7IcffojU1FT4+/sjNzcX3bp1Q5MmTeDp6YkFCxZURhuJTHP3lZaL9fAvCrKpF4oCrIe/8XKyREREVCPc9zyxBw8eRFxcHLRaLdq3b4/evXtbu21WwXlia4HiPa86DLBERETVklWXnTWlZ8+e6Nmz5/0eTmQVGq3A8WQgx/sZtLu9At5uztKyyJETGWCJiIhqsPsKsQcOHMCBAweQkpICrVZr8NjXX39tlYYRlWX3mSTM23EO6owUTHfajDhZOhTOjmhWzwP+sSvZE0tERFSDlbsmdt68eejbty8OHDiAO3fu4N69ewY3oqqw+0wSxm+I0wdYP1k6UkUdROUOw8EbMqQkJZqetYCIiIhqhHLXxAYGBmLp0qUYOXJkZbXJqlgTW/NotAJdlhw0CrBLC4fhHrzgg0zMdfsRjzd2hMyzHntkiYiIqhFLs1u5e2Lz8/PRuXPnCjWOqCKOJ9xFUoYaarggC24GARYA7sILUTlPI0XUAeRegLOi9CckIiKiaqfcIfbll1/Gd999Z5UX//XXXzF48GAEBQVBJpNh27Zt+scKCgowY8YMtGrVCu7u7ggKCsKoUaNw69Ytq7w2VV8pWdK8r7lQYEXhUwYBVucevBDXdBLQ4x3Axd0WzSQiIqJKZNHArqlTp+r/rtVqsWbNGuzfvx+tW7eGs7Ozwb7Lly+3+MVVKhXatGmDF1980WihhJycHMTFxeG9995DmzZtcO/ePUyePBmPP/44Tp48afFrUM3j71nUs5oLBXJhuqe1Tt1ABlgiIqIayqIQ+8cffxjcb9u2LQDgzJkzBttlMlm5XnzAgAEYMGCAyceUSiX27dtnsO2zzz5Dx44dkZiYiIYNG5brtaiayFdJK2yZqmFVpQHOCnQM80GgUoHkDDVMFXTLAAQopf2IiIioZrIoxEZHR1d2OyySkZEBmUyGOnXqmN0nLy8PeXl5+vuZmZlV0DKyinwVEL1QWiq25GAs3YIGci849ngHcwdHYPyGOMgAgyCr+2/U3MERcHQo33+qiIiIqPood02srajVasycORMjRowodaTaokWLoFQq9bfg4OAqbCVVSIFaCrC6JWR102MVX5ErLxMoUKN/y0CseqE9ApSGpQQBSgVWvdAe/VsGVnnziYiIqOrc97Kz1iaTybB161YMHTrU6LGCggI888wzSExMxKFDh0oNsaZ6YoODgznFVnVRPLB6+Esrb8WuLLpfoodWoxU4nnAXKVlq+HtKJQTsgSUiIqq+Kn3Z2apSUFCAZ599FgkJCTh48GCZQVQul0Mul1dR68jq3H2loKoLsvvmSNsZYImIiKgYuw6xugB76dIlREdHw9eXE9bXCu6+Ug+sLsAC0v1iAVa35GxShlq/LVCpwNzBESwlICIiqgVsWhObnZ2N+Ph4xMfHAwASEhIQHx+PxMREFBYW4umnn8bJkyexceNGaDQaJCcnIzk5Gfn5+bZsNlU2VZpUQlBc7Ep9jaxuydniARYAkjPUGL8hDrvPJFVVS4mIiMhGbFoTe+jQIfTo0cNo++jRoxEVFYWwsDCTx0VHR6N79+4WvQaXna1myqiJ1fSciy6f/WkUYHV002v9PqMnSwuIiIiqoWpRE9u9e3eUlqHtZMwZVRFN9h3c2fYORNZtyDzroW7PuXD0qGtQI3tn2ztQZ/QAYPqiFgCSMtQ4nnAXkY1ZfkJERFRTVZsptqhm230mCb0/PYaf/85G9E0Z+p/pgS6f/SmVBugGe3n4IwuuUMOlzOfTLU1LRERENZNdD+yi2kFX4yoArMBTUCAf9+AF2b81rvp5X3tFIe1GDnLPnC7zOYsvTUtEREQ1D3tiyaY0WoF5O87pV93KhQL3/i0V0G2bt+McNFoBuPuiQ7MGCFQqYK7aVQZplgIuOUtERFSzMcSSTR1PuGt2kBZgWOMKAI4OMswdHAEARkGWS84SERHVHgyxZFOW1q4W349LzhIRERFrYsmmLK1dLblf/5aB6BMRwBW7iIiIaimGWLKpjmE+CFQqkJyhhqkJ1XTzvpqqcXV0kHEaLSIiolqK5QRUObJTgdRLph9LvSQ9Dta4EhER0f1hiCXry04FNj4NbBpuHGRTL0nbNz6tD7KscSUiIqLyYjkBWV9mEpCXAeRlS4F1+CbAr6kUYL9/BlBnSfvlpgMefgBY40pERETlwxBL1pWvAv7aDAS0AZL/BNQZUpDtuwDYM1MKuE4K4Km1UrAthjWuREREZCmWE5B1FaiBvEygUC0FWbmHFGR/fl0KsDJHoFF3QFnf1i0lIiKiaowhlqzL3RfoFQV4+EtB1jsM0BZKPbRCCzR8GBi0XNqPiIiI6D4xxJL16YKskwJIPFYUYGUOwL0EIOeurVtIRERE1RxDLFWOnLvArThAWwA4OAP1Wkq9s7rBXuam3yIiIiKyAEMsWZ9uFoLsFCnA1m8nBdig9kU1sgyyREREVAEMsWSSRisQeyUNP8ffROyVNGi0ptbTMns0kJshDeJq2keqgdXVyOoGe7l4AK51Kqv5REREVMNxii0ysvtMEubtOIekDLV+W6BSgbmDI8peeECVBsT+BwhsAyiURYO4ekUBB6Kk3tmANkCPWfo5YomIiIjKiz2xZGD3mSSM3xBnEGABIDlDjfEb4rD7TFLpT+CsAOReQJ1gw1kIis9aoGwAKIMr5w0QERFRrSATQpTne+JqJzMzE0qlEhkZGfDy8rJ1c+yaRivQZclBowCrI4O0FOzvM3rqV9LSaIXxKluFOdJ8saam0VKlSUHXxb0S3wkRERFVV5ZmN5YTkN7xhLtmAywACABJGWocT7iLyMa+91d2wPlhiYiIyApYTkB6KVnmAywAuEINb2QiJUttVHbgjUy4Qm152QERERFRBbAnlvT8PRVmH3OFGlOctsATOQhwehCTd9yArg7FG5mY7rQZWXDDisKnoIYC83acQ5+IAH3ZAREREZE1sSeW9DqG+SBQqYCp2KlAPryQg2CXLIT8sQzqjBQARQHWT5YOT+RAgXyDsgMiIiKiysAQS3qODjLMHRwBAEZBNh1eWFo4DCEhYUB2CqY7bUZj2U19gE0VdbC0cBjuoagAu6zyBCIiIqL7xRBLBvq3DMSqF9ojQGlYWhCgVGDhC90R/PQSyDzrwU+WjllO35kNsEDp5QlEREREFcGaWDLSv2Ug+kQEGE+d9W99a90+U6G4OgF5BRoIAF9pBhoEWN1UXB3DfGzzBoiIiKjGY08smeToIENkY18MaVsfkY19iwZoqdLgeOw/aFbPA4AUWF923AVvZOrvA8DcwREc1EVERESVhiGWLKdK0y8d6x/YEC4DPkC2sy/8ZOmY7rQZ3shEgFKBVS+0L3t5WiIiIqIKYDkBWaZYgIWHP9ArCl3dfaFp3Rx3tr0DkXUbPT2jUXfoQjh61LV1a4mIiKiGY08sWcZZAci99AFWt/KWo0dd1HtiEQLqN0Q9Pz84urjatp1ERERUK8iEEKLs3aovS9ffJQvkq4ACtemlY1VpUtB1ca/6dhEREVGNYWl2YzkBWc7F3XxINRVsiYiIiCqJTcsJfv31VwwePBhBQUGQyWTYtm2bweNCCERFRSEoKAiurq7o3r07zp49a5vGEhEREZHdsGmIValUaNOmDVauXGny8aVLl2L58uVYuXIlTpw4gYCAAPTp0wdZWVlV3FIiIiIisic2LScYMGAABgwYYPIxIQQ+/vhjvPvuu3jyyScBAN988w3q1auH7777DuPGjavKphIRERGRHbHb2QkSEhKQnJyMvn376rfJ5XJ069YNR44csWHLiIiIiMjW7HZgV3JyMgCgXr16Btvr1auH69evmz0uLy8PeXl5+vuZmZmV00AiIiIishm77YnVkckMly4VQhhtK27RokVQKpX6W3BwcGU3sXrKV0nTYpmiSpMeJyIiIrJTdhtiAwICABT1yOqkpKQY9c4WN2vWLGRkZOhvN27cqNR2Vkv5KiB6obQCV8kgq1uZK3ohgywRERHZLbsNsWFhYQgICMC+ffv02/Lz8xETE4POnTubPU4ul8PLy8vgRiUUqIG8TGkJ2eJBtvjSsnmZ0n5EREREdsimITY7Oxvx8fGIj48HIA3mio+PR2JiImQyGSZPnoyFCxdi69atOHPmDMaMGQM3NzeMGDHCls2u/tx9paVjPfyLgmzqhaIAW2JpWSIiIiJ7Y9NlZw8dOoQePXoYbR89ejTWrVsHIQTmzZuHL774Avfu3cPDDz+M//znP2jZsqXFr8FlZ0tRvOdVhwGWiIiIbMjS7GbTEFsVGGLLkHoB2Den6H6f+YBfuO3aQ0RERLWapdnNbmtiqQqo0oDYEqulxa40P2sBERERkZ1giK2tipcSePhLPbDFa2QZZImIiMiOMcTWBiXnhC0eYBVK4NFpUglBycFeDLJERERkpxhiazpTc8I6KwC5lxRgIYBjX0j7FZ+1QO4l7UdERERkhxhia7J8FZBxy3hOWBd3oNUwQFMAqDMN54TVBdke70j7EREREdkhhtiaStcDG/spEDnJsEwg8Riw5UXgVhzg6m08pZa7LwMsERER2TWG2Jqq+KpcxYNsxj/Aj2MBdQbg4AR0eYtzwhIREVG1wxBbU5VclSv2UyB8EJB6EdDkAY5yYMjngF9TW7eUiIiIqNwYYmuy4kE24x9g51tFAdavGRC/gTMQEBERUbXEEFvTZKcCqZeK7rv7Am1fkHpgC1SAzBEY9BGgbMCptIiIiKjaYoitSbJTgY1PA5uGFwXZ1EvAz69LAVadJQ34+utH48FeDLJERERUjTDE1iS56UB+tjRoa9NwIPEEcHg5kJ8D5GUDcg9AJgOykwwHe3FOWCIiIqpmZEIIYetGVKbMzEwolUpkZGTAy8vL1s2pfKmXpACbe0+aB9a3CZB2CXBwBtx8gCGrgfj1RcvNRk4ClEGcUouIiIjsgqXZjT2xNY1fU2D4JkBRB9AWACnnAEcXKcAO3wQ0fMhwVS4GWCIiIqqGnGzdAKoEfk2BfguBnycAEIDMAei7oGg6Ld2sBc4KBlgiIiKqltgTWxOlXgL2vivVv8r+/RHvfdd41gIGWCIiIqqmGGJrGl1NrDoDUCiBx1dKf+oGexUPskRERETVFENsTVIywA7fBIT3/7dGlkGWiIiIag6G2JrEtQ7g4lEUYHU1sPrBXkrpcdc6tmwlERERUYVxYFdN4uEHPP+jNF+sLsDq6IKsax1pPyIiIqJqjCG2OspXAQVqaXBWSTIHadosU0oGWyIiIqJqiuUE1U2+CoheaHqpWFWatD16obQfERERUQ3FEFvdFKihVWfi7u1/cO2H6Thx9hI0WlEUYLNTgLxMqaeWiIiIqIZiOUE1szshHx+f7YnROd/AT3YFqVem4SPXJ7CswW8Ids6SVuLqFWW61ICIiIiohmBPbDWy+0wSxm+Iw9+ZLlhaOAypog78ZOl4Rb0WFy9fxo0CTwZYIiIiqhUYYqsJjVZg3o5zEP/evwcvfKUZCAD6bdP+6QqNq49N2kdERERUlRhiq4N8FeLOX0ZSRlGdqzcy8bLjLjijEA7QQAAYkrsVcecvA5BCb+yVNPwcfxOxV9KkulkiIiKiGoI1sfbu39kI/G4lwxu9cQ9e8EYmpjttRpDsDurL7iBOK02d5SdLh9/xxTigfguz9yYZhN5ApQJzB0egf8tAW70TIiIiIqthT6y9K1ADeZnwKryL6U6bEYZbBgH2pqgLDRyxqvBxpIo60GbeRurWd6DOSDF4muQMNcZviMPuM0k2eiNERERE1sMQa+/cfYFeUajj3wDBLll43Wk7HKHRB9hboi6WFg7DNQRhneso/JUhRybcoIaLwdPoignm7TjH0gIiIiKq9lhOYMc0WoHjCXeRkqVGUJNJaHhvMfKuXIEDtAYBNh1eAICBD7dE1P6noYYLcqEwej4BIClDjeMJdxHZmDMYEBERUfXFEGundp9Jwrwd5wzqWjt5Por3fZNwOzMP6gINvtIMxD146etd8wq1uPdvoC1NShYXQiAiIqLqjSHWDunmgy3+pb83MjE0dysS1TloWV8JZ0cH/J/TMaR2nIn2LZrA0UGG2CtpZp+zOH9P415aIiIiouqENbF2puR8sAD0sxHUlaUjVdTBm2lPoo5/A4TKs/HQlU/hmHsXANAxzAeBSgVkZp5bBmmWgo5hnEuWiIiIqje7DrGFhYWYPXs2wsLC4OrqikaNGmH+/PnQarW2blqlOZ5w12g+2OlOm+H3b4BdUjgMR7Pq4lSTSdISs9kpwIEoQJUGRwcZ5g6OAACjIKu7P3dwBBwdzMVcIiIiourBrkPskiVLsHr1aqxcuRLnz5/H0qVLsWzZMnz22We2blqlKVmvqoYLsuCGVFEHSwuH6WtebxW4SUvMevgDci/AWSoR6N8yEKteaI8ApWHJQIBSgVUvtOc8sURERFQj2HVNbGxsLIYMGYJBgwYBAEJDQ/H999/j5MmTNm5Z5QlQaOGNTH1YzYUCKwqfggL5AABXqJELhVTX+u/0W3BWAC7u+ufo3zIQfSIC9DMb+HtKJQTsgSUiIqKawq57Yrt06YIDBw7g4sWLAIA///wTv//+OwYOHGjjllWSfBUeurYaUW4/wgeZ+s266bKmO23GVKctCPNCUV2ru69BgNVxdJAhsrEvhrStj8jGvgywREREVKPYdU/sjBkzkJGRgebNm8PR0REajQYLFizAc889Z/aYvLw85OXl6e9nZmaa3deu5KuAjFtwyM9EZL1CTL+xGUsLh+EuvBCGW3jdaTt8ZFmQCWB2vzCGUiIiIqrV7DrEbt68GRs2bMB3332HBx54APHx8Zg8eTKCgoIwevRok8csWrQI8+bNq+KWVlC+CoheCORlApGT4B/7KXoiEa63f8Q3uY9gqfMaOEODv5wegN+ghej1YIStW0xERERkUzIhhN2uQRocHIyZM2diwoQJ+m0ffPABNmzYgL///tvkMaZ6YoODg5GRkQEvr7IXArAJVZo0w0B2ijRQK3ISEPspRPo/yE/+G6IwD1qFEvIXNsPRv5mtW0tERERUaTIzM6FUKsvMbnZdE5uTkwMHB8MmOjo6ljrFllwuh5eXl8HN7ukGaOmmzIr9FAgfBNmdi5DLCqBQuMLt6dUMsERERET/susQO3jwYCxYsAA7d+7EtWvXsHXrVixfvhxPPPGErZtmfcWDbMY/wM63AE0e4CgH/JoB8RukHlsiIiIisu9ygqysLLz33nvYunUrUlJSEBQUhOeeew5z5syBi4uLRc9haZe03Ug8Bvw4tijADvoIuLCzqNSgV5QUeImIiIhqIEuzm12HWGuoViE29RKwaTigzijqgVU20NfIMsgSERFRTVcjamJrtHyVYXmAKg34/SNAWwjIPYGh/5ECrK5GNnKS0epcRERERLWVXU+xVWMVn1KrV5S07UAUkHsPqP+gFGQv7QMenQb8uswwyCqDTC5uQERERFSbsCfWFgrUUoDNTpHCa36O1MOqUAIOjkBhnvS4s3vRYC+5FwMsERER0b9YE2srJeeG7fAScPQ/gDrTuO5VlSaVEDDAEhERUQ3Hmlh7V3Ju2EOLTAdY3b4MsERERER6DLG25O4LRE403BY5kTMPEBEREZWBIdaWVGlA7ErDbbEruagBERERURkYYm2lZE1sn/lFpQUHohhkiYiIiErBEGsLJQNsryjAL9ywRpZBloiIiMgshlhbcFZIU2aVHMRVfLAXFzUgIiIiMotTbNlKvkqaL9bUIC5OqUVERES1lKXZjSt22YqLu/mQytkJiIiIiErFcgIiIiIiqnYYYomIiIio2mGIJSIiIqJqhyGWiIiIiKodhlgiIiIiqnYYYomIiIio2mGItbZ8lfmVtlRp0uNEREREVCEMsdaUrwKiF5peMla31Gz0QgZZIiIiogpiiLWmAjWQlwlkpxgGWV2AzU6RHi9Q27KVRERERNUeQ6w1ufsCvaIAD/+iIJt6oSjAevhLj3NFLiIiIqIKYYi1tpJBdt8cBlgiIiIiK2OIrQzuvkDkRMNtkRMZYImIiIishCHWGkrOSKBKA2JXSn/X5APaQum+uVkLiIiIiKhcGGIrquSMBMUHcSm8ALknkHlLupmatYCIiIiIys3J1g2o9orPSLBnFgABqDOlAAsZUJgHBLSStusGe7E2loiIiKhC2BNbUcUHcuWkAclnAEcXADJAnSFt77cI6LdY+rvcC3BW2LrVRERERNUae2KtQRdkD0RJ9/OzpVrYkjMS9IqSAqyLu23aSURERFRDsCfWWnQzEjg4/dsTC+MZCdx9GWCJiIiIrIAh1lqKz0igwxkJiIiIiCoFQ6w1FJ+RwMMf6DPfcNUuBlkiIiIiq2KIraiSAbZXFOAXbrz8LIMsERERkdUwxFaUs0KacaDkIK7isxZwRgIiIiIiq7L7EHvz5k288MIL8PX1hZubG9q2bYtTp07ZullFXNyBHu+YnvtVF2R7vMMBXURERERWZNdTbN27dw+PPPIIevTogf/973/w9/fHlStXUKdOHVs3zZCLu/mQykUNiIiIiKzOrkPskiVLEBwcjLVr1+q3hYaG2q5BRERERGQX7LqcYPv27ejQoQOeeeYZ+Pv7o127dvjyyy9LPSYvLw+ZmZkGNyIiIiKqWew6xF69ehWrVq1C06ZNsWfPHrz22muYNGkS1q9fb/aYRYsWQalU6m/BwcFV2GIiIiIiqgoyIYSwdSPMcXFxQYcOHXDkyBH9tkmTJuHEiROIjY01eUxeXh7y8vL09zMzMxEcHIyMjAx4eXlVepuJiIiI6P5lZmZCqVSWmd3suic2MDAQERERBttatGiBxMREs8fI5XJ4eXkZ3IiIiIioZrHrEPvII4/gwoULBtsuXryIkJAQG7WIiIiIiOyBXYfYKVOm4OjRo1i4cCEuX76M7777DmvWrMGECRNs3TQiIiIisiG7DrEPPfQQtm7diu+//x4tW7bE+++/j48//hjPP/+8rZtGRERERDZk1wO7rCEjIwN16tTBjRs3WB9LREREZOd0g/LT09OhVCrN7mfXix1YQ1ZWFgBwqi0iIiKiaiQrK6vUEFvje2K1Wi1u3boFT09PyGQym7RB9z8K9gYb4nkxjefFNJ4X83huTON5MY3nxTSeF/Oq+twIIZCVlYWgoCA4OJivfK3xPbEODg5o0KCBrZsBAJzyywyeF9N4XkzjeTGP58Y0nhfTeF5M43kxryrPTWk9sDp2PbCLiIiIiMgUhlgiIiIiqnYYYquAXC7H3LlzIZfLbd0Uu8LzYhrPi2k8L+bx3JjG82Iaz4tpPC/m2eu5qfEDu4iIiIio5mFPLBERERFVOwyxRERERFTtMMQSERERUbXDEFtJoqKiIJPJDG4BAQG2bpZN/Prrrxg8eDCCgoIgk8mwbds2g8eFEIiKikJQUBBcXV3RvXt3nD171jaNrUJlnZcxY8YYXUOdOnWyTWOryKJFi/DQQw/B09MT/v7+GDp0KC5cuGCwT229Xiw5N7Xxmlm1ahVat26tn78yMjIS//vf//SP19brpazzUhuvFVMWLVoEmUyGyZMn67fV1mumOFPnxR6vGYbYSvTAAw8gKSlJf/vrr79s3SSbUKlUaNOmDVauXGny8aVLl2L58uVYuXIlTpw4gYCAAPTp00e/ZHBNVdZ5AYD+/fsbXEO7du2qwhZWvZiYGEyYMAFHjx7Fvn37UFhYiL59+0KlUun3qa3XiyXnBqh910yDBg2wePFinDx5EidPnkTPnj0xZMgQfeiorddLWecFqH3XSkknTpzAmjVr0Lp1a4PttfWa0TF3XgA7vGYEVYq5c+eKNm3a2LoZdgeA2Lp1q/6+VqsVAQEBYvHixfptarVaKJVKsXr1ahu00DZKnhchhBg9erQYMmSITdpjL1JSUgQAERMTI4Tg9VJcyXMjBK8ZHW9vb/HVV1/xeilBd16E4LWSlZUlmjZtKvbt2ye6desm3nzzTSEEP2PMnRch7POaYU9sJbp06RKCgoIQFhaG4cOH4+rVq7Zukt1JSEhAcnIy+vbtq98ml8vRrVs3HDlyxIYtsw+HDh2Cv78/mjVrhldeeQUpKSm2blKVysjIAAD4+PgA4PVSXMlzo1ObrxmNRoNNmzZBpVIhMjKS18u/Sp4Xndp8rUyYMAGDBg1C7969DbbX9mvG3HnRsbdrxsmmr16DPfzww1i/fj2aNWuG27dv44MPPkDnzp1x9uxZ+Pr62rp5diM5ORkAUK9ePYPt9erVw/Xr123RJLsxYMAAPPPMMwgJCUFCQgLee+899OzZE6dOnbK7CacrgxACU6dORZcuXdCyZUsAvF50TJ0boPZeM3/99RciIyOhVqvh4eGBrVu3IiIiQh86auv1Yu68ALX3WgGATZs2IS4uDidOnDB6rDZ/xpR2XgD7vGYYYivJgAED9H9v1aoVIiMj0bhxY3zzzTeYOnWqDVtmn2QymcF9IYTRttpm2LBh+r+3bNkSHTp0QEhICHbu3Iknn3zShi2rGhMnTsTp06fx+++/Gz1W268Xc+emtl4z4eHhiI+PR3p6OrZs2YLRo0cjJiZG/3htvV7MnZeIiIhae63cuHEDb775Jvbu3QuFQmF2v9p2zVhyXuzxmmE5QRVxd3dHq1atcOnSJVs3xa7oZmzQ/e9XJyUlxeh/wrVdYGAgQkJCasU19MYbb2D79u2Ijo5GgwYN9Nt5vZg/N6bUlmvGxcUFTZo0QYcOHbBo0SK0adMGn3zySa2/XsydF1Nqy7Vy6tQppKSk4MEHH4STkxOcnJwQExODTz/9FE5OTvrrorZdM2WdF41GY3SMPVwzDLFVJC8vD+fPn0dgYKCtm2JXwsLCEBAQgH379um35efnIyYmBp07d7Zhy+xPWloabty4UaOvISEEJk6ciJ9++gkHDx5EWFiYweO1+Xop69yYUhuuGVOEEMjLy6vV14spuvNiSm25Vnr16oW//voL8fHx+luHDh3w/PPPIz4+Ho0aNaqV10xZ58XR0dHoGLu4Zmw1oqyme+utt8ShQ4fE1atXxdGjR8Vjjz0mPD09xbVr12zdtCqXlZUl/vjjD/HHH38IAGL58uXijz/+ENevXxdCCLF48WKhVCrFTz/9JP766y/x3HPPicDAQJGZmWnjlleu0s5LVlaWeOutt8SRI0dEQkKCiI6OFpGRkaJ+/fo1+ryMHz9eKJVKcejQIZGUlKS/5eTk6PeprddLWeemtl4zs2bNEr/++qtISEgQp0+fFu+8845wcHAQe/fuFULU3uultPNSW68Vc0qOwq+t10xJxc+LvV4zDLGVZNiwYSIwMFA4OzuLoKAg8eSTT4qzZ8/aulk2ER0dLQAY3UaPHi2EkKY0mTt3rggICBByuVw8+uij4q+//rJto6tAaeclJydH9O3bV/j5+QlnZ2fRsGFDMXr0aJGYmGjrZlcqU+cDgFi7dq1+n9p6vZR1bmrrNTN27FgREhIiXFxchJ+fn+jVq5c+wApRe6+X0s5Lbb1WzCkZYmvrNVNS8fNir9eMTAghqq7fl4iIiIio4lgTS0RERETVDkMsEREREVU7DLFEREREVO0wxBIRERFRtcMQS0RERETVDkMsEREREVU7DLFEREREVO0wxBIRERFRtcMQS0QW6d69OyZPnqy/Hxoaio8//thm7amOSp7D6ubQoUOQyWRIT0+v0PPIZDJs27atyl9XJz8/H02aNMHhw4et8nym5OXloWHDhjh16lSlvQZRbccQS0T35cSJE3j11Vdt3QysW7cOderUsWkbxowZg6FDh9q0DRVx7do1yGQyxMfHV8nrJSUlYcCAAVZ9zqioKLRt29aifdesWYOQkBA88sgjVm1DcXK5HG+//TZmzJhRaa9BVNsxxBLRffHz84Obm5utm2E1Go0GWq3W1s0wkp+fb+smWF1AQADkcrnNXv+zzz7Dyy+/XOmv8/zzz+O3337D+fPnK/21iGojhlgiMqJSqTBq1Ch4eHggMDAQH330kdE+JcsJli9fjlatWsHd3R3BwcF4/fXXkZ2drX9c12P6yy+/IDw8HG5ubnj66aehUqnwzTffIDQ0FN7e3njjjTeg0Wj0x+Xn52P69OmoX78+3N3d8fDDD+PQoUMApK+ZX3zxRWRkZEAmk0EmkyEqKqrM40q2JyIiAnK5HNevXzd6nxqNBi+99BLCwsLg6uqK8PBwfPLJJ/rHo6Ki8M033+Dnn3/Wt6H465RUWFiIiRMnok6dOvD19cXs2bMhhDA4rx988AHGjBkDpVKJV155BQBw5MgRPProo3B1dUVwcDAmTZoElUqlP27Dhg3o0KEDPD09ERAQgBEjRiAlJUX/+L179/D888/Dz88Prq6uaNq0KdauXQsACAsLAwC0a9cOMpkM3bt3N9t+ADh16hQ6dOgANzc3dO7cGRcuXDB4fMeOHXjwwQehUCjQqFEjzJs3D4WFhfrHS5YTHDlyBG3btoVCoUCHDh2wbds2kz3D5l533bp1mDdvHv7880/9z2DdunUm2x4XF4fLly9j0KBBBtv/+ecfDB8+HD4+PnB3d0eHDh1w7NgxAEW9vF9//TUaNmwIDw8PjB8/HhqNBkuXLkVAQAD8/f2xYMECg+f09fVF586d8f3335d6PonoPgkiohLGjx8vGjRoIPbu3StOnz4tHnvsMeHh4SHefPNN/T4hISFixYoV+vsrVqwQBw8eFFevXhUHDhwQ4eHhYvz48frH165dK5ydnUWfPn1EXFyciImJEb6+vqJv377i2WefFWfPnhU7duwQLi4uYtOmTfrjRowYITp37ix+/fVXcfnyZbFs2TIhl8vFxYsXRV5envj444+Fl5eXSEpKEklJSSIrK6vM44q3p3PnzuLw4cPi77//FtnZ2UbnIj8/X8yZM0ccP35cXL16VWzYsEG4ubmJzZs3CyGEyMrKEs8++6zo37+/vg15eXkmz2u3bt305/Hvv//WP9eaNWsMzquXl5dYtmyZuHTpkrh06ZI4ffq08PDwECtWrBAXL14Uhw8fFu3atRNjxozRH/d///d/YteuXeLKlSsiNjZWdOrUSQwYMED/+IQJE0Tbtm3FiRMnREJCgti3b5/Yvn27EEKI48ePCwBi//79IikpSaSlpZlsf3R0tAAgHn74YXHo0CFx9uxZ0bVrV9G5c2f9Prt37xZeXl5i3bp14sqVK2Lv3r0iNDRUREVF6fcBILZu3SqEECIzM1P4+PiIF154QZw9e1bs2rVLNGvWTAAQf/zxh0Wvm5OTI9566y3xwAMP6H8GOTk5Jt/DihUrRPPmzQ22ZWVliUaNGomuXbuK3377TVy6dEls3rxZHDlyRAghxNy5c4WHh4d4+umnxdmzZ8X27duFi4uL6Nevn3jjjTfE33//Lb7++msBQMTGxho89/Tp00X37t1NtoWIKoYhlogMZGVlGQXJtLQ04erqWmqILemHH34Qvr6++vtr164VAMTly5f128aNGyfc3Nz0wVMIIfr16yfGjRsnhBDi8uXLQiaTiZs3bxo8d69evcSsWbP0z6tUKg0et/Q4ACI+Pr6Us2Ha66+/Lp566in9/dGjR4shQ4aUeVy3bt1EixYthFar1W+bMWOGaNGihf5+SEiIGDp0qMFxI0eOFK+++qrBtt9++004ODiI3Nxck6+lC6a6czt48GDx4osvmtw3ISHBIDSaowuT+/fv12/buXOnAKBvR9euXcXChQsNjvv2229FYGCg/n7xELtq1Srh6+tr8D6+/PJLkyG2tNedO3euaNOmTantF0KIN998U/Ts2dNg2xdffCE8PT3Nhve5c+cKNzc3kZmZqd/Wr18/ERoaKjQajX5beHi4WLRokcGxn3zyiQgNDS2zXURUfk5V3/dLRPbsypUryM/PR2RkpH6bj48PwsPDSz0uOjoaCxcuxLlz55CZmYnCwkKo1WqoVCq4u7sDANzc3NC4cWP9MfXq1UNoaCg8PDwMtum+Bo+Li4MQAs2aNTN4rby8PPj6+ppti6XHubi4oHXr1qW+LwBYvXo1vvrqK1y/fh25ubnIz8+3eBBRSZ06dYJMJtPfj4yMxEcffQSNRgNHR0cAQIcOHQyOOXXqFC5fvoyNGzfqtwkhoNVqkZCQgBYtWuCPP/5AVFQU4uPjcffuXX19b2JiIiIiIjB+/Hg89dRTiIuLQ9++fTF06FB07tz5vt5D8XMWGBgIAEhJSdGPxj9x4oTBV+sajQZqtRo5OTlGddQXLlxA69atoVAo9Ns6duxY7te1VG5ursFrAUB8fDzatWsHHx8fs8eFhobC09NTf79evXpwdHSEg4ODwbbiJRwA4OrqipycHIvbR0SWY4glIgOiWH2mpa5fv46BAwfitddew/vvvw8fHx/8/vvveOmll1BQUKDfz9nZ2eA4mUxmcpsugGm1Wjg6OuLUqVP6gKdTPPiWZOlxrq6uBoHSlB9++AFTpkzBRx99hMjISHh6emLZsmX6esnKoAv9OlqtFuPGjcOkSZOM9m3YsCFUKhX69u2Lvn37YsOGDfDz80NiYiL69eunHxg2YMAAXL9+HTt37sT+/fvRq1cvTJgwAR9++GG521f8Z6Y7f8V/ZvPmzcOTTz5pdFzJ8AhI11vJn4G5a7C017VU3bp18ddffxlsc3V1LfO48l67Onfv3oWfn1+52khElmGIJSIDTZo0gbOzM44eParv4bp37x4uXryIbt26mTzm5MmTKCwsxEcffaTvmfrhhx8q3JZ27dpBo9EgJSUFXbt2NbmPi4uLwUAwS4+z1G+//YbOnTvj9ddf12+7cuVKmW0w5+jRo0b3mzZtahS2i2vfvj3Onj2LJk2amHz8r7/+wp07d7B48WIEBwcDkH4mJfn5+WHMmDEYM2YMunbtimnTpuHDDz+Ei4sLAFj8HkrTvn17XLhwwWxbS2revDk2btyIvLw8/YwFptpeFkt/Bu3atcOqVasMwnPr1q3x1Vdf4e7du6X2xt6PM2fOoF27dlZ9TiKScHYCIjLg4eGBl156CdOmTcOBAwdw5swZjBkzxuBr05IaN26MwsJCfPbZZ7h69Sq+/fZbrF69usJtadasGZ5//nmMGjUKP/30ExISEnDixAksWbIEu3btAiB9zZudnY0DBw7gzp07yMnJseg4SzVp0gQnT57Enj17cPHiRbz33ns4ceKEwT6hoaE4ffo0Lly4gDt37hj0Ppd048YNTJ06FRcuXMD333+Pzz77DG+++WapbZgxYwZiY2MxYcIExMfH49KlS9i+fTveeOMNAFJvrIuLi/78b9++He+//77Bc8yZMwc///wzLl++jLNnz+KXX35BixYtAAD+/v5wdXXF7t27cfv2bWRkZJTrHJV8nfXr1yMqKgpnz57F+fPnsXnzZsyePdvk/iNGjIBWq8Wrr76K8+fPY8+ePfre4bJ6yYsLDQ1FQkIC4uPjcefOHeTl5Zncr0ePHlCpVDh79qx+23PPPYeAgAAMHToUhw8fxtWrV7FlyxbExsaW452b9ttvv6Fv374Vfh4iMsYQS0RGli1bhkcffRSPP/44evfujS5duuDBBx80u3/btm2xfPlyLFmyBC1btsTGjRuxaNEiq7Rl7dq1GDVqFN566y2Eh4fj8ccfx7Fjx/Q9jp07d8Zrr72GYcOGwc/PD0uXLrXoOEu99tprePLJJzFs2DA8/PDDSEtLM+iVBYBXXnkF4eHh6NChA/z8/EpdCWrUqFHIzc1Fx44dMWHCBLzxxhtlLhrRunVrxMTE4NKlS+jatSvatWuH9957T18X6ufnh3Xr1uG///0vIiIisHjxYqMyARcXF8yaNQutW7fGo48+CkdHR2zatAkA4OTkhE8//RRffPEFgoKCMGTIkHKdo+L69euHX375Bfv27cNDDz2ETp06Yfny5QgJCTG5v5eXF3bs2IH4+Hi0bdsW7777LubMmQPAdPmBOU899RT69++PHj16wM/Pz+y0Vr6+vnjyyScN6otdXFywd+9e+Pv7Y+DAgWjVqhUWL15cau+4JWJjY5GRkYGnn366Qs9DRKbJxP0UwBEREVWSjRs36uf/taRetbz++usv9O7dG5cvXzYYrGVtzzzzDNq1a4d33nmn0l6DqDZjTSwREdnU+vXr0ahRI9SvXx9//vknZsyYgWeffbZSAiwAtGrVCkuXLsW1a9fQqlWrSnmNvLw8tGnTBlOmTKmU5yci9sQSEZGNLV26FJ9//jmSk5MRGBiIoUOHYsGCBTVqWWMisj6GWCIiIiKqdjiwi4iIiIiqHYZYIiIiIqp2GGKJiIiIqNphiCUiIiKiaochloiIiIiqHYZYIiIiIqp2GGKJiIiIqNphiCUiIiKiaochloiIiIiqnf8H3GMoAoMjiVwAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 700x450 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(7, 4.5))\n",
    "for label, marker, alpha in ((\"measured\", \"o\", 1.0), (\"imputed\", \"x\", 0.7)):\n",
    "    subset = summary[summary[\"provenance\"] == label]\n",
    "    ax.scatter(subset[\"dbh_cm\"], subset[\"height_m\"], marker=marker, alpha=alpha, label=label)\n",
    "ax.set_xlabel(\"diameter at breast height (cm)\")\n",
    "ax.set_ylabel(\"height (m)\")\n",
    "ax.set_title(\"Measured heights and the Näslund curve fitted through them\")\n",
    "ax.legend()\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c10",
   "metadata": {},
   "source": [
    "The imputed points lie exactly on the fitted curve, which is what an\n",
    "interpolated value *is* -- it carries no residual scatter. The measured points\n",
    "scatter around it. Anything that treats the two alike will understate the\n",
    "variance it is working with."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c11",
   "metadata": {},
   "source": [
    "## 3. Competition indices\n",
    "\n",
    "`competition_indices` takes a plot (or a stand, or a plain list of trees), the\n",
    "indices you want, and a rule for choosing competitors. Here are three indices\n",
    "with contrasting behaviour:\n",
    "\n",
    "- **`Heg`** (Hegyi 1974) -- distance-dependent size ratio; higher means more\n",
    "  competition.\n",
    "- **`BAL`** (Wykoff et al. 1982) -- basal area of larger trees, distance\n",
    "  independent.\n",
    "- **`drg`** (Hamilton 1986) -- the subject's diameter over the plot's quadratic\n",
    "  mean diameter. Note this one runs the *other* way: a high value means the\n",
    "  subject dominates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c12",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-26T12:01:50.240637Z",
     "iopub.status.busy": "2026-07-26T12:01:50.240637Z",
     "iopub.status.idle": "2026-07-26T12:01:50.262474Z",
     "shell.execute_reply": "2026-07-26T12:01:50.262474Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  uid  dbh_cm  n_comp    Heg    BAL   drg\n",
      "p1t37     4.7      19 24.554 39.106 0.252\n",
      "p1t31     4.9      29 31.089 39.046 0.263\n",
      "p1t02     5.6      13 23.026 38.968 0.300\n",
      "p1t20     5.9      19 26.137 38.881 0.316\n",
      "p1t28     6.1      24 20.130 38.788 0.327\n",
      "   ...\n",
      "  uid  dbh_cm  n_comp   Heg   BAL   drg\n",
      "p1t29    24.4      24 5.038 9.669 1.308\n",
      "p1t25    26.3      12 3.517 7.939 1.410\n",
      "p1t11    26.4      21 3.837 6.197 1.415\n",
      "p1t44    34.6      10 3.183 3.204 1.854\n",
      "p1t34    35.8      11 1.814 0.000 1.919\n"
     ]
    }
   ],
   "source": [
    "plot = stand.plots[0]\n",
    "results = competition_indices(\n",
    "    plot,\n",
    "    indices=[\"Heg\", \"BAL\", \"drg\"],\n",
    "    selector=FixedRadius(8.0, min_size_ratio=0.3),\n",
    ")\n",
    "\n",
    "table = pd.DataFrame(\n",
    "    {\n",
    "        \"uid\": [r.tree.uid for r in results],\n",
    "        \"dbh_cm\": [float(r.tree.diameter_cm) for r in results],\n",
    "        \"n_comp\": [r.n_competitors for r in results],\n",
    "        **{name: [r.indices[name] for r in results] for name in (\"Heg\", \"BAL\", \"drg\")},\n",
    "    }\n",
    ").sort_values(\"dbh_cm\")\n",
    "\n",
    "print(table.head(5).round(3).to_string(index=False))\n",
    "print(\"   ...\")\n",
    "print(table.tail(5).round(3).to_string(index=False))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c13",
   "metadata": {},
   "source": [
    "Small trees carry high `Heg` and `BAL` and low `drg`; large trees the reverse.\n",
    "That is the sanity check you want -- a suppressed stem is under more competition\n",
    "than a dominant one.\n",
    "\n",
    "`n_comp` counts the competitors the *selector* kept, and it drives `Heg` only.\n",
    "`BAL` and `drg` are plot descriptors: they are computed over every other tree on\n",
    "the plot whatever the selector does, which is what makes `BAL` comparable\n",
    "between runs that used different selection rules.\n",
    "\n",
    "In the correlation below, `drg` sits at exactly 1.000 against diameter. That is\n",
    "not a finding: `drg` *is* the subject's diameter divided by a plot-level\n",
    "constant, so within one plot it is diameter on a different scale. It earns its\n",
    "place as a predictor across plots of differing size, not within one.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c14",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-26T12:01:50.262474Z",
     "iopub.status.busy": "2026-07-26T12:01:50.262474Z",
     "iopub.status.idle": "2026-07-26T12:01:50.273416Z",
     "shell.execute_reply": "2026-07-26T12:01:50.273416Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "        dbh_cm    Heg    BAL    drg\n",
      "dbh_cm   1.000 -0.870 -0.924  1.000\n",
      "Heg     -0.870  1.000  0.692 -0.870\n",
      "BAL     -0.924  0.692  1.000 -0.924\n",
      "drg      1.000 -0.870 -0.924  1.000\n"
     ]
    }
   ],
   "source": [
    "print(table[[\"dbh_cm\", \"Heg\", \"BAL\", \"drg\"]].corr().round(3).to_string())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c15",
   "metadata": {},
   "source": [
    "## 4. The selection rule is part of the index\n",
    "\n",
    "An index value is only interpretable next to the rule that chose its\n",
    "competitors. The same `Heg` formula on the same plot gives materially different\n",
    "numbers under different rules, so the rule has to be reported alongside the\n",
    "index."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c16",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-26T12:01:50.273416Z",
     "iopub.status.busy": "2026-07-26T12:01:50.273416Z",
     "iopub.status.idle": "2026-07-26T12:01:50.318436Z",
     "shell.execute_reply": "2026-07-26T12:01:50.318436Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "              selector  mean_competitors  mean_Heg         citation\n",
      "      FixedRadius(8 m)            18.533    10.012 (plain geometry)\n",
      " MeanHeightRadius(0.4)            10.578     6.891        Sims 2009\n",
      "MeanHeightRadius(0.25)             4.444     4.191        Sims 2009\n",
      "   LeeGadowRadius(k=2)             9.067     6.283         Lee 1997\n",
      "   LeeGadowRadius(k=3)            18.133     9.868         Lee 1997\n",
      "      BitterlichBAF(2)            13.000     6.066  Bitterlich 1952\n",
      "  NearestNeighbours(4)             4.000     3.209 (plain geometry)\n",
      "    SearchCone(80 deg)            31.289     8.765    Pretzsch 2009\n"
     ]
    }
   ],
   "source": [
    "selectors = {\n",
    "    \"FixedRadius(8 m)\": FixedRadius(8.0),\n",
    "    \"MeanHeightRadius(0.4)\": MeanHeightRadius(0.4),\n",
    "    \"MeanHeightRadius(0.25)\": MeanHeightRadius(0.25),\n",
    "    \"LeeGadowRadius(k=2)\": LeeGadowRadius(2.0),\n",
    "    \"LeeGadowRadius(k=3)\": LeeGadowRadius(3.0),\n",
    "    \"BitterlichBAF(2)\": BitterlichBAF(2.0),\n",
    "    \"NearestNeighbours(4)\": NearestNeighbours(4),\n",
    "    \"SearchCone(80 deg)\": SearchCone(80.0),\n",
    "}\n",
    "\n",
    "rows = []\n",
    "for label, selector in selectors.items():\n",
    "    res = competition_indices(plot, indices=[\"Heg\"], selector=selector)\n",
    "    src = selector_source(selector)\n",
    "    rows.append(\n",
    "        {\n",
    "            \"selector\": label,\n",
    "            \"mean_competitors\": sum(r.n_competitors for r in res) / len(res),\n",
    "            \"mean_Heg\": sum(r.indices[\"Heg\"] for r in res) / len(res),\n",
    "            \"citation\": f\"{src.author.split(',')[0]} {src.year}\" if src else \"(plain geometry)\",\n",
    "        }\n",
    "    )\n",
    "\n",
    "print(pd.DataFrame(rows).round(3).to_string(index=False))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c17",
   "metadata": {},
   "source": [
    "Note how far apart these are. Every number in these rules is a *parameter*\n",
    "with a literature default, not a constant: `MeanHeightRadius(0.4)` uses the\n",
    "fraction Sims et al. (2009) tested, and dropping it to 0.25 shrinks the zone and\n",
    "the index with it.\n",
    "\n",
    "`SearchCone` looks indiscriminate here, and at this scale it is: an 80 degree\n",
    "cone reaches `h / tan(50 deg)` horizontally, so a 20 m neighbour competes out to\n",
    "almost 17 m -- further than this 10 m plot is wide. The cone discriminates on\n",
    "larger neighbourhoods, or at a narrower angle.\n",
    "\n",
    "`SearchCone` and `BitterlichBAF` are the two rules whose reach depends on the\n",
    "*neighbour's* size rather than the subject's -- the cone on its height, the\n",
    "angle gauge on its diameter. That is what a variable-radius rule means: under\n",
    "`BitterlichBAF(2)` a 60 cm neighbour competes out to 21 m while a 10 cm one\n",
    "stops at 3.5 m. Neither defines a single competition zone, so neither reports a\n",
    "`zone_radius_m`, and neither is edge-corrected -- see section 5.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c18",
   "metadata": {},
   "source": [
    "### The two screens\n",
    "\n",
    "The comparison this set is drawn from applies two screens together for its\n",
    "influence-zone approaches: a neighbour must be at least 30% of the subject's\n",
    "diameter, *and* must not stand in the shadow of a nearer competitor within 30\n",
    "degrees. Reproducing those approaches needs both."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "c19",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-26T12:01:50.319764Z",
     "iopub.status.busy": "2026-07-26T12:01:50.319764Z",
     "iopub.status.idle": "2026-07-26T12:01:50.349866Z",
     "shell.execute_reply": "2026-07-26T12:01:50.349866Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "no screens               mean competitors 18.53   mean Heg 10.012\n",
      "size only (0.3)          mean competitors 16.98   mean Heg  9.856\n",
      "shadow only (30 deg)     mean competitors  9.11   mean Heg  6.849\n",
      "both                     mean competitors  8.76   mean Heg  6.888\n"
     ]
    }
   ],
   "source": [
    "screens = {\n",
    "    \"no screens\": FixedRadius(8.0),\n",
    "    \"size only (0.3)\": FixedRadius(8.0, min_size_ratio=0.3),\n",
    "    \"shadow only (30 deg)\": FixedRadius(8.0, elimination_angle_deg=30.0),\n",
    "    \"both\": FixedRadius(8.0, min_size_ratio=0.3, elimination_angle_deg=30.0),\n",
    "}\n",
    "for label, selector in screens.items():\n",
    "    res = competition_indices(plot, indices=[\"Heg\"], selector=selector)\n",
    "    mean_n = sum(r.n_competitors for r in res) / len(res)\n",
    "    mean_h = sum(r.indices[\"Heg\"] for r in res) / len(res)\n",
    "    print(f\"{label:24s} mean competitors {mean_n:5.2f}   mean Heg {mean_h:6.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c20",
   "metadata": {},
   "source": [
    "## 5. Edge correction\n",
    "\n",
    "A tree near the plot boundary has part of its competition zone outside the plot,\n",
    "so its distance-dependent indices come out too low. Each result reports\n",
    "`observed_zone_fraction`, the share of the zone that was actually inside, and by\n",
    "default the spatial indices are divided by it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c21",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-26T12:01:50.354904Z",
     "iopub.status.busy": "2026-07-26T12:01:50.354904Z",
     "iopub.status.idle": "2026-07-26T12:01:50.376664Z",
     "shell.execute_reply": "2026-07-26T12:01:50.376664Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Trees nearest the plot centre (zone fully observed):\n",
      "  uid  dist_from_centre_m  observed_fraction  Heg_raw  Heg_corrected\n",
      "p1t05                0.73              1.000    5.119          5.119\n",
      "p1t28                1.43              1.000   20.130         20.130\n",
      "p1t14                2.41              0.986    5.349          5.427\n",
      "\n",
      "Trees nearest the boundary (zone partly outside):\n",
      "  uid  dist_from_centre_m  observed_fraction  Heg_raw  Heg_corrected\n",
      "p1t02                9.43              0.455   10.486         23.026\n",
      "p1t18                9.45              0.454    8.708         19.171\n",
      "p1t01                9.96              0.417    4.036          9.682\n"
     ]
    }
   ],
   "source": [
    "raw = competition_indices(\n",
    "    plot, indices=[\"Heg\"], selector=FixedRadius(8.0), edge_correction=None\n",
    ")\n",
    "corrected = competition_indices(plot, indices=[\"Heg\"], selector=FixedRadius(8.0))\n",
    "\n",
    "edge = pd.DataFrame(\n",
    "    {\n",
    "        \"uid\": [r.tree.uid for r in raw],\n",
    "        \"dist_from_centre_m\": [\n",
    "            round(((r.tree.position.X - plot.position.X) ** 2\n",
    "                   + (r.tree.position.Y - plot.position.Y) ** 2) ** 0.5, 2)\n",
    "            for r in raw\n",
    "        ],\n",
    "        \"observed_fraction\": [round(r.observed_zone_fraction, 3) for r in raw],\n",
    "        \"Heg_raw\": [round(r.indices[\"Heg\"], 3) for r in raw],\n",
    "        \"Heg_corrected\": [round(r2.indices[\"Heg\"], 3) for r2 in corrected],\n",
    "    }\n",
    ").sort_values(\"dist_from_centre_m\")\n",
    "\n",
    "print(\"Trees nearest the plot centre (zone fully observed):\")\n",
    "print(edge.head(3).to_string(index=False))\n",
    "print(\"\\nTrees nearest the boundary (zone partly outside):\")\n",
    "print(edge.tail(3).to_string(index=False))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c22",
   "metadata": {},
   "source": [
    "Only the boundary trees change; a tree whose whole zone is inside the plot has\n",
    "`observed_fraction == 1.0` and is left alone.\n",
    "\n",
    "Three things are deliberately *not* corrected:\n",
    "\n",
    "- **Non-spatial indices.** `BAL` and its relatives are plot-level sums and do\n",
    "  not depend on where in the plot the subject sits. They also ignore the\n",
    "  selector entirely: `BAL` is the basal area per hectare in trees larger than\n",
    "  the subject **over the whole plot**, so it does not shrink when the search\n",
    "  radius does.\n",
    "- **`SBAr` and `Almdg`.** A ratio and a weight-normalised mean of areas do not\n",
    "  scale with the observed share -- halving the competitor set leaves both\n",
    "  unchanged -- so dividing them would invent competition. `INDEX_REGISTRY[name].additive`\n",
    "  says which indices the correction applies to.\n",
    "- **Selections with no zone fixed in advance** (`SearchCone`,\n",
    "  `NearestNeighbours`, `BitterlichBAF`). Their reach is an outcome of the data,\n",
    "  and a neighbourhood truncated by the boundary reaches less far -- so using it\n",
    "  to size the correction would under-correct by exactly the amount at issue.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c23",
   "metadata": {},
   "source": [
    "## 6. Where the numbers come from\n",
    "\n",
    "Every index carries the citation of the paper that proposed it -- not the review\n",
    "they were collected from. That means a result can always answer \"according to\n",
    "whom?\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "c24",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-26T12:01:50.379600Z",
     "iopub.status.busy": "2026-07-26T12:01:50.379600Z",
     "iopub.status.idle": "2026-07-26T12:01:50.389345Z",
     "shell.execute_reply": "2026-07-26T12:01:50.388840Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " index  spatial                                       author  year\n",
      " BA-gj    False                Steneker, G.A. & Jarvis, J.M.  1963\n",
      "   BAL    False  Wykoff, W.R., Crookston, N.L. & Stage, A.R.  1982\n",
      "   Sdr    False                                Lorimer, C.G.  1983\n",
      "   drg    False                               Hamilton, D.A.  1986\n",
      "   BAr    False                     Corona, P. & Ferrara, A.  1989\n",
      "  BALr    False                                Vanclay, J.K.  1991\n",
      "BALMOD    False                 Schröder, J. & von Gadow, K.  1999\n",
      "    Sl     True                               Staebler, G.R.  1951\n",
      "   SOr     True                                Gerrard, D.J.  1969\n",
      "  SOdr     True                                  Bella, I.E.  1971\n",
      "   Heg     True                                    Hegyi, F.  1974\n",
      " SAng1     True                                    Lin, J.Y.  1974\n",
      " Almdg     True                                Alemdag, I.S.  1978\n",
      " Sdrl1     True                                Lorimer, C.G.  1983\n",
      " Sdrl2     True                      Martin, G.L. & Ek, A.R.  1984\n",
      "  SBAr     True Daniels, R.F., Burkhart, H.E. & Clason, T.R.  1986\n",
      " SAng2     True               Rouvinen, S. & Kuuluvainen, T.  1997\n",
      "SdrAng     True               Rouvinen, S. & Kuuluvainen, T.  1997\n"
     ]
    }
   ],
   "source": [
    "catalogue = pd.DataFrame(\n",
    "    [\n",
    "        {\n",
    "            \"index\": name,\n",
    "            \"spatial\": entry.spatial,\n",
    "            \"author\": index_source(name).author,\n",
    "            \"year\": index_source(name).year,\n",
    "        }\n",
    "        for name, entry in INDEX_REGISTRY.items()\n",
    "    ]\n",
    ").sort_values([\"spatial\", \"year\"])\n",
    "print(catalogue.to_string(index=False))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c25",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "- `Stand.impute(\"height_m\")` fits Näslund's curve to the stand's own measured\n",
    "  pairs and fills the gaps, keeping measured and modelled values distinct and\n",
    "  recording which model produced each modelled one.\n",
    "- `competition_indices(...)` computes any of the eighteen indices; pass\n",
    "  `indices=` to choose and `selector=` to say how competitors are found.\n",
    "- The selection rule is part of the result. Report it alongside the index, and\n",
    "  treat its parameters as parameters -- the defaults are one study's choices.\n",
    "- Distance-dependent indices are biased low at the plot edge;\n",
    "  `observed_zone_fraction` reports the exposure and the default corrects for it.\n",
    "\n",
    "### Next steps\n",
    "\n",
    "- The influence-zone indices (`Sl`, `SOr`, `SOdr`) additionally need a crown\n",
    "  radius. Supply one with `stand.impute(\"crown_radius_m\", <callable>)` or the\n",
    "  `crown_radius=` argument; a caller-supplied value is recorded as uncited.\n",
    "- `SearchCone(apex=\"crown_base\")` needs `crown_base_height_m`, which can be\n",
    "  imputed the same way."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pyforestry",
   "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.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}