07 - Counterfactuals example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "a26d9f5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('fivethirtyeight')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56c9ae74",
   "metadata": {},
   "source": [
    "# Counterfactuals I\n",
    "\n",
    "An example from Peters et al. (2017). \n",
    "\n",
    "Different mechanisms leading to the same observational distribution.\n",
    "\n",
    "Let:\n",
    "\n",
    "* $N_1 \\sim Ber(.5)$ \n",
    "\n",
    "* $N_2 \\sim Ber(.5)$ \n",
    "\n",
    "* $N_3 \\sim U(\\{0, 1, 2\\})$ \n",
    "\n",
    "We have two SCMs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "id": "f0af8f25",
   "metadata": {},
   "outputs": [],
   "source": [
    "class SCM:\n",
    "    \n",
    "    def __init__(self, scm_type):\n",
    "        self.n_1 = stats.bernoulli(.5)\n",
    "        self.n_2 = stats.bernoulli(.5)\n",
    "        self.n_3 = stats.randint(0, 3)\n",
    "        \n",
    "        if scm_type not in (1, 2):\n",
    "            raise AttributeError(f'SCM {scm_type} not defined.')\n",
    "        \n",
    "        self.scm_type = scm_type\n",
    "        \n",
    "        \n",
    "    def get_n_3_term(self, n_3):\n",
    "        \n",
    "        if self.scm_type == 1:\n",
    "            n_3_term = n_3\n",
    "        elif self.scm_type == 2:\n",
    "            n_3_term = 2 - n_3\n",
    "            \n",
    "        return n_3_term\n",
    "\n",
    "        \n",
    "    def sample(self, n=100):\n",
    "        samples = []\n",
    "        x_1 = self.n_1.rvs(n)\n",
    "        x_2 = self.n_2.rvs(n)\n",
    "        n_3 = self.n_3.rvs(n)\n",
    "        \n",
    "        n_3_term = self.get_n_3_term(n_3)\n",
    "   \n",
    "        x_3 = ((n_3 > 0) * x_1 + (n_3 == 0) * x_2) + n_3_term * (x_1 == x_2)\n",
    "        \n",
    "        return x_3\n",
    "    \n",
    "    \n",
    "    def intervene(self, intervention_dict, n=100):\n",
    "        \n",
    "        x_1_intervened = x_2_intervened = False\n",
    "        \n",
    "        for k, v in intervention_dict.items():\n",
    "            \n",
    "            if k == 'x_1':\n",
    "                x_1 = np.ones(n) * v\n",
    "                x_1_intervened = True\n",
    "                \n",
    "            if k == 'x_2':\n",
    "                x_2 = np.ones(n) * v\n",
    "                x_2_intervened = True\n",
    "                \n",
    "            if k not in ['x_1', 'x_2']:\n",
    "                raise Warning(f'Unknown paramenter `{k}`.')\n",
    "                \n",
    "        if not x_1_intervened:\n",
    "            x_1 = self.n_1.rvs(n)\n",
    "        \n",
    "        if not x_2_intervened:\n",
    "            x_2 = self.n_2.rvs(n)\n",
    "                \n",
    "        n_3 = self.n_3.rvs(n)\n",
    "        \n",
    "        n_3_term = self.get_n_3_term(n_3)\n",
    "   \n",
    "        x_3 = ((n_3 > 0) * x_1 + (n_3 == 0) * x_2) + n_3_term * (x_1 == x_2)\n",
    "    \n",
    "        return x_3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "id": "68ea7d3b",
   "metadata": {},
   "outputs": [],
   "source": [
    "scm_1 = SCM(1)\n",
    "scm_2 = SCM(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "id": "df60ab41",
   "metadata": {},
   "outputs": [],
   "source": [
    "scm_1_sample = scm_1.sample(10000)\n",
    "scm_2_sample = scm_2.sample(10000)\n",
    "\n",
    "scm_1_int_x_1_1 = scm_1.intervene({'x_1': 1}, n=10000)\n",
    "scm_2_int_x_1_1 = scm_2.intervene({'x_1': 1}, n=10000)\n",
    "\n",
    "scm_1_int_x_2_1 = scm_1.intervene({'x_2': 1}, n=10000)\n",
    "scm_2_int_x_2_1 = scm_2.intervene({'x_2': 1}, n=10000)\n",
    "\n",
    "scm_1_int_x_1_0 = scm_1.intervene({'x_1': 0}, n=10000)\n",
    "scm_2_int_x_1_0 = scm_2.intervene({'x_1': 0}, n=10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "id": "a5ab21d8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA+EAAAF9CAYAAABrgmn3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAABxyUlEQVR4nO3de1xVdb7/8dcOUBQvECKagqYSCmoohWR5Q9Px+GtIx1vHn1OU4e1M6kyGZKOT9YgUycgQ7/3MciYFp2jyWHrChrxAM0ektB+HRxcvZHDANgqCGu7fHx7Wry0g173ZbN/Px6PHI9b67LW/3+/6br/rs9Z3rWUym80WRERERERERMTm7mjpAoiIiIiIiIjcLpSEi4iIiIiIiNiJknARERERERERO1ESLiIiIiIiImInSsJFRERERERE7ERJuIiIiIiIiIidKAkXERERERERsRMl4f8jLy+vpYvQom73+oPaANQGt3v9QW0AaoPbjfa32gDUBrd7/UFtAGoDe9ZfSbiIiIiIiIiInSgJFxEREREREbETJeEiIiIiIiIidqIkXERERERERMROXFu6ACIiYns///wzZWVldca5u7tTUlJihxI5roa0gYeHB66uGkpFRJxFfcdLZ3S7HwM0tP5NOQbQkYM4pcwNGxv8mfLyci60a2eD0lgbtmCezb9D5Jd+/vlnLl26hKenJyaT6Zaxbdu2xd3d3U4lc0z1bQOLxYLZbKZjx45KxEUcSNzxiw3+zIViN+4sbfjnGiJ2SCebbl+ariHjpTO63Y8BGlL/ph4DaDq6iIiTKysru20PKGzJZDLh6el5214xERFxNhovpb6aegzQ4CQ8ISEBT09Pli5daiyzWCzExcXRv39/unXrxqRJk/j666+tPnflyhWWLl1Knz59uOuuu5g5cyb5+flWMWazmejoaPz9/fH39yc6Ohqz2dyoiomIyP+nAwrbULuKiDgX/bsu9dWUvtKga+dffPEFO3bsIDg42Gp5YmIiSUlJJCUlERAQwJo1a5g8eTJffPEFHTt2BCA2NpZ9+/axbds2vLy8WL58OTNmzOCzzz7DxcUFgDlz5nDu3Dn27NmDyWTimWeeYe7cubz33nuNrqCIiIiI3D4iju5q8GfKy8tpZ+tb0obodjQRuaHeV8JLSkp4+umnWb9+PZ6ensZyi8VCcnIyixcvJjIykqCgIJKTkyktLSUlJcX47M6dO1m1ahVjxowhJCSETZs2cfLkSQ4dOgRAbm4uBw8e5PXXX2fYsGGEhYWxbt06Pv74Y/Ly8pq10iIiIiIiIiItod5JeFWSPWrUKKvlp0+fpqCggIiICGNZu3btGD58OJmZmQBkZ2dz7do1q5iePXsSGBhoxGRlZdGhQweGDRtmxISHh+Ph4WHEiIiIiIiIiLRm9ZqOvmPHDr799ls2bdpUbV1BQQEAPj4+Vst9fHw4f/48AIWFhbi4uODt7V0tprCw0Ijx9va2mltvMpno0qWLEVOT5rxKfrtfcXem+peXl9v1cw3h6O3s6OWzNWesv7u7O23btq13fEVFhQ1L0zo0pA0uXrxY4zgVEBDQnEUSERERJ1FnEp6Xl8eqVav493//d9q0aVNr3M03plssljpvVr85pqb4urbTXAc5eXl5t/UBk7PVvzGvGrPL/WA49oG5s/WDhnLW+peUlFR75UZtr/D5+eefbfK6raa8mueHH37gtdde49NPPyU/P5+OHTsyYMAAVq5cyX333QdATk4O69ev5/DhwxQXF+Pr60t4eDjPPPMMAwcOBOCRRx4hIyOD3//+96xYscLqO6ZOncrBgwd5/PHHWb16dYNe0dKpUyf8/PwaXT8RERFnZTabuf/++/nkk0+4++67W7o49fLb3/6WsLAw/u3f/s1m31HnkVZWVhbFxcU88MADxrLKykqOHDnC9u3bOXbsGHDjSnbPnj2NmKKiIuPqeNeuXamsrKS4uJguXbpYxQwfPtyIKSoqskq6LRYLxcXF1a6yi4jI7eHMmTNEREQQFhbGhg0b6N69O+fPn2ffvn3GyYJ3332XZ555hpkzZ7J9+3Z69OjBmTNn2LFjBxs3buTNN98EbiTqfn5+nDx50uo79u3bR05ODm5uboSEhNi7iiIi4uAa8+75pmjsieumnrTu168fUP+T1omJiXWWKSEhgYcffthhEvDDhw+zfv16Tpw4wfnz50lKSmLWrFlWMTExMUyaNInZs2fTuXNnm5SjznvCJ02axJEjR8jIyDD+GzJkCL/5zW/IyMigX79++Pr6kp6ebnymoqKCo0ePGvd3h4SE4ObmZhWTn59Pbm6uERMWFkZpaSlZWVlGTFZWFmVlZVb3iYuIyO1j06ZNuLi4sHPnTsLDw+nVqxfh4eGsWrWKkJAQsrKyeOaZZ1i1ahVJSUmEh4fj5+fHgw8+yObNm3nppZcA+O677ygpKWHWrFmcOnXK2P6VK1dYvnw5jz/+ONeuXVMSLiIirdKZM2cYOXIkP/zwAxs2bCArK4t33nmHIUOGWJ20HjNmDG3atGH79u384x//IDk5GYCNGzca22quk9aXL1/m7bffZvbs2c1X0SYqKysjKCiIV199tdYZsMHBwfTu3Zvdu3fbrBx1Xgn39PS0eho6QPv27fHy8iIoKAiA+fPnk5CQQEBAAP369WPt2rV4eHgwdepUADp37szs2bNZsWIFPj4+xivKgoODGT16NACBgYGMGzeOJUuWkJiYiMViYcmSJUyYMMEpp4eKiEjdzGYzV69e5cyZMzWeRX/hhRe47777WLhwYY2f9/LyAm48ILRNmzbMnDmTV199lYsXL9KpUyfefPNNunTpQv/+/XFzcyMoKAiLxWLTOomIiDS3X560rnr9c9WJa8DqpPUvx8yqE9c//fQT8P9PWs+fP593333XiPvlSev4+Ph6JeEHDhzgjjvuMMoA8P777xMdHc0//vEP/P39gRtXnj/++GM++eQTunbt2uS2uJXx48czfvx4ABYsWFBr3MSJE0lNTeXpp5+2STnq/XT0W1m0aBELFixg6dKljBkzhh9//JG9e/ca7wgHeOWVV/hf/+t/ERUVxa9+9Ss8PDz4y1/+YnQSgC1btjBw4ECmTJnCb37zGwYOHFjjw+BEROT2MHfuXDp27MjQoUMZOXIkL7zwAjk5OQB88803ZGVlER0dXed2srOzGTBgAL1798bb25uvv/6aH374gddff53Vq1dz4sQJBgwYYDzAbubMmfTq1Yvf/va3Nq2fiIhIc/jlSeuaNOak9dmzZ7l48cZU/JpOWtflyJEjhISEWD3fq+qV1mvXrgVg/fr1pKamkpqa2qAEPCEhgR49etzyvyNHjtR7ezcLDQ3ln//8p80e2tyop+989NFHVn+bTCZiY2OJjY2t9TPu7u7Ex8cTHx9fa4yXlxebN29uTJFERMQJDR48mOPHj5OZmUl6ejppaWls2LCBN99803h42pAhQ+rcTnZ2tnHWfvDgwZw6dYrNmzcTGRnJ0KFDefHFF7n33nuN+AULFvD444/z5z//2Sb1EhERaU5z584lIyODoUOHMmjQIEaOHMn06dMZPHiwcdJ627ZtdW6nppPWfn5+vP7663zwwQe8//771U5aHz16lFGjRvH2229bbevs2bP4+vpaLTOZTKxYsYLp06dz9913k5CQQFpaGn379jVibrXNKk8++SSTJ0++ZV26d+9eZ31r061bN65du8aPP/5ok/vZm+VKuIiIiK24uLgwfPhwli9fzpEjR4z7tC5fvgyAh4dHndvIyckxkuzBgwfz9ttvc+DAAVauXAnAiRMnrKbWjRw5kg4dOjR/ZURERGyg6qT1Rx99xIQJEzhw4ACjR49m165dnDhxAmj8Ses//vGPxknrEydOVDtp/cv7yX+poqKixreNREREMHToUF5++WXeeusthg4darX+Vtus4uXlRZ8+fW75X1PeelT1WVtdCVcSLiIirYbFYuHKlSt06dLFmAp3+PDhGmOrkvTvv/+en376yeqg4vjx4zz33HP4+Pjw/fffYzab9VA2ERFp1RztpLW3tzdms7na8s8++4yvvvoKi8VS41uw6nMi3NbT0avukf/lm72aU/O/DFZERKQZREdHc8899zBq1Ch8fX05c+YMb7zxBhcvXmTx4sUMHDiQ8ePH89xzz1FeXk54eDgmk4kTJ06wY8cOYmJieOCBBzhx4gSurq5G0h4ZGcno0aON+9+q1gcHB7dkdUVERJpNbSetp0yZUi328uXLtG/fntOnT1c7aZ2YmMjLL7/cqJPWgwcPZteuXVbLvvzyS2bPns3q1av55JNPWLVqFXv37m1w/Ww9Hf3UqVN0797dZg+KUxIuIiIO6d577yUtLY3k5GRKS0vp0aMHI0aMICMjg969ewOwc+dONm7cSFJSEkuXLsXNzY3evXszfvx4Y9pddnY2gYGBxpQ4V1dXvL29je+5eb2IiEhr0lwnrXNycpr1pHVERAQrV67kwoUL3HnnnZw5c4bp06ezcOFCZs+eTWhoKA8++CAZGRmMGDGiQXX28vIyytUQpaWlfPvttwBcv36dc+fOkZOTg5eXl9VV+aNHjzJ27NgGb7++lISLiNyGYod0qnF5bfdvtYSFCxfW+hTXKm3btmXRokUsWrSo1piVK1ca0+gas15ERMSRNddJ65ycnGY9aR0cHExoaCipqalMnTqVqVOnMmHCBGJiYgAICgri0UcfZdWqVRw4cKAZW6R2x48f55FHHjH+jouLIy4ujscee4x169YBN46FPvroI1JTU21WDiXhIiIiN4mMjOSrr77i8uXLBAUF8X/+z/8hLCyspYslIiItoLYT146iuU5aL1++nJdeeqnW9Y05aR0TE8OyZct48sknycrKqrb+rbfeatD2mmrEiBE13qcON5JvuHHCIjQ0lPvvv99m5VASLiIicpMPPvigpYsgIiLi8Oo6aT1u3DjmzJlDfn4+/v7+zbJNW3Nzc2PNmjU2/Q4l4SIiIiIiItJg9TlpPW/evGbfpi098cQTNv8OvaJMRERERERExE6UhIuIiIiIiIjYiZJwERERERERETtREi4iIiIiIiJiJ0rCRUREREREROxESbiIiIiIiIiInSgJFxEREREREbETJeEiIiJSp4SEBDw9PVm6dKmxzGKxEBcXR//+/enWrRuTJk3i66+/tvrclStXWLp0KX369OGuu+5i5syZ5OfnW8WYzWaio6Px9/fH39+f6OhozGazPaolIiJid64tXQBHsfm0G3eWXmzpYlQTO6RTSxdBRERuc1988QU7duwgODjYanliYiJJSUkkJSUREBDAmjVrmDx5Ml988QUdO3YEIDY2ln379rFt2za8vLxYvnw5M2bM4LPPPsPFxQWAOXPmcO7cOfbs2YPJZOKZZ55h7ty5vPfee3avq4iIiK0pCf8fk078lXbt2rV0MaobMq+lSyAiIrexkpISnn76adavX8+aNWuM5RaLheTkZBYvXkxkZCQAycnJBAQEkJKSQlRUFCUlJezcuZOkpCTGjBkDwKZNmxg0aBCHDh1i7Nix5ObmcvDgQfbv38+wYcMAWLduHRMnTiQvL4+AgAD7V1pERMSGlISLiNyG2vz1rRqX3/Hzz7i6Nv/QcHVyVKM/+8MPP/Daa6/x6aefkp+fT8eOHRkwYAArV67kvvvuAyAnJ4f169dz+PBhiouL8fX1JTw8nGeeeYaBAwcC8Mgjj5CRkcHvf/97VqxYYfUdU6dO5eDBgzz++OOsXr268RV1QlVJ9qhRo6yS8NOnT1NQUEBERISxrF27dgwfPpzMzEyioqLIzs7m2rVrVjE9e/YkMDCQzMxMxo4dS1ZWFh06dDAScIDw8HA8PDzIzMxUEi4iIk5HSbiIiDisM2fOEBERQVhYGBs2bKB79+6cP3+effv2GScL3n33XZ555hlmzpzJ9u3b6dGjB2fOnGHHjh1s3LiRN998E7iRqPv5+XHy5Emr79i3bx85OTm4ubkREhJi7yo6tB07dvDtt9+yadOmausKCgoA8PHxsVru4+PD+fPnASgsLMTFxQVvb+9qMYWFhUaMt7c3JpPJWG8ymejSpYsRU5u8vLyGV8qG22nNnKkNysvL7fq5+nL0Nnb08tnDxYsXadu2bUsXo0VVVFQ06/bMZjMPPfQQf/vb3+jdu3ezbtsWKioqeOqpp7j//vuZN6/uGckXL16scayq6wSyknAREXFYmzZtwsXFhZ07dxr3D/fq1Yvw8HAAsrKyeOaZZ1i1ahULFy40Pufn58eDDz7ITz/9BMB3331HSUkJ8+fP59133zXirly5wvLly3n88ceJj49XEv4LeXl5rFq1in//93+nTZs2tcb9MnmGG9PUb152s5tjaoqvz3aa4yq5prw7XxtcaMTtheXl5Ta/LdGR29jZ+kBj5OXl0alTJ9zd3autq232mK00dvZYU2eO9evXD3d393rPHEtMTKyzTElJSYwfP57+/fs3qk7N7fDhw6xfv54TJ05w/vx5kpKSmDVrFnAjAXd3d+f5559n0qRJREVF0blz51tur1OnTvj5+TW4HHo6uoiIOCyz2czVq1c5c+ZMjetfeOEF7rvvPqsE/Je8vLwAyM7Opk2bNsycOZOzZ89y8eKNB3G++eabdOnShf79++Pm5kZQUJBtKtIKZWVlUVxczAMPPIC3tzfe3t4cPnyYrVu34u3tzZ133glQ7QpAUVGRcXW8a9euVFZWUlxcfMuYoqIiLBaLsd5isVBcXFztKruIiNTszJkzjBw5kh9++IENGzaQlZXFO++8w5AhQ6xmjo0ZM4Y2bdqwfft2/vGPf5CcnAzAxo0bjW0118yxy5cv8/bbbzN79uzmq2gTlZWVERQUxKuvvlrribfg4GB69+7N7t27bVYOJeEiIuKw5s6dS8eOHRk6dCgjR47khRdeICcnB4BvvvmGrKwsoqOj69xOdnY2AwYMoHfv3nh7e/P111/zww8/8Prrr7N69WpOnDjBgAEDaNu2Lfn5+UyaNIlhw4bx4IMPkpaWZutqOqRJkyZx5MgRMjIyjP+GDBnCb37zGzIyMujXrx++vr6kp6cbn6moqODo0aPG/d0hISG4ublZxeTn55Obm2vEhIWFUVpaSlZWlhGTlZVFWVmZ1X3iIiJSu1/OHAsPDzdmja1atYqQkBCrmWNJSUmEh4cbs8Y2b97MSy+9BPz/mWOzZs3i1KlTxvZ/OXPs2rVr9UrCDxw4wB133GHMXgN4//336dq1q9XJ9ZiYGEJCQuq8Bak5jB8/nhUrVhAZGckdd9SeCk+cOJHU1FSblUNJuIiIOKzBgwdz/PhxPvroIyZMmMCBAwcYPXo0u3bt4sSJEwAMGTKkzu1kZ2cbBwyDBw/m1KlT/PGPfyQyMpKhQ4dy4sQJ7r33XgBcXV2Ji4sjMzOT999/n9jYWC5fvmyzOjoqT09PgoKCrP5r3749Xl5eBAUFYTKZmD9/Pq+//jppaWmcOnWKBQsW4OHhwdSpUwHo3Lkzs2fPZsWKFRw6dIgTJ04wd+5cgoODGT16NACBgYGMGzeOJUuW8MUXX5CVlcWSJUuYMGHCbT89VkSkvhxx5tiRI0cICQmxurUoMjKSoKAg1q5dC8D69etJTU0lNTWVrl271ru+CQkJ9OjR45b/HTlypN7bu1loaCj//Oc/bfasCN0TLiIiDs3FxYXhw4czfPhwli1bxv3338/u3buNRM/Dw6PObeTk5Biv0Ro8eDBvv/0233zzDf/85z8BOHHiBL/+9a8B8PX1pVevXsCNB4h17tyZ4uJi2rdvb4vqtWqLFi2ivLycpUuXYjabCQ0NZe/evcY7wgFeeeUVXFxciIqKoqKigpEjR7Jx40bjHn+ALVu2EBMTw5QpU4AbVyB++SR2ERG5tblz55KRkcHQoUMZNGgQI0eOZPr06QwePNiYObZt27Y6t1PTzDE/Pz9ef/11PvjgA95//31j5ti5c+eYO3cuRUVFuLq6EhMTY4ylAGfPnsXX19dq+yaTiRUrVjB9+nTuvvtuEhISSEtLo2/fvkbMzJkzOXr0KKNGjeLtt9+usZxPPvkkkydPvmVdunfvXmd9a9OtWzeuXbvGjz/+yN13393o7dSmzivhW7ZsYfjw4fj5+eHn58fDDz/Mxx9/bKy3WCzExcXRv39/unXrxqRJk/j666+ttnHlyhWWLl1Knz59uOuuu5g5cyb5+flWMWazmejoaPz9/fH39yc6Ohqz2dw8tRQREadgsVi4cuUKXbp0Mc7CHz58uMbYqqvX33//PT/99JPVlfDjx4/z3HPP4ePjw/fff4/ZbK5xat3x48f5+eef6dmzp03q09p89NFHxMfHG3+bTCZiY2PJzc2loKCAffv2Vbs64u7uTnx8PN999x3nz5/nvffeq9aeXl5ebN68mbNnz3L27Fk2b96Mp6enPaokIuIUHHHmWNWDzm4WERHB0KFDefnll3nrrbcYOnSo1foFCxZY3aNeEy8vL/r06XPL/5rysMWqz9rqSnidSfhdd93Fiy++yGeffUZ6ejojR45k1qxZfPXVVwAkJiaSlJTE6tWr+fTTT/Hx8WHy5MlcunTJ2EZsbCwffvgh27ZtY9++fVy6dIkZM2ZQWVlpxMyZM4ecnBz27NlDSkoKOTk5zJ071wZVFhGR1iA6Opq1a9fyxRdfcObMGT7//HP+9V//lYsXL7J48WKGDh3K+PHjee6553j33Xf55ptv+Pbbb/nrX//Ko48+ahx0nDhxAldXVyM5jIyM5JtvvmHBggVW64ODg62+/8KFC8ybN4/169fX+ZRuERGRllY1c2z58uUcOXLEeLhYVWJc35ljVUl21cyxAwcOsHLlSuDGmFmVpHfr1o3BgwcD1jPHqnh7e9d4UfWzzz7jq6++wmKx1PgAzpEjR9KhQ4dbltPW09Gr3q7SpUuXRm/jVuqcjj5p0iSrv//4xz+ybds2vvjiC4KDg0lOTmbx4sXGNL/k5GQCAgJISUkhKiqKkpISdu7cSVJSEmPGjAFuPDhg0KBBHDp0iLFjx5Kbm8vBgwfZv3+/8RCWdevWMXHiRL0yQUTEBmp7/UltZ61bwr333ktaWhrJycmUlpbSo0cPRowYQUZGhvGu0Z07d7Jx40aSkpJYunQpbm5u9O7dm/Hjxxtn/LOzswkMDDTq5erqavXe6pvXw40ZXLNmzWLJkiV6OJiIiLQ6tc0cq7rt55cuX75M+/btOX36dLWZY4mJibz88suNmjk2ePBgdu3aZRX35ZdfMnv2bFavXs0nn3zCqlWr2Lt3b4PrZ+vp6KdOnaJ79+4Nuk+9IRp0T3hlZSXvv/8+ZWVlhIWFcfr0aQoKCoiIiDBi2rVrx/Dhw8nMzCQqKors7GyuXbtmFdOzZ08CAwPJzMxk7NixZGVl0aFDB6sDnfDwcDw8PMjMzFQSLiJyG1q4cGGtD5Cp0rZtWxYtWsSiRYtqjVm5cqVxBr8+6y0WCwsXLmTkyJHMnDmz4QUXERGxo+joaO655x5GjRqFr68vZ86c4Y033jBmjg0cONCYOVZeXk54eDgmk4kTJ06wY8cOYmJieOCBB8jJyak2c2z06NHGQ9saOnMsIiKClStXcuHCBe68807OnDnD9OnTWbhwIbNnzyY0NJQHH3yQjIwMRowY0aA6e3l5GeVqiNLSUr799lsArl+/zrlz58jJycHLy8vqqvzRo0cZO3Zsg7dfX/VKwk+ePMn48eOpqKjAw8ODd955h+DgYDIzMwGqTSPw8fHh/PnzwI33h7q4uFhddaiKqXoMfWFhId7e3lY7zWQy0aVLlzofVZ+Xl1efKtSLreb8N0Vz1s+RvsvWGrsv7dEHHL2dHb18tuaM9Xd3d6dt27b1jq+oqLBhaRxfVlYWe/fuJSgoiL/97W/AjafCDhgwoMb4ixcv1jhW6QSyiIhzqG32mKNorpljOTk5zTpzLDg4mNDQUFJTU5k6dSpTp05lwoQJxMTEABAUFMSjjz7KqlWrOHDggC2byHD8+HEeeeQR4++4uDji4uJ47LHHWLduHXDjOOijjz6y6SvK6pWEBwQEkJGRQUlJCWlpacyfP984MAGq3StnsVjqvH/u5pia4uuzneY6yMmHJt28byv2Oohztmn/FxqxL8vLy+3SBxy5nZ2tHzSUs9a/pKSk3lPMHWk6eksZNmyYcS9YfXTq1Ak/Pz8blkhERKR2zTVzbPny5cb7wmtS08yxBQsW3HLmWExMDMuWLePJJ58kKyur2vq33nrrluVubiNGjKj14d9VFyF27txJaGgo999/v83KUa/3hLdp04Y+ffowZMgQVq5cyaBBg9iwYYPxyPmbrwAUFRUZV8e7du1KZWWl1U36NcUUFRVhsViM9RaLheLi4hpv1hcREREREZGWc+zYMfbu3ctHH33EQw89xEMPPcTJkyetYsaNG8ecOXOqvRnrViIjI3niiSc4cOAAQUFBNSbvtuTm5mbz12Q26j3h169f5+rVq/Tq1QtfX1/S09ONR8tXVFRw9OhRVq1aBUBISAhubm6kp6czbdo0APLz88nNzTWmLISFhVFaWkpWVpaxLCsri7KyMj0QR0RERERExME88MAD9Zo5Nm/evAZt94MPPmhskZrFE088YfPvqDMJ/9Of/sT48ePp0aMHpaWlpKSk8Pnnn7N7925MJhPz588nISGBgIAA+vXrx9q1a/Hw8GDq1KkAdO7cmdmzZ7NixQp8fHzw8vJi+fLlBAcHM3r0aAACAwMZN24cS5YsITExEYvFwpIlS5gwYYJTTg0VERERERGR21OdSXhBQQHR0dEUFhbSqVMngoODSUlJMZ4Wt2jRIsrLy1m6dClms5nQ0FD27t1Lx44djW288soruLi4EBUVRUVFBSNHjmTjxo24uLgYMVu2bCEmJsZ4bP7EiRNtPg1ARERERERExJ7qTMKTk5Nvud5kMhEbG0tsbGytMe7u7sTHxxMfH19rjJeXF5s3b66rOCIiIiIiIiKtVr0ezCYiIq2Xq6srZWVlVg+/lKazWCyUlZXh6tqox6uIiIjIbUpHDiIiTs7Dw4MrV65w8eLFOmMvXrxIp06d7FAqx9WQNmjoO9hFRERElISLiNwG2rZtW69ksbCw8LZ/57XaQETk9lQ1c6x9+/aYTKaWLo44MIvFwuXLlxs9G05JuIiIiIiI3PYaMnPMGd3us+EaWv+mzIZTEi4iIiIiIkL9Z445o9t9Jpg9668Hs4mIiIiIiIjYiZJwERERERERETtREi4iIiIiIiJiJ0rCRUREREREROxESbiIiIiIiIiInSgJFxEREREREbETJeEiIiIiIiIidqIkXERERERERMROlISLiIiIiIiI2ImScBERERERERE7URIuIiIiIiIiYidKwkVERERERETsREm4iIiIiIiIiJ0oCRcRERERERGxEyXhIiIiIiIiInaiJFxERERERETETpSEi4iIiIiIiNiJknARERERERERO1ESLiIiIiIiImInSsJFRERERERE7ERJuIiIiIiIiIid1JmEv/baa4wZMwY/Pz/69u3LjBkzOHXqlFWMxWIhLi6O/v37061bNyZNmsTXX39tFXPlyhWWLl1Knz59uOuuu5g5cyb5+flWMWazmejoaPz9/fH39yc6Ohqz2dz0WoqIiIiIiIg4gDqT8M8//5ynnnqKjz/+mLS0NFxdXXn00Uf56aefjJjExESSkpJYvXo1n376KT4+PkyePJlLly4ZMbGxsXz44Yds27aNffv2cenSJWbMmEFlZaURM2fOHHJyctizZw8pKSnk5OQwd+7cZq6yiIiIiIiISMtwrStg7969Vn9v2rQJf39/jh07xsSJE7FYLCQnJ7N48WIiIyMBSE5OJiAggJSUFKKioigpKWHnzp0kJSUxZswYYzuDBg3i0KFDjB07ltzcXA4ePMj+/fsZNmwYAOvWrWPixInk5eUREBDQ3HUXERERERERsasG3xNeWlrK9evX8fT0BOD06dMUFBQQERFhxLRr147hw4eTmZkJQHZ2NteuXbOK6dmzJ4GBgUZMVlYWHTp0MBJwgPDwcDw8PIwYERERERERkdaszivhN1u2bBmDBg0iLCwMgIKCAgB8fHys4nx8fDh//jwAhYWFuLi44O3tXS2msLDQiPH29sZkMhnrTSYTXbp0MWJqkpeX19Aq1Kq8vLzZttVcmrN+jvRdttbYfWmPPuDo7ezo5bO1273+oDaA5mkDzeASERGRmjQoCX/++ec5duwY+/fvx8XFxWrdL5NnuPGwtpuX3ezmmJri69pOcx3k5HPjCr6jsddBnLNN+b/QiH1ZXl5ulz7gyO3sbP2goW73+oPaANQGIiIiYlv1no4eGxtLamoqaWlp9O7d21ju6+sLUO1qdVFRkXF1vGvXrlRWVlJcXHzLmKKiIiwWi7HeYrFQXFxc7Sq7iIiIiIiISGtUryQ8JiaGlJQU0tLSuOeee6zW9erVC19fX9LT041lFRUVHD161Li/OyQkBDc3N6uY/Px8cnNzjZiwsDBKS0vJysoyYrKysigrK7O6T1xERERERESktapzOvqzzz7Le++9xzvvvIOnp6dxD7iHhwcdOnTAZDIxf/58EhISCAgIoF+/fqxduxYPDw+mTp0KQOfOnZk9ezYrVqzAx8cHLy8vli9fTnBwMKNHjwYgMDCQcePGsWTJEhITE7FYLCxZsoQJEyZoWqCIiIiIiIg4hTqT8K1btwIYrx+rEhMTQ2xsLACLFi2ivLycpUuXYjabCQ0NZe/evXTs2NGIf+WVV3BxcSEqKoqKigpGjhzJxo0bre4t37JlCzExMUyZMgWAiRMnsmbNmqbXUkRERERERMQB1JmEm83mOjdiMpmIjY01kvKauLu7Ex8fT3x8fK0xXl5ebN68uc7vExEREREREWmNGvyecBERERERERFpHCXhIiIiIiIiInaiJFxERERq9NprrzFmzBj8/Pzo27cvM2bM4NSpU1YxFouFuLg4+vfvT7du3Zg0aRJff/21VcyVK1dYunQpffr04a677mLmzJnk5+dbxZjNZqKjo/H398ff35/o6Oh63RInIiLS2igJFxERkRp9/vnnPPXUU3z88cekpaXh6urKo48+yk8//WTEJCYmkpSUxOrVq/n000/x8fFh8uTJXLp0yYiJjY3lww8/ZNu2bezbt49Lly4xY8YMKisrjZg5c+aQk5PDnj17SElJIScnh7lz59q1viIiIvZQ54PZRERE5Pa0d+9eq783bdqEv78/x44dY+LEiVgsFpKTk1m8eLHxFpXk5GQCAgJISUkhKiqKkpISdu7cSVJSEmPGjDG2M2jQIA4dOsTYsWPJzc3l4MGD7N+/n2HDhgGwbt06Jk6cSF5enl5VKiIiTkVJuIi0CpkbNtr8O8rLy7nQrl2DPjNswTwblUbE8ZSWlnL9+nU8PT0BOH36NAUFBURERBgx7dq1Y/jw4WRmZhIVFUV2djbXrl2ziunZsyeBgYFkZmYyduxYsrKy6NChg5GAA4SHh+Ph4UFmZqaScBERcSpKwkVERKReli1bxqBBgwgLCwOgoKAAAB8fH6s4Hx8fzp8/D0BhYSEuLi54e3tXiyksLDRivL29MZlMxnqTyUSXLl2MmJrk5eU1vVLNuJ3WzJnaoLy83K6fqy9Hb2NHL589qA3UBs1V/7pOHisJFxERkTo9//zzHDt2jP379+Pi4mK17pfJM9x4WNvNy252c0xN8XVtpzmukGu6u/O1QUNnNMGNBLxdIz7XEI7cxs7WBxpDbaA2sGf99WA2ERERuaXY2FhSU1NJS0ujd+/exnJfX1+Aaleri4qKjKvjXbt2pbKykuLi4lvGFBUVYbFYjPUWi4Xi4uJqV9lFRERaOyXhIiIiUquYmBhSUlJIS0vjnnvusVrXq1cvfH19SU9PN5ZVVFRw9OhR4/7ukJAQ3NzcrGLy8/PJzc01YsLCwigtLSUrK8uIycrKoqyszOo+cREREWeg6egiIiJSo2effZb33nuPd955B09PT+MecA8PDzp06IDJZGL+/PkkJCQQEBBAv379WLt2LR4eHkydOhWAzp07M3v2bFasWIGPjw9eXl4sX76c4OBgRo8eDUBgYCDjxo1jyZIlJCYmYrFYWLJkCRMmTLitp0aKiIhzUhIuIiIiNdq6dSuA8fqxKjExMcTGxgKwaNEiysvLWbp0KWazmdDQUPbu3UvHjh2N+FdeeQUXFxeioqKoqKhg5MiRbNy40ere8i1bthATE8OUKVMAmDhxImvWrLF1FUVEROxOSbiIiIjUyGw21xljMpmIjY01kvKauLu7Ex8fT3x8fK0xXl5ebN68uTHFFBERaVV0T7iIiIiIiIiInehKuIiIiIiIOKQ2f33LLt/TrfgCbb66s0GfuTo5ykalcW722qcNNvAhu32VroSLiIiIiIiI2ImuhIuIiMhta/NpN+4svdjSxahR7JBOLV0EkRaXcf6KXb6nvPwO2l1t2HfpBYqNY6992lB3DrTfd+lKuIiIiIiIiIid6Eq4iIiI3LYmnfgr7dq1a+li1GzIvJYugYiI2ICuhIuIiIiIiIjYiZJwERERERERETtREi4iIiIiIiJiJ0rCRUREREREROxESbiIiIiIiIiInSgJFxEREREREbGTeiXhhw8fZubMmQwYMABPT0/effddq/UWi4W4uDj69+9Pt27dmDRpEl9//bVVzJUrV1i6dCl9+vThrrvuYubMmeTn51vFmM1moqOj8ff3x9/fn+joaMxmc9NqKCIiIiIiIuIg6pWEl5WVERQUxKuvvlrjuzQTExNJSkpi9erVfPrpp/j4+DB58mQuXbpkxMTGxvLhhx+ybds29u3bx6VLl5gxYwaVlZVGzJw5c8jJyWHPnj2kpKSQk5PD3Llzm6GaIiIiIiIiIi3PtT5B48ePZ/z48QAsWLDAap3FYiE5OZnFixcTGRkJQHJyMgEBAaSkpBAVFUVJSQk7d+4kKSmJMWPGALBp0yYGDRrEoUOHGDt2LLm5uRw8eJD9+/czbNgwANatW8fEiRPJy8sjICCg2SotIiIiIiIi0hKafE/46dOnKSgoICIiwljWrl07hg8fTmZmJgDZ2dlcu3bNKqZnz54EBgYaMVlZWXTo0MFIwAHCw8Px8PAwYkRERERERERas3pdCb+VgoICAHx8fKyW+/j4cP78eQAKCwtxcXHB29u7WkxhYaER4+3tjclkMtabTCa6dOlixNQkLy+vqVUwlJeXN9u2mktz1s+RvsvWGrsv7dEHHL2dHbV89vp9NvR7HLW9msIZ69RQzdEGmsElIiIiNWlyEl7ll8kz3JimfvOym90cU1N8XdtproOcfKjxfveWZq+DOGeb8n+hEfuyvLzcLn3AkdvZkftBY/ZpQzWmDzhqezWWI/cBe1EbiIiIiC01eTq6r68vQLWr1UVFRcbV8a5du1JZWUlxcfEtY4qKirBYLMZ6i8VCcXFxtavsIiIiIiIiIq1Rk5PwXr164evrS3p6urGsoqKCo0ePGvd3h4SE4ObmZhWTn59Pbm6uERMWFkZpaSlZWVlGTFZWFmVlZVb3iYuIiIiIiIi0VvWajl5aWsq3334LwPXr1zl37hw5OTl4eXnh5+fH/PnzSUhIICAggH79+rF27Vo8PDyYOnUqAJ07d2b27NmsWLECHx8fvLy8WL58OcHBwYwePRqAwMBAxo0bx5IlS0hMTMRisbBkyRImTJigaYEiIiIiIiLiFOqVhB8/fpxHHnnE+DsuLo64uDgee+wxkpOTWbRoEeXl5SxduhSz2UxoaCh79+6lY8eOxmdeeeUVXFxciIqKoqKigpEjR7Jx40ZcXFyMmC1bthATE8OUKVMAmDhxImvWrGmuuoqIiINo89e3WroItRv4UEuXQERERJxYvZLwESNGYDaba11vMpmIjY0lNja21hh3d3fi4+OJj4+vNcbLy4vNmzfXp0giItKKZZy/0tJFqNWdA1u6BCIiIuLMmnxPuIiIiIiIiIjUj5JwERERERERETtREi4iIiIiIiJiJ0rCRUREREREROxESbiIiIiIiIiInSgJFxEREREREbETJeEiIiIiIiIidqIkXERERERERMROlISLiIiIiIiI2ImScBERERERERE7URIuIiIiIiIiYidKwkVERERERETsREm4iIiIiIiIiJ0oCRcRERERERGxEyXhIiIiIiIiInaiJFxERERERETETpSEi4iIiIiIiNiJknARERERERERO1ESLiIiIiIiImInSsJFRERERERE7ERJuIiIiIiIiIidKAkXERERERERsRMl4SIiIiIiIiJ2oiRcRERERERExE6UhIuIiIiIiIjYiZJwERERERERETtxyCR869atDB48GF9fX0aNGsWRI0daukgiIiJiYxr/RUTkduBwSfjevXtZtmwZf/jDH/j73/9OWFgY06ZN4+zZsy1dNBEREbERjf8iInK7MJnNZktLF+KXxo4dS3BwMG+88YaxbOjQoURGRrJy5coWLJmIiIjYisZ/ERG5XTjUlfCrV6+SnZ1NRESE1fKIiAgyMzNbqFQiIiJiSxr/RUTkduJQSXhxcTGVlZX4+PhYLffx8aGwsLCFSiUiIiK2pPFfRERuJw6VhFcxmUxWf1sslmrLRERExLlo/BcRkduBQyXh3t7euLi4VDvrXVRUVO3suIiIiDgHjf8iInI7cagkvE2bNoSEhJCenm61PD09nWHDhrVQqURERMSWNP6LiMjtxLWlC3CzhQsXMnfuXEJDQxk2bBjbt2/nxx9/JCoqqqWLJiIiIjai8V9ERG4XDnUlHGDKlCnExcURHx/PiBEjOHbsGLt378bf37/R29y6dSuDBw/G19eXUaNGceTIkVvGnzx5kn/5l3+hW7duDBgwgNWrV2OxONSb3BqsIW1w+vRpPD09q/138OBBO5a4+Rw+fJiZM2cyYMAAPD09effdd+v8jLP1gYa2gbP1gddee40xY8bg5+dH3759mTFjBqdOnarzc87UDxrTBs7WD7Zs2cLw4cPx8/PDz8+Phx9+mI8//viWn3GmPuDobDH+g44BQMcAOgbQMYCOAW7vYwBHHP8d7ko4wJw5c5gzZ06zbGvv3r0sW7aMhIQEwsPD2bp1K9OmTePYsWP4+flVi7948SKTJ09m+PDhfPrpp+Tl5bFw4ULat2/P7373u2Ypk701tA2qpKamMnDgQONvLy8vexS32ZWVlREUFMRjjz3GvHnz6ox3xj7Q0Dao4ix94PPPP+epp55i6NChWCwWXnnlFR599FEyMzNrrZOz9YPGtEEVZ+kHd911Fy+++CJ9+/bl+vXr/PnPf2bWrFkcOnTIqn5VnK0PtAbNOf6DjgFAxwA6BtAxgI4BdAzgiOO/yWw2t85TOvU0duxYgoODeeONN4xlQ4cOJTIykpUrV1aL37ZtG3/605/4r//6L9q1awdAfHw827dv59SpU63yKa0NbYPTp09z7733kp6ezpAhQ+xZVJvr0aMHa9asYdasWbXGOGMf+KX6tIEz9wGA0tJS/P39effdd5k4cWKNMc7eD+rTBs7eDwB69+7NypUra5zy7Ox94HagYwAdA/ySjgF0DAA6BgAdA0DLj/8ONx29OV29epXs7GwiIiKslkdERJCZmVnjZ7KysnjggQeMBocbA9j58+c5ffq0TctrC41pgyqzZ8+mX79+TJgwgQ8++MCWxXQoztYHmsJZ+0BpaSnXr1/H09Oz1hhn7wf1aYMqztgPKisrSU1NpaysjLCwsBpjnL0PODsdA+gYoDGcrQ80hbP2AR0D3N7HAI4y/jt1El5cXExlZWW115v4+PhUew1KlcLCwhrjq9a1No1pgw4dOvDSSy/x1ltvsWfPHkaOHElUVBTvvfeePYrc4pytDzSGs/eBZcuWMWjQoFr/8QXn7wf1aQNn7AcnT56kR48edO3alSVLlvDOO+8QHBxcY6yz9wFnp2MAHQM0hrP1gcZw9j6gY4Db8xjA0cZ/h7wnvLndPGXAYrHcchpBTfE1LW9NGtIG3t7eVvc7DBkyhAsXLpCYmMiMGTNsWk5H4Yx9oCGcuQ88//zzHDt2jP379+Pi4nLLWGftB/VtA2fsBwEBAWRkZFBSUkJaWhrz58/nb3/7G0FBQTXGO2sfuJ3oGEDHAA3ljH2gIZy5D+gY4PY9BnC08d+pr4R7e3vj4uJS7YxFUVFRtbMbVbp27VpjPFDrZxxZY9qgJqGhoXz77bfNXTyH5Gx9oLk4Qx+IjY0lNTWVtLQ0evfufctYZ+0HDWmDmrT2ftCmTRv69OnDkCFDWLlyJYMGDWLDhg01xjprH7hd6BhAxwCN4Wx9oLk4Qx/QMcDtfQzgaOO/Uyfhbdq0ISQkhPT0dKvl6enpDBs2rMbPhIWFcfToUSoqKqziu3fvTq9evWxaXltoTBvU5Msvv8TX17e5i+eQnK0PNJfW3gdiYmJISUkhLS2Ne+65p854Z+wHDW2DmrT2fnCz69evc/Xq1RrXOWMfuJ3oGEDHAI3hbH2gubT2PqBjAB0D3Kylx3+nTsIBFi5cyK5du3j77bfJzc0lJiaGH3/80XgS3osvvsivf/1rI37q1Km0a9eOBQsWcOrUKdLS0nj99ddZsGBBq51+0tA22LVrF3v27CE3N5e8vDzWr1/P1q1biY6ObqkqNElpaSk5OTnk5ORw/fp1zp07R05ODmfPngVujz7Q0DZwtj7w7LPPsmvXLrZu3YqnpycFBQUUFBRQWlpqxDh7P2hMGzhbP/jTn/7EkSNHOH36NCdPnuTFF1/k888/Z9q0aYDz94HbkY4BdAygYwAdA+gYQMcAjjj+O/094VOmTOHChQvEx8dTUFDAgAED2L17N/7+/gD8+OOPfPfdd0Z8586d+etf/8qzzz7LmDFj8PT0ZOHChfzbv/1bS1WhyRraBgBr167l7NmzuLi40LdvX958881Wef8HwPHjx3nkkUeMv+Pi4oiLi+Oxxx4jOTn5tugDDW0DcK4+sHXrVgAiIyOtlsfExBAbGws4/78FjWkDcK5+UFBQQHR0NIWFhXTq1Ing4GBSUlIYO3Ys4Px94HakYwAdA+gYQMcAOgbQMYAjjv9O/55wEREREREREUfh9NPRRURERERERByFknARERERERERO1ESLiIiIiIiImInSsJFRERERERE7ERJuIiIiIiIiIidKAkXERERERERsRMl4SIiIiIiIiJ2oiRcRERERERExE6UhIuIiIiIiIjYiZJwERERERERETtREi4iIiIiIiJiJ0rCRUREREREROxESbiIiIiIiIiInSgJFxEREREREbETJeEiIiIiIiIidqIkXERERERERMROlISLiIiIiIiI2ImS8P+Rl5fX0kVoUbd7/UFtAGqD273+oDYAtcHtxhn2d2uvQ2svP7T+OrT28kPrr0NrLz+0/jrYs/xKwkVERERERETsREm4iIiIiIiIiJ0oCRcRERERERGxEyXhIiIiIiIiInbi2tIFEBFxRD///DNlZWUtXQy7c3d3p6SkpKWL0aIa0gYeHh64umooFRFxFi01/jvD+Nva69DQ8jflGEBHDv+j22dptPnqzpYuRjVXJ0e1dBFEbjs///wzly5dwtPTE5PJ1NLFsau2bdvi7u7e0sVoUfVtA4vFgtlspmPHjkrEW7HNp924s/RiSxejRrFDOrV0EURuKy05/jvD+Nva69CQ8jf1GEBHDf/jnyV30O7qlZYuRjXDWroArVTmho0N/kx5eTkX2rWzQWmsDVswz+bfIU1TVlZ2Wybg0jAmkwlPT08uXrxI586dW7o4IiLSRBr/pb6aegygJFxEpAYagKU+1E9av0kn/ko7O5yAbZQhOmkrYm/6d13qqyl9RQ9mExEREREREbETJeEiIiIiIiIidqIkXERERERERMROlISLiIiIiIiI2ImScBERsTmz2UxAQADfffddSxel3n7729/y5ptvtnQxREREWi2N/zXT09FFROoh7rh93yPclPcD//DDD7z22mt8+umn5Ofn07FjRwYMGMDKlSu57777AMjJyWH9+vUcPnyY4uJifH19CQ8PZ968eQwdOhSARx55hIyMDH7/+9+zYsUKq++YOnUqBw8e5PHHHycxMbHOMiUkJPDwww9z9913N7pezenw4cOsX7+eEydOcP78eZKSkpg1a5ZVTExMDJMmTWL27Nl6BZmIyG3MHscAP//8M66uV4HGHwM0Zfx/5plnGDhwIODc4z/A1q1beeONNygoKKB///7ExcUxfPhwY709xn9dCRcRcSJnzpxh5MiR/PDDD2zYsIGsrCzeeecdhgwZgqvrjfOu7777LmPGjKFNmzZs376df/zjHyQnJwOwZcsWY1s5OTn4+flx8uRJq+/Yt28fOTk5uLm5ERISUmeZLl++zNtvv83s2bObr6JNVFZWRlBQEK+++mqtr6cKDg6md+/e7N69286lExERaZimjv8bN240tuXM4//evXtZtmwZf/jDH/j73/9OWFgY06ZN4+zZs0aMPcZ/JeEiIk5k06ZNuLi4sHPnTsLDw+nVqxfh4eGsWrWKkJAQsrKyeOaZZ1i1ahVJSUmEh4fj5+fHgw8+yObNm40z3t999x0lJSXMmjWLU6dOGdu/cuUKy5cv5/HHH+fatWv1GoQPHDjAHXfcQXh4uLHs/fffp2vXrpw5c8ZYFhMTQ0hICIWFhc3XILUYP348K1asIDIykjvuqH0onDhxIqmpqTYvj4iISFM0dfx/6aWXAOcf/5OSkvjXf/1XHn/8cQIDA4mPj8fX15ft27dbxdl6/FcSLiLiRMxmM1evXrUa3H7phRde4L777mPhwoU1rvfy8gIgOzubNm3aMHPmTM6ePcvFizem4r355pt06dKF/v374+bmRlBQUJ1lOnLkCCEhIZhMJmNZZGQkQUFBrF27FoD169eTmppKamoqXbt2rXd9ExIS6NGjxy3/O3LkSL23d7PQ0FD++c9/Ul5e3uhtiIiI2JrG/7rH/6tXr5KdnU1ERITV8oiICDIzM62W2Xr8r/Oe8C1btvDWW28Zl+j79+/Ps88+y4QJEwCwWCy8+uqr7NixA7PZTGhoKGvXrmXAgAHGNq5cucILL7xAamoqFRUVjBw50mi4Kmazmeeee479+/cD8Ktf/Yo1a9bg6enZnPUVEXFqc+fOJSMjg6FDhzJo0CBGjhzJ9OnTGTx4MN988w1ZWVls27atzu1kZ2czYMAAevfujbe3N19//TV+fn68/vrrfPDBB7z//vsMGDCAtm3bAjBz5kyOHj3KqFGjePvtt622dfbsWXx9fa2WmUwmVqxYwfTp07n77rtJSEggLS2Nvn37AnDu3Dnmzp1LUVERrq6uxMTE8Otf/7paOZ988kkmT558y7p07969zvrWplu3bly7do0ff/zRoe5nExER+SVHGP83b95stS1HG/+Li4uprKzEx8fHarmPj0+1q/C2Hv/rTMLvuusuXnzxRfr27cv169f585//zKxZszh06BADBw4kMTGRpKQkkpKSCAgIYM2aNUyePJkvvviCjh07AhAbG8u+ffvYtm0bXl5eLF++nBkzZvDZZ5/h4uICwJw5czh37hx79uzBZDLxzDPPMHfuXN57771mr7SIiLMaPHgwx48fJzMzk/T0dNLS0tiwYQNvvvkm7u7uAAwZMqTO7WRnZxtTzQYPHsypU6fYvHkzkZGRDB06lBdffJF7773XiF+wYAGPP/44f/7zn6ttq6Kiosaz2xEREQwdOpSXX36Zv/zlL8YD4QBcXV2Ji4tj8ODB/Pd//zejR49m3LhxtG/f3mobXl5extl7W6i6X1xXwkVExJFp/K+/X16ZhxsXlW9eZuvxv87p6JMmTeLhhx+mT58+9OvXjz/+8Y906NCBL774AovFQnJyMosXLzamFiQnJ1NaWkpKSgoAJSUl7Ny5k1WrVjFmzBhCQkLYtGkTJ0+e5NChQwDk5uZy8OBBXn/9dYYNG0ZYWBjr1q3j448/Ji8vzyYVFxFxVi4uLgwfPpzly5dz5MgR4+Eily9fBsDDw6PObeTk5BiD7ODBg3n77bc5cOAAK1euBODEiRNW94ONHDmSDh061Lgtb29vzGZzteWfffYZX331FRaLpdpZ6W7dujF48GDgxhnqzp07U1xcXG0btp6O/tNPPwHQpUuXRm9DRETEHjT+33r89/b2xsXFpdpV76KiomrlsPX436B7wisrK0lNTaWsrIywsDBOnz5NQUGB1bz6du3aMXz4cGNefXZ2NteuXbOK6dmzJ4GBgUZMVlYWHTp0YNiwYUZMeHg4Hh4e1ebni4hI/VksFq5cuUKXLl2M+7cOHz5cY2zVIP3999/z008/WZ0JP378OM899xw+Pj58//33mM3mej2Uperzubm5Vsu+/PJLZs+ezerVq5k0aRKrVq2q9fPHjx/n559/pmfPntXWPfnkk2RkZNzyv/qc+a/NqVOn6N69e4PuUxMREWlpGv+rj/9t2rQhJCSE9PR0q+Xp6elWeSjYfvyv13vCT548yfjx46moqMDDw4N33nmH4OBgI0GuaV79+fPnASgsLMTFxQVvb+9qMVVnIQoLC/H29raaBmAymejSpUudT8lrzivljjjd0J4zAZxp1kFj96U9+oCjt7Ojl8/W8vLycHd3N+51qvLzzz/btRwVFRUN/szChQsJCAhgxIgRdO3albNnz5KUlERJSQkLFiwgKCiIsWPHsnTpUi5evEhYWBgmk4mcnBzeffdd/vCHPzBs2DC++OILXF1d6dOnDxUVFUyYMIGTJ0/i5eVFRUWFsb5v375W5bx69SqVlZXVyv7QQw+xcuVKfvjhB+68807Onj3LtGnTmDt3LtOmTWPgwIFERETwH//xHzz44INWn71w4QLR0dG89tprXLlypVqd27VrV+srxm7VnmVlZXz33XcAXL9+ne+//55//OMfeHp60rNnTyP+888/Z/To0bfcHxcvXqxxrAoICKizXCIiIk0VHR3NPffcw6hRo/D19eXMmTO88cYbXLx4kcWLFzNw4EDGjx/Pc889R3l5OeHh4ZhMJk6cOMGOHTuIiYnhgQce4MSJE7i6uhpJe2RkJKNHjzamfVetDw4Orle5IiIiWLlyJRcuXODOO+/kzJkzTJ8+nYULFzJ79mxCQ0N58MEHycjIYMSIEVafvXDhAvPmzWP9+vXVpopD46ejL1y4kLlz5xIaGsqwYcPYvn07P/74I1FRUVZxR48eZezYsQ3efn3VKwkPCAggIyODkpIS0tLSmD9/Pn/729+M9fWZV3+zm2Nqiq/PdprrICcf6nUgZ2/2OojLy8tzqgPGC43Yl+Xl5XbpA47czs7WDxqqqv4lJSXG/VNVXF2v2rUsN39/fQwdOpS0tDS2bNlCaWkpPXr0YMSIESQkJNC7d28Adu3axcaNG9m8eTPPP/88bm5u9O7dm/HjxxtngU+dOkVgYKDVgzF/OdWspvVw4wyzi4tLtbIPGTKE0NBQPvroI6ZOncqsWbP41a9+xfLly431jz76KK+++ioHDhwwPnflyhWeeuop/vCHP1QbnJvqiy++4JFHHjH+jo+PJz4+nscee4x169bh7u5ORUUF+/fvJzU19Zb7o1OnTvj5+TVr+UREROrr3nvvJS0tzbgtuGr8z8jIMMb/nTt3snHjRpKSkli6dKnV+F91xTg7O5vAwEBjzHN1dbW6kHrz+roEBwcTGhpKamoqU6dOZerUqUyYMIGYmBgAgoKCePTRR1m1alW18X/WrFksWbKk2hXqppoyZQoXLlwgPj6egoICBgwYwO7du/H39zdOuFdUVPDRRx/Z9BVl9UrC27RpQ58+fYAbB0v/+Z//yYYNG3j22WeBG1eyfzlN4Jfz6rt27UplZSXFxcVWc+qLiooYPny4EVNUVGSVdFssFoqLi6tdZRcRaQmxQzq1dBHqtHDhwlpfPVKlbdu2LFq0iEWLFtW4vqKigpUrVxr3ftWkrvU1iYmJYdmyZTz55JNkZWVVW//WW29Z/W2xWFiwYAEjR45k5syZDfqu+hgxYkSN96nB/79qvnPnTkJDQ7n//vub/ftFRKT1sMcxQEVFRaNOwEPzjP9Q9/juDOM/3Hgg+Jw5c2pdb4/xv1HvCb9+/TpXr16lV69e+Pr6Ws2rr6io4OjRo8ZZi5CQENzc3Kxi8vPzyc3NNWLCwsIoLS212jFZWVmUlZU1+9kPERFpfpGRkTzxxBMcOHCAoKCgagPtuHHjmDNnDvn5+fXa3rFjx9i7dy8fffQRDz30EA899BAnT560RdFr5ebmxpo1a+z6nSIiIq3JL8f/IUOGaPyvpzqvhP/pT39i/Pjx9OjRw3jq+eeff87u3bsxmUzMnz+fhIQEAgIC6NevH2vXrsXDw4OpU6cC0LlzZ2bPns2KFSvw8fExXlEWHBzM6NGjAQgMDGTcuHEsWbKExMRELBYLS5YsYcKECbf11FgRkdbigw8+qDNm3rx59d7eAw88YDyZtKU88cQTLfr9IiIiju6X439tV/M1/ldXZxJeUFBAdHQ0hYWFdOrUieDgYFJSUowb1RctWkR5eTlLly7FbDYTGhrK3r17jXeEA7zyyiu4uLgQFRVFRUUFI0eOZOPGjcY7wgG2bNlCTEwMU6ZMAWDixIm6AiEiIiIiIiJOpc4kPDk5+ZbrTSYTsbGxxMbG1hrj7u5uPPSmNl5eXmzevLmu4oiIiIiIiIi0Wo26J1xEREREREREGk5JuIiIiIiIiIidKAkXERERERERsRMl4SIiIiIiIiJ2oiRcRERERERExE6UhIuIiIiIiIjYiZJwERERqVNCQgKenp4sXbrUWGaxWIiLi6N///5069aNSZMm8fXXX1t97sqVKyxdupQ+ffpw1113MXPmTPLz861izGYz0dHR+Pv74+/vT3R0NGaz2R7VEhERsTsl4SIiInJLX3zxBTt27CA4ONhqeWJiIklJSaxevZpPP/0UHx8fJk+ezKVLl4yY2NhYPvzwQ7Zt28a+ffu4dOkSM2bMoLKy0oiZM2cOOTk57Nmzh5SUFHJycpg7d67d6iciImJPSsJFRMTmzGYzAQEBfPfddy1dlHr77W9/y5tvvtnSxWhxJSUlPP3006xfvx5PT09jucViITk5mcWLFxMZGUlQUBDJycmUlpaSkpJifHbnzp2sWrWKMWPGEBISwqZNmzh58iSHDh0CIDc3l4MHD/L6668zbNgwwsLCWLduHR9//DF5eXktUGMREWkuGv9r5mrTrYuIOIk2f33Lrt93dXJUoz/7ww8/8Nprr/Hpp5+Sn59Px44dGTBgACtXruS+++4DICcnh/Xr13P48GGKi4vx9fUlPDycefPmMXToUAAeeeQRMjIy+P3vf8+KFSusvmPq1KkcPHiQxx9/nMTExDrLlJCQwMMPP8zdd9/d6Ho1p8OHD7N+/XpOnDjB+fPnSUpKYtasWVYxMTExTJo0idmzZ9O5c+cWKmnLq0qyR40axZo1a4zlp0+fpqCggIiICGNZu3btGD58OJmZmURFRZGdnc21a9esYnr27ElgYCCZmZmMHTuWrKwsOnTowLBhw4yY8PBwPDw8yMzMJCAgoNayNVeSXl5e3izbaW4NqV9rP2HR2ssPrb8Orb380PQ6uLu707Zt22rL23+4s0nbrY/2v/j/y4/MbtQ2zp8/T2JiIp999hk//PADHTp0oH///ixfvtwY27/66is2bNjA0aNHuXDhAl27diUsLIyFCxcSFBQEwG9+8xsOHz7MM888w/PPP2/1HY899hjp6en87//9v1m7dm21MlRUVFj9vXr1aiIiIujevXu1dS3lrbfeYsOGDRQWFhIYGMiqVasIDw8HbpR/yZIlTJ48menTp9OpU6dbbuvixYsUFhZWW36rsQuUhIuIOJUzZ84QERFBWFgYGzZsoHv37pw/f559+/bh6nrjn/x3332XZ555hpkzZ7J9+3Z69OjBmTNn2LFjB1u2bCE5ORm4kaj7+flx8uRJq+/Yt28fOTk5uLm5ERISUmeZLl++zNtvv81f/vKXZq9vY5WVlREUFMRjjz3GvHnzaowJDg6md+/e7N69m6efftrOJXQMO3bs4Ntvv2XTpk3V1hUUFADg4+NjtdzHx4fz588DUFhYiIuLC97e3tViqg5aCgsL8fb2xmQyGetNJhNdunSp8cDml+o6yKmPfG6cPHBE9a1fXl5es7RFS2nt5YfWX4fWXn5onjqUlJTg7u5ebXnV+GlLP//8s/E9NZWhLmfOnOHhhx8mLCyM5ORkq/G/ffv2uLu7W43/b731ltX4v337duPq71dffYWfnx+5ublWZdm3bx8nT57Ezc2N0NDQauWsqKiwWnb58mV27drFX/7yl0bVyRb27t3LH//4RxISEggPD2fr1q3MmjWLY8eO4ePjg7u7O0OGDKF379588MEHdY7/nTp1ws/Pr8Hl0HR0EREnsmnTJlxcXNi5cyfh4eH06tWL8PBwVq1aRUhICFlZWTzzzDOsWrWKpKQkwsPD8fPz48EHH2Tz5s3GFe/vvvuOkpISZs2axalTp4ztX7lyheXLl/P4449z7dq1eiXhBw4c4I477jDOMgO8//77dO3alTNnzhjLYmJiCAkJqTPxag7jx49nxYoVREZGcscdtQ+FEydOJDU11eblcUR5eXmsWrWKLVu20KZNm1rjfpk8w41p6jcvu9nNMTXF12c7IiJyQ1PH/5deeglw/vE/KSmJf/3Xf+Xxxx8nMDCQ+Ph4fH192b59u1Wcrcd/JeEiIk7EbDZz9epVq8Htl1544QXuu+8+Fi5cWON6Ly8vALKzs2nTpg0zZ87k7NmzXLx4EYA333yTLl260L9/f9zc3Iypa7dy5MgRQkJCrBKqqnuIq6ayrV+/ntTUVFJTU+natWu965uQkECPHj1u+d+RI0fqvb2bhYaG8s9//tNhpyvbUlZWFsXFxTzwwAN4e3vj7e3N4cOH2bp1K97e3tx5550A1Q6aioqKjKvjXbt2pbKykuLi4lvGFBUVYbFYjPUWi4Xi4uJqV9lFRKRmGv/rHv+vXr1Kdna21S1SABEREWRmZlots/X4r+noIiJOZO7cuWRkZDB06FAGDRrEyJEjmT59OoMHD+abb74hKyuLbdu21bmd7OxsBgwYQO/evfH29ubrr7/Gz8+P119/nQ8++ID333+fAQMG0LZtW86dO8fcuXMpKirC1dWVmJgYfv3rXxvbOnv2LL6+vlbbN5lMrFixgunTp3P33XeTkJBAWloaffv2Bahzm1WefPJJJk+efMu6dO/evT5NV6Nu3bpx7do1fvzxR4e5n91eJk2axJAhQ6yWLVy4kL59+/L73/+efv364evrS3p6unGvYUVFBUePHmXVqlUAhISE4ObmRnp6OtOmTQMgPz+f3Nxc4x7wsLAwSktLycrKMpZlZWVRVlZmdZ+4iIjUzhHG/9///vf85je/MbblaON/cXExlZWVNd5GdfMJZVuP/0rCRUScyODBgzl+/DiZmZmkp6eTlpbGhg0bePPNN437sW5OrGqSnZ1tTDUbPHgwp06dYvPmzURGRjJ06FBefPFF7r33XuDGvXJxcXEMHjyY//7v/2b06NGMGzeO9u1vPGamoqKixrPbERERDB06lJdffpm//OUvRiJXn21W8fLyMs7e20LVvcK345VwT09Pq6ehA7Rv3x4vLy/jCsj8+fNJSEggICCAfv36sXbtWjw8PJg6dSoAnTt3Zvbs2axYsQIfHx+8vLxYvnw5wcHBjB49GoDAwEDGjRvHkiVLSExMxGKxsGTJEiZMmNDq71EVEbEXRxj/R40axcSJEx1+/K/PbVS2Hv+VhIuIOBkXFxeGDx/O8OHDWbZsGffffz+7d+82EiMPD486t5GTk0NkZCRwYxB+++23+eabb/jnP/8JwIkTJ4wz0926daNbt27AjbPJnTt3pri42Bgwvb29MZvN1b7js88+46uvvsJisVQ7K13XNqskJCTw2muv3bIue/bsYfjw4XXWuSY//fQTAF26dGnU553dokWLKC8vZ+nSpZjNZkJDQ9m7dy8dO3Y0Yl555RVcXFyIioqioqKCkSNHsnHjRlxcXIyYLVu2EBMTw5QpU4Ab9+L98knsIiJSN43/1m4e/729vXFxcbnlbVRVbD3+KwkXEXFiFouFK1eu0KVLF+Pq5eHDh41k55cuX77MHXfcwffff89PP/1kdSY8MTGRl19+GR8fH77//nvMZnOND2U5fvw4P//8Mz179jSWDR48mF27dlnFffnll8yePZvVq1fzySefsGrVKvbu3VtjHWraZhVbT0c/deoU3bt3b9B9as7so48+svrbZDIRGxtLbGxsrZ9xd3cnPj6e+Pj4WmO8vLzYvHlzs5VTROR219Dxv3379k4//rdp04aQkBDS09N59NFHjeXp6enVprzbevxXEi4i4iSio6O55557GDVqFL6+vpw5c4Y33niDixcvsnjxYgYOHMj48eN57rnnKC8vJzw8HJPJxIkTJ9ixYwcxMTEMGTKEEydO4OrqagzakZGRjB492pj2VbU+ODjY6vsvXLjAvHnzWL9+vdW0roiICFauXMmFCxe48847OXPmDNOnT2fhwoXMnj2b0NBQHnzwQTIyMhgxYkS9tlmlsdPRSktL+fbbbwG4fv06586dIycnBy8vL6uz4UePHmXs2LEN3r6IiIi9NMf4/8ADDzR5/H/ttdccfvxfuHAhc+fOJTQ0lGHDhrF9+3Z+/PFHoqKirOJsPf4rCRcRcRL33nsvaWlpJCcnU1paSo8ePRgxYgQZGRn07t0bgJ07d7Jx40aSkpJYunQpbm5u9O7dm/Hjxxv3imVnZxMYGGjcQ+bq6mr1nueb18ONV5fMmjWLJUuWVHuYVnBwMKGhoaSmpjJ16lSmTp3KhAkTiImJASAoKIhHH32UVatWceDAgXpts6mOHz/OI488YvwdFxdHXFwcjz32GOvWrQNu3Mv20Ucf3bavKBMRkdbBUcb/+++/36pcjjj+T5kyhQsXLhAfH09BQQEDBgxg9+7d+Pv7U1FRAdhn/DeZzWZL3WHO71DCOuMGfEcybME8u3xPXl6eUz0AJ3PDxgZ/pry83C59wF77tDGcrR80VFX9S0pK6Ny5c0sXp0VUVFRYDa71YbFYmDNnDv369at1WvLBgwdZtmwZmZmZVvcCN2WbtlLVBlu2bGHfvn389a9/vWX87dxfnIGjjv9Q//Gitf/b3drLD62/Dq29/NA8dWjJf88bM/62tJvH6prq4MzjPzS+z+g94SIi0iTHjh1j7969fPTRRzz00EM89NBDnDx50ipm3LhxzJkzh/z8/Gbbpq25ubnp4WAiIiK1uHmsHjt2rMb/etJ0dBERaZIHHnjAeIrorcybV/9ZIPXdpi098cQTLfr9IiIijuzmsbq2q/ka/6vTlXARERERERERO1ESLiIiIiIiImInSsJFRERERERE7ERJuIiIiIiIiIidKAkXERERERERsRMl4SIiNbBYLC1dBGkF1E9ERJyL/l2X+mpKX9ErykREbuLh4YHZbMbT0xOTydTSxREHZbFYMJvNdOzYsaWLIiLSJN0+S6PNV3e2dDFqdHVylN2+S+O/1FdTjwGUhIuI3MTV1ZWOHTty8eLFli6K3V28eJFOnTq1dDFaVEPaoGPHjri6aigVEXEGLTn+O8P429rr0NDyN+UYQEcOIiI1cHV1pXPnzi1dDLsrLCzEz8+vpYvRotQGIiK3r5Ya/51h7GntdbBn+XVPuIiIiIiIiIid6Eq4iIiIiMht7J8ld9Du6pWWLkaNhrV0AURsQFfCRUREREREROykziT8tddeY8yYMfj5+dG3b19mzJjBqVOnrGIsFgtxcXH079+fbt26MWnSJL7++murmCtXrrB06VL69OnDXXfdxcyZM8nPz7eKMZvNREdH4+/vj7+/P9HR0ZjN5qbXUkRERERERMQB1JmEf/755zz11FN8/PHHpKWl4erqyqOPPspPP/1kxCQmJpKUlMTq1av59NNP8fHxYfLkyVy6dMmIiY2N5cMPP2Tbtm3s27ePS5cuMWPGDCorK42YOXPmkJOTw549e0hJSSEnJ4e5c+c2c5VFREREREREWkad94Tv3bvX6u9Nmzbh7+/PsWPHmDhxIhaLheTkZBYvXkxkZCQAycnJBAQEkJKSQlRUFCUlJezcuZOkpCTGjBljbGfQoEEcOnSIsWPHkpuby8GDB9m/fz/Dht24+2PdunVMnDiRvLw8AgICmrvuIiIiIiIiInbV4HvCS0tLuX79Op6engCcPn2agoICIiIijJh27doxfPhwMjMzAcjOzubatWtWMT179iQwMNCIycrKokOHDkYCDhAeHo6Hh4cRIyIiIiIiItKaNfjp6MuWLWPQoEGEhYUBUFBQAICPj49VnI+PD+fPnwduvHPNxcUFb2/vajGFhYVGjLe3NyaTyVhvMpno0qWLEVOTvLy8hlahVuXl5c22rebSnPVzpO+ytcbuS3v0AUdvZ0cvn63d7vUHtQE0TxtoBpeIiIjUpEFJ+PPPP8+xY8fYv38/Li4uVut+mTzDjYe13bzsZjfH1BRf13aa6yAnnxtX8B2NvQ7inG3K/4VG7Mvy8nK79AFHbmdn6wcNdbvXH9QGoDYQERER26r3dPTY2FhSU1NJS0ujd+/exnJfX1+Aaleri4qKjKvjXbt2pbKykuLi4lvGFBUVYbFYjPUWi4Xi4uJqV9lFREREREREWqN6JeExMTGkpKSQlpbGPffcY7WuV69e+Pr6kp6ebiyrqKjg6NGjxv3dISEhuLm5WcXk5+eTm5trxISFhVFaWkpWVpYRk5WVRVlZmdV94iIiIiIiIiKtVZ3T0Z999lnee+893nnnHTw9PY17wD08POjQoQMmk4n58+eTkJBAQEAA/fr1Y+3atXh4eDB16lQAOnfuzOzZs1mxYgU+Pj54eXmxfPlygoODGT16NACBgYGMGzeOJUuWkJiYiMViYcmSJUyYMEHTAkVERERERMQp1JmEb926FcB4/ViVmJgYYmNjAVi0aBHl5eUsXboUs9lMaGgoe/fupWPHjkb8K6+8gouLC1FRUVRUVDBy5Eg2btxodW/5li1biImJYcqUKQBMnDiRNWvWNL2WIiIiIiIiIg6gziTcbDbXuRGTyURsbKyRlNfE3d2d+Ph44uPja43x8vJi8+bNdX6fiIiIiIiISGvU4PeEi4iIiIiIiEjjKAkXERERERERsRMl4SIiIiIiIiJ2oiRcRERERERExE6UhIuIiIiIiIjYiZJwERERERERETup8xVlIiKOoM1f37L5d3QrvkCbr+5s0GeuTo6yUWlERERExBkpCReRViHj/BWbf0d5+R20u9qw7xlmo7KIiIiIiHPSdHQRERERERERO9GVcBERsbu44xdbugi1mtqhpUsgIiIizkxXwkVERERERETsRFfCRUTE7iKO7mrpItTu4bEtXQIRERFxYroSLiIiIiIiImInSsJFRERERERE7ERJuIiIiIiIiIidKAkXERERERERsRMl4SIiIiIiIiJ2oiRcRERERERExE6UhIuIiIiIiIjYiZJwERERERERETtREi4iIiI1eu211xgzZgx+fn707duXGTNmcOrUKasYi8VCXFwc/fv3p1u3bkyaNImvv/7aKubKlSssXbqUPn36cNdddzFz5kzy8/OtYsxmM9HR0fj7++Pv7090dDRms9nWVRQREbE7JeEiIiJSo88//5ynnnqKjz/+mLS0NFxdXXn00Uf56aefjJjExESSkpJYvXo1n376KT4+PkyePJlLly4ZMbGxsXz44Yds27aNffv2cenSJWbMmEFlZaURM2fOHHJyctizZw8pKSnk5OQwd+5cu9ZXRETEHlxbugAiIiLimPbu3Wv196ZNm/D39+fYsWNMnDgRi8VCcnIyixcvJjIyEoDk5GQCAgJISUkhKiqKkpISdu7cSVJSEmPGjDG2M2jQIA4dOsTYsWPJzc3l4MGD7N+/n2HDhgGwbt06Jk6cSF5eHgEBAfatuIiIiA3pSriIiIjUS2lpKdevX8fT0xOA06dPU1BQQEREhBHTrl07hg8fTmZmJgDZ2dlcu3bNKqZnz54EBgYaMVlZWXTo0MFIwAHCw8Px8PAwYkRERJyFroSLiIhIvSxbtoxBgwYRFhYGQEFBAQA+Pj5WcT4+Ppw/fx6AwsJCXFxc8Pb2rhZTWFhoxHh7e2MymYz1JpOJLl26GDE1ycvLa3qlgPLy8mbZTnNrSP2aqy1aSmsvP7T+Ouh30PJae/mh9dehucpf1wwuJeEiIiJSp+eff55jx46xf/9+XFxcrNb9MnmGGw9ru3nZzW6OqSm+ru00xzT1fG5cvXdE9a1fa5+y39rLD62/DvodtLzWXn5o/XWwZ/k1HV1ERERuKTY2ltTUVNLS0ujdu7ex3NfXF6Da1eqioiLj6njXrl2prKykuLj4ljFFRUVYLBZjvcViobi4uNpVdhERkdZOSbiIiIjUKiYmhpSUFNLS0rjnnnus1vXq1QtfX1/S09ONZRUVFRw9etS4vzskJAQ3NzermPz8fHJzc42YsLAwSktLycrKMmKysrIoKyuzuk9cRETEGWg6uoiIiNTo2Wef5b333uOdd97B09PTuAfcw8ODDh06YDKZmD9/PgkJCQQEBNCvXz/Wrl2Lh4cHU6dOBaBz587Mnj2bFStW4OPjg5eXF8uXLyc4OJjRo0cDEBgYyLhx41iyZAmJiYlYLBaWLFnChAkTWvXURhERkZooCRcREZEabd26FcB4/ViVmJgYYmNjAVi0aBHl5eUsXboUs9lMaGgoe/fupWPHjkb8K6+8gouLC1FRUVRUVDBy5Eg2btxodW/5li1biImJYcqUKQBMnDiRNWvW2LqKIiIidqckXERERGpkNpvrjDGZTMTGxhpJeU3c3d2Jj48nPj6+1hgvLy82b97cmGKKiIi0KronXERERERERMROlISLiIiIiIiI2Immo4uIiIiIiIhdxB2/2NJFqNHUDvb7rnol4YcPH2b9+vWcOHGC8+fPk5SUxKxZs4z1FouFV199lR07dhgPZVm7di0DBgwwYq5cucILL7xAamqq8VCWhIQEevToYcSYzWaee+459u/fD8CvfvUr1qxZg6enZzNVV0RERERERFpKxNFdLV2Emj081m5fVa/p6GVlZQQFBfHqq6/Srl27ausTExNJSkpi9erVfPrpp/j4+DB58mQuXbpkxMTGxvLhhx+ybds29u3bx6VLl5gxYwaVlZVGzJw5c8jJyWHPnj2kpKSQk5PD3Llzm6GaIiIiIiIiIi2vXlfCx48fz/jx4wFYsGCB1TqLxUJycjKLFy82XmGSnJxMQEAAKSkpREVFUVJSws6dO0lKSmLMmDEAbNq0iUGDBnHo0CHGjh1Lbm4uBw8eZP/+/QwbNgyAdevWMXHiRPLy8vSeUBEREREREWn1mvxgttOnT1NQUEBERISxrF27dgwfPpzMzEwAsrOzuXbtmlVMz549CQwMNGKysrLo0KGDkYADhIeH4+HhYcSIiIiIiIiItGZNfjBbQUEBAD4+PlbLfXx8OH/+PACFhYW4uLjg7e1dLaawsNCI8fb2xmQyGetNJhNdunQxYmqSl5fX1CoYysvLm21bzaU56+dI32Vrjd2X9ugDjt7Ojlo+e/0+G/o9jtpeTWGPOjniv7e/1BxtoBlcIiIiUpNmezr6L5NnuDFN/eZlN7s5pqb4urbTXAc5+VDj/e4tzV4Hcc425f9CI/ZleXm5XfqAI7ezI/eDxuzThmpMH3DU9mose/UBe+zPpnC2/SoiIiKOo8nT0X19fQGqXa0uKioyro537dqVyspKiouLbxlTVFSExWIx1lssFoqLi6tdZRcRERERERFpjZqchPfq1QtfX1/S09ONZRUVFRw9etS4vzskJAQ3NzermPz8fHJzc42YsLAwSktLycrKMmKysrIoKyuzuk9cREREREREpLWq13T00tJSvv32WwCuX7/OuXPnyMnJwcvLCz8/P+bPn09CQgIBAQH069ePtWvX4uHhwdSpUwHo3Lkzs2fPZsWKFfj4+ODl5cXy5csJDg5m9OjRAAQGBjJu3DiWLFlCYmIiFouFJUuWMGHCBE0LFBEREREREadQryT8+PHjPPLII8bfcXFxxMXF8dhjj5GcnMyiRYsoLy9n6dKlmM1mQkND2bt3Lx07djQ+88orr+Di4kJUVBQVFRWMHDmSjRs34uLiYsRs2bKFmJgYpkyZAsDEiRNZs2ZNc9VVREREREREpEXVKwkfMWIEZrO51vUmk4nY2FhiY2NrjXF3dyc+Pp74+PhaY7y8vNi8eXN9iiQiIiIiIiLS6jT5nnARERERERERqR8l4SIiIiIiIiJ2oiRcRERERERExE6UhIuIiIiIiIjYiZJwERERERERETtREi4iIiIiIiJiJ0rCRUREREREROxESbiIiIiIiIiInSgJFxEREREREbETJeEiIiIiIiIidqIkXERERERERMROlISLiIiIiIiI2ImScBERERERERE7URIuIiIiIiIiYidKwkVERERERETsREm4iIiIiIiIiJ0oCRcRERERERGxEyXhIiIiIiIiInaiJFxERERERETETpSEi4iIiIiIiNiJknARERERERERO1ESLiIiIiIiImInSsJFRERERERE7ERJuIiIiIiIiIidKAkXERERERERsRMl4SIiIiIiIiJ2oiRcRERERERExE6UhIuIiIiIiIjYiZJwERERERERETtREi4iIiIiIiJiJ0rCRUREREREROxESbiIiIiIiIiInSgJFxEREREREbETJeEiIiIiIiIidqIkXERERERERMROHDIJ37p1K4MHD8bX15dRo0Zx5MiRli6SiIiI2JjGfxERuR04XBK+d+9eli1bxh/+8Af+/ve/ExYWxrRp0zh79mxLF01ERERsROO/iIjcLkxms9nS0oX4pbFjxxIcHMwbb7xhLBs6dCiRkZGsXLmyBUsmIiIitqLxX0REbhcOdSX86tWrZGdnExERYbU8IiKCzMzMFiqViIiI2JLGfxERuZ04VBJeXFxMZWUlPj4+Vst9fHwoLCxsoVKJiIiILWn8FxGR24lDJeFVTCaT1d8Wi6XaMhEREXEuGv9FROR24FBJuLe3Ny4uLtXOehcVFVU7Oy4iIiLOQeO/iIjcThwqCW/Tpg0hISGkp6dbLU9PT2fYsGEtVCoRERGxJY3/IiJyO3Ft6QLcbOHChcydO5fQ0FCGDRvG9u3b+fHHH4mKimrpoomIiIiNaPwXEZHbhUNdCQeYMmUKcXFxxMfHM2LECI4dO8bu3bvx9/dv9Da3bt3K4MGD8fX1ZdSoURw5cuSW8SdPnuRf/uVf6NatGwMGDGD16tVYLA71JrcGa0gbnD59Gk9Pz2r/HTx40I4lbj6HDx9m5syZDBgwAE9PT9599906P+NsfaChbeBsfeC1115jzJgx+Pn50bdvX2bMmMGpU6fq/Jwz9YPGtIGz9YMtW7YwfPhw/Pz88PPz4+GHH+bjjz++5WecqQ84uuYY/20x3n/++eeMGjUKX19f7r33XrZv396o+tVXQ+qQkZHBY489RmBgIN27d2f48OHs3LmzWkxNv+P/+q//cog61PffGXvuh4aUPy4ursbye3p68t///d+AffeBrY557Nn+Da2Do/0ObHXM5cj7wNF+B7Y87muu/eBwV8IB5syZw5w5c5plW3v37mXZsmUkJCQQHh7O1q1bmTZtGseOHcPPz69a/MWLF5k8eTLDhw/n008/JS8vj4ULF9K+fXt+97vfNUuZ7K2hbVAlNTWVgQMHGn97eXnZo7jNrqysjKCgIB577DHmzZtXZ7wz9oGGtkEVZ+kDn3/+OU899RRDhw7FYrHwyiuv8Oijj5KZmVlrnZytHzSmDao4Sz+46667ePHFF+nbty/Xr1/nz3/+M7NmzeLQoUNW9avibH2gNWjK+G+L8f77779n+vTpzJo1i82bN3Ps2DH+8Ic/4O3tTWRkZJPq2hx1yMrKIjg4mEWLFtGtWzf+4z/+g8WLF+Pu7s60adOsYo8dO2b12+3SpUuzl78xdahyq39n7LkfGlr+3/3udzz55JNWy5588klMJlO15xnYYx/Y4pjH3r+DhtbB0X4HtjjmcvR94Gi/A1sd9zXnfjCZzWanPqU/duxYgoODeeONN4xlQ4cOJTIykpUrV1aL37ZtG3/605/4r//6L9q1awdAfHw827dv59SpU63yKa0NbYPTp09z7733kp6ezpAhQ+xZVJvr0aMHa9asYdasWbXGOGMf+KX6tIEz9wGA0tJS/P39effdd5k4cWKNMc7eD+rTBs7eDwB69+7NypUra5zy7Ox9wNnYYrxfuXIlH374If/5n/9pfO53v/sd//f//l8OHDjQ4nWoyRNPPEFlZaVxJTAjI4NHHnmEb775Bm9v72Yv881sccxhz/3Q1H1w7tw5Bg8ezKZNm4wE0N77oEpzHfPY+3fQ0DrUpKV/B1Wa65irte0DR/odQPMd9zXnfnC46ejN6erVq2RnZxMREWG1PCIigszMzBo/k5WVxQMPPGA0Ptz4B/n8+fOcPn3apuW1hca0QZXZs2fTr18/JkyYwAcffGDLYjoUZ+sDTeGsfaC0tJTr16/j6elZa4yz94P6tEEVZ+wHlZWVpKamUlZWRlhYWI0xzt4HnImtxvusrKxq2xw7dizHjx/n2rVrLV6Hmly6dKnG3/Xo0aMJDAzk17/+NX//+9+bWtwa2eqYw177oTn2wc6dO+ncuTO//vWvq62zxz5oKEf7HTSXlvwdNJYj/Aaai6P9DprruK8594NTJ+HFxcVUVlZWmwbh4+NT7TUoVQoLC2uMr1rX2jSmDTp06MBLL73EW2+9xZ49exg5ciRRUVG899579ihyi3O2PtAYzt4Hli1bxqBBg2pNvsD5+0F92sAZ+8HJkyfp0aMHXbt2ZcmSJbzzzjsEBwfXGOvsfcCZ2Gq8ry3m559/pri4uLmKDzSuDjfbv38/n332GU888YSxrFu3brz22mvs3LmTnTt3EhAQQGRkJIcPH27O4gO2O+aw135o6j64fv067777LjNnzqRt27bGcnvug4ZytN9Bc2jp30FDOdJvoDk44u+guY77mnM/OOQ94c3t5imDFovlltMIa4qvaXlr0pA28Pb2trrfcciQIVy4cIHExERmzJhh03I6CmfsAw3hzH3g+eef59ixY+zfvx8XF5dbxjprP6hvGzhjPwgICCAjI4OSkhLS0tKYP38+f/vb3wgKCqox3ln7gLOyxXhv7z7Q0DpUOXbsGE8//TSrV68mNDTUWB4QEEBAQIDxd1hYGGfOnGH9+vU8+OCDzVfwX7DFMYc990Nj98Enn3zCuXPn+O1vf2u1vCX2QUM44u+gsRzpd1BfjvgbaApH+x0093Ffc+0Hp74S7u3tjYuLS7Wzl0VFRdXOYlTp2rVrjfFArZ9xZI1pg5qEhoby7bffNnfxHJKz9YHm4gx9IDY2ltTUVNLS0ujdu/ctY521HzSkDWrS2vtBmzZt6NOnD0OGDGHlypUMGjSIDRs21BjrrH3AGdlqvK8txtXVlTvvvLO5ig80bbw+evQo06ZNIzY2lqeeeqrO77LV79hWxxz22g9NLf+OHTsYNmwYAwYMqDPWUf4tdbTfQVM4yu+gObTUb6A5ONLvoLmP+5pzPzh1Et6mTRtCQkJIT0+3Wp6ens6wYcNq/ExYWBhHjx6loqLCKr579+706tXLpuW1hca0QU2+/PJLfH19m7t4DsnZ+kBzae19ICYmhpSUFNLS0rjnnnvqjHfGftDQNqhJa+8HN7t+/TpXr16tcZ0z9gFnZavxPiwsjEOHDlXb5pAhQ3Bzc2vxOsCNVwlNmzaN5557jgULFtTru2z1O7bVMYe99kNTyn/+/Hk++eSTalf/auMo/5Y62u+gsRzpd9AcWuo30FSO9DuwxXFfc+4Hp07CARYuXMiuXbt4++23yc3NJSYmhh9//NF4Eu6LL75o9dCAqVOn0q5dOxYsWMCpU6dIS0vj9ddfZ8GCBQ433aO+GtoGu3btYs+ePeTm5pKXl8f69evZunUr0dHRLVWFJiktLSUnJ4ecnByuX7/OuXPnyMnJ4ezZs8Dt0Qca2gbO1geeffZZdu3axdatW/H09KSgoICCggJKS0uNGGfvB41pA2frB3/60584cuQIp0+f5uTJk7z44ot8/vnnxpNbnb0PODtbjPdRUVH88MMPLFu2jNzcXN5++2127drFv/3bvzlEHTIyMpg2bRpRUVFMnz7d+F1XXb0B2LBhA3/729/45ptv+Prrr3nxxRf56KOPePrppx2iDvX5d8ae+6Gh5a/yzjvv4OHhweTJk6uts+c+sMUxj71/Bw2tg6P9DmxxzOXo+6CKo/wObHXc15z7wenvCZ8yZQoXLlwgPj6egoICBgwYwO7du/H39wfgxx9/5LvvvjPiO3fuzF//+leeffZZxowZg6enJwsXLrRZJ7eHhrYBwNq1azl79iwuLi707duXN998s9XeA3r8+HEeeeQR4++4uDji4uJ47LHHSE5Ovi36QEPbAJyrD2zduhWg2jscY2JiiI2NBZz/34LGtAE4Vz8oKCggOjqawsJCOnXqRHBwMCkpKYwdOxZw/j7g7Gwx3vfu3Zvdu3fz/PPPs337drp168bq1att8l7extRh165dXL58mfXr17N+/XpjuZ+fH19++SUA165d449//CPnz5/H3d3d2Ob48eMdog5Q978z9twPjSm/xWJh586dTJs2jfbt21fbpj33gS2Oeez9O2hoHRztd2CLYy5H3wfgWL8DWx33Ned+cPr3hIuIiIiIiIg4Cqefji4iIiIiIiLiKJSEi4iIiIiIiNiJknARERERERERO1ESLiIiIiIiImInSsJFRERERERE7ERJuIiIiIiIiIidKAkXERERERERsRMl4SIiIiIiIiJ28v8Aee3n9BQUyMcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x432 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(15, 6))\n",
    "plt.subplot(221)\n",
    "plt.hist(scm_1_sample, label='$SCM_1$', alpha=.5)\n",
    "plt.hist(scm_2_sample, label='$SCM_2$', alpha=.5)\n",
    "plt.legend()\n",
    "\n",
    "plt.subplot(222)\n",
    "plt.hist(scm_1_int_x_1_1, label='$SCM_1(x_1 = 1)$', alpha=.5)\n",
    "plt.hist(scm_2_int_x_1_1, label='$SCM_2(x_1 = 1)$', alpha=.5)\n",
    "plt.legend()\n",
    "\n",
    "plt.subplot(223)\n",
    "plt.hist(scm_1_int_x_2_1, label='$SCM_1(x_2 = 1)$', alpha=.5)\n",
    "plt.hist(scm_2_int_x_2_1, label='$SCM_2(x_2 = 1)$', alpha=.5)\n",
    "plt.legend()\n",
    "\n",
    "plt.subplot(224)\n",
    "plt.hist(scm_1_int_x_1_0, label='$SCM_1(x_2 = 0)$', alpha=.5)\n",
    "plt.hist(scm_2_int_x_1_0, label='$SCM_2(x_2 = 0)$', alpha=.5)\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c57b4cc",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:econml-dowhy-py38]",
   "language": "python",
   "name": "conda-env-econml-dowhy-py38-py"
  },
  "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}