{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "# Python: Group Average Treatment Effects (GATEs) for IRM models\n",
    "\n",
    "In this simple example, we illustrate how the [DoubleML](https://docs.doubleml.org/stable/index.html) package can be used to estimate group average treatment effects in the [DoubleMLIRM](https://docs.doubleml.org/stable/guide/models.html#interactive-regression-model-irm) model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import doubleml as dml\n",
    "\n",
    "from doubleml.datasets import make_heterogeneous_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data\n",
    "\n",
    "We define a data generating process to create synthetic data to compare the estimates to the true effect. The data generating process is based on the Monte Carlo simulation from [Oprescu et al. (2019)](http://proceedings.mlr.press/v97/oprescu19a.html).\n",
    "\n",
    "The documentation of the data generating process can be found [here](https://docs.doubleml.org/dev/api/api.html#dataset-generators). In this example the true effect depends only the first covariate $X_0$ and takes the following form\n",
    "\n",
    "$$\n",
    "\\theta_0(X) = \\exp(2X_0) + 3\\sin(4X_0).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          y    d       X_0       X_1       X_2       X_3       X_4       X_5  \\\n",
      "0  6.114530  1.0  0.925248  0.180575  0.567945  0.915488  0.033946  0.697420   \n",
      "1  5.580922  1.0  0.474214  0.862043  0.844549  0.319100  0.828915  0.037008   \n",
      "2  1.278434  0.0  0.696289  0.339875  0.724767  0.065356  0.315290  0.539491   \n",
      "3  1.794805  0.0  0.615863  0.232959  0.024401  0.870099  0.021269  0.874702   \n",
      "4  6.178169  1.0  0.350712  0.767188  0.401931  0.479876  0.627505  0.873677   \n",
      "\n",
      "        X_6       X_7       X_8       X_9  \n",
      "0  0.297349  0.924396  0.971058  0.944266  \n",
      "1  0.596270  0.230009  0.120567  0.076953  \n",
      "2  0.790723  0.318753  0.625891  0.885978  \n",
      "3  0.528937  0.939068  0.798783  0.997934  \n",
      "4  0.984083  0.768273  0.417767  0.421357  \n"
     ]
    }
   ],
   "source": [
    "np.random.seed(42)\n",
    "data_dict = make_heterogeneous_data(\n",
    "    n_obs=500,\n",
    "    p=10,\n",
    "    support_size=5,\n",
    "    n_x=1,\n",
    "    binary_treatment=True,\n",
    ")\n",
    "data = data_dict['data']\n",
    "print(data.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The generated dictionary also contains the true individual effects saved in the key `effects`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[4.770944   5.4235839  5.07202564 5.30917769 4.97441062]\n"
     ]
    }
   ],
   "source": [
    "ite = data_dict['effects']\n",
    "print(ite[:5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal is to estimate the average treatment effect for different groups based on the covariate $X_0$. The groups can be specified as [DataFrame](https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.html) with boolean columns. We consider the following three groups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   Group 1  Group 2  Group 3\n",
      "0    False    False     True\n",
      "1    False     True    False\n",
      "2    False     True    False\n",
      "3    False     True    False\n",
      "4    False     True    False\n"
     ]
    }
   ],
   "source": [
    "groups = pd.DataFrame(\n",
    "    np.column_stack((data['X_0'] <= 0.3,\n",
    "                     (data['X_0'] > 0.3) & (data['X_0'] <= 0.7),\n",
    "                     data['X_0'] > 0.7)),\n",
    "    columns=['Group 1', 'Group 2', 'Group 3'])\n",
    "print(groups.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The true effects (still including sampling uncertainty) are given by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2.906716732639898, 5.223485956098176, 4.827938162750831]\n"
     ]
    }
   ],
   "source": [
    "true_effects = [ite[groups[group]].mean() for group in groups.columns]\n",
    "print(true_effects)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interactive Regression Model (IRM)\n",
    "The first step is to fit a [DoubleML IRM Model](https://docs.doubleml.org/stable/guide/models.html#interactive-regression-model-irm) to the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_dml_base = dml.DoubleMLData(\n",
    "    data,\n",
    "    y_col='y',\n",
    "    d_cols='d'\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training IRM Model\n",
      "       coef   std err          t  P>|t|     2.5 %    97.5 %\n",
      "d  4.482012  0.086889  51.583034    0.0  4.311712  4.652312\n"
     ]
    }
   ],
   "source": [
    "# First stage estimation\n",
    "from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor\n",
    "ml_g = RandomForestRegressor(n_estimators=500)\n",
    "ml_m = RandomForestClassifier(n_estimators=500)\n",
    "\n",
    "np.random.seed(42)\n",
    "\n",
    "dml_irm = dml.DoubleMLIRM(data_dml_base,\n",
    "                          ml_g=ml_g,\n",
    "                          ml_m=ml_m,\n",
    "                          n_folds=5)\n",
    "print(\"Training IRM Model\")\n",
    "dml_irm.fit()\n",
    "\n",
    "print(dml_irm.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Group Average Treatment Effects (GATEs)\n",
    "To calculate GATEs just call the ``gate()`` method and supply the DataFrame with the group definitions and the ``level`` (with default of ``0.95``). Remark that for straightforward interpretation of the GATEs the groups should be mutually exclusive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            2.5 %    effect    97.5 %\n",
      "Group 1  2.701265  3.016315  3.331365\n",
      "Group 2  5.096550  5.314651  5.532751\n",
      "Group 3  4.412004  4.668981  4.925957\n"
     ]
    }
   ],
   "source": [
    "gate = dml_irm.gate(groups=groups)\n",
    "print(gate.confint(level=0.95))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The confidence intervals above are point-wise, but by setting the option ``joint`` and providing a number of bootstrap repetitions ``n_rep_boot``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            2.5 %    effect    97.5 %\n",
      "Group 1  2.441219  3.016315  3.591411\n",
      "Group 2  4.916528  5.314651  5.712774\n",
      "Group 3  4.199893  4.668981  5.138068\n"
     ]
    }
   ],
   "source": [
    "ci = gate.confint(level=0.95, joint=True, n_rep_boot=1000)\n",
    "print(ci)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let us plot the estimates together with the true effect within each group."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA04AAAKXCAYAAAC15hrSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABUfElEQVR4nO3deVyVZf7/8fdhR2VxQUFFxVAU9yUTydQZUtIxzbZhbEwzp0V/6oy22DKilmhlZZlmNWmNmm1kU5qlJu7mlolrLiguIOUCIokE9+8Pvpw8At7nKHCO8no+Hvejc1/3dd/355x8nHx3Xfd1LIZhGAIAAAAAlMrN2QUAAAAAgKsjOAEAAACACYITAAAAAJggOAEAAACACYITAAAAAJggOAEAAACACYITAAAAAJggOAEAAACACYITAAAAAJggOAEAAACACYITAMClpaSkaMSIEWratKmqVKmiKlWqKDIyUsOHD9eOHTtKPOfJJ5+UxWLR/fffb9NusVjs2pKSknT48OEr9pkyZUpFvH0AgIuwGIZhOLsIAABK8vXXX+v++++Xh4eHBg4cqDZt2sjNzU179+5VYmKijhw5opSUFDVs2NB6jmEYatCggTw8PHTy5EmdPHlSfn5+kqR58+bZXP/DDz/UsmXL9N///tem/fbbb9dvv/2msLAwxcXFqXfv3sVqa9eunVq0aFEO7xoA4Io8nF0AAAAlOXjwoP7617+qYcOGWrFihUJCQmyOT506VTNnzpSbm+3kiaSkJB07dkzff/+9evXqpcTERD344IOSpAceeMCm78aNG7Vs2bJi7ZJ0+PBhSVL79u1LPA4AqFyYqgcAcEkvvfSSzp8/rzlz5hQLTZLk4eGhkSNHKjQ01KZ9/vz5ioyMVI8ePRQTE6P58+eXe61btmxRr169VKtWLfn6+iosLEwPPfRQud8XAFBxGHECALikr7/+WuHh4brlllvsPic3N1eff/65xowZI0mKi4vTkCFDlJ6eruDg4KuqIycnR7/++mux9sDAQHl4eCgjI0M9e/ZUUFCQnn76aQUGBurw4cNKTEy8qvsBAFwTI04AAJeTlZWlEydOqGXLlsWOnT17Vr/++qt1++2336zHvv76a509e1Z//etfJUn9+/eXp6enFi5ceNW1jB8/XkFBQcW2LVu2SJLWr1+vM2fOaP78+Ro7dqwefvhhvfDCC9q9e/dV3xMA4HoITgAAl5OVlSVJqlatWrFj3bt3twkwb731lvXY/Pnz1bFjR4WHh0uS/Pz81KdPn2uarvePf/xDy5YtK7ZFRkZKKhx5kgpDW15e3lXfBwDg2piqBwBwOUWr4GVnZxc7Nnv2bJ07d04nT560WbTh7NmzWrJkiUaMGKEDBw5Y26Ojo/X555/r559/VtOmTR2upUmTJoqJiSn1eLdu3XT33XdrwoQJeu2119S9e3f1799ff/vb3+Tt7e3w/QAArokRJwCAywkICFBISIh27txZ7Ngtt9yimJgYRUdH27R/+umnys3N1bRp09SkSRPr9q9//UuSym2RCIvFos8++0wbNmzQiBEjdPz4cT300EPq0KFDicEPAHB9IjgBAFxSnz59dODAAW3atMmu/vPnz1fLli316aefFttiYmK0YMGCcq23c+fOevHFF7VlyxbNnz9fu3btuqZnqwAAroWpegAAl/Tkk09qwYIFeuihh7RixQrVqVPH5vilv99+9OhRrV69WhMmTNA999xT7FoXL17UwIED9cMPPzi0Sp89zpw5o8DAQFksFmtb27ZtJRWu8gcAuDEQnAAALqlJkyZasGCB4uLiFBERoYEDB6pNmzYyDEMpKSlasGCB3NzcVL9+fS1YsECGYejOO+8s8Vq9e/eWh4eH5s+f73Bw2rZtm+bNm1es/aabblJUVJQ++OADzZw5U3fddZduuukmnTt3Tu+++678/f3Vu3fvq3rvAADXYzEu/V92AAC4mIMHD2ratGlatmyZjh07JovFooYNG6p79+569NFH1aZNG7Vu3VqZmZk6cuRIqdfp0aOHdu/erePHj8vDo/D/G44YMUJvvfWWSvpP4eHDhxUWFlbq9R588EHNnTtXP/74o15++WWtW7dOJ0+eVEBAgDp16qT4+Hh16NDh2j8AAIBLIDgBAAAAgAkWhwAAAAAAEwQnAAAAADBBcAIAAAAAEwQnAAAAADBBcAIAAAAAEwQnAAAAADBR6X4At6CgQCdOnJCfn5/Nr7wDAAAAqFwMw9C5c+dUt25dubldeUyp0gWnEydOKDQ01NllAAAAAHARR48eVf369a/Yp9IFJz8/P0mFH46/v7+TqwEAAADgLFlZWQoNDbVmhCupdMGpaHqev78/wQkAAACAXY/wsDgEAAAAAJggOAEAAACACYITAAAAAJiodM84AQCAG19+fr7y8vKcXQYAJ/P09JS7u3uZXIvgBAAAbhiGYSg9PV1nz551dikAXERgYKCCg4Ov+TdcCU4AAOCGURSaateurSpVqvBj90AlZhiGcnJylJGRIUkKCQm5pusRnAAAwA0hPz/fGppq1qzp7HIAuABfX19JUkZGhmrXrn1N0/ZYHAIAANwQip5pqlKlipMrAeBKir4TrvW5R4ITAAC4oTA9D8Clyuo7geAEAAAAACYITgAAALDb3r171blzZ/n4+Kht27altgE3GoITAACAE1kslitu8fHxFVZL9+7dS6zh0UcftfYZP368qlatqn379mnFihWltl0ri8WiRYsWlcm1gLLAqnoAAACXy8+X1qyR0tKkkBCpa1epjH5E83JpaWnW1x9//LH+/e9/a9++fda2atWqWV8bhqH8/Hx5eJTfX+GGDRumiRMn2rRduuDGwYMH1adPHzVs2PCKbcCNhhEnAACASyUmSo0aST16SH/7W+E/GzUqbC8HwcHB1i0gIEAWi8W6v3fvXvn5+embb75Rhw4d5O3trbVr12rw4MHq37+/zXVGjx6t7t27W/cLCgqUkJCgsLAw+fr6qk2bNvrss89M66lSpYpNTcHBwfL395dUOAq0detWTZw40ToaVlKbJB09elT33XefAgMDVaNGDfXr10+HDx+2udf777+vFi1ayNvbWyEhIRoxYoQkqVGjRpKku+66SxaLxboPOBPBCQAAoEhionTPPdKxY7btx48XtpdTeDLz9NNPa8qUKdqzZ49at25t1zkJCQn68MMP9fbbb2vXrl365z//qQceeECrVq266jrS0tLUokULjRkzRmlpaRo7dmyJbXl5eerVq5f8/Py0Zs0arVu3TtWqVVNsbKwuXrwoSZo1a5aGDx+uf/zjH0pOTtb//vc/hYeHS5I2b94sSZozZ47S0tKs+4AzMVUPAABAKpyeN2qUZBjFjxmGZLFIo0dL/fqV27S90kycOFG333673f1zc3M1efJkLV++XFFRUZKkxo0ba+3atZo9e7a6detW6rkzZ87Ue++9Z9M2e/ZsDRw4UMHBwfLw8FC1atUUHBwsqXAq4eVt8+bNU0FBgd577z3rUtBz5sxRYGCgkpKS1LNnT73wwgsaM2aMRo0aZb3PzTffLEkKCgqSJAUGBlqvCTgbwQkAAEAqfKbp8pGmSxmGdPRoYb9LpsRVhI4dOzrU/8CBA8rJySkWti5evKh27dpd8dyBAwfq2WeftWmrU6eOQ/f/6aefdODAAfn5+dm0X7hwQQcPHlRGRoZOnDihP//5zw5dF3AmghMAAIBUuBBEWfYrQ1WrVrXZd3Nzk3HZyFheXp71dXZ2tiRp8eLFqlevnk0/b2/vK94rICDAOmXuamVnZ6tDhw6aP39+sWNBQUFyc+NpEVx/CE4AAABS4ep5ZdmvHAUFBWnnzp02bdu3b5enp6ckKTIyUt7e3kpNTb3itLzy0r59e3388ceqXbu2dWGJyzVq1EgrVqxQjx49Sjzu6emp/Pz88iwTcAhxHwAAQCpccrx+/cJnmUpisUihoYX9nOxPf/qTtmzZog8//FD79+/X+PHjbYKUn5+fxo4dq3/+85/64IMPdPDgQW3btk1vvvmmPvjggyteOycnR+np6TbbmTNnHKpv4MCBqlWrlvr166c1a9YoJSVFSUlJGjlypI7933TI+Ph4TZs2TW+88Yb2799vra9IUbC6mvsD5YHgBAAAIBUu+DB9euHry8NT0f7rr1f4whAl6dWrl55//nk9+eSTuvnmm3Xu3DkNGjTIps+kSZP0/PPPKyEhQc2bN1dsbKwWL16ssLCwK1773XffVUhIiM0WFxfnUH1VqlTR6tWr1aBBAw0YMEDNmzfX0KFDdeHCBesI1IMPPqjXX39dM2fOVIsWLfSXv/xF+/fvt15j2rRpWrZsmUJDQ02fywIqgsW4fILsDS4rK0sBAQHKzMwsdegYAABcfy5cuKCUlBSFhYXJx8fn6i+UmFi4ut6lC0WEhhaGpgEDrrlOABXrSt8NjmQDnnECAAC41IABhUuOr1lTuBBESEjh9DwXGGkC4DwEJwBApZNz8XdF/vtbSdLuib1UxYv/HOIy7u4VvuQ4ANfGM04AAAAAYILgBAAAAAAmCE4AAAAAYILgBAAAAAAmCE4AAAAAYILgBAAAAAAmCE4AAACXybn4uxo9vViNnl6snIu/O7scAC6A4AQAAHADmjt3rgIDA51dhkPKquZFixYpPDxc7u7uGj16dKltFenvf/+7Jk+efFXnxsfHq23btlfsc/jwYVksFm3fvt2hazdq1Eivv/56mdZSkXbv3q369evr/Pnz5X4vghMAAMBl8gsM6+tNKadt9svD4MGDZbFYim2xsbF2nV/SX37vv/9+/fzzz+VQra2KDmglfU4Wi0ULFy609nnkkUd0zz336OjRo5o0aVKpbdciKSlJFotFZ8+eNe37008/acmSJRo5cuRV3Wvs2LFasWKFdX/w4MHq37//VV3rcps3b9Y//vGPq67FHo6Esx9//FH33nuv6tSpIx8fHzVp0kTDhg2z/lm+PCBGRkaqc+fOevXVVx2q6WoQnAAAAC6xdGeaYl5dZd0fPGezbp36vZbuTCvX+8bGxiotLc1m++ijj676er6+vqpdu3YZVug65syZU+yzKgoS2dnZysjIUK9evVS3bl35+fmV2FaR3nzzTd17772qVq3aVZ1frVo11axZs4yrKhQUFKQqVaq4RC1ff/21OnfurNzcXM2fP1979uzRvHnzFBAQoOeff77U84YMGaJZs2bp99/Ld1otwQkAAOD/LN2ZpsfmbdPJrFyb9vTMC3ps3rZyDU/e3t4KDg622apXry5JMgxD8fHxatCggby9vVW3bl3r6EX37t115MgR/fOf/7SOvkjFR4KKpli9//77atCggapVq6bHH39c+fn5eumllxQcHKzatWvrxRdftKnr1VdfVatWrVS1alWFhobq8ccfV3Z2tqTCUZchQ4YoMzPTeu/4+HhJUm5ursaOHat69eqpatWquuWWW5SUlGRz7blz56pBgwaqUqWK7rrrLp06dcquzyowMLDYZ+Xj46OkpCRrKPrTn/4ki8VSapskrV27Vl27dpWvr69CQ0M1cuRImylfubm5euqppxQaGipvb2+Fh4frP//5jw4fPqwePXpIkqpXry6LxaLBgweXWGt+fr4+++wz9e3b19o2Y8YMtWzZ0rq/aNEiWSwWvf3229a2mJgYPffcc5Jsp8fFx8frgw8+0Jdffmn9zC/9XA8dOqQePXqoSpUqatOmjTZs2HDFz/Ly0aDU1FT169dP1apVk7+/v+677z6dPHnSevzyqXpFo1+vvPKKQkJCVLNmTQ0fPlx5eXmSSv/zebmcnBwNGTJEvXv31v/+9z/FxMQoLCxMt9xyi1555RXNnj271Pdw++236/Tp01q1alWpfcoCwQkAAECF0/MmfLVbJU3KK2qb8NXucp+2V5LPP/9cr732mmbPnq39+/dr0aJFatWqlSQpMTFR9evX18SJE62jL6U5ePCgvvnmGy1dulQfffSR/vOf/6hPnz46duyYVq1apalTp+q5557TDz/8YD3Hzc1Nb7zxhnbt2qUPPvhA33//vZ588klJUpcuXfT666/L39/feu+xY8dKkkaMGKENGzZo4cKF2rFjh+69917FxsZq//79kqQffvhBQ4cO1YgRI7R9+3b16NFDL7zwwjV9Tl26dNG+ffusn1laWlqpbQcPHlRsbKzuvvtu7dixQx9//LHWrl2rESNGWK83aNAgffTRR3rjjTe0Z88ezZ49W9WqVVNoaKg+//xzSdK+ffuUlpam6dOnl1jTjh07lJmZqY4dO1rbunXrpt27d+uXX36RJK1atUq1atWyBqC8vDxt2LBB3bt3L3a9sWPH6r777rMZoezSpYv1+LPPPquxY8dq+/btatq0qeLi4uweiSkoKFC/fv2sIWTZsmU6dOiQ7r///iuet3LlSh08eFArV67UBx98oLlz52ru3LmS7P/z+e233+rXX3+1/tm63JWmg3p5ealt27Zas2aNXe/zanmU69UBAACuE5tSTist80Kpxw1JaZkXtCnltKJuKvupSl9//XWxqVzPPPOMnnnmGaWmpio4OFgxMTHy9PRUgwYN1KlTJ0lSjRo15O7uLj8/PwUHB1/xHgUFBXr//ffl5+enyMhI9ejRQ/v27dOSJUvk5uamiIgITZ06VStXrtQtt9wiSTYLKTRq1EgvvPCCHn30Uc2cOVNeXl4KCAiQxWKxuXdqaqrmzJmj1NRU1a1bV1LhX/iXLl2qOXPmaPLkyZo+fbpiY2Otf1Fu2rSp1q9fr6VLl5p+VnFxcXJ3d7dp2717txo0aGCdnlijRg1rTSW1JSQkaODAgdb316RJE73xxhvq1q2bZs2apdTUVH3yySdatmyZYmJiJEmNGze23q9GjRrWa1/pL/VHjhyRu7u7zbTJli1bqkaNGlq1apXuueceJSUlacyYMdbwtWnTJuXl5dkEoiLVqlWTr6+vcnNzS/z3PXbsWPXp00eSNGHCBLVo0UIHDhxQs2bNSq2xyIoVK5ScnKyUlBSFhoZKkj788EO1aNFCmzdv1s0331ziedWrV9eMGTPk7u6uZs2aqU+fPlqxYoWGDRtm95/PokBtT50lqVu3ro4cOXJV59qL4AQAACAp41zpoelq+jmqR48emjVrlk1b0V/O7733Xr3++utq3LixYmNj1bt3b/Xt21ceHo79Va5Ro0Y2z/fUqVNH7u7ucnNzs2nLyMiw7i9fvlwJCQnau3evsrKy9Pvvv+vChQvKyckp9dmY5ORk5efnq2nTpjbtubm51udj9uzZo7vuusvmeFRUlF3B6bXXXrOGmSJFAc1eP/30k3bs2KH58+db2wzDUEFBgVJSUpScnCx3d3d169bNoete7rfffpO3t7fNFDWLxaLbbrtNSUlJiomJ0e7du/X444/rpZde0t69e7Vq1SrdfPPNDj17VKR169bW1yEhIZKkjIwMuwLJnj17FBoaag1NUuHiC4GBgdqzZ0+pwalFixY2QTYkJETJyckO1W0Y1zaS6+vrq5ycnGu6hhmCEwAAgKTafj5l2s9RVatWVXh4eInHQkNDtW/fPi1fvlzLli3T448/rpdfflmrVq2Sp6en3fe4vK/FYimxraCgQFLhCmZ/+ctf9Nhjj+nFF19UjRo1tHbtWg0dOlQXL14s9S/22dnZcnd319atW4uNDF3tAgmXCg4OLvWzsld2drYeeeSREle6a9CggQ4cOHBN1y9Sq1Yt5eTk6OLFi/Ly8rK2d+/eXe+8847WrFmjdu3ayd/f3xqmVq1addWB7dJ/n0VhrejfZ3m50p8hexWF7L179yoqKsrhGk6fPq2bbrrJ4fMcwTNOAAAAkjqF1VBIgI9KfnRdskgKCfBRp7AaFVmWla+vr/r27as33nhDSUlJ2rBhg/X/6nt5eSk/P7/M77l161YVFBRo2rRp6ty5s5o2baoTJ07Y9Cnp3u3atVN+fr4yMjIUHh5usxVN12revLnNs1SStHHjxjJ/D6Vp3769du/eXay+8PBweXl5qVWrViooKCh1wYGiEGT2uRctpLB7926b9qLnnD799FPrs0zdu3fX8uXLtW7duhKfb7r03uXx77t58+Y6evSojh49am3bvXu3zp49q8jIyKu+rj319uzZU7Vq1dJLL71U4nGzZd937typdu3aXW2JdiE4AQAASHJ3s2h838K/HF4enor2x/eNlLtbadHq2uTm5io9Pd1m+/XXXyUVrj73n//8Rzt37tShQ4c0b948+fr6qmHDhpIKp+CtXr1ax48ft55TFsLDw5WXl6c333xThw4d0n//+1+bld+K7p2dna0VK1bo119/VU5Ojpo2baqBAwdq0KBBSkxMVEpKijZt2qSEhAQtXrxYkjRy5EgtXbpUr7zyivbv368ZM2bYNU1PKvxL9OWflaM/gPrUU09p/fr11sUp9u/fry+//NK6OESjRo304IMP6qGHHtKiRYuUkpKipKQkffLJJ5Kkhg0bymKx6Ouvv9Yvv/xiXWnwckFBQWrfvr3Wrl1r0966dWtVr15dCxYssAlOixYtUm5urqKjo0utvVGjRtqxY4f27dunX3/91bqC3bWKiYlRq1atNHDgQG3btk2bNm3SoEGD1K1bN5vFLRxlz5/PqlWr6r333tPixYt15513avny5Tp8+LC2bNmiJ598Uo8++mip1z98+LCOHz9ebPpmWSM4AQAA/J/YliGa9UB71fb3tmkPDvDRrAfaK7ZlSLnde+nSpQoJCbHZbr31VkmFK4q9++67io6OVuvWrbV8+XJ99dVX1ueFJk6cqMOHD+umm25SUFBQmdXUpk0bvfrqq5o6dapatmyp+fPnKyEhwaZPly5d9Oijj+r+++9XUFCQdcRgzpw5GjRokMaMGaOIiAj1799fmzdvVoMGDSRJnTt31rvvvqvp06erTZs2+u6776zLb5sZMmRIsc/qzTffdOi9tW7dWqtWrdLPP/+srl27ql27dvr3v/9t86zUrFmzdM899+jxxx9Xs2bNNGzYMGtAq1evniZMmKCnn35aderUsVmN73IPP/ywzbNUUuF0tq5du8pisVj/Pbdu3Vr+/v7q2LGjqlatWur1hg0bpoiICHXs2FFBQUFat26dQ++9NBaLRV9++aWqV6+u2267TTExMWrcuLE+/vjja7quvX8++/Xrp/Xr18vT01N/+9vf1KxZM8XFxSkzM/OKKy5+9NFH6tmzp/V/JJQXi3GtT2JdZ7KyshQQEKDMzEz5+/s7uxwAgBPkXPxdkf/+VpK0e2IvVfHikd8bwYULF5SSkqKwsDD5+Fzbc0jnLuSpVfx3kqS5Q25W1yZB5TbShBvfb7/9poiICH388cdX9fxOeQoJCdGkSZP08MMPO7uUq3Lx4kU1adJECxYsKHWU7krfDY5kA/5LAQAAcJlLQ1KnsBqEJlwTX19fffjhh2U6jfJa5eTkaN26dTp58qRatGjh7HKuWmpqqp555pkrTm0sKwQnAACAy1Tx8tDhKX2cXQZuIFda7MEZ3nnnHU2aNEmjR492uVEwRxQt6FERCE4AAABAJTN69GibHzeGORaHAAAAAAATBCcAAHBDqWTrXgEwUVbfCQQnAABwQ/D09JRU+NA7ABQp+k4o+o64WjzjBAAAbgju7u4KDAxURkaGJKlKlSqyWFgND6isDMNQTk6OMjIyFBgYKHd392u6HsEJAADcMIKDgyXJGp4AIDAw0PrdcC0ITgAA4IZhsVgUEhKi2rVrKy8vz9nlAHAyT0/Pax5pKkJwAgAANxx3d/cy+8sSAEhOXhwiPj5eFovFZmvWrFmp/efOnVusv4+PTwVWDAAAAKAycvqIU4sWLbR8+XLrvofHlUvy9/fXvn37rPs89AkAAACgvDk9OHl4eDj0sJbFYimTh7sAAAAAwF5O/x2n/fv3q27dumrcuLEGDhyo1NTUK/bPzs5Ww4YNFRoaqn79+mnXrl0VVCkAAACAysqpwemWW27R3LlztXTpUs2aNUspKSnq2rWrzp07V2L/iIgIvf/++/ryyy81b948FRQUqEuXLjp27Fip98jNzVVWVpbNBgAAAACOsBiGYTi7iCJnz55Vw4YN9eqrr2ro0KGm/fPy8tS8eXPFxcVp0qRJJfaJj4/XhAkTirVnZmbK39//mmsGAFx/ci7+rsh/fytJ2j2xl6p4OX3mOgDACbKyshQQEGBXNnD6VL1LBQYGqmnTpjpw4IBd/T09PdWuXbsr9h83bpwyMzOt29GjR8uqXAAAAACVhEsFp+zsbB08eFAhISF29c/Pz1dycvIV+3t7e8vf399mAwAAAABHODU4jR07VqtWrdLhw4e1fv163XXXXXJ3d1dcXJwkadCgQRo3bpy1/8SJE/Xdd9/p0KFD2rZtmx544AEdOXJEDz/8sLPeAgAAAIBKwKmTuo8dO6a4uDidOnVKQUFBuvXWW7Vx40YFBQVJklJTU+Xm9ke2O3PmjIYNG6b09HRVr15dHTp00Pr16xUZGemstwAAAACgEnCpxSEqgiMPgAEAbkwsDgEAkK7jxSEAAAAAwBURnAAAAADABMEJAAAAAEwQnAAAAADABMEJAAAAAEwQnAAAAADABMEJAAAAAEwQnAAAAADABMEJAAAAAEwQnAAAAADABMEJAAAAAEwQnAAAAADABMEJAAAAAEwQnAAAAADABMEJAAAAAEwQnAAAAADABMEJAAAAAEwQnAAAAADABMEJAFDp5BcY1tebUk7b7AMAUBKCEwCgUlm6M00xr66y7g+es1m3Tv1eS3emObEqAICrIzgBACqNpTvT9Ni8bTqZlWvTnp55QY/N20Z4AgCUiuAEAKgU8gsMTfhqt0qalFfUNuGr3UzbAwCUiOAEAKgUNqWcVlrmhVKPG5LSMi9oU8rpiisKAHDdIDgBACqFjHOlh6ar6QcAqFwITgCASqG2n0+Z9gMAVC4EJwBApdAprIZCAnxkKeW4RVJIgI86hdWoyLIAANcJghMAoFJwd7NofN9ISSoWnor2x/eNlLtbadEKAFCZEZwAAJVGbMsQzXqgvWr7e9u0Bwf4aNYD7RXbMsRJlQEAXJ2HswsAAKAixbYMUXR4LbWK/06SNHfIzeraJIiRJgDAFTHiBACodC4NSZ3CahCaAACmCE4AAAAAYILgBAAAAAAmCE4AAAAAYILgBACofPLz/3i9eo3tPgAAJSA4AQAql8REqXnkH/u975AaNSpsBwCgFAQnAEDlkZgo3XOPdOK4bfvx44XthCcAQCkITgCAyiE/Xxo1SjKM4seK2kaPZtoeAKBEBCcAQOWwZo107Fjpxw1DOnq0sB8AAJchOAEAKoe0tLLtBwCoVAhOAIDKISSkbPsBACoVghMAoHLo2lWqX1+yWEo+brFIoaGF/QAAuAzBCQBQObi7S9OnF76+PDwV7b/+emE/AAAuQ3ACAFQeAwZIn30mhdS1ba9fv7B9wADn1AUAcHkEJwBA5TJggLRn9x/7S76RUlIITQCAKyI4AQAqn0un493Wlel5AABTBCcAAAAAMEFwAgAAAAATBCcAAAAAMEFwAgAAAAATBCcAAAAAMEFwAgAAAAATBCcAAAAAMEFwAgAAAFAhci7+rkZPL1ajpxcr5+Lvzi7HIQQnAAAAADBBcAIAAAAAEwQnAAAAADBBcAIAAAAAEwQnAAAAADBBcAIAAAAAEwQnAAAAADBBcAIAAAAAEwQnAAAAADBBcAIAAAAAEwQnAAAAADBBcAIAAAAAEwQnAAAAADBBcAIAAAAAEwQnAAAAADBBcAIAAAAAE04NTvHx8bJYLDZbs2bNrnjOp59+qmbNmsnHx0etWrXSkiVLKqhaAAAAAJWV00ecWrRoobS0NOu2du3aUvuuX79ecXFxGjp0qH788Uf1799f/fv3186dOyuwYgAAAACVjdODk4eHh4KDg61brVq1Su07ffp0xcbG6oknnlDz5s01adIktW/fXjNmzKjAigEAAABUNk4PTvv371fdunXVuHFjDRw4UKmpqaX23bBhg2JiYmzaevXqpQ0bNpR3mQAAAAAqMQ9n3vyWW27R3LlzFRERobS0NE2YMEFdu3bVzp075efnV6x/enq66tSpY9NWp04dpaenl3qP3Nxc5ebmWvezsrLK7g0AAAAAqBScGpzuuOMO6+vWrVvrlltuUcOGDfXJJ59o6NChZXKPhIQETZgwoUyuBQAAAKBycvpUvUsFBgaqadOmOnDgQInHg4ODdfLkSZu2kydPKjg4uNRrjhs3TpmZmdbt6NGjZVozAAAAgBufSwWn7OxsHTx4UCEhISUej4qK0ooVK2zali1bpqioqFKv6e3tLX9/f5sNAAAAABzh1OA0duxYrVq1SocPH9b69et11113yd3dXXFxcZKkQYMGady4cdb+o0aN0tKlSzVt2jTt3btX8fHx2rJli0aMGOGstwAAAACgEnDqM07Hjh1TXFycTp06paCgIN16663auHGjgoKCJEmpqalyc/sj23Xp0kULFizQc889p2eeeUZNmjTRokWL1LJlS2e9BQAAAACVgFOD08KFC694PCkpqVjbvffeq3vvvbecKgIAAACA4lzqGScAAAAAcEUEJwAAAAAwQXACAAAAABMEJwAAAAAwQXACAAAAABMEJwAAAAAwQXACAAAAABNO/R0nAACcoYqXhw5P6ePsMgAA1xFGnAAAAADABMEJAAAAAEwQnAAAAADABMEJAAAAAEwQnAAAAADABMEJAAAAAEwQnAAAAADABMEJAAAAAEwQnAAAAABUjPz8P16vXmO77+IITgAAAADKX2Ki1Dzyj/3ed0iNGhW2XwcITgAAAADKV2KidM890onjtu3Hjxe2XwfhieAEAAAAoPzk50ujRkmGUfxYUdvo0S4/bY/gBAAAAKD8rFkjHTtW+nHDkI4eLeznwghOAAAAAMpPWlrZ9nMSghMAAACA8hMSUrb9nITgBAAAAKD8dO0q1a8vWSwlH7dYpNDQwn4ujOAEAAAAoPy4u0vTpxe+vjw8Fe2//nphPxdGcAIAAABQvgYMkD77TAqpa9tev35h+4ABzqnLAQQnAAAAAOVvwABpz+4/9pd8I6WkXBehSSI4AQAAAKgol07Hu62ry0/PuxTBCQAAAABMEJwAAAAAwATBCQAAAABMEJwAAAAAwATBCQAAAABMEJwAAAAAwATBCQAAAABMEJwAAAAAwATBCQAAAABMEJwAAAAAwATBCQAAAABMEJwAAAAAwATBCQAAAABMEJwAAAAAwATBCQAAAECFyC8wrK83pZy22Xd1BCcAAAAA5W7pzjTFvLrKuj94zmbdOvV7Ld2Z5sSq7EdwAgAAAFCulu5M02PztulkVq5Ne3rmBT02b9t1EZ4ITgAAAADKTX6BoQlf7VZJk/KK2iZ8tdvlp+0RnAAAAACUm00pp5WWeaHU44aktMwL2pRyuuKKugoEJwAAAADlJuNc6aHpavo5C8EJAAAAQLmp7edTpv2cheAEAAAAoNx0CquhkAAfWUo5bpEUEuCjTmE1KrIshxGcAAAAAJQbdzeLxveNlKRi4alof3zfSLm7lRatXAPBCQAAAEC5im0ZolkPtFdtf2+b9uAAH816oL1iW4Y4qTL7eTi7AAAAAAA3vtiWIYoOr6VW8d9JkuYOuVldmwS5/EhTEUacAAAAAFSIS0NSp7Aa101okghOAAAAAGCK4AQAAAAAJghOAAAAAGCC4AQAAAAAJghOAAAAAGCC4AQAAAAAJghOAAAAAGCC4AQAAAAAJghOAAAAAGCC4AQAAAAAJghOAAAAAGCC4AQAAAAAJghOAAAAAGCC4AQAAAAAJghOAAAAAGCC4AQAAAAAJlwmOE2ZMkUWi0WjR48utc/cuXNlsVhsNh8fn4orEgAAAECl5OHsAiRp8+bNmj17tlq3bm3a19/fX/v27bPuWyyW8iwNAAAAAJw/4pSdna2BAwfq3XffVfXq1U37WywWBQcHW7c6depUQJUAAAAAKjOnB6fhw4erT58+iomJsat/dna2GjZsqNDQUPXr10+7du26Yv/c3FxlZWXZbAAAAADgCKcGp4ULF2rbtm1KSEiwq39ERITef/99ffnll5o3b54KCgrUpUsXHTt2rNRzEhISFBAQYN1CQ0PLqnwAAAAAlYTTgtPRo0c1atQozZ8/3+4FHqKiojRo0CC1bdtW3bp1U2JiooKCgjR79uxSzxk3bpwyMzOt29GjR8vqLQAAAACoJJy2OMTWrVuVkZGh9u3bW9vy8/O1evVqzZgxQ7m5uXJ3d7/iNTw9PdWuXTsdOHCg1D7e3t7y9vYus7oBAAAAVD5OC05//vOflZycbNM2ZMgQNWvWTE899ZRpaJIKg1ZycrJ69+5dXmUCAAAAgPOCk5+fn1q2bGnTVrVqVdWsWdPaPmjQINWrV8/6DNTEiRPVuXNnhYeH6+zZs3r55Zd15MgRPfzwwxVePwAAAIDKwyV+x6k0qampcnP74zGsM2fOaNiwYUpPT1f16tXVoUMHrV+/XpGRkU6sEgAAAMCNzmIYhuHsIipSVlaWAgIClJmZKX9/f2eXAwAAAFQaORd/V+S/v5Uk7Z7YS1W8nDuO40g2sLvSN954w65+I0eOtPeSAAAAAHBdsDs4vfbaa6Z9LBYLwQkAAADADcfu4JSSklKedQAAAACAy3LaD+ACAAAAwPXC7uD0/fffKzIyUllZWcWOZWZmqkWLFlq9enWZFgcAAAAArsDu4PT6669r2LBhJa42ERAQoEceecSu56AAAAAA4Hpjd3D66aefFBsbW+rxnj17auvWrWVSFAAAAAC4EruD08mTJ+Xp6VnqcQ8PD/3yyy9lUhQAAAAAuBK7g1O9evW0c+fOUo/v2LFDISEhZVIUAAAAALgSu4NT79699fzzz+vChQvFjv32228aP368/vKXv5RpcQAAAADgCuz+HafnnntOiYmJatq0qUaMGKGIiAhJ0t69e/XWW28pPz9fzz77bLkVCgAAAADOYndwqlOnjtavX6/HHntM48aNk2EYkiSLxaJevXrprbfeUp06dcqtUAAAAABwFruDkyQ1bNhQS5Ys0ZkzZ3TgwAEZhqEmTZqoevXq5VUfAAAAADid3c84Xap69eq6+eabdfDgQXl5eZV1TQAAAADgUq4qOBV55JFHdPLkybKqBQAAAABc0jUFp6LnnAAAAADgRnZNwQkAAAAAKoNrCk7ffPON6tWrV1a1AAAAAIBLcmhVvcs1a9ZMy5cvV35+vm6++WaFhISUVV0AAAAA4DKuOjh9/vnnGjp0qJo2baq8vDzt27dPb731loYMGVKW9QEAAACA09k9VS87O9tmf8KECdq0aZM2bdqkH3/8UZ9++qmeffbZMi8QAAAAAJzN7uDUoUMHffnll9Z9Dw8PZWRkWPdPnjzJbzoBAAAAuCHZPVXv22+/1fDhwzV37ly99dZbmj59uu6//37l5+fr999/l5ubm+bOnVuOpQIAAACAc9gdnBo1aqTFixfro48+Urdu3TRy5EgdOHBABw4cUH5+vpo1ayYfH5/yrBUAAAAAnMLh5cjj4uK0efNm/fTTT+revbsKCgrUtm1bQhMAAACAG5ZDq+otWbJEe/bsUZs2bfTee+9p1apVGjhwoO644w5NnDhRvr6+5VUnAAAAADiN3SNOY8aM0ZAhQ7R582Y98sgjmjRpkrp166Zt27bJx8dH7dq10zfffFOetQIAAACAU9gdnObOnaslS5Zo4cKF2rx5s/773/9Kkry8vDRp0iQlJiZq8uTJ5VYoAAAAADiL3cGpatWqSklJkSQdPXq02DNNkZGRWrNmTdlWBwAAAAAuwO7glJCQoEGDBqlu3brq1q2bJk2aVJ51AQAAAIDLsHtxiIEDByo2NlaHDh1SkyZNFBgYWI5lAQAAAIDrcGhVvZo1a6pmzZrlVQsAAAAAuCS7p+plZGTY7G/fvl0PPvigoqOjdc899ygpKamsawMAAAAAl2B3cAoJCbGGp/Xr16tTp046cuSIoqOjlZWVpdtvv12rV68ut0IBAAAAwFnsnqpnGIb1dXx8vP7+97/rP//5j7Vt9OjRmjBhglasWFG2FQIAAACAk9k94nSpnTt3atiwYTZtw4YN044dO8qkKAAAAABwJQ4tDnHu3Dn5+PjIx8dH3t7eNsd8fHyUk5NTpsUBAAAAgCtwaMSpadOmql69ug4fPqwtW7bYHNu1a5fq1q1bpsUBAAAAgCuwe8Rp5cqVNvshISE2+ykpKfrHP/5RNlUBAAAAgAuxOzh169btisdHjRp1zcUAAAAAgCty6BmnIqmpqUpLS5Obm5saN27Mj+ICAAAAuKE5FJxmzpypqVOn6tixYzbtUVFRmj59ujp06FCmxQEAAAC4cVTx8tDhKX2cXcZVsXtxiFdeeUUvvviinnjiCc2ePVsRERGKj4/X4sWL1bhxY912223FFowAAAAAgBuBxbj0l22vICwsTDNnztQdd9whSfr555/VpUsXpaeny8PDQ6NGjdKePXv03XfflWvB1yorK0sBAQHKzMyUv7+/s8sBAAAA4CSOZAO7R5wyMjLUvHlz636TJk2UmZmpX375RZL00EMPacOGDVdZMgAAAAC4LruDU9OmTbVs2TLr/sqVK+Xl5aXg4GBJhT+Aa7FYyr5CAAAAAHAyuxeHGDdunB544AEtX75cPj4+SkxM1MiRI61hKSkpSS1btiy3QgEAAADAWex+xkmSvvnmG82bN0+5ubnq1auXhg0bZj126tQpSXL5pcl5xgkAAACA5Fg2cCg43QgITgAAAACkclocAgAAAAAqK4ITAAAAAJggOAEAAACACYITAAAAAJggOAEAAACACbt+x2nAgAF2XzAxMfGqiwEAAAAAV2TXiFNAQIB18/f314oVK7Rlyxbr8a1bt2rFihUKCAgot0IBAAAAwFnsGnGaM2eO9fVTTz2l++67T2+//bbc3d0lSfn5+Xr88cf5XSQAAAAANySHfwA3KChIa9euVUREhE37vn371KVLF506dapMCyxr/AAuAAAAAKmcfwD3999/1969e4u17927VwUFBY5eDgAAAABcnl1T9S41ZMgQDR06VAcPHlSnTp0kST/88IOmTJmiIUOGlHmBAAAAAOBsDgenV155RcHBwZo2bZrS0tIkSSEhIXriiSc0ZsyYMi8QAAAAAJzN4WecLpWVlSVJ19WzQjzjBAAAAEByLBs4POJ0KYIHAAAAgMrA4cUhTp48qb///e+qW7euPDw85O7ubrMBAAAAwI3G4RGnwYMHKzU1Vc8//7xCQkJksVjKoy4AAAAAcBkOB6e1a9dqzZo1atu2bTmUAwAAAACux+GpeqGhobqG9SQAAAAA4LrjcHB6/fXX9fTTT+vw4cPlUA4AAAAAuB6Hp+rdf//9ysnJ0U033aQqVarI09PT5vjp06fLrDgAAAAAcAUOB6fXX3+9HMqQpkyZonHjxmnUqFFXvMenn36q559/XocPH1aTJk00depU9e7du1xqAgAAAADpKoLTgw8+WOZFbN68WbNnz1br1q2v2G/9+vWKi4tTQkKC/vKXv2jBggXq37+/tm3bppYtW5Z5XQAAAAAgXcUzTpe6cOGCsrKybDZHZWdna+DAgXr33XdVvXr1K/adPn26YmNj9cQTT6h58+aaNGmS2rdvrxkzZlztWwAAAAAAUw4Hp/Pnz2vEiBGqXbu2qlatqurVq9tsjho+fLj69OmjmJgY074bNmwo1q9Xr17asGFDqefk5uZec7gDAAAAULk5HJyefPJJff/995o1a5a8vb313nvvacKECapbt64+/PBDh661cOFCbdu2TQkJCXb1T09PV506dWza6tSpo/T09FLPSUhIUEBAgHULDQ11qEYAAAAAcDg4ffXVV5o5c6buvvtueXh4qGvXrnruuec0efJkzZ8/3+7rHD16VKNGjdL8+fPl4+PjaBl2GzdunDIzM63b0aNHy+1eAAAAAG5MDi8Ocfr0aTVu3FiS5O/vb11+/NZbb9Vjjz1m93W2bt2qjIwMtW/f3tqWn5+v1atXa8aMGcrNzZW7u7vNOcHBwTp58qRN28mTJxUcHFzqfby9veXt7W13XQAAAABwOYdHnBo3bqyUlBRJUrNmzfTJJ59IKhyJCgwMtPs6f/7zn5WcnKzt27dbt44dO2rgwIHavn17sdAkSVFRUVqxYoVN27JlyxQVFeXo2wAAAAAAuzk84jRkyBD99NNP6tatm55++mn17dtXM2bMUF5enl599VW7r+Pn51dsCfGqVauqZs2a1vZBgwapXr161megRo0apW7dumnatGnq06ePFi5cqC1btuidd95x9G0AAAAAgN0cDk7//Oc/ra9jYmK0d+9ebd26VeHh4aa/w+So1NRUubn9MSjWpUsXLViwQM8995yeeeYZNWnSRIsWLeI3nAAAAACUK4thGIazi6hIWVlZCggIUGZmpvz9/Z1dDgAAAAAncSQbXNMP4AIAAABAZUBwAgAAAAATBCcAAAAAMEFwAgAAAAATDgcnd3d3ZWRkFGs/depUib+9BAAAAADXO4eDU2mL8OXm5srLy+uaCwIAAAAAV2P37zi98cYbkiSLxaL33ntP1apVsx7Lz8/X6tWr1axZs7KvEAAAAACczO7g9Nprr0kqHHF6++23bableXl5qVGjRnr77bfLvkIAAAAAcDK7g1NKSookqUePHkpMTFT16tXLrSgAAAAAcCV2B6ciK1euLI86AAAAAMBlObw4xN13362pU6cWa3/ppZd07733lklRAAAAAOBKHA5Oq1evVu/evYu133HHHVq9enWZFAUAAAAArsTh4JSdnV3isuOenp7Kysoqk6IAAAAAwJU4HJxatWqljz/+uFj7woULFRkZWSZFAQAAAIArcXhxiOeff14DBgzQwYMH9ac//UmStGLFCn300Uf69NNPy7xAAAAAAHA2h4NT3759tWjRIk2ePFmfffaZfH191bp1ay1fvlzdunUrjxoBAAAAwKkshmEYzi6iImVlZSkgIECZmZny9/d3djkAAAAAnMSRbODwM06SdPbsWb333nt65plndPr0aUnStm3bdPz48au5HAAAAAC4NIen6u3YsUMxMTEKCAjQ4cOH9fDDD6tGjRpKTExUamqqPvzww/KoEwAAAACcxuERp3/9618aPHiw9u/fLx8fH2t77969+R0nAAAAADckh4PT5s2b9cgjjxRrr1evntLT08ukKAAAAABwJQ4HJ29v7xJ/6Pbnn39WUFBQmRQFAAAAAK7E4eB05513auLEicrLy5MkWSwWpaam6qmnntLdd99d5gUCAAAAgLM5HJymTZum7Oxs1a5dW7/99pu6deum8PBw+fn56cUXXyyPGgEAAADAqRxeVS8gIEDLli3TunXr9NNPPyk7O1vt27dXTExMedQHAAAAAE5nV3CqUaOGfv75Z9WqVUsPPfSQpk+frujoaEVHR5d3fQAAAADgdHZN1bt48aJ1QYgPPvhAFy5cKNeiAAAAAMCV2DXiFBUVpf79+6tDhw4yDEMjR46Ur69viX3ff//9Mi0QAAAAAJzNruA0b948vfbaazp48KAkKTMzk1EnAAAAAJWGxTAMw5ETwsLCtGXLFtWsWbO8aipXWVlZCggIUGZmpvz9/Z1dDgAAAAAncSQb2PWMU40aNfTrr79Kknr06CEvL69rrxIAAAAArhMsDgEAAAAAJlgcAgAAAABMOLw4hMViYXEIAAAAAJUKi0MAAAAAqJQcyQZ2jThdKiUl5aoLAwAAAIDrkV2LQ0hS7969lZmZad2fMmWKzp49a90/deqUIiMjy7Q4AAAAAHAFdgenb7/9Vrm5udb9yZMn6/Tp09b933//Xfv27Svb6gAAAADABdgdnC5/FMrBR6MAAAAA4Lpld3ACAAAAgMrK7uBksVhksViKtQEAAADAjc7uVfUMw9DgwYPl7e0tSbpw4YIeffRRVa1aVZJsnn8CAAAAgBuJ3cHpwQcftNl/4IEHivUZNGjQtVcEAAAAAC7G7uA0Z86c8qwDAAAAAFwWi0MAAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCk5PkXPxdjZ5erEZPL1bOxd+dXQ4AAACAKyA4AQAAAIAJghMAAAAAmCA4AQAAAIAJghMAAAAAmCA4AQAAAIAJghMAAAAAmCA4AQAAAIAJghMAAAAAmCA4AQAAAIAJpwanWbNmqXXr1vL395e/v7+ioqL0zTfflNp/7ty5slgsNpuPj08FVgwAAACgMvJw5s3r16+vKVOmqEmTJjIMQx988IH69eunH3/8US1atCjxHH9/f+3bt8+6b7FYKqpcAAAAAJWUU4NT3759bfZffPFFzZo1Sxs3biw1OFksFgUHB1dEeQAAAAAgyYWeccrPz9fChQt1/vx5RUVFldovOztbDRs2VGhoqPr166ddu3ZVYJUAAAAAKiOnjjhJUnJysqKionThwgVVq1ZNX3zxhSIjI0vsGxERoffff1+tW7dWZmamXnnlFXXp0kW7du1S/fr1SzwnNzdXubm51v2srKxyeR8AAAAAblxOH3GKiIjQ9u3b9cMPP+ixxx7Tgw8+qN27d5fYNyoqSoMGDVLbtm3VrVs3JSYmKigoSLNnzy71+gkJCQoICLBuoaGh5fVWAAAAANygnB6cvLy8FB4erg4dOighIUFt2rTR9OnT7TrX09NT7dq104EDB0rtM27cOGVmZlq3o0ePllXpAAAAACoJpwenyxUUFNhMrbuS/Px8JScnKyQkpNQ+3t7e1uXOizYAAAAAcIRTn3EaN26c7rjjDjVo0EDnzp3TggULlJSUpG+//VaSNGjQINWrV08JCQmSpIkTJ6pz584KDw/X2bNn9fLLL+vIkSN6+OGHnfk2AAAAANzgnBqcMjIyNGjQIKWlpSkgIECtW7fWt99+q9tvv12SlJqaKje3PwbFzpw5o2HDhik9PV3Vq1dXhw4dtH79+lIXkwAAAACAsmAxDMNwdhEVKSsrSwEBAcrMzHTqtL2ci78r8t+FI2u7J/ZSFS+nL3AIAAAAVCqOZAOXe8YJAAAAAFwNwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcAAAAAMAEwQkAAAAATBCcnCS/wLC+3pRy2mYfAAAAgGshODnB0p1pinl1lXV/8JzNunXq91q6M82JVQEAAAAoDcGpgi3dmabH5m3Tyaxcm/b0zAt6bN42whMAAADggpwanGbNmqXWrVvL399f/v7+ioqK0jfffHPFcz799FM1a9ZMPj4+atWqlZYsWVJB1V67/AJDE77arZIm5RW1TfhqN9P2AAAAABfj1OBUv359TZkyRVu3btWWLVv0pz/9Sf369dOuXbtK7L9+/XrFxcVp6NCh+vHHH9W/f3/1799fO3furODKr86mlNNKy7xQ6nFDUlrmBW1KOV1xRQEAAAAwZTEMw6WGN2rUqKGXX35ZQ4cOLXbs/vvv1/nz5/X1119b2zp37qy2bdvq7bfftuv6WVlZCggIUGZmpvz9/cusbnt8uf24Ri3cbtpv+l/bql/beuVfEAAAAFCJOZINXOYZp/z8fC1cuFDnz59XVFRUiX02bNigmJgYm7ZevXppw4YNFVHiNavt51Om/QAAAABUDA9nF5CcnKyoqChduHBB1apV0xdffKHIyMgS+6anp6tOnTo2bXXq1FF6enqp18/NzVVu7h8LMWRlZZVN4VehU1gNhQT4KD3zQonPOVkkBQf4qFNYjYouDQAAAMAVOH3EKSIiQtu3b9cPP/ygxx57TA8++KB2795dZtdPSEhQQECAdQsNDS2zazvK3c2i8X0LQ6HlsmNF++P7Rsrd7fKjAAAAAJzJ6cHJy8tL4eHh6tChgxISEtSmTRtNnz69xL7BwcE6efKkTdvJkycVHBxc6vXHjRunzMxM63b06NEyrd9RsS1DNOuB9qrt723THhzgo1kPtFdsyxAnVQYAAACgNE4PTpcrKCiwmVp3qaioKK1YscKmbdmyZaU+EyVJ3t7e1uXOizZni20ZouX/6mbdnzvkZq196k+EJgAAAMBFOfUZp3HjxumOO+5QgwYNdO7cOS1YsEBJSUn69ttvJUmDBg1SvXr1lJCQIEkaNWqUunXrpmnTpqlPnz5auHChtmzZonfeeceZb+OqXDodr1NYDabnAQAAAC7MqcEpIyNDgwYNUlpamgICAtS6dWt9++23uv322yVJqampcnP7Y1CsS5cuWrBggZ577jk988wzatKkiRYtWqSWLVs66y0AAAAAqARc7necypszf8fpUjkXf1fkvwtH1nZP7KUqXk5f4BAAAACoVK7L33ECAAAAAFdFcAIAAAAAEwQnAAAAADBBcAIAAAAAEwQnZ8nP/+P16jW2+wAAAABcCsHJGRITpeaRf+z3vkNq1KiwHQAAAIDLIThVtMRE6Z57pBPHbduPHy9sJzwBAAAALofgVJHy86VRo6SSfjqrqG30aKbtAQAAAC6G4FSR1qyRjh0r/bhhSEePFvYDAAAA4DIIThUpLa1s+wEAAACoEASnihQSUrb9AAAAAFQIglNF6tpVql9fslhKPm6xSKGhhf0AAAAAuAyCU0Vyd5emTy98fXl4Ktp//fXCfgAAAABcBsGpog0YIH32mRRS17a9fv3C9gEDnFMXAAAAgFIRnJxhwABpz+4/9pd8I6WkEJoAAAAAF0VwcpZLp+Pd1pXpeQAAAIALIzgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYIDgBAAAAgAmCEwAAAACYcGpwSkhI0M033yw/Pz/Vrl1b/fv31759+654zty5c2WxWGw2Hx+fCqoYAAAAQGXk1OC0atUqDR8+XBs3btSyZcuUl5ennj176vz581c8z9/fX2lpadbtyJEjFVQxAAAAgMrIw5k3X7p0qc3+3LlzVbt2bW3dulW33XZbqedZLBYFBweXd3kAAAAAIMnFnnHKzMyUJNWoUeOK/bKzs9WwYUOFhoaqX79+2rVrV6l9c3NzlZWVZbMBAAAAgCNcJjgVFBRo9OjRio6OVsuWLUvtFxERoffff19ffvml5s2bp4KCAnXp0kXHjh0rsX9CQoICAgKsW2hoaHm9BQAAAAA3KIthGIazi5Ckxx57TN98843Wrl2r+vXr231eXl6emjdvrri4OE2aNKnY8dzcXOXm5lr3s7KyFBoaqszMTPn7+5dJ7Vcj5+Lvivz3t5Kk3RN7qYqXU2dNAgAAAJVOVlaWAgIC7MoGLvG39REjRujrr7/W6tWrHQpNkuTp6al27drpwIEDJR739vaWt7d3WZQJAAAAoJJy6lQ9wzA0YsQIffHFF/r+++8VFhbm8DXy8/OVnJyskJCQcqgQAAAAAJw84jR8+HAtWLBAX375pfz8/JSeni5JCggIkK+vryRp0KBBqlevnhISEiRJEydOVOfOnRUeHq6zZ8/q5Zdf1pEjR/Twww877X0AAAAAuLE5NTjNmjVLktS9e3eb9jlz5mjw4MGSpNTUVLm5/TEwdubMGQ0bNkzp6emqXr26OnTooPXr1ysyMrKiygYAAABQybjM4hAVxZEHwMoTi0MAAAAAzuVINnCZ5cgBAAAAwFURnAAAAADABPPDnKSKl4cOT+nj7DIAAAAA2IERJwAAAAAwQXACAAAAABMEJwAAAAAwQXACAAAAABMEJwAAAAAwQXACAAAAABMEJwAAAAAwQXACAAAAABMEJwAAAAAwQXACAAAAABMEJwAAAAAwQXACAAAAABMEJwAAAAAwQXACAAAAABMEJwAAAAAwQXACAAAAABMEJwAAAAAwQXACAAAAABMEJwAAAAAwQXACAAAAABMEJwAAAAAwQXACAAAAABMEJwAAAAAwQXACAAAAABMEJwAAAAAw4eHsAiqaYRiSpKysLCdXAgAAAMCZijJBUUa4kkoXnM6dOydJCg0NdXIlAAAAAFzBuXPnFBAQcMU+FsOeeHUDKSgo0IkTJ+Tn5yeLxeLUWrKyshQaGqqjR4/K39/fqbUAQGXDdzAAOIcrff8ahqFz586pbt26cnO78lNMlW7Eyc3NTfXr13d2GTb8/f2d/ocGACorvoMBwDlc5fvXbKSpCItDAAAAAIAJghMAAAAAmCA4OZG3t7fGjx8vb29vZ5cCAJUO38EA4BzX6/dvpVscAgAAAAAcxYgTAAAAAJggOAEAAACACYITAAAAAJggOAEAAACACYLTFaSnp2vUqFEKDw+Xj4+P6tSpo+joaM2aNUs5OTnOLq9U77zzjrp37y5/f39ZLBadPXvW2SUBgEOux+/f06dP6//9v/+niIgI+fr6qkGDBho5cqQyMzOdXRoA2O16/P6VpEceeUQ33XSTfH19FRQUpH79+mnv3r1leg+PMr3aDeTQoUOKjo5WYGCgJk+erFatWsnb21vJycl65513VK9ePd15550lnpuXlydPT88KrvgPOTk5io2NVWxsrMaNG+e0OgDgalyv378nTpzQiRMn9MorrygyMlJHjhzRo48+qhMnTuizzz5zSk0A4Ijr9ftXkjp06KCBAweqQYMGOn36tOLj49WzZ0+lpKTI3d29bG5ioES9evUy6tevb2RnZ5d4vKCgwPpakjFz5kyjb9++RpUqVYzx48cbhmEYM2fONBo3bmx4enoaTZs2NT788EPrOSkpKYYk48cff7S2nTlzxpBkrFy50jAMw1i5cqUhyfj666+NVq1aGd7e3sYtt9xiJCcn2/Ueis4/c+aMQ+8dAJzpRvj+LfLJJ58YXl5eRl5enkPnAYAz3Ejfvz/99JMhyThw4IBD510JU/VKcOrUKX333XcaPny4qlatWmIfi8Visx8fH6+77rpLycnJeuihh/TFF19o1KhRGjNmjHbu3KlHHnlEQ4YM0cqVKx2u54knntC0adO0efNmBQUFqW/fvsrLy7uq9wYAruxG+/7NzMyUv7+/PDyY4AHAtd1I37/nz5/XnDlzFBYWptDQUIfvXaoyi2A3kI0bNxqSjMTERJv2mjVrGlWrVjWqVq1qPPnkk9Z2Scbo0aNt+nbp0sUYNmyYTdu9995r9O7d2zAMxxL3woULrX1OnTpl+Pr6Gh9//LHp+2DECcD15kb5/jUMw/jll1+MBg0aGM8884xd/QHAmW6E79+33nrLqFq1qiHJiIiIKNPRJsNgxMkhmzZt0vbt29WiRQvl5ubaHOvYsaPN/p49exQdHW3TFh0drT179jh836ioKOvrGjVqKCIi4qquAwDXq+vt+zcrK0t9+vRRZGSk4uPjHb4vALiK6+n7d+DAgfrxxx+1atUqNW3aVPfdd58uXLjg8L1Lw9yBEoSHh8tisWjfvn027Y0bN5Yk+fr6FjuntCHN0ri5FWZWwzCsbUy/A1DZ3Qjfv+fOnVNsbKz8/Pz0xRdfOPVhaQCw143w/RsQEKCAgAA1adJEnTt3VvXq1fXFF18oLi6uTK7PiFMJatasqdtvv10zZszQ+fPnr+oazZs317p162za1q1bp8jISElSUFCQJCktLc16fPv27SVea+PGjdbXZ86c0c8//6zmzZtfVV0A4Mqu9+/frKws9ezZU15eXvrf//4nHx+fq3oPAFDRrvfv38sZhiHDMIqNkl0LRpxKMXPmTEVHR6tjx46Kj49X69at5ebmps2bN2vv3r3q0KHDFc9/4okndN9996ldu3aKiYnRV199pcTERC1fvlxSYWrv3LmzpkyZorCwMGVkZOi5554r8VoTJ05UzZo1VadOHT377LOqVauW+vfvX+q909PTlZ6ergMHDkiSkpOT5efnpwYNGqhGjRpX94EAQAW5Xr9/i0JTTk6O5s2bp6ysLGVlZUkq/MtCmS2HCwDl5Hr9/j106JA+/vhj9ezZU0FBQTp27JimTJkiX19f9e7d+5o+Extl+sTUDebEiRPGiBEjjLCwMMPT09OoVq2a0alTJ+Pll182zp8/b+0nyfjiiy+KnX+l5RgNwzB2795tREVFGb6+vkbbtm2N7777rsSH47766iujRYsWhpeXl9GpUyfjp59+umLd48ePNyQV2+bMmXOtHwkAVIjr8fu36JyStpSUlLL4WACg3F2P37/Hjx837rjjDqN27dqGp6enUb9+feNvf/ubsXfv3jL5TIpYDOOSSYZwKUlJSerRo4fOnDmjwMBAZ5cDAJUG378A4Byu/P3LM04AAAAAYILgBAAAAAAmmKoHAAAAACYYcQIAAAAAEwQnAAAAADBBcAIAAAAAEwQnAAAAADBBcAIAAAAAEwQnAIDLS09P16hRoxQeHi4fHx/VqVNH0dHRmjVrlnJycpxdHgCgEvBwdgEAAFzJoUOHFB0drcDAQE2ePFmtWrWSt7e3kpOT9c4776hevXq68847i52Xl5cnT09PJ1QMALgRMeIEAHBpjz/+uDw8PLRlyxbdd999at68uRo3bqx+/fpp8eLF6tu3ryTJYrFo1qxZuvPOO1W1alW9+OKLkqRZs2bppptukpeXlyIiIvTf//7Xeu3Dhw/LYrFo+/bt1razZ8/KYrEoKSlJkpSUlCSLxaLFixerdevW8vHxUefOnbVz507rOUeOHFHfvn1VvXp1Va1aVS1atNCSJUvK/8MBAFQYghMAwGWdOnVK3333nYYPH66qVauW2MdisVhfx8fH66677lJycrIeeughffHFFxo1apTGjBmjnTt36pFHHtGQIUO0cuVKh2t54oknNG3aNG3evFlBQUHq27ev8vLyJEnDhw9Xbm6uVq9ereTkZE2dOlXVqlW7ujcNAHBJTNUDALisAwcOyDAMRURE2LTXqlVLFy5ckFQYWqZOnSpJ+tvf/qYhQ4ZY+8XFxWnw4MF6/PHHJUn/+te/tHHjRr3yyivq0aOHQ7WMHz9et99+uyTpgw8+UP369fXFF1/ovvvuU2pqqu6++261atVKktS4ceOre8MAAJfFiBMA4LqzadMmbd++XS1atFBubq61vWPHjjb99uzZo+joaJu26Oho7dmzx+F7RkVFWV/XqFFDERER1uuMHDlSL7zwgqKjozV+/Hjt2LHD4esDAFwbwQkA4LLCw8NlsVi0b98+m/bGjRsrPDxcvr6+Nu2lTecrjZtb4X8GDcOwthVNv3PEww8/rEOHDunvf/+7kpOT1bFjR7355psOXwcA4LoITgAAl1WzZk3dfvvtmjFjhs6fP+/w+c2bN9e6dets2tatW6fIyEhJUlBQkCQpLS3NevzShSIutXHjRuvrM2fO6Oeff1bz5s2tbaGhoXr00UeVmJioMWPG6N1333W4XgCA6+IZJwCAS5s5c6aio6PVsWNHxcfHq3Xr1nJzc9PmzZu1d+9edejQodRzn3jiCd13331q166dYmJi9NVXXykxMVHLly+XJPn6+qpz586aMmWKwsLClJGRoeeee67Ea02cOFE1a9ZUnTp19Oyzz6pWrVrq37+/JGn06NG644471LRpU505c0YrV660CVUAgOsfwQkA4NJuuukm/fjjj5o8ebLGjRunY8eOydvbW5GRkRo7dqx14YeS9O/fX9OnT9crr7yiUaNGKSwsTHPmzFH37t2tfd5//30NHTpUHTp0UEREhF566SX17Nmz2LWmTJmiUaNGaf/+/Wrbtq2++uoreXl5SZLy8/M1fPhwHTt2TP7+/oqNjdVrr71W5p8FAMB5LMalE7sBAICNpKQk9ejRQ2fOnFFgYKCzywEAOAnPOAEAAACACYITAAAAAJhgqh4AAAAAmGDECQAAAABMEJwAAAAAwATBCQAAAABMEJwAAAAAwATBCQAAAABMEJwAAAAAwATBCQAAAABMEJwAAAAAwATBCQAAAABM/H/my7cAwunEkwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1000x750 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = 10., 7.5\n",
    "\n",
    "errors = np.full((2, ci.shape[0]), np.nan)\n",
    "errors[0, :] = ci['effect'] - ci['2.5 %']\n",
    "errors[1, :] = ci['97.5 %'] - ci['effect']\n",
    "\n",
    "plt.errorbar(ci.index, ci.effect, fmt='o', yerr=errors, label='Estimated Effect (with joint CI)')\n",
    "\n",
    "#add true effect\n",
    "ax = plt.subplot(1, 1, 1)\n",
    "ax.scatter(x=['Group 1', 'Group 2', 'Group 3'], y=true_effects, c='red', label='True Effect')\n",
    "\n",
    "plt.title('GATEs')\n",
    "plt.xlabel('Groups')\n",
    "plt.legend()\n",
    "_ = plt.ylabel('Effect and 95%-CI')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also possible to supply disjoint groups as a single vector (still as a data frame). Remark the slightly different name."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  Group\n",
      "0     3\n",
      "1     2\n",
      "2     2\n",
      "3     2\n",
      "4     2\n"
     ]
    }
   ],
   "source": [
    "groups =  pd.DataFrame(columns=['Group'], index=range(data['X_0'].shape[0]), dtype=str)\n",
    "for i, x_i in enumerate(data['X_0']):\n",
    "    if x_i <= 0.3:\n",
    "         groups['Group'][i] = '1'\n",
    "    elif (x_i > 0.3) & (x_i <= 0.7):\n",
    "         groups['Group'][i] = '2'\n",
    "    else:\n",
    "         groups['Group'][i] = '3'\n",
    "\n",
    "print(groups.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This time lets consider pointwise confidence intervals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            2.5 %    effect    97.5 %\n",
      "Group_1  2.701265  3.016315  3.331365\n",
      "Group_2  5.096550  5.314651  5.532751\n",
      "Group_3  4.412004  4.668981  4.925957\n"
     ]
    }
   ],
   "source": [
    "gate = dml_irm.gate(groups=groups)\n",
    "ci = gate.confint()\n",
    "print(ci)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The coefficients of the best linear predictor can be seen via the summary (the values can be accessed through the underlying model ``.blp_model``)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "             coef   std err          t          P>|t|    [0.025    0.975]\n",
      "Group_1  3.016315  0.144241  20.911662   3.948154e-70  2.732918  3.299712\n",
      "Group_2  5.314651  0.117072  45.396300  6.151047e-179  5.084633  5.544669\n",
      "Group_3  4.668981  0.138851  33.625766  3.991444e-130  4.396173  4.941788\n"
     ]
    }
   ],
   "source": [
    "print(gate.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remark that the confidence intervals in the summary are slightly smaller, since they are not based on the White's heteroskedasticity robus standard errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA04AAAKXCAYAAAC15hrSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA9RklEQVR4nO3de5TVdb3/8deGgRlFZxBFQCWEQC4qqZiJ6FHPwVTMJLty9KB4yUqXWGZF2c/byUvlMfulZPUr6wjZsbCW91uiJZqQWuA1DQSUy0lxBrygwv79wXJyBPzO4Ax7YB6PtfZasz/f797zHlZru559L7tULpfLAQAAYL06VXoAAACA9k44AQAAFBBOAAAABYQTAABAAeEEAABQQDgBAAAUEE4AAAAFhBMAAEAB4QQAAFBAOAEAABQQTgC0a3Pnzs1pp52WXXbZJVtuuWW23HLLDBs2LKeeemr++te/rvM1X/nKV1IqlfLpT3+6yXqpVGrWY/r06Zk3b9677nPxxRdvjD8fgHaiVC6Xy5UeAgDW5cYbb8ynP/3pVFVV5ZhjjskHPvCBdOrUKU888USmTZuWZ599NnPnzk2/fv0aX1Mul/O+970vVVVVWbJkSZYsWZKtt946SXLNNdc0ef9f/OIXueOOO/Lf//3fTdYPOeSQvPrqq+nfv3/GjRuXMWPGrDXbnnvumV133bUN/moA2qOqSg8AAOvyzDPP5DOf+Uz69euXu+66K3369Gmy/ZJLLsmVV16ZTp2anjwxffr0LFy4ML///e9z6KGHZtq0aTnuuOOSJMcee2yTfR944IHccccda60nybx585Ike+211zq3A9CxOFUPgHbp29/+dl5++eX87Gc/WyuakqSqqiqnn356+vbt22R9ypQpGTZsWA4++OCMHj06U6ZMafNZZ82alUMPPTTbbbddtthii/Tv3z8nnHBCm/9eADYeR5wAaJduvPHGDBw4MB/60Iea/ZqVK1fmN7/5Tc4888wkybhx4zJhwoQsXrw4vXv33qA5XnnllfzjH/9Ya7179+6pqqrK0qVL8+EPfzg9e/bM1772tXTv3j3z5s3LtGnTNuj3AdA+OeIEQLvT0NCQ559/Prvtttta21566aX84x//aHy8+uqrjdtuvPHGvPTSS/nMZz6TJBk7dmy6dOmSa6+9doNnOeecc9KzZ8+1HrNmzUqSzJgxI8uWLcuUKVPy5S9/OSeddFL+8z//M4899tgG/04A2h/hBEC709DQkCTZaqut1tp20EEHNQmYK664onHblClTsvfee2fgwIFJkq233jpHHHHEezpd77Of/WzuuOOOtR7Dhg1LsubIU7Im2t54440N/j0AtG9O1QOg3XnrLngrVqxYa9tVV12V5cuXZ8mSJU1u2vDSSy/l5ptvzmmnnZann366cX3UqFH5zW9+k6eeeiq77LJLi2cZNGhQRo8evd7tBx54YD7+8Y/nvPPOy2WXXZaDDjooY8eOzb//+7+nurq6xb8PgPbJEScA2p26urr06dMnc+bMWWvbhz70oYwePTqjRo1qsn7ddddl5cqVufTSSzNo0KDGx5e+9KUkabObRJRKpfz617/O/fffn9NOOy3PPfdcTjjhhIwYMWKd4QfApkk4AdAuHXHEEXn66afz4IMPNmv/KVOmZLfddst111231mP06NGZOnVqm86777775lvf+lZmzZqVKVOm5NFHH31P11YB0L44VQ+AdukrX/lKpk6dmhNOOCF33XVXevXq1WT727+/fcGCBbn33ntz3nnn5ROf+MRa7/X666/nmGOOyZ/+9KcW3aWvOZYtW5bu3bunVCo1ru2xxx5J1tzlD4DNg3ACoF0aNGhQpk6dmnHjxmXw4ME55phj8oEPfCDlcjlz587N1KlT06lTp+y0006ZOnVqyuVyPvrRj67zvcaMGZOqqqpMmTKlxeH00EMP5Zprrllr/f3vf39GjhyZn//857nyyivzsY99LO9///uzfPny/PjHP05tbW3GjBmzQX87AO1Pqfz2/8sOANqZZ555JpdeemnuuOOOLFy4MKVSKf369ctBBx2Uz33uc/nABz6Q4cOHp76+Ps8+++x63+fggw/OY489lueeey5VVWv+f8PTTjstV1xxRdb1n8J58+alf//+632/4447LldffXUefvjhfOc738l9992XJUuWpK6uLvvss0/OPffcjBgx4r3/AwDQLggnAACAAm4OAQAAUEA4AQAAFBBOAAAABYQTAABAAeEEAABQQDgBAAAU6HBfgLt69eo8//zz2XrrrZt8yzsAANCxlMvlLF++PDvssEM6dXr3Y0odLpyef/759O3bt9JjAAAA7cSCBQuy0047ves+HS6ctt566yRr/nFqa2srPA0AAFApDQ0N6du3b2MjvJsOF05vnZ5XW1srnAAAgGZdwuPmEAAAAAWEEwAAQAHhBAAAUEA4AQAAFBBOAAAABYQTAABAAeEEAABQQDgBAAAUEE4AAAAFhBMAAEAB4QQAAFBAOAEAABQQTgAAAAWEEwAAQAHhBAAAUEA4AQAAFBBOAAAABYQTAABAAeEEAABQQDgBAAAUEE4AAAAFhBMAAEAB4QQAAFBAOAHQ4bzy+pvZ+Ws3Zeev3ZRXXn+z0uMAsAkQTgAAAAUqGk7nnntuSqVSk8eQIUPWu//VV1+91v41NTUbcWIAAKAjqqr0ALvuumvuvPPOxudVVe8+Um1tbZ588snG56VSqc1mAwAASNpBOFVVVaV3797N3r9UKrVofwAAgPeq4tc4/e1vf8sOO+yQAQMG5Jhjjsn8+fPfdf8VK1akX79+6du3b4466qg8+uijG2lSAACgo6poOH3oQx/K1VdfnVtvvTWTJ0/O3Llzc8ABB2T58uXr3H/w4MH56U9/mt/97ne55pprsnr16uy3335ZuHDhen/HypUr09DQ0OQBAADQEqVyuVyu9BBveemll9KvX7/813/9V0488cTC/d94440MHTo048aNywUXXLDOfc4999ycd955a63X19entrb2Pc8MwKbnldffzLD/c1uS5LHzD82WXSt+5joAFdDQ0JC6urpmtUHFT9V7u+7du2eXXXbJ008/3az9u3Tpkj333PNd9580aVLq6+sbHwsWLGitcQEAgA6iXYXTihUr8swzz6RPnz7N2n/VqlWZPXv2u+5fXV2d2traJg8AAICWqGg4ffnLX84999yTefPmZcaMGfnYxz6Wzp07Z9y4cUmS8ePHZ9KkSY37n3/++bn99tvz97//PQ899FCOPfbYPPvssznppJMq9ScAAAAdQEVP6l64cGHGjRuXF154IT179sz++++fBx54ID179kySzJ8/P506/bPtli1blpNPPjmLFy/ONttskxEjRmTGjBkZNmxYpf4EAACgA2hXN4fYGFpyARgAmyc3hwAg2YRvDgEAANAeCScAAIACwgkAAKCAcAIAACggnAAAAAoIJwAAgALCCQAAoIBwAgAAKCCcAAAACggnAACAAsIJAACggHACAAAoIJwAAAAKCCcAAIACwgkAAKCAcAIAACggnAAAAAoIJwAAgALCCQAAoIBwAgAAKCCcAAAACggnAACAAsIJAACggHACoMNZtbrc+PODc19s8hwA1kU4AdCh3DpnUUb/1z2Nz4//2czsf8nvc+ucRRWcCoD2TjgB0GHcOmdRPn/NQ1nSsLLJ+uL61/L5ax4STwCsl3ACoENYtbqc8254LOs6Ke+ttfNueMxpewCsk3ACoEN4cO6LWVT/2nq3l5Msqn8tD859ceMNBcAmQzgB0CEsXb7+aNqQ/QDoWIQTAB3C9lvXtOp+AHQswgmADmGf/j3Sp64mpfVsLyXpU1eTffr32JhjAbCJEE4AdAidO5VyzpHDkmSteHrr+TlHDkvnTutLKwA6MuEEQIdx2G59MvnYvbJ9bXWT9d51NZl87F45bLc+FZoMgPauqtIDAMDGdNhufTJq4HbZ/dzbkyRXT/hgDhjU05EmAN6VI04AdDhvj6R9+vcQTQAUEk4AAAAFhBMAAEAB4QQAAFBAOAEAABQQTgAAAAWEEwAdz6pV//z53j80fQ4A6yCcAOhYpk1Lhg775/Mxhyc777xmHQDWQzgB0HFMm5Z84hPJ8881XX/uuTXr4gmA9RBOAHQMq1YlEycm5fLa295aO+MMp+0BsE7CCYCO4Q9/SBYuXP/2cjlZsGDNfgDwDsIJgI5h0aLW3Q+ADkU4AdAx9OnTuvsB0KEIJwA6hgMOSHbaKSmV1r29VEr69l2zHwC8g3ACoGPo3Dm5/PI1P78znt56/r3vrdkPAN5BOAHQcRx9dPLrXyd9dmi6vtNOa9aPProycwHQ7gknADqWo49OHn/sn89vviWZO1c0AfCuhBMAHc/bT8f7lwOcngdAIeEEAABQQDgBAAAUEE4AAAAFhBMAAEAB4QQAAFBAOAEAABQQTgAAAAWEEwAAQAHhBAAAUEA4AQAAFBBOAAAABYQTAABAAeEEAABQQDgBAAAUEE4AAAAFhBMAAECBqkoPAAAb25ZdqzLv4iMqPQYAmxBHnAAAAAoIJwAAgALCCQAAoIBwAgAAKCCcAAAACggnAACAAsIJAACggHACAAAoIJwAAAAKCCcAAIACwgkAAKCAcAIAACggnAAAAApUNJzOPffclEqlJo8hQ4a862uuu+66DBkyJDU1Ndl9991z8803b6RpAQCAjqriR5x23XXXLFq0qPHxxz/+cb37zpgxI+PGjcuJJ56Yhx9+OGPHjs3YsWMzZ86cjTgxAADQ0VQ8nKqqqtK7d+/Gx3bbbbfefS+//PIcdthhOeusszJ06NBccMEF2WuvvfKDH/xgI04MAAB0NBUPp7/97W/ZYYcdMmDAgBxzzDGZP3/+eve9//77M3r06CZrhx56aO6///71vmblypVpaGho8gAAAGiJiobThz70oVx99dW59dZbM3ny5MydOzcHHHBAli9fvs79Fy9enF69ejVZ69WrVxYvXrze33HRRRelrq6u8dG3b99W/RsAAIDNX0XD6fDDD88nP/nJDB8+PIceemhuvvnmvPTSS/mf//mfVvsdkyZNSn19feNjwYIFrfbeAABAx1BV6QHernv37tlll13y9NNPr3N77969s2TJkiZrS5YsSe/evdf7ntXV1amurm7VOQEAgI6l4tc4vd2KFSvyzDPPpE+fPuvcPnLkyNx1111N1u64446MHDlyY4wHAAB0UBUNpy9/+cu55557Mm/evMyYMSMf+9jH0rlz54wbNy5JMn78+EyaNKlx/4kTJ+bWW2/NpZdemieeeCLnnntuZs2aldNOO61SfwIAANABVPRUvYULF2bcuHF54YUX0rNnz+y///554IEH0rNnzyTJ/Pnz06nTP9tuv/32y9SpU3P22Wfn61//egYNGpTf/va32W233Sr1JwAAAB1AqVwulys9xMbU0NCQurq61NfXp7a2ttLjAAAAFdKSNmhX1zgBAAC0R8IJAACggHACAAAoIJwAAAAKCCcAAIACwgkAAKCAcAIAACggnAAAgI3ildffzM5fuyk7f+2mvPL6m5Uep0WEEwAAQAHhBAAAUEA4AQAAFBBOAAAABYQTAABAAeEEAABQQDgBAAAUEE4AAAAFhBMAAEAB4QQAAFBAOAEAABQQTgAAAAWEEwAAsHGsWvXPn+/9Q9Pn7ZxwAgAA2t60acnQYf98PubwZOed16xvAoQTAADQtqZNSz7xieT555quP/fcmvVNIJ6EEwAA0HZWrUomTkzK5bW3vbV2xhnt/rQ94QQAALSdP/whWbhw/dvL5WTBgjX7tWPCCQAAaDuLFrXufhUinAAAgLbTp0/r7lchwgkAAGg7BxyQ7LRTUiqte3uplPTtu2a/dkw4AQAAbadz5+Tyy9f8/M54euv59763Zr92TDgBAABt6+ijk1//OumzQ9P1nXZas3700ZWZqwWEEwAA0PaOPjp5/LF/Pr/5lmTu3E0imhLhBAAAbCxvPx3vXw5o96fnvZ1wAgAAKCCcAAAACggnAACAAsIJAACggHACAAAoIJwAAAAKCCcAAIACwgkAAKCAcAIAACggnAAAAAoIJwAAgALCCQAAoIBwAgAAKCCcAAAACggnAACAAsIJAACggHACAAAoIJwAAAAKCCcAAGCjWLW63Pjzg3NfbPK8vRNOAABAm7t1zqKM/q97Gp8f/7OZ2f+S3+fWOYsqOFXzCScAAKBN3TpnUT5/zUNZ0rCyyfri+tfy+Wse2iTiSTgBAABtZtXqcs674bGs66S8t9bOu+Gxdn/annACAADazINzX8yi+tfWu72cZFH9a3lw7osbb6gNIJwAAIA2s3T5+qNpQ/arFOEEAAC0me23rmnV/SpFOAEAAG1mn/490qeuJqX1bC8l6VNXk33699iYY7WYcAIAANpM506lnHPksCRZK57een7OkcPSudP60qp9EE4AAECbOmy3Ppl87F7Zvra6yXrvuppMPnavHLZbnwpN1nxVlR4AAADY/B22W5+MGrhddj/39iTJ1RM+mAMG9Wz3R5re4ogTAACwUbw9kvbp32OTiaZEOAEAABQSTgAAAAWEEwAAQAHhBAAAUEA4AQAAFBBOAAAABYQTAABAAeEEAABQQDgBAAAUEE4AAAAFhBMAAEAB4QQAAFBAOAEAABQQTgAAAAWEEwAAQAHhBAAAUEA4AQAAFBBOAAAABYQTAABAAeEEAABQoN2E08UXX5xSqZQzzjhjvftcffXVKZVKTR41NTUbb0gAAKBDqqr0AEkyc+bMXHXVVRk+fHjhvrW1tXnyyScbn5dKpbYcDQAAoPJHnFasWJFjjjkmP/7xj7PNNtsU7l8qldK7d+/GR69evTbClAAAQEdW8XA69dRTc8QRR2T06NHN2n/FihXp169f+vbtm6OOOiqPPvrou+6/cuXKNDQ0NHkAAAC0REXD6dprr81DDz2Uiy66qFn7Dx48OD/96U/zu9/9Ltdcc01Wr16d/fbbLwsXLlzvay666KLU1dU1Pvr27dta4wMAAB1ExcJpwYIFmThxYqZMmdLsGzyMHDky48ePzx577JEDDzww06ZNS8+ePXPVVVet9zWTJk1KfX1942PBggWt9ScAAAAdRMVuDvHnP/85S5cuzV577dW4tmrVqtx77735wQ9+kJUrV6Zz587v+h5dunTJnnvumaeffnq9+1RXV6e6urrV5gYAADqeioXTv/3bv2X27NlN1iZMmJAhQ4bkq1/9amE0JWtCa/bs2RkzZkxbjQkAAFC5cNp6662z2267NVnr1q1btt1228b18ePHZ8cdd2y8Bur888/Pvvvum4EDB+all17Kd77znTz77LM56aSTNvr8AABAx9EuvsdpfebPn59Onf55GdayZcty8sknZ/Hixdlmm20yYsSIzJgxI8OGDavglAAAwOauVC6Xy5UeYmNqaGhIXV1d6uvrU1tbW+lxAACgw3jl9Tcz7P/cliR57PxDs2XXyh7HaUkbVPx7nAAAANo74QQAAFBAOAEAABRo9kmF3//+95u13+mnn77BwwAAALRHzQ6nyy67rHCfUqkknAAAgM1Os8Np7ty5bTkHAABAu+UaJwAAgALNDqff//73GTZsWBoaGtbaVl9fn1133TX33ntvqw4HAADQHjQ7nL73ve/l5JNPXucXQ9XV1eWUU05p1nVQAAAAm5pSuVwuN2fHfv365dZbb83QoUPXuf2JJ57Ihz/84cyfP79VB2xtLfl2YAAAYPPVkjZo9hGnJUuWpEuXLuvdXlVVlf/93/9t/pQAAACbiGaH04477pg5c+asd/tf//rX9OnTp1WGAgAAaE+aHU5jxozJN7/5zbz22mtrbXv11Vdzzjnn5CMf+UirDgcAANAeNPsapyVLlmSvvfZK586dc9ppp2Xw4MFJ1lzbdMUVV2TVqlV56KGH0qtXrzYd+L1yjRMAAJC0rA2a/QW4vXr1yowZM/L5z38+kyZNylu9VSqVcuihh+aKK65o99EEAACwIZodTsmaO+vdfPPNWbZsWZ5++umUy+UMGjQo22yzTVvNBwAAUHHNvsbp7bbZZpt88IMfzDPPPJOuXbu29kwAAADtygaF01tOOeWULFmypLVmAQAAaJfeUzg1874SAAAAm7T3FE4AAAAdwXsKp1tuuSU77rhja80CAADQLrXornrvNGTIkNx5551ZtWpVPvjBD6ZPnz6tNRcAAEC7scHh9Jvf/CYnnnhidtlll7zxxht58sknc8UVV2TChAmtOR8AAEDFNftUvRUrVjR5ft555+XBBx/Mgw8+mIcffjjXXXddvvGNb7T6gAAAAJXW7HAaMWJEfve73zU+r6qqytKlSxufL1myxHc6AQAAm6VSuZn3FJ83b15OPfXUdO3aNVdccUWeeeaZfOYzn8mqVavy5ptvplOnTrn66qszZsyYtp75PWloaEhdXV3q6+tTW1tb6XEAAIAKaUkbNPsap5133jk33XRTfvnLX+bAAw/M6aefnqeffjpPP/10Vq1alSFDhqSmpuY9Dw8AANDetPh25OPGjcvMmTPzl7/8JQcddFBWr16dPfbYQzQBAACbrRbdVe/mm2/O448/ng984AP5yU9+knvuuSfHHHNMDj/88Jx//vnZYost2mpOAACAimn2EaczzzwzEyZMyMyZM3PKKafkggsuyIEHHpiHHnooNTU12XPPPXPLLbe05awAAAAV0eybQ2y77ba5/fbbM2LEiLz44ovZd99989RTTzVuf+yxx3LKKafkD3/4Q5sN2xrcHAIAAEha1gbNPuLUrVu3zJ07N0myYMGCta5pGjZsWLuPJgAAgA3R7HC66KKLMn78+Oywww458MADc8EFF7TlXAAAAO1Gs0/VS5IXXnghf//73zNo0KB07969DcdqO07VAwAAkjb6HqdkzXVO22677XsaDgAAYFPT7FP1li5d2uT5I488kuOOOy6jRo3KJz7xiUyfPr21ZwMAAGgXmh1Offr0aYynGTNmZJ999smzzz6bUaNGpaGhIYccckjuvffeNhsUAACgUpp9jVOnTp2yePHibL/99vnwhz+cvn375v/9v//XuP2MM87I7Nmzc9ddd7XZsK3BNU4AAEDSRrcjf7s5c+bk5JNPbrJ28skn569//euGvB0AAEC71qKbQyxfvjw1NTWpqalJdXV1k201NTV55ZVXWnU4AACA9qBFR5x22WWXbLPNNpk3b15mzZrVZNujjz6aHXbYoVWHAwAAaA+afcTp7rvvbvK8T58+TZ7PnTs3n/3sZ1tnKgAAgHakRV+AuzlwcwgAACBpwy/Afcv8+fOzaNGidOrUKQMGDPCluAAAwGatRdc4XXnllenXr1/69++f/fbbL/vuu2+233777L///vnzn//cVjMCAABUVLPD6bvf/W6+9a1v5ayzzspVV12VwYMH59xzz81NN92UAQMG5F/+5V/WumEEAADA5qDZ1zj1798/V155ZQ4//PAkyVNPPZX99tsvixcvTlVVVSZOnJjHH388t99+e5sO/F65xgkAAEja6Atwly5dmqFDhzY+HzRoUOrr6/O///u/SZITTjgh999//waODAAA0H41O5x22WWX3HHHHY3P77777nTt2jW9e/dOsuYLcEulUutPCAAAUGHNvqvepEmTcuyxx+bOO+9MTU1Npk2bltNPP70xlqZPn57ddtutzQYFAAColBZ9j9Mtt9ySa665JitXrsyhhx6ak08+uXHbCy+8kCTt/tbkrnECAACSlrWBL8AFAAA6pDa5OQQAAEBHJZwAAAAKCCcAAIACwgkAAKCAcAIAACjQrO9xOvroo5v9htOmTdvgYQAAANqjZh1xqqura3zU1tbmrrvuyqxZsxq3//nPf85dd92Vurq6NhsUAACgUpp1xOlnP/tZ489f/epX86lPfSo//OEP07lz5yTJqlWr8oUvfMH3IgEAAJulFn8Bbs+ePfPHP/4xgwcPbrL+5JNPZr/99ssLL7zQqgO2Nl+ACwAAJG38BbhvvvlmnnjiibXWn3jiiaxevbqlbwcAANDuNetUvbebMGFCTjzxxDzzzDPZZ599kiR/+tOfcvHFF2fChAmtPiAAAECltTicvvvd76Z379659NJLs2jRoiRJnz59ctZZZ+XMM89s9QEBAAAqrcXXOL1dQ0NDkmxS1wq5xgkAAEha1gYtPuL0dsIDAADoCFp8c4glS5bkP/7jP7LDDjukqqoqnTt3bvIAAADY3LT4iNPxxx+f+fPn55vf/Gb69OmTUqnUFnMBAAC0Gy0Opz/+8Y/5wx/+kD322KMNxgEAAGh/WnyqXt++ffMe7icBAACwyWlxOH3ve9/L1772tcybN68NxgEAAGh/Wnyq3qc//em88soref/7358tt9wyXbp0abL9xRdfbLXhAAAA2oMWh9P3vve9NhgDAACg/WpxOB133HFtMQcAAEC79Z6+APe1117L66+/3mTNl+ICAACbmxbfHOLll1/Oaaedlu233z7dunXLNtts0+QBAACwuWlxOH3lK1/J73//+0yePDnV1dX5yU9+kvPOOy877LBDfvGLX7TFjAAAABXV4lP1brjhhvziF7/IQQcdlAkTJuSAAw7IwIED069fv0yZMiXHHHNMW8wJAABQMS0+4vTiiy9mwIABSdZcz/TW7cf333//3Hvvva07HQAAQDvQ4nAaMGBA5s6dmyQZMmRI/ud//ifJmiNR3bt3b9XhAAAA2oMWh9OECRPyl7/8JUnyta99LVdccUVqamryxS9+MWeddVarDwgAAFBpLQ6nL37xizn99NOTJKNHj84TTzyRqVOn5uGHH87EiRM3eJCLL744pVIpZ5xxxrvud91112XIkCGpqanJ7rvvnptvvnmDfycAAEBztDic3qlfv345+uijM3z48A1+j5kzZ+aqq64qfI8ZM2Zk3LhxOfHEE/Pwww9n7NixGTt2bObMmbPBvxsAAKDIew6n92rFihU55phj8uMf/7jwe6Auv/zyHHbYYTnrrLMydOjQXHDBBdlrr73ygx/8YCNNCwAAdEQVD6dTTz01RxxxREaPHl247/3337/Wfoceemjuv//+thoPAACg5d/j1JquvfbaPPTQQ5k5c2az9l+8eHF69erVZK1Xr15ZvHjxel+zcuXKrFy5svF5Q0PDhg0LAAB0WBU74rRgwYJMnDgxU6ZMSU1NTZv9nosuuih1dXWNj759+7bZ7wIAADZPLQ6nzp07Z+nSpWutv/DCC+ncuXOz3+fPf/5zli5dmr322itVVVWpqqrKPffck+9///upqqrKqlWr1npN7969s2TJkiZrS5YsSe/evdf7eyZNmpT6+vrGx4IFC5o9IwAAQLIBp+qVy+V1rq9cuTJdu3Zt9vv827/9W2bPnt1kbcKECRkyZEi++tWvrjPCRo4cmbvuuqvJLcvvuOOOjBw5cr2/p7q6OtXV1c2eCwAA4J2aHU7f//73kySlUik/+clPstVWWzVuW7VqVe69994MGTKk2b946623zm677dZkrVu3btl2220b18ePH58dd9wxF110UZJk4sSJOfDAA3PppZfmiCOOyLXXXptZs2blRz/6UbN/LwAAQEs1O5wuu+yyJGuOOP3whz9sckSoa9eu2XnnnfPDH/6wVYebP39+OnX659mE++23X6ZOnZqzzz47X//61zNo0KD89re/XSvAAAAAWlOpvL5z79bj4IMPzrRp0wq/c6m9amhoSF1dXerr61NbW1vpcQAAgAppSRu0+Bqnu+++e4MHAwAA2BS1+K56H//4x3PJJZestf7tb387n/zkJ1tlKAAAgPakxeF07733ZsyYMWutH3744bn33ntbZSgAAID2pMXhtGLFinXedrxLly5paGholaEAAADakxaH0+67755f/epXa61fe+21GTZsWKsMBQAA0J60+OYQ3/zmN3P00UfnmWeeyb/+678mSe6666788pe/zHXXXdfqAwIAAFRai8PpyCOPzG9/+9tceOGF+fWvf50tttgiw4cPz5133pkDDzywLWYEAACoqBZ/j9Omzvc4AQAAScvaoMXXOCXJSy+9lJ/85Cf5+te/nhdffDFJ8tBDD+W5557bkLcDAABo11p8qt5f//rXjB49OnV1dZk3b15OOumk9OjRI9OmTcv8+fPzi1/8oi3mBAAAqJgWH3H60pe+lOOPPz5/+9vfUlNT07g+ZswY3+MEAABsllocTjNnzswpp5yy1vqOO+6YxYsXt8pQAAAA7UmLw6m6unqdX3T71FNPpWfPnq0yFAAAQHvS4nD66Ec/mvPPPz9vvPFGkqRUKmX+/Pn56le/mo9//OOtPiAAAECltTicLr300qxYsSLbb799Xn311Rx44IEZOHBgtt5663zrW99qixkBAAAqqsV31aurq8sdd9yR++67L3/5y1+yYsWK7LXXXhk9enRbzAcAAFBxzQqnHj165Kmnnsp2222XE044IZdffnlGjRqVUaNGtfV8AAAAFdesU/Vef/31xhtC/PznP89rr73WpkMBAAC0J8064jRy5MiMHTs2I0aMSLlczumnn54ttthinfv+9Kc/bdUBAQAAKq1Z4XTNNdfksssuyzPPPJMkqa+vd9QJAADoMErlcrnckhf0798/s2bNyrbbbttWM7WphoaG1NXVpb6+PrW1tZUeBwAAqJCWtEGzrnHq0aNH/vGPfyRJDj744HTt2vW9TwkAALCJcHMIAACAAm4OAQAAUKDFN4colUpuDgEAAHQobg4BAAB0SC1pg2YdcXq7uXPnbvBgAAAAm6Jm3RwiScaMGZP6+vrG5xdffHFeeumlxucvvPBChg0b1qrDAQAAtAfNDqfbbrstK1eubHx+4YUX5sUXX2x8/uabb+bJJ59s3ekAAADagWaH0zsvhWrhpVEAAACbrGaHEwAAQEfV7HAqlUoplUprrQEAAGzumn1XvXK5nOOPPz7V1dVJktdeey2f+9zn0q1btyRpcv0TAADA5qTZ4XTcccc1eX7ssceutc/48ePf+0QAAADtTLPD6Wc/+1lbzgEAANBuuTkEAABAAeEEAABQQDgBAAAUEE4AAAAFhBMAAEAB4QQAAFBAOAEAABQQTgAAAAWEEwAAQAHhBAAAUEA4AQAAFBBOAAAABYQTAABAAeEEAABQQDgBAAAUEE4AAAAFhBMAAEAB4QQAAFBAOAEAABQQTgAAAAWEEwAAQAHhBAAAUEA4AQAAFBBOAAAABYQTAABAAeEEAABQQDgBAAAUEE4AAAAFhBMAAEAB4QQAAFBAOAEAABQQTgAAAAWEEwAAQAHhBAAAUEA4AQAAFBBOAAAABYQTAABAAeEEAABQQDgBAAAUEE4AAAAFhBMAAEAB4QQAAFBAOAEAABQQTgAAAAWEEwAAQAHhBAAAUEA4AQAAFKhoOE2ePDnDhw9PbW1tamtrM3LkyNxyyy3r3f/qq69OqVRq8qipqdmIEwMAAB1RVSV/+U477ZSLL744gwYNSrlczs9//vMcddRRefjhh7Prrruu8zW1tbV58sknG5+XSqWNNS4AANBBVTScjjzyyCbPv/Wtb2Xy5Ml54IEH1htOpVIpvXv33hjjAQAAJGlH1zitWrUq1157bV5++eWMHDlyvfutWLEi/fr1S9++fXPUUUfl0Ucffdf3XblyZRoaGpo8AAAAWqLi4TR79uxstdVWqa6uzuc+97lcf/31GTZs2Dr3HTx4cH7605/md7/7Xa655pqsXr06++23XxYuXLje97/oootSV1fX+Ojbt29b/SkAAMBmqlQul8uVHOD111/P/PnzU19fn1//+tf5yU9+knvuuWe98fR2b7zxRoYOHZpx48blggsuWOc+K1euzMqVKxufNzQ0pG/fvqmvr09tbW2r/R0AAMCmpaGhIXV1dc1qg4pe45QkXbt2zcCBA5MkI0aMyMyZM3P55ZfnqquuKnxtly5dsueee+bpp59e7z7V1dWprq5utXkBAICOp+Kn6r3T6tWrmxwhejerVq3K7Nmz06dPnzaeCgAA6MgqesRp0qRJOfzww/O+970vy5cvz9SpUzN9+vTcdtttSZLx48dnxx13zEUXXZQkOf/887Pvvvtm4MCBeemll/Kd73wnzz77bE466aRK/hkAAMBmrqLhtHTp0owfPz6LFi1KXV1dhg8fnttuuy2HHHJIkmT+/Pnp1OmfB8WWLVuWk08+OYsXL84222yTESNGZMaMGc26HgoAAGBDVfzmEBtbSy4AAwAANl8taYN2d40TAABAeyOcAAAACggnAACAAsIJAACggHACAAAoIJwAAAAKCCcAAIACwgkAAKCAcAIAACggnAAAAAoIJwAAgALCCQAAoIBwAgAAKCCcAAAACggnAACAAsIJAACggHACAAAoIJwAAAAKCCcAAIACwgkAAKCAcAIAACggnAAAAAoIJwAAgALCCQAAoIBwAgAAKCCcKuSV19/Mzl+7KTt/7aa88vqblR4HAAB4F8IJAACggHACAAAoIJwAAAAKCCcAAIACwgkAAKCAcAIAACggnAAAAAoIJwAAgALCCQAAoIBwAgAAKCCcAAAACggnAACAAsIJAACggHACAAAoIJwAAAAKCCcAAIACwgkAAKCAcAIAACggnAAAAAoIJwAAgALCCQAAoIBwAgAAKCCcAAAACggnAACAAsIJAACggHACAAAoIJwAAAAKCCcAAIACwgkAAKCAcAIAACggnAAAAAoIJwAAgALCCQAAoIBwAgAAKCCcAAAACggnAACAAsIJAACggHACAAAoIJwAAAAKCCcAAIACwgkAAKCAcAIAACggnAAAAAoIJwAAgALCCQAAoIBwAgAAKCCcAAAACggnAACAAsIJAACggHACAAAoIJwAAAAKCCcAAIACwqlCVq0uN/784NwXmzwHAADaF+FUAbfOWZTR/3VP4/PjfzYz+1/y+9w6Z1EFpwIAANZHOG1kt85ZlM9f81CWNKxssr64/rV8/pqHxBMAALRDwmkjWrW6nPNueCzrOinvrbXzbnjMaXsAANDOVDScJk+enOHDh6e2tja1tbUZOXJkbrnllnd9zXXXXZchQ4akpqYmu+++e26++eaNNO179+DcF7Oo/rX1bi8nWVT/Wh6c++LGGwoAAChU0XDaaaedcvHFF+fPf/5zZs2alX/913/NUUcdlUcffXSd+8+YMSPjxo3LiSeemIcffjhjx47N2LFjM2fOnI08+YZZunz90bQh+wEAABtHqVwut6vzwnr06JHvfOc7OfHEE9fa9ulPfzovv/xybrzxxsa1fffdN3vssUd++MMfNuv9GxoaUldXl/r6+tTW1rba3M1x/zMvZNyPHyjc75cn75uR7992I0wEAAAdV0vaoN1c47Rq1apce+21efnllzNy5Mh17nP//fdn9OjRTdYOPfTQ3H///et935UrV6ahoaHJo1L26d8jfepqUlrP9lKSPnU12ad/j405FgAAUKDi4TR79uxstdVWqa6uzuc+97lcf/31GTZs2Dr3Xbx4cXr16tVkrVevXlm8ePF63/+iiy5KXV1d46Nv376tOn9LdO5UyjlHrvnb3hlPbz0/58hh6dxpfWkFAABUQsXDafDgwXnkkUfypz/9KZ///Odz3HHH5bHHHmu19580aVLq6+sbHwsWLGi1994Qh+3WJ5OP3Svb11Y3We9dV5PJx+6Vw3brU6HJAACA9amq9ABdu3bNwIEDkyQjRozIzJkzc/nll+eqq65aa9/evXtnyZIlTdaWLFmS3r17r/f9q6urU11dvd7tlXDYbn0yauB22f3c25MkV0/4YA4Y1NORJgAAaKcqfsTpnVavXp2VK1euc9vIkSNz1113NVm744471ntNVHv29kjap38P0QQAAO1YRY84TZo0KYcffnje9773Zfny5Zk6dWqmT5+e2267LUkyfvz47LjjjrnooouSJBMnTsyBBx6YSy+9NEcccUSuvfbazJo1Kz/60Y8q+WcAAACbuYqG09KlSzN+/PgsWrQodXV1GT58eG677bYccsghSZL58+enU6d/HhTbb7/9MnXq1Jx99tn5+te/nkGDBuW3v/1tdtttt0r9CQAAQAfQ7r7Hqa1V8nuc3u6V19/MsP+z5sjaY+cfmi27VvxyMwAA6FA2ye9xAgAAaK+EEwAAQAHhBAAAUEA4AQAAFBBOAAAABYRTpaxa9c+f7/1D0+cAAEC7IpwqYdq0ZOiwfz4fc3iy885r1gEAgHZHOG1s06Yln/hE8vxzTdefe27NungCAIB2RzhtTKtWJRMnJuv6zuG31s44w2l7AADQzginjekPf0gWLlz/9nI5WbBgzX4AAEC7IZw2pkWLWnc/AABgoxBOG1OfPq27HwAAsFEIp43pgAOSnXZKSqV1by+Vkr591+wHAAC0G8JpY+rcObn88jU/vzOe3nr+ve+t2Q8AAGg3hNPGdvTRya9/nfTZoen6TjutWT/66MrMBQAArJdwqoSjj04ef+yfz2++JZk7VzQBAEA7JZwq5e2n4/3LAU7PAwCAdkw4AQAAFBBOAAAABYQTAABAAeEEAABQQDgBAAAUEE4AAAAFhBMAAEAB4QQAAFBAOAEAABQQTgAAAAWEEwAAQAHhBAAAUEA4AQAAFBBOAAAABYQTAABAAeEEAABQQDgBAAAUEE4AAAAFhBMAAEAB4QQAAFBAOAEAABQQTgAAAAWEEwAAQAHhBAAAUEA4AQAAFBBOAAAABYQTAABAAeEEAABQQDgBAAAUEE4AAAAFhBMAAECBqkoP0FFt2bUq8y4+otJjAAAAzeCIEwAAQAHhBAAAUEA4AQAAFBBOAAAABYQTAABAAeEEAABQQDgBAAAUEE4AAAAFhBMAAEAB4QQAAFBAOAEAABQQTgAAAAWEEwAAQAHhBAAAUEA4AQAAFBBOAAAABYQTAABAAeEEAABQQDgBAAAUEE4AAAAFhBMAAEAB4QQAAFBAOAEAABQQTgAAAAWEEwAAQIGqSg+wsZXL5SRJQ0NDhScBAAAq6a0meKsR3k2HC6fly5cnSfr27VvhSQAAgPZg+fLlqaure9d9SuXm5NVmZPXq1Xn++eez9dZbp1QqVXSWhoaG9O3bNwsWLEhtbW1FZwHoaHwGA1RGe/r8LZfLWb58eXbYYYd06vTuVzF1uCNOnTp1yk477VTpMZqora2t+P9oADoqn8EAldFePn+LjjS9xc0hAAAACggnAACAAsKpgqqrq3POOeekurq60qMAdDg+gwEqY1P9/O1wN4cAAABoKUecAAAACggnAACAAsIJAACggHACAAAoIJzWY/HixZk4cWIGDhyYmpqa9OrVK6NGjcrkyZPzyiuvVHq8d/Xaa6/l+OOPz+67756qqqqMHTu20iMBNNum/Pk7ffr0HHXUUenTp0+6deuWPfbYI1OmTKn0WADNtil/Bj/55JM5+OCD06tXr9TU1GTAgAE5++yz88Ybb7TK+1e1yrtsZv7+979n1KhR6d69ey688MLsvvvuqa6uzuzZs/OjH/0oO+64Yz760Y+u9bo33ngjXbp0qcDETa1atSpbbLFFTj/99PzmN7+p9DgAzbapf/7OmDEjw4cPz1e/+tX06tUrN954Y8aPH5+6urp85CMfqfR4AO9qU/8M7tKlS8aPH5+99tor3bt3z1/+8pecfPLJWb16dS688ML3/gvKrOXQQw8t77TTTuUVK1asc/vq1avL5XK5nKR85ZVXlo888sjylltuWT7nnHPK5XK5fOWVV5YHDBhQ7tKlS3mXXXYp/+IXv2h87dy5c8tJyg8//HDj2rJly8pJynfffXe5XC6X77777nKS8o033ljefffdy9XV1eUPfehD5dmzZ7f4bznuuOPKRx11VItfB1AJm9Pn71vGjBlTnjBhwga/HmBj2Rw/g7/4xS+W999//w1+/ds5Ve8dXnjhhdx+++059dRT061bt3XuUyqVGn8+99xz87GPfSyzZ8/OCSeckOuvvz4TJ07MmWeemTlz5uSUU07JhAkTcvfdd7d4lrPOOiuXXnppZs6cmZ49e+bII49stUONAO3N5vr5W19fnx49emzQawE2ls3xM/jpp5/OrbfemgMPPLDFr12nVsmvzcgDDzxQTlKeNm1ak/Vtt9223K1bt3K3bt3KX/nKV8rl8praPuOMM5rst99++5VPPvnkJmuf/OQny2PGjCmXyy2r7WuvvbZxnxdeeKG8xRZblH/1q1+16O9xxAnYVGxun7/lcrn8q1/9qty1a9fynDlzWvxagI1pc/oMHjlyZLm6urqcpPzZz362vGrVqma/9t044tRMDz74YB555JHsuuuuWblyZeP63nvv3WS/xx9/PKNGjWqyNmrUqDz++OMt/p0jR45s/LlHjx4ZPHjwBr0PwKZsU/38vfvuuzNhwoT8+Mc/zq677triGQDag03xM/hXv/pVHnrooUydOjU33XRTvvvd77Z4hnVxc4h3GDhwYEqlUp588skm6wMGDEiSbLHFFk3W13coc306dVrTquVyuXHN6XcAm9fn7z333JMjjzwyl112WcaPH98mvwOgNW1On8F9+/ZNkgwbNiyrVq3KZz/72Zx55pnp3Lnze3pfR5zeYdttt80hhxySH/zgB3n55Zdb/PqhQ4fmvvvua7J23333ZdiwYUmSnj17JkkWLVrUuP2RRx5Z53s98MADjT8vW7YsTz31VIYOHdrimQA2BZvL5+/06dNzxBFH5JJLLslnP/vZlvwJABWzuXwGv9Pq1avzxhtvZPXq1Rv0+rdzxGkdrrzyyowaNSp77713zj333AwfPjydOnXKzJkz88QTT2TEiBHrfe1ZZ52VT33qU9lzzz0zevTo3HDDDZk2bVruvPPOJGtqfd99983FF1+c/v37Z+nSpTn77LPX+V7nn39+tt122/Tq1Svf+MY3st122zX7O5kee+yxvP7663nxxRezfPnyxv9h7rHHHi35pwDYqDb1z9+77747H/nIRzJx4sR8/OMfz+LFi5MkXbt2dYMIoN3b1D+Dp0yZki5dujTeRn3WrFmZNGlSPv3pT7fO7dJb5UqpzdDzzz9fPu2008r9+/cvd+nSpbzVVluV99lnn/J3vvOd8ssvv1wul9dcGHf99dev9dp3uxVjuVwuP/bYY+WRI0eWt9hii/Iee+xRvv3229d5YdwNN9xQ3nXXXctdu3Yt77PPPuW//OUvzZ6/X79+5SRrPQDau0358/e4445b52fvgQce+F7+SQA2mk35M/jaa68t77XXXuWtttqq3K1bt/KwYcPKF154YfnVV199T/8mbymVy2870ZB2Yfr06Tn44IOzbNmydO/evdLjAHQYPn8BKqe9fwa7xgkAAKCAcNoEHX744dlqq63W+bjwwgsrPR7AZsvnL0DlVPoz2Kl6m6Dnnnsur7766jq39ejRwwXIAG3E5y9A5VT6M1g4AQAAFHCqHgAAQAHhBAAAUEA4AQAAFBBOAAAABYQTAO3e4sWLM3HixAwcODA1NTXp1atXRo0alcmTJ+eVV16p9HgAdABVlR4AAN7N3//+94waNSrdu3fPhRdemN133z3V1dWZPXt2fvSjH2XHHXfMRz/60bVe98Ybb6RLly4VmBiAzZEjTgC0a1/4whdSVVWVWbNm5VOf+lSGDh2aAQMG5KijjspNN92UI488MklSKpUyefLkfPSjH023bt3yrW99K0kyefLkvP/970/Xrl0zePDg/Pd//3fje8+bNy+lUimPPPJI49pLL72UUqmU6dOnJ0mmT5+eUqmUm266KcOHD09NTU323XffzJkzp/E1zz77bI488shss8026datW3bdddfcfPPNbf+PA8BGI5wAaLdeeOGF3H777Tn11FPTrVu3de5TKpUafz733HPzsY99LLNnz84JJ5yQ66+/PhMnTsyZZ56ZOXPm5JRTTsmECRNy9913t3iWs846K5deemlmzpyZnj175sgjj8wbb7yRJDn11FOzcuXK3HvvvZk9e3YuueSSbLXVVhv2RwPQLjlVD4B26+mnn065XM7gwYObrG+33XZ57bXXkqyJlksuuSRJ8u///u+ZMGFC437jxo3L8ccfny984QtJki996Ut54IEH8t3vfjcHH3xwi2Y555xzcsghhyRJfv7zn2ennXbK9ddfn0996lOZP39+Pv7xj2f33XdPkgwYMGDD/mAA2i1HnADY5Dz44IN55JFHsuuuu2blypWN63vvvXeT/R5//PGMGjWqydqoUaPy+OOPt/h3jhw5svHnHj16ZPDgwY3vc/rpp+c///M/M2rUqJxzzjn561//2uL3B6B9E04AtFsDBw5MqVTKk08+2WR9wIABGThwYLbYYosm6+s7nW99OnVa85/BcrncuPbW6XctcdJJJ+Xvf/97/uM//iOzZ8/O3nvvnf/7f/9vi98HgPZLOAHQbm277bY55JBD8oMf/CAvv/xyi18/dOjQ3HfffU3W7rvvvgwbNixJ0rNnzyTJokWLGre//UYRb/fAAw80/rxs2bI89dRTGTp0aONa375987nPfS7Tpk3LmWeemR//+MctnheA9ss1TgC0a1deeWVGjRqVvffeO+eee26GDx+eTp06ZebMmXniiScyYsSI9b72rLPOyqc+9ansueeeGT16dG644YZMmzYtd955Z5Jkiy22yL777puLL744/fv3z9KlS3P22Wev873OP//8bLvttunVq1e+8Y1vZLvttsvYsWOTJGeccUYOP/zw7LLLLlm2bFnuvvvuJlEFwKZPOAHQrr3//e/Pww8/nAsvvDCTJk3KwoULU11dnWHDhuXLX/5y440f1mXs2LG5/PLL893vfjcTJ05M//7987Of/SwHHXRQ4z4//elPc+KJJ2bEiBEZPHhwvv3tb+fDH/7wWu918cUXZ+LEifnb3/6WPfbYIzfccEO6du2aJFm1alVOPfXULFy4MLW1tTnssMNy2WWXtfq/BQCVUyq//cRuAKCJ6dOn5+CDD86yZcvSvXv3So8DQIW4xgkAAKCAcAIAACjgVD0AAIACjjgBAAAUEE4AAAAFhBMAAEAB4QQAAFBAOAEAABQQTgAAAAWEEwAAQAHhBAAAUEA4AQAAFPj/KxPGJWf4DBYAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1000x750 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "errors = np.full((2, ci.shape[0]), np.nan)\n",
    "errors[0, :] = ci['effect'] - ci['2.5 %']\n",
    "errors[1, :] = ci['97.5 %'] - ci['effect']\n",
    "\n",
    "#add true effect\n",
    "ax = plt.subplot(1, 1, 1)\n",
    "ax.scatter(x=['Group_1', 'Group_2', 'Group_3'], y=true_effects, c='red', label='True Effect')\n",
    "\n",
    "plt.errorbar(ci.index, ci.effect, fmt='o', yerr=errors)\n",
    "plt.title('GATEs')\n",
    "plt.xlabel('Groups')\n",
    "_ =  plt.ylabel('Effect and 95%-CI')"
   ]
  }
 ],
 "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.11.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}