{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a2745dfd",
   "metadata": {},
   "source": [
    "# Fit your own xG — a reproduction of Lab 1 (Soccer analytics: maths and statistics)\n",
    "\n",
    "**Teaching artifact by Ahmed Beddi, on public data.** Shots from the 2018 FIFA World Cup, StatsBomb open data (competition 43, season 3), exported to `wc2018_shots.csv`. Data provided by StatsBomb under its open-data user agreement (https://github.com/statsbomb/open-data).\n",
    "\n",
    "This notebook reproduces the headline numbers of the course's Lab 1 (Sessions 3–5): a two-feature logistic expected-goals (xG) model built from shot geometry, judged against StatsBomb's own xG, and the penalty-spot gap that shows where geometry stops working. The interactive xG explorer on the course page uses the same formulas and the coefficients fitted here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2c3d4d5f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:48:13.741488Z",
     "iopub.status.busy": "2026-11-03T14:48:13.741245Z",
     "iopub.status.idle": "2026-11-03T14:48:15.423778Z",
     "shell.execute_reply": "2026-11-03T14:48:15.422568Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1706 shots; 183 recorded as goals\n"
     ]
    }
   ],
   "source": [
    "import numpy as np, pandas as pd, statsmodels.api as sm\n",
    "import matplotlib.pyplot as plt\n",
    "shots = pd.read_csv('wc2018_shots.csv')\n",
    "own = pd.read_csv('wc2018_own_goals.csv')\n",
    "matches = pd.read_csv('wc2018_matches.csv')\n",
    "print(len(shots), 'shots;', int(shots.goal.sum()), 'recorded as goals')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5004637b",
   "metadata": {},
   "source": [
    "## Lab 0's mystery: 183 goals in the shot file, 169 in the results\n",
    "Shoot-out kicks (period 5) are recorded as shots; own goals are not shots at all."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "af4abcc1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:48:15.425613Z",
     "iopub.status.busy": "2026-11-03T14:48:15.425421Z",
     "iopub.status.idle": "2026-11-03T14:48:15.431918Z",
     "shell.execute_reply": "2026-11-03T14:48:15.431004Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "in-game shots: 1667 | in-game goals from shots: 157\n",
      "own goals: 12 | total from results: 169\n",
      "base rate (in-game): 0.094\n"
     ]
    }
   ],
   "source": [
    "in_game = shots[shots.period < 5]\n",
    "print('in-game shots:', len(in_game), '| in-game goals from shots:', int(in_game.goal.sum()))\n",
    "print('own goals:', len(own), '| total from results:', int((matches.home_score + matches.away_score).sum()))\n",
    "print('base rate (in-game):', round(in_game.goal.mean(), 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5b6ebcf",
   "metadata": {},
   "source": [
    "## Geometry: distance and the angle the goal gives you\n",
    "StatsBomb coordinates: the pitch is 120 × 80 yards, the attacking team shoots towards x = 120, and the posts sit at y = 36 and y = 44."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3af2a2ee",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:48:15.433447Z",
     "iopub.status.busy": "2026-11-03T14:48:15.433248Z",
     "iopub.status.idle": "2026-11-03T14:48:15.444963Z",
     "shell.execute_reply": "2026-11-03T14:48:15.443576Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "median shot distance: 18.4 yd; median goal distance: 11.0 yd\n"
     ]
    }
   ],
   "source": [
    "def geometry(df):\n",
    "    d = np.hypot(120 - df.x, 40 - df.y)\n",
    "    ax, ay = 120 - df.x, 36 - df.y\n",
    "    bx, by = 120 - df.x, 44 - df.y\n",
    "    cos = (ax * bx + ay * by) / (np.hypot(ax, ay) * np.hypot(bx, by))\n",
    "    return d, np.arccos(np.clip(cos, -1, 1))\n",
    "in_game = in_game.assign(dist=geometry(in_game)[0], angle=geometry(in_game)[1])\n",
    "print('median shot distance:', round(in_game.dist.median(), 1), 'yd; median goal distance:', round(in_game[in_game.goal == 1].dist.median(), 1), 'yd')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7ce39d3",
   "metadata": {},
   "source": [
    "## Fit the model on open-play shots\n",
    "$z = b_0 + b_1\\,\\text{angle} + b_2\\,\\text{distance}$, $\\;xG = 1/(1+e^{-z})$. Unpenalised maximum likelihood (statsmodels `Logit`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ff42693f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:48:15.446507Z",
     "iopub.status.busy": "2026-11-03T14:48:15.446344Z",
     "iopub.status.idle": "2026-11-03T14:48:15.457836Z",
     "shell.execute_reply": "2026-11-03T14:48:15.456898Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "const   -2.0265\n",
      "angle    1.2662\n",
      "dist    -0.0689\n",
      "dtype: float64\n",
      "open-play shots: 1556 | goals: 129 | sum of fitted xG: 129.0\n",
      "Brier, our model: 0.0692\n",
      "Brier, StatsBomb xG: 0.0638\n"
     ]
    }
   ],
   "source": [
    "op = in_game[in_game.shot_type == 'Open Play']\n",
    "X = sm.add_constant(op[['angle', 'dist']])\n",
    "fit = sm.Logit(op.goal, X).fit(disp=0)\n",
    "print(fit.params.round(4))\n",
    "p = fit.predict(X)\n",
    "print('open-play shots:', len(op), '| goals:', int(op.goal.sum()), '| sum of fitted xG:', round(p.sum(), 1))\n",
    "print('Brier, our model:', round(((p - op.goal) ** 2).mean(), 4))\n",
    "print('Brier, StatsBomb xG:', round(((op.shot_statsbomb_xg - op.goal) ** 2).mean(), 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c58314ef",
   "metadata": {},
   "source": [
    "## Calibration by distance band\n",
    "Fitted mean xG against the actual conversion rate, and the illustrative teaching coefficients used in Session 3 ($b_0=-0.6$, $b_1=2.2$, $b_2=-0.1$) for comparison. The teaching model has the right shape and the wrong level: that is the point of Session 5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a7809ce7",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:48:15.459653Z",
     "iopub.status.busy": "2026-11-03T14:48:15.459471Z",
     "iopub.status.idle": "2026-11-03T14:48:15.627344Z",
     "shell.execute_reply": "2026-11-03T14:48:15.625834Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "           shots  goals  actual  fitted  teaching\n",
      "dist                                             \n",
      "[0, 6)        53     20   0.377   0.384     0.869\n",
      "[6, 12)      343     55   0.160   0.158     0.526\n",
      "[12, 18)     375     30   0.080   0.077     0.247\n",
      "[18, 24)     318     11   0.035   0.045     0.120\n",
      "[24, 30)     294     11   0.037   0.029     0.065\n",
      "[30, 120)    173      2   0.012   0.017     0.030\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxYAAAGGCAYAAADmRxfNAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjIsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvgI3uAAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAX+ZJREFUeF7t3Qm8TeX7///LPGXIlFllyBSSeSpCIkkaVCqVlKGBUkSUBvRRmgfKEMqYDJUhJJlKIYrMQ6aQzFPs/+N9/75r//eZz7EP9vB6Ph77wTln77XXutc+rGvd13XdaXw+n88AAAAAIAhpg3gtAAAAADgEFgAAAACCRmABAAAAIGgEFgAAAACCRmABAAAAIGgEFgAAAACCRmABAAAAIGgEFgAAAACCRmABAAAAIGgEFgAAAACCRmABAAAAIGgEFgAAAACCRmABIKTt3r3bjh07luQ+njx50v75558kn5cYbWPHjh3m8/mC2k6k0FhoTAAASA4CCwAhaciQIVagQAErWLCgZc+e3W644QZbt25dnH39+++/rWnTppYzZ07LmzevlS9f3ubNmxfneQpQGjdubBkyZLDmzZvH+PmWLVusRYsWli1bNitSpIjlyJHDXnjhBfvvv/8SHJvvv//ebWvKlClxnvPII4+4n+kYYnv++ectc+bM9u+//wY97tq+3ufgwYOJbuuXX35xz5s7d26S77lz5067/fbbLXfu3G4sNCa1a9dO1mtxbvQ5f+65587txedBau/Prl27bN++fedhTwGEGgILACHns88+s0cffdTatGlj8+fPt7Fjx7oL3htvvNEOHTrk398jR45YnTp1bM2aNfbxxx/b9OnTrWjRonbTTTfZggUL/M+bOXOmVapUybZu3WqXXXaZnT171v+zEydOWLNmzey3336zkSNH2o8//mjPPvus9e/f33r37p3g2NSqVcsyZswY7wX37Nmz3azHd999F+f1c+bMsWuuucZy5coV9LjrPRILfgLpeUnNxCxbtsyNk8a8S5cuNmHCBHvrrbfcjFGjRo3c10h9Z86cifGZvNhSe3/0+/X4449f7MMCcAGkvwDvAQAp8uabb7qLEV3UejQToce3335rd911l/u2ApBNmzbZqlWrrFy5cu57mtmoWrWqvfbaa+650qNHD3dh/NFHH1nDhg1j7MuSJUtcYPL111+79xQFK9u2bbNhw4bZgAED4t33TJky2XXXXRcnsNi4caMLYFq3bu1mTnQxnyZNGrcNBUWaPejVq1ecbSrlSNtMyunTp10wpFmcxCh40gWiZiqSQ/vZsWNHF/BoNqZw4cL+lynI04xOu3btXHCnGZ1zpX3XjE1Kaf907ArmkpLc90ju85J7bhLjXawn93x4Dh8+nOS5PlcKGLNmzZrky5Nz/Ml5TkKSux8AQh8zFgBCzoMPPmg9e/aMsV8lS5Z0F+jbt2/3f18X/g0aNPAHFaKLc71esxSqEZBBgwbZmDFj4r1A81KSChUqFOP99LUCgcDZjdiaNGliq1evtj179vh/pBkJXfx26tTJpX9oJsTzww8/uAtMvc4zceJEu/rqq13KkVK57r333jhpI0pNeeONN9wsQp48eeyhhx5K8JxpX1q1auVSw3TRqIAgcJYnIV999ZWbsdBYBQYVoovh119/3a688ko3oyOaHVJwEzs9TTM9gbMx77zzjjsuBVw1atRw+5QvXz43PsePH0/ys6fz3bJlS7fNSy65xAWNsWeClEakmZZff/3VKlSoYFmyZHHnTwGcxjuQZrk6dOjgZq60L/pcKUAN1L59e5c2p2Dz2muvdc/TORg6dGiS++tRYKrx0TaUGqfzq6/1ef3zzz8T3Y5qhXSO9XwFcTqX3bp1838W33vvPfez9evXx9iOPqulSpVyM30J0cyV0vF0TrRPGleNX+xxUiDXr18/K1asmBtPje/SpUvjbFbBepkyZdxzlLao4PTo0aPueQr4tZ/6HdDnXH/v27ev+9mpU6fs6aefdil32g/9qZ8lNasGIMT5ACAMTJw4URXVvoULF/r3NmfOnL7evXvH2ftFixa5586fPz/Oz6pWrepr2rSp//v//vuvL3fu3L57773Xd/LkSff9bdu2+YoXL+67++67Ex2Z33//3b3P2LFj/c+78847fTVr1vSdOHHClzlzZt8bb7zh/1nXrl19OXLk8J0+fdp9b+jQoe71bdu29X311Ve+Dz74wFeoUCHfFVdc4Tt69Kj/dZdccomvWLFiviZNmvg+/fRT3/fff+9+9tFHH7nX6xhE71miRAlf/vz5fe+//75v6tSpbttXXXWVe953332X4PE8//zzvmzZsvnOnj2b6DF7tL/a5tq1a2M8/8UXX/SlT5/e/71Bgwa5rytUqOB77rnnfF9//bXv5Zdf9mXKlMnXvHnzRN9r7969vrx58/pKly7tjnvSpEm+W265xZcmTRrflClT/K996qmn3DnUsb/yyivuPZ599llf2rRpfY8++qj/eWfOnPFVq1bNV7BgQTfWep5+ruMYOXKk/3kaM53/ypUr+0aMGOH78ssvfXXr1nXvu3Tp0kT32bNlyxa33auvvtp34403us/v8OHD3bHoc7t7927/c3WMzzzzjPta+1i9enX3OenTp49v+vTp7lj03ho30Ws1pn379o2xLzq/es/JkycnuI/6fdHnsn///u73Q+OVLl06X79+/WLsz+WXX+5r0aKFO3Z9TjVmhQsXjvG51Pvr/bp06eKbNm2a+6znypXLjbGO4+DBg27/9XmuV6+e+/uaNWvc++ic6fOmz4c+zzpWHeObb76ZrPEFEJoILACEPF2k5cuXz3fzzTf79/Xw4cPuoubdd9+Ns/+bNm1yPxszZkycn8UOLGTZsmW+MmXKuPfQxaQuvO644w7foUOHkhwbXWx16NDBPU8X5bpIfOGFF9zXjRo18jVr1sy/jUqVKvluvfVW97UCB72fApFAf/75p7vQe/XVV/3fVmBRo0YNd7EWKHZg8d5777nj/umnn2I8TxfjSQUWt912m9u/5EpJYKH3jn0R/Mknn8QJFGN74oknfFmzZvXt2rXL/yON8Q033OCCCG88dJGq99AFcCAFORrLjRs3um+PGjXKlyFDBhcQBlIgoWDMo68VlHivk/3797tgyDu3SfECi4oVK/r+++8//9N1LFmyZHFBpicwsFDQqYAwdlDcoEEDX506dfyvadmypRuDQO3atXOfqVOnTiW4e7rob9y4cYyfz5gxI8aYaH/Kly8fY7+9wN7brx07drjfk27dusXYln4eO1DT71SbNm1iPK9s2bK+Vq1axfieAj19/gGEL1KhAIQ05Zh7HZuU+uTxipbTpUsXZ/+97yknPylKyVD6z4YNG1z3I9VXqJbjm2++sfHjxyf1cpfWpPQnWblypUtjUhqN6GdKf9K+7t+/36WEeGlQShPZu3evS8sJVLp0aVcH4m3Tc/3111vatIn/k600JaUCVatWLcY277vvviSPQ+OQ0vz/lLj//vvj7FP69Olt4cKFCW5GXaxuueUWl4bkUTqcxsyrZQn8ftu2bWNsSylxSvHxUng0pldddZX7fCgdyXvUrFnTdRdTgwDP5Zdf7lK/PErV0fdUP+O55557XKOAwMe0adNi7INS2wI/ozoWdTFL6Lg1Jvq8169f39WA/PHHH26bSo/SMXt0bPpaNUKitLJJkya5cU7sPCrtTjU0r776quuGJqqb8WqUPKofCtxvpYSJd/yLFy92+xf786v9Llu2bJzPb2zajxkzZrgUOy+9UTVO+vwDCF8EFgBClgpCdWGpWgkVYis336PccOW+64IwNq/mQe1Sk/L222+7/G9ddKrOQPnrqjXQBZOKlteuXZvoJhQo6AJPF0cq1lYdgC5URT9TTv/PP//sOi0pf9wLLFR0LsWLF4+zfV3QBl5ESnIu+nXRF98xJ2cclCev4Op8iV23oToU1Tls3rw5wbfUGCU0PhI4Rqo9iV0UruNWMOa9h7anmhgda+BDtSsSWLMQX1GzzkFgzY2CDdUqBD5i70Ps4xbVLSR23Apq9RlS7YGCXNUtxK6TUctkjd/o0aPdNvXZVRCugCMxCqLV/ODll1+2K664wj3U/Sx2vUvs4/c+f97xp/TzG5t+zxRA6b01Hqp1UU0H66YA4Y3AAkBI0p1mFaHqIl9BhS4EY9NF54oVK+J8XzMHEt9FaWzatgqCdac1kIpltQ+6q5oYdZvSXXTdwdVdaM0seBdhFStWdBd/mklYtGiRu4grUaKE25y+L/Et6qdZj8C79Mmlbca3XkBy1hDQHWQVsid0J11Byx133BFnjZDYxe2a+YhP7H1QkKVZnMSOU8eT0PhI4Gu177ELkPVa7Z/3vPz587tzrTVN4nuohXBK6OJ41KhRMR7ebJUnofOR0HErkFUwrSBanz0FC3/99Zf7XiDNbGiGZty4ce64P//8c1ccr0AkMZdeeqnbX82Wafv6/A4cODDJgCS2YD+/Oj7tu/ZDBe5169a1l156yTp37pySUwAgxBBYAAhJmjHQhc/UqVOtevXq8e6jAg8FBrpI9OiCVetRKHVDd0GToo47usCN3Y3Gm/XQzxOju9VVqlRxgYWCh8ALSwUcunBTYKGfK+XEowt5BSA6vkCa4dDsiZd6khLapgItXYgG0h3wpOgutmZK1KknvpXOdTd58uTJLs1F1AFIYs9yqJ1ufGLvgzo7KZVGwVdCNK56XeyUNt2dV4chdUDyKN1s1qxZMTblpSV576H1Q5RaJLow9h6aDdC4J6eVbUrFPm7dkdc6Jwkdtz4nChTUgUqfJc2ASXwzZ+ocpcBFY6nPTGLdwkRB1k8//eRSwDR++jzqfW677bYUL4CosZTYn1999vQZiP35DQxAdXw6Ts046fdLKVAjRoxwKWLJWcQRQAgL3/IQAJGqe/furgC0Z8+eviVLlsR4BHYhUnccddipUqWK6wSlwk8VUuu148aNi3d4Yhdvq+uNnv/AAw/4Vq9e7fv77799s2fPdsWrKugN7N6TEO1nkSJF3Ha8rjcedRu69NJLXeGv3iuQio5VAKvnHDlyxLdhwwbXsUrFvepM5VHxdnzdr2IXb+/Zs8c9t1atWm6ctE0VSV922WVJFm+LimczZszojv3zzz/3rVq1ynUaat26tXv9wIED/fukbavIV0W4ixcv9m3dutV1LVL3oNhdofTaK6+80hUJq3OVzpW6Dum1gYXNsa1cudIVX6vwV+Oh7lcq1tf2BwwY4H+6xlH7rQJ8bfv48eO+mTNnuv1T0bNHr1dxszo86djUBezHH390BfgqhPa6dal4W12sYtP39LPk8Iq39fnUuBw4cMC3fft23+233+66H61YscK/mcDibe2PxkvF+3qNxlnjrnEoUKBAnLdWUb/2S58ZrztYQjTWRYsW9ZUrV8591jVOc+bMcd9T5ypP4P54/vrrL7df6pLl0bGoC5Q6dB07dsyNqTpa6fdGHaE8amCgzlD6fKshgsZd50G/t/qc6rX6bOhzqs8agPBFYAEgpOjCUxcwCT0CO0OJLj7VGtR7vi6whgwZkuAxxdcVSs/XhU7ge6p7jrpFJce8efPcaxVcxLZz5073M10Y6kIxkC7svHanutjU83SRq+0FSm5gIQoEdGGobWmbuhhVl6GkukJ5fvnlF9ciNXAsdMGnICY2tRhV0OS9l9rHvvbaa/EGFt98841rL+odp4KK2EFYfNR9KvA91K2pR48e/iBAFFjoOL2uT957qO2vLogDLV++3L23t8/6U18HdopKzcBCbX8VqHn7pM+n9jNQ7Av5J5980j1fnxn9Wb9+fd/jjz8eb2DhBZf33Xdfcj6qLnDRWOm4NZb6U78TgYFscgMLfZ7VTSxwLNXaWGMc6IcffnABiBfAi4J373fO24/atWsnK5AHELrSaNdCeEIFQJTRP0leYWh8VFTqpeEEUuHtgQMHrHLlyommtChVQ91uYm9DKSrahnK+lRKUnPoMj1I71GHHW5gtNhXqKi1KXYXio1QupTCpBkBdnWKn5Gg8VKyuYuFASuHR/qp2I7BjlApxly9f7hYqU92AFg3Ugm1aNE77mBxKD1Ohs/ZJKWUJFY8rbUqpNSrAVQH1wYMHXXqOV0uihf2eeeYZl6qkmgcVsmvslbalOoHk0PGoo5aOV+lRep9AXbt2tbFjx9quXbvc+dVYakyUthVfJy2dLy3sp/Ot4mpt01sdXZQGp1qRokWLxngfFeh7RedJUU2Kzrf2S6lGSv3RGKjGI/bK5fp86HuBx6X904rwGid9HnU+VG8R+zOkFCil26lxgOp7kkPHppQobVPjFDstK7798T7j+jzEXvVdny19VvQZ0JjH16lN+650LtWWeOOqVDjth+pj9HlJqj4EQOgjsAAAnDeBgUV8F5ypITCwCBWBgYXqV84XrUquAETBXWBwBAAXQ/JuFwEAgJDRvXt3V8Su2SwV1RNUAAgFBBYAgPNGqT9KhTqfF74NGjRwbVRDiVKbdNxakO98UNrUzTffbK1atYrTKhkALhZSoQAAAAAEjXUsAAAAAASNwAIAAABA0KK2xkKrgO7cudO1zTufub8AAABAuFIbeLWMVsvy+Fp4B4rawEJBRewe5QAAAADi0lo+RYoUSXRoojaw8Bb40SDFXqwIAAAAgLnFSXUzPvbimPGJ2sDCS39SUEFgAQAAACQsOaUDFG8DAAAACBqBBQAAAICgEVgAAAAACFrU1lgAAADgwjhz5oydPn36wrwZUiRDhgyWLl26FL0mIQQWAAAAOC+0BsLu3bvt33//PS/bR+rIlSuXFShQIFkF2okhsAAAAMB54QUV+fPnt6xZswZ94YrUpcDv2LFj9vfff7sNFyxYMKg3ILAAAABAqlP6k2YqFFTkyZMn1beP1JElSxa3IQUXOlfBpEVRvA0AAIBU59VUaKYCoc07R8HWwRBYAAAA4LxJzsJquLhS6xwRWAAAAAAIGjUWiGNK/wqpOiote65O1e0BAAAg9DBjAQAAACBoBBYAAAAAgkZgAQAAAISw5557zmbMmBHCe/j/UGMBAACAC+rydu0v6PttGfHJBXu/devW2bXXXmt//fWX5cyZM1Xed+7cuVasWLFU2db5xIwFAAAAEIsW+Fu9erVt3rw53rHRz3///fcYP9cq1lu2bLEjR47Ypk2b/Ctab9++3Q4dOhRjO/r5iRMnYnxv37599uuvv7rthCMCCwAAACDAyy+/7GYbGjZsaCVLlrTmzZvbyZMn/c+YMmWKW6W6bt267uf169d3QcSiRYvs3nvvdc9r2rSpvfbaa+7vt9xyi33xxRcxxrh8+fK2dOlS973169e7WY5ChQpZkyZNrGDBgvb111+H3TkhsAAAAAD+j1KYvvzyS5s0aZILFnbs2GELFy60iRMnumdoBqNNmzY2aNAg++eff2z37t1uxepnnnnGGjVqZD/++KN7nlKi3nrrrWSN69ixY61UqVK2d+9e956dO3e2p556KlmvDSXUWAAAAAD/p0iRIrZ8+XL/ePh8PjeDsHLlSjcbMXr0aKtUqZI9+OCD7jn58uWzmTNnWtasWc95DF944QX/azUzon3YsGGDHT161LJly3bO273QCCwAAACAAEOGDLGvvvrKFixY4L+w16yE6IK/QoWYiwnnyJEjqPHTNjW78d1337kajSuvvNJt77///gtquxcaqVAAAADA/1EthFKRNDuhomulOtWuXds/PpqhUApUSqRJk8aOHz/uf4kKv71ARVTDoRSr8ePHu1kKpUaFIwILAAAA4P/8+eefbkZCgUWuXLlch6aff/7ZPz433HCDS31SLYbnxRdftF69erkv06VL5/4M7OykVrGB61DMmjXLBReiWQl1iOrQoYNVrFjRvf77778Py/NBKhQAAADwf9TBSd2cWrdubVdddZXNnj07xnoUN910k912221Ws2ZNu+OOO1xHp3nz5rlOUVK8eHErXLiwtW3b1u6++25r3769CxpatGhhVapUcd2gNIORO3du9/z06dNby5Yt7cknn3SF3yreVmpUOCKwAAAAwAV1IResSyld/GuGQilRKpweOXKkS4nKkiWLf1OjRo2yqVOn2uLFi61GjRo2ePBg19VJMmTI4NrO6vXea5o1a2ZLlixxNRSqn1Dw8tlnn7lWtTJmzBgbPny4my1R4KI2t4Gv79atm5UrVy6lh3LBpfGp1D0KaZES9Sc+ePBg0AU3kWZK/5gFScFq2XN1sJsAAABhRou/qW7giiuusMyZM4fZ3keXxM5VSq6ZQ2rGQgelnLU1a9a4aLFx48Zuqighyl1TT+GtW7damTJlrFWrVm46CQAAAMCFFTLF24qGtNKgpnp+++03u//++13RjFfYEpsq9DUlpEKZVatWWceOHa1OnTpxlkYHAAAAcP6FzO399957z/bv328rVqyw7Nmz27Zt29ysxeTJk+3222+PMxLq9at2X1oKPW3atHbkyBGX2/b555/bQw89dP5HDgAAAIBfyMxYqGBFPYMVVIjacqmSftiwYfGeLq1KqPQnBRVyySWXuCp8pVPFR8/XrEjgAwAAAEDqCInAQulOKu6pVq1ajKOqWrVqgu227rnnHtcDWAGJev8OGjTI1VqoXVd8+vfv74q1vUfRokVTZwQBAAAAhEYq1N69e11wkTdv3hinRKlOO3fujPc0KehQT2Av7UlF3kOHDnV9g+PTs2dPV7/h0cwFwQUAAACQOkKixkKzCHL48OEYR6WL/0svvTTeI33++edtwoQJbmVCrY64cOFCa9eunal7rhYiiS1TpkzuAQAAACD1hUQqlBYP0eqDsVcZ3LhxoxUpUiTeo/7ggw/cqojXXXed5cmTxy008vTTT7sicAAAAAAXVkgEFqLaCK1s6Pnvv//cioOBNROBa/lpVUOtYxFIX+v7AAAAAC6skEiFkmeeecZq165tDz/8sDVq1Mgtba7UqEceecT9fMaMGa5LlJZCv/baa61Lly6uZmLlypWu3uKHH35wgcmQIUMu7AgCAAAACJ0ZCy12t2DBAjt9+rRbo0JF2IsWLXJpTqKUqPr161v+/Pnd13379rV3333X1q1bZ2+88YYrAP/qq6/cwnoAAAAALqw0vsD8oiiiwnAVjR88eNBy5MgRRUeetCn9KyT9pBRo2XN1Cp4NAAAigdYW03ICV1xxhWXOnDnsDkmZM8qU0c1tLYmg2l+tmZY+fcgk/KSaxM5VSq6ZI29kAAAAENJGN72w64m1nbE9ReOxfPlyq1OnjgsklKav4KhkyZL28ccfu+UO9u3bZ3v27LHy5cu77R4/ftz++OMPq1KlilsC4Vzdeuut1rBhQ3viiSfOdRMXFYEFAAAAEGDy5Ml29dVX29KlS/3f1RIHqvOVKVOmuNT9VatWua+Vmq+aX935D2Z5A2XSKEgJVwQWAAAAgJlbsHnEiBEuoNDfP/30U6tVq5arBd6yZYuVKlXK1qxZ4+qC//nnH/fzyy+/3DURkuHDh7uARLMdsmPHDpszZ47rXNq4cWMrUaJEjHHevn27zZo1y3U+vf3228P+HIRMu1kAAADgYjp79qyrq9i6datLddLf//rrL7dL6kiqoELrrint6ejRo+7nf/75p39mY+7cubZ27Vr3fAUMWsR57NixNnPmTKtYsaLreurRc5VKpWBk6tSpdv3119u///57MQ8/aMxYAAAAAGZuPTSto9a9e3c3Q6G/x6blDzQDoVQo7+eaoVDwMGrUKJcKdfLkSWvfvr0NHjzY2rVr5zah+ozOnTu7RZ2zZctmjz32mHXs2NEGDhzofv7tt99as2bNrE2bNmF7LggsAAAAgFS0YsUK+/vvv6106dL+2QylQR05csRWr15tl112ma1fvz5GkfZNN93kisTDGYEFAAAAkIpUW6G0Ks1QBFL9hQq0vdUeFGAE8tZrC1cEFgAAAEAqUptaWbJkSbxdory6jU2bNrlZDVEgotqOcEbxNgAAAJACWbJksb179/pnHvS1KP1JKlWq5Na9ePHFF/3PUTF37dq13YyFt+hev379XEAhQ4YMsd27d6dgL0IPgQUAAACQAo0aNbK0adNasWLF7JNPPnFBhIKGa665xrp27epW51YHqPHjx7vaCgUaN998szVt2tRy5szp3kmvU5tatatVx6hJkya5xfHCGalQAAAAuKBSuhL2habg4NSpUzHedtGiRf7iatVCqM2sirSvvPJKF2TMnz/fli1bZrly5XKvU5ChtrQrV650rWv1tWYqPAom1L72559/dt+qV6+eS4W65JJLLvThphoCCwAAACBAoUKF4oyHCq8DZc+e3QUDHs1S1KxZM8Zz1Hq2evXqCY6t2s5q/QqPgpRwRioUAAAAgKARWAAAAAAIGoEFAAAAgKARWAAAAAAIGoEFAAAAgKAF1RVq3bp1tn//fqtVq5YNHz7ctdgqV66c69PrrTgIAAAAIPKd84zF22+/bWXLlrU5c+a4UdICH7Nnz7bBgwdbhQoVbNy4cZE/egAAAACccwostJDHM888Yy+//LL17NnTP5T33nuvrV+/3nr16mXt27e3HTt2MMwAAABAFDinwEJLkN955532/PPPW7p06WIMU5o0aaxHjx5uefJp06ZFwRACAAAAOKcai40bN1qzZs1ijN7TTz9tWbJk8X+vfPny7nkAAAAAIt85BRZXXHGFS3kKpLqKQH/++ac1aNAg8kcQAAAAwLmlQjVv3txGjhxpK1asiHcIP/vsM1u7dq01atSIIQYAAAASsXXrVtu7d2+Cz9A19/Tp00N+DM9pxuLWW291hdpqM9uxY0erV6+eFShQwA3KlClTbPz48fb+++9biRIlQn4AAAAAcGFN6R8z0+V8a9lzdYre4uOPP3aZN6VLl07R687VY489ZjVr1rS+ffvGuwntz8KFC92SDqHsnAILGTp0qKuzePPNN93sxT///GMFCxa0a665xhYvXmzVq1cP5eMGAAAA4vXKK69Y3rx5L1hgkZQPPvjATp8+ndTTLrpzDiykVatW7iGnTp2yjBkzXvQDAgAAAM6FLt4feughd8Nca7bNnz/f3nnnHbcppfkPGTLENm/e7G6kd+3a1bJnz+5/mwULFtjUqVNdnbEWjFZjo3z58vl/vmnTJhcgqE75qquusieeeMKKFCkSYzeV9TN58mS33e7du1upUqXcz7/66ivbsGGD+96ePXvcsg9vvPGGvfvuu/b7779bnTp13PYyZMjg396ECRNcJlHatGntySeftJkzZ9q1115rN95447kMTbKcU7tZz7///usGWWIHFTpoTdkAAAAA4UAX4XXr1nXXtVoI2svA+eWXX1yq0n///ecydhRw6GJeN9Zl4sSJdsMNN9jRo0ftpptuckFGYAfV3377zSpXrmy7d+92z1u+fLnbtr72KAiYMWOG+/5PP/3kUrGOHTvmfqzXz5071/39yJEjNnr0aOvSpYvlyJHDihUrZi+++KL169fPv63+/fvbww8/bCVLlrRKlSq5v2u5CO+6/XwJasZCB1alShU3XfTUU0+5kyGKqjp06OAW0NOgAwAAAKFO67M9+uij7tq2SZMm1rp1a7fLWhBaMwLexfs999zjSgDGjBljDz74oFWrVs1mzZpl119/vXu+/tSshRaLLly4sFv7TTXKanAknTp1cuu+/f33365OWbS9YcOGub/fdtttbk04LUp93XXXxTtsel8FMaJASNffuvY+dOiQ2/9Ro0a57YgaL2kpiPMtqMBC0zeqUH/kkUfctI2mijQlo4hLtRcPPPDA+d5/AAAA4LxatGiRW25h4MCB/vfRdfDKlSvd18WLF7dcuXK5a+Bff/3VzViIZiQUWKj+eMSIEf7Xpk+f3gYNGhRjn3Wz3qPtFS1aNMaMRmyBz9csy/Dhw91T/vjjDztx4oS1aNHC/5IyZcpckKZKQaVCScOGDd0BaKpFB6hF8datW0dQAQAAgLCnuovjx4+7mQAFCt5DtQrexb2aGmnmQTUUoswdj8/nc1k+2bJlS3QsFGwE0uxJYgKfH/hcBRX62ssk8gTWX5wvQc1YyMGDB930jmYsnnvuORctKefrrbfesssuu+x87TcAAABw3iggEF2QawkFpUYp/chz+PBhf/G2MnVUUK1UJNm1a5f/eWnSpHHF2suWLXM35D0//PCDXX311XbppZem6jGoNuTMmTP2448/+tOodu7c6Yq/z7egZixUMa98rdWrV7silAEDBtiqVatcVKe8Mg0YAAAAEE7y5Mnj0p/27dvndlu1xC+88IK7WNd1rkoB8ufPb9OmTXM/V2taZfDs37/f1U2oY1Sgbt262euvv25z5sxxMx9Ki1Jxtjo6pTbd2NeMidac00zK2LFj7e6773apWudbUDMWmmrRwOnhTbdokFU8ogFTu6369euf72MAAAAAUo1mHlRDrGJqBRdaEFodmm6//XbX+VQX6b179/bXMajO+K677nLtZVV7oa5Mau/qpSip7vjAgQOu6FuBhwq1leWj2ghRmlTmzJlj7L86PnnpS/qZl0qla269v2ZCPEqL0vM97733nmtVO27cOPczdY169tlnU218EpLG583zRBlFizlz5nSpXIEnApbqq2GmdLVLAAAQ/nQDWms+qOg59kVzuNBlcuAFvOgaMqFrR11X6voyMYEpVOeLajoUlGTKlMn/FmpL+9prr1nbtm3jvG1i5yol18xBzVgAAAAAkSp2UCGJ3ZBOKqiQwEX1zpf777/fpWz17dvXrb2hNrcKepR+dT4FVWMBAAAAILR8+umnrsXs448/7uortJ7G7NmzXevb84kZCwAAACCCqNPU4MGDL/gRMWMBAAAAIGhBz1ho0RBVnKtdlgo/At1yyy3nPZcLAAAAwMUXdGDRvHlzmzt3rtWrV8+tOBjo6NGjF/8IAQAAcNFEaQPSsJJa5yiowGLTpk2uEGTBggX+PrwAAACAtwaD1n/IkiULAxLCdI7EO2fnKqjAImvWrK44xFvcAwAAABAtDqeF3LQgnOi6Mb72rbh4NFOhoELnSOfKW9DvXAUVWCj16cYbb7QZM2ZY06ZNz3UfAAAAEIG8NHkvuEBoUlARu6ThXAQVWOzdu9dOnTrlljcvX7685c6dO8Y+tG/f3lq3bp3s/dK2tPy5CsErVarkghYtW56Y+fPn25IlS9wS6nfeeaddcsklyX4/AAAAnD+aoShYsKDlz5/fNfxB6FH6U7AzFZ6gAouMGTNaxYoV3SM+hQoVivf78VGht7pIqW6jVq1a9s4777g/x44dG2++l6ZuOnfubBMnTrTrrrvOli5d6lYXXLx4sRUpUiTZ7wsAAIDzSxeuqXXxitAVVGChZct1MZ8aPv74Y9u2bZstX77cTcfs3LnTqlSp4lrZtm3bNs5bfPPNNzZt2jRbsWKFC2DOnj1r1atXtxdffNE++eST1NglAAAAAMkUVGDh+fPPP23lypVxprgUGJQtWzZZuzJ8+HDr2LGjCypEwcJ9991nw4YNizew+Oijj6xTp07+WRGlTKlDVUItbk+ePOkenkOHDiVrvwAAAAAkLejAokePHjZw4EDLnDmzu7jPkyeP7dixw/LmzWsDBgxIdmCxYcMGq1atWow9rlq1qpuxiM9vv/1mjzzyiD311FP23XffuRqLZ5991m666aZ4n9+/f3976aWXkh4RAAAAACmWeGV0ElS8PWjQIPv+++/t119/dRf3Smdat26dFS9e3Jo0aZLUJpx///3XrdqtYCSQtrdr164421Dak1KlFFSoKOjNN990C/SpRkOzFvHp2bOnHTx40P/Yvn17svYNAAAAQNKCmrHQRb/WsVDxtIqp1UpMqUglSpSwBg0a2IQJE9zFf1LUyUmzHYcPH47xVKUrealRsen5Wj9j8ODB7kcKYvbt2+cCncaNG8d5fqZMmdwDAAAAQOoLasaiTJkyLhj4448/3MyBUpm0CrdodmDz5s3J2uP06dO73rlKhwq0cePGeDs8Kai4/PLLrXbt2jGer5oOdZUCAAAAcGEFFVio3exdd93lL67WuhMPPfSQdejQwT777LNkp0JJy5YtbeTIkf6jP3PmjI0ZM8Z936NZEY9a0SoFK9C8efNcsAMAAADgwgoqFUqGDh1qCxcudH9/+umn7cCBA7ZmzRp77733rHnz5sk+mm7dulmNGjXs/vvvdwXYKtres2ePPfbYY24bChpatWrl2sxqpkJtZbUonwIZpT7NnTvXpV55+wIAAADgwglqxkI0a6F6ClENw//+9z+bPn26u+BPiZIlS9qPP/7o0pzeeustV7itxe685cXVfrZ+/fr+r5UKpYLx48eP2xtvvGFHjhxxi+TF7iwFAAAA4PxL4wvMLzpHkydPdg+tZzFz5ky3yJ3azypdKVSpMFwL/KkWJEeOHKG6mxfFlP4VUvV9W/ZcnarbAwAAwIWRkmvmoGcs+vTpY/fcc4+ridAq2N6fjRo1opAaAAAAiBJBBRaKYF5//XX7+eefXaG10qKka9eu1rBhQ/v888+jZBgBAACA6BZUYKHF8LTORIUKcVNnKlasaPv374/u0QUAAACiRFCBhQqu1QVK61gE0srYWgH76quvjpJhBAAAAKJbUIGFCrRfeOEFl/Y0YMAAO336tA0bNsytxK1uTffee290jy4AAAAQJYIKLKR37972wQcfuMXqihQpYhMnTrQ6deq41rFqPwsAAAAg8gW9QJ7cdttt7gEAAAAgOqVKYKEai40bN7raikBaGVt1GAAAAAAiW9CBRbt27WzkyJF26aWXxkl90hoXBBYAAABA5AsqsNBMhYKKqVOnWosWLSJ/tAAAAADEK6ji7QwZMriZiptvvpnhBQAAAKJYUIFFqVKlXB3FiBEj4tRXAAAAAIgeQQUWMnz4cHvqqacsd+7cVrBgwRiPDz/8MHpGEgAAAIhiQdVYnDp1ymrVqmVNmza1O++809Knj7k5zWYACC1T+ldItR1q2XN1qm0LAACEt6ACizVr1rjgYtSoUZYxY8bwHgkAAAAA5yyoVKgKFSpYgQIFbOfOnee8AwAAAADCX1AzFkeOHHEdoZo3b+7SobJmzRpjRPS9OnXqhP8oAQAAAEhUUIHFyZMnbd26dXbllVe6P2OrUqUKww8AAABEgaACi/z589u0adOiYJgAAAAAJCaoGgsAAAAAEAILAAAAAEEjsAAAAAAQNAILAAAAAEFLcWCxdetWt35F165d7ZtvvnEtZwEAAABEtxQHFoULF7aBAweaz+ezZ5991vLkyWPXXXedvfLKK7Z06VI7c+ZMdI8oAAAAEIVS3G42ffr0bkE8PUSrbn/33Xc2e/Zse//99+3EiRPWoEEDa9y4sXuULFkyCocVAAAAiC4pnrGIrVChQnb//ffbqFGjbNeuXbZgwQKrV6+eW9+icuXKVrp06egaUQAAACAKpXjGIimqv/BqME6dOuXSowAAAABEtqBnLBKTMWNGN3sBAAAAILKd18ACAAAAQHQgsAAAAAAQNAILAAAAAEFLleLt9evX2x9//OFazQZSV6irrroq6J0EAAAAENqCDix69OjhFszLkSOH5c6dO8bRPv/88wQWAAAAQBQIKrDYvXu3vf766/bOO+/YY489ZhkyZIiCIQMAAAAQW1CBxdGjRy179uzWqVMnS5cuHaMLAAAARKmgirdLlCjhHps3b47S4QMAAAAgQc1YHDp0yJo3b26NGjWyli1bxqmxaNKkidWqVYuRBgAAACJcUIGFukAtW7bMypYta+vWrYszVBUqVIjw4QMAAAAgQQUW+fPnt2+//ZaRBAAAAKJcUIGFRzMXs2bNcutZ5MuXz+rWrWtXXnlllA8tAAAAED2CDixWrFhht956qx05csStWbF//37btGmTvfTSS9azZ8/oGUkAAAAgigXVFUratGljrVu3tp07d9rChQtt7dq1Nm/ePLdo3o8//hjFQwsAAABEj6ACiy1btti2bdtswIABljFjRv+o1alTx9q2bWszZ86MnpEEAAAAolhQqVBaFO/s2bPuEdt///2X4kXzfv75ZzfLUaRIEdeqNmfOnMk6NUrDGjdunJUvX95q1qyZrNcAAAAASD1BzVgULVrUypQpY4899pgdPHjQv1dfffWVjR492po1a5bsPX355ZetYcOGtmDBAlefUb16dTcbkhy9e/e29u3bu+ACAAAAwIUX1IyFjB071u644w7XelaBxoEDB+zkyZP2v//9zwUHybFmzRoXTCio0IJ6Z86ccUGJAobPPvss0U1oHY1PP/3UqlWrlpy3AgAAAHAeBB1YaMZCnaEWL15sGzZssLx587qAQoFGcg0fPtyuv/56/yrdSrHq3r27W9X7ww8/tGzZssW7KaVbPfLII9anTx/79ddfE307BTt6eLRqOAAAAIDUcU6pULFrKhQIaO2Kdu3a2c0335yioEK0/kXs2Q3NQJw6dSrRdKjBgwe713ft2jXJ0ejfv7+r2fAeml0BAAAAkDrOKbC45pprXKrTrl273N8TeiSVxuTRdjTTEUgBQIYMGVwb2/hs3rzZ1WUMHTrU0qdPeuJFa2qoDsR7bN++PXVGEAAAAMC5pUI9/PDDVqNGDcuePbsr3E5IpUqVkjXECiLU2SmQ0pZOnz6dYGeojh07uoLtqlWrJus9MmXK5B4AAAAAUl/St/rj8cQTT/i/++ijjwa9V2ovu3Hjxhjb8b7Wz2KbOnWqzZo1y6VLaSE++eOPP9wsxFtvvWVPPfVU0PsEAAAAIPnOKbAIdPjwYfv666/dCtwnTpywJ5980l3kd+jQwe67775k7UmLFi1cfYbSlLy1K0aNGuUChwIFCsTZRubMme2uu+5yxeKef/75x73/qlWrkvWekWR009StF8l+XfLWDwEAAAA8QQUWPp/P6tev71KeFFi8+eabNn36dGvVqpULLNQxKjltYG+55RYrXbq0NW3a1M2GqP2stjVp0iS3n0ePHrWyZctap06drEePHm7xPD0CqeWtZje8gm4AAAAAF05QgYW6OSllaenSpW6PFQio9avSozT7oJSl5AQWadOmtTlz5rjXKjAoXLiwe+2NN97otqvaCNV0lCpVKsGRqVy5cpwC8FB1ebv2qbprr6Tq1gAAAICUCyqw0IyFOjLpsWfPHlu5cqV/JkEpTEqTSi4Vgic026DtT5gwIdFN9erVK7lvBQAAACCVBRVYlCxZ0nLlymUPPPCA7d271ypWrGhXXHGFW6xOK3IPGDAglXcXAAAAQCg6p3UsPFoYb/LkyW62Qq1hvXUrvvzyS7eStgqsAQAAAES+oGYsRIXbav0aqF+/fq5uAgAAAEB0COrqf8WKFXbVVVfFGSmCCgAAACC6BBVYqMZiy5Yttnr16ugaNQAAAAAxBJUKlTFjRuvZs6d17tzZunfv7jo7BSpRokS8K2cDAAAAiCxBBRa7d++2gQMHuhHRAnWx9e/f35566qnIGjEAAAAAcQQVWBQrVsyOHz/OsAIAAABRLqgai0A7d+60efPmubazp06divJhBQAAAKJL0IHF+vXrrUaNGla4cGFr2LChHTp0yKVA3X333Xb27NnoGk0AAAAgSgUVWPh8PrcInoq0f/vtN7vkkkvcMLZu3drmz5/vFsoDAAAAEPmCCizUalYzFlpx++qrr/aPVoUKFVwx96JFiyJ/BAEAAAAEF1hoIbw0adLEuyCe6ixy587NEAMAAABRIKgZi+LFi7vOUL169bIzZ874h+uHH36wMWPG2I033hgFQwgAAAAgqMBCxo8fb1OmTLH8+fPbsWPHrFKlSta0aVN79dVXrVq1aowwAAAAEAWCWsdCypUrZ6tWrbIFCxbYhg0bLE+ePFazZk0rWLBgFAwfAAAAAAkqsNizZ4+9//771q5dO7v++uvdAwAAAED0CSoVKn369DZ9+nQrWbKkW8Ni1KhRLh0KAAAAQHQJKrBQ2tOvv/7qUqG0SJ6KuJUC9eijj9qSJUuiayQBAACAKBZUYOEpX768W21769atbgYjS5Ys1qBBA3v33XejeGgBAACA6BFUjUUgrVvxzTff2BdffGEzZsywihUrWpUqVaJnJAEAAIAoFlRg4fP5bOHChTZ69GjXdlY1F23btnUrbmsWAwAAAEB0CCqw2LZtm+sEpXUrPvnkE2vRooVlyJAhOkYOAAAAgF9QgYUWxVNwUahQIYYUAAAAiGJBBRYq0tbj+++/t99//91OnDgRYyivu+46q1q1ahQPLwAAABAdggos5L777nM1FlrLokCBAjFGrXDhwgQWAAAAQBQIKrDYvn27CyomTpxorVu3joLhAgAAABCfoNaxUFeonDlz2m233cboAgAAAFEsqMCiWLFiLtVJq28DAAAAiF5BpUIdOHDASpcubU2aNHFtZ3Pnzh1jJO+44w73MwAAAACRLajA4uzZs3b69GlSoQAAAIAoF1RgkSdPHhs6dGiUDyEAAACAoAILz759+2z8+PG2fv16y5cvnzVo0MBq1arF6AIAAABRIqjibZk3b56VKlXK3nvvPdu5c6fNnTvXLYzXpUuXKBlCAAAAAEHNWKjdrBbI69Wrlz3zzDP+0dywYYPVrVvXbr75ZmvatCmjDAAAAES4oAILBRAHDx60rl27xhgmrcJ9zz332Pfff09gAQAXyJT+FVL1nVr2XJ2q2wMARLagUqFy5MhhJ06csGPHjsUZpf3797ufAwAAAIh8QQUWl112mdWuXdvatGnjZi/k6NGjrt5i4sSJ1rJly8gfQQAAAADBpULJ559/bu3atXMF3JkyZbKTJ09agQIFbNSoUVa+fHmGGAAAAIgCQdVYSOHChW327Nm2ZcsWN2uRN29eK1u2rAsyAAAAAESHoAML+fvvv+3yyy93D9m1a5cVLFgwOkYQAAAAQHA1FjJhwgQrVqyYCy48qrm444477MyZMwwxAAAAEAWCKt5WR6iOHTvagAEDLH/+/P7hGjp0qC1atMi+/PLLKBhCAAAAAEEFFqqpUHDRuXPnGCNZunRpa926tQsuAAAAAES+oGos1P1JgcWRI0fs0ksvjTFae/futcqVK6doBLXY3pw5c2znzp1Wr149q1SpUqKvV6rVtGnTXIBz1VVXWfPmzS1t2qBiJQAAAADnIKircHWAaty4sd166622bNkyF2QoKHjllVds+vTpdttttyV7l7Zt22bVqlWzPn362MyZM936GIMGDUrw9cePH7eaNWta9+7dbfHixXbfffdZ/fr17dSpU8l+TwAAAACpI6gZC/nss8+sU6dOLijwlChRwtVXaG2L5OrRo4d7/tSpUy1dunQ2b948F7S0atXKbS+2Dz/80A4fPmyrV692rW21MF+FChVs5MiR9sgjjyT3bQEAAACkgqADi3z58rnOUAcOHLCNGze6WQx1iUpJStKhQ4fcSt1Kg1JQIQ0aNLAaNWrYiBEj7OWXX45zqApCBg4c6F8vI1u2bO572of4aOE+PTx6TwAAAACpI+jAwqMai6pVq57TXm3evNlOnz5t1atXj/F6fb1u3bp4t9miRYsY31e72x9++MGefPLJeJ/fv39/e+mll85p/wAAAAAkLqgai9SiBfWyZ88eZ7VuzYaoZiMpqqu4/fbbXeqUCrjj07NnT1cc7j22b9+e1GYBAAAAJFOqzVgEI0eOHHbs2DHz+XyWJk0a/6ZUQ5EzZ85EN63XPPzwwy4VS/UZCVEdhh4AAAAAUl9IzFgUKVLEtY7dunVrjCNUvUThwoUTPequXbvaggUL7Ntvv7VcuXIl+lwAAAAA50dIBBYq9q5YsaKNGjXKf5SagVCwcMsttyR45K+++qqNGzfOvvvuOxecAAAAALg4QiIVSvr27Wv33HOP6wpVsmRJe+ONN6x8+fJ20003uZ+/88477nurVq1yqVPDhg2z3r172+OPP27z5893DylYsKA1a9bs4owmAAAAEKVCJrDQYnrjx4+3jz/+2CZNmmR169Z1MxJe21rNamgl7syZM7tTpbSpRo0a2Zo1a9zDo+cQWAAAAAAXVsgEFqK0p4RSn7S6tx4eWscCAAAAoSMkaiwAAAAAhDcCCwAAAABBI7AAAAAAELSQqrEAIsHl7dqn6mFsGfFJqm4PAADgfCCwAELc6KZFU3UPs1+X+Gr2AAAA54JUKAAAAABBI7AAAAAAEDQCCwAAAABBI7AAAAAAEDQCCwAAAABBI7AAAAAAEDQCCwAAAABBI7AAAAAAEDQCCwAAAABBY+VtAABC3JT+FVJ1D1v2XJ2q2wMAYcYCAAAAQNAILAAAAAAEjcACAAAAQNAILAAAAAAEjcACAAAAQNAILAAAAAAEjcACAAAAQNAILAAAAAAEjcACAAAAQNAILAAAAAAEjcACAAAAQNDSB70FAIgSl7drn6pH+srumam6vezX5UzV7QEAkBLMWAAAAAAIGoEFAAAAgKARWAAAAAAIGjUWAAAAF9GU/hVS9d1b9lydqtsDkosZCwAAAABBI7AAAAAAEDQCCwAAAABBI7AAAAAAEDSKtwEASGWjmxZN1S2y+CGAcMCMBQAAAICgEVgAAAAACBqpUAAAmNnl7dqn2ji8kmpbQmpIzXMrW0Z8khq7BUQcZiwAAAAABI0ZCwAAgBSgOB/JFW2rqjNjAQAAACBoBBYAAAAAgkYqFAAAAGBGmluQQiqwOHLkiM2ZM8e2bt1qtWvXtqpVqyZ6eMePH3fP37hxo9WoUcNq1qwZ5HAAAAAAOBchE1js2bPHmjRpYidOnLAyZcpYr169rGvXrtavX794j+uff/6xm266yfbt22dXX3219enTxx566CEbPHjwuYwDAAAAwkxqtxKmVXRwQiaweOmllyxnzpz2888/W8aMGd2fderUsVatWtk111wT5ygHDBhgZ8+etdWrV1uWLFls1apVbtaidevWVrdu3eBGBQAAAECKhERgcfr0afvss8/siy++cEGFVKtWzRo3bmzDhg2zd999N85Bed9XUCGatWjRooV9+umn8QYWJ0+edA/PwYMH3V8PHTqUogFLDWdPnUqNzfgd/+9sqm4v3Ykzqbq9izHGF1M0nV/ObXD43Q0tqfm7y7kNLal5boXzGzo4t+ef93+9z+dL+s18IWDLli3aU9+ePXti7M2LL77oa9KkSZw93Ldvn3v+xo0bY/xs0KBBvtq1a8d7RH379nWv4cEY8BngM8BngM8AnwE+A3wG+AzwGUjZZ2D79u1JRg0hMWOxa9cuFwHlyZMnRiSUL18+27lzZ5zoyHt+3rx5k/V86dmzp3Xr1s3/fKVRqU5D75kmTZqkI7Awp2izaNGitn37dsuRI0eYHw1i4/xGLs5t5OLcRjbOb+SKtnOrmYrDhw9boUKFkjypIRFYeBe6x44ds+zZs/t3WgehuovYvOcfPXo0xkVyQs+XTJkyuUegXLlyJTlAkUZjR2ARuTi/kYtzG7k4t5GN8xu5ouncxnc9Hp+QWCCvSJEibt/UNjaQvi5cuHCc/S5YsKClTZs22c8HAAAAcH6FRGChiK9+/fo2atQo/9FqNuKrr76yW265Jc4IZMiQwZo2bRrj+SrMnjRpUrzPBwAAAHB+hUQqlPTu3dt1dcqWLZtVrFjR3nvvPcufP7/dcccd7ucjRoyw559/3n799VcrUKCA+3ujRo1cOlP16tVtyJAhrqNU27Ztz++IhSmlgfXt2zdOOhgiA+c3cnFuIxfnNrJxfiMX5zZhaVTeHSqnfvbs2S6g2LZtm1tF+9VXX7XcuXO73dPPtPjduHHj/HUYP/zwg7311lu2adMm155Wz1cwAgAAAODCCqnAAgAAAEB4CokaCwAAAADhjcACAAAAQNAILAAAAAAELWS6QiH1rV+/3hYtWmQbNmxwf9cq5sWLF7fSpUtbpUqVrFmzZpY+PR+BcKR2zHPnzrU///zTf36zZs3qzq0ezZs3968Pg/CzfPlyW7ZsmTu3euzfv99KlCjhzq0aVTRs2DD8DgrOP//843539W+yzq3Wa1KTEp3bq666ym699Va79NJLGa0wdPbsWVu8eLGtXLnS/7t75MgR/7/L9erVc7+/CE9//fWXzZ8/339ut27d6lai1u9u+fLl3XIHWbJkCc+DS0UUb0cY/cP2xRdf2KBBg2zFihV2+eWXW6lSpaxkyZJuYcHt27e7/9D0s0suucS6dOliHTp04D+yMLF27VrX/UxrtigoLFu2rDu3uug8fvy4O7d6ji5W1Kr5qaeecu2YEfp0/tQ2+5133rHNmze7/6y8c5snTx7X/c773dXvtc6t2mvzH1l4WLJkib322mv27bffuhVsy5Qp487vFVdcYQcOHHDn9vfff7e9e/faAw88YE888YR7DkKfzp86VOr3d9++ff5zq3+X1ULfu7mnGwYKLPS727p1a0uXLl3oHxzs66+/toEDB9qPP/7orqN0A0DntlixYrZr1y73u/vbb7/ZmTNn7NFHH7XOnTu7gCNaEVhEGF1QarFBXXDcf//9LoqOz+nTp23mzJn2wQcf2PXXX2/PPvtshI1EZNIdkTRp0rhze/PNNye4Lsm6dets+PDh9s0337i7Zwh9+s+rW7du7tzq91ezi/HxFg999913rWfPntayZcvQPzi4WeKrr77anV+twZQ2bfyZyFqr6ZNPPnF3Q/WZQOh7++23bfTo0e7ctmnTxvLlyxfvTmvmcfz48e7mwcSJExP8/xmhQ8GCljHQedX5rVGjRrw7p5u6WgLho48+cudf/z5HKwKLCPPff/+5/7AS+k8roSBDq5kj9J06dcotBJlcCiA5t+EhpedWOL/hI6Xnl3MbPlJ6brV8mC5YSUUODyk9v9H+u0tgEcEUQSv1ybvzqdxeTcfrQlNRd44cOSL46KOD/nOKbzp9x44dduLECTddi8igmSelQ2nxUE3HI7JoBlnpE5rVQHhTapTq35STrxWaCxcubOXKlbPMmTOH94HBXVf98ccf7twePHjQ/c7q/9loTn2Kja5QEUqrl2vqvUePHu4IJ0yY4HJ577rrLpdOo9xerWaO8PTZZ5+54mz9p9W0aVOX4xloxIgR9uKLL4bnwUU53Qy48cYbXUGvCnmVs921a1erXLmy3Xvvva5matasWVE+SuFJM8onT56M9/H++++74EI/1w0DhJ9///3X2rVr5wL/OnXquFx71ctce+21VrRoUXvppZfcbAXCk/5f1b+/Cv7vvPNOVytzww03uP+LlZoc+//haEVgEaEef/xxlxeoIu7du3fbgw8+aC+88IKLsA8dOmQPPfSQu0g5fPhwhI5A5FIOtgru+/bt62oolM9ZpUoVVxSK8KffS81CjRw50hUH6kbAuHHj3F0y1Vfod7tjx45coIQhnTfdtY7vMW3aNOvevbv72SuvvBKGR4cWLVq4+jb9u6z0mT179rhMAf0/qzoM1TT27t2bgQpDY8eOdYGiggk1WNB1lAq3dSPg559/dpkgTZo0cbNVUc+HiHTJJZf45s+f745t0qRJvrJly8Y4zrNnz/oKFSrkmzt3bkQefyR78cUXfc2aNYtxiEOGDPFlyZLFN2LECPf9V155xde2bdtIHoaItG/fPt3O9G3YsMEd35kzZ3yFCxf2vfTSS/7j/euvv9xztm/fHpFjEMnefPNNX5o0aXwPPfSQ75dffvGtWLHC/2jQoIGvW7duvmXLlvl27twZycMQkdauXevLmjWr78iRIwke38yZM30FChSIyOOPdI0aNfK98847CR6m/q0uU6aMu96KdixiEKF0F1sFRKL8zthT6+ospLuimppHeNF51Z3sQI888oibaleLWd1JQXhSXYwaL3j1T/q7WpNqjRKP2leq6FN3yhBelNJWtWpV15JSrWWHDh3qr6nQedbvtdJmEH6Usqj1R/T7mRD9G610Kf2/S+F2eNH51flLiP6tVp2FnhftSIWKUEqn0BoVKiDTf2T6wL/++usufULFRy+//LK7MGGNg/CjmorPP//cpT7pQtSj7yv3XrUVygVF+NFNgCuvvNKtVaIURp1HBYqjRo1yC23Jl19+admzZ3drWSD8aJE0rSOkBUpr165tvXr1ivF7jPCkgFApT2+88Ua8aYpap0bnWu3gCSrCjxYl1To0f//9d7w7P3nyZLcgsc5vtKMrVIRSAKF83Q8//NAt5qI7KVrtNW/evC6w0OJ4yuHWGhYIP//73/9crq76oevuZyDdCVXxb4MGDdwFKcKL/nNS0bbyeHVDQE0WnnzySVu9erVrurB06VJ33h977LHwOjDEocUs9furXG1dbOrvOtcIT/PmzfOvnK5/fy+77DI3w6yGDCrML1CggHsOHYTCjwJD1dBokbzGjRu7G0CandIMhVZbV23Ne++957IHoh2BRYRT4Zj+QdM/bIq01a1CFye68KT1XXjTbIVmnZRCEZvOte6e0W42POliRJ3d9J+X0hZVCDpmzBi3IrfunHFDIHLoJtCnn37qmmvo0alTp8g5uCikmUX9rmolZrUk1bpDmom87rrrXCMGVtsObwsWLLDp06f7283qmqp06dJ23333ucARRmABAAAAIHgUbwNhZuHChfbdd98ludcVK1a0Vq1aJfk8hA4VdepuZ3Lcc889rKoehjRDoXRFFfFqscPY+fZqVapZSK2DgPCjdqNaoLR8+fJutlEpjR9//LFrFX377bdTnB/G1ARn1apV7gi0rlBsqn1UrY3Sz6MZgQUQZpTjqQJBLYCYGHWpILAIL0pt0wJbys1WX/TE6Nwm9RyEFp3fli1buroZXXQqRWb48OEuvc2jxUxVmE9gEX5Gjx7tcuyVpqoOX8q9b9u2rUtN1bpSb775pk2dOtWlIiO86BzedNNNtnLlShdgqCmOFqotW7as/0Deffddt84FgQWAsKLuT+rqpc5BdKCILFmyZHGLHeq8Dh48OLIODm517Y0bN9rOnTvdTMWQIUNcdygV42vRS4Qv1TOqoUKfPn3cArSqm1EBt+qkVG+hGz09e/Z0DwKL8KPzqm58mpHav3+/DRw40M04Tpo0yRo1ahR+B3Qe0W4WCDOaqVBXqLvvvpue2RFGFx/jx4+3YcOGuTubiCxKY1SRp2ak8uTJ4y4ydZ7VwU83ChC+dCdb3Rafe+45d367devmAkilLOr3Wu68804XZHhrTCF8qGi7Y8eO/lbf6rip4EIF+ZplxP+PVKgopR75aj+rtmlaTA/hRf/AlStXzo4dOxZeO44kqZOXOrnpogSRJVeuXHHWrGjSpIlrQao0C9VdqAYD4Ud3s7XmjNLdNPOorouamShevLj/YDSrkSlTJtaxCEP63Y0dEGqGSp2gtG6YujDi/yGwiFKaztMvhVrjXXPNNfb444+73G6ED7UvRGTSFDsizwMPPOACCN3ZVi6+d1NH6W+qnVKQsWfPHnfXG+FFM8lan0KzT1qkVGtGTZw40X8Q+j/3rbfeculRqq9BeNFMo86trpl0reS169eaQ19//bX7U23BYaRCRSsVHOnuivpsK60m8K4KACD1qXZGC5Nqld41a9bEeINSpUq5VCn1xGdl5vCjRgrKt//+++/djGNsSo1SGpQKuBF+tHjlM8884+qhNPMUSOsKadZR3dwIGlnHIiqoQ4UCCKXNFClSxHLnzh0Vxw0AwIV29uxZf12FZ9OmTa4LmFKhgEhG8XYEmzJlimtZqHxP3Q3TVK0KBjVd269fP7dEPQAASD1esXYgdYciqEA0SOOjUiwiaaEWtbzT9Kv63auntv5RU/6upmqV/qTiXwUfTN0BAADEr2vXrrZ69eokh6dXr14uNSqaEVhEKK26rF+EBx98MN4j1GIvqqv49ddfYyzwgshw8803u0BSxaDqTKIFtxAZVCg4YsQId151fnXTAABw/qhZilZPVyfNxKgNvK6/ohmBRYRS67slS5ZY+fLlEzzCMmXKuBWcmzdvHqGjEL2+/PJL12Xmzz//tLVr19rPP/9MbU2EWLVqlesuo/Oq89ujRw9XVAgAOD++/fZbu+uuu9z/pdG+snZSCCwi1A033ODqKrSIS3ypTupaoYBix44dbjEfAOFJvdXVkQYAcP6oDfSMGTNs6dKl/naziIvAIkLpTma9evVcn/T4aixmz57t7np27tw5QkcgsqlWRm0Llc+ZI0cOW7Rokb3//vu2bds2Nw2rFX3VAQyR6eDBg661ISJPo0aN3NpCSnPTv+Fev3yEv9GjR7sZR6UxqrEKxdzhRSXJWqsiY8aM1KYmgq5QEUpTdevWrXMF3MuXL3cXna+++qp98803rjuFaisIKsKT1h+pUaOGK7xPly6dy7lXf3y1ONTFyMqVK10HMM1GIfxoBdcuXbq45gq6AJk/f36Mg1BQoVVg9Scij2abt2zZ4hpvKF0VkUO/t/q/VwusqUOj0msQPpT9oWCQhjeJY8YCCDPvvvuuW6Tnjz/+cCkw1apVsxYtWlifPn3ckeiuStOmTd2F6eDBg8Ps6NC6dWvbvn27dejQwU256y6nZqZ69+7tBscLLP79919mLSKcfpe5iIlMKgTWzSDVQwKRhBmLCKOLjY8//jjZR6XVIn/66adkPx8Xny46GzRo4M+r1+KHN910k3/HdCHSsmXLZLXGQ2jRbJQK74cOHWrt27d3f86ZM8cFk1r5VRciCH/qyqcLS49mlV988UX3b3fg9wkqIoNuAik/P1C2bNkIKhCRCCwizP79++2VV15J9lFNnTrVBRcIH9dee6398MMP/guQunXrug5QgdQRTO2EEV40G5E+fXqrUKGCf8dr167tzufcuXPtnnvucTm+CE///feftWnTxjXMyJs3rwsmlMqo1EalxfTt29f9PnOOw48CQqUtxvdQwe+YMWPcz1QPB0QyUqEizMaNG6106dLJbi165MgR959b7LspCF266Ljlllts7969Lg+7YMGC9vDDD9t9991nVapUcRcoSp9RLi9t8cKPAsLHH388TgtZFewrxU3F+gosSYUKP+rS179/fzcrpaJsBYr//POPS2287bbbXH2NuvnpOQmtQYTQpNq2ypUruxsDegQ6c+aM696m76u2Yvfu3aF5EDhnummg80rTFCOwiDQnTpxwedkpoWLuokWLpuQluMiOHz9uw4cPtyFDhtiaNWv8dzgvvfRS10lGFyaqsUD40czE7bffbiVLloyTpqgZDQWVBBbh6Y477nAF2S+//LI7AP3+vvTSSzEaLSgFLkuWLC79DeHjwIEDrj5q8+bN9tFHH7nGCx6lMz722GO2fv368DkgpIhX+6a1w9RE5dZbb3VNVaJRzLAaYU93wbRCJCKbLjw6derkHirw1OyFCrkVWCC8NWzY0M08qmV0bGoxqzVoJkyYYFmzZg3vA41C6igT2M1Lfz958mSMIm3V2ah7EMKL/u3VTQHd8Ln33nvd7KKaZ6jlOyKfivB1U9dbuHTWrFlRG1iQCgWEEd3x0n9Uybnw0AWLLk61rgWAi08BoYrwtYaQggx1+7rkkktc3YXSUXW3W2tYjBw50t3xRHjSjZ6uXbu62oo33njDChcubB07dmTGAlGBwAIII+r+pCl1dX1Kiu56X3/99a6LFMKD8nQnTZqUrJ1VTj4rboeffv36uTQn5dq/8MILLjVKi+Lp3KsAWOdVAQgdocKf7loroNCsstYcIhUqvGmGccSIETZt2jTXjVFfFypUyNVFKYVRv8egxgIIKwos9I+Z7oAlRYX5q1evJrAII7qw1B1szUjFLgCNTYGjCrkR/lTArY5BBQoUsJo1axJURJBjx465To1asFQzUQhPO3fudB36lG6uGjjVpqplsGanlAI1efJk13CB2igCCyCsqGXhhg0bkr3PuvDUlDzCg7rHqLOX8rMHDhwYHjsNABFOq6Wr/mn8+PHx3vTRTbxatWrZ7Nmz3c2BaEYqFACEENXFaDX1cePGxVj4EOFP51RrkiRFgWVgVyGEj++++86WLVvmugM1b97c0qb9/5cL0yyjuoB99tln4XNAcLS20Ouvv27NmjVLcESUUaBz3qVLl6geNRbIA4AQorVHpk+f7u6OIbLoolPBhWYdE3soNQrh59lnn7VWrVq5HHy1FlaNm1JoPFp7RusMIfyo/XfshWgDaQ2a3377zdVbRDtmLAAgBKiLl/Kxk9syWH3z1XJW3YUQHtSOVOuQaNYicHV1hD8V82pxS+XbV61a1QURyrn/5ZdfXDChGQz9XbNRystHeFm4cKE1btzYdXBTO2GvxmLfvn3u9/nNN990/3br5kG0/5vMjAUAhAD9h3TnnXcme0+0arPaWSJ8aI2Sp59+2p1nFeojcmj9gmLFirmgQtSAQSus60JUC6X99NNPkXOwUahOnTr2/fffu+YpSndSYHHZZZe5gFHd3dTNTTPN0R5UCAvkAUCIUH1FcovtteI6wk+fPn3crMWpU6dcVxlEBgUVmonQHey8efO6g1LLYOXl58+f37Ui1QUowlf16tVdO3C1D9Z5VpBRsGBBfo9jIRUKAEKA8nNT2qqwc+fOVrly5RDYeyRFM1I1atRwK/Qm5e+//3ZrHuguKcKH2pAuX77cevbs6dY1CKT1D/Q9pcuQChU5tP6MWgkXKVLErVUC2s0CAHDeqajzq6++cqkTSZk5c6YNGDDA5s2bl9RTEUJOnz5t77//vmsbrZS32KZOnWrDhg1znwOEn08++cS++OILK126tHXr1s2taVG3bl3btm2bCxgVPGo2MtoxYwEAwHmmwELF9nokRYW/WiyPwAIIDQoIH3/8cVczo8Bx5cqVVrZsWTdToZkoLZCnhRDV7U21F9GMwAIAgPNswoQJKWohrNztxHrmI3So3kn59sldGE1F3SrkV4E3wkO9evXc76PS3KRFixYufXXLli2ulkZ1FwoytLBp27Ztw+OgzhOKtwEAOM+0rgEi06xZs1xXqOQGFs8995xLhyKwCB8q1i5Xrpx/h7WA5YkTJ1xQIfqzTJky7nnRjsACAAAgCGPHjnXtSJNj69atyXkaQkiVKlXcrGPLli3dXqm9rOorPFrQdPHixW6RxGhHKhQAAMA5Ul69cu5TQout5ciRIyUvwUWkGSkFEjfccIONGzcuxp4oUNQMRpYsWdwiiGnTRvcSccxYAAAAnKOSJUu6ByKX0pxWrFgR70KH6gammovnn38+6oMKYcYCAAAAQNCie74GAAAAiMfZs2etY8eO7s/kmDt3rn3++efJeWrEYsYCAAAAiEVrVqRPn951e9KfSfnwww9dzYX+jFZJjxIAAAAQpfLmzZvsI3/ssceS/dxIRGABAAAAxJIuXTr766+/UjQu2bJlS9HzIw2pUAAAAACCRvE2AAAAgKARWAAAAAAIGoEFAAAAgKARWAAAAAAIGoEFAAAAgKARWAAAAAAIGoEFAAAAgKARWAAAotqhQ4fsv//+C8kxOHjwoJ09ezYk9w0AYiOwAIAQMWXKFFuwYIHbm++//97uueeekL3gTYljx47ZG2+84S7gL6QNGzbYhx9+mOhbrl692vLmzWt79+69YLu2adMm++CDD5L1fnXr1rX3338/Wc8FgIuNwAIAQsSIESPsm2++cXtz4sQJO3z4cLLuVn/55Zf2xBNPhMhRxKXjeOaZZ+zAgQMXdB9///1369evX6LvOWjQIGvUqJEVLFjwgu3bmjVrrG/fvsl6v/vvv9/to8/nS9bzAeBiIrAAgBDUtGlTmzZtmmXMmDHJvdu+fbv9/PPPST4PMZ08edImT55sDzzwQKoPjQKB1Ehhuvfee23Hjh22ePHiVN9HAEhtBBYAcBGcOnXKunXrZgUKFLBLL700zp11BRUVK1b079kvv/xiderUscyZM1uhQoXs7bffdj/r0aOH9erVy5YtW2b58uVzF8qiVCrdic+ZM6dL9dHzvLveS5cutRIlStjcuXOtUqVKli1bNmvVqlWMGQXNMjz66KPuTn727Nnttttus507d/r3R7MP9913n9u23rdTp07uNYlZsWKFO6asWbNavXr1bNWqVf6nnz592l544QW7+uqrXTBVuXJlf1qYdOjQwV577TV78skn3ftpDIYMGRLj7aZPn+5enyVLFmvYsGGS6U0rV660I0eO2C233OK288orr9gNN9wQY5sK2vR+v/32m/v+sGHDrHbt2u49ihcv7r72fPzxx9a6dWvr3bu35c6d24YOHep+9O2337rj1muuv/5627NnT4z30H7Ur1/f/Vyfh//973/+n+s4q1WrZvPnz78In1IASCEfAOCC69q1q++KK67wzZkzx7d27VrfAw884MuZM6evR48ebl8mTpzoK1KkiH+/9PdOnTr5duzY4Zs9e7YvT548vlWrVvn279/v69Onj++aa67xbdmyxXf06FH3vVy5cvm6d+/uW7dunW/mzJm+7Nmz+8aMGeO29+OPP/oyZcrka9eunXvvSZMm+QoUKODr1q2b//2aNWvmq1mzpu+nn37yLVmyxNekSRNf5cqVfWfPnnXPqVu3rq9169bu9StWrPBVrVrVd//998c7jrt371Yej69cuXK+WbNmudc8+uijvrx58/oOHjzoXvP222/7LrvsMt+UKVN8mzZt8j3xxBPuGM+cOeN+3qZNGzcGw4cPd6/X+KVLl863fv169/Ply5f70qdP7xs0aJBvw4YNvpEjR7pj0iMhGo/ixYv7f/z777+7/fzzzz/93+vfv7+vfPny7uvFixf7MmbM6Pvwww99mzdv9g0ZMsSXJk0atz/y1ltv+bJmzer2Xfuzd+9e32+//eb2a+DAgW6/Ro0a5StYsKA7dk+JEiV87du39/3111++uXPnup8tW7bM/3ON68MPP3zBP6MAkFIEFgBwgZ06dcpd2Osi2nPy5Elfvnz54g0sTpw44UubNq27qPacPn3a/3dd0CoICHT48OEYXzdo0MDXpUsX9z0FFrqA3rVrl/85CkJq1arlvlbAogtmXeB79u3b55s3b5670NfrFQQpiPF888037jVHjhyJM5peYKGLao+2o8Dqo48+8n8vcJ9Xr17t9tG7aFdgoUAmcLwUPHnB0kMPPeRr1apVjPfWBX5igcVLL73ku+GGG2K8RuOoQM2joOLNN9+Mdx91zhT8jBgxwv1c50HHFKhDhw6+Fi1axPieAjgvsNA4eMFKfOdW+vXr56tfv36ccQWAUEMqFABcYFu2bHH5/TVq1PC/s9J/AlOfAmXKlMkV+z700EMuLea5556ztWvXJrrXeg+lDt14440u7enHH390aT+edOnSubQbz1VXXeVPZVLRs1JwrrjiCv/P8+TJ49J40qZN6zopnTlzxpo3b24NGjRwD72XUq0S2y/tu0fbqV69uitkFrV7nTFjhnXs2NGNi95LAve5SJEiMcbr8ssv9++z3jdwPCXw/eLz119/WeHChWP86MEHH7TRo0e77y1fvtzWrVtnbdu29T9HKWdKK7vuuutcKpTawQbu42WXXRZje0ntl8ZBBeadO3e2KlWqWPfu3d34B9JxKyULAEIdgQUAXGCqkxAFF4Fifx2oT58+tnXrVlf38Mcff7iLUNVIxOfXX3+1a6+91jZu3Oi6Rc2ePdt1FwqUJk2aBI9agYxqHhKiIEg1BLrA9h6qK1BgcOWVVyb4OtWVBFLnK28sunTp4uonSpcubW+99ZYLXmJLbJ+1nZSMpyh4CqwbkTZt2tiuXbtsyZIlNmbMGLv55ptdjYWoa1eTJk3c8SvQ0ziXL18+weOV5OyXAsVt27a5MVCL3KpVq7q6DI/2MTCoAoBQRWABABeYLhJV9Bx48bh//35XoB2ff//91z7//HNXSN2+fXt/Ybcu5D2BC7yp0Ffv8emnn7pZBV3sJzXDEeiaa65x+xPYaUozCyoSV8ChoGX37t1Wrlw5NyOihzfbkiNHjgRHM3B/dZf/hx9+cEXaouDn6aeftq5du1qtWrXiXPAnRdsJHE+ZNWtWoi9TEKMZiUDafxVgf/XVV25dEc0SebSPt956q5thUHF4rly5XLCXmKT2SwXvOrcKXvReKr7XzEvgsaxfv97tKwCEuvShvoMAEGl05113vHUhffToUXfn/L333kvwrrTalipVRmtc3H777S4AUYel/v37u6FRGpPu8GsxON1RV4rR5s2b7eWXX3YXtrpI1uxFqVKlkjWUSvHR3XN1gtJ+Ku1Jswi6O58hQwYXROjOvt5Ld9uVzqOf6+68vpeQCRMmuHQpbV/HqzQkHY/oYlodlpSepSBJQVFKqMPWJ5984tqzar8VtGhmJzHqIKUUIwVJgWlhusC/8847LX369K7tr0f7qPH46KOPLH/+/K7rk449MU899ZTrXqXxuuOOO1xKmtdhyqNZJQUxeo5+ptkSzQR5NDMSGOAAQKhixgIALgJduOsCdc6cOe4iWheS+p53518zDl4bVKUd6UJZAcnAgQNdDr5aqzZu3NjtuS5+X331Vfc9BRhqS/vFF1+4lrNa4VlfqwbCq0HQbIkunANpViPwIlqBgtq/KijRttRaVqlBnuHDh7v9nTRpkrug12t1lz2+dCUFHLpoHjt2rJs5UVChWQ8dky7eRfupWYB3333XXXxrbPQateIVzWJoJiWQVz8iClI0w6LAZ/Dgwa51q/axZcuWCZ5dBRaqw9DzAqm+Q9vR+hbe/omOV4HUqFGjbNy4cS6YUZtdL2DTrIJa/AbSLJP2S/uj/dJ2R44c6dr7ilr5arV1jdHrr7/uAkbNWigFSzRemrEIrPMAgFCVRtXkobpzAACcTwqalOKlGR0vKFK9Q8mSJd06G5oNuphUd6JZLQVvABDqSIUCAEQtpWJpgTqlQ2l2wSt2V+rZxQ4qRLNV2hcACAfMWAAAYOYK0726FKUiqfsTACD5CCwAAAAABI3ibQAAAABBI7AAAAAAEDQCCwAAAABBI7AAAAAAEDQCCwAAAABBI7AAAAAAEDQCCwAAAABBI7AAAAAAEDQCCwAAAABB+/8AIkIDcU6LPL0AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sig = lambda z: 1 / (1 + np.exp(-z))\n",
    "op = op.assign(fitted=p, teaching=sig(-0.6 + 2.2 * op.angle - 0.1 * op.dist))\n",
    "bands = pd.cut(op.dist, [0, 6, 12, 18, 24, 30, 120], right=False)\n",
    "tab = op.groupby(bands, observed=True).agg(shots=('goal', 'size'), goals=('goal', 'sum'), actual=('goal', 'mean'), fitted=('fitted', 'mean'), teaching=('teaching', 'mean')).round(3)\n",
    "print(tab)\n",
    "ax = tab[['actual', 'fitted', 'teaching']].plot.bar(figsize=(8, 4), color=['#1f5f63', '#a3471c', '#b0862b'])\n",
    "ax.set_ylabel('conversion / mean xG'); ax.set_xlabel('distance band (yards)'); ax.set_title('2018 World Cup open-play shots')\n",
    "plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ef90be8",
   "metadata": {},
   "source": [
    "## Where geometry stops working: the penalty spot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "66e3eae5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:48:15.629607Z",
     "iopub.status.busy": "2026-11-03T14:48:15.629359Z",
     "iopub.status.idle": "2026-11-03T14:48:15.640219Z",
     "shell.execute_reply": "2026-11-03T14:48:15.638593Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "in-game penalties: 29 | scored: 22 | actual rate: 0.76 | geometry model predicts: 0.12\n"
     ]
    }
   ],
   "source": [
    "pens = in_game[in_game.shot_type == 'Penalty']\n",
    "pp = fit.predict(sm.add_constant(pens[['angle', 'dist']], has_constant='add'))\n",
    "print('in-game penalties:', len(pens), '| scored:', int(pens.goal.sum()), '| actual rate:', round(pens.goal.mean(), 2), '| geometry model predicts:', round(pp.mean(), 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ee03ddb",
   "metadata": {},
   "source": [
    "## Limits\n",
    "In-sample scores flatter every model; the course's Session 6 refits on the group stage and tests on the knockout rounds. A geometry-only model does not know about the goalkeeper, defenders, body part or shot type, which is why the penalty gap appears. One tournament is a small sample: 1,556 open-play shots and 129 goals."
   ]
  }
 ],
 "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.13.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}