diff --git a/DESU_regression.ipynb b/DESU_regression.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..f811d72947094736331d646980e18bb598bee5f5
--- /dev/null
+++ b/DESU_regression.ipynb
@@ -0,0 +1,903 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 65,
+ "id": "ef2a766b-a3a1-431f-8dac-6c57c1157312",
+ "metadata": {
+ "tags": []
+ },
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import scipy.stats as stt\n",
+ "import statsmodels.api as sm\n",
+ "import statsmodels.formula.api as smf\n",
+ "from patsy import dmatrices\n",
+ "import pandas as pd\n",
+ "\n",
+ "import matplotlib.pyplot as plt\n",
+ "%matplotlib inline\n",
+ "import seaborn as sb\n",
+ "from statsmodels.graphics.factorplots import interaction_plot\n",
+ "\n",
+ "sb.set_style('whitegrid')\n",
+ "sb.set(font_scale=1.5)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "102b1593-da1a-485c-a8a3-0129619a69ac",
+ "metadata": {},
+ "source": [
+ "## Synthetic example\n",
+ "\n",
+ "We start with an example where we control the dependency from 2 predictor variables ($x1$ and $x2$) to a response variable ($y$).\n",
+ "\n",
+ "For further reference on the models and , see the docs about [ordinary least square](https://www.statsmodels.org/dev/examples/notebooks/generated/ols.html) estimation, as well as [patsy](https://patsy.readthedocs.io/en/latest/).\n",
+ "\n",
+ "### Exercise\n",
+ "- Generate $y$ such that it has a linear dependency with each of the variable, test a linear model on it.\n",
+ "- Add a non-linear dependence to $y$ with the product "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 57,
+ "id": "3576be4a-1c48-4f7d-abbc-f18064c930e6",
+ "metadata": {
+ "tags": []
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "<div>\n",
+ "<style scoped>\n",
+ " .dataframe tbody tr th:only-of-type {\n",
+ " vertical-align: middle;\n",
+ " }\n",
+ "\n",
+ " .dataframe tbody tr th {\n",
+ " vertical-align: top;\n",
+ " }\n",
+ "\n",
+ " .dataframe thead th {\n",
+ " text-align: right;\n",
+ " }\n",
+ "</style>\n",
+ "<table border=\"1\" class=\"dataframe\">\n",
+ " <thead>\n",
+ " <tr style=\"text-align: right;\">\n",
+ " <th></th>\n",
+ " <th>x1</th>\n",
+ " <th>x2</th>\n",
+ " <th>y</th>\n",
+ " </tr>\n",
+ " </thead>\n",
+ " <tbody>\n",
+ " <tr>\n",
+ " <th>0</th>\n",
+ " <td>1.951076</td>\n",
+ " <td>1.654086</td>\n",
+ " <td>0.061812</td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <th>1</th>\n",
+ " <td>-0.024425</td>\n",
+ " <td>1.342077</td>\n",
+ " <td>0.091460</td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <th>2</th>\n",
+ " <td>0.638307</td>\n",
+ " <td>-0.854780</td>\n",
+ " <td>1.663423</td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <th>3</th>\n",
+ " <td>0.958588</td>\n",
+ " <td>1.135066</td>\n",
+ " <td>0.393107</td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <th>4</th>\n",
+ " <td>-0.944212</td>\n",
+ " <td>-0.223619</td>\n",
+ " <td>-0.959402</td>\n",
+ " </tr>\n",
+ " </tbody>\n",
+ "</table>\n",
+ "</div>"
+ ],
+ "text/plain": [
+ " x1 x2 y\n",
+ "0 1.951076 1.654086 0.061812\n",
+ "1 -0.024425 1.342077 0.091460\n",
+ "2 0.638307 -0.854780 1.663423\n",
+ "3 0.958588 1.135066 0.393107\n",
+ "4 -0.944212 -0.223619 -0.959402"
+ ]
+ },
+ "execution_count": 57,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "# number of samples\n",
+ "n = 50\n",
+ "\n",
+ "# exogeneous variables\n",
+ "x1 = stt.norm.rvs(size=n)\n",
+ "x2 = stt.norm.rvs(size=n)\n",
+ "\n",
+ "# error\n",
+ "e = stt.norm.rvs(size=n)\n",
+ "\n",
+ "# endogenous variable scaled by factor a\n",
+ "a1 = -0.5\n",
+ "a2 = 0.7\n",
+ "y = a1 * x1 + a2 * x2 + e\n",
+ "\n",
+ "# build dataframe\n",
+ "df = pd.DataFrame(np.column_stack((x1,x2,y)), columns=['x1','x2','y'])\n",
+ "df.head()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "6345745a-3046-48fe-b3b8-54e09c507dc9",
+ "metadata": {
+ "tags": []
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAHMCAYAAAAarpbgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABGgklEQVR4nO3de3RU1f3//9dMboQwIQlXiRFQQ2IIfCopVy9YpYK2VVCLgFX5VMUC0t8SRWsV+6k31FqpUoVaq3IRxQvIEixWLloVQRovEJAASiAYIGISCAm5zZzfH3yTGpJMZiaTOefMPB9rdVVmdmbezk7klX32eW+HYRiGAAAA4DOn2QUAAADYDQEKAADATwQoAAAAPxGgAAAA/ESAAgAA8BMBCgAAwE8EKAAAAD8RoAAAAPxEgAIAAPATAQoAAMBP0WYXYJb169frww8/1Pbt23Xo0CGVlpYqOjpaqampGj58uCZPnqzU1FSzywQAABbkiNSz8K6//np9+umniomJUbdu3ZSSkqLS0lIdPHhQHo9H8fHx+utf/6rzzz/f7FIBAIDFRGyAeuutt9SjRw/l5OQoNja24fH9+/fr97//vbZs2aLk5GStX79eHTt2NLFSAABgNREboLw5cuSIzjvvPEnSc889p5EjR5pcEQAAsJKI3QPlTdeuXZWUlKSysjJVVVUF/DqGYcjjIZ/6yul08HmZjDkwH3NgPubAfGbOgdPpkMPhaHUcAaoZX3/9tcrKyuR0OpWVlRXw63g8hkpKKoJYWfiKjnYqOTlBx45Vqq7OY3Y5EYk5MB9zYD7mwHxmz0FKSoKioghQPjMMQyUlJcrNzdUTTzwhSfr1r3+ttLQ0kysDAABWE/EBauXKlbrrrrsaPXbmmWfqiSee0C9+8Ys2v350NK22fBEV5Wz0/wg95sB8zIH5mAPz2WUOIj5AdenSRYMGDZJhGDp06JAOHz6sgoICvf322xo8eLB69uwZ8Gs7nQ4lJycEsdrwl5gYb3YJEY85MB9zYD7mwHxWnwPuwjtFYWGhHn30Ua1du1Y9evTQ6tWr5XK5Anott9ujY8dOBLnC8BQV5VRiYryOHTsht5t9B2ZgDszHHJiPOTCf2XOQmBjv0+pXxK9AnSotLU1PP/20rrzySu3evVtLlizR1KlTA349NiH6x+328JmZjDkwH3NgPubAfFafA2tfYDRJVFSULrjgAklSXl6eydUAAACrIUC1oK6uTpLk8Vg3/QIAAHMQoJpRU1Oj999/X5La1AcKAACEp4gMUNu2bdNf/vIXFRQUNHlu7969mjp1qvbv36+OHTtq/PjxoS8QAABYWkRuIq+srNT8+fM1f/58paSk6LTTTlN0dLS+++47FRUVSZKSkpL0l7/8RT169DC5WsA8Ho+hXYVlKquoVlJCnPqlJcnpbL1DLwCEu4gMUJmZmbrvvvv06aefateuXdq3b5+qqqrUqVMn5eTk6IILLtC1116rlJQUs0sFTJObX6yla3ertLy64bFkV5wmjUpXTkZ3EysDAPPRB6odud0ezsLzUf3ZR6WlFZa+bTWc/XAONm8/pGdWtHwH6vRx2YSodsDPgfmYA/OZPQcnz8JrfYdTRO6BAtAyj8fQ0rW7vY55Ze1uTqsHENEIUAAayd9f2uiyXXNKyqu1q7AsNAUBgAURoAA0Una8xrdxFd5DFgCEMwIUgEaSOsX6Ni4hrp0rAQDrIkABaCTjjGQlu7yHoxTXyZYGABCpCFAAGnE6HZo0Kt3rmImj0ukHBSCiEaAANJGT0V3Tx2U3WYlKccXRwgAAFKGNNAG0Lieju85N70YncgBoBgEKQIucTocyeyebXQYAWA6X8AAAAPxEgAIAAPATAQoAAMBPBCgAAAA/EaAAAAD8RIACAADwEwEKAADATwQoAAAAP9FIE5bk8Rh0wAYAWBYBCpaTm1+spWt3q7S8uuGxZFecJo1K5ww2AIAlcAkPlrJlZ7GeWZHXKDxJUml5tZ5Zkafc/GKTKgMA4L8IULAMt8fQy+/mex3zytrd8niMEFUEAEDzCFCwjB3ffK+SU1aeTlVSXq1dhWWhKQgAgBYQoGAZJceqfBpXVuE9ZAEA0N4IULCMlMQOPo1LSohr50oAAPCOAAXLyDqzi1Jc3sNRiutkSwMAAMxEgIJlRDkdum50htcxE0el0w8KAGA6AhQsZXBmd00fl63kU1aiUlxxmj4umz5QAABLoJEmLCcno7vOTe9GJ3IAgGURoGBJTqdDmb2TzS4DAIBmEaAAcfYeAMA/BChEPM7eAwD4i03kiGi5+Zy9BwDwHwEKEcvjMbR07W6vYzh7DwDQHAIUItauwrImK0+n4uw9AEBzCFCIWL6eqcfZewCAUxGgELF8PVOPs/cAAKciQCFi9UtLatLx/FScvQcAaA4BChHL6XRo0qh0r2M4ew8A0BwCFCJaTgZn7wEA/EcjTUQ8zt4DAPiLAAWIs/cAAP7hEh4AAICfCFAAAAB+IkABAAD4iQAFAADgJwIUAACAnwhQAAAAfiJAAQAA+IkABQAA4CcCFAAAgJ8irhO5YRj6/PPPtX79euXm5uqbb77R8ePH5XK5lJWVpbFjx+oXv/iFHA6O8QAAAM2LuAC1adMmTZ48ueHPaWlpSk1N1bfffquPP/5YH3/8sVavXq158+YpNjbWvEIBAIBlRVyAMgxDp59+um688Ub97Gc/U5cuXRqee+uttzR79my9//77evrpp3XnnXeaWCkAALCqiNsDNXDgQK1Zs0Y33HBDo/AkSWPHjtX06dMlSa+//ro8Ho8ZJQIAAIuLuADVqVMnxcTEtPj8hRdeKEkqKytTSUlJqMoCAAA2EnEBqjXV1dUN/9yhQwcTKwEAAFZFgDrF6tWrJUmZmZnq1KmTydUAAAArirhN5N5s375dr776qiRpypQpQXnN6Ggyqi+iopyN/h+hxxyYjzkwH3NgPrvMgcMwDMPsIqzgyJEj+uUvf6mioiL99Kc/1V//+tc2v6ZhGPSTAgAgDBGgJJWXl+uGG27Qjh071L9/fy1atCgol+/cbo+OHTsRhArDX1SUU4mJ8Tp27ITcbu5+NANzYD7mwHzMgfnMnoPExHifVr8i/hJeRUWFbr75Zu3YsUPp6en6xz/+EdS9T3V1/AD6w+328JmZjDkwH3NgPubAfFafA2tfYGxnJ06c0K233qovvvhCffr00Ysvvqjk5GSzywIAABYXsQGqurpa06ZN05YtW5SamqqFCxeqW7duZpcFAABsICIDVG1trWbMmKGNGzeqZ8+eWrhwoXr27Gl2WQAAwCYiLkC53W7deeed+uCDD9StWzctXLhQaWlpZpcFAABsJOI2kf/zn//UmjVrJEmxsbG65557Whw7e/ZsZWVlhao0AABgExEXoGpqahr++dtvv9W3337b4tjy8vJQlAQAAGwm4gLUVVddpauuusrsMgAAgI1F3B4oAACAtiJAAQAA+IkABQAA4CcCFAAAgJ8ibhM5AOvyeAztKixTWUW1khLi1C8tSU6nw+yyAKAJAhQAS8jNL9bStbtVWl7d8FiyK06TRqUrJ6O7iZUBQFNcwgNgutz8Yj2zIq9ReJKk0vJqPbMiT7n5xSZVBgDNI0ABMJXHY2jp2t1ex7yydrc8HiNEFQFA6whQAEy1q7CsycrTqUrKq7WrsCw0BQGADwhQAExVVuE9PPk7DgBCgQAFwFRJCXFBHQcAoUCAAmCqfmlJSnZ5D0cprpMtDQDAKghQAEzldDo0aVS61zETR6XTDwqApRCgAJguJ6O7po/LbrISleKK0/Rx2fSBAmA5NNIEYAk5Gd11bno3OpEDsAUCFADLcDodyuydbHYZANAqAhTaBWeaAQDCGQEKQceZZgCAcMcmcgQVZ5oBACIBAQpBw5lmAIBIQYBC0HCmGdB+PB5DO/eVatOOQ9q5r5RfRACTsQcKQcOZZkD7YF8hYD2sQCFoONMMCD72FQLWRIBC0HCmGRBc7CsErIsAhaDhTDMguNhXCFgXAQpBxZlmQPCwrxCwLjaRI+g40wwIDvYVAtZFgEK74EwzoO3q9xV6u4zHvkLAHFzCAwCLYl8hYF0EKACwMPYVAtbEJTwAsDj2FQLWQ4ACABtgXyFgLVzCAwAA8BMBCgAAwE8EKAAAAD+xBwpAxPJ4DDZmAwgIAQpARMrNL9bStbsbNalMdsVp0qh0WgMAaBWX8ABEnNz8Yj2zIq9Jh+/S8mo9syJPufnFJlUGwC4IUAAiisdjaOna3V7HvLJ2tzweI0QVAbAjAhSAiLKrsMzr2XKSVFJerV2FZaEpCIAtEaAARJSyCu/hyd9xACITAQpARElKiGt9kB/jAEQmAhSAiNIvLanJwbynSnGdbGkAAC0hQAGIKE6nQ5NGpXsdM3FUOv2gAHhFgAIQcXIyumv6uOwmK1EprjhNH5dNHygAraKRJoCIlJPRXeemd6MTOYCAEKAARCyn06HM3slmlwHAhriEBwAA4CcCFAAAgJ8IUAAAAH6KyD1Q3333nTZu3Kht27YpLy9PX331laqqqtS/f38tX77c7PIAAIDFRWSAWr16tebMmWN2GQAAwKYiMkB16tRJI0aMUHZ2trKzs1VQUKAnn3zS7LIAAIBNBBygtmzZosGDBwezlpC55pprdM011zT8mct2AADAHwFvIr/++ut1+eWX66WXXlJZWVkQSwIAALC2Nt2F98033+ixxx7ThRdeqFmzZmnLli3BqgsAAMCyAg5Q//rXv3TzzTerS5cuqqmp0apVq3TDDTewKtWOPB5DO/eVatOOQ9q5r1Qej2F2SQAARCSHYRht+lu4rq5O69ev17Jly/TJJ5/I4/HI4XAoNjZWl156qcaPH2/5vVLLly/XPffcE/Q2Bm63R8eOnQjKa23ZWayX381XSXl1w2MprjhdNzpDgzPtf/BpVJRTiYnxOnbshNxuj9nlRCTmwHzMgfmYA/OZPQeJifGKimp9fanNAeqHioqK9Nprr2nFihU6fPjwyTdwONS3b1+NHz9eY8eOVVJSUrDeLmjaK0AZhiGHo+0Hk27cWqQ5C1u+PHrPjYM1YmCvNr8PrMXtMbTjm+9VcqxKKYkdlHVmF0Vx0C0AWEJQA1Q9j8ej999/X6+//rr+/e9/y+12N1qVmjBhgnJycoL9tgGz8gqUx2No5ryPGq08nSolMU5P3na+rU+RN/s3Dn95PIby95eq7HiNkjrFKuOM5KB+/masONptDsIRc2A+5sB8Zs+BrytQ7dIHyul06uKLL1Z0dLTKysr0xRdfyDAMVVdX6+2339aqVav0P//zP/r973+vgQMHtkcJllFX17bJ37mv1Gt4kqSSY9XasbckLE6Vd7s9bf7M2ltufrGWrt2t0h/MS7IrTpNGpSsno+3hJje/WM+syGvyeEl5tea9sVXTx2UH5X1aYoc5CHfMgfmYA/NZfQ6CfhZecXGxnn32WV188cW69dZb9fnnn8swDOXk5Ojee+/VRRddJIfDoS+++EKTJk3S5s2bg11CWCmr8B6e/B2HtqkPN6WnhNrS8mo9syJPufnFbXp9j8fQ0rW7vY55Ze1ubiAAAJMFZQXKMAx98MEHeu211xou2RmGoU6dOunKK6/UhAkTlJ6eLulk/6jCwkL93//9nz7++GM99dRTWrp0aTDKCEtJCXFBHYfA+Rpuzk3vFvDlvF2FZU3C2alKyqu1q7AsLFYcAcCu2hSgDh06pDfeeENvvvmmDh06pPrtVFlZWZo4caJ+/vOfKz4+vsnXpaWl6amnntLw4cOVn5/flhLCXr+0JCW74rz+pZriilO/tKTQFRWhQhFuWHEEAHsIOEDdeuut+uijj+TxeGQYhuLj43XZZZdpwoQJPu1r6tSpk7p166aDBw8GWkJEcDodmjQqvdk9MfUmjkq39QZyuwhFuGHFEQDsIeAA9cEHH0iSzjrrLF177bUaN26cXC6XX68xevRoUxpuHjx4UGPHjm34c01NjSQpPz9fQ4cObXj85ptv1i233BLq8prIyeiu6eOym2xcTnHFaWKQNi6jdaEIN6w4AoA9BBygLr/8ck2YMEFDhgwJ+M3vvvvugL+2Ldxud7PBra6urtHjVVVVoSuqFTkZ3XVuejftKixTWUW1khJO/iXKylPohCLcsOIIAPbQLn2gcJLb7VFJSUVQX9PjMcIyREVHO5WcnKDS0gpL37baUouBesFqMdBcq4T2XnG0yxyEM+bAfMyB+cyeg5SUBPP6QKF9tHf/IbQuVJdTWXEEAGsjQNlESysf9f2H2ru5Iv4rVOHG6XTQqgAALIoAZQOh6D8E/xBuAHuoq/No/WcHVFx2Qt2T4nXxoNMVHR30HtKIQAQoG6C5IgD477X1u/XulkL9cKfvsg17NHpwmsZfnG5eYQgLBCgboLkiAPjntfW7tebTwiaPG4YaHidEoS1Yx7QBmisCLfN4DO3cV6pNOw5p575SzgmE6uo8endL0/D0Q+9uKeQuO7QJK1A2QHNFoHncmYrmrP/sgFpr0GMYJ8ddOuSM0BSFsMMKlA3UN1f0huaKiDT1d6ae+otF/Z2pufnFJlUGsxWXnQjqOKA5BCibqO8/lOxqfJkuxRVHCwNEHF/vTOVyXmTqntT0EPu2jAOaE/AlvKKiIr/Gx8XFyeVyKTY2NtC3jHg0V4xM4dp9vi24MxXeXDzodC3bsMfrZTyH4+Q4IFABB6hLLrkkoK9LS0vThRdeqOuvv169e/cO9O0jFv2HIgt7fJrHnanwJjraqdGD05q9C6/e6MFp9INCmwT83WMYRkD/279/v15++WVdeeWV+uc//xnMfxcgrLDHp2Xcmfpf3IXYvPEXp2vMkDQ5TlmsdTikMUPoA4W2C3gFat26ddq2bZv+8Ic/yOl0auLEiRoyZIh69OghwzBUXFysTz/9VK+++qrcbrcefPBBpaWladu2bVq0aJH27Nmju+++W1lZWaxEAaeg+7x3dr8zNViXZVmh9G78xem66sKz6ESOduEwjNZu9mze/v37dfXVV+v000/XCy+8oOTk5i8rlZaW6te//rWKioq0fPlypaamqqamRjfccIO+/PJLTZo0SbNnz27Tv4RVud0elZRUmF2GLZh9+rbV7NxXqsdf+bzVcXdNPDdol3TtNgctnQ9Zz6o3V3gLPUP79/R5Duz67291dvs5CEdmz0FKSoKioloP2QHH8Pnz5+v48eN68MEHWwxPkpScnKwHHnhAR48e1YIFCyRJsbGxuuOOO2QYhjZt2hRoCUDYYo9P6+x4Z2prl2W37PTtsix3IQLmC/gS3saNG9WxY0dlZ2e3OnbAgAHq2LGjPvroo4bHBg0apJiYGB08eDDQEoCwxR4f39jpzlRfQs/L/8rXJUP7tPpa3IUImC/gAFVSUqLoaN+/3DAMff/99w1/joqKUseOHVVVVRVoCUDYsvsen1Cyy52pPoWeY9Xa8c33Or2L9/5ErFAC5gv4El5KSoqqqqq0ZcuWVsdu2bJFJ06caHSpr7a2VseOHfN6+Q+IVHSfDz++hpmSY63/UhnpK5TceQgrCHgF6rzzztPy5ct133336fnnn1daWlqz4w4cOKD77rtPDodD559/fsPjBQUFMgxDvXr1CrQEIKzV7/E5dcNxiitOE7nLynZ8DTMpiR1aHRPJK5TceQirCDhA3XbbbXr33Xe1f/9+/eIXv9Bll12mIUOGqHv37nI4HCouLtbmzZu1Zs0anThxQgkJCZo2bVrD169evVqSNGTIkLb/WwBhyk57fOCdT6EnMU5ZZ3bRsaOVXl+rfoXS21144bhC2dKdh/Wb8K168wDCU8BtDCTp888/14wZM3TkyBE5Tu1W9v8YhqGuXbvqqaeeUk5OTsPj77zzjr777jv95Cc/0RlnhOdp2LQx8J3Zt62COQiF1loPzLhmoC4d3tfnOWhuNSZcVyg9HkOz5m9sddXt8akj2hQc+Tkwn9lz4Gsbg4BXoCTp3HPP1T//+U8tXrxY7777rvbs2SO32y3p5Cbxs88+W6NHj9avfvUrJSYmNvrayy+/vC1vDQC209pl2cGZ/oWeSFqh5M5DWE2bApQkuVwuTZs2TdOmTVNtba2OHj0qwzCUlJSkmJiYYNQIAGEj2KHHanchttfh19x5CKtpc4D6oZiYGHXt2jWYLwkAYcdqoSdY2nODd6TfeQjrCeqBQG63WyUlJSopKWm4lAcACH/tffh1/SZ8b8L1zkNYU5tXoE6cOKFXX31Vq1atUn5+fqM9UJmZmfr5z3+ua6+9VvHx3hvDAQDsKRSHX0fqnYewrjYFqG+++UZTp07V/v37derNfHV1dcrLy9P27dv1yiuvaMGCBerbt2+bigVgjvba1xIqdq/f6kK1wZveaLCSgAPU8ePHddNNN+ngwYOKjo7WT3/6U40YMUI9e/aUJB06dEiffPKJ/vWvf2nfvn266aab9PbbbyshISFoxQNof3ZvXGj3+u0glBu8I+nOQ1hbwAFq4cKFOnjwoLp3766//e1vOuecc5qM+eUvf6mdO3dqypQpOnjwoBYtWqSpU6e2qWAAoWP3xoV2r98uQr3BO1w34cNeAt5Evm7dOjkcDj3wwAPNhqd6mZmZevDBB2UYht57771A3w5AiPm6r8Wq55DZvX47YYM3IlHAAWrfvn2KjY3VRRdd1OrYCy+8UHFxcdq3b1+gbwcgxPzZ12JFdq/fTjj8GpEo4ABVV1fnc6NMh8OhmJgY1dXVBfp2AELM7o0LrV6/x2No575SbdpxSDv3ldp+Jax+g/epK1EprjgulSIsBbwHqmfPntq/f792796t9HTvv3ns2rVLx48fV+/evQN9OwAhZvfGhVauv6WN7b8anaFLh9v3bmU2eCOSBLwCNWzYMBmGoT/+8Y+qrm75N7jq6mr98Y9/lMPh0PDhwwN9OwAhZvd9LVat31vDyXlvbNXGrUUhrSfY6jd4D8vqqczeyYQnhK2AA9TNN9+s2NhY5ebm6oorrtDrr7+uAwcOqLa2VrW1tSosLNTrr7+uK664Qrm5uYqJidFNN90UzNoBtCO772uxYv2+bGz/+8o821/OAyKBwzi1A6Yf3nnnHd11112qq6uTw9H8f4QMw1B0dLQef/xxXX755QEXakdut0clJRVml2EL0dFOJScnqLS0QnV1HrPLiUgtzUFzl5t+2LjQ6k0qW6s/lHbuK9Xjr3ze6rh7fjVI6acntX9BaIL/FpnP7DlISUlQVFTr60tt6kR++eWXq3fv3po7d64+/vjjJt3InU6nLrjgAt1+++1eWx0AsC5v+1rs0KTSSvtyfN7YfrymnSsB0FZtPguvf//+ev7551VeXq7t27erpKREkpSSkqL+/fvL5XK1uUgA5mqucaGdmlRapfGizxvbO8W2cyUA2qrNAaqey+XSsGHDgvVyACwsFIfHhqP6je3e+lN1TYpXxhnJ7IMCLC7gTeQAIhdNKgPjy8b2W67MJnQCNkCAAuA3qzeptDJvDSdnXDNQIwb2MqkyAP7w6RJesDaAOxwO7dixIyivBcA8Vm5SaQctbWyPjY0yuzQAPvIpQLWh0wGAMOTLXh4rN9m0AqtsbAcQGJ8C1KJFi5p9/MCBA3r00UdVWVmp0aNHa9iwYerZs6ck6fDhw9q0aZPeffdddezYUb/73e+UmpoavMoBmKZ+L09zd+HVs3KTTQBoq4AbaX7//fcaO3asoqOj9dxzz7V4Ht6ePXs0ZcoUud1urVixQikpKW0q2E5opOk7sxunIbA5sFKTSjMFq5koPwfmYw7MZ/YctHsjzfnz5+vIkSP6+9//7vUw4bPPPlsPPPCAbr75Zs2fP1/33ntvoG8JwGKs1KTSLHZoJgog+AK+C+/9999XXFyczj///FbHnn/++erQoYM2bNgQ6NsBsKhIPjzW28HAz6zIU25+sUmVAWhvAQeo4uJiRUX5fsdIVFSUvvvuu0DfDgAsxddmojTEBMJTwAEqMTFRlZWVystreRNpvby8PFVUVHCsC4CwQTNRILIFHKCGDh0qwzA0e/ZslZaWtjiurKxMs2fPlsPh0NChQwN9OwCwFJqJApEt4E3k06dP13vvvaedO3fq8ssv18SJEzV06FD16NFDDodDhw4d0ubNm/Xqq6+qpKREcXFxmjZtWjBrb7NNmzbpxRdf1JdffqnKykr16tVLY8aM0ZQpU9SxY0ezywNgYTQTBSJbwG0MJOnf//63Zs6cqePHj8vhaH7jqGEYSkhI0JNPPqmRI0cGXGiwLV68WA8//LAMw1DPnj2VkpKiPXv2qKamRmeddZaWLl2qpKSkNr0HbQx8Z/Ztq2AO/OXxGJo1f2OrzUQfnzrC5431zIH5mAPzmT0HvrYxaNNZeBdeeKFWr16tCRMmKDExUYZhNPpfYmKiJkyYoFWrVlkqPOXl5emRRx6RJD3wwAN6//33tWLFCq1du1b9+/fX119/rdmzZ5tcJQAr8+VgYJqJAuGrTStQpyosLFRJSYkkKSUlRWlpacF66aCaNm2a1q1bp7Fjx+qxxx5r9FxBQYEuu+wyeTwerVy5UpmZmQG/DytQvjP7Nw4wB4EKZjNR5sB8zIH5zJ6Ddm+k2Zy0tDTLhqZ6FRUV+vDDDyVJ48ePb/J8nz59NGzYMG3cuFFr1qxpU4ACEP7CuZlosDqsA+EoqAHKDr766ivV1NQoNjZWAwcObHZMTk6ONm7cqC+//DLE1QGwo3A8GJgO64B3QQlQHo9HBQUFOnr0qOrq6ryOHTx4cDDeMmB79+6VJPXq1UsxMTHNjjnjjDMajQWASFLfYf1U9R3Wp4/LDjhEWX1Vy+MxtG3PERUePCpXfIzl6oN1tClAFRcX68knn9S7776rqqqqVsc7HA7t2LGjLW/ZZkePHpUkde7cucUx9c/Vj22L6Og27dOPGPXXm3257oz2wRyYzwpz4PEYeqW1DuvrdmvwOT38DhZbdhbr5XfzVXLKfrHrRmdocKb5q1pbdhbr5X/lq+SYNeuLFFb4OfBFwAHq8OHDGj9+vIqLi+XrPvQg7lcPWHX1yR+MllafJCk2NrbR2EA5nQ4lJye06TUiTWJivNklRDzmwHxmzsG2PUcaBZzmlByrVlFplQac3dXn1924tUjz3tja9LXKqzXvja2658bBGjGwl9/1BovV64tEVv9vUcAB6q9//asOHz6shIQE3X777brkkkvUvXt3v87HM0Nc3MmmdrW1tS2OqampaTQ2UB6PoWPHKtv0GpEiKsqpxMR4HTt2Qm43d76YgTkwnxXmoPCgbyvvhQeP6vQuvv0F5/EY+tvypuHkh/62YqsyUhNNuVxm9foijdk/B4mJ8e17F96///1vORwOPfzwwxozZkygLxNyvlye8+Uyn6+4DdY/breHz8xkzIH5zJwDV3zLq/OnjvO1xp37Sn1a1dqxt8SUzfhWry9SWf2/RQFfYCwpKVFUVJRGjRoVzHraXZ8+fSRJRUVFLa5C7d+/v9FYAIgU/dKSlOzyvvqe4jq5+dtXVj830Or1wZoCDlBdunRRhw4dFB1tr04IWVlZiomJUU1NjbZubX7JNjc3V5L0ox/9KISVAYD52qPDutXPDbR6fbCmgAPU8OHDVVFRoYKCgiCW0/4SEhJ0/vnnS5Jee+21Js8XFBRo06ZNkmSrS5MAECw5Gd01fVx2k5WoFFdcQC0M2mNVK5isXh+sKeAA9Zvf/Ebx8fF64okngllPSEybNk0Oh0MrV67UsmXLGu4OLC4u1syZM+XxeDRq1Ci6kAOIWDkZ3fWnqSN018RzNeWKLN018Vw9PnVEQP2frH5uoNXrgzW16Sy8zZs367e//a2ysrJ06623auDAgerYsWMw62s3L730kh599FEZhqHTTjtNycnJ2rNnj2pqatS3b18tXbpUKSkpbXoPzsLzndlnH4E5sIJwn4NgnhvYHnLzi/XK2t1N+lRZpb5IYfbPga9n4QUcoM455xy/v8YKjTR/6JNPPtELL7ygrVu3qrKyUr169dKYMWM0ZcoUJSS0vX8TAcp3Zv/AgDmwgkiYA6t3Inc6HSoqraITuYnM/jlo98OErdAUs62GDx+u4cOHm10GAEQMq58b6HQ6NODsrjq9S3zYhlgER8ABatGiRcGsAwAAwDYCDlBDhgwJZh0AAAC2Ya8mTgBszer7XwDAV0ELUIZhqLS0VFVVVerViwMXATTW3B1Yya44TeIOJwA21OYAtX37ds2fP18bN27UiRMnmtxpd/ToUf35z3+WJN13332KjY1t61sCsJnc/GI9syKvyeOl5dV6ZkVeQM0Z0b5YLQw9PnN7aVOAeuutt3Tfffeprq6uxTGdO3fWgQMH9Mknn+jiiy/WRRdd1Ja3BGAzHo+hpWt3ex3zytrdOje9G39ZWASrhaHHZ24/AXci//rrrzV79mzV1dXp+uuv15tvvqnk5OZvTb3yyitlGIbWrVsXcKEA7GlXYVmjvxSaU1JerV2FZaEpCF7VrxaeOmf1q4W5+cUmVRa++MztKeAA9eKLL6q2tlbXXXed7r33XvXv319RUVHNjh02bJgk6Ysvvgj07QDYFCfd24evq4Uej/37AFoFn7l9BRygNm3aJIfDoVtuuaXVsT169FB8fLyKiooCfTsANsVJ9/bBamHo8ZnbV8ABqri4WPHx8erZs6dP4+Pi4lRdzW+YQKThpHv7YLUw9PjM7SvgABUbG6va2lqfjnSpqqpSeXm5OnXqFOjbAbApTrq3D1YLQ4/P3L4CDlCpqamqq6tTQUFBq2M/+OADud1unX322YG+HQAby8norunjspusRKW44mhhYCGsFoYen7l9BdzG4IILLlB+fr4WLVqkP/zhDy2OKy0t1Z/+9Cc5HA6NHDky0LcDWkTvFHvIyeiuc9O7MVcWVr9a2FzPrnqsFgYXn7l9BRygJk+erKVLl+rVV19Vly5dNHny5EbPV1VV6b333tPcuXNVVFSk5ORkTZw4sa31Ao3QO8VenE6HMns33+4E1lC/Wnjqz1WKK04Tg/RzxS89jYXiM0fwOQxfNjG1YMOGDfrtb3+ruro6RUdHyzAMud1unXnmmSosLGzYIxUbG6u//e1vGj58eDBrtzy326OSkgqzy7CF6GinkpMTVFpaobo6j09f01J363pcGvJPIHOA4LLSHLRXyLH6Lz1mzgHB8iSzfw5SUhIUFdX6DqeA90BJ0k9+8hO9/PLL6t+/v2pra1VXVyfDMPT111+rpqZGhmEoKytLS5YsibjwhPZF7xSgfdWvFg7L6qnM3slBC080jGxZe3zmaD9tPgtv4MCBeuONN7Rz507l5uaquLhYHo9HXbt21aBBgzRgwIBg1Ak04k/vFC4ZAebjSB+EmzYHqHqZmZnKzMwM1ssBXtE7BbAXfulBuGnTJTzALPROAeyFX3oQbghQsCV6pwD2wi89CDcEKNgS3a0Be+GXHoQbAhRsKxy6W3s8hnbuK9WmHYe0c18pdw0ibPFLD8JN0DaRA2awc3drq/fDAYKNhpEIJwQo2J4du1u31AS0vh+OXVbQAH/Z+Zce4IcIUECI0Q8Hkc6Ov/QAp2IPFBBi/vTDAQBYEwEKCDH64QCA/RGggBCjHw4A2B97oIAQq++H4+0yHv1wAKB5Ho9hiZsQCFBAiNX3w2nuLrx69MMBgKas1P6FS3iACcKhCSgAhFJ9+5dTV+/r27/k5heHtB5WoACT0A8HAHxjxfYvBCjARPTDAQJXV+fR+s8OqLjshLonxeviQacrOpoLK+HIn/YvofpvKgEKAGA7r63frXe3FMr4wfGRyzbs0ejBaRp/sfcz92A/Vmz/QoACANjKa+t3a82nhU0eNww1PE6ICi9WbP/CWicAwDbq6jx6d0vT8PRD724pVF2dJ0QVIRTq2794E+r2LwQoAIBtrP/sQKPLds0xjJPjED7q2794E+r2LwQoAIBtFJedCOo42IfV2r+wBwoAYBvdk+KDOg72YqX2LwQoAIBtXDzodC3bsMfrZTyH4+Q4hCertH/hEh4AwDaio50aPTjN65jRg9PoB4V2xwoUAMBW6lsUnNoHyuEQfaAQMgQoAIDtjL84XVddeBadyGEaAhQAwJaio526dMgZZpeBCEWAAgAgBDwewxJ3jyE4CFAAALSz3PxiLV27u9GBuMmuOE0alR7y/kUIDi4WAwDQjnLzi/XMirxG4UmSSsur9cyKPOXmF5tUGdqCAAUAQDvxeAwtXbvb65hX1u6Wx9PK+TSwHAIUAADtZFdhWZOVp1OVlFdrV2FZaApC0BCgAABoJ2UV3sOTv+NgHRG3ibyqqkofffSRtm3bpry8POXl5amsrEyS9NlnnykhIcHcAgEAYSMpIa71QX6Mg3VEXIDau3evpk+fbnYZAHzErd+ws35pSUp2xXm9jJfiOvl9DXuJuAAVHR2tgQMHasCAAcrOzlaXLl00ZcoUs8sC0Axu/YbdOZ0OTRqVrmdW5LU4ZuKodH4psKGIC1Dp6el6/fXXG/584MABE6sB0JL6W79PVX/r9/Rx2YQo2EJORndNH5fd5JeBFFecJvLLgG1FXIACYH2+3vp9bno3fnOHLeRkdNe56d24HB1GCFAALMefW78zeyeHqCqgbZxOB9+vYYQ2BgAsh1u/AVgdK1DtLDqajOqLqChno/9H6FlpDrokdvB5XDj9jFlpDiIVc2A+u8wBAaodOZ0OJSfTV8ofiYnxZpcQ8awwB0M7d1SXt3fo+6NVLY7pmhSvof9zuqLCcA+JFeYg0jEH5rP6HNgmQN1///1atmyZ3183ZMgQLV68uB0qap3HY+jYsUpT3ttuoqKcSkyM17FjJ+R2e8wuJyJZbQ4m/bSf5r2xtcXnJ45K17Gj4fXzZbU5iETMgfnMnoPExHifVr9sE6BcLpe6du3q99d17ty5HarxXV0dP4D+cLs9fGYms8ocnHt2V6+3fp97dldL1NkerDIHkYw5MJ/V58A2AWrWrFmaNWuW2WUACCFu/QZgVbYJUEA44XgS33Hrt3/43gJCgwAFhBjHk6A9eDyGVm3cq/f+c0AVVXUNj/O9BbQPa98jCISZ+uNJTm0SWX88SW5+sUmVwc5y84v1/z39od76qKBReJL43gLaS0SuQI0bN05FRUWSJMMwGh6/+OKLG/550KBBmj9/fshrQ/jieBK0h5bODDwV31tAcEVkgDp69KjKysqaPP7Dx44fPx66ghAROJ4EweZLKK/H9xYQXBEZoNavX292CYhAHE+CYPMllP8Q31tA8LAHCgiRpIS4oI4D/A1EfG8BwUOAAkKkX1qSkl3e/wJLcZ287RzwhT+BiO8tILgIUECIOJ0OTRqV7nXMxFHpbPKFz3wJ5fX43gKCiwAFhFBORndNH5fd5C+9FFecpo/LplcP/OJLKO/UIZrvLaAdROQmcsBMHE+CYKoP5ac2Z03oEK2f/jhNPx/Rh+8toB0QoAATcDwJgolQDoQeAQoAwoBdQzln98GuCFAAAFNwLiTsjE3kAICQ41xI2B0BCgAQUr6eC+nxGF7HAGYiQAEAQsqfcyEBqyJAAQBCinMhEQ7YRA4AJmju7rNIwbmQCAcEKAAIsZbuPvvV6AxdOryviZWFRv0RNN4u43F2H6yOS3gAIpbHY2jnvlJt2nFIO/eVhmTTsre7z+a9sVUbtxa1ew1m41xIhANWoABEJDN6EPly99nfV+bpiWkj2uX9raSlI2hSXHGaSB8o2AABCkDEqV8FOlV9D6L2OnzXl7vPjpSdUP7+UqWfnhT09w+11rqMcwQN7IwABSCi+NqD6Nz0bkH/i9znu8+O1wT1fc3g6wqf0+lQv7SkhhC1q7CMEAVbIEABiCj+9CAK9tlyPt991ik2qO8bav6s8HGcC+yKTeQAIoqZPYjq7z7zpmtSvDLOsN+hwPX86TLOcS6wMwIUgIhiZg8iX+4+u+XKbFtfvvJ1hW/n/lKOc4GtEaAARBRfVoHaswdR/d1np9aQ4orTjGsGasTAXu3yvqHi68rdzn2lHOcCW2MPFICIUr8K1NwenXrt3YOopbvPYmOj2u09Q8XnlTsfP16Oc4FVsQIFIOJ4WwVqrxYGp3I6HcrsnaxhWT2V2TvZ1pftfsjXFb7MNN/2eXGcC6yKFSgAEYkeRO3D1xW+zN7JHOcCW2MFCkDECtdVILP5ssLHcS6wO1agAABB58sKH8e5wM4IUACAdlG/wucNl1JhVwQoAICpfAlagNWwBwoAAMBPBCgAAAA/EaAAAAD8xB4oAIBleTwGG8xhSQQoAIAl5eYXN2lxkOyK0yRaHMACuIQHALCc3PxiPbMir0mn8tLyaj2zIk+5+cUmVQacRIACAFiKx2No6drdXse8sna3PB4jRBUBTRGgAACWsquwzOsZeZJUUl6tXYVloSkIaAZ7oAAAllJW4T08+TvOytgkb18EKACApSQlxLU+yI9xVsUmeXvjEh4AwFL6pSUp2eU9HKW4Tq7W2BWb5O2PAAUAsBSn06FJo9K9jpk4Kt22l7rYJB8eCFAAAMvJyeiu6eOym6xEpbjiNH1ctq0vcbFJPjywBwoAYEk5Gd11bnq3sNtkHUmb5MMZAQoAYFlOp0OZvZPNLiOoImWTfLjjEh4AACEUCZvkIwEBCgCAEAr3TfKRggAFAECIhfMm+UjBHigAAEwQrpvkIwUBCgAAk4TjJvlIwSU8AAAAP0XcClRBQYHee+89bd68Wfn5+SotLVVcXJz69u2rSy+9VNddd50SEhLMLhMAAFhYRAUot9ut0aNHN/y5W7duysjI0JEjR7Rt2zZt27ZNr7/+ul566SWlpqaaWCkAALCyiApQhmGoU6dOmjBhgq666iqdddZZDc998cUXuuOOO7R//37dfvvteu2110ysFAAAWFlEBaioqCitW7dOSUlJTZ770Y9+pD/96U+aOHGivvzyS3311Vc655xzQl8kAACwvIjaRO5wOJoNT/UGDRokl8slSdq7d2+IqgIAAHYTUQGqNW63W3V1dZKkDh06mFwNAACwKgLUD6xbt04nTpxQdHS0fvSjH5ldDgAAsKiI2gPlzfHjx/XYY49Jkq6++mqlpKQE5XWjo8movoiKcjb6f4Qec2A+5sB8dp0Dj8dQ/v5SlR2vUVKnWGWckWzbjuZ2mQOHYRiG2UWYze12a+rUqfrggw+UmpqqlStXNuyFagvDMORw2PMbGABgDxu3Fum5t7bp+6NVDY916dxBU8YO0IiBvUysLLzZJkDdf//9WrZsmd9fN2TIEC1evLjF5w3D0L333qs333xTnTt31pIlS9SvX7+2lNrA7fbo2LETQXmtcBcV5VRiYryOHTsht9tjdjkRiTkwH3NgPrvNwZadxZr3xtYWn59xzUANzrTXwcRmz0FiYrxPq1+2uYTncrnUtWtXv7+uc+fOXp9/6KGH9OabbyohIUHPP/980MJTvbo66/8AWonb7eEzMxlzYD7mwHx2mAOPx9CSd/O9jnn53Xz9z5ldbHk5z+pzYJsANWvWLM2aNSuor/nYY49pyZIlio+P13PPPaeBAwcG9fUBwM48HkO7CstUVlGtpIQ49UtLsuVfxOFqV2GZSsurvY4pKa/WrsIyDixuB7YJUME2d+5cvfDCC4qLi9P8+fP14x//2OySAMAycvOLtXTt7kZ/QSe74jRpVLpyMux1SShclVV4D0/+joN/rL3FvZ0sWLBACxYsUExMjObNm6fhw4ebXRIAWEZufrGeWZHXZHWjtLxaz6zIU25+sUmV4YeSEuKCOg7+ibgAtWjRIs2dO1fR0dGaO3euRo4caXZJAGAZHo+hpWt3ex3zytrd8nhscf9RWOuXlqRkl/dwlOI6eekVwRdRl/AOHz6sRx55RJKUkJCgF154QS+88EKzY6+++mpdc801oSwPAEzHvhr7cDodmjQqXc+syGtxzMRR6exbaycRFaBqa2tV37Xh6NGj+uyzz1ocO2LEiFCVBQCWwb4ae8nJ6K7p47Kb7FdLccVpIvvV2lVEBajTTz9d+fneb/kEgEjGvhr7ycnornPTu3HHZIhFVIACAHhXv6/G22U89tVYj9Pp4JJqiEXcJnIAQMvq99V4w74agAAFADhF/b6aU+/wSnHFafq4bPbVAOISHgCgGeyrAbwjQAEAmsW+GqBlXMIDAADwEwEKAADATwQoAAAAPxGgAAAA/ESAAgAA8BMBCgAAwE8EKAAAAD8RoAAAAPxEgAIAAPCTwzAMw+wiwpVhGPJ4+Hh9FRXllNvtMbuMiMYcmI85MB9zYD4z58DpdMjhaP3IIgIUAACAn7iEBwAA4CcCFAAAgJ8IUAAAAH4iQAEAAPiJAAUAAOAnAhQAAICfCFAAAAB+IkABAAD4iQAFAADgJwIUAACAnwhQAAAAfiJAAQAA+IkABQAA4CcCFAAAgJ+izS4AaElBQYHee+89bd68Wfn5+SotLVVcXJz69u2rSy+9VNddd50SEhLMLjNsVVVV6aOPPtK2bduUl5envLw8lZWVSZI+++wzPvsg2rRpk1588UV9+eWXqqysVK9evTRmzBhNmTJFHTt2NLu8sPbdd99p48aNDd/nX331laqqqtS/f38tX77c7PLCnmEY+vzzz7V+/Xrl5ubqm2++0fHjx+VyuZSVlaWxY8fqF7/4hRwOh9mlNuEwDMMwuwjgVG63W1lZWQ1/7tatm3r06KEjR47o0KFDkqQzzjhDL730klJTU80qM6x99dVXGjt2bLPPEaCCZ/HixXr44YdlGIZ69uyplJQU7dmzRzU1NTrrrLO0dOlSJSUlmV1m2HrppZc0Z86cJo8ToELjk08+0eTJkxv+nJaWpsTERH377bcNv7BddNFFmjdvnmJjY80psgWsQMGSDMNQp06dNGHCBF111VU666yzGp774osvdMcdd2j//v26/fbb9dprr5lYafiKjo7WwIEDNWDAAGVnZ6tLly6aMmWK2WWFlby8PD3yyCOSpAceeEDjx4+Xw+HQ4cOHNXXqVG3fvl2zZ8/WvHnzTK40fHXq1EkjRoxQdna2srOzVVBQoCeffNLssiKGYRg6/fTTdeONN+pnP/uZunTp0vDcW2+9pdmzZ+v999/X008/rTvvvNPESptiBQqWZBiGjh492uJv3p999pkmTpwo6eQP2TnnnBPC6iLTgQMHdMkll0hiBSpYpk2bpnXr1mns2LF67LHHGj1XUFCgyy67TB6PRytXrlRmZqZJVUaW5cuX65577mEFKkSOHz+uuLg4xcTENPv8ggULNHfuXCUlJemTTz6R02mdrdvWqQT4AYfD4fWyxaBBg+RyuSRJe/fuDVFVQPBUVFToww8/lCSNHz++yfN9+vTRsGHDJElr1qwJaW1AqHTq1KnF8CRJF154oSSprKxMJSUloSrLJwQo2JLb7VZdXZ0kqUOHDiZXA/jvq6++Uk1NjWJjYzVw4MBmx+Tk5EiSvvzyy1CWBlhGdXV1wz9b7b/1BCjY0rp163TixAlFR0frRz/6kdnlAH6rXznt1atXi7+Bn3HGGY3GApFm9erVkqTMzEx16tTJ5GoaI0DBdo4fP96wX+Tqq69WSkqKyRUB/jt69KgkqXPnzi2OqX+ufiwQSbZv365XX31Vkix5AwsBCrbidrs1c+ZMHThwQKmpqZo1a5bZJQEBqb804W3/R/1t2z+8jAFEgiNHjui2225TbW2tfvrTn+pnP/uZ2SU1QRsDBN3999+vZcuW+f11Q4YM0eLFi1t83jAMzZ49Wx988IE6d+6sBQsWNGwkx3+11+eP4IqLi5Mk1dbWtjimpqam0VggEpSXl+uWW25RUVGR+vfvr0cffdTskppFgELQuVwude3a1e+v83YpQ5Ieeughvfnmm0pISNDzzz+vfv36BVpiWGuvzx/B5cvlOV8u8wHhpKKiQjfffLN27Nih9PR0/eMf/7Dc3qd6BCgE3axZs4J+ae2xxx7TkiVLFB8fr+eee67Fu5bQPp8/gq9Pnz6SpKKiItXW1jZ7KW///v2NxgLh7MSJE7r11lv1xRdfqE+fPnrxxReVnJxsdlktYg8ULG/u3Ll64YUXFBcXp/nz5+vHP/6x2SUBbZaVlaWYmBjV1NRo69atzY7Jzc2VJO40Rdirrq7WtGnTtGXLFqWmpmrhwoXq1q2b2WV5RYCCpS1YsEALFixQTEyM5s2bp+HDh5tdEhAUCQkJOv/88yWp2eOICgoKtGnTJknSmDFjQlobEEq1tbWaMWOGNm7cqJ49e2rhwoXq2bOn2WW1igAFy1q0aJHmzp2r6OhozZ07VyNHjjS7JCCopk2bJofDoZUrV2rZsmWqP1mruLhYM2fOlMfj0ahRozjGBWHL7Xbrzjvv1AcffKBu3bpp4cKFSktLM7ssn3AWHizp8OHDGjlypAzDUOfOnRsdJnyqq6++Wtdcc00Iq4sc48aNU1FRkaT/nk8oqdExO4MGDdL8+fPNKC8svPTSS3r00UdlGIZOO+00JScna8+ePaqpqVHfvn21dOlSep21o4MHD2rs2LENf66pqVFlZaWio6MbbV6++eabdcstt5hQYXhbtWqV7rjjDklSamqqevTo0eLY2bNnKysrK1SltYpN5LCk2traht/Gjx49qs8++6zFsSNGjAhVWRHn6NGjKisra/L4Dx87fvx46AoKQ5MnT1ZGRoZeeOEFbd26Vd9//7169eqlMWPGaMqUKRza3M7cbnez3+N1dXWNHq+qqgpdURGkvlWHJH377bf69ttvWxxbXl4eipJ8xgoUAACAn9gDBQAA4CcCFAAAgJ8IUAAAAH4iQAEAAPiJAAUAAOAnAhQAAICfCFAAAAB+IkABAAD4iQAFAADgJwIUAACAnzgLDwBs7tixY/r000+1fft27dixQ9u3b9d3330nSZozZ46uuuoqkysEwg8BCgBsbu3atbrnnnvMLgOIKAQoAAgD3bp10znnnKOsrCz1799fM2bMMLskIKwRoADA5q644gou0wEhxiZyAKb67rvv9MQTT+iKK65QTk6OBgwYoEsuuUT33nuv9uzZ02T8c889p4yMDGVnZ2vr1q3NvuYHH3ygzMxMZWRk6O23327y/ObNm/Xb3/5WF1xwgbKzszV06FDdeOONevPNN+V2u73Wu2XLFv3mN7/R0KFDNXDgQI0ePVpz585VRUWFli9froyMDF188cV+fQbvvPOOMjIylJGRoffee6/ZMV999ZUGDBigjIwMLViwoNFz0dH8LgyEmsMwDMPsIgBEpg0bNmjmzJmqrKyUJMXExCgmJqbRnx966CGNHTu24WsMw9Cvf/1rbdy4UWlpaXrrrbfUqVOnhueLi4t15ZVXqqSkRGPHjtVjjz3W6D3nzJmjl156SZLkcDjkcrlUUVHREJyGDRumZ555ptFr1lu8eLEefvhh1f9n0+VyqaqqSrW1tTrrrLM0fvx4zZkzR6mpqVq/fr1fn8Xvf/97vfnmm0pKStJbb72l0047reG5yspKXXXVVdq7d6+GDBmihQsXyun0/vtvRkZGw78vq1NA8LECBcAUW7du1YwZM1RZWalrr71W77zzjr788kt9/vnn2rBhgyZNmqTa2lrde++92rZtW8PXORwOPf744+rSpYsKCwv1hz/8oeE5wzB09913q6SkRL1799b999/f6D2XLFnSEJ6uvfZaffjhh9qyZYv+85//6J577lF0dLQ2bdqk2bNnN6n3s88+0yOPPCLDMHTeeedpzZo1+s9//qMvvvhCTz31lL7//ns9++yzAX8es2fP1plnnqmysjLNmjWr0UrYgw8+qL179yopKUlPPPFEq+EJQPvjpxCAKR544AHV1tZq2rRpeuCBB3TWWWcpKipKktSrVy/94Q9/0PXXX6+6ujrNnz+/0dd269ZNc+bMkcPh0KpVq7RixQpJ0t///ndt3LhRMTEx+vOf/6yEhISGr6mqqtK8efMkST//+c/1wAMPqFu3bpKkjh07avLkyfrd734n6eQltR+GNkl6+umn5fF4dPbZZ2vBggXq27evpJOXz8aMGaOnnnpKR48eDfjziI+P19y5cxUbG6stW7Y0/Du/8847Wr58uSTp4YcfVo8ePQJ+DwDBQ4ACEHI7d+7Utm3bFBMTo1//+tctjqu/dPfJJ5802Zs0cuRITZ48WdLJMLZy5Uo9/fTTkqTbb79dAwYMaDT+448/VllZmSTptttua/b9Jk2a1BCqVq9e3fB4WVmZNm3aJEm66aabFBsb2+Rrhw0bph//+Mct/rv4IjMzU3fddZck6dlnn9Vbb73VsIp23XXXadSoUW16fQDBw85DACGXm5srSfJ4PBozZkyL4+pDU2VlpcrKytSlS5dGz8+cObOhgWR98Dj//PObDWV5eXmSpNNOO61h9ehUUVFRGjZsmN5+++2G8dLJDdz1+54GDx7cYr1DhgzRf/7znxaf98X111+vjz/+WBs2bNDdd98tSerXr1/DPwOwBlagAIRccXGxpJMB6ciRIy3+r7S0tOFrTpw40eR1YmNj9eijjzb82eVy6dFHH5XD4Wgy9vvvv5ekVi+B9ezZs9F4SSopKWn4Z29fH6zLa4888oji4uIknQx1Tz75ZMOfAVgDK1AAQq5+ZenMM8/UP//5zza91muvvdbwz8ePH9dXX33VcBmuOc2Fq9bG+XqzcrBual65cqWqq6slnfyscnNzlZ6eHpTXBhAcrEABCLmuXbtKkgoLCxtaFgRiw4YNWrx4saSTt+0bhqHf/e53OnLkSJOx9Zf/Dh486PU1Dx06JElKSUlp8rXSf1fPmuPtOV9t375df/7znyU1bkXQXE8sAOYhQAEIuUGDBkmSamtrW2wc2Zri4uKG89+uuuoqvfzyy0pNTdX333+vu+++u8lqUHZ2tqSTAWnv3r3Nvqbb7dbmzZslqdEm9HPOOadhRerTTz9tsSZvz/misrJSM2fOVG1trYYNG6Y33nhDAwcOVFVVlWbOnKmampo2vT6A4CFAAQi5AQMGKCsrS5L0l7/8pdEeo+bU3z1Xz+Px6K677lJpaan69Omj2bNny+Vy6c9//rOio6P10Ucf6cUXX2z0Needd56SkpIkSX/961+bfZ9XX321YRXpZz/7WcPjSUlJGjp0qCTpxRdfbDbI1PeTaosHH3xQBQUFSkpK0uOPP67Y2NiGdgz5+flNmoICMA8BCkDIORwO/fGPf1RsbKyKior0y1/+UmvWrGm0Ufzw4cNauXKl/vd//1dPPPFEo69//vnn9cknnzT0e+rYsaMk6dxzz9X06dMlSU8++aS2b9/e8DUdOnRoOGB31apVuv/++xsu9Z04cUKLFy/WnDlzJEmXX355w4pVvRkzZsjhcGjXrl2aOnWqCgoKJEl1dXX617/+pRkzZqhz584BfyarV69u6Pf0yCOPNGxIP+OMMxqahS5ZskQbNmxo9utLSkoa/a9eZWVlo8eb24wPwH8c5QLANB9//LFmzpzZsMIUFRUll8ul6urqRn/R//KXv9RDDz0k6WQH8/ou5XfddZduuummRq/p8Xh044036tNPP1WfPn20YsWKhoAlNT3KJTExURUVFaqrq5MkDR06VM8++2yzR7m89NJLDSFLkhITE1VVVaWamhr169dPV199tebMmaO+fftqzZo1Pn8OBw4c0NixY1VeXq7rrruuSQd1Sbrrrru0cuVKJScna+XKlU3u+KvfL9Wa2267rSFIAggcK1AATHPeeefpvffe0x133KGcnBy5XC6Vl5fL4XDo7LPP1jXXXKP58+c3HK1y/Pjxhj1C5513XrP9npxOp/70pz8pKSlJBQUFevDBBxs9f88992jhwoUaPXq0unbtqsrKSiUkJGjo0KF65JFH9OKLLzYbniRp8uTJWrx4sUaOHKnOnTururpaqampmjp1qpYtW9YwLjEx0efPoK6uTnfccYfKy8u99nu6//771bt3b5WWluruu++Wx+Px+T0ABB8rUAAQJHfccYdWrVqlq6++Wo888ojZ5QBoR6xAAUAQ7N27t+GOwgsuuMDkagC0NwIUAPjoqaee0pIlS1RUVNRwCa2yslLvvPOObrjhBlVXV+vMM8/kzDogAnAJDwB8NG3aNK1bt06SFBMTo4SEBB07dqwhTPXo0UPPP/+8+vXrZ2aZAEKAo1wAwEeTJ09W9+7d9fnnn+u7777T0aNHlZCQoD59+uiiiy7Sr371q4ZeUwDCGytQAAAAfmIPFAAAgJ8IUAAAAH4iQAEAAPiJAAUAAOAnAhQAAICfCFAAAAB+IkABAAD4iQAFAADgJwIUAACAn/5/eLVo0lcXI1kAAAAASUVORK5CYII=",
+ "text/plain": [
+ "<Figure size 640x480 with 1 Axes>"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAHMCAYAAAAarpbgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABDzElEQVR4nO3deXjU1d3//9dkJYQJybAKoqAmQQgopLKJS5Uq0lqDuADW5XbBG6i9v2LBWsW22rpVRYsKWosIFJcqy6W4VBCtyiINIoQlgBIIBgiQTEiAbDOf3x/8JhKSTGbfPs/HdXldMnNm5ihnZl5zzvm8j8UwDEMAAADwWFy4OwAAABBtCFAAAABeIkABAAB4iQAFAADgJQIUAACAlwhQAAAAXiJAAQAAeIkABQAA4CUCFAAAgJcIUAAAAF5KCHcHwuXTTz/VF198oc2bN2v//v0qLy9XQkKCunfvrqFDh+q2225T9+7dw91NAAAQgSxmPQvv5ptv1tdff63ExER16tRJNptN5eXl2rdvn5xOp1JSUvTCCy9o+PDh4e4qAACIMKYNUEuWLFGXLl2Um5urpKSkhtv37Nmj3//+91q3bp0yMjL06aefqm3btmHsKQAAiDSmDVDuHDp0SBdeeKEk6ZVXXtEll1wS5h4BAIBIYto9UO507NhR6enpstvtqq6u9vl5DMOQ00k+NbO4OAtjwOQYA2AMRJe4OIssFkur7QhQzfjuu+9kt9sVFxenPn36+Pw8FotFR44cU329M4C9Q7RISIhTRkYqY8DEGANgDEQfmy1V8fEEKI8ZhqGysjLl5+fr6aefliTdfvvt6tGjR5h7BgAAIo3pA9TSpUs1bdq0RredddZZevrpp3X11Vf7/fzx8ZTaMivX3z1jwLwYA2AMxC7TB6gOHTpo4MCBMgxD+/fv14EDB1RUVKT33ntPF1xwgbp27erX86elpQSop4hWjAEwBsAYiD1chXeK4uJiPfHEE1q+fLm6dOmiZcuWyWq1+vx8R44cl8PBurcZxcfHKS0thTFgYowBMAaiT1paikczhqafgTpVjx499Le//U3XXHONduzYoQULFmjixIk+P5/D4WTjoMkxBsAYAGMg9rAo24z4+HhddNFFkqSCgoIw9wYAAEQaAlQL6uvrJUlOJ78YAABAYwSoZtTW1uqzzz6TJL/qQAEAgNhkygC1adMmPffccyoqKmpy365duzRx4kTt2bNHbdu21Q033BD6DgIAgIhmyk3kx44d06xZszRr1izZbDaddtppSkhI0MGDB1VSUiJJSk9P13PPPacuXbqEubcAgGjndBratrtc9qM1Sk9NVlaPdMXFtV7tGpHLlAGqd+/eeuihh/T1119r+/bt2r17t6qrq9WuXTvl5ubqoosu0o033iibzRburgIAotyqjSV6edFGlVXWNNyWYU3W+BGZys3uHMaewR/UgQqy8vKjXLpqUq4zsBgD5sUYwDc7D2nmOxtbvH/y6BxCVIQ5cRZe6zucTLkHCgCAYHM6Df3z40K3bd5YvkNOJ/MY0YgABQBAEGwvtjdatmtOWWWNthfbQ9MhBBQBCgCAILAfdR+evG2HyEKAAgAgCNJTkwPaDpGFAAUAQBBk9UiXzeo+HNmsJ0oaIPoQoAAACIK4OItuujLbbZtxIzKpBxWlCFAAAATJBb0764FbL2gyE2WzJlPCIMqZspAmAAChMqx/N2V3T9OWXWVUIo8hBCgAAIIsLs6i3mdmhLsbCCCW8AAAALxEgAIAAPASAQoAAMBLBCgAAAAvEaAAAAC8RIACAADwEgEKAADASwQoAAAAL1FIEwAABJ3TaWh7sT1mqrEToAAAQFDlF5Zq4fIdKq+sabgtw5qs8SMyo/Y8QJbwAABA0OQXlurFxQWNwpMklVfW6MXFBcovLA1Tz/xDgAIAAEHhdBpauHyH2zZvLN8hp9MIUY8ChwAFAACCYnuxvcnM06nKKmu0vdgemg4FEAEKAAAEhf2o+/DkbbtIQoACAABBkZ6aHNB2kYQABQAAgiKrR7oyrO7Dkc16oqRBtCFAAQCAoIiLs2j8iEy3bcaNyIzKelAEKAAAEDS52Z01eXROk5komzVZk0fnRG0dKAppAgCAoMrN7qwBmZ2oRA4AAOCNuDiLep+ZEe5uBAwBCgDQINbOKwOChQAFAJAUm+eVAcHCJnIAQMyeVwYECwEKAEwuls8rA4KFAAUAJhfL55UBwUKAAgCTi+XzyoBgIUABgMnF8nllQLAQoADA5GL5vDIgWAhQAGBysXxeGRAsBCgAQMyeVwYEC4U0AQCSYvO8MiBYCFAAgAaxdl4ZECws4QEAAHiJAAUAAOAlAhQAAICXCFAAAABeIkABAAB4iQAFAADgJQIUAACAlwhQAAAAXiJAAQAAeMl0lcgNw9A333yjTz/9VPn5+fr+++9VVVUlq9WqPn36KC8vT1dffbUsFo4uAAAAzTNdgFqzZo1uu+22hj/36NFD3bt31w8//KCvvvpKX331lZYtW6aZM2cqKSkpfB0FAAARy3QByjAMnX766br11lv185//XB06dGi4b8mSJZo+fbo+++wz/e1vf9Nvf/vbMPYUAABEKothGEa4OxFKVVVVSk5OVmJiYrP3z549WzNmzFB6erpWr16tuDj/tomVlx9Vfb3Tr+dAdEpIiFNGRipjwMQYA2AMRB+bLVXx8a1/95tuE3m7du1aDE+SdPHFF0uS7Ha7ysrKQtUtAAAQRUwXoFpTU1PT8O9t2rQJY08AAECkIkCdYtmyZZKk3r17q127dmHuDQAAiESm20TuzubNm/Xmm29KkiZMmBCQ5/RkHRWxyfV3zxgwL8YAGAOxy3SbyFty6NAhXX/99SopKdHPfvYzvfDCC+HuEgAAiFAEKEmVlZW65ZZbtGXLFvXt21fz5s0L2PLdkSPH5XBw5YUZxcfHKS0thTFgYowBMAaiT1paikczhqZfwjt69KjuvPNObdmyRZmZmfrHP/4R0L1PDoeTS1dNjjEAxgAYA7HH1Iuyx48f1913360NGzaoZ8+eeu2115SRkRHubgEAgAhn2gBVU1OjSZMmad26derevbtef/11derUKdzdAgAAUcCUAaqurk733HOPVq1apa5du+r1119X165dw90tAAAQJUwXoBwOh37729/q888/V6dOnfT666+rR48e4e4WAACIIqbbRP7hhx/qo48+kiQlJSXpgQceaLHt9OnT1adPn1B1DQAARAnTBaja2tqGf//hhx/0ww8/tNi2srIyFF0CAABRhjpQQcYJ3ObFKexgDIAxEH1stlSP6kCZbg8UAACAvwhQAAAAXiJAAQAAeIkABQAA4CXTXYUHALHO6TS0vdgu+9EapacmK6tHuuLiLOHuFhBTCFAAEEPyC0u1cPkOlVfWNNyWYU3W+BGZys3uHMaeAbGFJTwAiBH5haV6cXFBo/AkSeWVNXpxcYHyC0vD1DMg9hCgACAGOJ2GFi7f4bbNG8t3yOmk9B8QCAQoAIgB24vtTWaeTlVWWaPtxfbQdAiIcQQoAIgB9qPuw5O37QC4R4ACgBiQnpoc0HYA3CNAAUAMyOqRrgyr+3Bks54oaQDAfwQoAIgBcXEWjR+R6bbNuBGZ1IMCAoQABQAxIje7syaPzmkyE2WzJmvy6BzqQAEBRCFNAIghudmdNSCzE5XIgSAjQAFAjImLs6j3mRnh7gYQ0whQAADEMM5GDA4CFAAAMYqzEYOHTeQAAMQgzkYMLgIUAAAxhrMRg48ABSConE5D23aXa82W/dq2u5wPbCAEOBsx+NgDBSBo2H8BhAdnIwYfM1AAgoL9F0D4cDZi8BGgAAQc+y+A8OJsxOAjQAEIOPZfAOHF2YjBR4ACEHDsvwDCj7MRg4tN5AACjv0XQGTgbMTgIUABCDjX/gt3y3jsvwBCg7MRg4MlPAABx/4LALGOAAUgKNh/ASCWsYQHIGjYfwEgVhGgAAQV+y8AxCKW8AAAALxEgAIAAPASAQoAAMBL7IECACAGOJ0GF2yEEAEKAIAol19YqoXLdzQqXpthTdb4EZmUDAkSlvAAAIhi+YWlenFxQZPK/+WVNXpxcYHyC0vD1LPYRoACACBKOZ2GFi7f4bbNG8t3yOk0QtQj8yBAAQAQpbYX292eOSlJZZU12l5sD02HTIQABQBAlLIfdR+evG0HzxGgAACIUumpya038qIdPEeAAgAgSmX1SG9yYPepbNYTJQ0QWAQoAACiVFycReNHZLptM25EJvWggoAABQBAFMvN7qzJo3OazETZrMmaPDqHOlBBQiFNAACiXG52Zw3I7EQl8hAiQAEAEAPi4izqfWZGuLthGizhAQAAeIkABQAA4CUCFAAAgJdMuQfq4MGDWrVqlTZt2qSCggJt3bpV1dXV6tu3rxYtWhTu7gEAgAhnygC1bNkyPf744+HuBgAAiFKmDFDt2rXTsGHDlJOTo5ycHBUVFenZZ58Nd7cAAECU8DlArVu3ThdccEEg+xIy1113na677rqGP7NsBwAAvOHzJvKbb75Zo0aN0ty5c2W32wPYJQAAgMjm11V433//vZ588kldfPHFmjp1qtatWxeofgEAAEQsnwPUv//9b915553q0KGDamtr9f777+uWW25hVgoAAASN02lo2+5yrdmyX9t2l8vpNMLSD4thGH69cn19vT799FO99dZbWr16tZxOpywWi5KSknTFFVfohhtuiPi9UosWLdIDDzwQlDIGR44cl8PhDOhzIjrEx8cpLS2FMWBijAEwBgJr3bZS/fPjQpVV1jTcZrMm66Yrs3VB78AcmpyWlqL4+Nbnl/wOUCcrKSnR22+/rcWLF+vAgQMnXsBiUa9evXTDDTcoLy9P6enpgXq5gAlmgAIQexxOQ1u+P6yyI9WypbVRn7M6KJ5DW4GgWrWxRI+/3vJWoQduvUDD+ncLWX8CGqBcnE6nPvvsM/3rX//Sf/7zHzkcjkazUmPHjlVubm6gX9ZnzEAh0JxOQzt+qFBNvaHkBIsyu7fnVPQY4c0vYGYfwBgIDKfT0JSZXzZ6353KlpasZ3893O/PWk9noIJSByouLk6XXXaZEhISZLfbtWHDBhmGoZqaGr333nt6//33dd555+n3v/+9+vfvH4wuRAyHw6n6et40ZpJfWKqFy3eo/KQ3eoY1WeNHZCo3OzBTzAiP/MJSvbi4oMntZZU1mvnORk0endPs3zGfA2AM+Gfb7nK34UmSyo7UaMuuMvU+MyMkfQr4WXilpaV66aWXdNlll+nuu+/WN998I8MwlJubqwcffFCXXnqpLBaLNmzYoPHjx2vt2rWB7gIQNq4v2PJT3ujllTV6cXGB8gtLw9Qz+MvpNLRw+Q63bd5YviNsG1qBWGY/6j48edsuEAIyA2UYhj7//HO9/fbbDUt2hmGoXbt2uuaaazR27FhlZmZKOlE/qri4WH/84x/11Vdf6fnnn9fChQsD0Q0grDz9gh2Q2YnlvCi0vdjeJBifqqyyRtuL7SH7BQyYRXpqckDbBYJfAWr//v1655139O6772r//v1ybafq06ePxo0bp1/84hdKSUlp8rgePXro+eef19ChQ1VYWOhPF4CIwRdsbIvEX8CAWWT1SFeGNdntZ6zNmqysHukh65PPAeruu+/Wl19+KafTKcMwlJKSoquuukpjx471aF9Tu3bt1KlTJ+3bt8/XLgARhS/Y2BaJv4ABs4iLs2j8iMxm9yC6jBuRGdLZfZ8D1Oeffy5JOvvss3XjjTdq9OjRslqtXj3HlVdeGZaCm/v27VNeXl7Dn2trayVJhYWFGjx4cMPtd955p+66665Qdw9Rii/Y2BaJv4ABM8nN7qzJo3OaXKRjsyZrXBgu0vE5QI0aNUpjx47VoEGDfH7x+++/3+fH+sPhcDQb3Orr6xvdXl1dHbpOIerxBRvbIvEXMGA2udmdNSCzk7YX22U/WqP01BOfqeF43wWlDhR+VF5+lEtXTaSly9xdWrrMHdGjuTIVLf0CTkiIU0ZGaqPPAafTiIgPf4RGc2MAkc1mSw19JXI0xZvGfLz5gkV08jQEnfrlSY0w8yFARR8CVITgTWNOTqeh70oqVGdYlGgxdHY3KpGb0clfnms372d20oQIUNHH0wAVlErkgNnFxVl0bk8bH5yQRI2wSOfPsipLsuZFgAKAICvcU06NsAjlz7IqS7LmFvCjXAAAjdmraj1rR42wkPLn6CWObQIBCgCCLL1dkmftqBEWMv6cbRjN5yI6nYa27S7Xmi37tW13eUT2MVqwhAcAQZZ9RgY1wiKMP0cvefPYnLM7+N3XQGHJMbCYgQKAIHMV4XSHIpyh5c/RS9F4bBNLjoFHgAKAEHAdQ5FhbbxMZ7MmU8IgDPw5einajm2K5iXHSObzEl5JSYlX7ZOTk2W1WpWU5NleAACINZF0DIXZ+XP0UrQd2+TPcmUwRXsJCJ8D1OWXX+7T43r06KGLL75YN998s84880xfXx4AolJcnIVSBRHAn7MNo+1cxEhccoyF/Vg+L+EZhuHTP3v27NE///lPXXPNNfrwww8D+d8CAAHFFUuxzZ9l1Whako20JcdY2Y/l8wzUihUrtGnTJv3hD39QXFycxo0bp0GDBqlLly4yDEOlpaX6+uuv9eabb8rhcOjRRx9Vjx49tGnTJs2bN087d+7U/fffrz59+jATBSDixMIvZCn6l0mCzZ9l1WhZko2kJcdYqsrv81l4e/bs0ZgxY3T66adrzpw5yshofkq6vLxct99+u0pKSrRo0SJ1795dtbW1uuWWW/Ttt99q/Pjxmj59ul//EZGMYzzMizOwopfrF3JLPJ1hCPcYiJUQGM3CPQZcAjWm/bVtd7meeuObVttNGzcgbEvdnp6F5/MS3qxZs1RVVaVHH320xfAkSRkZGXrkkUdUUVGh2bNnS5KSkpJ03333yTAMrVmzxtcuAEDAxcoVS7GyTILAiJQlx0jcj+Urn5fwVq1apbZt2yonJ6fVtv369VPbtm315ZdfNtw2cOBAJSYmat++fb52AQACLlKvWPJGLC2TIHAiYckx0vZj+cPnAFVWVqaEBM8fbhiGDh8+3PDn+Ph4tW3bVtXV1b52AQACLhZ+IcdCCERwhPsq0Ejaj+Uvn5fwbDabqqurtW7dulbbrlu3TsePH2+01FdXV6cjR464Xf4DgFALxS/kYF/dFwshELEplqry+zwDdeGFF2rRokV66KGH9Oqrr6pHjx7Nttu7d68eeughWSwWDR8+vOH2oqIiGYahbt26+doFAAi4YP9CDsXG7lhaJkHsce3HOvV9YLMma1wUXeDgc4D69a9/rY8//lh79uzR1VdfrauuukqDBg1S586dZbFYVFpaqrVr1+qjjz7S8ePHlZqaqkmTJjU8ftmyZZKkQYMG+f9fAQABEswiiS1dCeXa2B2ozbyxtEyC2BQJ+7H85XMZA0n65ptvdM899+jQoUOyWJr/jzYMQx07dtTzzz+v3Nzchts/+OADHTx4UD/96U91xhln+NqFiBfuS1cRPpFy+TJ809xMkbe/kE8eA7W1Dk2dtarVUPPUxGEB+RKJlMvWzY7PgejjaRkDvwKUJFVWVmr+/Pn6+OOPtXPnTjkcDkknNomfc845uvLKK/WrX/1KaWlp/rxM1OJNY158cEY/f4tQnjwGCr47HPL6N4EIgfAPnwPRx9MA5fMSnovVatWkSZM0adIk1dXVqaKiQoZhKD09XYmJif4+PQCETSCvWArHxu5YWCaJJFR1x8n8DlAnS0xMVMeOHQP5lAAQE8K1sTvcl63HCqq641Q+lzFojsPhUFlZmcrKyhqW8gAAP27sdoeN3ZGJqu5ojt8zUMePH9ebb76p999/X4WFhY32QPXu3Vu/+MUvdOONNyolJcXvzgJAtArm1X0IHqq6oyV+Bajvv/9eEydO1J49e3TqXvT6+noVFBRo8+bNeuONNzR79mz16tXLr84CiB2u/STllTWqPFardqmJsrVrE9P7SmKl/o2ZUNUdLfE5QFVVVemOO+7Qvn37lJCQoJ/97GcaNmyYunbtKknav3+/Vq9erX//+9/avXu37rjjDr333ntKTU0NWOcBRKfm9pO4xPq+EjZ2RxequqMlPgeo119/Xfv27VPnzp318ssv69xzz23S5vrrr9e2bds0YcIE7du3T/PmzdPEiRP96jCA6NZafaJAF5WMRGzsjh5UdUdLfN5EvmLFClksFj3yyCPNhieX3r1769FHH5VhGPrkk098fTkAMcCT/SQubyzfEfAz4gBvsfkfLfE5QO3evVtJSUm69NJLW2178cUXKzk5Wbt37/b15QDEAE/2k7i49pUA4RRLh98isHwOUPX19R4XyrRYLEpMTFR9fb2vLwcgBni7T4R9JYgErs3/p85E2azJMb3UDPd83gPVtWtX7dmzRzt27FBmpvt0vn37dlVVVenMM8/09eUAxABv94mUlh0PUk8A77D5H6fyeQZqyJAhMgxDf/rTn1RT0/KvxJqaGv3pT3+SxWLR0KFDfX05ADHAk/0kJ1vy5S6KFCJiuDb/D+nTVb3PzCA8mZzPAerOO+9UUlKS8vPz9ctf/lL/+te/tHfvXtXV1amurk7FxcX617/+pV/+8pfKz89XYmKi7rjjjkD2HUCU8WQ/yanYTA4gElmMUytgeuGDDz7QtGnTVF9fL4ul+SRuGIYSEhL01FNPadSoUT53NFpxArd5cQp7y9zVgWrOtHEDovKyf8ZA6EXagb+Mgehjs6UqPr71+SW/KpGPGjVKZ555pmbMmKGvvvqqSTXyuLg4XXTRRbr33nvdljoAYC6u/SRLvvhe769u/epcNpPDExz4i1Dy+yy8vn376tVXX1VlZaU2b96ssrIySZLNZlPfvn1ltVr97iSA2BMXZ1GfnjaPAhRFCtGalgq0mqEwK8LD7wDlYrVaNWTIkEA9HQATcG0qd7eUR5FCtIYDfxEOPm8iBwB/UaQQgeDNgb9AoBCgAIQVRQrhLw78RTh4tIQXqA3gFotFW7ZsCchzAYgdFCmEPzjwF+HgUYDyo9IBAHjEVaQQ8BZ76RAOHgWoefPmNXv73r179cQTT+jYsWO68sorNWTIEHXt2lWSdODAAa1Zs0Yff/yx2rZtq9/97nfq3r174HoOAIB+3EvX3FV4LuylQ6D5XEjz8OHDysvLU0JCgl555ZUWz8PbuXOnJkyYIIfDocWLF8tms/nV4WhD8TTzCmUBvUgrHogTKKIYWs3VgbJZkzUujHWgGAPRJ+iFNGfNmqVDhw7p73//u9vDhM855xw98sgjuvPOOzVr1iw9+OCDvr4kgGZQPBA4gb10CCWfr8L77LPPlJycrOHDh7fadvjw4WrTpo1Wrlzp68sBaIareOCpez9cxQM5iBdmw4G/CBWfA1Rpaani4+M9bh8fH6+DBw/6+nIATuFp8UAO4gWAwPM5QKWlpenYsWMqKGh5055LQUGBjh49yrEuQABRPBAAwsfnADV48GAZhqHp06ervLy8xXZ2u13Tp0+XxWLR4MGDfX05AKegeCAAhI/Pm8gnT56sTz75RNu2bdOoUaM0btw4DR48WF26dJHFYtH+/fu1du1avfnmmyorK1NycrImTZoUyL77bc2aNXrttdf07bff6tixY+rWrZtGjhypCRMmqG3btuHuHuAWxQMBIHx8LmMgSf/5z380ZcoUVVVVyWJpfqOeYRhKTU3Vs88+q0suucTnjgba/Pnz9Ze//EWGYahr166y2WzauXOnamtrdfbZZ2vhwoVKT0/3+3W4dNW8gn35stNpaOqsVa0WD3xq4jA20oYJl7CDMRB9PC1j4NdZeBdffLGWLVumsWPHKi0tTYZhNPonLS1NY8eO1fvvvx9R4amgoECPPfaYJOmRRx7RZ599psWLF2v58uXq27evvvvuO02fPj3MvQTc4yBeAAgfv2agTlVcXKyysjJJks1mU48ePQL11AE1adIkrVixQnl5eXryyScb3VdUVKSrrrpKTqdTS5cuVe/evf16LX51mFeofnlGYvFAnMDsAxgD0SfohTSb06NHj4gNTS5Hjx7VF198IUm64YYbmtzfs2dPDRkyRKtWrdJHH33kd4ACgo3igYg2VM5HLAhogIoGW7duVW1trZKSktS/f/9m2+Tm5mrVqlX69ttvQ9w7wDccxItoQeV8xIqABCin06mioiJVVFSovr7ebdsLLrggEC/ps127dkmSunXrpsTExGbbnHHGGY3aIjj4FQqYi6ty/qlclfMnj86JyBDFZxWa41eAKi0t1bPPPquPP/5Y1dXVrba3WCzasmWLPy/pt4qKCklS+/btW2zjus/V1h+erKOa0bptpfrnx4UqO2Xfzk1XZuuC3pH3AeoL1989Y8C8GAM/cjoNvdFa5fwVO3TBuV0iKpz4+1nFGIhdPgeoAwcO6IYbblBpaak83YcewP3qPqupOfEmaGn2SZKSkpIatfVHWlqK388Ra1ZtLNHMdzY2ub2sskYz39moB269QMP6dwtDz4KDMQDGgLRp56FGIaQ5ZUdqVFJerX7ndAxRr9wL5GcVYyD2+BygXnjhBR04cECpqam69957dfnll6tz585enY8XDsnJJ4oK1tXVtdimtra2UVt/HDlyXA4HV164OJ2GXl7U9APpZC8v3qjs7mkR9SvUF/HxcUpLS2EMmBhj4EfF+zyb0S/eV6HTO4Q/bATqs4oxEH3S0lKCexXef/7zH1ksFv3lL3/RyJEjfX2akPNkec6TZT5PORxOLl09ybbd5R79Ct2yqyxmNkUzBsAYkKwpLc/6n9ouEv5fBfqzijEQe3xelC0rK1N8fLxGjBgRyP4EXc+ePSVJJSUlLc5C7dmzp1FbBA7ntwHmlNUjXRlW97P6NuuJDdqRgM8qtMbnANWhQwe1adNGCQnRVQmhT58+SkxMVG1trTZubH56Nj8/X5J0/vnnh7Bn5sD5bYA5RVvlfD6r0BqfA9TQoUN19OhRFRUVBbA7wZeamqrhw4dLkt5+++0m9xcVFWnNmjWSFFVLk9Ei2n6FAgic3OzOmjw6p8lngM2aHHElDPisQmt8DlD/+7//q5SUFD399NOB7E9ITJo0SRaLRUuXLtVbb73VcHVgaWmppkyZIqfTqREjRlCFPAii7VcogMDKze6sv04cpmnjBmjCL/to2rgBemrisIgKTxKfVWidX2fhrV27Vr/5zW/Up08f3X333erfv7/atm0byP4Fzdy5c/XEE0/IMAyddtppysjI0M6dO1VbW6tevXpp4cKFstlsfr8O5x81zwznt3EGFhgD0c/fzyrGQPTx9Cw8nwPUueee6/VjIqGQ5slWr16tOXPmaOPGjTp27Ji6deumkSNHasKECUpNTQ3Ia/CmaVmsV/flgxOMgdjgz2cVYyD6BP0w4UgoiumvoUOHaujQoeHuhmlxfhuAaMBnFZrjc4CaN29eIPsBAAAQNXwOUIMGDQpkPwAAAKJGdBVxAiJArO/dAgC0LmAByjAMlZeXq7q6Wt26xc5BsMDJmrsiJ8OarPExdPUgAKB1fgeozZs3a9asWVq1apWOHz/e5Eq7iooKPfPMM5Kkhx56SElJSf6+JBAW+YWlenFxQZPbyytr9OLigogrBGgGzAbCE4wTBINfAWrJkiV66KGHVF9f32Kb9u3ba+/evVq9erUuu+wyXXrppf68JBAWTqehhct3uG3zxvIdGpDZiQ/mEGE2EJ5gnCBYfK5E/t1332n69Omqr6/XzTffrHfffVcZGc1f5nnNNdfIMAytWLHC544C4bS92N7oA7g5ZZU12l5sD02HTM41G3jq34lrNjC/sDRMPUMkYZwgmHwOUK+99prq6up000036cEHH1Tfvn0VHx/fbNshQ4ZIkjZs2ODrywFhxcnskcPT2UCnM/pr1cF3jBMEm88Bas2aNbJYLLrrrrtabdulSxelpKSopKTE15cDwoqT2SMHs4HwBOMEweZzgCotLVVKSoq6du3qUfvk5GTV1PDrHNGJk9kjB7OB8ATjBMHmc4BKSkpSXV2dR0e6VFdXq7KyUu3atfP15YCw4mT2yMFsIDzBOEGw+Rygunfvrvr6ehUVFbXa9vPPP5fD4dA555zj68sBYZeb3VmTR+c0mYmyWZMpYRBCzAbCE4wTBJvPZQwuuugiFRYWat68efrDH/7QYrvy8nL99a9/lcVi0SWXXOLrywERITe7swZkdqKmTBi5ZgObq8nl4m42kJpA5uDvOAFa43OAuu2227Rw4UK9+eab6tChg2677bZG91dXV+uTTz7RjBkzVFJSooyMDI0bN87f/gJhx8nsoXdq6BmQ2UmTR+c0qe9jsyZrnJv6PtQEMhfXrLG34wTwhMXwZBNTC1auXKnf/OY3qq+vV0JCggzDkMPh0FlnnaXi4uKGPVJJSUl6+eWXNXTo0ED2PSqUlx9Vfb0z3N1AGCQkxCkjI5Ux4Cd3oceb2cCWKsm7BGMZljEQGcI568gYiD42W6ri41vf4eRXgJKkjRs36pFHHlFBQfMfTH369NEf//hH9e/f35+XiVq8acyLD07/BSr0OJ2Gps5a5faydps1WU9NHBbQL1bGABgD0cfTAOX3WXj9+/fXO++8o23btik/P1+lpaVyOp3q2LGjBg4cqH79+vn7EgBMKJDH53hTE4jlWQCe8DtAufTu3Vu9e/cO1NMBMLlAhh5qAgEINJ/LGABAMAUy9FATCECgEaAARKRAhh5qAgEINAIUgIgUyNBDJXkAgUaAAhCRAh16qCTvGafT0Lbd5VqzZb+27S6X0+nXhdpAzPK7jAHc49JV8+Ly5cBorg6UP4UQQ1kTKNrGAIVGAy/axgBCWAcK7vGmMS8+OAMnWo9fiaYxEI5Co2YQTWMAJ4SsDhQABBvH5wRXIGtuAWbBHigAMDlvam4BOIEABQAmR6FRwHsEKAAwOQqNAt5jDxQA+CBaN7Y3x1Vzq7XDlik0CvyIAAUAXoq1y/1dNbfcXYVHoVGgMZbwAMALrsv9T52tKa+s0YuLC5RfWBqmnvmHQqOAd5iBAgAPxfrl/rnZnTUgs1PMLE0CwUSAijGxtC8DiDTeXO4frXWrIr3mFp9xiBQEqBgSa/sygEjD5f7hxWccIgl7oGJErO7LACIJl/uHD59xiDQEqBjg6b4MTlUH/OO63N8dLvcPPD7jEIkIUDGAYxiA0HBd7u8Ol/sHHp9xiEQEqBjAvgwgdLjcP/T4jEMkYhN5DGBfBhBaXO4fWnzGIRIRoGIAxzAAoRfpl/vHEj7jEIlYwosB7MsAEMv4jEMkIkDFCPZlAIhlfMYh0rCEF0PYlwEglvEZh0hCgIox7MsAEMv4jEOkIEABIcIZXgAQOwhQQAhwhhcAxBY2kQNBtm4bZ3gBQKwhQAFB5HAa+ufHhW7bcIYXAEQfAhQQRFu+P6wyzvACgJhDgAKCqOxItUftOMMLAKKL6TaRV1dX68svv9SmTZtUUFCggoIC2e12SdL69euVmpoa3g4iptjS2njUjjO8ACC6mC5A7dq1S5MnTw53NxAl/C090OesDrJZk90u43GGFwBEH9MFqISEBPXv31/9+vVTTk6OOnTooAkTJoS7W4hAgSg9EB9n0U1XZmvmOxtbbMMZXgAQfSyGYZj68p+9e/fq8ssvlxScJbzy8qOqr3cG9DkRfPmFJ0oPtMSTs7cSEuKUkZGq8vKjWrt5f5MwZrMmaxx1oGLayWOAzwFzYgxEH5stVfHxrW8RN90MFNAap9PQwuU73LZ5Y/kODcjs5PHMEWd4AUBsIUABp9hebG9S9PJUrtID3pzJxRleABA7KGMAnMLTkgKUHgAA82IGKsg8WUdFZOngYemBDmltlJDQ8t+v6++eMWBejAEwBmIXASrI0tJSwt0FeGlw+7bq8N4WHa5ouQhmx/QUDT7vdMV7sIeJMQDGABgDsSdqAtTDDz+st956y+vHDRo0SPPnzw9Cjzxz5MhxORxceRFtxv8sq9XSA0cqjrl9jvj4OKWlpTAGTIwxAMZA9ElLS4mtq/CsVqs6duzo9ePat28fhN54zuFwculqFBpwTkdNHp3TYumBAed09PjvlTEAxgAYA7GHOlDUgYIb/lQip/4LGANgDEQf6kABARCrpQf8PaIGQHjxHg4/AhRgMoE4ogY/4osMocZ7ODIQoAATaemImvLKGr24uMCjI2rwI77IEGq8hyOHKQPU6NGjVVJSIkk6eQvYZZdd1vDvAwcO1KxZs0LeNyBYgnFEjZl58kU2uG/XMPQMsYr3cGQxZWWviooK2e122e12VVRUNNzuus1ut6uqqiqMPQQCz5sjauCep19kTqepr9FBgPEejiymnIH69NNPw90FIOQ4oiZwPP0iK9xTrmEd2oWoV4h1vIcjiylnoAAzSk9NDmg7M/P4i6yqNsg9gZnwHo4sBCjAJLJ6pCvD6v6D1WY9cRUZ3PP4i6xdUpB7AjPhPRxZCFCAScTFWTR+RKbbNuNGZLL51AOefpFlnxF7NcQQPryHIwsBCjCR3OzOmjw6p8mXv82azOXPXuCLDOHCezhymP4ol2CjfL95RfIRDhR/DIzm6kC5zkvMze4c0WMAoRGsMcB7OHg4ygVAi2L1iJpQy83urAGZnfgiQ8jxHg4/AhQQJE6noU07D6l4X4WsKYl8scYovsgAcyJAAUGQX1iqN5bvUBlHfABATGITORBgriM+yk4ptOg64iO/sDRMPQMABAoBCgggjvgAAHMgQAEBxFlVAGAOBCgggDirCgDMgU3kQAAF46wq6r0AQOQhQAEB5Driw90ynjdnVTVXqJGr+QAg/FjCMzmn09C23eVas2W/tu0uZ3OznwJ5xIfrar5TwxhX8wFA+DEDZWLMbgSH66yqU+tAnXzER2s8vZpvQGYnlvMAIAwIUCblmt04lWt2g0Mp/ZOb3VkXnNtFJeXVPlUi9+ZqPn+qYLO/CgB8Q4AyIWY3QiMuzqJ+53TU6R1SvD5ENBRX8zEDGVqEVSC2EKBMKFSzG/BdMK7mOxkzkKG1blupFnxcSFgFYgibyE2IWkWRz3U1nzveXM13Mqqlh9aqjSWa+c5GLgYAYgwByoSCPbsB/wXyar5TUS09dJxOQ68s2eS2DWEViE4EKBMK5uwGAsd1Nd+pf1c2a7JfS2zMQIZO4Z5yHa6odtuGsApEJ/ZAmZBrdqO5PTAuvs5uILBysztrQGangG4+ZgYydOxVtZ61I6wCUYcAZVKu2Y1Tr8LyplYRQiMuzhLQzfyBrpaOlqW3S/KsHWEViDoEKBMLxuwGIh8zkKGTfUaGOrRv43YZj7AKRCf2QJmca3ZjSJ+u6n1mBl+aJhGs/VVoLC7Oogl5/dy2IawC0cliGAaXfwRReflRr4soIjYkJMQpIyM1oscAxR2DyzUG/r16V5M6UCyXm0M0fA6gMZstVfHxrc8vsYQHmFig91eheRf07qzzzupAWAViCAEKAEKAsArEFvZAAQAAeIkABQAA4CUCFAAAgJfYAxUDuJIKAIDQIkBFufzC0ibVxDOsyRrP5dEAAAQNS3hRLL+wVC8uLmhyJEd5ZY1eXFyg/MLSMPUMAIDYRoCKUk6noYXLd7ht88byHXI6qZMKAECgEaCi1PZiu9vDYCWprLJG24vtoekQAAAmwh6oKGU/6j48edsumrGJHgAQagSoKJWemtx6Iy/aRSs20QMAwoElvCiV1SNdGVb34chmPTEbE6vYRA8ACBcCVJSKi7No/IhMt23GjciM2aUsNtEDAMKJABXFcrM7a/LonCYzUTZrsiaPzonpJSw20QMAwok9UFEuN7uzBmR2Mt0majbRAwDCiQAVA+LiLOp9Zka4uxFSbKIHAIQTS3iISmyiBwCEEwEKUcnsm+gBAOFFgELUMvMmegBAeLEHClHNrJvoAQDhRYBC1DPjJnoAQHixhAcAAOAl081AFRUV6ZNPPtHatWtVWFio8vJyJScnq1evXrriiit00003KTU1NdzdBAAAEcxiGIZpzrpwOBzq06dPw587deqkLl266NChQ9q/f78k6YwzztDcuXPVvXv3gLxmeflR1dc7A/JciC4JCXHKyEhlDJgYYwCMgehjs6UqPr71BTpTzUAZhqF27dpp7Nixuvbaa3X22Wc33Ldhwwbdd9992rNnj+699169/fbbYewpAACIZKaagTIMQxUVFUpPT2/2/vXr12vcuHGSpCVLlujcc8/1+zX51WFe/PIEYwCMgejj6QyUqTaRWyyWFsOTJA0cOFBWq1WStGvXrhD1CgAARBtTBajWOBwO1dfXS5LatGkT5t4AAIBIRYA6yYoVK3T8+HElJCTo/PPPD3d3AABAhDLVJnJ3qqqq9OSTT0qSxowZI5vNFpDn9WQdFbHJ9XfPGAgOp9NQ4Z5y2atqld4uSdlnZERcBXrGABgDsctUm8hb4nA4NHHiRH3++efq3r27li5d2rAXCkDkWbWxRK8s2aTDFdUNt3Vo30YT8vppWP9uYewZALOImgD18MMP66233vL6cYMGDdL8+fNbvN8wDD344IN699131b59ey1YsEBZWVn+dLWRI0eOy+Hgygszio+PU1paCmMgwNZtK9XMdza2eP891/XXBb0j4yBpxgAYA9EnLS0ltupAWa1WdezY0evHtW/f3u39f/7zn/Xuu+8qNTVVr776akDDkyQ5HE4uXTU5xkDgOJ2GFnxc6LbNPz8u1HlndYio5TzGABgDsSdqAtTUqVM1derUgD7nk08+qQULFiglJUWvvPKK+vfvH9DnDzSn09D2YrvsR2uUnpqsrB7pEfUlAQTb9mK7yitr3LYpq6zR9mI7B0wDCKqoCVCBNmPGDM2ZM0fJycmaNWuWfvKTn4S7S27lF5Zq4fIdjb48MqzJGj8iU7nZkbFcAQSb/aj78ORtOwDwlSkvC5g9e7Zmz56txMREzZw5U0OHDg13l9zKLyzVi4sLmvzyLq+s0YuLC5RfWBqmngGhlZ6aHNB2AOAr0wWoefPmacaMGUpISNCMGTN0ySWXhLtLbjmdhhYu3+G2zRvLd8jpjIprAQC/ZPVIV4bVfTiyWU8sbwNAMJlqCe/AgQN67LHHJEmpqamaM2eO5syZ02zbMWPG6Lrrrgtl95rFng/gR3FxFo0fkakXFxe02GbciEz2BgIIOlMFqLq6OrmqNlRUVGj9+vUtth02bFiouuUWez6AxnKzO2vy6JwmewJt1mSNY08ggBAxVYA6/fTTVVjo/hLoSMOeD6Cp3OzOGpDZiatSAYSNqQJUNHLt+XC3jMeeD5hRXJyFZWsAYWO6TeTRxrXnwx32fAAAEFoEqCjg2vNx6tVHNmuyJo/OYc8HAAAhxhJelGDPBwAAkYMAFUXY8wEAQGRgCQ8AAMBLBCgAAAAvEaAAAAC8RIACAADwEgEKAADASwQoAAAALxGgAAAAvESAAgAA8BIBCgAAwEsWwzCMcHciljkcznB3AWEUHx/HGDA5xgAYA9ElLs4ii6X1Y9IIUAAAAF5iCQ8AAMBLBCgAAAAvEaAAAAC8RIACAADwEgEKAADASwQoAAAALxGgAAAAvESAAgAA8BIBCgAAwEsEKAAAAC8RoAAAALxEgAIAAPASAQoAAMBLBCgAAAAvJYS7A2ZRVFSkTz75RGvXrlVhYaHKy8uVnJysXr166YorrtBNN92k1NTUcHcTQVZdXa0vv/xSmzZtUkFBgQoKCmS32yVJ69evZwzEiDVr1ui1117Tt99+q2PHjqlbt24aOXKkJkyYoLZt24a7ewiygwcPatWqVQ3v861bt6q6ulp9+/bVokWLwt09BIjFMAwj3J2IdQ6HQ3369Gn4c6dOndSlSxcdOnRI+/fvlySdccYZmjt3rrp37x6ubiIEtm7dqry8vGbvI0DFhvnz5+svf/mLDMNQ165dZbPZtHPnTtXW1urss8/WwoULlZ6eHu5uIojmzp2rxx9/vMntBKjYwgxUCBiGoXbt2mns2LG69tprdfbZZzfct2HDBt13333as2eP7r33Xr399tth7CmCLSEhQf3791e/fv2Uk5OjDh06aMKECeHuFgKkoKBAjz32mCTpkUce0Q033CCLxaIDBw5o4sSJ2rx5s6ZPn66ZM2eGuacIpnbt2mnYsGHKyclRTk6OioqK9Oyzz4a7WwgwZqBCwDAMVVRUtPirc/369Ro3bpwkacmSJTr33HND2DuE0969e3X55ZdLYgYqFkyaNEkrVqxQXl6ennzyyUb3FRUV6aqrrpLT6dTSpUvVu3fvMPUSobZo0SI98MADzEDFGDaRh4DFYnE7ZT9w4EBZrVZJ0q5du0LUKwCBdPToUX3xxReSpBtuuKHJ/T179tSQIUMkSR999FFI+wYg8AhQEcDhcKi+vl6S1KZNmzD3BoAvtm7dqtraWiUlJal///7NtsnNzZUkffvtt6HsGoAgIEBFgBUrVuj48eNKSEjQ+eefH+7uAPCBa/a4W7duSkxMbLbNGWec0agtgOhFgAqzqqqqhr0SY8aMkc1mC3OPAPiioqJCktS+ffsW27juc7UFEL0IUGHkcDg0ZcoU7d27V927d9fUqVPD3SUAPqqpqZGkFmefJCkpKalRWwDRizIGrXj44Yf11ltvef24QYMGaf78+S3ebxiGpk+frs8//1zt27fX7NmzGzaSI/IEaxwgdiQnJ0uS6urqWmxTW1vbqC2A6EWAaoXValXHjh29fpy7aXxJ+vOf/6x3331XqampevXVV5WVleVrFxECwRoHiB2eLM95sswHIDoQoFoxderUgC+tPfnkk1qwYIFSUlL0yiuvtHjFDiJHMMYBYkvPnj0lSSUlJaqrq2t2KW/Pnj2N2gKIXuyBCrEZM2Zozpw5Sk5O1qxZs/STn/wk3F0CEAB9+vRRYmKiamtrtXHjxmbb5OfnSxJX2wIxgAAVQrNnz9bs2bOVmJiomTNnaujQoeHuEoAASU1N1fDhwyWp2SOZioqKtGbNGknSyJEjQ9o3AIFHgAqRefPmacaMGUpISNCMGTN0ySWXhLtLAAJs0qRJslgsWrp0qd566y25TsoqLS3VlClT5HQ6NWLECI5xAWIAZ+GFwIEDB3TJJZfIMAy1b9++0WHCpxozZoyuu+66EPYOoTZ69GiVlJRI+vGcREmNjvsZOHCgZs2aFY7uwU9z587VE088IcMwdNpppykjI0M7d+5UbW2tevXqpYULF1LvLcbt27dPeXl5DX+ura3VsWPHlJCQoHbt2jXcfuedd+quu+4KQw8RCGwiD4G6urqGX6IVFRVav359i22HDRsWqm4hTCoqKmS325vcfvJtVVVVoesQAuq2225Tdna25syZo40bN+rw4cPq1q2bRo4cqQkTJnBgtAk4HI5m3+P19fWNbq+urg5dpxBwzEABAAB4iT1QAAAAXiJAAQAAeIkABQAA4CUCFAAAgJcIUAAAAF4iQAEAAHiJAAUAAOAlAhQAAICXCFAAAABeIkABAAB4ibPwACDKHThwQMuXL9fatWu1detWHThwQJLUsWNHnX/++br++us1dOjQMPcSiC2chQcAUWzfvn366U9/qpM/ylNSUmQYRqPDaseMGaNHH31U8fHx4egmEHOYgQKAKOZwOGQYhoYOHaq8vDwNHTpUXbp0kdPp1Pfff69nn31WK1as0LvvvqvOnTvr//2//xfuLgMxgRkoAIhilZWV2rNnj/r27dvs/YZh6K677tIXX3yhtm3bas2aNUpOTg5xL4HYwyZyAGF18OBBPf300/rlL3+p3Nxc9evXT5dffrkefPBB7dy5s0n7V155RdnZ2crJydHGjRubfc7PP/9cvXv3VnZ2tt57770m969du1a/+c1vdNFFFyknJ0eDBw/WrbfeqnfffVcOh8Ntf9etW6f//d//1eDBg9W/f39deeWVmjFjho4ePapFixYpOztbl112mVf/Dz744ANlZ2crOztbn3zySbNttm7dqn79+ik7O1uzZ89uuN1qtbYYniTJYrFozJgxkqRjx47pu+++86pvAJpHgAIQNitXrtQVV1yhv//97yosLFRNTY0SEhK0d+9evfPOO8rLy9OSJUsaPeauu+7SsGHDVFdXpylTpqiqqqrR/aWlpfrd734nwzCUl5enq6++utH9jz/+uG655RZ9/PHHOnjwoFJSUlRZWak1a9bo97//vW6//fYmz+kyf/583XzzzVq5cqXsdruSkpL0ww8/aPbs2br++ut15MgRn/4/jBo1qiHkPPTQQ9q3b1+j+48dO6Z7771XtbW1GjRokCZMmODV858849RaQATgGQIUgLDYuHGj7rnnHh07dkw33nijPvjgA3377bf65ptvtHLlSo0fP151dXV68MEHtWnTpobHWSwWPfXUU+rQoYOKi4v1hz/8oeE+wzB0//33q6ysTGeeeaYefvjhRq+5YMECzZ07V5J044036osvvtC6dev03//+Vw888IASEhK0Zs0aTZ8+vUl/169fr8cee0yGYejCCy/URx99pP/+97/asGGDnn/+eR0+fFgvvfSSz/8/pk+frrPOOkt2u11Tp05tFHQeffRR7dq1S+np6Xr66acVF+fdR/fXX38tSUpMTFSvXr187iOAkxgAEAZjxowxsrKyjOeee67FNo8++qiRlZVlTJw4scl9n332mZGdnW1kZWUZixYtMgzDMF5++WUjKyvL6Nu3r7Fx48ZG7Y8fP24MGjTIyMrKMqZMmdLs682bN8/IysoysrKymjz+1ltvNbKysoxRo0YZNTU1TR67evXqhsf+9Kc/bfW/vzlbt241cnJyjKysLGPmzJmGYRjGsmXLGp73k08+8fo59+zZY5x33nlGVlaWMXXqVJ/6BaApZqAAhNy2bdu0adMmJSYm6vbbb2+xXV5eniRp9erVTZaeLrnkEt12222SpEceeURLly7V3/72N0nSvffeq379+jVq/9VXX8lut0uSfv3rXzf7euPHj1enTp0kScuWLWu43W63a82aNZKkO+64Q0lJSU0eO2TIEP3kJz9p8b/FE71799a0adMkSS+99JKWLFnSMIt20003acSIEV49X3V1tf7v//5Px48fV3p6uu677z6/+gfgR5QxABBy+fn5kiSn06mRI0e22M4Vmo4dOya73a4OHTo0un/KlCn6+uuvtXnz5obgMXz48GZDWUFBgSTptNNOa3EZKz4+XkOGDNF7773X0F46sYHb+P8vWL7gggta7O+gQYP03//+t8X7PXHzzTfrq6++0sqVK3X//fdLkrKyshr+3VP19fW67777tHnzZiUmJuqZZ55Rly5d/OobgB8xAwUg5EpLSyWdCEiHDh1q8Z/y8vKGxxw/frzJ8yQlJemJJ55o+LPVatUTTzwhi8XSpO3hw4clqdUQ0bVr10btJamsrKzh3909PlAB5bHHHmvY+B0fH69nn33Wq9IDDodDU6dO1fLly5WQkKCnn35aw4cPD0jfAJzADBSAkHPNLJ111ln68MMP/Xqut99+u+Hfq6qqtHXr1oZluOY0F65aa2d4WC7P03atWbp0qWpqaiSd+H+Vn5+vzMxMjx7rCk8ffPCB4uPj9de//tXtLB8A3zADBSDkOnbsKEkqLi7WsWPHfH6elStXav78+ZKk7OxsGYah3/3udzp06FCTtq7lv1NLBJxq//79kiSbzdbksdKPs2fNcXefpzZv3qxnnnlG0on/JulE6YXmamKdyuFw6Le//a2WLVvWEJ5GjRrld58ANEWAAhByAwcOlCTV1dW1WDiyNaWlpXrggQckSddee63++c9/qnv37jp8+LDuv//+JrNBOTk5kk4EpF27djX7nA6HQ2vXrpWkRpvQzz333IYZKVdJgOa4u88Tx44d05QpU1RXV6chQ4bonXfeUf/+/VVdXa0pU6aotra2xce6wtPJM08///nP/eoPgJYRoACEXL9+/dSnTx9J0nPPPddoj1FzXFfPuTidTk2bNk3l5eXq2bOnpk+fLqvVqmeeeUYJCQn68ssv9dprrzV6zIUXXqj09HRJ0gsvvNDs67z55psNs0gnh4/09HQNHjxYkvTaa681G2Rc9aT88eijj6qoqEjp6el66qmnlJSUpGeeeUapqakqLCzUk08+2ezjHA6H7rvvPn3wwQcNe54IT0BwEaAAhJzFYtGf/vQnJSUlqaSkRNdff70++uijRhvFDxw4oKVLl+p//ud/9PTTTzd6/KuvvqrVq1c3XF3Wtm1bSdKAAQM0efJkSdKzzz6rzZs3NzymTZs2uueeeyRJ77//vh5++OGGpb7jx49r/vz5evzxxyWdqAzumrFyueeee2SxWLR9+3ZNnDhRRUVFkk5c7fbvf/9b99xzj9q3b+/z/5Nly5Zp0aJFkk5sIndtSD/jjDMaioUuWLBAK1eubPQ4h8OhadOm6cMPP2wITyzbAcHHYcIAwuarr77SlClTGmaY4uPjZbVaVVNT0yhMXX/99frzn/8s6UQFc1eV8mnTpumOO+5o9JxOp1O33nqrvv76a/Xs2VOLFy9uCFjSif1ErmrkFotFaWlpOnr0qOrr6yVJgwcP1ksvvaR27do16e/cuXMbQpYkpaWlqbq6WrW1tcrKytKYMWP0+OOPq1evXvroo488/v+wd+9e5eXlqbKyUjfddFOTCuqSNG3aNC1dulQZGRlaunRpQ8Bat26dfvWrX0k6UWm8tRD34IMPErCAACBAAQirI0eO6M0339Rnn32m7777TpWVlUpOTla3bt10/vnn6/LLL9eFF16o5ORkVVVVKS8vT8XFxbrwwgv1j3/8o9mr6vbv369rrrlGdrtd1157baPQI0lr1qzRwoULtX79etntdrVt21a9e/fWNddco7y8PMXHx7fY36+//lqvvvqqNmzYoOrqanXr1k0jR47UhAkT9Pbbb+vxxx/Xeeed1+jqQHfq6+t10003acOGDcrKytI777zTbMmCqqoqXXvttdq9e7eGDh2qOXPmKC4uTmvXrtUtt9zi0WtJJwLktdde63F7AM0jQAFAgNx33316//33NWbMGD322GPh7g6AIGIPFAAEwK5duxquKLzooovC3BsAwUaAAgAPPf/881qwYIFKSkrkdDolnSg98MEHH+iWW25RTU2NzjrrLK/PrAMQfVjCAwAPTZo0SStWrJB0YsN2amqqjhw50hCmunTpoldffVVZWVnh7CaAEOAoFwDw0G233abOnTvrm2++0cGDB1VRUaHU1FT17NlTl156qX71q1811JoCENuYgQIAAPASe6AAAAC8RIACAADwEgEKAADASwQoAAAALxGgAAAAvESAAgAA8BIBCgAAwEsEKAAAAC8RoAAAALz0/wFluKBfoH8kOgAAAABJRU5ErkJggg==",
+ "text/plain": [
+ "<Figure size 640x480 with 1 Axes>"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plt.figure()\n",
+ "plt.scatter(x1, y)\n",
+ "plt.xlabel('exog x1')\n",
+ "plt.ylabel('endog y')\n",
+ "\n",
+ "plt.figure()\n",
+ "plt.scatter(x2, y)\n",
+ "plt.xlabel('exog x2')\n",
+ "plt.ylabel('endog y')\n",
+ "\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "d30de122-7250-4c31-9e6d-f2e22370b242",
+ "metadata": {
+ "tags": []
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ " OLS Regression Results \n",
+ "==============================================================================\n",
+ "Dep. Variable: y R-squared: 0.484\n",
+ "Model: OLS Adj. R-squared: 0.462\n",
+ "Method: Least Squares F-statistic: 22.01\n",
+ "Date: Thu, 06 Jul 2023 Prob (F-statistic): 1.80e-07\n",
+ "Time: 00:36:28 Log-Likelihood: -62.306\n",
+ "No. Observations: 50 AIC: 130.6\n",
+ "Df Residuals: 47 BIC: 136.3\n",
+ "Df Model: 2 \n",
+ "Covariance Type: nonrobust \n",
+ "==============================================================================\n",
+ " coef std err t P>|t| [0.025 0.975]\n",
+ "------------------------------------------------------------------------------\n",
+ "Intercept -0.0413 0.125 -0.332 0.741 -0.292 0.209\n",
+ "x1 -0.4097 0.118 -3.468 0.001 -0.647 -0.172\n",
+ "x2 0.5906 0.136 4.335 0.000 0.316 0.865\n",
+ "==============================================================================\n",
+ "Omnibus: 1.032 Durbin-Watson: 2.003\n",
+ "Prob(Omnibus): 0.597 Jarque-Bera (JB): 0.986\n",
+ "Skew: 0.167 Prob(JB): 0.611\n",
+ "Kurtosis: 2.398 Cond. No. 1.47\n",
+ "==============================================================================\n",
+ "\n",
+ "Notes:\n",
+ "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
+ ]
+ }
+ ],
+ "source": [
+ "# define linear model\n",
+ "lm = smf.ols('y ~ x1 + x2', df)\n",
+ "\n",
+ "# fit model to data\n",
+ "lmf = lm.fit()\n",
+ "\n",
+ "# summary\n",
+ "print(lmf.summary())"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "776bee7d-176d-443d-90ff-256ab36bb5e2",
+ "metadata": {},
+ "source": [
+ "The construction of the model relies on the design matrix created by `patsy.dmatrix`."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "9237dd64-3e28-4a6e-8a26-fb8fa63d9842",
+ "metadata": {
+ "tags": []
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "<div>\n",
+ "<style scoped>\n",
+ " .dataframe tbody tr th:only-of-type {\n",
+ " vertical-align: middle;\n",
+ " }\n",
+ "\n",
+ " .dataframe tbody tr th {\n",
+ " vertical-align: top;\n",
+ " }\n",
+ "\n",
+ " .dataframe thead th {\n",
+ " text-align: right;\n",
+ " }\n",
+ "</style>\n",
+ "<table border=\"1\" class=\"dataframe\">\n",
+ " <thead>\n",
+ " <tr style=\"text-align: right;\">\n",
+ " <th></th>\n",
+ " <th>Intercept</th>\n",
+ " <th>x1</th>\n",
+ " <th>x2</th>\n",
+ " </tr>\n",
+ " </thead>\n",
+ " <tbody>\n",
+ " <tr>\n",
+ " <th>0</th>\n",
+ " <td>1.0</td>\n",
+ " <td>-0.598903</td>\n",
+ " <td>-1.107940</td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <th>1</th>\n",
+ " <td>1.0</td>\n",
+ " <td>1.484895</td>\n",
+ " <td>-0.308404</td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <th>2</th>\n",
+ " <td>1.0</td>\n",
+ " <td>-1.279121</td>\n",
+ " <td>-0.071924</td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <th>3</th>\n",
+ " <td>1.0</td>\n",
+ " <td>0.796557</td>\n",
+ " <td>-1.757192</td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <th>4</th>\n",
+ " <td>1.0</td>\n",
+ " <td>0.772527</td>\n",
+ " <td>0.049753</td>\n",
+ " </tr>\n",
+ " </tbody>\n",
+ "</table>\n",
+ "</div>"
+ ],
+ "text/plain": [
+ " Intercept x1 x2\n",
+ "0 1.0 -0.598903 -1.107940\n",
+ "1 1.0 1.484895 -0.308404\n",
+ "2 1.0 -1.279121 -0.071924\n",
+ "3 1.0 0.796557 -1.757192\n",
+ "4 1.0 0.772527 0.049753"
+ ]
+ },
+ "execution_count": 5,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "y_aff, X_aff = dmatrices('y ~ x1 + x2', data=df, return_type='dataframe')\n",
+ "X_aff.head()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "b9596e22-c8b4-43ea-88ca-894c1cb4bbff",
+ "metadata": {},
+ "source": [
+ "## Interactions between predictors\n",
+ "\n",
+ "Now we change the model to incorporate a non-linear interaction between $x1$ and $x2$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "id": "0ae0753d-9007-489a-a102-5a1a1f1b97c9",
+ "metadata": {
+ "tags": []
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "<div>\n",
+ "<style scoped>\n",
+ " .dataframe tbody tr th:only-of-type {\n",
+ " vertical-align: middle;\n",
+ " }\n",
+ "\n",
+ " .dataframe tbody tr th {\n",
+ " vertical-align: top;\n",
+ " }\n",
+ "\n",
+ " .dataframe thead th {\n",
+ " text-align: right;\n",
+ " }\n",
+ "</style>\n",
+ "<table border=\"1\" class=\"dataframe\">\n",
+ " <thead>\n",
+ " <tr style=\"text-align: right;\">\n",
+ " <th></th>\n",
+ " <th>x1</th>\n",
+ " <th>x2</th>\n",
+ " <th>y</th>\n",
+ " <th>y2</th>\n",
+ " </tr>\n",
+ " </thead>\n",
+ " <tbody>\n",
+ " <tr>\n",
+ " <th>0</th>\n",
+ " <td>-0.598903</td>\n",
+ " <td>-1.107940</td>\n",
+ " <td>-0.983083</td>\n",
+ " <td>-0.717664</td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <th>1</th>\n",
+ " <td>1.484895</td>\n",
+ " <td>-0.308404</td>\n",
+ " <td>-2.069647</td>\n",
+ " <td>-2.252826</td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <th>2</th>\n",
+ " <td>-1.279121</td>\n",
+ " <td>-0.071924</td>\n",
+ " <td>0.927597</td>\n",
+ " <td>0.964397</td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <th>3</th>\n",
+ " <td>0.796557</td>\n",
+ " <td>-1.757192</td>\n",
+ " <td>-0.702349</td>\n",
+ " <td>-1.262231</td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <th>4</th>\n",
+ " <td>0.772527</td>\n",
+ " <td>0.049753</td>\n",
+ " <td>0.529516</td>\n",
+ " <td>0.544891</td>\n",
+ " </tr>\n",
+ " </tbody>\n",
+ "</table>\n",
+ "</div>"
+ ],
+ "text/plain": [
+ " x1 x2 y y2\n",
+ "0 -0.598903 -1.107940 -0.983083 -0.717664\n",
+ "1 1.484895 -0.308404 -2.069647 -2.252826\n",
+ "2 -1.279121 -0.071924 0.927597 0.964397\n",
+ "3 0.796557 -1.757192 -0.702349 -1.262231\n",
+ "4 0.772527 0.049753 0.529516 0.544891"
+ ]
+ },
+ "execution_count": 6,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "# coefficient for the interaction, which is simply a second-order polynomial in x1, x2\n",
+ "a12 = 0.4\n",
+ "y2 = a1 * x1 + a2 * x2 + a12 * x1 * x2 + e\n",
+ "\n",
+ "# build dataframe\n",
+ "df['y2'] = y2\n",
+ "df.head()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "id": "fc12371e-4da7-4497-aebb-a55a7aae0daf",
+ "metadata": {
+ "tags": []
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "<div>\n",
+ "<style scoped>\n",
+ " .dataframe tbody tr th:only-of-type {\n",
+ " vertical-align: middle;\n",
+ " }\n",
+ "\n",
+ " .dataframe tbody tr th {\n",
+ " vertical-align: top;\n",
+ " }\n",
+ "\n",
+ " .dataframe thead th {\n",
+ " text-align: right;\n",
+ " }\n",
+ "</style>\n",
+ "<table border=\"1\" class=\"dataframe\">\n",
+ " <thead>\n",
+ " <tr style=\"text-align: right;\">\n",
+ " <th></th>\n",
+ " <th>Intercept</th>\n",
+ " <th>x1</th>\n",
+ " <th>x2</th>\n",
+ " <th>x1:x2</th>\n",
+ " </tr>\n",
+ " </thead>\n",
+ " <tbody>\n",
+ " <tr>\n",
+ " <th>0</th>\n",
+ " <td>1.0</td>\n",
+ " <td>-0.598903</td>\n",
+ " <td>-1.107940</td>\n",
+ " <td>0.663548</td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <th>1</th>\n",
+ " <td>1.0</td>\n",
+ " <td>1.484895</td>\n",
+ " <td>-0.308404</td>\n",
+ " <td>-0.457947</td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <th>2</th>\n",
+ " <td>1.0</td>\n",
+ " <td>-1.279121</td>\n",
+ " <td>-0.071924</td>\n",
+ " <td>0.092000</td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <th>3</th>\n",
+ " <td>1.0</td>\n",
+ " <td>0.796557</td>\n",
+ " <td>-1.757192</td>\n",
+ " <td>-1.399704</td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <th>4</th>\n",
+ " <td>1.0</td>\n",
+ " <td>0.772527</td>\n",
+ " <td>0.049753</td>\n",
+ " <td>0.038436</td>\n",
+ " </tr>\n",
+ " </tbody>\n",
+ "</table>\n",
+ "</div>"
+ ],
+ "text/plain": [
+ " Intercept x1 x2 x1:x2\n",
+ "0 1.0 -0.598903 -1.107940 0.663548\n",
+ "1 1.0 1.484895 -0.308404 -0.457947\n",
+ "2 1.0 -1.279121 -0.071924 0.092000\n",
+ "3 1.0 0.796557 -1.757192 -1.399704\n",
+ "4 1.0 0.772527 0.049753 0.038436"
+ ]
+ },
+ "execution_count": 7,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "y_aff, X_aff = dmatrices('y2 ~ x1 * x2', data=df, return_type='dataframe')\n",
+ "X_aff.head()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "id": "8c0ae1bb-4493-4a1f-bb2a-d467d216fc62",
+ "metadata": {
+ "tags": []
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ " OLS Regression Results \n",
+ "==============================================================================\n",
+ "Dep. Variable: y2 R-squared: 0.477\n",
+ "Model: OLS Adj. R-squared: 0.443\n",
+ "Method: Least Squares F-statistic: 14.01\n",
+ "Date: Thu, 06 Jul 2023 Prob (F-statistic): 1.28e-06\n",
+ "Time: 00:36:41 Log-Likelihood: -62.280\n",
+ "No. Observations: 50 AIC: 132.6\n",
+ "Df Residuals: 46 BIC: 140.2\n",
+ "Df Model: 3 \n",
+ "Covariance Type: nonrobust \n",
+ "==============================================================================\n",
+ " coef std err t P>|t| [0.025 0.975]\n",
+ "------------------------------------------------------------------------------\n",
+ "Intercept -0.0335 0.131 -0.256 0.799 -0.297 0.230\n",
+ "x1 -0.4150 0.122 -3.405 0.001 -0.660 -0.170\n",
+ "x2 0.5960 0.140 4.259 0.000 0.314 0.878\n",
+ "x1:x2 0.4273 0.126 3.389 0.001 0.174 0.681\n",
+ "==============================================================================\n",
+ "Omnibus: 0.909 Durbin-Watson: 1.989\n",
+ "Prob(Omnibus): 0.635 Jarque-Bera (JB): 0.917\n",
+ "Skew: 0.164 Prob(JB): 0.632\n",
+ "Kurtosis: 2.423 Cond. No. 1.71\n",
+ "==============================================================================\n",
+ "\n",
+ "Notes:\n",
+ "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
+ ]
+ }
+ ],
+ "source": [
+ "# define linear model\n",
+ "lm2 = smf.ols('y2 ~ x1 * x2', df)\n",
+ "\n",
+ "# fit model to data\n",
+ "lmf2 = lm2.fit()\n",
+ "\n",
+ "# summary\n",
+ "print(lmf2.summary())"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "3761c74f-1dc0-47a0-aeb3-ccb928875dad",
+ "metadata": {},
+ "source": [
+ "This model accurately captures the coefficients used to generate $y2$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "id": "dfc36822-9af7-4eab-9e25-baaa0b9b6eb9",
+ "metadata": {
+ "tags": []
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ " OLS Regression Results \n",
+ "==============================================================================\n",
+ "Dep. Variable: y2 R-squared: 0.347\n",
+ "Model: OLS Adj. R-squared: 0.319\n",
+ "Method: Least Squares F-statistic: 12.48\n",
+ "Date: Thu, 06 Jul 2023 Prob (F-statistic): 4.48e-05\n",
+ "Time: 00:36:52 Log-Likelihood: -67.852\n",
+ "No. Observations: 50 AIC: 141.7\n",
+ "Df Residuals: 47 BIC: 147.4\n",
+ "Df Model: 2 \n",
+ "Covariance Type: nonrobust \n",
+ "==============================================================================\n",
+ " coef std err t P>|t| [0.025 0.975]\n",
+ "------------------------------------------------------------------------------\n",
+ "Intercept -0.1561 0.139 -1.121 0.268 -0.436 0.124\n",
+ "x1 -0.3311 0.132 -2.509 0.016 -0.597 -0.066\n",
+ "x2 0.5106 0.152 3.354 0.002 0.204 0.817\n",
+ "==============================================================================\n",
+ "Omnibus: 0.529 Durbin-Watson: 2.231\n",
+ "Prob(Omnibus): 0.768 Jarque-Bera (JB): 0.634\n",
+ "Skew: -0.016 Prob(JB): 0.728\n",
+ "Kurtosis: 2.449 Cond. No. 1.47\n",
+ "==============================================================================\n",
+ "\n",
+ "Notes:\n",
+ "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
+ ]
+ }
+ ],
+ "source": [
+ "# model without interaction\n",
+ "lm = smf.ols('y2 ~ x1 + x2', df)\n",
+ "\n",
+ "# fit model to data\n",
+ "lmf = lm.fit()\n",
+ "\n",
+ "# summary\n",
+ "print(lmf.summary())"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "edf8ac33-7145-4bed-b9d9-c0362104debc",
+ "metadata": {},
+ "source": [
+ "Note that correctly estimating the interaction requires including the linear dependencies."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "id": "db858a08-c017-4c7b-9f2a-24d9fbe5dfbd",
+ "metadata": {
+ "tags": []
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ " OLS Regression Results \n",
+ "==============================================================================\n",
+ "Dep. Variable: y2 R-squared: 0.026\n",
+ "Model: OLS Adj. R-squared: 0.005\n",
+ "Method: Least Squares F-statistic: 1.257\n",
+ "Date: Thu, 06 Jul 2023 Prob (F-statistic): 0.268\n",
+ "Time: 00:37:28 Log-Likelihood: -77.858\n",
+ "No. Observations: 50 AIC: 159.7\n",
+ "Df Residuals: 48 BIC: 163.5\n",
+ "Df Model: 1 \n",
+ "Covariance Type: nonrobust \n",
+ "==============================================================================\n",
+ " coef std err t P>|t| [0.025 0.975]\n",
+ "------------------------------------------------------------------------------\n",
+ "Intercept -0.0803 0.173 -0.465 0.644 -0.428 0.267\n",
+ "x1:x2 0.1794 0.160 1.121 0.268 -0.142 0.501\n",
+ "==============================================================================\n",
+ "Omnibus: 0.414 Durbin-Watson: 2.034\n",
+ "Prob(Omnibus): 0.813 Jarque-Bera (JB): 0.567\n",
+ "Skew: 0.064 Prob(JB): 0.753\n",
+ "Kurtosis: 2.494 Cond. No. 1.35\n",
+ "==============================================================================\n",
+ "\n",
+ "Notes:\n",
+ "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
+ ]
+ }
+ ],
+ "source": [
+ "# model with only the interaction\n",
+ "lm3 = smf.ols('y2 ~ x1 : x2', df)\n",
+ "\n",
+ "# fit model to data\n",
+ "lmf3 = lm3.fit()\n",
+ "\n",
+ "# summary\n",
+ "print(lmf3.summary())"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "793b95f7-1df7-4b58-858d-309f2ab93717",
+ "metadata": {},
+ "source": [
+ "### EXERCISE\n",
+ "- Redo the fit of $y$ by the model with interactions (`lm2`)\n",
+ "- with the below data, check relationships and interactions between the cognitive scores.\n",
+ "- Adapt the synthetic model to the data and compare the perf with the model fitted on data.\n",
+ "\n",
+ "Dataset taken from R dataset repository ([link](https://github.com/vincentarelbundock/Rdatasets/)), use e.g.\n",
+ "`df = sm.datasets.get_rdataset(\"NeuroCog\", \"heplots\").data`"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "id": "e1625945-157d-4716-a6fd-5f287a4b17d4",
+ "metadata": {
+ "tags": []
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "<div>\n",
+ "<style scoped>\n",
+ " .dataframe tbody tr th:only-of-type {\n",
+ " vertical-align: middle;\n",
+ " }\n",
+ "\n",
+ " .dataframe tbody tr th {\n",
+ " vertical-align: top;\n",
+ " }\n",
+ "\n",
+ " .dataframe thead th {\n",
+ " text-align: right;\n",
+ " }\n",
+ "</style>\n",
+ "<table border=\"1\" class=\"dataframe\">\n",
+ " <thead>\n",
+ " <tr style=\"text-align: right;\">\n",
+ " <th></th>\n",
+ " <th>Unnamed: 0</th>\n",
+ " <th>Dx</th>\n",
+ " <th>Speed</th>\n",
+ " <th>Attention</th>\n",
+ " <th>Memory</th>\n",
+ " <th>Verbal</th>\n",
+ " <th>Visual</th>\n",
+ " <th>ProbSolv</th>\n",
+ " <th>SocialCog</th>\n",
+ " <th>Age</th>\n",
+ " <th>Sex</th>\n",
+ " </tr>\n",
+ " </thead>\n",
+ " <tbody>\n",
+ " <tr>\n",
+ " <th>0</th>\n",
+ " <td>14</td>\n",
+ " <td>Schizophrenia</td>\n",
+ " <td>19</td>\n",
+ " <td>9</td>\n",
+ " <td>19</td>\n",
+ " <td>33</td>\n",
+ " <td>24</td>\n",
+ " <td>39</td>\n",
+ " <td>28</td>\n",
+ " <td>44</td>\n",
+ " <td>Female</td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <th>1</th>\n",
+ " <td>15</td>\n",
+ " <td>Schizophrenia</td>\n",
+ " <td>8</td>\n",
+ " <td>25</td>\n",
+ " <td>15</td>\n",
+ " <td>28</td>\n",
+ " <td>24</td>\n",
+ " <td>40</td>\n",
+ " <td>37</td>\n",
+ " <td>26</td>\n",
+ " <td>Male</td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <th>2</th>\n",
+ " <td>16</td>\n",
+ " <td>Schizophrenia</td>\n",
+ " <td>14</td>\n",
+ " <td>23</td>\n",
+ " <td>15</td>\n",
+ " <td>20</td>\n",
+ " <td>13</td>\n",
+ " <td>32</td>\n",
+ " <td>24</td>\n",
+ " <td>55</td>\n",
+ " <td>Female</td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <th>3</th>\n",
+ " <td>17</td>\n",
+ " <td>Schizophrenia</td>\n",
+ " <td>7</td>\n",
+ " <td>18</td>\n",
+ " <td>14</td>\n",
+ " <td>34</td>\n",
+ " <td>16</td>\n",
+ " <td>31</td>\n",
+ " <td>36</td>\n",
+ " <td>53</td>\n",
+ " <td>Male</td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <th>4</th>\n",
+ " <td>18</td>\n",
+ " <td>Schizophrenia</td>\n",
+ " <td>21</td>\n",
+ " <td>9</td>\n",
+ " <td>35</td>\n",
+ " <td>28</td>\n",
+ " <td>29</td>\n",
+ " <td>45</td>\n",
+ " <td>28</td>\n",
+ " <td>51</td>\n",
+ " <td>Male</td>\n",
+ " </tr>\n",
+ " </tbody>\n",
+ "</table>\n",
+ "</div>"
+ ],
+ "text/plain": [
+ " Unnamed: 0 Dx Speed Attention Memory Verbal Visual \n",
+ "0 14 Schizophrenia 19 9 19 33 24 \\\n",
+ "1 15 Schizophrenia 8 25 15 28 24 \n",
+ "2 16 Schizophrenia 14 23 15 20 13 \n",
+ "3 17 Schizophrenia 7 18 14 34 16 \n",
+ "4 18 Schizophrenia 21 9 35 28 29 \n",
+ "\n",
+ " ProbSolv SocialCog Age Sex \n",
+ "0 39 28 44 Female \n",
+ "1 40 37 26 Male \n",
+ "2 32 24 55 Female \n",
+ "3 31 36 53 Male \n",
+ "4 45 28 51 Male "
+ ]
+ },
+ "execution_count": 12,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "# load NeuroCog data \n",
+ "df = pd.read_csv('NeuroCog_dataset.csv', sep=',')\n",
+ "\n",
+ "df.head()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "04a9875a-c70f-4d36-8046-d5651092109a",
+ "metadata": {},
+ "source": [
+ "## "
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "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.11.3"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/NeuroCog_dataset.csv b/NeuroCog_dataset.csv
new file mode 100644
index 0000000000000000000000000000000000000000..af4bb793715645760ce54a907de5d003ba2f46e6
--- /dev/null
+++ b/NeuroCog_dataset.csv
@@ -0,0 +1,243 @@
+,Dx,Speed,Attention,Memory,Verbal,Visual,ProbSolv,SocialCog,Age,Sex
+14,Schizophrenia,19,9,19,33,24,39,28,44,Female
+15,Schizophrenia,8,25,15,28,24,40,37,26,Male
+16,Schizophrenia,14,23,15,20,13,32,24,55,Female
+17,Schizophrenia,7,18,14,34,16,31,36,53,Male
+18,Schizophrenia,21,9,35,28,29,45,28,51,Male
+19,Schizophrenia,31,10,26,29,21,33,28,21,Male
+20,Schizophrenia,-1,8,3,20,12,29,28,53,Male
+21,Schizophrenia,17,20,27,30,32,29,44,56,Female
+22,Schizophrenia,7,30,26,26,27,30,39,48,Male
+23,Schizophrenia,37,15,17,33,21,33,24,46,Female
+24,Schizophrenia,30,27,28,34,19,30,32,48,Male
+25,Schizophrenia,26,20,22,33,18,39,36,31,Male
+26,Schizophrenia,32,23,22,29,21,31,26,45,Female
+27,Schizophrenia,36,16,38,27,23,36,28,27,Male
+28,Schizophrenia,29,42,31,28,23,37,18,44,Male
+29,Schizophrenia,21,35,42,33,35,32,27,30,Male
+30,Schizophrenia,48,32,24,37,23,42,24,43,Male
+31,Schizophrenia,36,16,30,36,37,50,23,41,Female
+32,Schizophrenia,36,37,31,32,26,30,36,36,Female
+33,Schizophrenia,46,28,33,33,21,37,32,48,Male
+34,Schizophrenia,29,28,31,28,35,45,39,52,Male
+35,Schizophrenia,27,42,35,25,19,39,49,45,Male
+36,Schizophrenia,33,21,38,28,37,48,29,27,Male
+37,Schizophrenia,32,41,38,33,24,35,39,53,Male
+38,Schizophrenia,30,36,35,42,32,35,32,33,Male
+39,Schizophrenia,32,50,38,33,21,41,27,51,Male
+40,Schizophrenia,39,31,40,28,26,37,39,45,Male
+41,Schizophrenia,45,35,36,29,19,46,32,20,Male
+42,Schizophrenia,27,21,35,34,34,39,57,59,Male
+43,Schizophrenia,28,50,47,37,18,37,33,53,Male
+44,Schizophrenia,39,34,47,30,21,40,44,49,Male
+45,Schizophrenia,22,28,37,30,46,50,40,20,Female
+46,Schizophrenia,29,29,42,36,40,50,29,35,Male
+47,Schizophrenia,47,32,39,39,24,39,34,53,Female
+48,Schizophrenia,42,31,51,30,38,50,10,40,Male
+49,Schizophrenia,25,43,42,36,46,42,24,50,Female
+50,Schizophrenia,39,34,44,29,27,41,42,57,Male
+51,Schizophrenia,32,36,37,29,35,50,40,28,Male
+52,Schizophrenia,51,32,38,40,29,37,41,45,Male
+53,Schizophrenia,34,50,46,30,41,37,32,51,Male
+54,Schizophrenia,25,38,44,40,38,32,60,31,Male
+55,Schizophrenia,43,35,49,48,43,36,29,53,Female
+56,Schizophrenia,43,17,45,40,40,50,45,40,Female
+57,Schizophrenia,38,32,38,29,48,52,42,43,Male
+58,Schizophrenia,44,31,42,33,37,61,32,28,Male
+59,Schizophrenia,39,48,45,37,24,37,53,44,Male
+60,Schizophrenia,30,46,47,44,34,40,50,30,Male
+61,Schizophrenia,37,48,53,46,38,35,39,63,Male
+62,Schizophrenia,35,40,49,44,44,46,56,25,Female
+63,Schizophrenia,46,57,52,67,34,39,27,50,Female
+64,Schizophrenia,39,53,47,48,48,41,52,25,Male
+65,Schizophrenia,41,39,44,42,57,55,50,40,Female
+66,Schizophrenia,34,43,56,51,48,65,37,25,Male
+67,Schizophrenia,48,54,51,48,46,55,31,31,Male
+68,Schizophrenia,44,51,51,53,51,52,51,29,Female
+69,Schizophrenia,48,55,46,56,57,45,52,37,Female
+70,Schizophrenia,49,63,49,42,55,50,50,30,Male
+71,Schizophrenia,55,43,63,51,55,65,41,26,Male
+79,Schizoaffective,25,11,35,30,24,45,32,49,Male
+80,Schizoaffective,26,29,28,27,26,33,32,52,Female
+81,Schizoaffective,22,21,22,32,21,31,24,64,Male
+82,Schizoaffective,14,1,20,23,15,36,37,37,Male
+83,Schizoaffective,23,30,28,34,23,29,32,51,Male
+84,Schizoaffective,32,10,26,33,19,55,24,35,Male
+85,Schizoaffective,30,31,13,25,19,35,28,56,Female
+86,Schizoaffective,13,19,15,27,26,30,44,49,Male
+87,Schizoaffective,26,20,22,32,18,37,36,46,Male
+88,Schizoaffective,27,44,24,33,26,31,36,53,Female
+89,Schizoaffective,22,23,33,27,12,45,60,37,Male
+90,Schizoaffective,33,19,20,37,43,41,19,40,Male
+91,Schizoaffective,31,38,33,26,32,35,31,40,Female
+92,Schizoaffective,32,33,32,42,26,30,37,50,Female
+93,Schizoaffective,39,31,40,40,34,35,27,46,Male
+94,Schizoaffective,34,38,38,28,24,39,53,47,Female
+95,Schizoaffective,34,30,45,32,35,40,37,39,Male
+96,Schizoaffective,42,33,47,30,32,35,39,56,Male
+97,Schizoaffective,43,26,42,33,38,45,29,37,Male
+98,Schizoaffective,24,31,33,40,41,33,59,42,Female
+99,Schizoaffective,46,21,36,42,37,41,40,28,Male
+100,Schizoaffective,44,42,44,33,32,37,31,35,Female
+101,Schizoaffective,31,34,36,36,24,55,53,33,Male
+102,Schizoaffective,34,41,40,40,40,29,45,51,Female
+103,Schizoaffective,35,47,51,48,26,31,37,57,Male
+104,Schizoaffective,27,34,35,29,55,52,41,34,Female
+105,Schizoaffective,23,41,42,48,46,29,50,52,Female
+106,Schizoaffective,29,35,38,48,41,41,56,50,Female
+107,Schizoaffective,30,47,37,37,37,52,57,27,Female
+108,Schizoaffective,37,45,42,34,44,58,40,37,Male
+109,Schizoaffective,37,38,51,40,37,48,51,52,Male
+110,Schizoaffective,38,51,61,34,43,37,39,40,Male
+111,Schizoaffective,37,45,54,51,34,37,51,43,Female
+112,Schizoaffective,39,49,44,27,43,42,65,39,Male
+113,Schizoaffective,44,42,58,59,35,32,43,44,Male
+114,Schizoaffective,41,53,47,39,30,50,63,32,Female
+115,Schizoaffective,52,58,47,48,35,40,65,43,Male
+116,Schizoaffective,45,53,58,48,41,48,53,35,Male
+117,Schizoaffective,52,54,68,53,49,58,48,32,Male
+119,Control,25,17,22,32,23,31,24,30,Male
+120,Control,27,14,17,36,16,30,15,55,Female
+121,Control,24,17,20,34,21,37,28,35,Female
+122,Control,2,2,10,29,16,35,27,51,Male
+123,Control,39,34,15,26,26,43,23,30,Male
+124,Control,21,17,26,29,26,52,51,37,Female
+125,Control,26,20,40,32,27,36,42,53,Female
+126,Control,24,23,38,39,46,33,23,37,Male
+127,Control,26,46,20,33,19,41,37,37,Male
+128,Control,34,39,29,34,32,34,24,39,Male
+129,Control,40,12,27,40,38,39,29,30,Male
+130,Control,44,25,21,37,32,43,29,43,Female
+131,Control,33,28,35,37,23,48,29,31,Male
+132,Control,37,34,27,30,18,41,43,49,Female
+133,Control,36,16,38,37,38,48,24,31,Male
+134,Control,38,33,37,32,30,42,28,48,Male
+135,Control,37,39,29,34,35,46,29,63,Male
+136,Control,37,37,38,34,34,45,32,48,Male
+137,Control,51,19,46,33,29,40,41,38,Male
+138,Control,33,63,29,39,26,37,39,63,Male
+139,Control,32,32,40,42,21,41,55,31,Male
+140,Control,28,30,42,34,30,35,64,46,Male
+141,Control,35,37,40,44,37,30,42,34,Male
+142,Control,37,32,35,37,40,45,43,53,Male
+143,Control,43,39,26,42,30,34,51,49,Female
+144,Control,35,28,48,28,38,50,50,42,Female
+145,Control,33,32,38,42,55,33,42,22,Female
+146,Control,43,50,40,36,35,45,26,57,Female
+147,Control,45,30,55,42,46,42,17,38,Female
+148,Control,41,39,40,36,34,46,42,34,Male
+149,Control,42,57,42,42,27,39,29,55,Male
+150,Control,35,45,44,40,30,41,44,38,Male
+151,Control,44,48,44,34,29,43,40,47,Male
+152,Control,23,35,35,46,40,46,53,60,Female
+153,Control,40,43,44,36,43,42,34,56,Male
+154,Control,42,39,47,39,29,42,47,56,Male
+155,Control,38,46,36,44,29,39,57,40,Female
+156,Control,36,36,47,51,37,40,43,52,Female
+157,Control,39,38,49,34,30,48,49,20,Male
+158,Control,41,47,39,40,43,39,45,30,Male
+159,Control,35,52,44,40,32,41,48,35,Male
+160,Control,39,38,55,44,34,58,27,35,Male
+161,Control,47,36,44,32,34,45,57,38,Male
+162,Control,56,36,44,39,40,45,37,48,Female
+163,Control,44,40,49,32,24,50,60,53,Male
+164,Control,40,38,37,40,43,58,44,35,Male
+165,Control,35,35,43,53,57,37,40,30,Male
+166,Control,46,37,49,51,43,43,33,55,Female
+167,Control,38,31,45,37,49,58,47,51,Male
+168,Control,36,44,49,56,38,32,47,64,Female
+169,Control,52,48,42,42,23,34,60,60,Female
+170,Control,39,46,42,46,32,43,57,43,Male
+171,Control,36,37,45,39,43,58,45,39,Female
+172,Control,40,43,26,42,51,48,50,38,Female
+173,Control,42,53,44,46,35,43,44,32,Female
+174,Control,40,32,49,44,46,46,50,42,Female
+175,Control,46,29,40,37,35,50,71,60,Female
+176,Control,48,54,45,42,35,43,45,52,Female
+177,Control,48,54,45,42,35,43,45,53,Female
+178,Control,46,40,44,39,43,48,50,36,Male
+179,Control,40,16,45,48,55,58,48,19,Male
+180,Control,53,44,36,42,38,55,41,18,Male
+181,Control,44,35,44,33,44,45,65,43,Male
+182,Control,41,46,33,53,35,43,65,56,Female
+183,Control,41,53,38,48,32,48,53,55,Male
+184,Control,54,58,44,36,41,43,41,27,Male
+185,Control,41,50,54,36,30,43,63,44,Male
+186,Control,40,45,42,40,35,42,72,39,Male
+187,Control,43,34,42,39,41,52,64,44,Female
+188,Control,48,28,45,48,43,52,52,29,Male
+189,Control,50,48,51,39,41,52,39,48,Male
+190,Control,57,36,45,46,38,50,45,21,Male
+191,Control,43,33,53,48,51,52,39,29,Male
+192,Control,36,41,51,46,43,55,50,42,Female
+193,Control,48,54,49,32,29,48,59,53,Male
+194,Control,47,40,53,53,40,45,49,55,Female
+195,Control,50,52,38,42,38,58,45,25,Male
+196,Control,37,44,52,46,40,50,58,42,Female
+197,Control,46,47,36,44,46,48,56,34,Male
+198,Control,46,53,36,48,52,35,56,55,Female
+199,Control,42,41,65,46,29,61,43,28,Male
+200,Control,46,45,49,42,37,50,59,37,Male
+201,Control,48,43,44,46,40,46,63,38,Female
+202,Control,47,44,36,48,38,52,63,52,Female
+203,Control,45,52,43,51,49,58,32,25,Female
+204,Control,50,49,54,59,48,45,24,44,Male
+205,Control,38,50,53,44,40,39,67,50,Male
+206,Control,52,41,51,44,40,52,53,37,Female
+207,Control,40,39,49,56,43,55,50,24,Female
+208,Control,45,49,49,56,30,45,59,47,Female
+209,Control,48,44,51,37,51,61,44,27,Male
+210,Control,49,49,56,59,35,45,47,53,Male
+211,Control,49,43,45,53,51,42,53,42,Female
+212,Control,47,42,53,44,54,61,37,27,Male
+213,Control,49,54,53,44,49,48,40,32,Male
+214,Control,39,52,51,40,38,61,59,48,Female
+215,Control,50,36,58,39,49,52,57,36,Male
+216,Control,49,50,38,39,55,48,64,66,Male
+217,Control,57,54,56,37,43,61,36,42,Male
+218,Control,53,60,58,32,37,48,60,56,Male
+219,Control,56,57,22,56,51,50,55,30,Male
+220,Control,52,35,58,53,49,55,44,35,Male
+221,Control,40,52,51,51,32,61,63,34,Female
+222,Control,50,67,47,48,48,40,49,55,Female
+223,Control,60,55,56,37,49,50,43,19,Male
+224,Control,57,51,60,42,37,61,49,37,Female
+225,Control,49,42,51,48,41,61,62,35,Female
+226,Control,54,47,62,42,48,61,44,31,Male
+227,Control,47,47,54,53,40,65,56,45,Male
+228,Control,49,56,51,53,41,61,48,45,Female
+229,Control,57,49,53,63,43,46,47,42,Female
+230,Control,62,46,63,48,40,52,48,23,Male
+231,Control,50,42,63,53,46,43,60,33,Female
+232,Control,52,47,47,67,49,52,51,30,Male
+233,Control,62,58,44,51,59,58,34,37,Female
+234,Control,50,54,52,56,43,50,58,34,Male
+235,Control,49,60,52,53,44,45,63,46,Female
+236,Control,53,55,47,56,59,61,39,31,Male
+237,Control,51,59,60,53,38,55,56,48,Male
+238,Control,45,51,51,59,48,43,71,50,Female
+239,Control,51,46,49,56,40,65,60,57,Male
+240,Control,53,51,51,59,44,61,51,62,Male
+241,Control,55,48,55,40,44,65,60,29,Female
+242,Control,58,55,58,53,43,45,58,29,Female
+243,Control,56,51,54,46,55,61,43,44,Male
+244,Control,49,42,47,59,65,48,63,47,Male
+245,Control,54,63,56,59,46,65,29,43,Male
+246,Control,53,46,63,46,52,55,59,30,Female
+247,Control,53,63,51,42,40,58,68,24,Male
+248,Control,58,49,53,51,44,55,67,33,Male
+249,Control,46,57,69,59,35,58,55,49,Male
+250,Control,54,53,65,56,41,61,57,32,Female
+251,Control,56,61,54,56,52,61,47,52,Female
+252,Control,50,56,53,78,48,65,39,32,Male
+253,Control,61,59,58,56,41,58,55,33,Female
+254,Control,73,48,49,63,46,65,49,33,Male
+255,Control,45,61,47,59,51,65,63,31,Male
+256,Control,45,57,61,63,63,65,42,29,Male
+257,Control,68,40,53,53,59,65,59,24,Female
+258,Control,51,56,56,63,52,61,55,33,Female
+259,Control,60,63,61,46,60,61,44,35,Male
+260,Control,74,54,56,59,57,58,41,29,Male
+261,Control,53,62,62,59,55,65,48,34,Male
+262,Control,74,62,50,53,52,65,52,24,Male
+263,Control,60,64,71,51,57,65,60,40,Male